21#ifndef OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
22#define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
24#include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
26#include <opm/material/common/MathToolbox.hpp>
27#include <opm/material/densead/Evaluation.hpp>
29#include <opm/output/data/Aquifer.hpp>
41template <
typename TypeTag>
48 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
49 using GridView = GetPropType<TypeTag, Properties::GridView>;
51 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
52 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
53 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
57 static constexpr int numEq = BlackoilIndices::numEq;
59 using Eval = DenseAd::Evaluation<Scalar, numEq>;
69 , cumulative_flux_(0.0)
70 , init_pressure_ (aquifer.numCells(), 0.0)
72 this->cell_to_aquifer_cell_idx_.resize(this->
simulator_.gridView().size(0), -1);
74 auto aquifer_on_process =
false;
75 for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
76 const auto* cell = aquifer.getCellPrt(idx);
79 const int compressed_idx = simulator.vanguard().compressedIndexForInterior(cell->global_index);
80 if (compressed_idx >= 0) {
81 this->cell_to_aquifer_cell_idx_[compressed_idx] = idx;
82 aquifer_on_process =
true;
86 if (aquifer_on_process) {
87 this->checkConnectsToReservoir();
94 result.flux_rate_ = 1.0;
95 result.cumulative_flux_ = 2.0;
96 result.init_pressure_ = {3.0, 4.0};
97 result.pressure_ = 5.0;
104 auto xaqPos = aquiferSoln.find(this->
aquiferID());
105 if (xaqPos == aquiferSoln.end())
108 if (this->connects_to_reservoir_) {
109 this->cumulative_flux_ = xaqPos->second.volume;
112 if (
const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
115 this->init_pressure_.resize(aqData->initPressure.size());
116 std::copy(aqData->initPressure.begin(),
117 aqData->initPressure.end(),
118 this->init_pressure_.begin());
121 this->solution_set_from_restart_ =
true;
129 this->pressure_ = this->calculateAquiferPressure();
130 this->flux_rate_ = this->calculateAquiferFluxRate();
131 this->cumulative_flux_ += this->flux_rate_ * this->
simulator_.timeStepSize();
136 data::AquiferData data;
138 data.pressure = this->pressure_;
139 data.fluxRate = this->flux_rate_;
140 data.volume = this->cumulative_flux_;
142 auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
143 aquNum->initPressure.resize(this->init_pressure_.size());
144 std::copy(this->init_pressure_.begin(),
145 this->init_pressure_.end(),
146 aquNum->initPressure.begin());
153 if (this->solution_set_from_restart_) {
157 this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
158 this->flux_rate_ = 0.;
159 this->cumulative_flux_ = 0.;
170 template<
class Serializer>
173 serializer(flux_rate_);
174 serializer(cumulative_flux_);
175 serializer(init_pressure_);
176 serializer(pressure_);
181 return this->flux_rate_ == rhs.flux_rate_ &&
182 this->cumulative_flux_ == rhs.cumulative_flux_ &&
183 this->init_pressure_ == rhs.init_pressure_ &&
184 this->pressure_ == rhs.pressure_;
189 return this->cumulative_flux_;
193 void checkConnectsToReservoir()
196 auto elemIt = std::find_if(this->
simulator_.gridView().template begin</*codim=*/0>(),
197 this->simulator_.gridView().template end</*codim=*/0>(),
198 [&elem_ctx,
this](
const auto& elem) ->
bool
200 elem_ctx.updateStencil(elem);
202 const auto cell_index = elem_ctx
203 .globalSpaceIndex(0, 0);
205 return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
208 assert ((elemIt != this->
simulator_.gridView().template end</*codim=*/0>())
209 &&
"Internal error locating numerical aquifer's connecting cell");
211 this->connects_to_reservoir_ =
212 elemIt->partitionType() == Dune::InteriorEntity;
215 Scalar calculateAquiferPressure()
const
217 auto capture = std::vector<Scalar>(this->init_pressure_.size(), 0.0);
218 return this->calculateAquiferPressure(capture);
221 Scalar calculateAquiferPressure(std::vector<Scalar>& cell_pressure)
const
223 Scalar sum_pressure_watervolume = 0.;
224 Scalar sum_watervolume = 0.;
227 const auto& gridView = this->
simulator_.gridView();
230 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
231 elem_ctx.updatePrimaryStencil(elem);
233 const std::size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
234 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
239 elem_ctx.updatePrimaryIntensiveQuantities(0);
240 const auto& iq0 = elem_ctx.intensiveQuantities(0, 0);
241 const auto& fs = iq0.fluidState();
246 const Scalar water_saturation = fs.saturation(this->
phaseIdx_()).value();
247 const Scalar porosity = iq0.porosity().value();
248 const Scalar volume = elem_ctx.dofTotalVolume(0, 0);
250 const Scalar water_pressure_reservoir = fs.pressure(this->
phaseIdx_()).value();
251 const Scalar water_volume = volume * porosity * water_saturation;
252 sum_pressure_watervolume += water_volume * water_pressure_reservoir;
253 sum_watervolume += water_volume;
255 cell_pressure[idx] = water_pressure_reservoir;
259 const auto& comm = this->
simulator_.vanguard().grid().comm();
260 comm.sum(&sum_pressure_watervolume, 1);
261 comm.sum(&sum_watervolume, 1);
264 comm.sum(cell_pressure.data(), cell_pressure.size());
266 return sum_pressure_watervolume / sum_watervolume;
269 template <
class ElemCtx>
270 Scalar getWaterFlux(
const ElemCtx& elem_ctx,
unsigned face_idx)
const
272 const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, 0);
273 const Scalar water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
277 Scalar calculateAquiferFluxRate()
const
279 Scalar aquifer_flux = 0.0;
281 if (! this->connects_to_reservoir_) {
286 const auto& gridView = this->
simulator_.gridView();
287 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
289 elem_ctx.updateStencil(elem);
291 const std::size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
292 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
298 const std::size_t num_interior_faces = elem_ctx.numInteriorFaces( 0);
300 const auto& stencil = elem_ctx.stencil(0);
303 for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
304 const auto& face = stencil.interiorFace(face_idx);
306 const std::size_t i = face.interiorIndex();
307 const std::size_t j = face.exteriorIndex();
310 const std::size_t J = stencil.globalSpaceIndex(j);
312 assert(stencil.globalSpaceIndex(i) == cell_index);
316 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
319 elem_ctx.updateAllIntensiveQuantities();
320 elem_ctx.updateAllExtensiveQuantities();
322 const Scalar water_flux = getWaterFlux(elem_ctx,face_idx);
323 const std::size_t up_id = water_flux >= 0.0 ? i : j;
324 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
325 const Scalar invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
326 const Scalar face_area = face.area();
327 aquifer_flux += water_flux * invB * face_area;
339 std::vector<Scalar> init_pressure_{};
341 bool solution_set_from_restart_ {
false};
342 bool connects_to_reservoir_ {
false};
345 std::vector<int> cell_to_aquifer_cell_idx_;
#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: AquiferInterface.hpp:35
const Simulator & simulator_
Definition: AquiferInterface.hpp:98
int phaseIdx_() const
Definition: AquiferInterface.hpp:88
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: AquiferInterface.hpp:39
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: AquiferInterface.hpp:37
int aquiferID() const
Definition: AquiferInterface.hpp:79
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: AquiferInterface.hpp:38
Definition: AquiferNumerical.hpp:43
void computeFaceAreaFraction(const std::vector< Scalar > &) override
Definition: AquiferNumerical.hpp:162
GetPropType< TypeTag, Properties::Indices > BlackoilIndices
Definition: AquiferNumerical.hpp:45
@ dimWorld
Definition: AquiferNumerical.hpp:55
MathToolbox< Eval > Toolbox
Definition: AquiferNumerical.hpp:60
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: AquiferNumerical.hpp:46
void serializeOp(Serializer &serializer)
Definition: AquiferNumerical.hpp:171
void initialSolutionApplied() override
Definition: AquiferNumerical.hpp:151
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: AquiferNumerical.hpp:53
GetPropType< TypeTag, Properties::ExtensiveQuantities > ExtensiveQuantities
Definition: AquiferNumerical.hpp:47
GetPropType< TypeTag, Properties::MaterialLaw > MaterialLaw
Definition: AquiferNumerical.hpp:51
AquiferNumerical(const SingleNumericalAquifer &aquifer, const Simulator &simulator)
Definition: AquiferNumerical.hpp:65
Scalar cumulativeFlux() const
Definition: AquiferNumerical.hpp:187
Scalar totalFaceArea() const override
Definition: AquiferNumerical.hpp:165
GetPropType< TypeTag, Properties::IntensiveQuantities > IntensiveQuantities
Definition: AquiferNumerical.hpp:50
@ numPhases
Definition: AquiferNumerical.hpp:56
void endTimeStep() override
Definition: AquiferNumerical.hpp:127
DenseAd::Evaluation< Scalar, numEq > Eval
Definition: AquiferNumerical.hpp:59
void addToSource(RateVector &, const unsigned, const unsigned) override
Definition: AquiferNumerical.hpp:125
data::AquiferData aquiferData() const override
Definition: AquiferNumerical.hpp:134
bool operator==(const AquiferNumerical &rhs) const
Definition: AquiferNumerical.hpp:179
void initFromRestart(const data::Aquifers &aquiferSoln) override
Definition: AquiferNumerical.hpp:102
static constexpr int numEq
Definition: AquiferNumerical.hpp:57
GetPropType< TypeTag, Properties::GridView > GridView
Definition: AquiferNumerical.hpp:49
void beginTimeStep() override
Definition: AquiferNumerical.hpp:124
static AquiferNumerical serializationTestObject(const Simulator &simulator)
Definition: AquiferNumerical.hpp:91
Definition: BlackoilPhases.hpp:27