22#ifndef OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
23#define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
25#include <opm/simulators/aquifers/AquiferInterface.hpp>
27#include <opm/common/utility/numeric/linearInterpolation.hpp>
28#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
30#include <opm/output/data/Aquifer.hpp>
32#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
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>
44#include <unordered_map>
49template <
typename TypeTag>
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>;
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>() };
67 static constexpr int numEq = BlackoilIndices::numEq;
68 using Scalar = double;
70 using Eval = DenseAd::Evaluation<double, numEq>;
72 using FluidState = BlackOilFluidState<Eval,
76 BlackoilIndices::gasEnabled,
79 enableSaltPrecipitation,
80 BlackoilIndices::numPhases>;
84 const std::vector<Aquancon::AquancCell>& connections,
85 const Simulator& ebosSimulator)
87 , connections_(connections)
96 void initFromRestart(
const data::Aquifers& aquiferSoln)
override
98 auto xaqPos = aquiferSoln.find(this->aquiferID());
99 if (xaqPos == aquiferSoln.end())
102 this->assignRestartData(xaqPos->second);
104 this->W_flux_ = xaqPos->second.volume;
105 this->pa0_ = xaqPos->second.initPressure;
106 this->solution_set_from_restart_ =
true;
109 void initialSolutionApplied()
override
114 void beginTimeStep()
override
116 ElementContext elemCtx(this->ebos_simulator_);
117 OPM_BEGIN_PARALLEL_TRY_CATCH();
119 for (
const auto& elem : elements(this->ebos_simulator_.gridView())) {
120 elemCtx.updatePrimaryStencil(elem);
122 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
123 const int idx = cellToConnectionIdx_[cellIdx];
127 elemCtx.updateIntensiveQuantities(0);
128 const auto& iq = elemCtx.intensiveQuantities(0, 0);
129 pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
132 OPM_END_PARALLEL_TRY_CATCH(
"AquiferAnalytical::beginTimeStep() failed: ",
133 this->ebos_simulator_.vanguard().grid().comm());
136 void addToSource(RateVector& rates,
137 const unsigned cellIdx,
138 const unsigned timeIdx)
override
140 const auto& model = this->ebos_simulator_.model();
142 const int idx = this->cellToConnectionIdx_[cellIdx];
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()");
152 this->updateCellPressure(this->pressure_current_, idx, *intQuantsPtr);
153 this->calculateInflowRate(idx, this->ebos_simulator_);
155 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
156 += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
158 if constexpr (enableEnergy) {
159 auto fs = intQuantsPtr->fluidState();
160 if (this->Ta0_.has_value() && this->Qai_[idx] > 0)
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);
172 rates[BlackoilIndices::contiEnergyEqIdx]
173 += this->Qai_[idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuantsPtr->pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
178 std::size_t size()
const
180 return this->connections_.size();
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;
190 Scalar gravity_()
const
192 return this->ebos_simulator_.problem().gravity()[2];
197 if (this->co2store_())
198 return FluidSystem::oilCompIdx;
200 return FluidSystem::waterCompIdx;
204 void initQuantities()
207 if (!this->solution_set_from_restart_) {
213 initializeConnections();
214 calculateAquiferCondition();
215 calculateAquiferConstants();
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});
222 void updateCellPressure(std::vector<Eval>& pressure_water,
224 const IntensiveQuantities& intQuants)
226 const auto& fs = intQuants.fluidState();
227 pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
230 void updateCellPressure(std::vector<Scalar>& pressure_water,
232 const IntensiveQuantities& intQuants)
234 const auto& fs = intQuants.fluidState();
235 pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
238 void initializeConnections()
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});
245 FaceDir::DirEnum faceDirection;
247 bool has_active_connection_on_proc =
false;
250 Scalar denom_face_areas{0};
251 this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(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< 0>();
258 std::advance(elemIt, cell_index);
261 if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
264 has_active_connection_on_proc =
true;
266 this->cellToConnectionIdx_[cell_index] = idx;
267 this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
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];
279 for (
const auto& intersection : intersections(gridView, elem)) {
281 if (!intersection.boundary())
284 int insideFaceIdx = intersection.indexInInside();
285 switch (insideFaceIdx) {
287 faceDirection = FaceDir::XMinus;
290 faceDirection = FaceDir::XPlus;
293 faceDirection = FaceDir::YMinus;
296 faceDirection = FaceDir::YPlus;
299 faceDirection = FaceDir::ZMinus;
302 faceDirection = FaceDir::ZPlus;
305 OPM_THROW(std::logic_error,
306 "Internal error in initialization of aquifer.");
310 if (faceDirection == this->connections_[idx].face_dir) {
311 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
315 denom_face_areas += this->faceArea_connected_.at(idx);
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) {
323 this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
325 : this->faceArea_connected_.at(idx) / denom_face_areas;
328 if (this->solution_set_from_restart_) {
329 this->rescaleProducedVolume(has_active_connection_on_proc);
333 void rescaleProducedVolume(
const bool has_active_connection_on_proc)
341 const auto this_area = has_active_connection_on_proc
342 ? std::accumulate(this->alphai_.begin(),
347 const auto tot_area = this->ebos_simulator_.vanguard()
348 .grid().comm().sum(this_area);
350 this->W_flux_ *= this_area / tot_area;
354 Scalar calculateReservoirEquilibrium()
357 std::vector<Scalar> pw_aquifer;
358 Scalar water_pressure_reservoir;
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);
365 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
366 const auto idx = this->cellToConnectionIdx_[cellIdx];
370 elemCtx.updatePrimaryIntensiveQuantities(0);
371 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
372 const auto& fs = iq0.fluidState();
374 water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
375 const auto water_density = fs.density(this->phaseIdx_());
378 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
380 pw_aquifer.push_back(this->alphai_[idx] *
381 (water_pressure_reservoir - water_density.value()*gdz));
385 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
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});
393 return vals[1] / vals[0];
396 const std::vector<Aquancon::AquancCell> connections_;
399 std::vector<Scalar> faceArea_connected_;
400 std::vector<int> cellToConnectionIdx_;
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_;
411 std::optional<Scalar> Ta0_{};
416 bool solution_set_from_restart_ {
false};
Definition: AquiferAnalytical.hpp:51
Definition: AquiferInterface.hpp:32
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27