diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index 8d6dc3773..e73b15503 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -84,6 +84,10 @@ public: data.type = data::AquiferType::CarterTracy; data.aquCT = std::make_shared(); + data.aquCT->timeConstant = this->aquct_data_.timeConstant(); + data.aquCT->influxConstant = this->aquct_data_.influxConstant(); + data.aquCT->waterDensity = this->aquct_data_.waterDensity(); + data.aquCT->waterViscosity = this->aquct_data_.waterViscosity(); data.aquCT->dimensionless_time = this->dimensionless_time_; data.aquCT->dimensionless_pressure = this->dimensionless_pressure_; @@ -92,10 +96,10 @@ public: protected: // Variables constants - const AquiferCT::AQUCT_data aquct_data_; + AquiferCT::AQUCT_data aquct_data_; + Scalar beta_; // Influx constant // TODO: it is possible it should be a AD variable - Scalar mu_w_{1}; // water viscosity Scalar fluxValue_{0}; // value of flux Scalar dimensionless_time_{0}; @@ -112,26 +116,31 @@ protected: { // We use the opm-common numeric linear interpolator this->dimensionless_pressure_ = - linearInterpolation(this->aquct_data_.td, - this->aquct_data_.pi, + linearInterpolation(this->aquct_data_.dimensionless_time, + this->aquct_data_.dimensionless_pressure, this->dimensionless_time_); const auto PItd = - linearInterpolation(this->aquct_data_.td, - this->aquct_data_.pi, td_plus_dt); + linearInterpolation(this->aquct_data_.dimensionless_time, + this->aquct_data_.dimensionless_pressure, + td_plus_dt); const auto PItdprime = - linearInterpolationDerivative(this->aquct_data_.td, - this->aquct_data_.pi, td_plus_dt); + linearInterpolationDerivative(this->aquct_data_.dimensionless_time, + this->aquct_data_.dimensionless_pressure, + td_plus_dt); return std::make_pair(PItd, PItdprime); } Scalar dpai(const int idx) const { - Scalar dp = this->pa0_ - + this->rhow_.at(idx).value() * this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth()) + const auto gdz = + this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth()); + + const auto dp = this->pa0_ + this->rhow_*gdz - this->pressure_previous_.at(idx); + return dp; } @@ -162,48 +171,29 @@ protected: inline void calculateAquiferConstants() override { - // We calculate the influx constant - beta_ = aquct_data_.c2 * aquct_data_.h * aquct_data_.theta * aquct_data_.phi_aq * aquct_data_.C_t - * aquct_data_.r_o * aquct_data_.r_o; - // We calculate the time constant - this->Tc_ = mu_w_ * aquct_data_.phi_aq * aquct_data_.C_t * aquct_data_.r_o * aquct_data_.r_o - / (aquct_data_.k_a * aquct_data_.c1); + this->Tc_ = this->aquct_data_.timeConstant(); + this->beta_ = this->aquct_data_.influxConstant(); } inline void calculateAquiferCondition() override { + if (! this->aquct_data_.initial_pressure.has_value()) { + this->aquct_data_.initial_pressure = + this->calculateReservoirEquilibrium(); - int pvttableIdx = aquct_data_.pvttableID - 1; - this->rhow_.resize(this->size(), 0.); - if (!aquct_data_.p0.first) { - this->pa0_ = this->calculateReservoirEquilibrium(); - } else { - this->pa0_ = aquct_data_.p0.second; + const auto& tables = this->ebos_simulator_.vanguard() + .eclState().getTableManager(); + + this->aquct_data_.finishInitialisation(tables); } - // use the thermodynamic state of the first active cell as a - // reference. there might be better ways to do this... - ElementContext elemCtx(this->ebos_simulator_); - auto elemIt = this->ebos_simulator_.gridView().template begin(); - elemCtx.updatePrimaryStencil(*elemIt); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - // Initialize a FluidState object first - FluidState fs_aquifer; - // We use the temperature of the first cell connected to the aquifer - // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs - fs_aquifer.assign(iq0.fluidState()); - Eval temperature_aq, pa0_mean, saltConcentration_aq; - temperature_aq = fs_aquifer.temperature(0); - saltConcentration_aq = fs_aquifer.saltConcentration(); - pa0_mean = this->pa0_; - Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean, saltConcentration_aq); - mu_w_ = mu_w_aquifer.value(); + this->pa0_ = this->aquct_data_.initial_pressure.value(); + this->rhow_ = this->aquct_data_.waterDensity(); } virtual Scalar aquiferDepth() const override { - return aquct_data_.d0; + return this->aquct_data_.datum_depth; } }; // class AquiferCarterTracy } // namespace Opm diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index 0b5faf21b..631f8e0ae 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -85,14 +85,16 @@ public: data.type = data::AquiferType::Fetkovich; data.aquFet = std::make_shared(); + data.aquFet->initVolume = this->aqufetp_data_.initial_watvolume; + data.aquFet->prodIndex = this->aqufetp_data_.prod_index; + data.aquFet->timeConstant = this->aqufetp_data_.timeConstant(); return data; } protected: // Aquifer Fetkovich Specific Variables - // TODO: using const reference here will cause segmentation fault, which is very strange - const Aquifetp::AQUFETP_data aqufetp_data_; + Aquifetp::AQUFETP_data aqufetp_data_; Scalar aquifer_pressure_; // aquifer void assignRestartData(const data::AquiferData& xaq) override @@ -109,52 +111,66 @@ protected: inline Eval dpai(int idx) { - const Eval dp = aquifer_pressure_ - this->pressure_current_.at(idx) - + this->rhow_[idx] * this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth()); - return dp; + const auto gdz = + this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth()); + + return this->aquifer_pressure_ + this->rhow_*gdz + - this->pressure_current_.at(idx); } // This function implements Eq 5.12 of the EclipseTechnicalDescription inline Scalar aquiferPressure() { Scalar Flux = this->W_flux_.value(); + const auto& comm = this->ebos_simulator_.vanguard().grid().comm(); comm.sum(&Flux, 1); - Scalar pa_ = this->pa0_ - Flux / (aqufetp_data_.C_t * aqufetp_data_.V0); - return pa_; + + const auto denom = + this->aqufetp_data_.total_compr * this->aqufetp_data_.initial_watvolume; + + return this->pa0_ - (Flux / denom); } inline void calculateAquiferConstants() override { - this->Tc_ = (aqufetp_data_.C_t * aqufetp_data_.V0) / aqufetp_data_.J; + this->Tc_ = this->aqufetp_data_.timeConstant(); } + // This function implements Eq 5.14 of the EclipseTechnicalDescription inline void calculateInflowRate(int idx, const Simulator& simulator) override { const Scalar td_Tc_ = simulator.timeStepSize() / this->Tc_; const Scalar coef = (1 - exp(-td_Tc_)) / td_Tc_; - this->Qai_.at(idx) = this->alphai_[idx] * aqufetp_data_.J * dpai(idx) * coef; + + this->Qai_.at(idx) = coef * this->alphai_[idx] * + this->aqufetp_data_.prod_index * dpai(idx); } inline void calculateAquiferCondition() override { - this->rhow_.resize(this->size(), 0.); - if (this->solution_set_from_restart_) { return; } - if (!aqufetp_data_.p0.first) { - this->pa0_ = this->calculateReservoirEquilibrium(); - } else { - this->pa0_ = aqufetp_data_.p0.second; + if (! this->aqufetp_data_.initial_pressure.has_value()) { + this->aqufetp_data_.initial_pressure = + this->calculateReservoirEquilibrium(); + + const auto& tables = this->ebos_simulator_.vanguard() + .eclState().getTableManager(); + + this->aqufetp_data_.finishInitialisation(tables); } - aquifer_pressure_ = this->pa0_; + + this->rhow_ = this->aqufetp_data_.waterDensity(); + this->pa0_ = this->aqufetp_data_.initial_pressure.value(); + this->aquifer_pressure_ = this->pa0_; } virtual Scalar aquiferDepth() const override { - return aqufetp_data_.d0; + return this->aqufetp_data_.datum_depth; } }; // Class AquiferFetkovich } // namespace Opm diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 88e739a77..247e1570e 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -146,7 +146,6 @@ public: // This is the pressure at td + dt this->updateCellPressure(this->pressure_current_, idx, intQuants); - this->updateCellDensity(idx, intQuants); this->calculateInflowRate(idx, context.simulator()); rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx] @@ -197,13 +196,6 @@ protected: pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); } - inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants) - { - const auto& fs = intQuants.fluidState(); - rhow_.at(idx) = fs.density(waterPhaseIdx); - } - - virtual void endTimeStep() = 0; const int aquiferID_{}; @@ -219,11 +211,11 @@ protected: std::vector pressure_previous_; std::vector pressure_current_; std::vector Qai_; - std::vector rhow_; std::vector alphai_; Scalar Tc_{}; // Time constant Scalar pa0_{}; // initial aquifer pressure + Scalar rhow_{}; Eval W_flux_; @@ -359,11 +351,13 @@ protected: const auto& fs = iq0.fluidState(); water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); - this->rhow_[idx] = fs.density(waterPhaseIdx); - pw_aquifer.push_back( - (water_pressure_reservoir - - this->rhow_[idx].value() * this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth())) - * this->alphai_[idx]); + const auto water_density = fs.density(waterPhaseIdx); + + const auto gdz = + this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth()); + + pw_aquifer.push_back(this->alphai_[idx] * + (water_pressure_reservoir - water_density.value()*gdz)); } // We take the average of the calculated equilibrium pressures.