21 #ifndef OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
22 #define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
24 #include <opm/output/data/Aquifer.hpp>
26 #include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
28 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
33 #include <unordered_map>
39 template <
typename TypeTag>
43 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
44 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
45 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
46 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
48 using GridView = GetPropType<TypeTag, Properties::GridView>;
49 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
51 enum { dimWorld = GridView::dimensionworld };
53 static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
55 static const int numEq = BlackoilIndices::numEq;
57 using Eval = DenseAd::Evaluation<double, numEq>;
58 using Toolbox = MathToolbox<Eval>;
62 const std::unordered_map<int, int>& cartesian_to_compressed,
63 const Simulator& ebos_simulator,
64 const int* global_cell)
66 , ebos_simulator_ (ebos_simulator)
68 , cumulative_flux_(0.0)
69 , global_cell_ (global_cell)
70 , init_pressure_ (aquifer.numCells(), 0.0)
72 this->cell_to_aquifer_cell_idx_.resize(this->ebos_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 auto search = cartesian_to_compressed.find(cell->global_index);
80 if (search != cartesian_to_compressed.end()) {
81 this->cell_to_aquifer_cell_idx_[search->second] = idx;
82 aquifer_on_process =
true;
86 if (aquifer_on_process) {
87 this->checkConnectsToReservoir();
91 void initFromRestart(
const data::Aquifers& aquiferSoln)
93 auto xaqPos = aquiferSoln.find(this->aquiferID());
94 if (xaqPos == aquiferSoln.end())
97 if (this->connects_to_reservoir_) {
98 this->cumulative_flux_ = xaqPos->second.volume;
101 if (
const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
104 this->init_pressure_ = aqData->initPressure;
107 this->solution_set_from_restart_ =
true;
112 this->pressure_ = this->calculateAquiferPressure();
113 this->flux_rate_ = this->calculateAquiferFluxRate();
114 this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
117 data::AquiferData aquiferData()
const
119 data::AquiferData data;
120 data.aquiferID = this->aquiferID();
121 data.pressure = this->pressure_;
122 data.fluxRate = this->flux_rate_;
123 data.volume = this->cumulative_flux_;
125 auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
126 aquNum->initPressure = this->init_pressure_;
131 void initialSolutionApplied()
133 if (this->solution_set_from_restart_) {
137 this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
138 this->flux_rate_ = 0.;
139 this->cumulative_flux_ = 0.;
142 int aquiferID()
const
144 return static_cast<int>(this->id_);
148 const std::size_t id_;
149 const Simulator& ebos_simulator_;
151 double cumulative_flux_;
152 const int* global_cell_;
153 std::vector<double> init_pressure_{};
155 bool solution_set_from_restart_ {
false};
156 bool connects_to_reservoir_ {
false};
159 std::vector<int> cell_to_aquifer_cell_idx_;
161 void checkConnectsToReservoir()
163 ElementContext elem_ctx(this->ebos_simulator_);
164 auto elemIt = std::find_if(this->ebos_simulator_.gridView().template begin</*codim=*/0>(),
165 this->ebos_simulator_.gridView().template end</*codim=*/0>(),
166 [&elem_ctx,
this](
const auto& elem) ->
bool
168 elem_ctx.updateStencil(elem);
170 const auto cell_index = elem_ctx
171 .globalSpaceIndex(0, 0);
173 return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
176 assert ((elemIt != this->ebos_simulator_.gridView().template end</*codim=*/0>())
177 &&
"Internal error locating numerical aquifer's connecting cell");
179 this->connects_to_reservoir_ =
180 elemIt->partitionType() == Dune::InteriorEntity;
183 double calculateAquiferPressure()
const
185 auto capture = std::vector<double>(this->init_pressure_.size(), 0.0);
186 return this->calculateAquiferPressure(capture);
189 double calculateAquiferPressure(std::vector<double>& cell_pressure)
const
191 double sum_pressure_watervolume = 0.;
192 double sum_watervolume = 0.;
194 ElementContext elem_ctx(this->ebos_simulator_);
195 const auto& gridView = this->ebos_simulator_.gridView();
196 auto elemIt = gridView.template begin<0>();
197 const auto& elemEndIt = gridView.template end<0>();
198 OPM_BEGIN_PARALLEL_TRY_CATCH();
200 for (; elemIt != elemEndIt; ++elemIt) {
201 const auto& elem = *elemIt;
202 if (elem.partitionType() != Dune::InteriorEntity) {
205 elem_ctx.updatePrimaryStencil(elem);
207 const size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
208 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
213 elem_ctx.updatePrimaryIntensiveQuantities(0);
214 const auto& iq0 = elem_ctx.intensiveQuantities(0, 0);
215 const auto& fs = iq0.fluidState();
220 const double water_saturation = fs.saturation(waterPhaseIdx).value();
221 const double porosity = iq0.porosity().value();
222 const double volume = elem_ctx.dofTotalVolume(0, 0);
224 const double water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
225 const double water_volume = volume * porosity * water_saturation;
226 sum_pressure_watervolume += water_volume * water_pressure_reservoir;
227 sum_watervolume += water_volume;
229 cell_pressure[idx] = water_pressure_reservoir;
231 OPM_END_PARALLEL_TRY_CATCH(
"AquiferNumerical::calculateAquiferPressure() failed: ", this->ebos_simulator_.vanguard().grid().comm());
232 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
233 comm.sum(&sum_pressure_watervolume, 1);
234 comm.sum(&sum_watervolume, 1);
237 comm.sum(cell_pressure.data(), cell_pressure.size());
239 return sum_pressure_watervolume / sum_watervolume;
242 double calculateAquiferFluxRate()
const
244 double aquifer_flux = 0.0;
246 if (! this->connects_to_reservoir_) {
250 ElementContext elem_ctx(this->ebos_simulator_);
251 const auto& gridView = this->ebos_simulator_.gridView();
252 auto elemIt = gridView.template begin<0>();
253 const auto& elemEndIt = gridView.template end<0>();
254 for (; elemIt != elemEndIt; ++elemIt) {
255 const auto& elem = *elemIt;
256 if (elem.partitionType() != Dune::InteriorEntity) {
260 elem_ctx.updateStencil(elem);
262 const std::size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
263 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
268 elem_ctx.updateAllIntensiveQuantities();
269 elem_ctx.updateAllExtensiveQuantities();
271 const std::size_t num_interior_faces = elem_ctx.numInteriorFaces( 0);
273 const auto& stencil = elem_ctx.stencil(0);
276 for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
277 const auto& face = stencil.interiorFace(face_idx);
279 const std::size_t i = face.interiorIndex();
280 const std::size_t j = face.exteriorIndex();
283 const std::size_t J = stencil.globalSpaceIndex(j);
285 assert(stencil.globalSpaceIndex(i) == cell_index);
289 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
292 const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, 0);
293 const double water_flux = Toolbox::value(exQuants.volumeFlux(waterPhaseIdx));
295 const std::size_t up_id = water_flux >= 0.0 ? i : j;
296 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
297 const double invB = Toolbox::value(intQuantsIn.fluidState().invB(waterPhaseIdx));
298 const double face_area = face.area();
299 aquifer_flux += water_flux * invB * face_area;
Definition: AquiferNumerical.hpp:41
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26