56 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
57 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
58 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
59 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
60 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
61 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
62 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
64 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
65 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
66 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
67 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
68 enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
70 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
72 static constexpr int numEq = BlackoilIndices::numEq;
73 using Scalar = double;
75 using Eval = DenseAd::Evaluation<double, numEq>;
77 using FluidState = BlackOilFluidState<Eval,
81 BlackoilIndices::gasEnabled,
84 enableSaltPrecipitation,
86 BlackoilIndices::numPhases>;
90 const std::vector<Aquancon::AquancCell>& connections,
91 const Simulator& simulator)
93 , connections_(connections)
95 this->initializeConnectionMappings();
102 void computeFaceAreaFraction(
const std::vector<double>& total_face_area)
override
104 assert (total_face_area.size() >=
static_cast<std::vector<double>::size_type
>(this->aquiferID()));
106 const auto tfa = total_face_area[this->aquiferID() - 1];
107 const auto eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
109 if (tfa < eps_sqrt) {
110 this->alphai_.assign(this->size(), Scalar{0});
113 std::transform(this->faceArea_connected_.begin(),
114 this->faceArea_connected_.end(),
115 this->alphai_.begin(),
116 [tfa](
const Scalar area)
122 this->area_fraction_ = this->totalFaceArea() / tfa;
125 double totalFaceArea()
const override
127 return this->total_face_area_;
130 void initFromRestart(
const data::Aquifers& aquiferSoln)
override
132 auto xaqPos = aquiferSoln.find(this->aquiferID());
133 if (xaqPos == aquiferSoln.end())
136 this->assignRestartData(xaqPos->second);
138 this->W_flux_ = xaqPos->second.volume * this->area_fraction_;
139 this->pa0_ = xaqPos->second.initPressure;
141 this->solution_set_from_restart_ =
true;
144 void initialSolutionApplied()
override
149 void beginTimeStep()
override
151 ElementContext elemCtx(this->simulator_);
152 OPM_BEGIN_PARALLEL_TRY_CATCH();
154 for (
const auto& elem : elements(this->simulator_.gridView())) {
155 elemCtx.updatePrimaryStencil(elem);
157 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
158 const int idx = cellToConnectionIdx_[cellIdx];
162 elemCtx.updateIntensiveQuantities(0);
163 const auto& iq = elemCtx.intensiveQuantities(0, 0);
164 pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
167 OPM_END_PARALLEL_TRY_CATCH(
"AquiferAnalytical::beginTimeStep() failed: ",
168 this->simulator_.vanguard().grid().comm());
171 void addToSource(RateVector& rates,
172 const unsigned cellIdx,
173 const unsigned timeIdx)
override
175 const auto& model = this->simulator_.model();
177 const int idx = this->cellToConnectionIdx_[cellIdx];
181 const auto& intQuants = model.intensiveQuantities(cellIdx, timeIdx);
184 this->updateCellPressure(this->pressure_current_, idx, intQuants);
185 this->calculateInflowRate(idx, this->simulator_);
187 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
188 += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
190 if constexpr (enableEnergy) {
191 auto fs = intQuants.fluidState();
192 if (this->Ta0_.has_value() && this->Qai_[idx] > 0)
194 fs.setTemperature(this->Ta0_.value());
195 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
196 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
197 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
198 paramCache.setRegionIndex(pvtRegionIdx);
199 paramCache.setMaxOilSat(this->simulator_.problem().maxOilSaturation(cellIdx));
200 paramCache.updatePhase(fs, this->phaseIdx_());
201 const auto& h = FluidSystem::enthalpy(fs, paramCache, this->phaseIdx_());
202 fs.setEnthalpy(this->phaseIdx_(), h);
204 rates[BlackoilIndices::contiEnergyEqIdx]
205 += this->Qai_[idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
210 std::size_t size()
const
212 return this->connections_.size();
215 template<
class Serializer>
216 void serializeOp(Serializer& serializer)
218 serializer(pressure_previous_);
219 serializer(pressure_current_);
227 return this->pressure_previous_ == rhs.pressure_previous_ &&
228 this->pressure_current_ == rhs.pressure_current_ &&
229 this->Qai_ == rhs.Qai_ &&
230 this->rhow_ == rhs.rhow_ &&
231 this->W_flux_ == rhs.W_flux_;
235 virtual void assignRestartData(
const data::AquiferData& xaq) = 0;
236 virtual void calculateInflowRate(
int idx,
const Simulator& simulator) = 0;
237 virtual void calculateAquiferCondition() = 0;
238 virtual void calculateAquiferConstants() = 0;
239 virtual Scalar aquiferDepth()
const = 0;
241 Scalar gravity_()
const
243 return this->simulator_.problem().gravity()[2];
248 if (this->co2store_or_h2store_())
249 return FluidSystem::oilCompIdx;
251 return FluidSystem::waterCompIdx;
254 void initQuantities()
257 if (! this->solution_set_from_restart_) {
261 this->initializeConnectionDepths();
262 this->calculateAquiferCondition();
263 this->calculateAquiferConstants();
265 this->pressure_previous_.resize(this->size(), Scalar{0});
266 this->pressure_current_.resize(this->size(), Scalar{0});
267 this->Qai_.resize(this->size(), Scalar{0});
270 void updateCellPressure(std::vector<Eval>& pressure_water,
272 const IntensiveQuantities& intQuants)
274 const auto& fs = intQuants.fluidState();
275 pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
278 void updateCellPressure(std::vector<Scalar>& pressure_water,
280 const IntensiveQuantities& intQuants)
282 const auto& fs = intQuants.fluidState();
283 pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
286 void initializeConnectionMappings()
288 this->alphai_.resize(this->size(), 1.0);
289 this->faceArea_connected_.resize(this->size(), Scalar{0});
292 this->total_face_area_ = Scalar{0};
293 this->cellToConnectionIdx_.resize(this->simulator_.gridView().size(0), -1);
294 const auto& gridView = this->simulator_.vanguard().gridView();
295 for (std::size_t idx = 0; idx < this->size(); ++idx) {
296 const auto global_index = this->connections_[idx].global_index;
297 const int cell_index = this->simulator_.vanguard().compressedIndex(global_index);
298 if (cell_index < 0) {
302 auto elemIt = gridView.template begin< 0>();
303 std::advance(elemIt, cell_index);
306 if (elemIt->partitionType() != Dune::InteriorEntity) {
310 this->cellToConnectionIdx_[cell_index] = idx;
314 FaceDir::DirEnum faceDirection;
317 const auto& elemMapper = this->simulator_.model().dofMapper();
318 for (
const auto& elem : elements(gridView)) {
319 const unsigned cell_index = elemMapper.index(elem);
320 const int idx = this->cellToConnectionIdx_[cell_index];
327 for (
const auto& intersection : intersections(gridView, elem)) {
329 if (! intersection.boundary()) {
333 switch (intersection.indexInInside()) {
335 faceDirection = FaceDir::XMinus;
338 faceDirection = FaceDir::XPlus;
341 faceDirection = FaceDir::YMinus;
344 faceDirection = FaceDir::YPlus;
347 faceDirection = FaceDir::ZMinus;
350 faceDirection = FaceDir::ZPlus;
353 OPM_THROW(std::logic_error,
354 "Internal error in initialization of aquifer.");
357 if (faceDirection == this->connections_[idx].face_dir) {
358 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
363 this->total_face_area_ += this->faceArea_connected_.at(idx);
367 void initializeConnectionDepths()
369 this->cell_depth_.resize(this->size(), this->aquiferDepth());
371 const auto& gridView = this->simulator_.vanguard().gridView();
372 for (std::size_t idx = 0; idx < this->size(); ++idx) {
373 const int cell_index = this->simulator_.vanguard()
374 .compressedIndex(this->connections_[idx].global_index);
375 if (cell_index < 0) {
379 auto elemIt = gridView.template begin< 0>();
380 std::advance(elemIt, cell_index);
383 if (elemIt->partitionType() != Dune::InteriorEntity) {
387 this->cell_depth_.at(idx) = this->simulator_.vanguard().cellCenterDepth(cell_index);
392 Scalar calculateReservoirEquilibrium()
395 std::vector<Scalar> pw_aquifer;
396 Scalar water_pressure_reservoir;
398 ElementContext elemCtx(this->simulator_);
399 const auto& gridView = this->simulator_.gridView();
400 for (
const auto& elem : elements(gridView)) {
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(this->phaseIdx_()).value();
413 const auto water_density = fs.density(this->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 = this->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];
434 const std::vector<Aquancon::AquancCell> connections_;
437 std::vector<Scalar> faceArea_connected_;
438 std::vector<int> cellToConnectionIdx_;
441 std::vector<Scalar> cell_depth_;
442 std::vector<Scalar> pressure_previous_;
443 std::vector<Eval> pressure_current_;
444 std::vector<Eval> Qai_;
445 std::vector<Scalar> alphai_;
449 std::optional<Scalar> Ta0_{};
452 Scalar total_face_area_{};
453 Scalar area_fraction_{Scalar{1}};
457 bool solution_set_from_restart_ {
false};
458 bool has_active_connection_on_proc_{
false};