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>
31#include <opm/simulators/aquifers/AquiferInterface.hpp>
32#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
37#include <unordered_map>
43template <
typename TypeTag>
47 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
48 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
49 using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
50 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
51 using GridView = GetPropType<TypeTag, Properties::GridView>;
52 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
53 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
54 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
56 enum { dimWorld = GridView::dimensionworld };
57 enum { numPhases = FluidSystem::numPhases };
58 static constexpr int numEq = BlackoilIndices::numEq;
60 using Eval = DenseAd::Evaluation<double, numEq>;
61 using Toolbox = MathToolbox<Eval>;
63 using typename AquiferInterface<TypeTag>::RateVector;
67 const Simulator& ebos_simulator)
70 , cumulative_flux_(0.0)
71 , init_pressure_ (aquifer.numCells(), 0.0)
73 this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(0), -1);
75 auto aquifer_on_process =
false;
76 for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
77 const auto* cell = aquifer.getCellPrt(idx);
80 const int compressed_idx = ebos_simulator.vanguard().compressedIndexForInterior(cell->global_index);
81 if (compressed_idx >= 0) {
82 this->cell_to_aquifer_cell_idx_[compressed_idx] = idx;
83 aquifer_on_process =
true;
87 if (aquifer_on_process) {
88 this->checkConnectsToReservoir();
92 static AquiferNumerical serializationTestObject(
const Simulator& ebos_simulator)
95 result.flux_rate_ = 1.0;
96 result.cumulative_flux_ = 2.0;
97 result.init_pressure_ = {3.0, 4.0};
98 result.pressure_ = 5.0;
103 void initFromRestart(
const data::Aquifers& aquiferSoln)
override
105 auto xaqPos = aquiferSoln.find(this->aquiferID());
106 if (xaqPos == aquiferSoln.end())
109 if (this->connects_to_reservoir_) {
110 this->cumulative_flux_ = xaqPos->second.volume;
113 if (
const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
116 this->init_pressure_ = aqData->initPressure;
119 this->solution_set_from_restart_ =
true;
122 void beginTimeStep()
override {}
123 void addToSource(RateVector&,
const unsigned,
const unsigned)
override {}
125 void endTimeStep()
override
127 this->pressure_ = this->calculateAquiferPressure();
128 this->flux_rate_ = this->calculateAquiferFluxRate();
129 this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
132 data::AquiferData aquiferData()
const override
134 data::AquiferData data;
135 data.aquiferID = this->aquiferID();
136 data.pressure = this->pressure_;
137 data.fluxRate = this->flux_rate_;
138 data.volume = this->cumulative_flux_;
140 auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
141 aquNum->initPressure = this->init_pressure_;
146 void initialSolutionApplied()
override
148 if (this->solution_set_from_restart_) {
152 this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
153 this->flux_rate_ = 0.;
154 this->cumulative_flux_ = 0.;
157 template<
class Serializer>
158 void serializeOp(Serializer& serializer)
160 serializer(flux_rate_);
161 serializer(cumulative_flux_);
162 serializer(init_pressure_);
163 serializer(pressure_);
168 return this->flux_rate_ == rhs.flux_rate_ &&
169 this->cumulative_flux_ == rhs.cumulative_flux_ &&
170 this->init_pressure_ == rhs.init_pressure_ &&
171 this->pressure_ == rhs.pressure_;
174 double cumulativeFlux()
const
176 return this->cumulative_flux_;
180 void checkConnectsToReservoir()
182 ElementContext elem_ctx(this->ebos_simulator_);
183 auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin</*codim=*/0>(),
184 this->ebos_simulator_.gridView().template end</*codim=*/0>(),
185 [&elem_ctx,
this](
const auto& elem) ->
bool
187 elem_ctx.updateStencil(elem);
189 const auto cell_index = elem_ctx
190 .globalSpaceIndex(0, 0);
192 return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
195 assert ((elemIt != this->ebos_simulator_.gridView().template end</*codim=*/0>())
196 &&
"Internal error locating numerical aquifer's connecting cell");
198 this->connects_to_reservoir_ =
199 elemIt->partitionType() == Dune::InteriorEntity;
202 double calculateAquiferPressure()
const
204 auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
205 return this->calculateAquiferPressure(capture);
208 double calculateAquiferPressure(std::vector<double>& cell_pressure)
const
210 double sum_pressure_watervolume = 0.;
211 double sum_watervolume = 0.;
213 ElementContext elem_ctx(this->ebos_simulator_);
214 const auto& gridView = this->ebos_simulator_.gridView();
215 OPM_BEGIN_PARALLEL_TRY_CATCH();
217 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
218 elem_ctx.updatePrimaryStencil(elem);
220 const size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
221 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
226 elem_ctx.updatePrimaryIntensiveQuantities(0);
227 const auto& iq0 = elem_ctx.intensiveQuantities(0, 0);
228 const auto& fs = iq0.fluidState();
233 const double water_saturation = fs.saturation(this->phaseIdx_()).value();
234 const double porosity = iq0.porosity().value();
235 const double volume = elem_ctx.dofTotalVolume(0, 0);
237 const double water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
238 const double water_volume = volume * porosity * water_saturation;
239 sum_pressure_watervolume += water_volume * water_pressure_reservoir;
240 sum_watervolume += water_volume;
242 cell_pressure[idx] = water_pressure_reservoir;
244 OPM_END_PARALLEL_TRY_CATCH(
"AquiferNumerical::calculateAquiferPressure() failed: ", this->ebos_simulator_.vanguard().grid().comm());
245 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
246 comm.sum(&sum_pressure_watervolume, 1);
247 comm.sum(&sum_watervolume, 1);
250 comm.sum(cell_pressure.data(), cell_pressure.size());
252 return sum_pressure_watervolume / sum_watervolume;
255 template <
class ElemCtx>
256 double getWaterFlux(
const ElemCtx& elem_ctx,
unsigned face_idx)
const
258 const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, 0);
259 const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
263 double calculateAquiferFluxRate()
const
265 double aquifer_flux = 0.0;
267 if (! this->connects_to_reservoir_) {
271 ElementContext elem_ctx(this->ebos_simulator_);
272 const auto& gridView = this->ebos_simulator_.gridView();
273 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
275 elem_ctx.updateStencil(elem);
277 const std::size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
278 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
284 const std::size_t num_interior_faces = elem_ctx.numInteriorFaces( 0);
286 const auto& stencil = elem_ctx.stencil(0);
289 for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
290 const auto& face = stencil.interiorFace(face_idx);
292 const std::size_t i = face.interiorIndex();
293 const std::size_t j = face.exteriorIndex();
296 const std::size_t J = stencil.globalSpaceIndex(j);
298 assert(stencil.globalSpaceIndex(i) == cell_index);
302 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
305 elem_ctx.updateAllIntensiveQuantities();
306 elem_ctx.updateAllExtensiveQuantities();
308 const double water_flux = getWaterFlux(elem_ctx,face_idx);
309 const std::size_t up_id = water_flux >= 0.0 ? i : j;
310 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
311 const double invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
312 const double face_area = face.area();
313 aquifer_flux += water_flux * invB * face_area;
324 double cumulative_flux_;
325 std::vector<double> init_pressure_{};
327 bool solution_set_from_restart_ {
false};
328 bool connects_to_reservoir_ {
false};
331 std::vector<int> cell_to_aquifer_cell_idx_;
Definition: AquiferInterface.hpp:35
Definition: AquiferNumerical.hpp:45
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27