22 #ifndef OPM_AQUIFERINTERFACE_HEADER_INCLUDED
23 #define OPM_AQUIFERINTERFACE_HEADER_INCLUDED
25 #include <opm/common/utility/numeric/linearInterpolation.hpp>
26 #include <opm/parser/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>() };
63 static const int numEq = BlackoilIndices::numEq;
64 typedef double Scalar;
66 typedef DenseAd::Evaluation<double, numEq> Eval;
68 typedef BlackOilFluidState<Eval,
72 BlackoilIndices::gasEnabled,
74 BlackoilIndices::numPhases>
79 const std::vector<Aquancon::AquancCell>& connections,
80 const Simulator& ebosSimulator)
82 , connections_(connections)
83 , ebos_simulator_(ebosSimulator)
92 void initFromRestart(
const data::Aquifers& aquiferSoln)
94 auto xaqPos = aquiferSoln.find(this->aquiferID());
95 if (xaqPos == aquiferSoln.end())
98 this->assignRestartData(xaqPos->second);
100 this->W_flux_ = xaqPos->second.volume;
101 this->pa0_ = xaqPos->second.initPressure;
102 this->solution_set_from_restart_ =
true;
105 void initialSolutionApplied()
112 ElementContext elemCtx(ebos_simulator_);
113 auto elemIt = ebos_simulator_.gridView().template begin<0>();
114 const auto& elemEndIt = ebos_simulator_.gridView().template end<0>();
115 OPM_BEGIN_PARALLEL_TRY_CATCH();
117 for (; elemIt != elemEndIt; ++elemIt) {
118 const auto& elem = *elemIt;
120 elemCtx.updatePrimaryStencil(elem);
122 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
123 const int idx = cellToConnectionIdx_[cellIdx];
127 elemCtx.updateIntensiveQuantities(0);
128 const auto& iq = elemCtx.intensiveQuantities(0, 0);
129 pressure_previous_[idx] = getValue(iq.fluidState().pressure(phaseIdx_()));
132 OPM_END_PARALLEL_TRY_CATCH(
"AquiferInterface::beginTimeStep() failed: ", ebos_simulator_.vanguard().grid().comm());
135 template <
class Context>
136 void addToSource(RateVector& rates,
137 const Context& context,
138 const unsigned spaceIdx,
139 const unsigned timeIdx)
141 const unsigned cellIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
143 const int idx = this->cellToConnectionIdx_[cellIdx];
150 const auto& intQuants = context.intensiveQuantities(spaceIdx, timeIdx);
153 this->updateCellPressure(this->pressure_current_, idx, intQuants);
154 this->calculateInflowRate(idx, context.simulator());
156 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
157 += this->Qai_[idx] / context.dofVolume(spaceIdx, timeIdx);
160 std::size_t size()
const {
161 return this->connections_.size();
164 int aquiferID()
const {
return this->aquiferID_; }
167 inline Scalar gravity_()
const
169 return ebos_simulator_.problem().gravity()[2];
172 inline bool co2store_()
const
174 return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
177 inline bool phaseIdx_()
const
180 return FluidSystem::oilPhaseIdx;
182 return FluidSystem::waterPhaseIdx;
185 inline bool compIdx_()
const
188 return FluidSystem::oilCompIdx;
190 return FluidSystem::waterCompIdx;
194 inline void initQuantities()
197 if (!this->solution_set_from_restart_) {
203 initializeConnections();
204 calculateAquiferCondition();
205 calculateAquiferConstants();
207 pressure_previous_.resize(this->connections_.size(), Scalar{0});
208 pressure_current_.resize(this->connections_.size(), Scalar{0});
209 Qai_.resize(this->connections_.size(), Scalar{0});
213 updateCellPressure(std::vector<Eval>& pressure_water,
const int idx,
const IntensiveQuantities& intQuants)
215 const auto& fs = intQuants.fluidState();
216 pressure_water.at(idx) = fs.pressure(phaseIdx_());
220 updateCellPressure(std::vector<Scalar>& pressure_water,
const int idx,
const IntensiveQuantities& intQuants)
222 const auto& fs = intQuants.fluidState();
223 pressure_water.at(idx) = fs.pressure(phaseIdx_()).value();
226 virtual void endTimeStep() = 0;
228 const int aquiferID_{};
229 const std::vector<Aquancon::AquancCell> connections_;
230 const Simulator& ebos_simulator_;
233 std::vector<Scalar> faceArea_connected_;
234 std::vector<int> cellToConnectionIdx_;
237 std::vector<Scalar> cell_depth_;
238 std::vector<Scalar> pressure_previous_;
239 std::vector<Eval> pressure_current_;
240 std::vector<Eval> Qai_;
241 std::vector<Scalar> alphai_;
249 bool solution_set_from_restart_ {
false};
251 void initializeConnections()
253 this->cell_depth_.resize(this->size(), this->aquiferDepth());
254 this->alphai_.resize(this->size(), 1.0);
255 this->faceArea_connected_.resize(this->size(), Scalar{0});
258 FaceDir::DirEnum faceDirection;
260 bool has_active_connection_on_proc =
false;
263 Scalar denom_face_areas{0};
264 this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(0), -1);
265 const auto& gridView = this->ebos_simulator_.vanguard().gridView();
266 for (std::size_t idx = 0; idx < this->size(); ++idx) {
267 const auto global_index = this->connections_[idx].global_index;
268 const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(global_index);
269 auto elemIt = gridView.template begin< 0>();
271 std::advance(elemIt, cell_index);
274 if ( cell_index < 0 || elemIt->partitionType() != Dune::InteriorEntity)
277 has_active_connection_on_proc =
true;
279 this->cellToConnectionIdx_[cell_index] = idx;
280 this->cell_depth_.at(idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
283 ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
284 auto elemIt = gridView.template begin< 0>();
285 const auto& elemEndIt = gridView.template end< 0>();
286 for (; elemIt != elemEndIt; ++elemIt) {
287 const auto& elem = *elemIt;
288 unsigned cell_index = elemMapper.index(elem);
289 int idx = this->cellToConnectionIdx_[cell_index];
295 auto isIt = gridView.ibegin(elem);
296 const auto& isEndIt = gridView.iend(elem);
297 for (; isIt != isEndIt; ++ isIt) {
299 const auto& intersection = *isIt;
302 if (!intersection.boundary())
305 int insideFaceIdx = intersection.indexInInside();
306 switch (insideFaceIdx) {
308 faceDirection = FaceDir::XMinus;
311 faceDirection = FaceDir::XPlus;
314 faceDirection = FaceDir::YMinus;
317 faceDirection = FaceDir::YPlus;
320 faceDirection = FaceDir::ZMinus;
323 faceDirection = FaceDir::ZPlus;
326 OPM_THROW(std::logic_error,
327 "Internal error in initialization of aquifer.");
331 if (faceDirection == this->connections_[idx].face_dir) {
332 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
336 denom_face_areas += this->faceArea_connected_.at(idx);
339 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
340 comm.sum(&denom_face_areas, 1);
341 const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
342 for (std::size_t idx = 0; idx < this->size(); ++idx) {
344 this->alphai_.at(idx) = (denom_face_areas < eps_sqrt)
346 : this->faceArea_connected_.at(idx) / denom_face_areas;
349 if (this->solution_set_from_restart_) {
350 this->rescaleProducedVolume(has_active_connection_on_proc);
354 void rescaleProducedVolume(
const bool has_active_connection_on_proc)
362 const auto this_area = has_active_connection_on_proc
363 ? std::accumulate(this->alphai_.begin(),
368 const auto tot_area = this->ebos_simulator_.vanguard()
369 .grid().comm().sum(this_area);
371 this->W_flux_ *= this_area / tot_area;
374 virtual void assignRestartData(
const data::AquiferData& xaq) = 0;
376 virtual void calculateInflowRate(
int idx,
const Simulator& simulator) = 0;
378 virtual void calculateAquiferCondition() = 0;
380 virtual void calculateAquiferConstants() = 0;
382 virtual Scalar aquiferDepth()
const = 0;
385 virtual Scalar calculateReservoirEquilibrium()
388 std::vector<Scalar> pw_aquifer;
389 Scalar water_pressure_reservoir;
391 ElementContext elemCtx(this->ebos_simulator_);
392 const auto& gridView = this->ebos_simulator_.gridView();
393 auto elemIt = gridView.template begin<0>();
394 const auto& elemEndIt = gridView.template end<0>();
395 for (; elemIt != elemEndIt; ++elemIt) {
396 const auto& elem = *elemIt;
397 elemCtx.updatePrimaryStencil(elem);
399 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
400 const auto idx = this->cellToConnectionIdx_[cellIdx];
404 elemCtx.updatePrimaryIntensiveQuantities(0);
405 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
406 const auto& fs = iq0.fluidState();
408 water_pressure_reservoir = fs.pressure(phaseIdx_()).value();
409 const auto water_density = fs.density(phaseIdx_());
412 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
414 pw_aquifer.push_back(this->alphai_[idx] *
415 (water_pressure_reservoir - water_density.value()*gdz));
419 const auto& comm = ebos_simulator_.vanguard().grid().comm();
422 vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
423 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
427 return vals[1] / vals[0];
Definition: AquiferInterface.hpp:49
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:26