From 578fa4b6c38e5ff2587568f05c7e855aae0135bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 22 Mar 2021 15:46:14 +0100 Subject: [PATCH 1/2] Report Dimensionless Time and Pressure for CT Aquifers This commit adds support for calculating and reporting the dimensionless time (simulator time divided by aquifer's time constant) and pressure (influence function evaluated at dimensionless time) values as part of the Carter-Tracy aquifer's 'aquiferData()' reporting function. These values are useful in their own right, e.g., for summary output through the keywords AAQTD and AAQPD, but they are also needed for ECLIPSE restart purposes. --- .../aquifers/AquiferCarterTracy.hpp | 40 ++++++++++++++----- 1 file changed, 30 insertions(+), 10 deletions(-) diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index e033da5fc..f68cb29ac 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -26,6 +26,7 @@ #include #include +#include #include namespace Opm @@ -80,7 +81,12 @@ public: } data.volume = this->W_flux_.value(); data.initPressure = this->pa0_; - data.type = Opm::data::AquiferType::CarterTracey; + data.type = Opm::data::AquiferType::CarterTracy; + + data.aquCT = std::make_shared(); + data.aquCT->dimensionless_time = this->dimensionless_time_; + data.aquCT->dimensionless_pressure = this->dimensionless_pressure_; + return data; } @@ -92,17 +98,27 @@ protected: Scalar mu_w_{1}; // water viscosity Scalar fluxValue_{0}; // value of flux + Scalar dimensionless_time_{0}; + Scalar dimensionless_pressure_{0}; + void assignRestartData(const data::AquiferData& /* xaq */) override { throw std::runtime_error {"Restart-based initialization not currently supported " "for Carter-Tracey analytic aquifers"}; } - inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) + inline void getInfluenceTableValues(const Scalar td_plus_dt, + Scalar& PItd, + Scalar& PItdprime) { // We use the opm-common numeric linear interpolator - pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td); - pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td); + this->dimensionless_pressure_ = + Opm::linearInterpolation(this->aquct_data_.td, + this->aquct_data_.pi, + this->dimensionless_time_); + + PItd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td_plus_dt); + PItdprime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td_plus_dt); } inline Scalar dpai(int idx) @@ -117,12 +133,16 @@ protected: inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator) { const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / this->Tc_; - const Scalar td = simulator.time() / this->Tc_; - Scalar PItdprime = 0.; - Scalar PItd = 0.; - getInfluenceTableValues(PItd, PItdprime, td_plus_dt); - a = 1.0 / this->Tc_ * ((beta_ * dpai(idx)) - (this->fluxValue_ * PItdprime)) / (PItd - td * PItdprime); - b = beta_ / (this->Tc_ * (PItd - td * PItdprime)); + this->dimensionless_time_ = simulator.time() / this->Tc_; + + auto PItd = Scalar{0}; + auto PItdprime = Scalar{0}; + this->getInfluenceTableValues(td_plus_dt, PItd, PItdprime); + + const auto denom = this->Tc_ * (PItd - this->dimensionless_time_*PItdprime); + + a = (this->beta_*dpai(idx) - this->fluxValue_*PItdprime) / denom; + b = this->beta_ / denom; } // This function implements Eq 5.7 of the EclipseTechnicalDescription From 7503cfd76ada9ef01a710518eaf2972cf8109e01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 24 Mar 2021 18:21:30 +0100 Subject: [PATCH 2/2] Address PR Review Comments --- .../aquifers/AquiferCarterTracy.hpp | 39 +++++++++++-------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index f68cb29ac..4ab6534d0 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -28,6 +28,7 @@ #include #include #include +#include namespace Opm { @@ -107,9 +108,8 @@ protected: "for Carter-Tracey analytic aquifers"}; } - inline void getInfluenceTableValues(const Scalar td_plus_dt, - Scalar& PItd, - Scalar& PItdprime) + std::pair + getInfluenceTableValues(const Scalar td_plus_dt) { // We use the opm-common numeric linear interpolator this->dimensionless_pressure_ = @@ -117,11 +117,18 @@ protected: this->aquct_data_.pi, this->dimensionless_time_); - PItd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td_plus_dt); - PItdprime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td_plus_dt); + const auto PItd = + Opm::linearInterpolation(this->aquct_data_.td, + this->aquct_data_.pi, td_plus_dt); + + const auto PItdprime = + Opm::linearInterpolationDerivative(this->aquct_data_.td, + this->aquct_data_.pi, td_plus_dt); + + return std::make_pair(PItd, PItdprime); } - inline Scalar dpai(int idx) + Scalar dpai(const int idx) const { Scalar dp = this->pa0_ + this->rhow_.at(idx).value() * this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth()) @@ -130,28 +137,28 @@ protected: } // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription - inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const Simulator& simulator) + std::pair + calculateEqnConstants(const int idx, const Simulator& simulator) { const Scalar td_plus_dt = (simulator.timeStepSize() + simulator.time()) / this->Tc_; this->dimensionless_time_ = simulator.time() / this->Tc_; - auto PItd = Scalar{0}; - auto PItdprime = Scalar{0}; - this->getInfluenceTableValues(td_plus_dt, PItd, PItdprime); + const auto [PItd, PItdprime] = this->getInfluenceTableValues(td_plus_dt); const auto denom = this->Tc_ * (PItd - this->dimensionless_time_*PItdprime); + const auto a = (this->beta_*dpai(idx) - this->fluxValue_*PItdprime) / denom; + const auto b = this->beta_ / denom; - a = (this->beta_*dpai(idx) - this->fluxValue_*PItdprime) / denom; - b = this->beta_ / denom; + return std::make_pair(a, b); } // This function implements Eq 5.7 of the EclipseTechnicalDescription inline void calculateInflowRate(int idx, const Simulator& simulator) override { - Scalar a, b; - calculateEqnConstants(a, b, idx, simulator); - this->Qai_.at(idx) - = this->alphai_.at(idx) * (a - b * (this->pressure_current_.at(idx) - this->pressure_previous_.at(idx))); + const auto [a, b] = this->calculateEqnConstants(idx, simulator); + + this->Qai_.at(idx) = this->alphai_.at(idx) * + (a - b*(this->pressure_current_.at(idx) - this->pressure_previous_.at(idx))); } inline void calculateAquiferConstants() override