relaxed tolerance for CNV if only small fraction of PV is violated

Previously we used relaxed tolerance once a certain number of Newton
steps was exceeded. Now we check for all cells violating CNV locally
and if their pore volume is less than a certaun fraction (default 3%)
we use the relaxed tolerance (default: 1e9)

Original idea originated from Norce.

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -677,6 +677,40 @@ namespace Opm {
             return pvSumLocal;
         }
 
+        double computeCnvErrorPv(const std::vector<Scalar>& B_avg, double dt)
+        {
+            double errorPV{};
+            const auto& ebosModel = ebosSimulator_.model();
+            const auto& ebosProblem = ebosSimulator_.problem();
+            const auto& ebosResid = ebosSimulator_.model().linearizer().residual();
+            const auto& gridView = ebosSimulator().gridView();
+            ElementContext elemCtx(ebosSimulator_);
+
+            for (const auto& elem: elements(gridView, Dune::Partitions::interiorBorder))
+            {
+                elemCtx.updatePrimaryStencil(elem);
+                elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+                const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+                const double pvValue = ebosProblem.referencePorosity(cell_idx, /*timeIdx=*/0) * ebosModel.dofTotalVolume( cell_idx );
+                const auto cellResidual = ebosResid[cell_idx];
+                bool cnvViolated = false;
+
+                for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx)
+                {
+                    using std::abs;
+                    Scalar CNV = cellResidual[eqIdx] * dt * B_avg[eqIdx] / pvValue;
+                    cnvViolated = cnvViolated || (abs(CNV) > param_.tolerance_cnv_);
+                }
+
+                if (cnvViolated)
+                {
+                    errorPV += pvValue;
+                }
+            }
+
+            return grid_.comm().sum(errorPV);
+        }
+
         ConvergenceReport getReservoirConvergence(const double dt,
                                                   const int iteration,
                                                   std::vector<Scalar>& B_avg,
@@ -684,9 +718,6 @@ namespace Opm {
         {
             typedef std::vector< Scalar > Vector;
 
-            const double tol_mb  = param_.tolerance_mb_;
-            const double tol_cnv = (iteration < param_.max_strict_iter_) ? param_.tolerance_cnv_ : param_.tolerance_cnv_relaxed_;
-
             const int numComp = numEq;
             Vector R_sum(numComp, 0.0 );
             Vector maxCoeff(numComp, std::numeric_limits< Scalar >::lowest() );
@@ -696,6 +727,12 @@ namespace Opm {
             const double pvSum = convergenceReduction(grid_.comm(), pvSumLocal,
                                                       R_sum, maxCoeff, B_avg);
 
+            auto cnvErrorPvFraction = computeCnvErrorPv(B_avg, dt);
+            cnvErrorPvFraction /= pvSum;
+
+            const double tol_mb  = param_.tolerance_mb_;
+            const double tol_cnv = (cnvErrorPvFraction < param_.relaxed_max_pv_fraction_) ? param_.tolerance_cnv_relaxed_ :  param_.tolerance_cnv_;
+
             // Finish computation
             std::vector<Scalar> CNV(numComp);
             std::vector<Scalar> mass_balance_residual(numComp);

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -32,6 +32,7 @@ NEW_TYPE_TAG(FlowModelParameters);
 NEW_PROP_TAG(DbhpMaxRel);
 NEW_PROP_TAG(DwellFractionMax);
 NEW_PROP_TAG(MaxResidualAllowed);
+NEW_PROP_TAG(RelaxedMaxPvFraction);
 NEW_PROP_TAG(ToleranceMb);
 NEW_PROP_TAG(ToleranceCnv);
 NEW_PROP_TAG(ToleranceCnvRelaxed);
@@ -62,6 +63,7 @@ NEW_PROP_TAG(MaxInnerIterWells);
 SET_SCALAR_PROP(FlowModelParameters, DbhpMaxRel, 1.0);
 SET_SCALAR_PROP(FlowModelParameters, DwellFractionMax, 0.2);
 SET_SCALAR_PROP(FlowModelParameters, MaxResidualAllowed, 1e7);
+SET_SCALAR_PROP(FlowModelParameters, RelaxedMaxPvFraction, .03);
 SET_SCALAR_PROP(FlowModelParameters, ToleranceMb, 1e-6);
 SET_SCALAR_PROP(FlowModelParameters, ToleranceCnv,1e-2);
 SET_SCALAR_PROP(FlowModelParameters, ToleranceCnvRelaxed, 1e9);
@@ -111,6 +113,9 @@ namespace Opm
         double dwell_fraction_max_;
         /// Absolute max limit for residuals.
         double max_residual_allowed_;
+        //// Max allowed pore volume faction where CNV is violated. Below the
+        //// relaxed tolerance tolerance_cnv_relaxed_ is used.
+        double relaxed_max_pv_fraction_;
         /// Relative mass balance tolerance (total mass balance error).
         double tolerance_mb_;
         /// Local convergence tolerance (max of local saturation errors).
@@ -189,6 +194,7 @@ namespace Opm
             dbhp_max_rel_=  EWOMS_GET_PARAM(TypeTag, Scalar, DbhpMaxRel);
             dwell_fraction_max_ = EWOMS_GET_PARAM(TypeTag, Scalar, DwellFractionMax);
             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_cnv_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnv);
             tolerance_cnv_relaxed_ = EWOMS_GET_PARAM(TypeTag, Scalar, ToleranceCnvRelaxed);
@@ -221,6 +227,8 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, DbhpMaxRel, "Maximum relative change of the bottom-hole pressure in a single iteration");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, DwellFractionMax, "Maximum absolute change of a well's volume fraction in a single iteration");
             EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxResidualAllowed, "Absolute maximum tolerated for residuals without cutting the time step size");
+            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, 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");