46 using typename Base::BlackoilIndices;
47 using typename Base::ElementContext;
48 using typename Base::Eval;
49 using typename Base::FluidState;
50 using typename Base::FluidSystem;
51 using typename Base::IntensiveQuantities;
52 using typename Base::RateVector;
53 using typename Base::Scalar;
54 using typename Base::Simulator;
55 using typename Base::ElementMapper;
58 const Simulator& simulator,
59 const AquiferCT::AQUCT_data& aquct_data)
60 :
Base(aquct_data.aquiferID, connections, simulator)
61 , aquct_data_(aquct_data)
68 result.pressure_previous_ = {1.0, 2.0, 3.0};
69 result.pressure_current_ = {4.0, 5.0};
70 result.Qai_ = {{6.0}};
73 result.fluxValue_ = 9.0;
74 result.dimensionless_time_ = 10.0;
75 result.dimensionless_pressure_ = 11.0;
80 void endTimeStep()
override
82 for (
const auto& q : this->Qai_) {
83 this->W_flux_ += q * this->simulator_.timeStepSize();
85 this->fluxValue_ = this->W_flux_.value();
86 const auto& comm = this->simulator_.vanguard().grid().comm();
87 comm.sum(&this->fluxValue_, 1);
90 data::AquiferData aquiferData()
const override
92 data::AquiferData data;
93 data.aquiferID = this->aquiferID();
95 data.pressure = this->pa0_;
97 for (
const auto& q : this->Qai_) {
98 data.fluxRate += q.value();
100 data.volume = this->W_flux_.value();
101 data.initPressure = this->pa0_;
103 auto* aquCT = data.typeData.template create<data::AquiferType::CarterTracy>();
105 aquCT->dimensionless_time = this->dimensionless_time_;
106 aquCT->dimensionless_pressure = this->dimensionless_pressure_;
107 aquCT->influxConstant = this->aquct_data_.influxConstant();
109 if (!this->co2store_or_h2store_()) {
110 aquCT->timeConstant = this->aquct_data_.timeConstant();
111 aquCT->waterDensity = this->aquct_data_.waterDensity();
112 aquCT->waterViscosity = this->aquct_data_.waterViscosity();
114 aquCT->waterDensity = this->rhow_;
115 aquCT->timeConstant = this->Tc_;
116 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
117 aquCT->waterViscosity = this->Tc_ * this->aquct_data_.permeability / x;
123 template<
class Serializer>
124 void serializeOp(Serializer& serializer)
126 serializer(
static_cast<Base&
>(*
this));
127 serializer(fluxValue_);
128 serializer(dimensionless_time_);
129 serializer(dimensionless_pressure_);
135 this->fluxValue_ == rhs.fluxValue_ &&
136 this->dimensionless_time_ == rhs.dimensionless_time_ &&
137 this->dimensionless_pressure_ == rhs.dimensionless_pressure_;
142 AquiferCT::AQUCT_data aquct_data_;
146 Scalar fluxValue_{0};
148 Scalar dimensionless_time_{0};
149 Scalar dimensionless_pressure_{0};
151 void assignRestartData(
const data::AquiferData& xaq)
override
153 this->fluxValue_ = xaq.volume;
154 this->rhow_ = this->aquct_data_.waterDensity();
157 std::pair<Scalar, Scalar>
158 getInfluenceTableValues(
const Scalar td_plus_dt)
161 this->dimensionless_pressure_ =
162 linearInterpolation(this->aquct_data_.dimensionless_time,
163 this->aquct_data_.dimensionless_pressure,
164 this->dimensionless_time_);
167 linearInterpolation(this->aquct_data_.dimensionless_time,
168 this->aquct_data_.dimensionless_pressure,
171 const auto PItdprime =
172 linearInterpolationDerivative(this->aquct_data_.dimensionless_time,
173 this->aquct_data_.dimensionless_pressure,
176 return std::make_pair(PItd, PItdprime);
179 Scalar dpai(
const int idx)
const
182 this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth());
184 const auto dp = this->pa0_ + this->rhow_*gdz
185 - this->pressure_previous_.at(idx);
191 std::pair<Scalar, Scalar>
192 calculateEqnConstants(
const int idx,
const Simulator& simulator)
194 const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / this->Tc_;
195 this->dimensionless_time_ = simulator.time() / this->Tc_;
197 const auto [PItd, PItdprime] = this->getInfluenceTableValues(td_plus_dt);
199 const auto denom = this->Tc_ * (PItd - this->dimensionless_time_*PItdprime);
200 const auto a = (this->beta_*dpai(idx) - this->fluxValue_*PItdprime) / denom;
201 const auto b = this->beta_ / denom;
203 return std::make_pair(a, b);
206 std::size_t pvtRegionIdx()
const
208 return this->aquct_data_.pvttableID - 1;
212 void calculateInflowRate(
int idx,
const Simulator& simulator)
override
214 const auto [a, b] = this->calculateEqnConstants(idx, simulator);
216 this->Qai_.at(idx) = this->alphai_.at(idx) *
217 (a - b*(this->pressure_current_.at(idx) - this->pressure_previous_.at(idx)));
220 void calculateAquiferConstants()
override
222 this->Tc_ = this->co2store_or_h2store_()
223 ? this->timeConstantCO2Store()
224 : this->aquct_data_.timeConstant();
226 this->beta_ = this->aquct_data_.influxConstant();
229 void calculateAquiferCondition()
override
231 if (this->solution_set_from_restart_) {
235 if (! this->aquct_data_.initial_pressure.has_value()) {
236 this->aquct_data_.initial_pressure =
237 this->calculateReservoirEquilibrium();
239 const auto& tables = this->simulator_.vanguard()
240 .eclState().getTableManager();
242 this->aquct_data_.finishInitialisation(tables);
245 this->pa0_ = this->aquct_data_.initial_pressure.value();
246 if (this->aquct_data_.initial_temperature.has_value()) {
247 this->Ta0_ = this->aquct_data_.initial_temperature.value();
250 this->rhow_ = this->co2store_or_h2store_()
251 ? this->waterDensityCO2Store()
252 : this->aquct_data_.waterDensity();
255 Scalar aquiferDepth()
const override
257 return this->aquct_data_.datum_depth;
261 Scalar timeConstantCO2Store()
const
263 const Scalar press = this->aquct_data_.initial_pressure.value();
264 const auto temp = this->reservoirTemperatureCO2Store();
266 auto waterViscosity = Scalar { 0 };
267 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
268 const auto rs = Scalar { 0 };
269 waterViscosity = FluidSystem::oilPvt()
270 .viscosity(pvtRegionIdx(), temp, press, rs);
272 else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
273 const auto salt = Scalar { 0 };
274 const auto rsw = Scalar { 0 };
275 waterViscosity = FluidSystem::waterPvt()
276 .viscosity(pvtRegionIdx(), temp, press, rsw, salt);
279 OPM_THROW(std::runtime_error,
"water or oil phase is needed to run CO2Store.");
282 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr
283 * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
285 return waterViscosity * x / this->aquct_data_.permeability;
288 Scalar waterDensityCO2Store()
const
290 const Scalar press = this->aquct_data_.initial_pressure.value();
291 const auto temp = this->reservoirTemperatureCO2Store();
293 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
294 const auto& pvt = FluidSystem::oilPvt();
295 const auto reg = this->pvtRegionIdx();
297 const auto rs = Scalar { 0 };
298 return pvt.inverseFormationVolumeFactor(reg, temp, press, rs)
299 * pvt.oilReferenceDensity(reg);
301 else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
302 const auto& pvt = FluidSystem::waterPvt();
303 const auto reg = this->pvtRegionIdx();
305 const auto salinity = Scalar { 0 };
306 const auto rsw = Scalar { 0 };
308 return pvt.inverseFormationVolumeFactor(reg, temp, press, rsw, salinity)
309 * pvt.waterReferenceDensity(reg);
312 OPM_THROW(std::runtime_error,
"water or oil phase is needed to run CO2Store.");
316 Scalar reservoirTemperatureCO2Store()
const
318 return this->aquct_data_.initial_temperature.has_value()
319 ? this->aquct_data_.initial_temperature.value()
320 : FluidSystem::reservoirTemperature();