Add an optional pore-volume floor in the CNV normalisation

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -43,6 +43,7 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     inj_mult_min_damp_factor_ = Parameters::Get<Parameters::InjMultMinDampFactor<Scalar>>();
     max_residual_allowed_ = Parameters::Get<Parameters::MaxResidualAllowed<Scalar>>();
     relaxed_max_pv_fraction_ = Parameters::Get<Parameters::RelaxedMaxPvFraction<Scalar>>();
+    cnv_pv_floor_fraction_ = Parameters::Get<Parameters::CnvPoreVolumeFloorFraction<Scalar>>();
     tolerance_mb_ = Parameters::Get<Parameters::ToleranceMb<Scalar>>();
     tolerance_mb_relaxed_ = std::max(tolerance_mb_, Parameters::Get<Parameters::ToleranceMbRelaxed<Scalar>>());
     tolerance_energy_balance_ = Parameters::Get<Parameters::ToleranceEnergyBalance<Scalar>>();
@@ -143,6 +144,12 @@ void BlackoilModelParameters<Scalar>::registerParameters()
         ("Minimum injection multiplier dampening factor (maximum dampening level)");
     Parameters::Register<Parameters::MaxResidualAllowed<Scalar>>
         ("Absolute maximum tolerated for residuals without cutting the time step size");
+    Parameters::Register<Parameters::CnvPoreVolumeFloorFraction<Scalar>>
+        ("Floor the per-cell pore volume used in the CNV normalisation at this fraction "
+         "of the average cell pore volume. The CNV criterion divides the residual by the "
+         "cell pore volume, so cells with vanishing PV (inner radial cells, salt-clogged "
+         "cells, degenerate geometry) are held to a divergently strict standard; the "
+         "floor bounds the criterion for such cells. 0 (default) disables.");
     Parameters::Register<Parameters::RelaxedMaxPvFraction<Scalar>>
         ("The fraction of the pore volume of the reservoir "
          "where the volumetric error (CNV) may be violated "

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -49,6 +49,9 @@ struct MaxResidualAllowed { static constexpr Scalar value = 1e7; };
 template<class Scalar>
 struct RelaxedMaxPvFraction { static constexpr Scalar value = 0.03; };
 
+template<class Scalar>
+struct CnvPoreVolumeFloorFraction { static constexpr Scalar value = 0.0; };
+
 template<class Scalar>
 struct ToleranceMb { static constexpr Scalar value = 1e-7; };
 
@@ -221,6 +224,11 @@ struct BlackoilModelParameters
     //// Max allowed pore volume faction where CNV is violated. Below the
     //// relaxed tolerance tolerance_cnv_relaxed_ is used.
     Scalar relaxed_max_pv_fraction_;
+
+    /// Floor (fraction of average cell PV) for the per-cell PV in the CNV
+    /// normalisation; 0 disables. Vanishing-PV cells otherwise make CNV
+    /// divergently strict.
+    Scalar cnv_pv_floor_fraction_;
     /// Relative mass balance tolerance (total mass balance error).
     Scalar tolerance_mb_;
     /// Relaxed mass balance tolerance (can be used when iter >= min_strict_mb_iter_).

opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp

@@ -590,6 +590,19 @@ localConvergenceData(std::vector<Scalar>& R_sum,
 
     const auto& residual = this->simulator_.model().linearizer().residual();
 
+    // Absolute PV floor for the CNV normalisation (0 disables); the average cell pore
+    // volume is cheap to form from static data.
+    Scalar cnv_pv_floor = 0.0;
+    if (this->param_.cnv_pv_floor_fraction_ > 0.0) {
+        Scalar pvTot = 0.0;
+        const unsigned n = model.numGridDof();
+        for (unsigned dof = 0; dof < n; ++dof) {
+            pvTot += problem.referencePorosity(dof, /*timeIdx=*/0) * model.dofTotalVolume(dof);
+        }
+        const Scalar avg = this->grid_.comm().sum(pvTot) / Scalar(this->global_nc_);
+        cnv_pv_floor = this->param_.cnv_pv_floor_fraction_ * avg;
+    }
+
     ElementContext elemCtx(this->simulator_);
     const auto& gridView = this->simulator().gridView();
     IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
@@ -610,7 +623,14 @@ localConvergenceData(std::vector<Scalar>& R_sum,
             numAquiferPvSumLocal += pvValue;
         }
 
-        this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvValue,
+        // CNV divides the residual by the cell pore volume, so a cell with vanishing PV
+        // (inner radial cells, salt-clogged cells, degenerate geometry) is held to a
+        // divergently strict standard although its capacity to accumulate mass error is
+        // limited by throughput, not storage. Optionally bound the criterion for such
+        // cells by flooring the normalisation volume.
+        const auto pvNorm = std::max(pvValue, cnv_pv_floor);
+
+        this->getMaxCoeff(cell_idx, intQuants, fs, residual, pvNorm,
                           B_avg, R_sum, maxCoeff, maxCoeffCell);
     }