21#ifndef OPM_AQUIFERCT_HEADER_INCLUDED
22#define OPM_AQUIFERCT_HEADER_INCLUDED
24#include <opm/simulators/aquifers/AquiferAnalytical.hpp>
26#include <opm/input/eclipse/EclipseState/Aquifer/AquiferCT.hpp>
28#include <opm/output/data/Aquifer.hpp>
38template <
typename TypeTag>
44 using typename Base::BlackoilIndices;
45 using typename Base::ElementContext;
46 using typename Base::Eval;
47 using typename Base::FluidState;
48 using typename Base::FluidSystem;
49 using typename Base::IntensiveQuantities;
50 using typename Base::RateVector;
51 using typename Base::Scalar;
52 using typename Base::Simulator;
53 using typename Base::ElementMapper;
56 const Simulator& ebosSimulator,
57 const AquiferCT::AQUCT_data& aquct_data)
58 :
Base(aquct_data.aquiferID, connections, ebosSimulator)
59 , aquct_data_(aquct_data)
66 result.pressure_previous_ = {1.0, 2.0, 3.0};
67 result.pressure_current_ = {4.0, 5.0};
68 result.Qai_ = {{6.0}};
71 result.fluxValue_ = 9.0;
72 result.dimensionless_time_ = 10.0;
73 result.dimensionless_pressure_ = 11.0;
78 void endTimeStep()
override
80 for (
const auto& q : this->Qai_) {
81 this->W_flux_ += q * this->ebos_simulator_.timeStepSize();
83 this->fluxValue_ = this->W_flux_.value();
84 const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
85 comm.sum(&this->fluxValue_, 1);
88 data::AquiferData aquiferData()
const override
90 data::AquiferData data;
91 data.aquiferID = this->aquiferID();
93 data.pressure = this->pa0_;
95 for (
const auto& q : this->Qai_) {
96 data.fluxRate += q.value();
98 data.volume = this->W_flux_.value();
99 data.initPressure = this->pa0_;
101 auto* aquCT = data.typeData.template create<data::AquiferType::CarterTracy>();
103 aquCT->dimensionless_time = this->dimensionless_time_;
104 aquCT->dimensionless_pressure = this->dimensionless_pressure_;
105 aquCT->influxConstant = this->aquct_data_.influxConstant();
107 if (!this->co2store_()) {
108 aquCT->timeConstant = this->aquct_data_.timeConstant();
109 aquCT->waterDensity = this->aquct_data_.waterDensity();
110 aquCT->waterViscosity = this->aquct_data_.waterViscosity();
112 aquCT->waterDensity = this->rhow_;
113 aquCT->timeConstant = this->Tc_;
114 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
115 aquCT->waterViscosity = this->Tc_ * this->aquct_data_.permeability / x;
121 template<
class Serializer>
122 void serializeOp(Serializer& serializer)
124 serializer(
static_cast<Base&
>(*
this));
125 serializer(fluxValue_);
126 serializer(dimensionless_time_);
127 serializer(dimensionless_pressure_);
133 this->fluxValue_ == rhs.fluxValue_ &&
134 this->dimensionless_time_ == rhs.dimensionless_time_ &&
135 this->dimensionless_pressure_ == rhs.dimensionless_pressure_;
140 AquiferCT::AQUCT_data aquct_data_;
144 Scalar fluxValue_{0};
146 Scalar dimensionless_time_{0};
147 Scalar dimensionless_pressure_{0};
149 void assignRestartData(
const data::AquiferData& xaq)
override
151 this->fluxValue_ = xaq.volume;
152 this->rhow_ = this->aquct_data_.waterDensity();
155 std::pair<Scalar, Scalar>
156 getInfluenceTableValues(
const Scalar td_plus_dt)
159 this->dimensionless_pressure_ =
160 linearInterpolation(this->aquct_data_.dimensionless_time,
161 this->aquct_data_.dimensionless_pressure,
162 this->dimensionless_time_);
165 linearInterpolation(this->aquct_data_.dimensionless_time,
166 this->aquct_data_.dimensionless_pressure,
169 const auto PItdprime =
170 linearInterpolationDerivative(this->aquct_data_.dimensionless_time,
171 this->aquct_data_.dimensionless_pressure,
174 return std::make_pair(PItd, PItdprime);
177 Scalar dpai(
const int idx)
const
180 this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth());
182 const auto dp = this->pa0_ + this->rhow_*gdz
183 - this->pressure_previous_.at(idx);
189 std::pair<Scalar, Scalar>
190 calculateEqnConstants(
const int idx,
const Simulator& simulator)
192 const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / this->Tc_;
193 this->dimensionless_time_ = simulator.time() / this->Tc_;
195 const auto [PItd, PItdprime] = this->getInfluenceTableValues(td_plus_dt);
197 const auto denom = this->Tc_ * (PItd - this->dimensionless_time_*PItdprime);
198 const auto a = (this->beta_*dpai(idx) - this->fluxValue_*PItdprime) / denom;
199 const auto b = this->beta_ / denom;
201 return std::make_pair(a, b);
204 std::size_t pvtRegionIdx()
const
206 return this->aquct_data_.pvttableID - 1;
210 inline void calculateInflowRate(
int idx,
const Simulator& simulator)
override
212 const auto [a, b] = this->calculateEqnConstants(idx, simulator);
214 this->Qai_.at(idx) = this->alphai_.at(idx) *
215 (a - b*(this->pressure_current_.at(idx) - this->pressure_previous_.at(idx)));
218 inline void calculateAquiferConstants()
override
220 if(this->co2store_()) {
221 const auto press = this->aquct_data_.initial_pressure.value();
222 Scalar temp = FluidSystem::reservoirTemperature();
223 if (this->aquct_data_.initial_temperature.has_value())
224 temp = this->aquct_data_.initial_temperature.value();
226 Scalar waterViscosity = 0.;
227 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
229 waterViscosity = FluidSystem::oilPvt().viscosity(pvtRegionIdx(), temp, press, rs);
230 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
233 waterViscosity = FluidSystem::waterPvt().viscosity(pvtRegionIdx(), temp, press, rsw, salt);
235 OPM_THROW(std::runtime_error,
"water or oil phase is needed to run CO2Store.");
237 const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
238 this->Tc_ = waterViscosity * x / this->aquct_data_.permeability;
240 this->Tc_ = this->aquct_data_.timeConstant();
242 this->beta_ = this->aquct_data_.influxConstant();
245 inline void calculateAquiferCondition()
override
247 if (this->solution_set_from_restart_) {
251 if (! this->aquct_data_.initial_pressure.has_value()) {
252 this->aquct_data_.initial_pressure =
253 this->calculateReservoirEquilibrium();
255 const auto& tables = this->ebos_simulator_.vanguard()
256 .eclState().getTableManager();
258 this->aquct_data_.finishInitialisation(tables);
261 this->pa0_ = this->aquct_data_.initial_pressure.value();
262 if (this->aquct_data_.initial_temperature.has_value())
263 this->Ta0_ = this->aquct_data_.initial_temperature.value();
265 if(this->co2store_()) {
266 const auto press = this->aquct_data_.initial_pressure.value();
268 Scalar temp = FluidSystem::reservoirTemperature();
269 if (this->aquct_data_.initial_temperature.has_value())
270 temp = this->aquct_data_.initial_temperature.value();
272 Scalar waterDensity = 0.;
273 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
275 waterDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx(), temp, press, rs)
276 * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIdx());
277 }
else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
278 Scalar salinity = 0.;
280 waterDensity = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx(), temp, press, rsw, salinity)
281 * FluidSystem::waterPvt().waterReferenceDensity(pvtRegionIdx());
283 OPM_THROW(std::runtime_error,
"water or oil phase is needed to run CO2Store.");
285 this->rhow_ = waterDensity;
287 this->rhow_ = this->aquct_data_.waterDensity();
291 virtual Scalar aquiferDepth()
const override
293 return this->aquct_data_.datum_depth;
Definition: AquiferAnalytical.hpp:55
Definition: AquiferCarterTracy.hpp:40
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27