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.updatePhase(fs, this->
phaseIdx_());
198 const auto& h = FluidSystem::enthalpy(fs, paramCache, this->
phaseIdx_());
201 rates[BlackoilIndices::contiEnergyEqIdx]
202 += this->
Qai_[idx] *fs.enthalpy(this->
phaseIdx_()) * FluidSystem::referenceDensity( this->
phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
212 template<
class Serializer>
240 return this->
simulator_.problem().gravity()[2];
246 return FluidSystem::oilCompIdx;
248 return FluidSystem::waterCompIdx;
264 this->
Qai_.resize(this->
size(), Scalar{0});
271 const auto& fs = intQuants.fluidState();
272 pressure_water.at(idx) = fs.pressure(this->
phaseIdx_());
279 const auto& fs = intQuants.fluidState();
280 pressure_water.at(idx) = fs.pressure(this->
phaseIdx_()).value();
291 const auto& gridView = this->
simulator_.vanguard().gridView();
292 for (std::size_t idx = 0; idx < this->
size(); ++idx) {
293 const auto global_index = this->
connections_[idx].global_index;
294 const int cell_index = this->
simulator_.vanguard().compressedIndex(global_index);
295 if (cell_index < 0) {
299 auto elemIt = gridView.template begin< 0>();
300 std::advance(elemIt, cell_index);
303 if (elemIt->partitionType() != Dune::InteriorEntity) {
311 FaceDir::DirEnum faceDirection;
314 const auto& elemMapper = this->
simulator_.model().dofMapper();
315 for (
const auto& elem : elements(gridView)) {
316 const unsigned cell_index = elemMapper.index(elem);
324 for (
const auto& intersection : intersections(gridView, elem)) {
326 if (! intersection.boundary()) {
330 switch (intersection.indexInInside()) {
332 faceDirection = FaceDir::XMinus;
335 faceDirection = FaceDir::XPlus;
338 faceDirection = FaceDir::YMinus;
341 faceDirection = FaceDir::YPlus;
344 faceDirection = FaceDir::ZMinus;
347 faceDirection = FaceDir::ZPlus;
350 OPM_THROW(std::logic_error,
351 "Internal error in initialization of aquifer.");
354 if (faceDirection == this->
connections_[idx].face_dir) {
368 const auto& gridView = this->
simulator_.vanguard().gridView();
369 for (std::size_t idx = 0; idx < this->
size(); ++idx) {
370 const int cell_index = this->
simulator_.vanguard()
372 if (cell_index < 0) {
376 auto elemIt = gridView.template begin< 0>();
377 std::advance(elemIt, cell_index);
380 if (elemIt->partitionType() != Dune::InteriorEntity) {
392 std::vector<Scalar> pw_aquifer;
393 Scalar water_pressure_reservoir;
396 const auto& gridView = this->
simulator_.gridView();
397 for (
const auto& elem : elements(gridView)) {
398 elemCtx.updatePrimaryStencil(elem);
400 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
405 elemCtx.updatePrimaryIntensiveQuantities(0);
406 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
407 const auto& fs = iq0.fluidState();
409 water_pressure_reservoir = fs.pressure(this->
phaseIdx_()).value();
410 const auto water_density = fs.density(this->
phaseIdx_());
415 pw_aquifer.push_back(this->
alphai_[idx] *
416 (water_pressure_reservoir - water_density.value()*gdz));
420 const auto& comm = this->
simulator_.vanguard().grid().comm();
423 vals[0] = std::accumulate(this->
alphai_.begin(), this->alphai_.end(),
Scalar{0});
424 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(),
Scalar{0});
428 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:275
Scalar total_face_area_
Definition: AquiferAnalytical.hpp:449
Scalar area_fraction_
Definition: AquiferAnalytical.hpp:450
Scalar gravity_() const
Definition: AquiferAnalytical.hpp:238
std::vector< int > cellToConnectionIdx_
Definition: AquiferAnalytical.hpp:435
Scalar Tc_
Definition: AquiferAnalytical.hpp:444
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:243
bool solution_set_from_restart_
Definition: AquiferAnalytical.hpp:454
void initializeConnectionMappings()
Definition: AquiferAnalytical.hpp:283
@ has_disgas_in_water
Definition: AquiferAnalytical.hpp:71
virtual Scalar aquiferDepth() const =0
std::size_t size() const
Definition: AquiferAnalytical.hpp:207
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:434
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:442
Scalar pa0_
Definition: AquiferAnalytical.hpp:445
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:439
AquiferAnalytical(const int aqID, const std::vector< Aquancon::AquancCell > &connections, const Simulator &simulator)
Definition: AquiferAnalytical.hpp:91
void initQuantities()
Definition: AquiferAnalytical.hpp:251
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: AquiferAnalytical.hpp:59
Eval W_flux_
Definition: AquiferAnalytical.hpp:452
@ enableTemperature
Definition: AquiferAnalytical.hpp:67
virtual void calculateAquiferConstants()=0
std::vector< Eval > pressure_current_
Definition: AquiferAnalytical.hpp:440
@ enableBrine
Definition: AquiferAnalytical.hpp:69
static constexpr int numEq
Definition: AquiferAnalytical.hpp:75
std::optional< Scalar > Ta0_
Definition: AquiferAnalytical.hpp:446
std::vector< Scalar > cell_depth_
Definition: AquiferAnalytical.hpp:438
@ enableVapwat
Definition: AquiferAnalytical.hpp:70
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: AquiferAnalytical.hpp:61
void serializeOp(Serializer &serializer)
Definition: AquiferAnalytical.hpp:213
bool has_active_connection_on_proc_
Definition: AquiferAnalytical.hpp:455
const std::vector< Aquancon::AquancCell > connections_
Definition: AquiferAnalytical.hpp:431
std::vector< Eval > Qai_
Definition: AquiferAnalytical.hpp:441
@ enableEnergy
Definition: AquiferAnalytical.hpp:68
Scalar calculateReservoirEquilibrium()
Definition: AquiferAnalytical.hpp:389
virtual void assignRestartData(const data::AquiferData &xaq)=0
void initializeConnectionDepths()
Definition: AquiferAnalytical.hpp:364
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:267
GetPropType< TypeTag, Properties::Indices > BlackoilIndices
Definition: AquiferAnalytical.hpp:62
bool operator==(const AquiferAnalytical &rhs) const
Definition: AquiferAnalytical.hpp:222
Scalar rhow_
Definition: AquiferAnalytical.hpp:447
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