Merge pull request #5173 from totto82/relaxed_mb

Add relaxed mb tolerance parameter

ebos/ebos.hh

@@ -199,6 +199,7 @@ public:
         EWOMS_HIDE_PARAM(TypeTag, DwellFractionMax);
         EWOMS_HIDE_PARAM(TypeTag, MaxResidualAllowed);
         EWOMS_HIDE_PARAM(TypeTag, ToleranceMb);
+        EWOMS_HIDE_PARAM(TypeTag, ToleranceMbRelaxed);
         EWOMS_HIDE_PARAM(TypeTag, ToleranceCnv);
         EWOMS_HIDE_PARAM(TypeTag, ToleranceCnvRelaxed);
         EWOMS_HIDE_PARAM(TypeTag, ToleranceWells);
@@ -212,6 +213,7 @@ public:
         EWOMS_HIDE_PARAM(TypeTag, MaxInnerIterWells);
         EWOMS_HIDE_PARAM(TypeTag, MaxSinglePrecisionDays);
         EWOMS_HIDE_PARAM(TypeTag, MinStrictCnvIter);
+        EWOMS_HIDE_PARAM(TypeTag, MinStrictMbIter);
         EWOMS_HIDE_PARAM(TypeTag, SolveWelleqInitially);
         EWOMS_HIDE_PARAM(TypeTag, UpdateEquationsScaling);
         EWOMS_HIDE_PARAM(TypeTag, UseUpdateStabilization);

opm/simulators/flow/BlackoilModel.hpp

@@ -335,6 +335,7 @@ namespace Opm {
         void initialLinearization(SimulatorReportSingle& report,
                                   const int iteration,
                                   const int minIter,
+                                  const int maxIter,
                                   const SimulatorTimerInterface& timer)
         {
             // -----------   Set up reports and timer   -----------
@@ -361,7 +362,7 @@ namespace Opm {
             perfTimer.start();
             // the step is not considered converged until at least minIter iterations is done
             {
-                auto convrep = getConvergence(timer, iteration, residual_norms);
+                auto convrep = getConvergence(timer, iteration, maxIter, residual_norms);
                 report.converged = convrep.converged() && iteration >= minIter;
                 ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
                 convergence_reports_.back().report.push_back(std::move(convrep));
@@ -426,7 +427,7 @@ namespace Opm {
             SimulatorReportSingle report;
             Dune::Timer perfTimer;
 
-            this->initialLinearization(report, iteration, nonlinear_solver.minIter(), timer);
+            this->initialLinearization(report, iteration, nonlinear_solver.minIter(), nonlinear_solver.maxIter(), timer);
 
             // -----------   If not converged, solve linear system and do Newton update  -----------
             if (!report.converged) {
@@ -886,6 +887,7 @@ namespace Opm {
         ConvergenceReport getReservoirConvergence(const double reportTime,
                                                   const double dt,
                                                   const int iteration,
+                                                  const int maxIter,
                                                   std::vector<Scalar>& B_avg,
                                                   std::vector<Scalar>& residual_norms)
         {
@@ -907,12 +909,38 @@ namespace Opm {
             auto cnvErrorPvFraction = computeCnvErrorPv(B_avg, dt);
             cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
 
-            const double tol_mb  = param_.tolerance_mb_;
+
-            // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0.
+            // For each iteration, we need to determine whether to use the relaxed tolerances.
-            // For each iteration, we need to determine whether to use the relaxed CNV tolerance.
+            // To disable the usage of relaxed tolerances, you can set the relaxed tolerances as the strict tolerances.
-            // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0.
+
-            const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ && iteration >= param_.min_strict_cnv_iter_;
+            // If min_strict_mb_iter = -1 (default) we use a relaxed tolerance for the mass balance for the last iterations
-            const double tol_cnv = use_relaxed ? param_.tolerance_cnv_relaxed_ :  param_.tolerance_cnv_;
+            // For positive values we use the relaxed tolerance after the given number of iterations
+            const bool relax_final_iteration_mb = (param_.min_strict_mb_iter_ < 0) && (iteration == maxIter);
+            const bool use_relaxed_mb = relax_final_iteration_mb || (param_.min_strict_mb_iter_ >= 0 && iteration >= param_.min_strict_mb_iter_);
+
+
+            // If min_strict_cnv_iter = -1 we use a relaxed tolerance for the cnv for the last iterations
+            // For positive values we use the relaxed tolerance after the given number of iterations
+            // We also use relaxed tolerances for cells with total poro volume less than relaxed_max_pv_fraction_
+            // Default value of relaxed_max_pv_fraction_ is 0.03
+            const bool relax_final_iteration_cnv = (param_.min_strict_cnv_iter_ < 0) && (iteration == maxIter);
+            const bool relax_iter_cnv = param_.min_strict_mb_iter_ >= 0 && iteration >= param_.min_strict_mb_iter_;
+            const bool relax_pv_fraction_cnv = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_;
+            const bool use_relaxed_cnv = relax_final_iteration_cnv || relax_pv_fraction_cnv || relax_iter_cnv;
+
+            if  (relax_final_iteration_mb || relax_final_iteration_cnv)  {
+                if ( terminal_output_ ) {
+                    std::string message = "Number of newton iterations reached its maximum try to continue with relaxed tolerances:";
+                    if (relax_final_iteration_mb)
+                        message += fmt::format("  MB: {:.1e}", param_.tolerance_mb_relaxed_);
+                    if (relax_final_iteration_cnv)
+                        message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_);
+
+                    OpmLog::debug(message);
+                }
+            }
+            const double tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ :  param_.tolerance_cnv_;
+            const double tol_mb  = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
 
             // Finish computation
             std::vector<Scalar> CNV(numComp);
@@ -995,9 +1023,11 @@ namespace Opm {
         /// residual mass balance (tol_cnv).
         /// \param[in]   timer       simulation timer
         /// \param[in]   iteration   current iteration number
+        /// \param[in]   maxIter     maximum number of iterations
         /// \param[out]  residual_norms   CNV residuals by phase
         ConvergenceReport getConvergence(const SimulatorTimerInterface& timer,
                                          const int iteration,
+                                         const int maxIter,
                                          std::vector<double>& residual_norms)
         {
             OPM_TIMEBLOCK(getConvergence);
@@ -1005,7 +1035,7 @@ namespace Opm {
             std::vector<Scalar> B_avg(numEq, 0.0);
             auto report = getReservoirConvergence(timer.simulationTimeElapsed(),
                                                   timer.currentStepLength(),
-                                                  iteration, B_avg, residual_norms);
+                                                  iteration, maxIter, B_avg, residual_norms);
             {
                 OPM_TIMEBLOCK(getWellConvergence);
                 report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged());

opm/simulators/flow/BlackoilModelNldd.hpp

@@ -193,7 +193,7 @@ public:
         SimulatorReportSingle report;
         Dune::Timer perfTimer;
 
-        model_.initialLinearization(report, iteration, nonlinear_solver.minIter(), timer);
+        model_.initialLinearization(report, iteration, nonlinear_solver.minIter(), nonlinear_solver.maxIter(), timer);
 
         if (report.converged) {
             return report;
@@ -610,19 +610,21 @@ private:
         auto cnvErrorPvFraction = computeCnvErrorPvLocal(domain, B_avg, dt);
         cnvErrorPvFraction /= (pvSum - numAquiferPvSum);
 
-        const double tol_mb  = model_.param().local_tolerance_scaling_mb_ *
-                               model_.param().tolerance_mb_;
         // Default value of relaxed_max_pv_fraction_ is 0.03 and min_strict_cnv_iter_ is 0.
         // For each iteration, we need to determine whether to use the relaxed CNV tolerance.
         // To disable the usage of relaxed CNV tolerance, you can set the relaxed_max_pv_fraction_ to be 0.
-        const bool use_relaxed = cnvErrorPvFraction < model_.param().relaxed_max_pv_fraction_ &&
+        const bool use_relaxed_cnv = cnvErrorPvFraction < model_.param().relaxed_max_pv_fraction_ &&
                                  iteration >= model_.param().min_strict_cnv_iter_;
         // Tighter bound for local convergence should increase the
         // likelyhood of: local convergence => global convergence
         const double tol_cnv = model_.param().local_tolerance_scaling_cnv_
-            * (use_relaxed ? model_.param().tolerance_cnv_relaxed_
+            * (use_relaxed_cnv ? model_.param().tolerance_cnv_relaxed_
                            : model_.param().tolerance_cnv_);
 
+        const bool use_relaxed_mb = iteration >= model_.param().min_strict_mb_iter_;
+        const double tol_mb  = model_.param().local_tolerance_scaling_mb_ 
+            * (use_relaxed_mb ? model_.param().tolerance_mb_relaxed_ : model_.param().tolerance_mb_);
+
         // Finish computation
         std::vector<Scalar> CNV(numComp);
         std::vector<Scalar> mass_balance_residual(numComp);

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -62,6 +62,10 @@ struct ToleranceMb {
     using type = UndefinedProperty;
 };
 template<class TypeTag, class MyTypeTag>
+struct ToleranceMbRelaxed {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
 struct ToleranceCnv {
     using type = UndefinedProperty;
 };
@@ -94,6 +98,10 @@ struct MinStrictCnvIter {
     using type = UndefinedProperty;
 };
 template<class TypeTag, class MyTypeTag>
+struct MinStrictMbIter {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
 struct SolveWelleqInitially {
     using type = UndefinedProperty;
 };
@@ -261,6 +269,11 @@ struct ToleranceMb<TypeTag, TTag::FlowModelParameters> {
     static constexpr type value = 1e-6;
 };
 template<class TypeTag>
+struct ToleranceMbRelaxed<TypeTag, TTag::FlowModelParameters> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1e-6;
+};
+template<class TypeTag>
 struct ToleranceCnv<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
     static constexpr type value = 1e-2;
@@ -298,6 +311,10 @@ struct MinStrictCnvIter<TypeTag, TTag::FlowModelParameters> {
     static constexpr int value = 0;
 };
 template<class TypeTag>
+struct MinStrictMbIter<TypeTag, TTag::FlowModelParameters> {
+    static constexpr int value = -1;
+};
+template<class TypeTag>
 struct SolveWelleqInitially<TypeTag, TTag::FlowModelParameters> {
     static constexpr bool value = true;
 };
@@ -480,6 +497,8 @@ namespace Opm
         double relaxed_max_pv_fraction_;
         /// Relative mass balance tolerance (total mass balance error).
         double tolerance_mb_;
+        /// Relaxed mass balance tolerance (can be used when iter >= min_strict_mb_iter_).
+        double tolerance_mb_relaxed_;
         /// Local convergence tolerance (max of local saturation errors).
         double tolerance_cnv_;
         /// Relaxed local convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViolatedPV < relaxed_max_pv_fraction_).
@@ -531,6 +550,9 @@ namespace Opm
         /// Minimum number of Newton iterations before we can use relaxed CNV convergence criterion
         int min_strict_cnv_iter_;
 
+        /// Minimum number of Newton iterations before we can use relaxed MB convergence criterion
+        int min_strict_mb_iter_;
+
         /// Solve well equation initially
         bool solve_welleq_initially_;
 
@@ -602,6 +624,7 @@ namespace Opm
             max_residual_allowed_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxResidualAllowed);
             relaxed_max_pv_fraction_ = EWOMS_GET_PARAM(TypeTag, Scalar, RelaxedMaxPvFraction);
             tolerance_mb_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceMb);
+            tolerance_mb_relaxed_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceMbRelaxed);
             tolerance_cnv_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnv);
             tolerance_cnv_relaxed_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnvRelaxed);
             tolerance_wells_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceWells);
@@ -621,6 +644,7 @@ namespace Opm
             max_inner_iter_wells_ = EWOMS_GET_PARAM(TypeTag, int, MaxInnerIterWells);
             maxSinglePrecisionTimeStep_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays) *24*60*60;
             min_strict_cnv_iter_ = EWOMS_GET_PARAM(TypeTag, int, MinStrictCnvIter);
+            min_strict_mb_iter_ = EWOMS_GET_PARAM(TypeTag, int, MinStrictMbIter);
             solve_welleq_initially_ = EWOMS_GET_PARAM(TypeTag, bool, SolveWelleqInitially);
             update_equations_scaling_ = EWOMS_GET_PARAM(TypeTag, bool, UpdateEquationsScaling);
             use_update_stabilization_ = EWOMS_GET_PARAM(TypeTag, bool, UseUpdateStabilization);
@@ -663,6 +687,7 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, RelaxedMaxPvFraction, "The fraction of the pore volume of the reservoir "
                                  "where the volumetric error (CNV) may be voilated during strict Newton iterations.");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceMb, "Tolerated mass balance error relative to total mass present");
+            EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceMbRelaxed, "Relaxed tolerated mass balance error that applies for iterations after the iterations with the strict tolerance");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceCnv, "Local convergence tolerance (Maximum of local saturation errors)");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceCnvRelaxed, "Relaxed local convergence tolerance that applies for iterations after the iterations with the strict tolerance");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, ToleranceWells, "Well convergence tolerance");
@@ -683,6 +708,8 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, RegularizationFactorWells, "Regularization factor for 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, MinStrictCnvIter, "Minimum number of Newton iterations before relaxed tolerances can be used for the CNV convergence criterion");
+            EWOMS_REGISTER_PARAM(TypeTag, int, MinStrictMbIter, "Minimum number of Newton iterations before relaxed tolerances can be used for the MB convergence criterion. "
+                                "Default -1 means that the relaxed tolerance is used when maximum number of Newton iterations are reached.");
             EWOMS_REGISTER_PARAM(TypeTag, bool, SolveWelleqInitially, "Fully solve the well equations before each iteration of the reservoir model");
             EWOMS_REGISTER_PARAM(TypeTag, bool, UpdateEquationsScaling, "Update scaling factors for mass balance equations during the run");
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseUpdateStabilization, "Try to detect and correct oscillations or stagnation during the Newton method");

regressionTests.cmake

@@ -299,7 +299,7 @@ add_test_compareECLFiles(CASENAME polymer_simple2D
                          SIMULATOR flow
                          ABS_TOL ${abs_tol}
                          REL_TOL ${coarse_rel_tol}
-                         TEST_ARGS --tolerance-mb=1.e-7)
+                         TEST_ARGS --tolerance-mb=1.e-7 --tolerance-mb-relaxed=1.e-7)
 
 add_test_compareECLFiles(CASENAME spe5
                          FILENAME SPE5CASE1