22#ifndef OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
23#define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
25#include <opm/common/ErrorMacros.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>
52template <
typename TypeTag>
56 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
57 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
59 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
61 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
66 enum {
enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
67 enum {
enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
68 enum {
enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
73 static constexpr int numEq = BlackoilIndices::numEq;
81 BlackoilIndices::gasEnabled,
86 BlackoilIndices::numPhases>;
90 const std::vector<Aquancon::AquancCell>& connections,
104 assert (total_face_area.size() >=
static_cast<typename std::vector<Scalar>::size_type
>(this->aquiferID()));
106 const auto tfa = total_face_area[this->
aquiferID() - 1];
107 const auto eps_sqrt = std::sqrt(std::numeric_limits<Scalar>::epsilon());
109 if (tfa < eps_sqrt) {
114 this->faceArea_connected_.end(),
115 this->alphai_.begin(),
132 auto xaqPos = aquiferSoln.find(this->
aquiferID());
133 if (xaqPos == aquiferSoln.end())
139 this->
pa0_ = xaqPos->second.initPressure;
154 for (
const auto& elem : elements(this->
simulator_.gridView())) {
155 elemCtx.updatePrimaryStencil(elem);
157 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
162 elemCtx.updateIntensiveQuantities(0);
163 const auto& iq = elemCtx.intensiveQuantities(0, 0);
172 const unsigned cellIdx,
173 const unsigned timeIdx)
override
181 const auto& intQuants = model.intensiveQuantities(cellIdx, timeIdx);
187 rates[BlackoilIndices::conti0EqIdx +
compIdx_()]
188 += this->
Qai_[idx] / model.dofTotalVolume(cellIdx);
191 auto fs = intQuants.fluidState();
192 if (this->
Ta0_.has_value() && this->
Qai_[idx] > 0)
194 fs.setTemperature(this->
Ta0_.value());
195 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
196 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
197 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
198 paramCache.setRegionIndex(pvtRegionIdx);
199 paramCache.setMaxOilSat(this->
simulator_.problem().maxOilSaturation(cellIdx));
200 paramCache.updatePhase(fs, this->
phaseIdx_());
201 const auto& h = FluidSystem::enthalpy(fs, paramCache, this->
phaseIdx_());
204 rates[BlackoilIndices::contiEnergyEqIdx]
205 += this->
Qai_[idx] *fs.enthalpy(this->
phaseIdx_()) * FluidSystem::referenceDensity( this->
phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
215 template<
class Serializer>
243 return this->
simulator_.problem().gravity()[2];
249 return FluidSystem::oilCompIdx;
251 return FluidSystem::waterCompIdx;
267 this->
Qai_.resize(this->
size(), Scalar{0});
274 const auto& fs = intQuants.fluidState();
275 pressure_water.at(idx) = fs.pressure(this->
phaseIdx_());
282 const auto& fs = intQuants.fluidState();
283 pressure_water.at(idx) = fs.pressure(this->
phaseIdx_()).value();
294 const auto& gridView = this->
simulator_.vanguard().gridView();
295 for (std::size_t idx = 0; idx < this->
size(); ++idx) {
296 const auto global_index = this->
connections_[idx].global_index;
297 const int cell_index = this->
simulator_.vanguard().compressedIndex(global_index);
298 if (cell_index < 0) {
302 auto elemIt = gridView.template begin< 0>();
303 std::advance(elemIt, cell_index);
306 if (elemIt->partitionType() != Dune::InteriorEntity) {
314 FaceDir::DirEnum faceDirection;
317 const auto& elemMapper = this->
simulator_.model().dofMapper();
318 for (
const auto& elem : elements(gridView)) {
319 const unsigned cell_index = elemMapper.index(elem);
327 for (
const auto& intersection : intersections(gridView, elem)) {
329 if (! intersection.boundary()) {
333 switch (intersection.indexInInside()) {
335 faceDirection = FaceDir::XMinus;
338 faceDirection = FaceDir::XPlus;
341 faceDirection = FaceDir::YMinus;
344 faceDirection = FaceDir::YPlus;
347 faceDirection = FaceDir::ZMinus;
350 faceDirection = FaceDir::ZPlus;
353 OPM_THROW(std::logic_error,
354 "Internal error in initialization of aquifer.");
357 if (faceDirection == this->
connections_[idx].face_dir) {
371 const auto& gridView = this->
simulator_.vanguard().gridView();
372 for (std::size_t idx = 0; idx < this->
size(); ++idx) {
373 const int cell_index = this->
simulator_.vanguard()
375 if (cell_index < 0) {
379 auto elemIt = gridView.template begin< 0>();
380 std::advance(elemIt, cell_index);
383 if (elemIt->partitionType() != Dune::InteriorEntity) {
395 std::vector<Scalar> pw_aquifer;
396 Scalar water_pressure_reservoir;
399 const auto& gridView = this->
simulator_.gridView();
400 for (
const auto& elem : elements(gridView)) {
401 elemCtx.updatePrimaryStencil(elem);
403 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
408 elemCtx.updatePrimaryIntensiveQuantities(0);
409 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
410 const auto& fs = iq0.fluidState();
412 water_pressure_reservoir = fs.pressure(this->
phaseIdx_()).value();
413 const auto water_density = fs.density(this->
phaseIdx_());
418 pw_aquifer.push_back(this->
alphai_[idx] *
419 (water_pressure_reservoir - water_density.value()*gdz));
423 const auto& comm = this->
simulator_.vanguard().grid().comm();
426 vals[0] = std::accumulate(this->
alphai_.begin(), this->alphai_.end(),
Scalar{0});
427 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(),
Scalar{0});
431 return vals[1] / vals[0];
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:172
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:138
Definition: AquiferAnalytical.hpp:54
void updateCellPressure(std::vector< Scalar > &pressure_water, const int idx, const IntensiveQuantities &intQuants)
Definition: AquiferAnalytical.hpp:278
Scalar total_face_area_
Definition: AquiferAnalytical.hpp:452
Scalar area_fraction_
Definition: AquiferAnalytical.hpp:453
Scalar gravity_() const
Definition: AquiferAnalytical.hpp:241
std::vector< int > cellToConnectionIdx_
Definition: AquiferAnalytical.hpp:438
Scalar Tc_
Definition: AquiferAnalytical.hpp:447
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: AquiferAnalytical.hpp:56
@ enableVapwat
Definition: AquiferAnalytical.hpp:68
void addToSource(RateVector &rates, const unsigned cellIdx, const unsigned timeIdx) override
Definition: AquiferAnalytical.hpp:171
@ enableSaltPrecipitation
Definition: AquiferAnalytical.hpp:71
virtual void calculateAquiferCondition()=0
Scalar totalFaceArea() const override
Definition: AquiferAnalytical.hpp:125
@ enableEnergy
Definition: AquiferAnalytical.hpp:66
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: AquiferAnalytical.hpp:62
void computeFaceAreaFraction(const std::vector< Scalar > &total_face_area) override
Definition: AquiferAnalytical.hpp:102
int compIdx_() const
Definition: AquiferAnalytical.hpp:246
bool solution_set_from_restart_
Definition: AquiferAnalytical.hpp:457
void initializeConnectionMappings()
Definition: AquiferAnalytical.hpp:286
virtual Scalar aquiferDepth() const =0
std::size_t size() const
Definition: AquiferAnalytical.hpp:210
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: AquiferAnalytical.hpp:58
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: AquiferAnalytical.hpp:61
void initFromRestart(const data::Aquifers &aquiferSoln) override
Definition: AquiferAnalytical.hpp:130
std::vector< Scalar > faceArea_connected_
Definition: AquiferAnalytical.hpp:437
void initialSolutionApplied() override
Definition: AquiferAnalytical.hpp:144
void beginTimeStep() override
Definition: AquiferAnalytical.hpp:149
std::vector< Scalar > alphai_
Definition: AquiferAnalytical.hpp:445
Scalar pa0_
Definition: AquiferAnalytical.hpp:448
DenseAd::Evaluation< Scalar, numEq > Eval
Definition: AquiferAnalytical.hpp:75
GetPropType< TypeTag, Properties::ElementMapper > ElementMapper
Definition: AquiferAnalytical.hpp:63
virtual void calculateInflowRate(int idx, const Simulator &simulator)=0
std::vector< Scalar > pressure_previous_
Definition: AquiferAnalytical.hpp:442
AquiferAnalytical(const int aqID, const std::vector< Aquancon::AquancCell > &connections, const Simulator &simulator)
Definition: AquiferAnalytical.hpp:89
void initQuantities()
Definition: AquiferAnalytical.hpp:254
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: AquiferAnalytical.hpp:57
Eval W_flux_
Definition: AquiferAnalytical.hpp:455
virtual void calculateAquiferConstants()=0
@ has_disgas_in_water
Definition: AquiferAnalytical.hpp:69
std::vector< Eval > pressure_current_
Definition: AquiferAnalytical.hpp:443
static constexpr int numEq
Definition: AquiferAnalytical.hpp:73
@ enableBrine
Definition: AquiferAnalytical.hpp:67
std::optional< Scalar > Ta0_
Definition: AquiferAnalytical.hpp:449
std::vector< Scalar > cell_depth_
Definition: AquiferAnalytical.hpp:441
@ enableTemperature
Definition: AquiferAnalytical.hpp:65
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: AquiferAnalytical.hpp:59
void serializeOp(Serializer &serializer)
Definition: AquiferAnalytical.hpp:216
bool has_active_connection_on_proc_
Definition: AquiferAnalytical.hpp:458
const std::vector< Aquancon::AquancCell > connections_
Definition: AquiferAnalytical.hpp:434
std::vector< Eval > Qai_
Definition: AquiferAnalytical.hpp:444
Scalar calculateReservoirEquilibrium()
Definition: AquiferAnalytical.hpp:392
virtual ~AquiferAnalytical()
Definition: AquiferAnalytical.hpp:99
virtual void assignRestartData(const data::AquiferData &xaq)=0
void initializeConnectionDepths()
Definition: AquiferAnalytical.hpp:367
BlackOilFluidState< Eval, FluidSystem, enableTemperature, enableEnergy, BlackoilIndices::gasEnabled, enableVapwat, enableBrine, enableSaltPrecipitation, has_disgas_in_water, BlackoilIndices::numPhases > FluidState
Definition: AquiferAnalytical.hpp:86
void updateCellPressure(std::vector< Eval > &pressure_water, const int idx, const IntensiveQuantities &intQuants)
Definition: AquiferAnalytical.hpp:270
GetPropType< TypeTag, Properties::Indices > BlackoilIndices
Definition: AquiferAnalytical.hpp:60
bool operator==(const AquiferAnalytical &rhs) const
Definition: AquiferAnalytical.hpp:225
Scalar rhow_
Definition: AquiferAnalytical.hpp:450
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: BlackoilPhases.hpp:27