Merge pull request #880 from totto82/appelyard_fixes

Fixes in the Appelyard in updateState and updateWellState

opm/autodiff/BlackoilModelBase.hpp

@@ -540,6 +540,7 @@ namespace Opm {
                              int nc) const;
 
         double dpMaxRel() const { return param_.dp_max_rel_; }
+        double dbhpMaxRel() const {return param_.dbhp_max_rel_; }
         double dsMax() const { return param_.ds_max_; }
         double drMaxRel() const { return param_.dr_max_rel_; }
         double maxResidualAllowed() const { return param_.max_residual_allowed_; }

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -1064,7 +1064,7 @@ typedef Eigen::Array<double,
                 ADB::V total_residual_v = total_residual.value();
                 const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
                 assert(dx.size() == total_residual_v.size());
-                asImpl().wellModel().updateWellState(dx.array(), dpMaxRel(), well_state);
+                asImpl().wellModel().updateWellState(dx.array(), dbhpMaxRel(), well_state);
             }
             // We have to update the well controls regardless whether there are local
             // wells active or not as parallel logging will take place that needs to
@@ -1141,6 +1141,7 @@ typedef Eigen::Array<double,
         const V null;
         assert(null.size() == 0);
         const V zero = V::Zero(nc);
+        const V ones = V::Constant(nc,1.0);
 
         // Extract parts of dx corresponding to each part.
         const V dp = subset(dx, Span(nc));
@@ -1165,7 +1166,6 @@ typedef Eigen::Array<double,
         const V p = (p_old - dp_limited).max(zero);
         std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
 
-
         // Saturation updates.
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
         const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
@@ -1213,7 +1213,6 @@ typedef Eigen::Array<double,
             so = so_old - step * dso;
         }
 
-        // Appleyard chop process.
         if (active_[Gas]) {
             auto ixg = sg < 0;
             for (int c = 0; c < nc; ++c) {
@@ -1255,45 +1254,23 @@ typedef Eigen::Array<double,
             }
         }
 
-        //const V sumSat = sw + so + sg;
-        //sw = sw / sumSat;
-        //so = so / sumSat;
-        //sg = sg / sumSat;
-
-        // Update the reservoir_state
-        if (active_[Water]) {
-            for (int c = 0; c < nc; ++c) {
-                reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
-            }
-        }
-
-        if (active_[Gas]) {
-            for (int c = 0; c < nc; ++c) {
-                reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
-            }
-        }
-
-        if (active_[ Oil ]) {
-            for (int c = 0; c < nc; ++c) {
-                reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
-            }
-        }
-
         // Update rs and rv
         const double drmaxrel = drMaxRel();
         V rs;
         if (has_disgas_) {
             const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
             const V drs = isRs_ * dxvar;
-            const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
+            const V drs_limited = sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
             rs = rs_old - drs_limited;
+            rs = rs.max(zero);
         }
         V rv;
         if (has_vapoil_) {
             const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
             const V drv = isRv_ * dxvar;
-            const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
+            const V drv_limited = sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
             rv = rv_old - drv_limited;
+            rv = rv.max(zero);
         }
 
         // Sg is used as primal variable for water only cells.
@@ -1318,11 +1295,16 @@ typedef Eigen::Array<double,
             for (int c = 0; c < nc; ++c) {
                 if (useSg[c]) {
                     rs[c] = rsSat[c];
+                    if (watOnly[c]) {
+                        so[c] = 0;
+                        sg[c] = 0;
+                        rs[c] = 0;
+                    }
                 } else {
                     hydroCarbonState[c] = HydroCarbonState::OilOnly;
                 }
             }
-
+            //rs = rs.min(rsSat);
         }
 
         // phase transitions so <-> rv
@@ -1345,14 +1327,37 @@ typedef Eigen::Array<double,
             for (int c = 0; c < nc; ++c) {
                 if (useSg[c]) {
                     rv[c] = rvSat[c];
+                    if (watOnly[c]) {
+                        so[c] = 0;
+                        sg[c] = 0;
+                        rv[c] = 0;
+                    }
                 } else {
                     hydroCarbonState[c] = HydroCarbonState::GasOnly;
                 }
             }
-
+            //rv = rv.min(rvSat);
         }
 
         // Update the reservoir_state
+        if (active_[Water]) {
+            for (int c = 0; c < nc; ++c) {
+                reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
+            }
+        }
+
+        if (active_[Gas]) {
+            for (int c = 0; c < nc; ++c) {
+                reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
+            }
+        }
+
+        if (active_[ Oil ]) {
+            for (int c = 0; c < nc; ++c) {
+                reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
+            }
+        }
+
         if (has_disgas_) {
             std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
         }
@@ -1361,8 +1366,7 @@ typedef Eigen::Array<double,
             std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
         }
 
-
+        asImpl().wellModel().updateWellState(dwells, dbhpMaxRel(), well_state);
-        asImpl().wellModel().updateWellState(dwells, dpMaxRel(), well_state);
 
         // Update phase conditions used for property calculations.
         updatePhaseCondFromPrimalVariable(reservoir_state);

opm/autodiff/BlackoilModelParameters.cpp

@@ -21,7 +21,6 @@
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/parser/eclipse/Units/Units.hpp>
 
-
 namespace Opm
 {
 
@@ -44,6 +43,8 @@ namespace Opm
         dp_max_rel_  = param.getDefault("dp_max_rel", dp_max_rel_);
         ds_max_      = param.getDefault("ds_max", ds_max_);
         dr_max_rel_  = param.getDefault("dr_max_rel", dr_max_rel_);
+        dbhp_max_rel_=  param.getDefault("dbhp_max_rel", dbhp_max_rel_);
+        dwell_fraction_max_ = param.getDefault("dwell_fraction_max", dwell_fraction_max_);
         max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_);
         tolerance_mb_    = param.getDefault("tolerance_mb", tolerance_mb_);
         tolerance_cnv_   = param.getDefault("tolerance_cnv", tolerance_cnv_);
@@ -64,9 +65,11 @@ namespace Opm
     void BlackoilModelParameters::reset()
     {
         // default values for the solver parameters
-        dp_max_rel_      = 1.0e9;
+        dp_max_rel_      = 1.0;
         ds_max_          = 0.2;
         dr_max_rel_      = 1.0e9;
+        dbhp_max_rel_    = 1.0;
+        dwell_fraction_max_ = 0.2;
         max_residual_allowed_ = 1e7;
         tolerance_mb_    = 1.0e-5;
         tolerance_cnv_   = 1.0e-2;

opm/autodiff/BlackoilModelParameters.hpp

@@ -36,6 +36,10 @@ namespace Opm
         double ds_max_;
         /// Max relative change in gas-oil or oil-gas ratio in single iteration.
         double dr_max_rel_;
+        /// Max relative change in bhp in single iteration.
+        double dbhp_max_rel_;
+        /// Max absolute change in well volume fraction in single iteration.
+        double dwell_fraction_max_;
         /// Absolute max limit for residuals.
         double max_residual_allowed_;
         /// Relative mass balance tolerance (total mass balance error).

opm/autodiff/BlackoilSolventModel_impl.hpp

@@ -401,6 +401,7 @@ namespace Opm {
         const V null;
         assert(null.size() == 0);
         const V zero = V::Zero(nc);
+        const V ones = V::Constant(nc,1.0);
 
         // Extract parts of dx corresponding to each part.
         const V dp = subset(dx, Span(nc));
@@ -459,6 +460,9 @@ namespace Opm {
             }
             if (has_solvent_){
                 maxVal = dss.abs().max(maxVal);
+                // solvent is not added note that the so calculated
+                // here is overwritten later
+                //dso = dso - dss;
             }
 
             maxVal = dso.abs().max(maxVal);
@@ -488,6 +492,7 @@ namespace Opm {
             so = so_old - step * dso;
         }
 
+        // solvent is not included in the adjustment for negative saturation
         auto ixg = sg < 0;
         for (int c = 0; c < nc; ++c) {
             if (ixg[c]) {
@@ -515,10 +520,18 @@ namespace Opm {
             }
         }
 
+        auto ixs = ss < 0;
+        for (int c = 0; c < nc; ++c) {
+            if (ixs[c]) {
+                ss[c] = 0;
+            }
+        }
+
+
         // The oil saturation is defined to
         // fill the rest of the pore space.
         // For convergence reasons oil saturations
-        // is included in the appelyard chopping.
+        // is included in the appelyard chopping
         so = V::Constant(nc,1.0) - sw - sg - ss;
 
         // Update rs and rv
@@ -527,15 +540,17 @@ namespace Opm {
         if (has_disgas_) {
             const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
             const V drs = Base::isRs_ * dxvar;
-            const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel);
+            const V drs_limited = sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
             rs = rs_old - drs_limited;
+            rs = rs.max(zero);
         }
         V rv;
         if (has_vapoil_) {
             const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
             const V drv = Base::isRv_ * dxvar;
-            const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel);
+            const V drv_limited = sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
             rv = rv_old - drv_limited;
+            rv = rv.max(zero);
         }
 
         // Sg is used as primal variable for water only cells.
@@ -560,11 +575,18 @@ namespace Opm {
             for (int c = 0; c < nc; ++c) {
                 if (useSg[c]) {
                     rs[c] = rsSat[c];
+                    if (watOnly[c]) {
+                        so[c] = 0;
+                        sg[c] = 0;
+                        ss[c] = 0;
+                        rs[c] = 0;
+                    }
+
                 } else {
                     hydroCarbonState[c] = HydroCarbonState::OilOnly;
                 }
             }
-
+            rs = rs.min(rsSat);
         }
 
         // phase transitions so <-> rv
@@ -587,10 +609,19 @@ namespace Opm {
             for (int c = 0; c < nc; ++c) {
                 if (useSg[c]) {
                     rv[c] = rvSat[c];
+                    if (watOnly[c]) {
+                        so[c] = 0;
+                        sg[c] = 0;
+                        ss[c] = 0;
+                        rv[c] = 0;
+                    }
+
                 } else {
                     hydroCarbonState[c] = HydroCarbonState::GasOnly;
                 }
             }
+            rv = rv.min(rvSat);
+
         }
 
         // Update the reservoir_state
@@ -623,7 +654,7 @@ namespace Opm {
             std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
         }
 
-        wellModel().updateWellState(dwells, dpMaxRel(), well_state);
+        wellModel().updateWellState(dwells, Base::dbhpMaxRel(), well_state);
 
         for( auto w = 0; w < wells().number_of_wells; ++w ) {
             if (wells().type[w] == INJECTOR) {

opm/autodiff/StandardWellsDense.hpp

@@ -942,8 +942,8 @@ enum WellVariablePositions {
                     const int np = wells().number_of_phases;
                     const int nw = wells().number_of_wells;
 
-                    double dFLimit = 0.2;
+                    double dFLimit = dWellFractionMax();
-                    double dBHPLimit = 2.0;
+                    double dBHPLimit = dbhpMaxRel();
                     std::vector<double> xvar_well_old = well_state.wellSolutions();
 
                     for (int w = 0; w < nw; ++w) {
@@ -1510,6 +1510,9 @@ enum WellVariablePositions {
             mutable BVector invDrw_;
             mutable BVector scaleAddRes_;
 
+            double dbhpMaxRel() const {return param_.dbhp_max_rel_; }
+            double dWellFractionMax() const {return param_.dwell_fraction_max_; }
+
             // protected methods
             EvalWell getBhp(const int wellIdx) const {
                 const WellControls* wc = wells().ctrls[wellIdx];