Merge pull request #139 from GitPaean/limited_variable_update

Limited variable update

examples/test_implicit_ad.cpp

@@ -105,7 +105,7 @@ try
     Opm::LinearSolverFactory linsolver(param);
     Opm::NewtonIterationBlackoilSimple fis_solver(linsolver);
 
-    Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(*g, props, geo, 0, *wells, fis_solver);
+    Opm::FullyImplicitBlackoilSolver<UnstructuredGrid> solver(param, *g, props, geo, 0, *wells, fis_solver);
 
     Opm::BlackoilState state;
     initStateBasic(*g, props0, param, 0.0, state);

opm/autodiff/FullyImplicitBlackoilSolver.hpp

@@ -31,6 +31,7 @@ struct Wells;
 
 namespace Opm {
 
+    namespace parameter { class ParameterGroup; }
     class DerivedGeology;
     class RockCompressibility;
     class NewtonIterationBlackoilInterface;
@@ -56,13 +57,15 @@ namespace Opm {
         /// Construct a solver. It will retain references to the
         /// arguments of this functions, and they are expected to
         /// remain in scope for the lifetime of the solver.
+        /// \param[in] param       parameters
         /// \param[in] grid             grid data structure
         /// \param[in] fluid            fluid properties
         /// \param[in] geo              rock properties
         /// \param[in] rock_comp_props  if non-null, rock compressibility properties
         /// \param[in] wells            well structure
         /// \param[in] linsolver        linear solver
-        FullyImplicitBlackoilSolver(const Grid&                     grid ,
+        FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param,
+                                    const Grid&                     grid ,
                                     const BlackoilPropsAdInterface& fluid,
                                     const DerivedGeology&           geo  ,
                                     const RockCompressibility*      rock_comp_props,
@@ -137,6 +140,9 @@ namespace Opm {
         HelperOps                       ops_;
         const WellOps                   wops_;
         const M                         grav_;
+        double                          dp_max_rel_;
+        double                          ds_max_;
+        double                          drs_max_rel_;
 
         std::vector<ReservoirResidualQuant> rq_;
         std::vector<PhasePresence> phaseCondition_;
@@ -266,6 +272,9 @@ namespace Opm {
         /// residual mass balance (tol_cnv).
         bool getConvergence(const double dt);
 
+        const double dpMaxRel() const { return dp_max_rel_; }
+        const double dsMax() const { return ds_max_; }
+        const double drsMaxRel() const { return drs_max_rel_; }
 
     };
 } // namespace Opm

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -35,6 +35,7 @@
 #include <opm/core/utility/Exceptions.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/core/well_controls.h>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 
 #include <cassert>
 #include <cmath>
@@ -212,7 +213,8 @@ namespace {
 
     template<class T>
     FullyImplicitBlackoilSolver<T>::
-    FullyImplicitBlackoilSolver(const Grid&                     grid ,
+    FullyImplicitBlackoilSolver(const parameter::ParameterGroup& param,
+                                const Grid&                     grid ,
                                 const BlackoilPropsAdInterface& fluid,
                                 const DerivedGeology&           geo  ,
                                 const RockCompressibility*      rock_comp_props,
@@ -230,17 +232,21 @@ namespace {
         , ops_   (grid)
         , wops_  (wells)
         , grav_  (gravityOperator(grid_, ops_, geo_))
+        , dp_max_rel_ ( 1.0e9 )
+        , ds_max_ ( 0.2 )
+        , drs_max_rel_ ( 1.0e9 )
         , rq_    (fluid.numPhases())
         , phaseCondition_(AutoDiffGrid::numCells(grid))
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
                         ADB::null(),
                         ADB::null() } )
     {
+        dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_);
+        ds_max_ = param.getDefault("ds_max", ds_max_);
+        drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_);
     }
 
 
-
-
     template<class T>
     void
     FullyImplicitBlackoilSolver<T>::
@@ -1230,7 +1236,7 @@ namespace {
         assert(varstart == dx.size());
 
         // Pressure update.
-        const double dpmaxrel = 0.8;
+        const double dpmaxrel = dpMaxRel();
         const V p_old = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
         const V absdpmax = dpmaxrel*p_old.abs();
         const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
@@ -1241,40 +1247,53 @@ namespace {
         // Saturation updates.
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
         const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
-        const double dsmax = 0.3;
+        const double dsmax = dsMax();
         V so = one;
         V sw;
+        V sg;
 
-        if (active_[ Water ]) {
+        {
+            V maxVal = zero;
+            V dso = zero;
+            if (active_[Water]){
+                maxVal = dsw.abs().max(maxVal);
+                dso = dso - dsw;
+            }
+
+            V dsg;
+            if (active_[Gas]){
+                dsg = isSg * dxvar - isRv * dsw;
+                maxVal = dsg.abs().max(maxVal);
+                dso = dso - dsg;
+            }
+
+            maxVal = dso.abs().max(maxVal);
+
+            V step = dsmax/maxVal;
+            step = step.min(1.);
+
+            if (active_[Water]) {
                 const int pos = pu.phase_pos[ Water ];
                 const V sw_old = s_old.col(pos);
-            const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax);
+                sw = sw_old - step * dsw;
-            sw = (sw_old - dsw_limited).unaryExpr(Chop01());
                 so -= sw;
-            for (int c = 0; c < nc; ++c) {
-                state.saturation()[c*np + pos] = sw[c];
-            }
             }
 
-        V sg;
             if (active_[Gas]) {
                 const int pos = pu.phase_pos[ Gas ];
                 const V sg_old = s_old.col(pos);
-            const V dsg = isSg * dxvar - isRv * dsw;
+                sg = sg_old - step * dsg;
-            const V dsg_limited = sign(dsg) * dsg.abs().min(dsmax);
-            sg = sg_old - dsg_limited;
                 so -= sg;
             }
+        }
 
-
+        const double drsmaxrel = drsMaxRel();
-
-        const double drsmax = 1e9;
         const double drvmax = 1e9;//% same as in Mrst
         V rs;
         if (disgas) {
             const V rs_old = Eigen::Map<const V>(&state.gasoilratio()[0], nc);
             const V drs = isRs * dxvar;
-            const V drs_limited = sign(drs) * drs.abs().min(drsmax);
+            const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel);
             rs = rs_old - drs_limited;
         }
         V rv;
@@ -1410,7 +1429,8 @@ namespace {
 
         // Bhp update.
         const V bhp_old = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
-        const V bhp = bhp_old - dbhp;
+        const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
+        const V bhp = bhp_old - dbhp_limited;
         std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
 
     }

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -200,7 +200,7 @@ namespace Opm
           wells_(wells_manager.c_wells()),
           gravity_(gravity),
           geo_(grid_, props_, gravity_),
-          solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver)
+          solver_(param, grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver)
           /*                   param.getDefault("nl_pressure_residual_tolerance", 0.0),
                                param.getDefault("nl_pressure_change_tolerance", 1.0),
                                param.getDefault("nl_pressure_maxiter", 10),