Merge pull request #2152 from GitPaean/wip_fixing_aquifer_2

fixing and improvement of the Aquifer models

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

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
 
@@ -145,9 +146,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_[idx] * Base::gravity_()*(Base::cell_depth_[idx] - aqufetp_data_.d0);
       return dp;
     }
 
@@ -166,14 +168,13 @@ 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;
       Base::rhow_.resize(Base::cell_idx_.size(),0.);
       if (!aqufetp_data_.p0)
       {
@@ -184,24 +185,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</*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;
-      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()
@@ -232,7 +215,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

opm/simulators/aquifers/AquiferInterface.hpp

@@ -202,7 +202,6 @@ namespace Opm
     std::vector<Eval> rhow_;
     std::vector<Scalar> alphai_;
 
-    Scalar mu_w_; //water viscosity
     Scalar Tc_; // Time constant
     Scalar pa0_; // initial aquifer pressure
 
@@ -210,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;