diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index c6f506321..45089c89a 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -677,6 +677,40 @@ namespace Opm { return pvSumLocal; } + double computeCnvErrorPv(const std::vector& 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& 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 CNV(numComp); std::vector mass_balance_residual(numComp); diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index 0e7bc8c51..a18646016 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/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");