Improvments in convergence for flow_ebos

- restrict pressure changes. Set default to 1.0 (this also effects flow)
- change default number of linear iterations to 150
- tell stabilized newton the residual occilates even if it occilates in
only one phase (this also effects flow)
- avoid problems realated to division on small numbers

Tested on SPE9, norne and Model 2 with significant improvments.

opm/autodiff/BlackoilModelEbos.hpp

@@ -532,11 +532,13 @@ namespace Opm {
             const int nc = numCells(grid_);
 
             for (int cell_idx = 0; cell_idx < nc; ++cell_idx) {
-                double dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)];
+                const double& dp = dx[cell_idx][flowPhaseToEbosCompIdx(0)];
                 //reservoir_state.pressure()[cell_idx] -= dp;
                 double& p = reservoir_state.pressure()[cell_idx];
-                p -= dp;
+                const double& dp_rel_max = dpMaxRel();
-                p = std::max(p, 1e5);
+                const int sign_dp = dp > 0 ? 1: -1;
+                p -= sign_dp * std::min(std::abs(dp), std::abs(p)*dp_rel_max);
+                p = std::max(p, 0.0);
 
                 // Saturation updates.
                 const double dsw = active_[Water] ? dx[cell_idx][flowPhaseToEbosCompIdx(1)] : 0.0;
@@ -623,7 +625,7 @@ namespace Opm {
                         reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::OilOnly; // sg --> rs
                         sg = 0;
                         so = 1.0 - sw - sg;
-                        double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
+                        const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
                         double& rs = reservoir_state.gasoilratio()[cell_idx];
                         rs = rsSat*(1-epsilon);
                     } else if (so <= 0.0 && has_vapoil_) {
@@ -632,7 +634,7 @@ namespace Opm {
                         sg = 1.0 - sw - so;
                         double& rv = reservoir_state.rv()[cell_idx];
                         // use gas pressure?
-                        double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
+                        const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
                         rv = rvSat*(1-epsilon);
                     }
                     break;
@@ -650,7 +652,7 @@ namespace Opm {
 
 
 
-                    double rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
+                    const double& rsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
                     if (rs > ( rsSat * (1+epsilon) ) ) {
                         reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
                         sg = epsilon;
@@ -668,7 +670,7 @@ namespace Opm {
                         //std::cout << "watonly rv -> sg" << cell_idx << std::endl;
                         break;
                     }
-                    double rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
+                    const double& rvSat = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), reservoir_state.temperature()[cell_idx], reservoir_state.pressure()[cell_idx]);
                     if (rv > rvSat * (1+epsilon) ) {
                         reservoir_state.hydroCarbonState()[cell_idx] = HydroCarbonState::GasAndOil;
                         so = epsilon;

opm/autodiff/BlackoilModelParameters.cpp

@@ -59,7 +59,7 @@ 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;
         max_residual_allowed_ = 1e7;

opm/autodiff/NewtonIterationBlackoilInterleaved.hpp

@@ -66,7 +66,7 @@ namespace Opm
         {
             newton_use_gmres_        = false;
             linear_solver_reduction_ = 1e-2;
-            linear_solver_maxiter_   = 75;
+            linear_solver_maxiter_   = 150;
             linear_solver_restart_   = 40;
             linear_solver_verbosity_ = 0;
             require_full_sparsity_pattern_ = false;

opm/autodiff/NonlinearSolver_impl.hpp

@@ -248,7 +248,7 @@ namespace Opm
             stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3));
         }
 
-        oscillate = (oscillatePhase > 1);
+        oscillate = (oscillatePhase > 0);
     }
 
 

opm/autodiff/StandardWellsDense.hpp

@@ -1502,7 +1502,8 @@ namespace Opm {
                         currentControlIdx += wells().comp_frac[np*wellIdx + i] * i;
                     }
 
-                    if (wellVolumeFractionScaled(wellIdx,currentControlIdx) == 0) {
+                    const double eps = 1e-6;
+                    if (wellVolumeFractionScaled(wellIdx,currentControlIdx) < eps) {
                         return qs;
                     }
                     return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx,currentControlIdx));