My Project
AquiferAnalytical.hpp
1/*
2 Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3 Copyright 2017 Statoil ASA.
4 Copyright 2017 IRIS
5
6 This file is part of the Open Porous Media project (OPM).
7
8 OPM is free software: you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12
13 OPM is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with OPM. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#ifndef OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
23#define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
24
25#include <opm/common/utility/numeric/linearInterpolation.hpp>
26
27#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
28
29#include <opm/material/common/MathToolbox.hpp>
30#include <opm/material/densead/Evaluation.hpp>
31#include <opm/material/densead/Math.hpp>
32#include <opm/material/fluidstates/BlackOilFluidState.hpp>
33
34#include <opm/models/blackoil/blackoilproperties.hh>
35#include <opm/models/utils/basicproperties.hh>
36
37#include <opm/output/data/Aquifer.hpp>
38
39#include <opm/simulators/aquifers/AquiferInterface.hpp>
40#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
41
42#include <algorithm>
43#include <cmath>
44#include <cstddef>
45#include <limits>
46#include <numeric>
47#include <optional>
48#include <unordered_map>
49#include <vector>
50
51namespace Opm
52{
53template <typename TypeTag>
54class AquiferAnalytical : public AquiferInterface<TypeTag>
55{
56public:
57 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
58 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
59 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
60 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
61 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
62 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
63 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
64
65 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
66 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
67 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
68 enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
69 enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
70
71 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
72
73 static constexpr int numEq = BlackoilIndices::numEq;
74 using Scalar = double;
75
76 using Eval = DenseAd::Evaluation<double, /*size=*/numEq>;
77
78 using FluidState = BlackOilFluidState<Eval,
79 FluidSystem,
80 enableTemperature,
81 enableEnergy,
82 BlackoilIndices::gasEnabled,
83 enableEvaporation,
84 enableBrine,
85 enableSaltPrecipitation,
86 has_disgas_in_water,
87 BlackoilIndices::numPhases>;
88
89 // Constructor
90 AquiferAnalytical(int aqID,
91 const std::vector<Aquancon::AquancCell>& connections,
92 const Simulator& ebosSimulator)
93 : AquiferInterface<TypeTag>(aqID, ebosSimulator)
94 , connections_(connections)
95 {
96 }
97
98 // Destructor
99 virtual ~AquiferAnalytical()
100 {
101 }
102
103 void initFromRestart(const data::Aquifers& aquiferSoln) override
104 {
105 auto xaqPos = aquiferSoln.find(this->aquiferID());
106 if (xaqPos == aquiferSoln.end())
107 return;
108
109 this->assignRestartData(xaqPos->second);
110
111 this->W_flux_ = xaqPos->second.volume;
112 this->pa0_ = xaqPos->second.initPressure;
113 this->solution_set_from_restart_ = true;
114 }
115
116 void initialSolutionApplied() override
117 {
118 initQuantities();
119 }
120
121 void beginTimeStep() override
122 {
123 ElementContext elemCtx(this->ebos_simulator_);
124 OPM_BEGIN_PARALLEL_TRY_CATCH();
125
126 for (const auto& elem : elements(this->ebos_simulator_.gridView())) {
127 elemCtx.updatePrimaryStencil(elem);
128
129 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
130 const int idx = cellToConnectionIdx_[cellIdx];
131 if (idx < 0)
132 continue;
133
134 elemCtx.updateIntensiveQuantities(0);
135 const auto& iq = elemCtx.intensiveQuantities(0, 0);
136 pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
137 }
138
139 OPM_END_PARALLEL_TRY_CATCH("AquiferAnalytical::beginTimeStep() failed: ",
140 this->ebos_simulator_.vanguard().grid().comm());
141 }
142
143 void addToSource(RateVector& rates,
144 const unsigned cellIdx,
145 const unsigned timeIdx) override
146 {
147 const auto& model = this->ebos_simulator_.model();
148
149 const int idx = this->cellToConnectionIdx_[cellIdx];
150 if (idx < 0)
151 return;
152
153 const auto& intQuants = model.intensiveQuantities(cellIdx, timeIdx);
154
155 // This is the pressure at td + dt
156 this->updateCellPressure(this->pressure_current_, idx, intQuants);
157 this->calculateInflowRate(idx, this->ebos_simulator_);
158
159 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
160 += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
161
162 if constexpr (enableEnergy) {
163 auto fs = intQuants.fluidState();
164 if (this->Ta0_.has_value() && this->Qai_[idx] > 0)
165 {
166 fs.setTemperature(this->Ta0_.value());
167 typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
168 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
169 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
170 paramCache.setRegionIndex(pvtRegionIdx);
171 paramCache.setMaxOilSat(this->ebos_simulator_.problem().maxOilSaturation(cellIdx));
172 paramCache.updatePhase(fs, this->phaseIdx_());
173 const auto& h = FluidSystem::enthalpy(fs, paramCache, this->phaseIdx_());
174 fs.setEnthalpy(this->phaseIdx_(), h);
175 }
176 rates[BlackoilIndices::contiEnergyEqIdx]
177 += this->Qai_[idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
178
179 }
180 }
181
182 std::size_t size() const
183 {
184 return this->connections_.size();
185 }
186
187 template<class Serializer>
188 void serializeOp(Serializer& serializer)
189 {
190 serializer(pressure_previous_);
191 serializer(pressure_current_);
192 serializer(Qai_);
193 serializer(rhow_);
194 serializer(W_flux_);
195 }
196
197 bool operator==(const AquiferAnalytical& rhs) const
198 {
199 return this->pressure_previous_ == rhs.pressure_previous_ &&
200 this->pressure_current_ == rhs.pressure_current_ &&
201 this->Qai_ == rhs.Qai_ &&
202 this->rhow_ == rhs.rhow_ &&
203 this->W_flux_ == rhs.W_flux_;
204 }
205
206protected:
207 virtual void assignRestartData(const data::AquiferData& xaq) = 0;
208 virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
209 virtual void calculateAquiferCondition() = 0;
210 virtual void calculateAquiferConstants() = 0;
211 virtual Scalar aquiferDepth() const = 0;
212
213 Scalar gravity_() const
214 {
215 return this->ebos_simulator_.problem().gravity()[2];
216 }
217
218 int compIdx_() const
219 {
220 if (this->co2store_())
221 return FluidSystem::oilCompIdx;
222
223 return FluidSystem::waterCompIdx;
224 }
225
226 void initQuantities()
227 {
228 // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
229 if (!this->solution_set_from_restart_) {
230 W_flux_ = Scalar{0};
231 }
232
233 // We next get our connections to the aquifer and initialize these quantities using the initialize_connections
234 // function
235 initializeConnections();
236 calculateAquiferCondition();
237 calculateAquiferConstants();
238
239 pressure_previous_.resize(this->connections_.size(), Scalar{0});
240 pressure_current_.resize(this->connections_.size(), Scalar{0});
241 Qai_.resize(this->connections_.size(), Scalar{0});
242 }
243
244 void updateCellPressure(std::vector<Eval>& pressure_water,
245 const int idx,
246 const IntensiveQuantities& intQuants)
247 {
248 const auto& fs = intQuants.fluidState();
249 pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
250 }
251
252 void updateCellPressure(std::vector<Scalar>& pressure_water,
253 const int idx,
254 const IntensiveQuantities& intQuants)
255 {
256 const auto& fs = intQuants.fluidState();
257 pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
258 }
259
260 void initializeConnections()
261 {
262 this->cell_depth_.resize(this->size(), this->aquiferDepth());
263 this->alphai_.resize(this->size(), 1.0);
264 this->faceArea_connected_.resize(this->size(), Scalar{0});
265
266 // Translate the C face tag into the enum used by opm-parser's TransMult class
267 FaceDir::DirEnum faceDirection;
268
269 bool has_active_connection_on_proc = false;
270
271 // denom_face_areas is the sum of the areas connected to an aquifer
272 Scalar denom_face_areas{0};
273 this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
274 const auto& gridView = this->ebos_simulator_.vanguard().gridView();
275 for (std::size_t idx = 0; idx < this->size(); ++idx) {
276 const auto global_index = this->connections_[idx].global_index;
277 const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
278 auto elemIt = gridView.template begin</*codim=*/ 0>();
279 if (cell_index > 0)
280 std::advance(elemIt, cell_index);
281
282 //the global_index is not part of this grid
283 if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
284 continue;
285
286 has_active_connection_on_proc = true;
287
288 this->cellToConnectionIdx_[cell_index] = idx;
289 this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
290 }
291 // get areas for all connections
292 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
293 for (const auto& elem : elements(gridView)) {
294 unsigned cell_index = elemMapper.index(elem);
295 int idx = this->cellToConnectionIdx_[cell_index];
296
297 // only deal with connections given by the aquifer
298 if( idx < 0)
299 continue;
300
301 for (const auto& intersection : intersections(gridView, elem)) {
302 // only deal with grid boundaries
303 if (!intersection.boundary())
304 continue;
305
306 int insideFaceIdx = intersection.indexInInside();
307 switch (insideFaceIdx) {
308 case 0:
309 faceDirection = FaceDir::XMinus;
310 break;
311 case 1:
312 faceDirection = FaceDir::XPlus;
313 break;
314 case 2:
315 faceDirection = FaceDir::YMinus;
316 break;
317 case 3:
318 faceDirection = FaceDir::YPlus;
319 break;
320 case 4:
321 faceDirection = FaceDir::ZMinus;
322 break;
323 case 5:
324 faceDirection = FaceDir::ZPlus;
325 break;
326 default:
327 OPM_THROW(std::logic_error,
328 "Internal error in initialization of aquifer.");
329 }
330
331
332 if (faceDirection == this->connections_[idx].face_dir) {
333 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
334 break;
335 }
336 }
337 denom_face_areas += this->faceArea_connected_.at(idx);
338 }
339
340 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
341 comm.sum(&denom_face_areas, 1);
342 const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
343 for (std::size_t idx = 0; idx < this->size(); ++idx) {
344 // Protect against division by zero NaNs.
345 this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
346 ? Scalar{0}
347 : this->faceArea_connected_.at(idx) / denom_face_areas;
348 }
349
350 if (this->solution_set_from_restart_) {
351 this->rescaleProducedVolume(has_active_connection_on_proc);
352 }
353 }
354
355 void rescaleProducedVolume(const bool has_active_connection_on_proc)
356 {
357 // Needed in parallel restart to approximate influence of aquifer
358 // being "owned" by a subset of the parallel processes. If the
359 // aquifer is fully owned by a single process--i.e., if all cells
360 // connecting to the aquifer are on a single process--then this_area
361 // is tot_area on that process and zero elsewhere.
362
363 const auto this_area = has_active_connection_on_proc
364 ? std::accumulate(this->alphai_.begin(),
365 this->alphai_.end(),
366 Scalar{0})
367 : Scalar{0};
368
369 const auto tot_area = this->ebos_simulator_.vanguard()
370 .grid().comm().sum(this_area);
371
372 this->W_flux_ *= this_area / tot_area;
373 }
374
375 // This function is for calculating the aquifer properties from equilibrium state with the reservoir
376 Scalar calculateReservoirEquilibrium()
377 {
378 // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
379 std::vector<Scalar> pw_aquifer;
380 Scalar water_pressure_reservoir;
381
382 ElementContext elemCtx(this->ebos_simulator_);
383 const auto& gridView = this->ebos_simulator_.gridView();
384 for (const auto& elem : elements(gridView)) {
385 elemCtx.updatePrimaryStencil(elem);
386
387 const auto cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
388 const auto idx = this->cellToConnectionIdx_[cellIdx];
389 if (idx < 0)
390 continue;
391
392 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
393 const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
394 const auto& fs = iq0.fluidState();
395
396 water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
397 const auto water_density = fs.density(this->phaseIdx_());
398
399 const auto gdz =
400 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
401
402 pw_aquifer.push_back(this->alphai_[idx] *
403 (water_pressure_reservoir - water_density.value()*gdz));
404 }
405
406 // We take the average of the calculated equilibrium pressures.
407 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
408
409 Scalar vals[2];
410 vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
411 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
412
413 comm.sum(vals, 2);
414
415 return vals[1] / vals[0];
416 }
417
418 const std::vector<Aquancon::AquancCell> connections_;
419
420 // Grid variables
421 std::vector<Scalar> faceArea_connected_;
422 std::vector<int> cellToConnectionIdx_;
423
424 // Quantities at each grid id
425 std::vector<Scalar> cell_depth_;
426 std::vector<Scalar> pressure_previous_;
427 std::vector<Eval> pressure_current_;
428 std::vector<Eval> Qai_;
429 std::vector<Scalar> alphai_;
430
431 Scalar Tc_{}; // Time constant
432 Scalar pa0_{}; // initial aquifer pressure
433 std::optional<Scalar> Ta0_{}; // initial aquifer temperature
434 Scalar rhow_{};
435
436 Eval W_flux_;
437
438 bool solution_set_from_restart_ {false};
439};
440
441} // namespace Opm
442
443#endif
Definition: AquiferAnalytical.hpp:55
Definition: AquiferInterface.hpp:35
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27