AquiferAnalytical.hpp
Go to the documentation of this file.
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 <dune/grid/common/partitionset.hh>
26
27#include <opm/common/ErrorMacros.hpp>
28
29#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
30
31#include <opm/material/common/MathToolbox.hpp>
32#include <opm/material/densead/Evaluation.hpp>
33#include <opm/material/densead/Math.hpp>
34#include <opm/material/fluidstates/BlackOilFluidState.hpp>
35
38
39#include <opm/output/data/Aquifer.hpp>
40
43
44#include <algorithm>
45#include <cmath>
46#include <cstddef>
47#include <limits>
48#include <numeric>
49#include <optional>
50#include <vector>
51
52namespace Opm
53{
54template <typename TypeTag>
55class AquiferAnalytical : public AquiferInterface<TypeTag>
56{
57public:
66
67 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
68 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
69 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
70 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
71 enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
72
73 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
74
75 static constexpr int numEq = BlackoilIndices::numEq;
76
77 using Eval = DenseAd::Evaluation<Scalar, /*size=*/numEq>;
78
79 using FluidState = BlackOilFluidState<Eval,
83 BlackoilIndices::gasEnabled,
88 BlackoilIndices::numPhases>;
89
90 // Constructor
91 AquiferAnalytical(const int aqID,
92 const std::vector<Aquancon::AquancCell>& connections,
93 const Simulator& simulator)
94 : AquiferInterface<TypeTag>(aqID, simulator)
95 , connections_(connections)
96 {
98 }
99
100 void computeFaceAreaFraction(const std::vector<Scalar>& total_face_area) override
101 {
102 assert (total_face_area.size() >= static_cast<typename std::vector<Scalar>::size_type>(this->aquiferID()));
103
104 const auto tfa = total_face_area[this->aquiferID() - 1];
105 const auto eps_sqrt = std::sqrt(std::numeric_limits<Scalar>::epsilon());
106
107 if (tfa < eps_sqrt) {
108 this->alphai_.assign(this->size(), Scalar{0});
109 }
110 else {
111 std::transform(this->faceArea_connected_.begin(),
112 this->faceArea_connected_.end(),
113 this->alphai_.begin(),
114 [tfa](const Scalar area)
115 {
116 return area / tfa;
117 });
118 }
119
120 this->area_fraction_ = this->totalFaceArea() / tfa;
121 }
122
123 Scalar totalFaceArea() const override
124 {
125 return this->total_face_area_;
126 }
127
128 void initFromRestart(const data::Aquifers& aquiferSoln) override
129 {
130 auto xaqPos = aquiferSoln.find(this->aquiferID());
131 if (xaqPos == aquiferSoln.end())
132 return;
133
134 this->assignRestartData(xaqPos->second);
135
136 this->W_flux_ = xaqPos->second.volume * this->area_fraction_;
137 this->pa0_ = xaqPos->second.initPressure;
138
139 this->solution_set_from_restart_ = true;
140 }
141
143 {
145 }
146
147 void beginTimeStep() override
148 {
149 ElementContext elemCtx(this->simulator_);
151
152 for (const auto& elem : elements(this->simulator_.gridView())) {
153 elemCtx.updatePrimaryStencil(elem);
154
155 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
156 const int idx = cellToConnectionIdx_[cellIdx];
157 if (idx < 0)
158 continue;
159
160 elemCtx.updateIntensiveQuantities(0);
161 const auto& iq = elemCtx.intensiveQuantities(0, 0);
162 pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
163 }
164
165 OPM_END_PARALLEL_TRY_CATCH("AquiferAnalytical::beginTimeStep() failed: ",
166 this->simulator_.vanguard().grid().comm());
167 }
168
170 const unsigned cellIdx,
171 const unsigned timeIdx) override
172 {
173 const auto& model = this->simulator_.model();
174
175 const int idx = this->cellToConnectionIdx_[cellIdx];
176 if (idx < 0)
177 return;
178
179 const auto& intQuants = model.intensiveQuantities(cellIdx, timeIdx);
180
181 // This is the pressure at td + dt
182 this->updateCellPressure(this->pressure_current_, idx, intQuants);
183 this->calculateInflowRate(idx, this->simulator_);
184
185 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
186 += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
187
188 if constexpr (enableEnergy) {
189 auto fs = intQuants.fluidState();
190 if (this->Ta0_.has_value() && this->Qai_[idx] > 0)
191 {
192 fs.setTemperature(this->Ta0_.value());
193 typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
194 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
195 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
196 paramCache.setRegionIndex(pvtRegionIdx);
197 paramCache.setMaxOilSat(this->simulator_.problem().maxOilSaturation(cellIdx));
198 paramCache.updatePhase(fs, this->phaseIdx_());
199 const auto& h = FluidSystem::enthalpy(fs, paramCache, this->phaseIdx_());
200 fs.setEnthalpy(this->phaseIdx_(), h);
201 }
202 rates[BlackoilIndices::contiEnergyEqIdx]
203 += this->Qai_[idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
204
205 }
206 }
207
208 std::size_t size() const
209 {
210 return this->connections_.size();
211 }
212
213 template<class Serializer>
214 void serializeOp(Serializer& serializer)
215 {
216 serializer(pressure_previous_);
217 serializer(pressure_current_);
218 serializer(Qai_);
219 serializer(rhow_);
220 serializer(W_flux_);
221 }
222
223 bool operator==(const AquiferAnalytical& rhs) const
224 {
225 return this->pressure_previous_ == rhs.pressure_previous_ &&
227 this->Qai_ == rhs.Qai_ &&
228 this->rhow_ == rhs.rhow_ &&
229 this->W_flux_ == rhs.W_flux_;
230 }
231
232protected:
233 virtual void assignRestartData(const data::AquiferData& xaq) = 0;
234 virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0;
235 virtual void calculateAquiferCondition() = 0;
236 virtual void calculateAquiferConstants() = 0;
237 virtual Scalar aquiferDepth() const = 0;
238
240 {
241 return this->simulator_.problem().gravity()[2];
242 }
243
244 int compIdx_() const
245 {
246 if (this->co2store_or_h2store_())
247 return FluidSystem::oilCompIdx;
248
249 return FluidSystem::waterCompIdx;
250 }
251
253 {
254 // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
255 if (! this->solution_set_from_restart_) {
256 W_flux_ = Scalar{0};
257 }
258
262
263 this->pressure_previous_.resize(this->size(), Scalar{0});
264 this->pressure_current_.resize(this->size(), Scalar{0});
265 this->Qai_.resize(this->size(), Scalar{0});
266 }
267
268 void updateCellPressure(std::vector<Eval>& pressure_water,
269 const int idx,
270 const IntensiveQuantities& intQuants)
271 {
272 const auto& fs = intQuants.fluidState();
273 pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
274 }
275
276 void updateCellPressure(std::vector<Scalar>& pressure_water,
277 const int idx,
278 const IntensiveQuantities& intQuants)
279 {
280 const auto& fs = intQuants.fluidState();
281 pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
282 }
283
285 {
286 this->alphai_.resize(this->size(), 1.0);
287 this->faceArea_connected_.resize(this->size(), Scalar{0});
288
289 // total_face_area_ is the sum of the areas connected to an aquifer
290 this->total_face_area_ = Scalar{0};
291 this->cellToConnectionIdx_.resize(this->simulator_.gridView().size(/*codim=*/0), -1);
292 const auto& gridView = this->simulator_.vanguard().gridView();
293 for (std::size_t idx = 0; idx < this->size(); ++idx) {
294 const auto global_index = this->connections_[idx].global_index;
295 const int cell_index = this->simulator_.vanguard().compressedIndex(global_index);
296 if (cell_index < 0) {
297 continue;
298 }
299
300 auto elemIt = gridView.template begin</*codim=*/ 0>();
301 std::advance(elemIt, cell_index);
302
303 // The global_index is not part of this grid
304 if (elemIt->partitionType() != Dune::InteriorEntity) {
305 continue;
306 }
307
308 this->cellToConnectionIdx_[cell_index] = idx;
309 }
310
311 // Translate the C face tag into the enum used by opm-parser's TransMult class
312 FaceDir::DirEnum faceDirection;
313
314 // Get areas for all connections
315 const auto& elemMapper = this->simulator_.model().dofMapper();
316 for (const auto& elem : elements(gridView)) {
317 const unsigned cell_index = elemMapper.index(elem);
318 const int idx = this->cellToConnectionIdx_[cell_index];
319
320 // Only deal with connections given by the aquifer
321 if (idx < 0) {
322 continue;
323 }
324
325 for (const auto& intersection : intersections(gridView, elem)) {
326 // Only deal with grid boundaries
327 if (! intersection.boundary()) {
328 continue;
329 }
330
331 switch (intersection.indexInInside()) {
332 case 0:
333 faceDirection = FaceDir::XMinus;
334 break;
335 case 1:
336 faceDirection = FaceDir::XPlus;
337 break;
338 case 2:
339 faceDirection = FaceDir::YMinus;
340 break;
341 case 3:
342 faceDirection = FaceDir::YPlus;
343 break;
344 case 4:
345 faceDirection = FaceDir::ZMinus;
346 break;
347 case 5:
348 faceDirection = FaceDir::ZPlus;
349 break;
350 default:
351 OPM_THROW(std::logic_error,
352 "Internal error in initialization of aquifer.");
353 }
354
355 if (faceDirection == this->connections_[idx].face_dir) {
356 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
357 break;
358 }
359 }
360
361 this->total_face_area_ += this->faceArea_connected_.at(idx);
362 }
363 }
364
366 {
367 this->cell_depth_.resize(this->size(), this->aquiferDepth());
368
369 const auto& gridView = this->simulator_.vanguard().gridView();
370 for (std::size_t idx = 0; idx < this->size(); ++idx) {
371 const int cell_index = this->simulator_.vanguard()
372 .compressedIndex(this->connections_[idx].global_index);
373 if (cell_index < 0) {
374 continue;
375 }
376
377 auto elemIt = gridView.template begin</*codim=*/ 0>();
378 std::advance(elemIt, cell_index);
379
380 // The global_index is not part of this grid
381 if (elemIt->partitionType() != Dune::InteriorEntity) {
382 continue;
383 }
384
385 this->cell_depth_.at(idx) = this->simulator_.vanguard().cellCenterDepth(cell_index);
386 }
387 }
388
389 // This function is for calculating the aquifer properties from equilibrium state with the reservoir
391 {
392 // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices
393 std::vector<Scalar> pw_aquifer;
394 Scalar water_pressure_reservoir;
395
396 ElementContext elemCtx(this->simulator_);
397 const auto& gridView = this->simulator_.gridView();
398 for (const auto& elem : elements(gridView)) {
399 elemCtx.updatePrimaryStencil(elem);
400
401 const auto cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
402 const auto idx = this->cellToConnectionIdx_[cellIdx];
403 if (idx < 0)
404 continue;
405
406 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
407 const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
408 const auto& fs = iq0.fluidState();
409
410 water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
411 const auto water_density = fs.density(this->phaseIdx_());
412
413 const auto gdz =
414 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
415
416 pw_aquifer.push_back(this->alphai_[idx] *
417 (water_pressure_reservoir - water_density.value()*gdz));
418 }
419
420 // We take the average of the calculated equilibrium pressures.
421 const auto& comm = this->simulator_.vanguard().grid().comm();
422
423 Scalar vals[2];
424 vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
425 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
426
427 comm.sum(vals, 2);
428
429 return vals[1] / vals[0];
430 }
431
432 const std::vector<Aquancon::AquancCell> connections_;
433
434 // Grid variables
435 std::vector<Scalar> faceArea_connected_;
436 std::vector<int> cellToConnectionIdx_;
437
438 // Quantities at each grid id
439 std::vector<Scalar> cell_depth_;
440 std::vector<Scalar> pressure_previous_;
441 std::vector<Eval> pressure_current_;
442 std::vector<Eval> Qai_;
443 std::vector<Scalar> alphai_;
444
445 Scalar Tc_{}; // Time constant
446 Scalar pa0_{}; // initial aquifer pressure
447 std::optional<Scalar> Ta0_{}; // initial aquifer temperature
449
452
454
457};
458
459} // namespace Opm
460
461#endif
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:192
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
Defines a type tags and some fundamental properties all models.
Declares the properties required by the black oil model.
Definition: AquiferAnalytical.hpp:56
void updateCellPressure(std::vector< Scalar > &pressure_water, const int idx, const IntensiveQuantities &intQuants)
Definition: AquiferAnalytical.hpp:276
Scalar total_face_area_
Definition: AquiferAnalytical.hpp:450
Scalar area_fraction_
Definition: AquiferAnalytical.hpp:451
Scalar gravity_() const
Definition: AquiferAnalytical.hpp:239
std::vector< int > cellToConnectionIdx_
Definition: AquiferAnalytical.hpp:436
Scalar Tc_
Definition: AquiferAnalytical.hpp:445
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: AquiferAnalytical.hpp:58
void addToSource(RateVector &rates, const unsigned cellIdx, const unsigned timeIdx) override
Definition: AquiferAnalytical.hpp:169
virtual void calculateAquiferCondition()=0
Scalar totalFaceArea() const override
Definition: AquiferAnalytical.hpp:123
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: AquiferAnalytical.hpp:64
void computeFaceAreaFraction(const std::vector< Scalar > &total_face_area) override
Definition: AquiferAnalytical.hpp:100
int compIdx_() const
Definition: AquiferAnalytical.hpp:244
bool solution_set_from_restart_
Definition: AquiferAnalytical.hpp:455
void initializeConnectionMappings()
Definition: AquiferAnalytical.hpp:284
@ has_disgas_in_water
Definition: AquiferAnalytical.hpp:71
virtual Scalar aquiferDepth() const =0
std::size_t size() const
Definition: AquiferAnalytical.hpp:208
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: AquiferAnalytical.hpp:60
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: AquiferAnalytical.hpp:63
void initFromRestart(const data::Aquifers &aquiferSoln) override
Definition: AquiferAnalytical.hpp:128
std::vector< Scalar > faceArea_connected_
Definition: AquiferAnalytical.hpp:435
void initialSolutionApplied() override
Definition: AquiferAnalytical.hpp:142
@ enableSaltPrecipitation
Definition: AquiferAnalytical.hpp:73
void beginTimeStep() override
Definition: AquiferAnalytical.hpp:147
std::vector< Scalar > alphai_
Definition: AquiferAnalytical.hpp:443
Scalar pa0_
Definition: AquiferAnalytical.hpp:446
DenseAd::Evaluation< Scalar, numEq > Eval
Definition: AquiferAnalytical.hpp:77
GetPropType< TypeTag, Properties::ElementMapper > ElementMapper
Definition: AquiferAnalytical.hpp:65
virtual void calculateInflowRate(int idx, const Simulator &simulator)=0
std::vector< Scalar > pressure_previous_
Definition: AquiferAnalytical.hpp:440
AquiferAnalytical(const int aqID, const std::vector< Aquancon::AquancCell > &connections, const Simulator &simulator)
Definition: AquiferAnalytical.hpp:91
void initQuantities()
Definition: AquiferAnalytical.hpp:252
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: AquiferAnalytical.hpp:59
Eval W_flux_
Definition: AquiferAnalytical.hpp:453
@ enableTemperature
Definition: AquiferAnalytical.hpp:67
virtual void calculateAquiferConstants()=0
std::vector< Eval > pressure_current_
Definition: AquiferAnalytical.hpp:441
@ enableBrine
Definition: AquiferAnalytical.hpp:69
static constexpr int numEq
Definition: AquiferAnalytical.hpp:75
std::optional< Scalar > Ta0_
Definition: AquiferAnalytical.hpp:447
std::vector< Scalar > cell_depth_
Definition: AquiferAnalytical.hpp:439
@ enableVapwat
Definition: AquiferAnalytical.hpp:70
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: AquiferAnalytical.hpp:61
void serializeOp(Serializer &serializer)
Definition: AquiferAnalytical.hpp:214
bool has_active_connection_on_proc_
Definition: AquiferAnalytical.hpp:456
const std::vector< Aquancon::AquancCell > connections_
Definition: AquiferAnalytical.hpp:432
std::vector< Eval > Qai_
Definition: AquiferAnalytical.hpp:442
@ enableEnergy
Definition: AquiferAnalytical.hpp:68
Scalar calculateReservoirEquilibrium()
Definition: AquiferAnalytical.hpp:390
virtual void assignRestartData(const data::AquiferData &xaq)=0
void initializeConnectionDepths()
Definition: AquiferAnalytical.hpp:365
BlackOilFluidState< Eval, FluidSystem, enableTemperature, enableEnergy, BlackoilIndices::gasEnabled, enableVapwat, enableBrine, enableSaltPrecipitation, has_disgas_in_water, BlackoilIndices::numPhases > FluidState
Definition: AquiferAnalytical.hpp:88
void updateCellPressure(std::vector< Eval > &pressure_water, const int idx, const IntensiveQuantities &intQuants)
Definition: AquiferAnalytical.hpp:268
GetPropType< TypeTag, Properties::Indices > BlackoilIndices
Definition: AquiferAnalytical.hpp:62
bool operator==(const AquiferAnalytical &rhs) const
Definition: AquiferAnalytical.hpp:223
Scalar rhow_
Definition: AquiferAnalytical.hpp:448
Definition: AquiferInterface.hpp:35
const Simulator & simulator_
Definition: AquiferInterface.hpp:98
int phaseIdx_() const
Definition: AquiferInterface.hpp:88
bool co2store_or_h2store_() const
Definition: AquiferInterface.hpp:82
int aquiferID() const
Definition: AquiferInterface.hpp:79
Definition: blackoilboundaryratevector.hh:39
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233