57 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
58 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
59 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
60 using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
61 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
62 using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
63 using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
65 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
66 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
67 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
68 enum { enableEvaporation = getPropValue<TypeTag, Properties::EnableEvaporation>() };
69 enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
71 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
73 static constexpr int numEq = BlackoilIndices::numEq;
74 using Scalar = double;
76 using Eval = DenseAd::Evaluation<double, numEq>;
78 using FluidState = BlackOilFluidState<Eval,
82 BlackoilIndices::gasEnabled,
85 enableSaltPrecipitation,
87 BlackoilIndices::numPhases>;
91 const std::vector<Aquancon::AquancCell>& connections,
92 const Simulator& ebosSimulator)
96 this->initializeConnectionMappings();
103 void computeFaceAreaFraction(
const std::vector<double>&
total_face_area)
override
108 const auto eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
111 this->alphai_.assign(this->size(), Scalar{0});
114 std::transform(this->faceArea_connected_.begin(),
115 this->faceArea_connected_.end(),
116 this->alphai_.begin(),
123 this->area_fraction_ = this->totalFaceArea() /
tfa;
126 double totalFaceArea()
const override
128 return this->total_face_area_;
131 void initFromRestart(
const data::Aquifers&
aquiferSoln)
override
137 this->assignRestartData(
xaqPos->second);
139 this->W_flux_ =
xaqPos->second.volume * this->area_fraction_;
140 this->pa0_ =
xaqPos->second.initPressure;
142 this->solution_set_from_restart_ =
true;
145 void initialSolutionApplied()
override
150 void beginTimeStep()
override
152 ElementContext
elemCtx(this->ebos_simulator_);
153 OPM_BEGIN_PARALLEL_TRY_CATCH();
155 for (
const auto&
elem :
elements(this->ebos_simulator_.gridView())) {
159 const int idx = cellToConnectionIdx_[
cellIdx];
163 elemCtx.updateIntensiveQuantities(0);
164 const auto&
iq =
elemCtx.intensiveQuantities(0, 0);
168 OPM_END_PARALLEL_TRY_CATCH(
"AquiferAnalytical::beginTimeStep() failed: ",
169 this->ebos_simulator_.vanguard().grid().comm());
172 void addToSource(RateVector& rates,
174 const unsigned timeIdx)
override
176 const auto& model = this->ebos_simulator_.model();
178 const int idx = this->cellToConnectionIdx_[
cellIdx];
185 this->updateCellPressure(this->pressure_current_,
idx,
intQuants);
186 this->calculateInflowRate(
idx, this->ebos_simulator_);
188 rates[BlackoilIndices::conti0EqIdx + compIdx_()]
189 += this->Qai_[
idx] / model.dofTotalVolume(
cellIdx);
191 if constexpr (enableEnergy) {
193 if (this->Ta0_.has_value() && this->Qai_[
idx] > 0)
195 fs.setTemperature(this->Ta0_.value());
196 typedef typename std::decay<
decltype(
fs)>::type::Scalar
FsScalar;
198 const unsigned pvtRegionIdx =
intQuants.pvtRegionIndex();
200 paramCache.setMaxOilSat(this->ebos_simulator_.problem().maxOilSaturation(
cellIdx));
202 const auto&
h = FluidSystem::enthalpy(
fs,
paramCache, this->phaseIdx_());
203 fs.setEnthalpy(this->phaseIdx_(),
h);
205 rates[BlackoilIndices::contiEnergyEqIdx]
206 += this->Qai_[
idx] *
fs.enthalpy(this->phaseIdx_()) * FluidSystem::referenceDensity( this->phaseIdx_(),
intQuants.pvtRegionIndex()) / model.dofTotalVolume(
cellIdx);
211 std::size_t size()
const
213 return this->connections_.size();
216 template<
class Serializer>
228 return this->pressure_previous_ == rhs.pressure_previous_ &&
229 this->pressure_current_ == rhs.pressure_current_ &&
230 this->Qai_ == rhs.Qai_ &&
231 this->rhow_ == rhs.rhow_ &&
232 this->W_flux_ == rhs.W_flux_;
236 virtual void assignRestartData(
const data::AquiferData&
xaq) = 0;
237 virtual void calculateInflowRate(
int idx,
const Simulator&
simulator) = 0;
238 virtual void calculateAquiferCondition() = 0;
239 virtual void calculateAquiferConstants() = 0;
240 virtual Scalar aquiferDepth()
const = 0;
242 Scalar gravity_()
const
244 return this->ebos_simulator_.problem().gravity()[2];
249 if (this->co2store_or_h2store_())
250 return FluidSystem::oilCompIdx;
252 return FluidSystem::waterCompIdx;
255 void initQuantities()
258 if (! this->solution_set_from_restart_) {
262 this->initializeConnectionDepths();
263 this->calculateAquiferCondition();
264 this->calculateAquiferConstants();
266 this->pressure_previous_.resize(this->size(), Scalar{0});
267 this->pressure_current_.resize(this->size(), Scalar{0});
268 this->Qai_.resize(this->size(), Scalar{0});
287 void initializeConnectionMappings()
289 this->alphai_.resize(this->size(), 1.0);
290 this->faceArea_connected_.resize(this->size(), Scalar{0});
293 this->total_face_area_ = Scalar{0};
294 this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(0), -1);
295 const auto&
gridView = this->ebos_simulator_.vanguard().gridView();
296 for (std::size_t
idx = 0;
idx < this->size(); ++
idx) {
298 const int cell_index = this->ebos_simulator_.vanguard().compressedIndex(
global_index);
299 if (cell_index < 0) {
304 std::advance(
elemIt, cell_index);
307 if (
elemIt->partitionType() != Dune::InteriorEntity) {
311 this->cellToConnectionIdx_[cell_index] =
idx;
318 const auto&
elemMapper = this->ebos_simulator_.model().dofMapper();
321 const int idx = this->cellToConnectionIdx_[cell_index];
355 "Internal error in initialization of aquifer.");
359 this->faceArea_connected_[
idx] = this->connections_[
idx].influx_coeff;
364 this->total_face_area_ += this->faceArea_connected_.at(
idx);
368 void initializeConnectionDepths()
370 this->cell_depth_.resize(this->size(), this->aquiferDepth());
372 const auto&
gridView = this->ebos_simulator_.vanguard().gridView();
373 for (std::size_t
idx = 0;
idx < this->size(); ++
idx) {
374 const int cell_index = this->ebos_simulator_.vanguard()
376 if (cell_index < 0) {
381 std::advance(
elemIt, cell_index);
384 if (
elemIt->partitionType() != Dune::InteriorEntity) {
388 this->cell_depth_.at(
idx) = this->ebos_simulator_.vanguard().cellCenterDepth(cell_index);
393 Scalar calculateReservoirEquilibrium()
399 ElementContext
elemCtx(this->ebos_simulator_);
400 const auto&
gridView = this->ebos_simulator_.gridView();
405 const auto idx = this->cellToConnectionIdx_[
cellIdx];
409 elemCtx.updatePrimaryIntensiveQuantities(0);
410 const auto&
iq0 =
elemCtx.intensiveQuantities(0, 0);
411 const auto&
fs =
iq0.fluidState();
417 this->gravity_() * (this->cell_depth_[
idx] - this->aquiferDepth());
424 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
427 vals[0] = std::accumulate(this->alphai_.begin(),
this->alphai_.end(), Scalar{0});
435 const std::vector<Aquancon::AquancCell> connections_;
438 std::vector<Scalar> faceArea_connected_;
439 std::vector<int> cellToConnectionIdx_;
442 std::vector<Scalar> cell_depth_;
443 std::vector<Scalar> pressure_previous_;
444 std::vector<Eval> pressure_current_;
445 std::vector<Eval> Qai_;
446 std::vector<Scalar> alphai_;
450 std::optional<Scalar> Ta0_{};
453 Scalar total_face_area_{};
454 Scalar area_fraction_{Scalar{1}};
458 bool solution_set_from_restart_ {
false};
459 bool has_active_connection_on_proc_{
false};