22 #ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED
23 #define OPM_AQUIFERINTERFACE_HEADER_INCLUDED
25 #include <opm/common/utility/numeric/linearInterpolation.hpp>
26 #include <opm/input/eclipse/EclipseState/Aquifer/Aquancon.hpp>
28 #include <opm/output/data/Aquifer.hpp>
30 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
32 #include <opm/material/common/MathToolbox.hpp>
33 #include <opm/material/densead/Evaluation.hpp>
34 #include <opm/material/densead/Math.hpp>
35 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
42 #include <unordered_map>
47 template <
typename TypeTag>
51 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
52 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
53 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
54 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
55 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
56 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
57 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
59 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
60 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
61 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
62 enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
63 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
65 static const int numEq = BlackoilIndices::numEq;
66 typedef double Scalar;
68 typedef DenseAd::Evaluation<double, numEq> Eval;
70 typedef BlackOilFluidState<Eval,
74 BlackoilIndices::gasEnabled,
77 enableSaltPrecipitation,
78 BlackoilIndices::numPhases>
83 const std::vector<Aquancon::AquancCell>& connections,
84 const Simulator& ebosSimulator)
86 , connections_(connections)
87 , ebos_simulator_(ebosSimulator)
96 void initFromRestart(
const data::Aquifers& aquiferSoln)
98 auto xaqPos = aquiferSoln.find(this->aquiferID());
99 if (xaqPos == aquiferSoln.end())
102 this->assignRestartData(xaqPos->second);
104 this->W_flux_ = xaqPos->second.volume;
105 this->pa0_ = xaqPos->second.initPressure;
106 this->solution_set_from_restart_ =
true;
109 void initialSolutionApplied()
116 ElementContext elemCtx(ebos_simulator_);
117 auto elemIt = ebos_simulator_.gridView().template begin<0>();
118 const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
119 OPM_BEGIN_PARALLEL_TRY_CATCH();
121 for (; elemIt != elemEndIt; ++elemIt) {
122 const auto& elem = *elemIt;
124 elemCtx.updatePrimaryStencil(elem);
126 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
127 const int idx = cellToConnectionIdx_[cellIdx];
131 elemCtx.updateIntensiveQuantities(0);
132 const auto& iq = elemCtx.intensiveQuantities(0, 0);
133 pressure_previous_[idx] = getValue(iq.fluidState().pressure(phaseIdx_()));
136 OPM_END_PARALLEL_TRY_CATCH(
"AquiferInterface::beginTimeStep() failed: ", ebos_simulator_.vanguard().grid().comm());
139 template <
class Context>
140 void addToSource(RateVector& rates,
141 const Context& context,
142 const unsigned spaceIdx,
143 const unsigned timeIdx)
145 const unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
147 const int idx = this->cellToConnectionIdx_[cellIdx];
154 const auto& intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
157 this->updateCellPressure(this->pressure_current_, idx, intQuants);
158 this->calculateInflowRate(idx, context.simulator());
160 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
161 += this->Qai_[idx] / context.dofVolume(spaceIdx, timeIdx);
164 std::size_t size()
const {
165 return this->connections_.size();
168 int aquiferID()
const {
return this->aquiferID_; }
171 inline Scalar gravity_()
const
173 return ebos_simulator_.problem().gravity()[2];
176 inline bool co2store_()
const
178 return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
181 inline int phaseIdx_()
const
184 return FluidSystem::oilPhaseIdx;
186 return FluidSystem::waterPhaseIdx;
189 inline int compIdx_()
const
192 return FluidSystem::oilCompIdx;
194 return FluidSystem::waterCompIdx;
198 inline void initQuantities()
201 if (!this->solution_set_from_restart_) {
207 initializeConnections();
208 calculateAquiferCondition();
209 calculateAquiferConstants();
211 pressure_previous_.resize(this->connections_.size(), Scalar{0});
212 pressure_current_.resize(this->connections_.size(), Scalar{0});
213 Qai_.resize(this->connections_.size(), Scalar{0});
217 updateCellPressure(std::vector<Eval>& pressure_water,
const int idx,
const IntensiveQuantities& intQuants)
219 const auto& fs = intQuants.fluidState();
220 pressure_water.at(idx) = fs.pressure(phaseIdx_());
224 updateCellPressure(std::vector<Scalar>& pressure_water,
const int idx,
const IntensiveQuantities& intQuants)
226 const auto& fs = intQuants.fluidState();
227 pressure_water.at(idx) = fs.pressure(phaseIdx_()).value();
230 virtual void endTimeStep() = 0;
232 const int aquiferID_{};
233 const std::vector<Aquancon::AquancCell> connections_;
234 const Simulator& ebos_simulator_;
237 std::vector<Scalar> faceArea_connected_;
238 std::vector<int> cellToConnectionIdx_;
241 std::vector<Scalar> cell_depth_;
242 std::vector<Scalar> pressure_previous_;
243 std::vector<Eval> pressure_current_;
244 std::vector<Eval> Qai_;
245 std::vector<Scalar> alphai_;
253 bool solution_set_from_restart_ {
false};
255 void initializeConnections()
257 this->cell_depth_.resize(this->size(), this->aquiferDepth());
258 this->alphai_.resize(this->size(), 1.0);
259 this->faceArea_connected_.resize(this->size(), Scalar{0});
262 FaceDir::DirEnum faceDirection;
264 bool has_active_connection_on_proc =
false;
267 Scalar denom_face_areas{0};
268 this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(0), -1);
269 const auto& gridView = this->ebos_simulator_.vanguard().gridView();
270 for (std::size_t idx = 0; idx < this->size(); ++idx) {
271 const auto global_index = this->connections_[idx].global_index;
272 const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
273 auto elemIt = gridView.template begin< 0>();
275 std::advance(elemIt, cell_index);
278 if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
281 has_active_connection_on_proc =
true;
283 this->cellToConnectionIdx_[cell_index] = idx;
284 this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
287 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
288 auto elemIt = gridView.template begin< 0>();
289 const auto& elemEndIt = gridView.template end< 0>();
290 for (; elemIt != elemEndIt; ++elemIt) {
291 const auto& elem = *elemIt;
292 unsigned cell_index = elemMapper.index(elem);
293 int idx = this->cellToConnectionIdx_[cell_index];
299 auto isIt = gridView.ibegin(elem);
300 const auto& isEndIt = gridView.iend(elem);
301 for (; isIt != isEndIt; ++ isIt) {
303 const auto& intersection = *isIt;
306 if (!intersection.boundary())
309 int insideFaceIdx = intersection.indexInInside();
310 switch (insideFaceIdx) {
312 faceDirection = FaceDir::XMinus;
315 faceDirection = FaceDir::XPlus;
318 faceDirection = FaceDir::YMinus;
321 faceDirection = FaceDir::YPlus;
324 faceDirection = FaceDir::ZMinus;
327 faceDirection = FaceDir::ZPlus;
330 OPM_THROW(std::logic_error,
331 "Internal error in initialization of aquifer.");
335 if (faceDirection == this->connections_[idx].face_dir) {
336 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
340 denom_face_areas += this->faceArea_connected_.at(idx);
343 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
344 comm.sum(&denom_face_areas, 1);
345 const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
346 for (std::size_t idx = 0; idx < this->size(); ++idx) {
348 this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
350 : this->faceArea_connected_.at(idx) / denom_face_areas;
353 if (this->solution_set_from_restart_) {
354 this->rescaleProducedVolume(has_active_connection_on_proc);
358 void rescaleProducedVolume(
const bool has_active_connection_on_proc)
366 const auto this_area = has_active_connection_on_proc
367 ? std::accumulate(this->alphai_.begin(),
372 const auto tot_area = this->ebos_simulator_.vanguard()
373 .grid().comm().sum(this_area);
375 this->W_flux_ *= this_area / tot_area;
378 virtual void assignRestartData(
const data::AquiferData& xaq) = 0;
380 virtual void calculateInflowRate(
int idx,
const Simulator& simulator) = 0;
382 virtual void calculateAquiferCondition() = 0;
384 virtual void calculateAquiferConstants() = 0;
386 virtual Scalar aquiferDepth()
const = 0;
389 virtual Scalar calculateReservoirEquilibrium()
392 std::vector<Scalar> pw_aquifer;
393 Scalar water_pressure_reservoir;
395 ElementContext elemCtx(this->ebos_simulator_);
396 const auto& gridView = this->ebos_simulator_.gridView();
397 auto elemIt = gridView.template begin<0>();
398 const auto& elemEndIt = gridView.template end<0>();
399 for (; elemIt != elemEndIt; ++elemIt) {
400 const auto& elem = *elemIt;
401 elemCtx.updatePrimaryStencil(elem);
403 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
404 const auto idx = this->cellToConnectionIdx_[cellIdx];
408 elemCtx.updatePrimaryIntensiveQuantities(0);
409 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
410 const auto& fs = iq0.fluidState();
412 water_pressure_reservoir = fs.pressure(phaseIdx_()).value();
413 const auto water_density = fs.density(phaseIdx_());
416 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
418 pw_aquifer.push_back(this->alphai_[idx] *
419 (water_pressure_reservoir - water_density.value()*gdz));
423 const auto& comm = ebos_simulator_.vanguard().grid().comm();
426 vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
427 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
431 return vals[1] / vals[0];
Definition: AquiferInterface.hpp:49
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27