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();
100 void computeFaceAreaFraction(
const std::vector<Scalar>&
total_face_area)
override
105 const auto eps_sqrt = std::sqrt(std::numeric_limits<Scalar>::epsilon());
108 this->alphai_.assign(this->size(), Scalar{0});
111 std::transform(this->faceArea_connected_.begin(),
112 this->faceArea_connected_.end(),
113 this->alphai_.begin(),
114 [
tfa](
const Scalar area)
120 this->area_fraction_ = this->totalFaceArea() /
tfa;
123 Scalar totalFaceArea()
const override
125 return this->total_face_area_;
128 void initFromRestart(
const data::Aquifers&
aquiferSoln)
override
134 this->assignRestartData(
xaqPos->second);
136 this->W_flux_ =
xaqPos->second.volume * this->area_fraction_;
137 this->pa0_ =
xaqPos->second.initPressure;
139 this->solution_set_from_restart_ =
true;
142 void initialSolutionApplied()
override
147 void beginTimeStep()
override
149 ElementContext elemCtx(this->simulator_);
150 OPM_BEGIN_PARALLEL_TRY_CATCH();
152 for (
const auto&
elem :
elements(this->simulator_.gridView())) {
153 elemCtx.updatePrimaryStencil(
elem);
155 const int cellIdx = elemCtx.globalSpaceIndex(0, 0);
156 const int idx = cellToConnectionIdx_[
cellIdx];
160 elemCtx.updateIntensiveQuantities(0);
161 const auto&
iq = elemCtx.intensiveQuantities(0, 0);
165 OPM_END_PARALLEL_TRY_CATCH(
"AquiferAnalytical::beginTimeStep() failed: ",
166 this->simulator_.vanguard().grid().comm());
169 void addToSource(RateVector& rates,
171 const unsigned timeIdx)
override
173 const auto& model = this->simulator_.model();
175 const int idx = this->cellToConnectionIdx_[
cellIdx];
179 const auto& intQuants = model.intensiveQuantities(
cellIdx,
timeIdx);
182 this->updateCellPressure(this->pressure_current_,
idx, intQuants);
183 this->calculateInflowRate(
idx, this->simulator_);
185 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
186 += this->Qai_[
idx] / model.dofTotalVolume(
cellIdx);
188 if constexpr (enableEnergy) {
189 auto fs = intQuants.fluidState();
190 if (this->Ta0_.has_value() && this->Qai_[
idx] > 0)
192 fs.setTemperature(this->Ta0_.value());
193 typedef typename std::decay<
decltype(fs)>::type::Scalar
FsScalar;
195 const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
198 paramCache.updatePhase(fs, this->phaseIdx_());
199 const auto&
h = FluidSystem::enthalpy(fs,
paramCache, this->phaseIdx_());
200 fs.setEnthalpy(this->phaseIdx_(),
h);
202 rates[BlackoilIndices::contiEnergyEqIdx]
203 += this->Qai_[
idx] *fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(), intQuants.pvtRegionIndex()) / model.dofTotalVolume(
cellIdx);
208 std::size_t size()
const
210 return this->connections_.size();
213 template<
class Serializer>
225 return this->pressure_previous_ == rhs.pressure_previous_ &&
226 this->pressure_current_ == rhs.pressure_current_ &&
227 this->Qai_ == rhs.Qai_ &&
228 this->rhow_ == rhs.rhow_ &&
229 this->W_flux_ == rhs.W_flux_;
233 virtual void assignRestartData(
const data::AquiferData&
xaq) = 0;
234 virtual void calculateInflowRate(
int idx,
const Simulator& simulator) = 0;
235 virtual void calculateAquiferCondition() = 0;
236 virtual void calculateAquiferConstants() = 0;
237 virtual Scalar aquiferDepth()
const = 0;
239 Scalar gravity_()
const
241 return this->simulator_.problem().gravity()[2];
246 if (this->co2store_or_h2store_())
247 return FluidSystem::oilCompIdx;
249 return FluidSystem::waterCompIdx;
252 void initQuantities()
255 if (! this->solution_set_from_restart_) {
259 this->initializeConnectionDepths();
260 this->calculateAquiferCondition();
261 this->calculateAquiferConstants();
263 this->pressure_previous_.resize(this->size(), Scalar{0});
264 this->pressure_current_.resize(this->size(), Scalar{0});
265 this->Qai_.resize(this->size(), Scalar{0});
270 const IntensiveQuantities& intQuants)
272 const auto& fs = intQuants.fluidState();
278 const IntensiveQuantities& intQuants)
280 const auto& fs = intQuants.fluidState();
284 void initializeConnectionMappings()
286 this->alphai_.resize(this->size(), 1.0);
287 this->faceArea_connected_.resize(this->size(), Scalar{0});
290 this->total_face_area_ = Scalar{0};
291 this->cellToConnectionIdx_.resize(this->simulator_.gridView().size(0), -1);
292 const auto& gridView = this->simulator_.vanguard().gridView();
293 for (std::size_t
idx = 0;
idx < this->size(); ++
idx) {
295 const int cell_index = this->simulator_.vanguard().compressedIndex(
global_index);
296 if (cell_index < 0) {
300 auto elemIt = gridView.template begin< 0>();
301 std::advance(
elemIt, cell_index);
304 if (
elemIt->partitionType() != Dune::InteriorEntity) {
308 this->cellToConnectionIdx_[cell_index] =
idx;
315 const auto&
elemMapper = this->simulator_.model().dofMapper();
318 const int idx = this->cellToConnectionIdx_[cell_index];
327 if (! intersection.boundary()) {
331 switch (intersection.indexInInside()) {
352 "Internal error in initialization of aquifer.");
356 this->faceArea_connected_[
idx] = this->connections_[
idx].influx_coeff;
361 this->total_face_area_ += this->faceArea_connected_.at(
idx);
365 void initializeConnectionDepths()
367 this->cell_depth_.resize(this->size(), this->aquiferDepth());
369 const auto& gridView = this->simulator_.vanguard().gridView();
370 for (std::size_t
idx = 0;
idx < this->size(); ++
idx) {
371 const int cell_index = this->simulator_.vanguard()
373 if (cell_index < 0) {
377 auto elemIt = gridView.template begin< 0>();
378 std::advance(
elemIt, cell_index);
381 if (
elemIt->partitionType() != Dune::InteriorEntity) {
385 this->cell_depth_.at(
idx) = this->simulator_.vanguard().cellCenterDepth(cell_index);
390 Scalar calculateReservoirEquilibrium()
396 ElementContext elemCtx(this->simulator_);
397 const auto& gridView = this->simulator_.gridView();
399 elemCtx.updatePrimaryStencil(
elem);
401 const auto cellIdx = elemCtx.globalSpaceIndex(0, 0);
402 const auto idx = this->cellToConnectionIdx_[
cellIdx];
406 elemCtx.updatePrimaryIntensiveQuantities(0);
407 const auto&
iq0 = elemCtx.intensiveQuantities(0, 0);
408 const auto& fs =
iq0.fluidState();
414 this->gravity_() * (this->cell_depth_[
idx] - this->aquiferDepth());
421 const auto& comm = this->simulator_.vanguard().grid().comm();
424 vals[0] = std::accumulate(this->alphai_.begin(),
this->alphai_.end(), Scalar{0});
432 const std::vector<Aquancon::AquancCell> connections_;
435 std::vector<Scalar> faceArea_connected_;
436 std::vector<int> cellToConnectionIdx_;
439 std::vector<Scalar> cell_depth_;
440 std::vector<Scalar> pressure_previous_;
441 std::vector<Eval> pressure_current_;
442 std::vector<Eval> Qai_;
443 std::vector<Scalar> alphai_;
447 std::optional<Scalar> Ta0_{};
450 Scalar total_face_area_{};
451 Scalar area_fraction_{Scalar{1}};
455 bool solution_set_from_restart_ {
false};
456 bool has_active_connection_on_proc_{
false};