From 21d69a7af0ac3737818d53e5a887cbce4193df27 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Nov 2019 13:24:16 +0100 Subject: [PATCH 1/4] making dp is Eval for AquiferFetkovich --- opm/simulators/aquifers/AquiferFetkovich.hpp | 13 +++++++------ opm/simulators/aquifers/AquiferInterface.hpp | 2 +- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index b01c71bbf..a0d0e07e7 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -145,9 +145,10 @@ namespace Opm } } - inline Scalar dpai(int idx) + inline Eval dpai(int idx) { - Scalar dp = aquifer_pressure_ + Base::rhow_.at(idx).value()*Base::gravity_()*(Base::cell_depth_.at(idx) - aqufetp_data_.d0) - Base::pressure_current_.at(idx).value() ; + const Eval dp = aquifer_pressure_ - Base::pressure_current_.at(idx) + + Base::rhow_.at(idx) * Base::gravity_()*(Base::cell_depth_.at(idx) - aqufetp_data_.d0); return dp; } @@ -166,14 +167,14 @@ namespace Opm // This function implements Eq 5.14 of the EclipseTechnicalDescription inline void calculateInflowRate(int idx, const Simulator& simulator) { - Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_ ; - Scalar exp_ = (1 - exp(-td_Tc_)) / td_Tc_; - Base::Qai_.at(idx) = Base::alphai_.at(idx) * aqufetp_data_.J * dpai(idx) * exp_; + const Scalar td_Tc_ = simulator.timeStepSize() / Base::Tc_ ; + const Scalar coef = (1 - exp(-td_Tc_)) / td_Tc_; + Base::Qai_.at(idx) = Base::alphai_[idx] * aqufetp_data_.J * dpai(idx) * coef; } inline void calculateAquiferCondition() { - int pvttableIdx = aqufetp_data_.pvttableID - 1; + const int pvttableIdx = aqufetp_data_.pvttableID - 1; Base::rhow_.resize(Base::cell_idx_.size(),0.); if (!aqufetp_data_.p0) { diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 5a8e8cf94..8aaa96596 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -210,7 +210,7 @@ namespace Opm virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0; - virtual Scalar dpai(int idx) = 0; + // virtual Scalar dpai(int idx) = 0; virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0; From 383748b1616cfdab8932e0604c39f12e3154ed4e Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Nov 2019 13:24:34 +0100 Subject: [PATCH 2/4] some cleaning up for the AquiferFetkovich --- opm/simulators/aquifers/AquiferFetkovich.hpp | 19 +------------------ opm/simulators/aquifers/AquiferInterface.hpp | 1 + 2 files changed, 2 insertions(+), 18 deletions(-) diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index a0d0e07e7..3edd0389f 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -64,6 +64,7 @@ namespace Opm protected: // Aquifer Fetkovich Specific Variables + // TODO: using const reference here will cause segmentation fault, which is very strange const Aquifetp::AQUFETP_data aqufetp_data_; Scalar aquifer_pressure_; // aquifer @@ -185,24 +186,6 @@ namespace Opm Base::pa0_ = *(aqufetp_data_.p0); } aquifer_pressure_ = Base::pa0_ ; - - // use the thermodynamic state of the first active cell as a - // reference. there might be better ways to do this... - ElementContext elemCtx(Base::ebos_simulator_); - auto elemIt = Base::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; - temperature_aq = fs_aquifer.temperature(0); - pa0_mean = Base::pa0_; - Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); - Base::mu_w_ = mu_w_aquifer.value(); } inline Scalar calculateReservoirEquilibrium() diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 8aaa96596..85d13e406 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -202,6 +202,7 @@ namespace Opm std::vector rhow_; std::vector alphai_; + // TODO: Fetkovich does not use mu_w_, it should be moved to CT class Scalar mu_w_; //water viscosity Scalar Tc_; // Time constant Scalar pa0_; // initial aquifer pressure From eea833ced2a9da1f98931b544f61659c7e1a54a6 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Nov 2019 13:25:22 +0100 Subject: [PATCH 3/4] correcting the way of calculate initial aquifer pressure for AquiferFetkovich. --- opm/simulators/aquifers/AquiferFetkovich.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index 3edd0389f..ddcb92aec 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -216,7 +216,8 @@ namespace Opm } // We take the average of the calculated equilibrium pressures. - Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.)/pw_aquifer.size(); + const Scalar sum_alpha = std::accumulate(this->alphai_.begin(), this->alphai_.end(), 0.); + const Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.) / sum_alpha; return aquifer_pres_avg; } }; //Class AquiferFetkovich From 0c6adbbba0a7925221d690dc72a1552ac73b8811 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 14 Nov 2019 22:53:39 +0100 Subject: [PATCH 4/4] some small cleaning up for Aquifer models there should be no functional change. --- opm/simulators/aquifers/AquiferCarterTracy.hpp | 7 +++++-- opm/simulators/aquifers/AquiferFetkovich.hpp | 3 +-- opm/simulators/aquifers/AquiferInterface.hpp | 4 ---- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/opm/simulators/aquifers/AquiferCarterTracy.hpp b/opm/simulators/aquifers/AquiferCarterTracy.hpp index 61de9108d..714576cab 100644 --- a/opm/simulators/aquifers/AquiferCarterTracy.hpp +++ b/opm/simulators/aquifers/AquiferCarterTracy.hpp @@ -63,6 +63,8 @@ namespace Opm // Variables constants const AquiferCT::AQUCT_data aquct_data_; Scalar beta_; // Influx constant + // TODO: it is possible it should be a AD variable + Scalar mu_w_; // water viscosity // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer inline void initializeConnections(const Aquancon::AquanconOutput& connection) @@ -184,7 +186,7 @@ namespace Opm * aquct_data_.C_t * aquct_data_.r_o * aquct_data_.r_o; // We calculate the time constant - Base::Tc_ = Base::mu_w_ * aquct_data_.phi_aq + Base::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 ); @@ -220,11 +222,12 @@ namespace Opm temperature_aq = fs_aquifer.temperature(0); pa0_mean = Base::pa0_; Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); - Base::mu_w_ = mu_w_aquifer.value(); + mu_w_ = mu_w_aquifer.value(); } // This function is for calculating the aquifer properties from equilibrium state with the reservoir + // TODO: this function can be moved to the Inteface class, since it is the same for both Aquifer models inline Scalar calculateReservoirEquilibrium() { // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices diff --git a/opm/simulators/aquifers/AquiferFetkovich.hpp b/opm/simulators/aquifers/AquiferFetkovich.hpp index ddcb92aec..4cd873b72 100644 --- a/opm/simulators/aquifers/AquiferFetkovich.hpp +++ b/opm/simulators/aquifers/AquiferFetkovich.hpp @@ -149,7 +149,7 @@ namespace Opm inline Eval dpai(int idx) { const Eval dp = aquifer_pressure_ - Base::pressure_current_.at(idx) - + Base::rhow_.at(idx) * Base::gravity_()*(Base::cell_depth_.at(idx) - aqufetp_data_.d0); + + Base::rhow_[idx] * Base::gravity_()*(Base::cell_depth_[idx] - aqufetp_data_.d0); return dp; } @@ -175,7 +175,6 @@ namespace Opm inline void calculateAquiferCondition() { - const int pvttableIdx = aqufetp_data_.pvttableID - 1; Base::rhow_.resize(Base::cell_idx_.size(),0.); if (!aqufetp_data_.p0) { diff --git a/opm/simulators/aquifers/AquiferInterface.hpp b/opm/simulators/aquifers/AquiferInterface.hpp index 85d13e406..3a9fb115c 100644 --- a/opm/simulators/aquifers/AquiferInterface.hpp +++ b/opm/simulators/aquifers/AquiferInterface.hpp @@ -202,8 +202,6 @@ namespace Opm std::vector rhow_; std::vector alphai_; - // TODO: Fetkovich does not use mu_w_, it should be moved to CT class - Scalar mu_w_; //water viscosity Scalar Tc_; // Time constant Scalar pa0_; // initial aquifer pressure @@ -211,8 +209,6 @@ namespace Opm virtual void initializeConnections(const Aquancon::AquanconOutput& connection) =0; - // virtual Scalar dpai(int idx) = 0; - virtual void calculateInflowRate(int idx, const Simulator& simulator) = 0; virtual void calculateAquiferCondition() = 0;