75 static constexpr int numEq = BlackoilIndices::numEq;
77 using Eval = DenseAd::Evaluation<Scalar, numEq>;
79 using FluidState = BlackOilFluidState<Eval,
83 BlackoilIndices::gasEnabled,
86 enableSaltPrecipitation,
88 BlackoilIndices::numPhases>;
92 const std::vector<Aquancon::AquancCell>& connections,
93 const Simulator& simulator)
95 , connections_(connections)
97 this->initializeConnectionMappings();
104 void computeFaceAreaFraction(
const std::vector<Scalar>& total_face_area)
override
106 assert (total_face_area.size() >=
static_cast<typename std::vector<Scalar>::size_type
>(this->aquiferID()));
108 const auto tfa = total_face_area[this->aquiferID() - 1];
109 const auto eps_sqrt = std::sqrt(std::numeric_limits<Scalar>::epsilon());
111 if (tfa < eps_sqrt) {
112 this->alphai_.assign(this->size(), Scalar{0});
115 std::transform(this->faceArea_connected_.begin(),
116 this->faceArea_connected_.end(),
117 this->alphai_.begin(),
118 [tfa](
const Scalar area)
124 this->area_fraction_ = this->totalFaceArea() / tfa;
127 Scalar totalFaceArea()
const override
129 return this->total_face_area_;
132 void initFromRestart(
const data::Aquifers& aquiferSoln)
override
134 auto xaqPos = aquiferSoln.find(this->aquiferID());
135 if (xaqPos == aquiferSoln.end())
138 this->assignRestartData(xaqPos->second);
140 this->W_flux_ = xaqPos->second.volume * this->area_fraction_;
141 this->pa0_ = xaqPos->second.initPressure;
143 this->solution_set_from_restart_ =
true;
146 void initialSolutionApplied()
override
151 void beginTimeStep()
override
153 ElementContext elemCtx(this->simulator_);
154 OPM_BEGIN_PARALLEL_TRY_CATCH();
156 for (
const auto& elem : elements(this->simulator_.gridView())) {
157 elemCtx.updatePrimaryStencil(elem);
159 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
160 const int idx = cellToConnectionIdx_[cellIdx];
164 elemCtx.updateIntensiveQuantities(0);
165 const auto& iq = elemCtx.intensiveQuantities(0, 0);
166 pressure_previous_[idx] = getValue(iq.fluidState().pressure(this->phaseIdx_()));
169 OPM_END_PARALLEL_TRY_CATCH(
"AquiferAnalytical::beginTimeStep() failed: ",
170 this->simulator_.vanguard().grid().comm());
173 void addToSource(RateVector& rates,
174 const unsigned cellIdx,
175 const unsigned timeIdx)
override
177 const auto& model = this->simulator_.model();
179 const int idx = this->cellToConnectionIdx_[cellIdx];
183 const auto& intQuants = model.intensiveQuantities(cellIdx, timeIdx);
186 this->updateCellPressure(this->pressure_current_, idx, intQuants);
187 this->calculateInflowRate(idx, this->simulator_);
189 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
190 += this->Qai_[idx] / model.dofTotalVolume(cellIdx);
192 if constexpr (enableEnergy) {
193 auto fs = intQuants.fluidState();
194 if (this->Ta0_.has_value() && this->Qai_[idx] > 0)
196 fs.setTemperature(this->Ta0_.value());
197 typedef typename std::decay<
decltype(fs)>::type::Scalar FsScalar;
198 typename FluidSystem::template ParameterCache<FsScalar> paramCache;
199 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
200 paramCache.setRegionIndex(pvtRegionIdx);
201 paramCache.setMaxOilSat(this->simulator_.problem().maxOilSaturation(cellIdx));
202 paramCache.updatePhase(fs, this->phaseIdx_());
203 const auto& h = FluidSystem::enthalpy(fs, paramCache, this->phaseIdx_());
204 fs.setEnthalpy(this->phaseIdx_(), h);
206 rates[BlackoilIndices::contiEnergyEqIdx]
207 += this->Qai_[idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(cellIdx);
212 std::size_t size()
const
214 return this->connections_.size();
217 template<
class Serializer>
218 void serializeOp(Serializer& serializer)
220 serializer(pressure_previous_);
221 serializer(pressure_current_);
229 return this->pressure_previous_ == rhs.pressure_previous_ &&
230 this->pressure_current_ == rhs.pressure_current_ &&
231 this->Qai_ == rhs.Qai_ &&
232 this->rhow_ == rhs.rhow_ &&
233 this->W_flux_ == rhs.W_flux_;
237 virtual void assignRestartData(
const data::AquiferData& xaq) = 0;
238 virtual void calculateInflowRate(
int idx,
const Simulator& simulator) = 0;
239 virtual void calculateAquiferCondition() = 0;
240 virtual void calculateAquiferConstants() = 0;
241 virtual Scalar aquiferDepth()
const = 0;
243 Scalar gravity_()
const
245 return this->simulator_.problem().gravity()[2];
250 if (this->co2store_or_h2store_())
251 return FluidSystem::oilCompIdx;
253 return FluidSystem::waterCompIdx;
256 void initQuantities()
259 if (! this->solution_set_from_restart_) {
263 this->initializeConnectionDepths();
264 this->calculateAquiferCondition();
265 this->calculateAquiferConstants();
267 this->pressure_previous_.resize(this->size(), Scalar{0});
268 this->pressure_current_.resize(this->size(), Scalar{0});
269 this->Qai_.resize(this->size(), Scalar{0});
272 void updateCellPressure(std::vector<Eval>& pressure_water,
274 const IntensiveQuantities& intQuants)
276 const auto& fs = intQuants.fluidState();
277 pressure_water.at(idx) = fs.pressure(this->phaseIdx_());
280 void updateCellPressure(std::vector<Scalar>& pressure_water,
282 const IntensiveQuantities& intQuants)
284 const auto& fs = intQuants.fluidState();
285 pressure_water.at(idx) = fs.pressure(this->phaseIdx_()).value();
288 void initializeConnectionMappings()
290 this->alphai_.resize(this->size(), 1.0);
291 this->faceArea_connected_.resize(this->size(), Scalar{0});
294 this->total_face_area_ = Scalar{0};
295 this->cellToConnectionIdx_.resize(this->simulator_.gridView().size(0), -1);
296 const auto& gridView = this->simulator_.vanguard().gridView();
297 for (std::size_t idx = 0; idx < this->size(); ++idx) {
298 const auto global_index = this->connections_[idx].global_index;
299 const int cell_index = this->simulator_.vanguard().compressedIndex(global_index);
300 if (cell_index < 0) {
304 auto elemIt = gridView.template begin< 0>();
305 std::advance(elemIt, cell_index);
308 if (elemIt->partitionType() != Dune::InteriorEntity) {
312 this->cellToConnectionIdx_[cell_index] = idx;
316 FaceDir::DirEnum faceDirection;
319 const auto& elemMapper = this->simulator_.model().dofMapper();
320 for (
const auto& elem : elements(gridView)) {
321 const unsigned cell_index = elemMapper.index(elem);
322 const int idx = this->cellToConnectionIdx_[cell_index];
329 for (
const auto& intersection : intersections(gridView, elem)) {
331 if (! intersection.boundary()) {
335 switch (intersection.indexInInside()) {
337 faceDirection = FaceDir::XMinus;
340 faceDirection = FaceDir::XPlus;
343 faceDirection = FaceDir::YMinus;
346 faceDirection = FaceDir::YPlus;
349 faceDirection = FaceDir::ZMinus;
352 faceDirection = FaceDir::ZPlus;
355 OPM_THROW(std::logic_error,
356 "Internal error in initialization of aquifer.");
359 if (faceDirection == this->connections_[idx].face_dir) {
360 this->faceArea_connected_[idx] = this->connections_[idx].influx_coeff;
365 this->total_face_area_ += this->faceArea_connected_.at(idx);
369 void initializeConnectionDepths()
371 this->cell_depth_.resize(this->size(), this->aquiferDepth());
373 const auto& gridView = this->simulator_.vanguard().gridView();
374 for (std::size_t idx = 0; idx < this->size(); ++idx) {
375 const int cell_index = this->simulator_.vanguard()
376 .compressedIndex(this->connections_[idx].global_index);
377 if (cell_index < 0) {
381 auto elemIt = gridView.template begin< 0>();
382 std::advance(elemIt, cell_index);
385 if (elemIt->partitionType() != Dune::InteriorEntity) {
389 this->cell_depth_.at(idx) = this->simulator_.vanguard().cellCenterDepth(cell_index);
394 Scalar calculateReservoirEquilibrium()
397 std::vector<Scalar> pw_aquifer;
398 Scalar water_pressure_reservoir;
400 ElementContext elemCtx(this->simulator_);
401 const auto& gridView = this->simulator_.gridView();
402 for (
const auto& elem : elements(gridView)) {
403 elemCtx.updatePrimaryStencil(elem);
405 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
406 const auto idx = this->cellToConnectionIdx_[cellIdx];
410 elemCtx.updatePrimaryIntensiveQuantities(0);
411 const auto& iq0 = elemCtx.intensiveQuantities(0, 0);
412 const auto& fs = iq0.fluidState();
414 water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
415 const auto water_density = fs.density(this->phaseIdx_());
418 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
420 pw_aquifer.push_back(this->alphai_[idx] *
421 (water_pressure_reservoir - water_density.value()*gdz));
425 const auto& comm = this->simulator_.vanguard().grid().comm();
428 vals[0] = std::accumulate(this->alphai_.begin(), this->alphai_.end(), Scalar{0});
429 vals[1] = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), Scalar{0});
433 return vals[1] / vals[0];
436 const std::vector<Aquancon::AquancCell> connections_;
439 std::vector<Scalar> faceArea_connected_;
440 std::vector<int> cellToConnectionIdx_;
443 std::vector<Scalar> cell_depth_;
444 std::vector<Scalar> pressure_previous_;
445 std::vector<Eval> pressure_current_;
446 std::vector<Eval> Qai_;
447 std::vector<Scalar> alphai_;
451 std::optional<Scalar> Ta0_{};
454 Scalar total_face_area_{};
455 Scalar area_fraction_{Scalar{1}};
459 bool solution_set_from_restart_ {
false};
460 bool has_active_connection_on_proc_{
false};