My Project
Loading...
Searching...
No Matches
AquiferNumerical.hpp
1/*
2 Copyright (C) 2020 Equinor ASA
3 Copyright (C) 2020 SINTEF Digital
4
5 This file is part of the Open Porous Media project (OPM).
6
7 OPM is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11
12 OPM is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with OPM. If not, see <http://www.gnu.org/licenses/>.
19*/
20
21#ifndef OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
22#define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
23
24#include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
25
26#include <opm/material/common/MathToolbox.hpp>
27#include <opm/material/densead/Evaluation.hpp>
28
29#include <opm/output/data/Aquifer.hpp>
30
31#include <opm/simulators/aquifers/AquiferInterface.hpp>
32#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
33
34#include <algorithm>
35#include <cassert>
36#include <cstddef>
37#include <vector>
38
39namespace Opm
40{
41template <typename TypeTag>
42class AquiferNumerical : public AquiferInterface<TypeTag>
43{
44public:
45 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
46 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
47 using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
48 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
49 using GridView = GetPropType<TypeTag, Properties::GridView>;
50 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
51 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
52 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
53
54 enum { dimWorld = GridView::dimensionworld };
55 enum { numPhases = FluidSystem::numPhases };
56 static constexpr int numEq = BlackoilIndices::numEq;
57
58 using Eval = DenseAd::Evaluation<double, numEq>;
59 using Toolbox = MathToolbox<Eval>;
60
61 using typename AquiferInterface<TypeTag>::RateVector;
62
63 // Constructor
64 AquiferNumerical(const SingleNumericalAquifer& aquifer,
65 const Simulator& simulator)
66 : AquiferInterface<TypeTag>(aquifer.id(), simulator)
67 , flux_rate_ (0.0)
68 , cumulative_flux_(0.0)
69 , init_pressure_ (aquifer.numCells(), 0.0)
70 {
71 this->cell_to_aquifer_cell_idx_.resize(this->simulator_.gridView().size(/*codim=*/0), -1);
72
73 auto aquifer_on_process = false;
74 for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
75 const auto* cell = aquifer.getCellPrt(idx);
76
77 // Due to parallelisation, the cell might not exist in the current process
78 const int compressed_idx = simulator.vanguard().compressedIndexForInterior(cell->global_index);
79 if (compressed_idx >= 0) {
80 this->cell_to_aquifer_cell_idx_[compressed_idx] = idx;
81 aquifer_on_process = true;
82 }
83 }
84
85 if (aquifer_on_process) {
86 this->checkConnectsToReservoir();
87 }
88 }
89
90 static AquiferNumerical serializationTestObject(const Simulator& simulator)
91 {
92 AquiferNumerical result({}, simulator);
93 result.flux_rate_ = 1.0;
94 result.cumulative_flux_ = 2.0;
95 result.init_pressure_ = {3.0, 4.0};
96 result.pressure_ = 5.0;
97
98 return result;
99 }
100
101 void initFromRestart(const data::Aquifers& aquiferSoln) override
102 {
103 auto xaqPos = aquiferSoln.find(this->aquiferID());
104 if (xaqPos == aquiferSoln.end())
105 return;
106
107 if (this->connects_to_reservoir_) {
108 this->cumulative_flux_ = xaqPos->second.volume;
109 }
110
111 if (const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
112 aqData != nullptr)
113 {
114 this->init_pressure_ = aqData->initPressure;
115 }
116
117 this->solution_set_from_restart_ = true;
118 }
119
120 void beginTimeStep() override {}
121 void addToSource(RateVector&, const unsigned, const unsigned) override {}
122
123 void endTimeStep() override
124 {
125 this->pressure_ = this->calculateAquiferPressure();
126 this->flux_rate_ = this->calculateAquiferFluxRate();
127 this->cumulative_flux_ += this->flux_rate_ * this->simulator_.timeStepSize();
128 }
129
130 data::AquiferData aquiferData() const override
131 {
132 data::AquiferData data;
133 data.aquiferID = this->aquiferID();
134 data.pressure = this->pressure_;
135 data.fluxRate = this->flux_rate_;
136 data.volume = this->cumulative_flux_;
137
138 auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
139 aquNum->initPressure = this->init_pressure_;
140
141 return data;
142 }
143
144 void initialSolutionApplied() override
145 {
146 if (this->solution_set_from_restart_) {
147 return;
148 }
149
150 this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
151 this->flux_rate_ = 0.;
152 this->cumulative_flux_ = 0.;
153 }
154
155 void computeFaceAreaFraction(const std::vector<double>& /*total_face_area*/) override
156 {}
157
158 double totalFaceArea() const override
159 {
160 return 1.0;
161 }
162
163 template<class Serializer>
164 void serializeOp(Serializer& serializer)
165 {
166 serializer(flux_rate_);
167 serializer(cumulative_flux_);
168 serializer(init_pressure_);
169 serializer(pressure_);
170 }
171
172 bool operator==(const AquiferNumerical& rhs) const
173 {
174 return this->flux_rate_ == rhs.flux_rate_ &&
175 this->cumulative_flux_ == rhs.cumulative_flux_ &&
176 this->init_pressure_ == rhs.init_pressure_ &&
177 this->pressure_ == rhs.pressure_;
178 }
179
180 double cumulativeFlux() const
181 {
182 return this->cumulative_flux_;
183 }
184
185private:
186 void checkConnectsToReservoir()
187 {
188 ElementContext elem_ctx(this->simulator_);
189 auto elemIt = std::find_if(this->simulator_.gridView().template begin</*codim=*/0>(),
190 this->simulator_.gridView().template end</*codim=*/0>(),
191 [&elem_ctx, this](const auto& elem) -> bool
192 {
193 elem_ctx.updateStencil(elem);
194
195 const auto cell_index = elem_ctx
196 .globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
197
198 return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
199 });
200
201 assert ((elemIt != this->simulator_.gridView().template end</*codim=*/0>())
202 && "Internal error locating numerical aquifer's connecting cell");
203
204 this->connects_to_reservoir_ =
205 elemIt->partitionType() == Dune::InteriorEntity;
206 }
207
208 double calculateAquiferPressure() const
209 {
210 auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
211 return this->calculateAquiferPressure(capture);
212 }
213
214 double calculateAquiferPressure(std::vector<double>& cell_pressure) const
215 {
216 double sum_pressure_watervolume = 0.;
217 double sum_watervolume = 0.;
218
219 ElementContext elem_ctx(this->simulator_);
220 const auto& gridView = this->simulator_.gridView();
221 OPM_BEGIN_PARALLEL_TRY_CATCH();
222
223 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
224 elem_ctx.updatePrimaryStencil(elem);
225
226 const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
227 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
228 if (idx < 0) {
229 continue;
230 }
231
232 elem_ctx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
233 const auto& iq0 = elem_ctx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
234 const auto& fs = iq0.fluidState();
235
236 // TODO: the porosity of the cells are still wrong for numerical aquifer cells
237 // Because the dofVolume still based on the grid information.
238 // The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
239 const double water_saturation = fs.saturation(this->phaseIdx_()).value();
240 const double porosity = iq0.porosity().value();
241 const double volume = elem_ctx.dofTotalVolume(0, 0);
242 // TODO: not sure we should use water pressure here
243 const double water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
244 const double water_volume = volume * porosity * water_saturation;
245 sum_pressure_watervolume += water_volume * water_pressure_reservoir;
246 sum_watervolume += water_volume;
247
248 cell_pressure[idx] = water_pressure_reservoir;
249 }
250 OPM_END_PARALLEL_TRY_CATCH("AquiferNumerical::calculateAquiferPressure() failed: ",
251 this->simulator_.vanguard().grid().comm());
252 const auto& comm = this->simulator_.vanguard().grid().comm();
253 comm.sum(&sum_pressure_watervolume, 1);
254 comm.sum(&sum_watervolume, 1);
255
256 // Ensure all processes have same notion of the aquifer cells' pressure values.
257 comm.sum(cell_pressure.data(), cell_pressure.size());
258
259 return sum_pressure_watervolume / sum_watervolume;
260 }
261
262 template <class ElemCtx>
263 double getWaterFlux(const ElemCtx& elem_ctx, unsigned face_idx) const
264 {
265 const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
266 const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
267 return water_flux;
268 }
269
270 double calculateAquiferFluxRate() const
271 {
272 double aquifer_flux = 0.0;
273
274 if (! this->connects_to_reservoir_) {
275 return aquifer_flux;
276 }
277
278 ElementContext elem_ctx(this->simulator_);
279 const auto& gridView = this->simulator_.gridView();
280 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
281 // elem_ctx.updatePrimaryStencil(elem);
282 elem_ctx.updateStencil(elem);
283
284 const std::size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
285 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
286 // we only need the first aquifer cell
287 if (idx != 0) {
288 continue;
289 }
290
291 const std::size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
292 // const auto &problem = elem_ctx.problem();
293 const auto& stencil = elem_ctx.stencil(0);
294 // const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
295
296 for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
297 const auto& face = stencil.interiorFace(face_idx);
298 // dof index
299 const std::size_t i = face.interiorIndex();
300 const std::size_t j = face.exteriorIndex();
301 // compressed index
302 // const std::size_t I = stencil.globalSpaceIndex(i);
303 const std::size_t J = stencil.globalSpaceIndex(j);
304
305 assert(stencil.globalSpaceIndex(i) == cell_index);
306
307 // we do not consider the flux within aquifer cells
308 // we only need the flux to the connections
309 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
310 continue;
311 }
312 elem_ctx.updateAllIntensiveQuantities();
313 elem_ctx.updateAllExtensiveQuantities();
314
315 const double water_flux = getWaterFlux(elem_ctx,face_idx);
316 const std::size_t up_id = water_flux >= 0.0 ? i : j;
317 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
318 const double invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
319 const double face_area = face.area();
320 aquifer_flux += water_flux * invB * face_area;
321 }
322
323 // we only need to handle the first aquifer cell, we can exit loop here
324 break;
325 }
326
327 return aquifer_flux;
328 }
329
330 double flux_rate_; // aquifer influx rate
331 double cumulative_flux_; // cumulative aquifer influx
332 std::vector<double> init_pressure_{};
333 double pressure_; // aquifer pressure
334 bool solution_set_from_restart_ {false};
335 bool connects_to_reservoir_ {false};
336
337 // TODO: maybe unordered_map can also do the work to save memory?
338 std::vector<int> cell_to_aquifer_cell_idx_;
339};
340
341} // namespace Opm
342
343#endif
Definition AquiferInterface.hpp:35
Definition AquiferNumerical.hpp:43
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition BlackoilPhases.hpp:27