Merge pull request #3254 from bska/chase-aquifer-redefinition

Use Aquifer's Notion of Water Properties

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -84,6 +84,10 @@ public:
         data.type = data::AquiferType::CarterTracy;
 
         data.aquCT = std::make_shared<data::CarterTracyData>();
+        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,
+            linearInterpolation(this->aquct_data_.dimensionless_time,
-                                this->aquct_data_.pi,
+                                this->aquct_data_.dimensionless_pressure,
                                 this->dimensionless_time_);
 
         const auto PItd =
-            linearInterpolation(this->aquct_data_.td,
+            linearInterpolation(this->aquct_data_.dimensionless_time,
-                                this->aquct_data_.pi, td_plus_dt);
+                                this->aquct_data_.dimensionless_pressure,
+                                td_plus_dt);
 
         const auto PItdprime =
-            linearInterpolationDerivative(this->aquct_data_.td,
+            linearInterpolationDerivative(this->aquct_data_.dimensionless_time,
-                                          this->aquct_data_.pi, td_plus_dt);
+                                          this->aquct_data_.dimensionless_pressure,
+                                          td_plus_dt);
 
         return std::make_pair(PItd, PItdprime);
     }
 
     Scalar dpai(const int idx) const
     {
-        Scalar dp = this->pa0_
+        const auto gdz =
-            + this->rhow_.at(idx).value() * this->gravity_() * (this->cell_depth_.at(idx) - this->aquiferDepth())
+            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
+        this->Tc_ = this->aquct_data_.timeConstant();
-        beta_ = aquct_data_.c2 * aquct_data_.h * aquct_data_.theta * aquct_data_.phi_aq * aquct_data_.C_t
+        this->beta_ = this->aquct_data_.influxConstant();
-            * 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);
     }
 
     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;
+            const auto& tables = this->ebos_simulator_.vanguard()
-        this->rhow_.resize(this->size(), 0.);
+                .eclState().getTableManager();
-        if (!aquct_data_.p0.first) {
+
-            this->pa0_ = this->calculateReservoirEquilibrium();
+            this->aquct_data_.finishInitialisation(tables);
-        } else {
-            this->pa0_ = aquct_data_.p0.second;
         }
 
-        // use the thermodynamic state of the first active cell as a
+        this->pa0_ = this->aquct_data_.initial_pressure.value();
-        // reference. there might be better ways to do this...
+        this->rhow_ = this->aquct_data_.waterDensity();
-        ElementContext elemCtx(this->ebos_simulator_);
-        auto elemIt = this->ebos_simulator_.gridView().template begin</*codim=*/0>();
-        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();
     }
 
     virtual Scalar aquiferDepth() const override
     {
-        return aquct_data_.d0;
+        return this->aquct_data_.datum_depth;
     }
 }; // class AquiferCarterTracy
 } // namespace Opm

opm/simulators/aquifers/AquiferFetkovich.hpp

@@ -85,14 +85,16 @@ public:
         data.type = data::AquiferType::Fetkovich;
 
         data.aquFet = std::make_shared<data::FetkovichData>();
+        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
+    Aquifetp::AQUFETP_data aqufetp_data_;
-    const 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)
+        const auto gdz =
-            + this->rhow_[idx] * this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
+            this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
-        return dp;
+
+        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) {
+        if (! this->aqufetp_data_.initial_pressure.has_value()) {
-            this->pa0_ = this->calculateReservoirEquilibrium();
+            this->aqufetp_data_.initial_pressure =
-        } else {
+                this->calculateReservoirEquilibrium();
-            this->pa0_ = aqufetp_data_.p0.second;
+
+            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

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<Scalar> pressure_previous_;
     std::vector<Eval> pressure_current_;
     std::vector<Eval> Qai_;
-    std::vector<Eval> rhow_;
     std::vector<Scalar> 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);
+            const auto water_density = fs.density(waterPhaseIdx);
-            pw_aquifer.push_back(
+
-                (water_pressure_reservoir
+            const auto gdz =
-                 - this->rhow_[idx].value() * this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth()))
+                this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());
-                * this->alphai_[idx]);
+
+            pw_aquifer.push_back(this->alphai_[idx] *
+                (water_pressure_reservoir - water_density.value()*gdz));
         }
 
         // We take the average of the calculated equilibrium pressures.