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