Merge pull request #2611 from totto82/improveConvergenceMSW

Improve convergence msw

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -52,6 +52,10 @@ NEW_PROP_TAG(TolerancePressureMsWells);
 NEW_PROP_TAG(MaxPressureChangeMsWells);
 NEW_PROP_TAG(UseInnerIterationsMsWells);
 NEW_PROP_TAG(MaxInnerIterMsWells);
+NEW_PROP_TAG(StrictInnerIterMsWells);
+NEW_PROP_TAG(RelaxedFlowTolInnerIterMsw);
+NEW_PROP_TAG(RelaxedPressureTolInnerIterMsw);
+NEW_PROP_TAG(RegularizationFactorMsw);
 
 SET_SCALAR_PROP(FlowModelParameters, DbhpMaxRel, 1.0);
 SET_SCALAR_PROP(FlowModelParameters, DwellFractionMax, 0.2);
@@ -73,8 +77,13 @@ SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01*1e5);
 SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 10*1e5);
 SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true);
 SET_INT_PROP(FlowModelParameters, MaxInnerIterMsWells, 100);
+SET_INT_PROP(FlowModelParameters, StrictInnerIterMsWells, 40);
+SET_SCALAR_PROP(FlowModelParameters, RegularizationFactorMsw, 1);
 SET_BOOL_PROP(FlowModelParameters, EnableWellOperabilityCheck, true);
 
+SET_SCALAR_PROP(FlowModelParameters, RelaxedFlowTolInnerIterMsw, 1);
+SET_SCALAR_PROP(FlowModelParameters, RelaxedPressureTolInnerIterMsw, 0.5e5);
+
 // if openMP is available, determine the number threads per process automatically.
 #if _OPENMP
 SET_INT_PROP(FlowModelParameters, ThreadsPerProcess, -1);
@@ -111,6 +120,12 @@ namespace Opm
         double tolerance_well_control_;
         /// Tolerance for the pressure equations for multisegment wells
         double tolerance_pressure_ms_wells_;
+
+        /// Relaxed tolerance for the inner iteration for the MSW flow solution
+        double relaxed_inner_tolerance_flow_ms_well_;
+        /// Relaxed tolerance for the inner iteration for the MSW pressure solution
+        double relaxed_inner_tolerance_pressure_ms_well_;
+
         /// Maximum pressure change over an iteratio for ms wells
         double max_pressure_change_ms_wells_;
 
@@ -120,6 +135,12 @@ namespace Opm
         /// Maximum inner iteration number for ms wells
         int max_inner_iter_ms_wells_;
 
+        /// Strict inner iteration number for ms wells
+        int strict_inner_iter_ms_wells_;
+
+        /// Regularization factor for ms wells
+        int regularization_factor_ms_wells_;
+
         /// Maximum iteration number of the well equation solution
         int max_welleq_iter_;
 
@@ -166,9 +187,13 @@ namespace Opm
             max_welleq_iter_ = EWOMS_GET_PARAM(TypeTag, int, MaxWelleqIter);
             use_multisegment_well_ = EWOMS_GET_PARAM(TypeTag, bool, UseMultisegmentWell);
             tolerance_pressure_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, TolerancePressureMsWells);
+            relaxed_inner_tolerance_flow_ms_well_ = EWOMS_GET_PARAM(TypeTag, Scalar, RelaxedFlowTolInnerIterMsw);
+            relaxed_inner_tolerance_pressure_ms_well_ = EWOMS_GET_PARAM(TypeTag, Scalar, RelaxedPressureTolInnerIterMsw);
             max_pressure_change_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxPressureChangeMsWells);
             use_inner_iterations_ms_wells_ = EWOMS_GET_PARAM(TypeTag, bool, UseInnerIterationsMsWells);
             max_inner_iter_ms_wells_ = EWOMS_GET_PARAM(TypeTag, int, MaxInnerIterMsWells);
+            strict_inner_iter_ms_wells_ = EWOMS_GET_PARAM(TypeTag, int, StrictInnerIterMsWells);
+            regularization_factor_ms_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, RegularizationFactorMsw);
             maxSinglePrecisionTimeStep_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays) *24*60*60;
             max_strict_iter_ = EWOMS_GET_PARAM(TypeTag, int, MaxStrictIter);
             solve_welleq_initially_ = EWOMS_GET_PARAM(TypeTag, bool, SolveWelleqInitially);
@@ -192,9 +217,13 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, int, MaxWelleqIter, "Maximum number of iterations to determine solution the  well equations");
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseMultisegmentWell, "Use the well model for multi-segment wells instead of the one for single-segment wells");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, TolerancePressureMsWells, "Tolerance for the pressure equations for multi-segment wells");
+            EWOMS_REGISTER_PARAM(TypeTag, Scalar, RelaxedFlowTolInnerIterMsw, "Relaxed tolerance for the inner iteration for the MSW flow solution");
+            EWOMS_REGISTER_PARAM(TypeTag, Scalar, RelaxedPressureTolInnerIterMsw, "Relaxed tolerance for the inner iteration for the MSW pressure solution");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxPressureChangeMsWells, "Maximum relative pressure change for a single iteration of the multi-segment well model");
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseInnerIterationsMsWells, "Use nested iterations for multi-segment wells");
             EWOMS_REGISTER_PARAM(TypeTag, int, MaxInnerIterMsWells, "Maximum number of inner iterations for multi-segment wells");
+            EWOMS_REGISTER_PARAM(TypeTag, int, StrictInnerIterMsWells, "Number of inner iterations for multi-segment wells with strict tolerance");
+            EWOMS_REGISTER_PARAM(TypeTag, Scalar, RegularizationFactorMsw, "Regularization factor for ms wells");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays, "Maximum time step size where single precision floating point arithmetic can be used solving for the linear systems of equations");
             EWOMS_REGISTER_PARAM(TypeTag, int, MaxStrictIter, "Maximum number of Newton iterations before relaxed tolerances are used for the CNV convergence criterion");
             EWOMS_REGISTER_PARAM(TypeTag, bool, SolveWelleqInitially, "Fully solve the well equations before each iteration of the reservoir model");

opm/simulators/wells/MultisegmentWell.hpp

@@ -128,7 +128,7 @@ namespace Opm
                                                Opm::DeferredLogger& deferred_logger) const override;
 
         /// check whether the well equations get converged for this well
-        virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const override;
+        virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger, const bool relax_tolerance = false) const override;
 
         /// Ax = Ax - C D^-1 B x
         virtual void apply(const BVector& x, BVector& Ax) const override;

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -536,7 +536,7 @@ namespace Opm
     template <typename TypeTag>
     ConvergenceReport
     MultisegmentWell<TypeTag>::
-    getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const
+    getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger, const bool relax_tolerance) const
     {
         assert(int(B_avg.size()) == num_components_);
 
@@ -577,11 +577,14 @@ namespace Opm
             if (eq_idx < num_components_) { // phase or component mass equations
                 const double flux_residual = maximum_residual[eq_idx];
                 // TODO: the report can not handle the segment number yet.
+
                 if (std::isnan(flux_residual)) {
                     report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::NotANumber, eq_idx, name()});
                 } else if (flux_residual > param_.max_residual_allowed_) {
                     report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::TooLarge, eq_idx, name()});
-                } else if (flux_residual > param_.tolerance_wells_) {
+                } else if (!relax_tolerance && flux_residual > param_.tolerance_wells_) {
+                    report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()});
+                } else if (flux_residual > param_.relaxed_inner_tolerance_flow_ms_well_) {
                     report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, name()});
                 }
             } else { // pressure equation
@@ -591,7 +594,9 @@ namespace Opm
                     report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, name()});
                 } else if (std::isinf(pressure_residual)) {
                     report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::TooLarge, dummy_component, name()});
-                } else if (pressure_residual > param_.tolerance_pressure_ms_wells_) {
+                } else if (!relax_tolerance && pressure_residual > param_.tolerance_pressure_ms_wells_) {
+                    report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()});
+                } else if (pressure_residual > param_.relaxed_inner_tolerance_pressure_ms_well_) {
                     report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::Normal, dummy_component, name()});
                 }
             }
@@ -1119,8 +1124,8 @@ namespace Opm
             // update the segment pressure
             {
                 const int sign = dwells[seg][SPres] > 0.? 1 : -1;
-                const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change);
+                //const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change);
-                // const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
+                const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
                 primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5);
             }
 
@@ -2369,17 +2374,20 @@ namespace Opm
         int it = 0;
         // relaxation factor
         double relaxation_factor = 1.;
-        const double min_relaxation_factor = 0.2;
+        const double min_relaxation_factor = 0.6;
         bool converged = false;
         int stagnate_count = 0;
+        bool relax_convergence = false;
         for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
 
             assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, deferred_logger);
 
             const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
 
+            if (it > param_.strict_inner_iter_ms_wells_)
+                relax_convergence = true;
 
-            const auto report = getWellConvergence(well_state, B_avg, deferred_logger);
+            const auto report = getWellConvergence(well_state, B_avg, deferred_logger, relax_convergence);
             if (report.converged()) {
                 converged = true;
                 break;
@@ -2395,19 +2403,22 @@ namespace Opm
             // TODO: maybe we should have more sophiscated strategy to recover the relaxation factor,
             // for example, to recover it to be bigger
 
-            if (!is_stagnate) {
-                stagnate_count = 0;
-            }
             if (is_oscillate || is_stagnate) {
                 // HACK!
-                if (is_stagnate && relaxation_factor == min_relaxation_factor) {
+                std::ostringstream sstr;
+                if (relaxation_factor == min_relaxation_factor) {
                     // Still stagnating, terminate iterations if 5 iterations pass.
                     ++stagnate_count;
-                    if (stagnate_count == 5) {
+                    if (stagnate_count == 6) {
-                        // break;
+                        sstr << " well " << name() << " observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
+                        const auto reportStag = getWellConvergence(well_state, B_avg, deferred_logger, true);
+                        if (reportStag.converged()) {
+                            converged = true;
+                            sstr << " well " << name() << " manages to get converged with relaxed tolerances in " << it << " inner iterations";
+                            deferred_logger.debug(sstr.str());
+                            return;
+                        }
                     }
-                } else {
-                    stagnate_count = 0;
                 }
 
                 // a factor value to reduce the relaxation_factor
@@ -2415,7 +2426,6 @@ namespace Opm
                 relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
 
                 // debug output
-                std::ostringstream sstr;
                 if (is_stagnate) {
                     sstr << " well " << name() << " observes stagnation in inner iteration " << it << "\n";
 
@@ -2433,7 +2443,9 @@ namespace Opm
         // TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state
         if (converged) {
             std::ostringstream sstr;
-            sstr << " well " << name() << " manage to get converged within " << it << " inner iterations";
+            sstr << " well " << name() << " manage to get converged within " << it << " inner iterations \n";
+            if (relax_convergence)
+                sstr << "A relaxed tolerance is used after "<< param_.strict_inner_iter_ms_wells_ << " iterations";
             deferred_logger.debug(sstr.str());
         } else {
             std::ostringstream sstr;
@@ -2495,9 +2507,14 @@ namespace Opm
             // TODO: without considering the efficiencty factor for now
             {
                 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
+
+                // Add a regularization_factor to increase the accumulation term
+                // This will make the system less stiff and help convergence for 
+                // difficult cases
+                const Scalar regularization_factor =  param_.regularization_factor_ms_wells_;
                 // for each component
                 for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
-                    const EvalWell accumulation_term = (segment_surface_volume * surfaceVolumeFraction(seg, comp_idx)
+                    const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * surfaceVolumeFraction(seg, comp_idx)
                                                      - segment_fluid_initial_[seg][comp_idx]) / dt;
 
                     resWell_[seg][comp_idx] += accumulation_term.value();
@@ -2818,7 +2835,7 @@ namespace Opm
             vol_ratio += mix[comp_idx] / b[comp_idx];
         }
 
-        // segment volume
+        // We increase the segment volume with a factor 10 to stabilize the system.
         const double volume = segmentSet()[seg_idx].volume();
 
         return volume / vol_ratio;

opm/simulators/wells/StandardWell.hpp

@@ -177,7 +177,8 @@ namespace Opm
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const WellState& well_state,
                                                      const std::vector<double>& B_avg,
-                                                     Opm::DeferredLogger& deferred_logger) const override;
+                                                     Opm::DeferredLogger& deferred_logger,
+                                                     const bool relax_tolerance = false) const override;
 
         /// Ax = Ax - C D^-1 B x
         virtual void apply(const BVector& x, BVector& Ax) const override;

opm/simulators/wells/StandardWell_impl.hpp

@@ -2087,7 +2087,8 @@ namespace Opm
     StandardWell<TypeTag>::
     getWellConvergence(const WellState& well_state,
                        const std::vector<double>& B_avg,
-                       Opm::DeferredLogger& deferred_logger) const
+                       Opm::DeferredLogger& deferred_logger,
+                       const bool /*relax_tolerance*/) const
     {
         // the following implementation assume that the polymer is always after the w-o-g phases
         // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells

opm/simulators/wells/WellInterface.hpp

@@ -157,7 +157,7 @@ namespace Opm
 
         virtual void initPrimaryVariablesEvaluation() const = 0;
 
-        virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const = 0;
+        virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger, const bool relax_tolerance = false) const = 0;
 
         virtual void solveEqAndUpdateWellState(WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0;