75 static constexpr int numEq = BlackoilIndices::numEq;
77 using Eval = DenseAd::Evaluation<Scalar, numEq>;
83 BlackoilIndices::gasEnabled,
86 enableSaltPrecipitation,
88 BlackoilIndices::numPhases>;
92 const std::vector<Aquancon::AquancCell>&
connections,
93 const Simulator& simulator)
97 this->initializeConnectionMappings();
104 void computeFaceAreaFraction(
const std::vector<Scalar>&
total_face_area)
override
109 const auto eps_sqrt = std::sqrt(std::numeric_limits<Scalar>::epsilon());
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
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())) {
160 const int idx = cellToConnectionIdx_[
cellIdx];
164 elemCtx.updateIntensiveQuantities(0);
165 const auto&
iq =
elemCtx.intensiveQuantities(0, 0);
169 OPM_END_PARALLEL_TRY_CATCH(
"AquiferAnalytical::beginTimeStep() failed: ",
170 this->simulator_.vanguard().grid().comm());
173 void addToSource(RateVector& rates,
175 const unsigned timeIdx)
override
177 const auto& model = this->simulator_.model();
179 const int idx = this->cellToConnectionIdx_[
cellIdx];
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) {
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;
199 const unsigned pvtRegionIdx =
intQuants.pvtRegionIndex();
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>
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});
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) {
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;
319 const auto&
elemMapper = this->simulator_.model().dofMapper();
322 const int idx = this->cellToConnectionIdx_[cell_index];
331 if (! intersection.boundary()) {
335 switch (intersection.indexInInside()) {
356 "Internal error in initialization of aquifer.");
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()
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()
400 ElementContext
elemCtx(this->simulator_);
401 const auto& gridView = this->simulator_.gridView();
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();
418 this->gravity_() * (this->cell_depth_[
idx] - this->aquiferDepth());
425 const auto& comm = this->simulator_.vanguard().grid().comm();
428 vals[0] = std::accumulate(this->alphai_.begin(),
this->alphai_.end(), Scalar{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};