Merge pull request #3531 from totto82/aquctAndCo2Store

support aquct and co2store

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -50,8 +50,6 @@ public:
     using typename Base::Simulator;
     using typename Base::ElementMapper;
 
-    using Base::waterCompIdx;
-    using Base::waterPhaseIdx;
     AquiferCarterTracy(const std::vector<Aquancon::AquancCell>& connections,
                        const Simulator& ebosSimulator,
                        const AquiferCT::AQUCT_data& aquct_data)
@@ -83,12 +81,21 @@ public:
         data.initPressure = this->pa0_;
 
         auto* aquCT = data.typeData.template create<data::AquiferType::CarterTracy>();
-        aquCT->timeConstant = this->aquct_data_.timeConstant();
-        aquCT->influxConstant = this->aquct_data_.influxConstant();
-        aquCT->waterDensity = this->aquct_data_.waterDensity();
-        aquCT->waterViscosity = this->aquct_data_.waterViscosity();
+
         aquCT->dimensionless_time = this->dimensionless_time_;
         aquCT->dimensionless_pressure = this->dimensionless_pressure_;
+        aquCT->influxConstant = this->aquct_data_.influxConstant();
+
+        if (!this->co2store_()) {
+            aquCT->timeConstant = this->aquct_data_.timeConstant();
+            aquCT->waterDensity = this->aquct_data_.waterDensity();
+            aquCT->waterViscosity = this->aquct_data_.waterViscosity();
+        } else {
+            aquCT->waterDensity = this->rhow_;
+            aquCT->timeConstant = this->Tc_;
+            const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
+            aquCT->waterViscosity = this->Tc_ *  this->aquct_data_.permeability / x;
+        }
 
         return data;
     }
@@ -159,6 +166,11 @@ protected:
         return std::make_pair(a, b);
     }
 
+    std::size_t pvtRegionIdx() const
+    {
+        return this->aquct_data_.pvttableID - 1;
+    }
+
     // This function implements Eq 5.7 of the EclipseTechnicalDescription
     inline void calculateInflowRate(int idx, const Simulator& simulator) override
     {
@@ -170,7 +182,19 @@ protected:
 
     inline void calculateAquiferConstants() override
     {
-        this->Tc_ = this->aquct_data_.timeConstant();
+        if(this->co2store_()) {
+             const auto press = this->aquct_data_.initial_pressure.value();
+             Scalar temp = FluidSystem::reservoirTemperature();
+             if (this->aquct_data_.initial_temperature.has_value())
+                 temp = this->aquct_data_.initial_temperature.value();
+
+             Scalar rs = 0.0; // no dissolved CO2
+             Scalar waterViscosity = FluidSystem::oilPvt().viscosity(pvtRegionIdx(), temp, press, rs);
+             const auto x = this->aquct_data_.porosity * this->aquct_data_.total_compr * this->aquct_data_.inner_radius * this->aquct_data_.inner_radius;
+             this->Tc_  = waterViscosity * x / this->aquct_data_.permeability;
+        } else {
+             this->Tc_ = this->aquct_data_.timeConstant();
+        }
         this->beta_ = this->aquct_data_.influxConstant();
     }
 
@@ -191,7 +215,20 @@ protected:
         }
 
         this->pa0_ = this->aquct_data_.initial_pressure.value();
-        this->rhow_ = this->aquct_data_.waterDensity();
+        if(this->co2store_()) {
+             const auto press = this->aquct_data_.initial_pressure.value();
+
+             Scalar temp = FluidSystem::reservoirTemperature();
+             if (this->aquct_data_.initial_temperature.has_value())
+                 temp = this->aquct_data_.initial_temperature.value();
+
+             Scalar rs = 0.0; // no dissolved CO2
+             Scalar waterDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx(), temp, press, rs)
+                     * FluidSystem::oilPvt().oilReferenceDensity(pvtRegionIdx());
+             this->rhow_  = waterDensity;
+        } else {
+            this->rhow_ = this->aquct_data_.waterDensity();
+        }
     }
 
     virtual Scalar aquiferDepth() const override

opm/simulators/aquifers/AquiferFetkovich.hpp

@@ -50,9 +50,6 @@ public:
     using typename Base::Simulator;
     using typename Base::ElementMapper;
 
-    using Base::waterCompIdx;
-    using Base::waterPhaseIdx;
-
     AquiferFetkovich(const std::vector<Aquancon::AquancCell>& connections,
                      const Simulator& ebosSimulator,
                      const Aquifetp::AQUFETP_data& aqufetp_data)

opm/simulators/aquifers/AquiferInterface.hpp

@@ -72,9 +72,6 @@ public:
                                BlackoilIndices::numPhases>
         FluidState;
 
-    static const auto waterCompIdx = FluidSystem::waterCompIdx;
-    static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
-
     // Constructor
     AquiferInterface(int aqID,
                      const std::vector<Aquancon::AquancCell>& connections,
@@ -125,7 +122,7 @@ public:
 
             elemCtx.updateIntensiveQuantities(0);
             const auto& iq = elemCtx.intensiveQuantities(0, 0);
-            pressure_previous_[idx] = getValue(iq.fluidState().pressure(waterPhaseIdx));
+            pressure_previous_[idx] = getValue(iq.fluidState().pressure(phaseIdx_()));
         }
     }
 
@@ -150,7 +147,7 @@ public:
         this->updateCellPressure(this->pressure_current_, idx, intQuants);
         this->calculateInflowRate(idx, context.simulator());
 
-        rates[BlackoilIndices::conti0EqIdx + FluidSystem::waterCompIdx]
+        rates[BlackoilIndices::conti0EqIdx + compIdx_()]
             += this->Qai_[idx] / context.dofVolume(spaceIdx, timeIdx);
     }
 
@@ -166,6 +163,28 @@ protected:
         return ebos_simulator_.problem().gravity()[2];
     }
 
+    inline bool co2store_() const
+    {
+        return ebos_simulator_.vanguard().eclState().runspec().co2Storage();
+    }
+
+    inline bool phaseIdx_() const
+    {
+        if(co2store_())
+            return FluidSystem::oilPhaseIdx;
+
+        return FluidSystem::waterPhaseIdx;
+    }
+
+    inline bool compIdx_() const
+    {
+        if(co2store_())
+            return FluidSystem::oilCompIdx;
+
+        return FluidSystem::waterCompIdx;
+    }
+
+
     inline void initQuantities()
     {
         // We reset the cumulative flux at the start of any simulation, so, W_flux = 0
@@ -188,14 +207,14 @@ protected:
     updateCellPressure(std::vector<Eval>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
     {
         const auto& fs = intQuants.fluidState();
-        pressure_water.at(idx) = fs.pressure(waterPhaseIdx);
+        pressure_water.at(idx) = fs.pressure(phaseIdx_());
     }
 
     inline void
     updateCellPressure(std::vector<Scalar>& pressure_water, const int idx, const IntensiveQuantities& intQuants)
     {
         const auto& fs = intQuants.fluidState();
-        pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value();
+        pressure_water.at(idx) = fs.pressure(phaseIdx_()).value();
     }
 
     virtual void endTimeStep() = 0;
@@ -380,8 +399,8 @@ protected:
             const auto& iq0 = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
             const auto& fs = iq0.fluidState();
 
-            water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
-            const auto water_density = fs.density(waterPhaseIdx);
+            water_pressure_reservoir = fs.pressure(phaseIdx_()).value();
+            const auto water_density = fs.density(phaseIdx_());
 
             const auto gdz =
                 this->gravity_() * (this->cell_depth_[idx] - this->aquiferDepth());