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