Use Aquifer's Notion of Water Properties

This commit switches to using the analytic aquifer's intrinsic water
properties (i.e., the mass density and the viscosity), and to get
the time constant from the *_data structure instead of calculating
this value with separate logic.  Note that this switches to using a
single density value for the aquifer instead of separate density
values for each aquifer connection.

If the aquifer's initial pressure is defaulted we still compute an
equilibrated initial pressure value.  We then use the *_data
structure's 'finishInitialisation()' member function to derive the
aforementioned PVT properties.

Finally, report these values in the aquifer type-specific sub
structures of data::AquiferData for restart output purposes.

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,
-                                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</*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();
+        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

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
-    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

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);
-            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.