From e0404845801eb9cdc001a14402a8700a7d87bccd Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 2 Jul 2020 11:52:44 +0200 Subject: [PATCH 1/4] 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 | 43 +++++++++++++++++-- .../flow/BlackoilModelParametersEbos.hpp | 8 ++++ 2 files changed, 48 insertions(+), 3 deletions(-) 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"); From f65eca331015e483f5ae8142d620c74c032d17ad Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 18 Aug 2020 16:55:23 +0200 Subject: [PATCH 2/4] Use RelaxedMaxPvFraction 0 for now. This should make the tests work. --- opm/simulators/flow/BlackoilModelParametersEbos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index a18646016..de22348b0 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -63,7 +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, RelaxedMaxPvFraction, .0); SET_SCALAR_PROP(FlowModelParameters, ToleranceMb, 1e-6); SET_SCALAR_PROP(FlowModelParameters, ToleranceCnv,1e-2); SET_SCALAR_PROP(FlowModelParameters, ToleranceCnvRelaxed, 1e9); From efe17780e2b8e9db69985d39e1a52d2c9d74ab87 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 19 Aug 2020 08:50:39 +0200 Subject: [PATCH 3/4] Reinforce maxStrictIter for the newton. That was removed before in lieu of the fraction of cells that violate CNV. This change should make the results as before unless somebody changes maxStrictIter or RelaxedMaxPvFraction --- opm/simulators/flow/BlackoilModelEbos.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 45089c89a..c7600631c 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -731,7 +731,8 @@ namespace Opm { 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_; + const bool use_relaxed = cnvErrorPvFraction < param_.relaxed_max_pv_fraction_ || iteration >= param_.max_strict_iter_; + const double tol_cnv = use_relaxed ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_; // Finish computation std::vector CNV(numComp); From 92b1708433f638a53f973244b1c3f4090b3df009 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Tue, 25 Aug 2020 14:07:27 +0200 Subject: [PATCH 4/4] Use reference to ebos residual to prevent copying. --- opm/simulators/flow/BlackoilModelEbos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index c7600631c..92f03e458 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -692,7 +692,7 @@ namespace Opm { 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]; + const auto& cellResidual = ebosResid[cell_idx]; bool cnvViolated = false; for (unsigned eqIdx = 0; eqIdx < cellResidual.size(); ++eqIdx)