22#ifndef OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
23#define OPM_AQUIFERANALYTICAL_HEADER_INCLUDED
25#include <dune/grid/common/partitionset.hh>
27#include <opm/common/ErrorMacros.hpp>
29#include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
31#include <opm/material/common/MathToolbox.hpp>
32#include <opm/material/densead/Evaluation.hpp>
33#include <opm/material/densead/Math.hpp>
34#include <opm/material/fluidstates/BlackOilFluidState.hpp>
39#include <opm/output/data/Aquifer.hpp>
54template <
typename TypeTag>
68 enum {
enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
69 enum {
enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
70 enum {
enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
75 static constexpr int numEq = BlackoilIndices::numEq;
83 BlackoilIndices::gasEnabled,
88 BlackoilIndices::numPhases>;
92 const std::vector<Aquancon::AquancCell>& connections,
102 assert (total_face_area.size() >=
static_cast<typename std::vector<Scalar>::size_type
>(this->aquiferID()));
104 const auto tfa = total_face_area[this->
aquiferID() - 1];
105 const auto eps_sqrt = std::sqrt(std::numeric_limits<Scalar>::epsilon());
107 if (tfa < eps_sqrt) {
112 this->faceArea_connected_.end(),
113 this->alphai_.begin(),
130 auto xaqPos = aquiferSoln.find(this->
aquiferID());
131 if (xaqPos == aquiferSoln.end())
137 this->
pa0_ = xaqPos->second.initPressure;
152 for (
const auto& elem : elements(this->
simulator_.gridView())) {
153 elemCtx.updatePrimaryStencil(elem);
155 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
160 elemCtx.updateIntensiveQuantities(0);
161 const auto& iq = elemCtx.intensiveQuantities(0, 0);
170 const unsigned cellIdx,
171 const unsigned timeIdx)
override
179 const auto& intQuants = model.intensiveQuantities(cellIdx, timeIdx);
185 rates[BlackoilIndices::conti0EqIdx +
compIdx_()]
186 += this->
Qai_[idx] / model.dofTotalVolume(cellIdx);
189 auto fs = intQuants.fluidState();
190 if (this->
Ta0_.has_value() && this->
Qai_[idx] > 0)
192 fs.setTemperature(this->
Ta0_.value());
193 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
194 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
195 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
196 paramCache.setRegionIndex(pvtRegionIdx);
197 paramCache.setMaxOilSat(this->
simulator_.problem().maxOilSaturation(cellIdx));
198 paramCache.updatePhase(fs, this->
phaseIdx_());
199 const auto& h = FluidSystem::enthalpy(fs, paramCache, this->
phaseIdx_());
202 rates[BlackoilIndices::contiEnergyEqIdx]
203 += this->
Qai_[idx] *fs.enthalpy(this->
phaseIdx_()) * FluidSystem::referenceDensity( this->
phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
213 template<
class Serializer>
241 return this->
simulator_.problem().gravity()[2];
247 return FluidSystem::oilCompIdx;
249 return FluidSystem::waterCompIdx;
265 this->
Qai_.resize(this->
size(), Scalar{0});
272 const auto& fs = intQuants.fluidState();
273 pressure_water.at(idx) = fs.pressure(this->
phaseIdx_());
280 const auto& fs = intQuants.fluidState();
281 pressure_water.at(idx) = fs.pressure(this->
phaseIdx_()).value();
292 const auto& gridView = this->
simulator_.vanguard().gridView();
293 for (std::size_t idx = 0; idx < this->
size(); ++idx) {
294 const auto global_index = this->
connections_[idx].global_index;
295 const int cell_index = this->
simulator_.vanguard().compressedIndex(global_index);
296 if (cell_index < 0) {
300 auto elemIt = gridView.template begin< 0>();
301 std::advance(elemIt, cell_index);
304 if (elemIt->partitionType() != Dune::InteriorEntity) {
312 FaceDir::DirEnum faceDirection;
315 const auto& elemMapper = this->
simulator_.model().dofMapper();
316 for (
const auto& elem : elements(gridView)) {
317 const unsigned cell_index = elemMapper.index(elem);
325 for (
const auto& intersection : intersections(gridView, elem)) {
327 if (! intersection.boundary()) {
331 switch (intersection.indexInInside()) {
333 faceDirection = FaceDir::XMinus;
336 faceDirection = FaceDir::XPlus;
339 faceDirection = FaceDir::YMinus;
342 faceDirection = FaceDir::YPlus;
345 faceDirection = FaceDir::ZMinus;
348 faceDirection = FaceDir::ZPlus;
351 OPM_THROW(std::logic_error,
352 "Internal error in initialization of aquifer.");
355 if (faceDirection == this->
connections_[idx].face_dir) {
369 const auto& gridView = this->
simulator_.vanguard().gridView();
370 for (std::size_t idx = 0; idx < this->
size(); ++idx) {
371 const int cell_index = this->
simulator_.vanguard()
373 if (cell_index < 0) {
377 auto elemIt = gridView.template begin< 0>();
378 std::advance(elemIt, cell_index);
381 if (elemIt->partitionType() != Dune::InteriorEntity) {
393 std::vector<Scalar> pw_aquifer;
394 Scalar water_pressure_reservoir;
397 const auto& gridView = this->
simulator_.gridView();
398 for (
const auto& elem : elements(gridView)) {
399 elemCtx.updatePrimaryStencil(elem);
401 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
406 elemCtx.updatePrimaryIntensiveQuantities(0);
407 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
408 const auto& fs = iq0.fluidState();
410 water_pressure_reservoir = fs.pressure(this->
phaseIdx_()).value();
411 const auto water_density = fs.density(this->
phaseIdx_());
416 pw_aquifer.push_back(this->
alphai_[idx] *
417 (water_pressure_reservoir - water_density.value()*gdz));
421 const auto& comm = this->
simulator_.vanguard().grid().comm();
424 vals[0] = std::accumulate(this->
alphai_.begin(), this->alphai_.end(),
Scalar{0});
425 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(),
Scalar{0});
429 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:192
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
Defines a type tags and some fundamental properties all models.
Declares the properties required by the black oil model.
Definition: AquiferAnalytical.hpp:56
void updateCellPressure(std::vector< Scalar > &pressure_water, const int idx, const IntensiveQuantities &intQuants)
Definition: AquiferAnalytical.hpp:276
Scalar total_face_area_
Definition: AquiferAnalytical.hpp:450
Scalar area_fraction_
Definition: AquiferAnalytical.hpp:451
Scalar gravity_() const
Definition: AquiferAnalytical.hpp:239
std::vector< int > cellToConnectionIdx_
Definition: AquiferAnalytical.hpp:436
Scalar Tc_
Definition: AquiferAnalytical.hpp:445
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: AquiferAnalytical.hpp:58
void addToSource(RateVector &rates, const unsigned cellIdx, const unsigned timeIdx) override
Definition: AquiferAnalytical.hpp:169
virtual void calculateAquiferCondition()=0
Scalar totalFaceArea() const override
Definition: AquiferAnalytical.hpp:123
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: AquiferAnalytical.hpp:64
void computeFaceAreaFraction(const std::vector< Scalar > &total_face_area) override
Definition: AquiferAnalytical.hpp:100
int compIdx_() const
Definition: AquiferAnalytical.hpp:244
bool solution_set_from_restart_
Definition: AquiferAnalytical.hpp:455
void initializeConnectionMappings()
Definition: AquiferAnalytical.hpp:284
@ has_disgas_in_water
Definition: AquiferAnalytical.hpp:71
virtual Scalar aquiferDepth() const =0
std::size_t size() const
Definition: AquiferAnalytical.hpp:208
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: AquiferAnalytical.hpp:60
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: AquiferAnalytical.hpp:63
void initFromRestart(const data::Aquifers &aquiferSoln) override
Definition: AquiferAnalytical.hpp:128
std::vector< Scalar > faceArea_connected_
Definition: AquiferAnalytical.hpp:435
void initialSolutionApplied() override
Definition: AquiferAnalytical.hpp:142
@ enableSaltPrecipitation
Definition: AquiferAnalytical.hpp:73
void beginTimeStep() override
Definition: AquiferAnalytical.hpp:147
std::vector< Scalar > alphai_
Definition: AquiferAnalytical.hpp:443
Scalar pa0_
Definition: AquiferAnalytical.hpp:446
DenseAd::Evaluation< Scalar, numEq > Eval
Definition: AquiferAnalytical.hpp:77
GetPropType< TypeTag, Properties::ElementMapper > ElementMapper
Definition: AquiferAnalytical.hpp:65
virtual void calculateInflowRate(int idx, const Simulator &simulator)=0
std::vector< Scalar > pressure_previous_
Definition: AquiferAnalytical.hpp:440
AquiferAnalytical(const int aqID, const std::vector< Aquancon::AquancCell > &connections, const Simulator &simulator)
Definition: AquiferAnalytical.hpp:91
void initQuantities()
Definition: AquiferAnalytical.hpp:252
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: AquiferAnalytical.hpp:59
Eval W_flux_
Definition: AquiferAnalytical.hpp:453
@ enableTemperature
Definition: AquiferAnalytical.hpp:67
virtual void calculateAquiferConstants()=0
std::vector< Eval > pressure_current_
Definition: AquiferAnalytical.hpp:441
@ enableBrine
Definition: AquiferAnalytical.hpp:69
static constexpr int numEq
Definition: AquiferAnalytical.hpp:75
std::optional< Scalar > Ta0_
Definition: AquiferAnalytical.hpp:447
std::vector< Scalar > cell_depth_
Definition: AquiferAnalytical.hpp:439
@ enableVapwat
Definition: AquiferAnalytical.hpp:70
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: AquiferAnalytical.hpp:61
void serializeOp(Serializer &serializer)
Definition: AquiferAnalytical.hpp:214
bool has_active_connection_on_proc_
Definition: AquiferAnalytical.hpp:456
const std::vector< Aquancon::AquancCell > connections_
Definition: AquiferAnalytical.hpp:432
std::vector< Eval > Qai_
Definition: AquiferAnalytical.hpp:442
@ enableEnergy
Definition: AquiferAnalytical.hpp:68
Scalar calculateReservoirEquilibrium()
Definition: AquiferAnalytical.hpp:390
virtual void assignRestartData(const data::AquiferData &xaq)=0
void initializeConnectionDepths()
Definition: AquiferAnalytical.hpp:365
BlackOilFluidState< Eval, FluidSystem, enableTemperature, enableEnergy, BlackoilIndices::gasEnabled, enableVapwat, enableBrine, enableSaltPrecipitation, has_disgas_in_water, BlackoilIndices::numPhases > FluidState
Definition: AquiferAnalytical.hpp:88
void updateCellPressure(std::vector< Eval > &pressure_water, const int idx, const IntensiveQuantities &intQuants)
Definition: AquiferAnalytical.hpp:268
GetPropType< TypeTag, Properties::Indices > BlackoilIndices
Definition: AquiferAnalytical.hpp:62
bool operator==(const AquiferAnalytical &rhs) const
Definition: AquiferAnalytical.hpp:223
Scalar rhow_
Definition: AquiferAnalytical.hpp:448
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: blackoilboundaryratevector.hh:39
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233