22#ifndef OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
23#define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
25#include <opm/common/utility/numeric/linearInterpolation.hpp>
27#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
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>
34#include <opm/models/blackoil/blackoilproperties.hh>
35#include <opm/models/utils/basicproperties.hh>
37#include <opm/output/data/Aquifer.hpp>
39#include <opm/simulators/aquifers/AquiferInterface.hpp>
40#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
48#include <unordered_map>
53template <
typename TypeTag>
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>;
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>() };
71 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
73 static constexpr int numEq = BlackoilIndices::numEq;
74 using Scalar = double;
76 using Eval = DenseAd::Evaluation<double, numEq>;
78 using FluidState = BlackOilFluidState<Eval,
82 BlackoilIndices::gasEnabled,
85 enableSaltPrecipitation,
87 BlackoilIndices::numPhases>;
91 const std::vector<Aquancon::AquancCell>& connections,
92 const Simulator& ebosSimulator)
94 , connections_(connections)
103 void initFromRestart(
const data::Aquifers& aquiferSoln)
override
105 auto xaqPos = aquiferSoln.find(this->aquiferID());
106 if (xaqPos == aquiferSoln.end())
109 this->assignRestartData(xaqPos->second);
111 this->W_flux_ = xaqPos->second.volume;
112 this->pa0_ = xaqPos->second.initPressure;
113 this->solution_set_from_restart_ =
true;
116 void initialSolutionApplied()
override
121 void beginTimeStep()
override
123 ElementContext elemCtx(this->ebos_simulator_);
124 OPM_BEGIN_PARALLEL_TRY_CATCH();
126 for (
const auto& elem : elements(this->ebos_simulator_.gridView())) {
127 elemCtx.updatePrimaryStencil(elem);
129 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
130 const int idx = cellToConnectionIdx_[cellIdx];
134 elemCtx.updateIntensiveQuantities(0);
135 const auto& iq = elemCtx.intensiveQuantities(0, 0);
136 pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
139 OPM_END_PARALLEL_TRY_CATCH(
"AquiferAnalytical::beginTimeStep() failed: ",
140 this->ebos_simulator_.vanguard().grid().comm());
143 void addToSource(RateVector& rates,
144 const unsigned cellIdx,
145 const unsigned timeIdx)
override
147 const auto& model = this->ebos_simulator_.model();
149 const int idx = this->cellToConnectionIdx_[cellIdx];
153 const auto& intQuants = model.intensiveQuantities(cellIdx, timeIdx);
156 this->updateCellPressure(this->pressure_current_, idx, intQuants);
157 this->calculateInflowRate(idx, this->ebos_simulator_);
159 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
160 += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
162 if constexpr (enableEnergy) {
163 auto fs = intQuants.fluidState();
164 if (this->Ta0_.has_value() && this->Qai_[idx] > 0)
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);
176 rates[BlackoilIndices::contiEnergyEqIdx]
177 += this->Qai_[idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
182 std::size_t size()
const
184 return this->connections_.size();
187 template<
class Serializer>
188 void serializeOp(Serializer& serializer)
190 serializer(pressure_previous_);
191 serializer(pressure_current_);
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_;
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;
213 Scalar gravity_()
const
215 return this->ebos_simulator_.problem().gravity()[2];
220 if (this->co2store_())
221 return FluidSystem::oilCompIdx;
223 return FluidSystem::waterCompIdx;
226 void initQuantities()
229 if (!this->solution_set_from_restart_) {
235 initializeConnections();
236 calculateAquiferCondition();
237 calculateAquiferConstants();
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});
244 void updateCellPressure(std::vector<Eval>& pressure_water,
246 const IntensiveQuantities& intQuants)
248 const auto& fs = intQuants.fluidState();
249 pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
252 void updateCellPressure(std::vector<Scalar>& pressure_water,
254 const IntensiveQuantities& intQuants)
256 const auto& fs = intQuants.fluidState();
257 pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
260 void initializeConnections()
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});
267 FaceDir::DirEnum faceDirection;
269 bool has_active_connection_on_proc =
false;
272 Scalar denom_face_areas{0};
273 this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(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< 0>();
280 std::advance(elemIt, cell_index);
283 if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
286 has_active_connection_on_proc =
true;
288 this->cellToConnectionIdx_[cell_index] = idx;
289 this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
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];
301 for (
const auto& intersection : intersections(gridView, elem)) {
303 if (!intersection.boundary())
306 int insideFaceIdx = intersection.indexInInside();
307 switch (insideFaceIdx) {
309 faceDirection = FaceDir::XMinus;
312 faceDirection = FaceDir::XPlus;
315 faceDirection = FaceDir::YMinus;
318 faceDirection = FaceDir::YPlus;
321 faceDirection = FaceDir::ZMinus;
324 faceDirection = FaceDir::ZPlus;
327 OPM_THROW(std::logic_error,
328 "Internal error in initialization of aquifer.");
332 if (faceDirection == this->connections_[idx].face_dir) {
333 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
337 denom_face_areas += this->faceArea_connected_.at(idx);
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) {
345 this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
347 : this->faceArea_connected_.at(idx) / denom_face_areas;
350 if (this->solution_set_from_restart_) {
351 this->rescaleProducedVolume(has_active_connection_on_proc);
355 void rescaleProducedVolume(
const bool has_active_connection_on_proc)
363 const auto this_area = has_active_connection_on_proc
364 ? std::accumulate(this->alphai_.begin(),
369 const auto tot_area = this->ebos_simulator_.vanguard()
370 .grid().comm().sum(this_area);
372 this->W_flux_ *= this_area / tot_area;
376 Scalar calculateReservoirEquilibrium()
379 std::vector<Scalar> pw_aquifer;
380 Scalar water_pressure_reservoir;
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);
387 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
388 const auto idx = this->cellToConnectionIdx_[cellIdx];
392 elemCtx.updatePrimaryIntensiveQuantities(0);
393 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
394 const auto& fs = iq0.fluidState();
396 water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
397 const auto water_density = fs.density(this->phaseIdx_());
400 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
402 pw_aquifer.push_back(this->alphai_[idx] *
403 (water_pressure_reservoir - water_density.value()*gdz));
407 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
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});
415 return vals[1] / vals[0];
418 const std::vector<Aquancon::AquancCell> connections_;
421 std::vector<Scalar> faceArea_connected_;
422 std::vector<int> cellToConnectionIdx_;
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_;
433 std::optional<Scalar> Ta0_{};
438 bool solution_set_from_restart_ {
false};
Definition: AquiferAnalytical.hpp:55
Definition: AquiferInterface.hpp:35
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27