From 3aaff5320606207c6373e6ceb6dfa341b78e4ce3 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Fri, 9 Jan 2026 14:55:44 +0100 Subject: [PATCH 01/11] TUNINGDP implementation. - If TRGDDP and/or TRGDDS is set and met -> use relaxed CNV tolerances - Runtime parameters for max dp and ds --- opm/simulators/flow/BlackoilModel.hpp | 15 +- .../flow/BlackoilModelParameters.cpp | 10 + .../flow/BlackoilModelParameters.hpp | 12 +- opm/simulators/flow/BlackoilModel_impl.hpp | 206 ++++++++++++++++-- .../flow/SimulatorFullyImplicitBlackoil.hpp | 35 +++ .../utils/UnsupportedFlowKeywords.cpp | 1 - 6 files changed, 262 insertions(+), 17 deletions(-) diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index 32fbff04a5..643b393ba0 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -75,6 +75,8 @@ public: using Scalar = GetPropType; using ModelParameters = BlackoilModelParameters; + static constexpr bool enableSaltPrecipitation = getPropValue(); + static constexpr int numEq = Indices::numEq; static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx; static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx; @@ -182,6 +184,11 @@ public: /// Apply an update to the primary variables. void updateSolution(const BVector& dx); + /// Get solution update vector as a PrimaryVarible + void prepareStoringSolutionUpdate(); + void storeSolutionUpdate(const BVector& dx); + std::pair getMaxSolutionUpdate(const std::vector& ixCells); + /// Return true if output to cout is wanted. bool terminalOutputEnabled() const { return terminal_output_; } @@ -206,7 +213,10 @@ public: /// \brief Compute pore-volume/cell count split among "converged", /// "relaxed converged", "unconverged" cells based on CNV point /// measures. - std::pair, std::vector> + std::tuple< + std::pair, std::vector>, + std::vector + > characteriseCnvPvSplit(const std::vector& B_avg, const double dt); /// \brief Compute the number of Newtons required by each cell in order to @@ -218,6 +228,7 @@ public: const int iteration); void updateTUNING(const Tuning& tuning); + void updateTUNINGDP(const TuningDp& tuning_dp); ConvergenceReport getReservoirConvergence(const double reportTime, @@ -344,6 +355,8 @@ protected: Scalar current_relaxation_; BVector dx_old_; + SolutionVector solUpd_; + std::vector convergence_reports_; ComponentName compNames_{}; diff --git a/opm/simulators/flow/BlackoilModelParameters.cpp b/opm/simulators/flow/BlackoilModelParameters.cpp index be75780752..0aa4715477 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -49,6 +49,8 @@ BlackoilModelParameters::BlackoilModelParameters() tolerance_cnv_relaxed_ = std::max(tolerance_cnv_, Parameters::Get>()); tolerance_cnv_energy_ = Parameters::Get>(); tolerance_cnv_energy_relaxed_ = std::max(tolerance_cnv_energy_, Parameters::Get>()); + tolerance_max_dp_ = Parameters::Get>(); + tolerance_max_ds_ = Parameters::Get>(); tolerance_wells_ = Parameters::Get>(); tolerance_well_control_ = Parameters::Get>(); max_welleq_iter_ = Parameters::Get(); @@ -159,6 +161,14 @@ void BlackoilModelParameters::registerParameters() Parameters::Register> ("Relaxed local energy convergence tolerance that applies for iterations " "after the iterations with the strict tolerance"); + Parameters::Register> + ("Tolerance for max pressure change during a Newton iteration. " + "A value greater than 0.0 allows for convergence according to " + "relaxed local convergence tolerance"); + Parameters::Register> + ("Tolerance for max saturation change during a Newton iteration. " + "A value greater than 0.0 allows for convergence according to " + "relaxed local convergence tolerance."); Parameters::Register> ("Well convergence tolerance"); Parameters::Register> diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index 55ed751d5e..a366130b56 100644 --- a/opm/simulators/flow/BlackoilModelParameters.hpp +++ b/opm/simulators/flow/BlackoilModelParameters.hpp @@ -56,7 +56,7 @@ template struct ToleranceMbRelaxed { static constexpr Scalar value = 1e-6; }; //TODO change to a simpler number with fewer digits -//converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water +//converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water template struct ToleranceEnergyBalance { static constexpr Scalar value = 1e-7*41.82; }; @@ -75,6 +75,12 @@ struct ToleranceCnvEnergy { static constexpr Scalar value = 1e-2*41.82; }; template struct ToleranceCnvEnergyRelaxed { static constexpr Scalar value = 1.0*41.82; }; +template +struct ToleranceMaxDp { static constexpr Scalar value = 0.0; }; + +template +struct ToleranceMaxDs { static constexpr Scalar value = 0.0; }; + template struct ToleranceWells { static constexpr Scalar value = 1e-4; }; @@ -212,6 +218,10 @@ public: Scalar tolerance_cnv_energy_; /// Relaxed local energy convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViolatedPV < relaxed_max_pv_fraction_). Scalar tolerance_cnv_energy_relaxed_; + /// Max pressure change during a Newton iteration (TUNINGDP item = TRGDDP) + Scalar tolerance_max_dp_; + /// Max saturation change during a Newton iteration (TUNINGDP item = TRGDDS) + Scalar tolerance_max_ds_; /// Well convergence tolerance. Scalar tolerance_wells_; /// Tolerance for the well control equations diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index c951c42f87..f5dd0eb3be 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -520,6 +520,11 @@ BlackoilModel:: updateSolution(const BVector& dx) { OPM_TIMEBLOCK(updateSolution); + // Prepare to store solution update for convergence check + if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0) { + prepareStoringSolutionUpdate(); + } + auto& newtonMethod = simulator_.model().newtonMethod(); SolutionVector& solution = simulator_.model().solution(/*timeIdx=*/0); @@ -535,6 +540,95 @@ updateSolution(const BVector& dx) OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities); simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); } + + // Store solution update + if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0) { + storeSolutionUpdate(dx); + } +} + +template +void +BlackoilModel:: +prepareStoringSolutionUpdate() +{ + // Init. solution update vector + unsigned nc = simulator_.model().numGridDof(); + solUpd_.resize(nc); + + const auto& elemMapper = simulator_.model().elementMapper(); + const auto& gridView = simulator_.gridView(); + for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { + // Copy solution vector to transfer primary variables meaning + unsigned globalElemIdx = elemMapper.index(elem); + solUpd_[globalElemIdx] = simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; + + // Ensure each element is zero + for (auto& elem : solUpd_[globalElemIdx]) { + elem = 0.0; + } + } +} + +template +void +BlackoilModel:: +storeSolutionUpdate(const BVector& dx) +{ + const auto& elemMapper = simulator_.model().elementMapper(); + const auto& gridView = simulator_.gridView(); + for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { + // Get cell vectors + unsigned globalElemIdx = elemMapper.index(elem); + auto& value = solUpd_[globalElemIdx]; + auto& update = dx[globalElemIdx]; + assert(value.size() == update.size()); + + // Transfer update from dx to solution update container (SolutionVector type) + for (std::size_t i = 0; i < value.size(); ++i) { + value[i] = update[i]; + } + } +} + +template +std::pair::Scalar, + typename BlackoilModel::Scalar> +BlackoilModel:: +getMaxSolutionUpdate(const std::vector& ixCells) +{ + static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0; + static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0; + + // Init output + Scalar dPMax = 0.0; + Scalar dSMax = 0.0; + + // Loop over solution update, get the correct variables and calculate max. + for (const auto& ix : ixCells) { + const auto& value = solUpd_[ix]; + for (std::size_t pvIdx = 0; pvIdx < value.size(); ++pvIdx) { + if (pvIdx == Indices::pressureSwitchIdx) { + dPMax = std::max(dPMax, value[pvIdx]); + } + else if ( (pvIdx == Indices::waterSwitchIdx + && value.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) + || (pvIdx == Indices::compositionSwitchIdx + && value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) + || (enableSolvent && pvIdx == Indices::solventSaturationIdx + && value.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) + || (enableBrine && enableSaltPrecipitation && pvIdx == Indices::saltConcentrationIdx + && value.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) ) { + dSMax = std::max(dSMax, value[pvIdx]); + } + } + } + + // Communicate max values + dPMax = this->grid_.comm().max(dPMax); + dSMax = this->grid_.comm().max(dSMax); + + return {dPMax, dSMax}; } template @@ -650,7 +744,10 @@ localConvergenceData(std::vector& R_sum, } template -std::pair, std::vector> +std::tuple< +std::pair, std::vector>, +std::vector +> BlackoilModel:: characteriseCnvPvSplit(const std::vector& B_avg, const double dt) { @@ -690,6 +787,8 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) ElementContext elemCtx(this->simulator()); + std::vector ixCells; + OPM_BEGIN_PARALLEL_TRY_CATCH(); for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { // Skip cells of numerical Aquifer @@ -710,6 +809,11 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) splitPV[ix] += static_cast(pvValue); ++cellCntPV[ix]; + + // For (later) dP and dS check, we need cell indices of [1] violations + if ((this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0) && ix == 1) { + ixCells.push_back(cell_idx); + } } OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::characteriseCnvPvSplit() failed: ", @@ -718,7 +822,7 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) this->grid_.comm().sum(splitPV .data(), splitPV .size()); this->grid_.comm().sum(cellCntPV.data(), cellCntPV.size()); - return cnvPvSplit; + return { cnvPvSplit, ixCells }; } template @@ -734,6 +838,16 @@ updateTUNING(const Tuning& tuning) this->param_.newton_min_iter_ = tuning.NEWTMN; } +template +void +BlackoilModel:: +updateTUNINGDP(const TuningDp& tuning_dp) +{ + // Note: If TUNINGDP is _not_ set, these should be 0.0 + this->param_.tolerance_max_dp_ = tuning_dp.TRGDDP; + this->param_.tolerance_max_ds_ = tuning_dp.TRGDDS; +} + template ConvergenceReport BlackoilModel:: @@ -765,8 +879,8 @@ getReservoirConvergence(const double reportTime, numAquiferPvSumLocal, R_sum, maxCoeff, B_avg); - report.setCnvPoreVolSplit(this->characteriseCnvPvSplit(B_avg, dt), - pvSum - numAquiferPvSum); + auto [cnvPvSplit, ixCells] = this->characteriseCnvPvSplit(B_avg, dt); + report.setCnvPoreVolSplit(cnvPvSplit, pvSum - numAquiferPvSum); // For each iteration, we need to determine whether to use the // relaxed tolerances. To disable the usage of relaxed @@ -778,9 +892,12 @@ getReservoirConvergence(const double reportTime, const bool relax_final_iteration_mb = this->param_.min_strict_mb_iter_ < 0 && iteration == maxIter; - const bool use_relaxed_mb = relax_final_iteration_mb || - (this->param_.min_strict_mb_iter_ >= 0 && - iteration >= this->param_.min_strict_mb_iter_); + const bool relax_iter_mb = + this->param_.min_strict_mb_iter_ >= 0 && + iteration >= this->param_.min_strict_mb_iter_; + + const bool use_relaxed_mb = relax_final_iteration_mb + || relax_iter_mb; // If min_strict_cnv_iter = -1 we use a relaxed tolerance for // the cnv for the last iterations. For positive values we use @@ -809,22 +926,70 @@ getReservoirConvergence(const double reportTime, this->param_.relaxed_max_pv_fraction_ * eligible; }(); + // If tolerances for dP and/or dS are met, we use + // relaxed cnv tolerances + Scalar dPMax = 0.0; + Scalar dSMax = 0.0; + bool use_dp_tol = iteration > 0 && + this->param_.tolerance_max_dp_ > 0.0; + bool use_ds_tol = iteration > 0 && + this->param_.tolerance_max_ds_ > 0.0; + bool relax_dp_ds_cnv = false; + if (use_dp_tol || use_ds_tol) { + std::tie(dPMax, dSMax) = getMaxSolutionUpdate(ixCells); + if (use_dp_tol && !use_ds_tol) { + relax_dp_ds_cnv = dPMax < this->param_.tolerance_max_dp_; + } + else if (!use_dp_tol && use_ds_tol) { + relax_dp_ds_cnv = dSMax < this->param_.tolerance_max_ds_; + } + else { + relax_dp_ds_cnv = dPMax < this->param_.tolerance_max_dp_ + && dSMax < this->param_.tolerance_max_ds_; + } + } + + // Determine if relaxed CNV tolerances should be used const bool use_relaxed_cnv = relax_final_iteration_cnv || relax_pv_fraction_cnv - || relax_iter_cnv; + || relax_iter_cnv + || relax_dp_ds_cnv; - if ((relax_final_iteration_mb || relax_final_iteration_cnv) && + if ((use_relaxed_cnv || use_relaxed_mb) && this->terminal_output_) { - std::string message = - "Number of newton iterations reached its maximum " - "try to continue with relaxed tolerances:"; + std::string message; + if (relax_final_iteration_mb || relax_final_iteration_cnv) { + message = "Number of newton iterations reached its maximum, "; + } + else if (relax_iter_cnv || relax_iter_mb) { + message = "Number of newton iterations using strict tolerances have been reached, "; + } + else if (relax_pv_fraction_cnv) { + message = "Sum of pore volumes below tolerance, "; + } + else if (relax_dp_ds_cnv) { + if (use_dp_tol && !use_ds_tol) { + message = fmt::format("Max pressure change is below tolerance ({:.2e}<{:.2e}), ", + dPMax, this->param_.tolerance_max_dp_); + } + else if (!use_dp_tol && use_ds_tol) { + message = fmt::format("Max saturation change is below tolerance ({:.2e}<{:.2e}), ", + dSMax, this->param_.tolerance_max_ds_); + } + else { + message = fmt::format("Both pressure and saturation changes are below tolerances ({:.2e}<{:.2e} " + "and {:.2e}<{:.2e}), ", + dPMax, this->param_.tolerance_max_dp_, dSMax, this->param_.tolerance_max_ds_); + } + } + message += "try to continue with relaxed tolerances:"; - if (relax_final_iteration_mb) { + if (use_relaxed_mb) { message += fmt::format(" MB: {:.1e}", param_.tolerance_mb_relaxed_); } - if (relax_final_iteration_cnv) { + if (use_relaxed_cnv) { message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_); } @@ -910,6 +1075,9 @@ getReservoirConvergence(const double reportTime, msg += ") "; } + msg += " DP "; + msg += " DS "; + OpmLog::debug(msg); } @@ -926,6 +1094,16 @@ getReservoirConvergence(const double reportTime, ss << std::setw(11) << CNV[compIdx]; } + if (iteration == 0) { + ss << std::string(5, ' ') << "-" << std::string(5, ' '); + ss << std::string(5, ' ') << "-" << std::string(5, ' '); + } + + else { + ss << std::setw(11) << dPMax; + ss << std::setw(11) << dSMax; + } + ss.precision(oprec); ss.flags(oflags); diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp index 7c8fe7b638..10fb5d26b0 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp @@ -355,6 +355,28 @@ public: } } + void updateTUNINGDP(const TuningDp& tuning_dp) + { + // NOTE: If TUNINGDP is _not_ set, these are 0.0 + modelParam_.tolerance_max_dp_ = tuning_dp.TRGDDP; + modelParam_.tolerance_max_ds_ = tuning_dp.TRGDDS; + + // Terminal warnings + if (terminalOutput_) { + // TUNINGDP values + OpmLog::debug(fmt::format("TUNINGDP values: TRGDDP: {:.2e}, TRGDDS: {:.2e}", + tuning_dp.TRGDDP, tuning_dp.TRGDDS)); + + // Warnings unsupported items + if (tuning_dp.TRGLCV_has_value) { + OpmLog::warning("TUNINGDP item 1 (TRGLCV) is not supported."); + } + if (tuning_dp.XXXLCV_has_value) { + OpmLog::warning("TUNINGDP item 2 (XXXLCV) is not supported."); + } + } + } + bool runStep(SimulatorTimer& timer) { if (schedule().exitStatus().has_value()) { @@ -449,6 +471,19 @@ public: result = max_next_tstep > 0; } + if (enableTUNING && events.hasEvent(ScheduleEvents::TUNINGDP_CHANGE)) { + // Unset the event to not trigger it again on the next sub step + schedule.clear_event(ScheduleEvents::TUNINGDP_CHANGE, reportStep); + + // Update TUNINGDP parameters + // NOTE: Need to update both solver (model) and simulator since solver is re-created each report + // step. + const auto& sched_state = schedule[reportStep]; + const auto& tuning_dp = sched_state.tuning_dp(); + solver_->model().updateTUNINGDP(tuning_dp); + this->updateTUNINGDP(tuning_dp); + } + const auto& wcycle = schedule[reportStep].wcycle.get(); if (wcycle.empty()) { return result; diff --git a/opm/simulators/utils/UnsupportedFlowKeywords.cpp b/opm/simulators/utils/UnsupportedFlowKeywords.cpp index c3771e9a4c..88bb983148 100644 --- a/opm/simulators/utils/UnsupportedFlowKeywords.cpp +++ b/opm/simulators/utils/UnsupportedFlowKeywords.cpp @@ -631,7 +631,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords() {"TRNHD", {true, std::nullopt}}, {"TRPLPORO", {true, std::nullopt}}, {"TRROCK", {true, std::nullopt}}, - {"TUNINGDP", {false, std::nullopt}}, {"TUNINGH", {false, std::nullopt}}, {"TUNINGL", {false, std::nullopt}}, {"TUNINGS", {false, std::nullopt}}, From 59148ef5b93130b7e14afd6d21c5969939639701 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Fri, 9 Jan 2026 17:20:39 +0100 Subject: [PATCH 02/11] Check max of absolute solution change --- opm/simulators/flow/BlackoilModel_impl.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index f5dd0eb3be..ca8bab6cf5 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -609,7 +609,7 @@ getMaxSolutionUpdate(const std::vector& ixCells) const auto& value = solUpd_[ix]; for (std::size_t pvIdx = 0; pvIdx < value.size(); ++pvIdx) { if (pvIdx == Indices::pressureSwitchIdx) { - dPMax = std::max(dPMax, value[pvIdx]); + dPMax = std::max(dPMax, std::abs(value[pvIdx])); } else if ( (pvIdx == Indices::waterSwitchIdx && value.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) @@ -619,7 +619,7 @@ getMaxSolutionUpdate(const std::vector& ixCells) && value.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) || (enableBrine && enableSaltPrecipitation && pvIdx == Indices::saltConcentrationIdx && value.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) ) { - dSMax = std::max(dSMax, value[pvIdx]); + dSMax = std::max(dSMax, std::abs(value[pvIdx])); } } } From 211a48717a3d3a5713c72e6c75a7dfdbbbcebfdb Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Mon, 12 Jan 2026 11:01:45 +0100 Subject: [PATCH 03/11] Use TUNINGDP regardless of --enable-tuning parameter --- opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp index 10fb5d26b0..0e5679a835 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp @@ -471,7 +471,7 @@ public: result = max_next_tstep > 0; } - if (enableTUNING && events.hasEvent(ScheduleEvents::TUNINGDP_CHANGE)) { + if (events.hasEvent(ScheduleEvents::TUNINGDP_CHANGE)) { // Unset the event to not trigger it again on the next sub step schedule.clear_event(ScheduleEvents::TUNINGDP_CHANGE, reportStep); From 2d763a0cddeb34936f1bf06890950a798dead4e7 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Mon, 12 Jan 2026 11:13:40 +0100 Subject: [PATCH 04/11] Update TUNINGDP checks with RS and RV . Also ensure that solution updates are > 0.0 in tolerance check. --- opm/simulators/flow/BlackoilModel.hpp | 3 +- .../flow/BlackoilModelParameters.cpp | 10 ++ .../flow/BlackoilModelParameters.hpp | 10 ++ opm/simulators/flow/BlackoilModel_impl.hpp | 99 ++++++++++++------- .../flow/SimulatorFullyImplicitBlackoil.hpp | 10 +- 5 files changed, 92 insertions(+), 40 deletions(-) diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index 643b393ba0..059c350792 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -187,7 +187,8 @@ public: /// Get solution update vector as a PrimaryVarible void prepareStoringSolutionUpdate(); void storeSolutionUpdate(const BVector& dx); - std::pair getMaxSolutionUpdate(const std::vector& ixCells); + std::tuple + getMaxSolutionUpdate(const std::vector& ixCells); /// Return true if output to cout is wanted. bool terminalOutputEnabled() const diff --git a/opm/simulators/flow/BlackoilModelParameters.cpp b/opm/simulators/flow/BlackoilModelParameters.cpp index 0aa4715477..9fdd1eb540 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -51,6 +51,8 @@ BlackoilModelParameters::BlackoilModelParameters() tolerance_cnv_energy_relaxed_ = std::max(tolerance_cnv_energy_, Parameters::Get>()); tolerance_max_dp_ = Parameters::Get>(); tolerance_max_ds_ = Parameters::Get>(); + tolerance_max_drs_ = Parameters::Get>(); + tolerance_max_drv_ = Parameters::Get>(); tolerance_wells_ = Parameters::Get>(); tolerance_well_control_ = Parameters::Get>(); max_welleq_iter_ = Parameters::Get(); @@ -169,6 +171,14 @@ void BlackoilModelParameters::registerParameters() ("Tolerance for max saturation change during a Newton iteration. " "A value greater than 0.0 allows for convergence according to " "relaxed local convergence tolerance."); + Parameters::Register> + ("Tolerance for max RS change during a Newton iteration. " + "A value greater than 0.0 allows for convergence according to " + "relaxed local convergence tolerance."); + Parameters::Register> + ("Tolerance for max RV change during a Newton iteration. " + "A value greater than 0.0 allows for convergence according to " + "relaxed local convergence tolerance."); Parameters::Register> ("Well convergence tolerance"); Parameters::Register> diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index a366130b56..9ac306abbc 100644 --- a/opm/simulators/flow/BlackoilModelParameters.hpp +++ b/opm/simulators/flow/BlackoilModelParameters.hpp @@ -81,6 +81,12 @@ struct ToleranceMaxDp { static constexpr Scalar value = 0.0; }; template struct ToleranceMaxDs { static constexpr Scalar value = 0.0; }; +template +struct ToleranceMaxDrs { static constexpr Scalar value = 0.0; }; + +template +struct ToleranceMaxDrv { static constexpr Scalar value = 0.0; }; + template struct ToleranceWells { static constexpr Scalar value = 1e-4; }; @@ -222,6 +228,10 @@ public: Scalar tolerance_max_dp_; /// Max saturation change during a Newton iteration (TUNINGDP item = TRGDDS) Scalar tolerance_max_ds_; + /// Max RS change during a Newton iteration (TUNINGDP item = TRGDDRS) + Scalar tolerance_max_drs_; + /// Max RV change during a Newton iteration (TUNINGDP item = TRGDDRV) + Scalar tolerance_max_drv_; /// Well convergence tolerance. Scalar tolerance_wells_; /// Tolerance for the well control equations diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index ca8bab6cf5..90ebd07667 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -521,7 +521,8 @@ updateSolution(const BVector& dx) { OPM_TIMEBLOCK(updateSolution); // Prepare to store solution update for convergence check - if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0) { + if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_) { prepareStoringSolutionUpdate(); } @@ -542,7 +543,8 @@ updateSolution(const BVector& dx) } // Store solution update - if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0) { + if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_) { storeSolutionUpdate(dx); } } @@ -592,8 +594,10 @@ storeSolutionUpdate(const BVector& dx) } template -std::pair::Scalar, - typename BlackoilModel::Scalar> +std::tuple::Scalar, + typename BlackoilModel::Scalar, + typename BlackoilModel::Scalar, + typename BlackoilModel::Scalar> BlackoilModel:: getMaxSolutionUpdate(const std::vector& ixCells) { @@ -603,6 +607,8 @@ getMaxSolutionUpdate(const std::vector& ixCells) // Init output Scalar dPMax = 0.0; Scalar dSMax = 0.0; + Scalar dRsMax = 0.0; + Scalar dRvMax = 0.0; // Loop over solution update, get the correct variables and calculate max. for (const auto& ix : ixCells) { @@ -621,14 +627,24 @@ getMaxSolutionUpdate(const std::vector& ixCells) && value.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) ) { dSMax = std::max(dSMax, std::abs(value[pvIdx])); } + else if (pvIdx == Indices::compositionSwitchIdx + && value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) { + dRsMax = std::max(dRsMax, std::abs(value[pvIdx])); + } + else if (pvIdx == Indices::compositionSwitchIdx + && value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) { + dRvMax = std::max(dRvMax, std::abs(value[pvIdx])); + } } } // Communicate max values dPMax = this->grid_.comm().max(dPMax); dSMax = this->grid_.comm().max(dSMax); + dRsMax = this->grid_.comm().max(dRsMax); + dRvMax = this->grid_.comm().max(dRvMax); - return {dPMax, dSMax}; + return { dPMax, dSMax, dRsMax, dRvMax }; } template @@ -810,8 +826,10 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) splitPV[ix] += static_cast(pvValue); ++cellCntPV[ix]; - // For (later) dP and dS check, we need cell indices of [1] violations - if ((this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0) && ix == 1) { + // For dP and dS check, we need cell indices of [1] violations + if ( ix == 1 && + (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0 ) ) { ixCells.push_back(cell_idx); } } @@ -843,9 +861,11 @@ void BlackoilModel:: updateTUNINGDP(const TuningDp& tuning_dp) { - // Note: If TUNINGDP is _not_ set, these should be 0.0 + // NOTE: If TUNINGDP item is _not_ set it should be 0.0 this->param_.tolerance_max_dp_ = tuning_dp.TRGDDP; this->param_.tolerance_max_ds_ = tuning_dp.TRGDDS; + this->param_.tolerance_max_drs_ = tuning_dp.TRGDDRS; + this->param_.tolerance_max_drv_ = tuning_dp.TRGDDRV; } template @@ -926,34 +946,35 @@ getReservoirConvergence(const double reportTime, this->param_.relaxed_max_pv_fraction_ * eligible; }(); - // If tolerances for dP and/or dS are met, we use - // relaxed cnv tolerances + // If tolerances for solution changes are met, we use the + // relaxed cnv tolerance. Note that all tolerances > 0.0 + // must be met to enable use of relaxed cnv tolerance. Scalar dPMax = 0.0; Scalar dSMax = 0.0; + Scalar dRsMax = 0.0; + Scalar dRvMax = 0.0; bool use_dp_tol = iteration > 0 && this->param_.tolerance_max_dp_ > 0.0; bool use_ds_tol = iteration > 0 && this->param_.tolerance_max_ds_ > 0.0; - bool relax_dp_ds_cnv = false; - if (use_dp_tol || use_ds_tol) { - std::tie(dPMax, dSMax) = getMaxSolutionUpdate(ixCells); - if (use_dp_tol && !use_ds_tol) { - relax_dp_ds_cnv = dPMax < this->param_.tolerance_max_dp_; - } - else if (!use_dp_tol && use_ds_tol) { - relax_dp_ds_cnv = dSMax < this->param_.tolerance_max_ds_; - } - else { - relax_dp_ds_cnv = dPMax < this->param_.tolerance_max_dp_ - && dSMax < this->param_.tolerance_max_ds_; - } + bool use_drs_tol = iteration > 0 && + this->param_.tolerance_max_drs_ > 0.0; + bool use_drv_tol = iteration > 0 && + this->param_.tolerance_max_drv_ > 0.0; + bool relax_dsol_cnv = false; + if (!ixCells.empty() && (use_dp_tol || use_ds_tol || use_drs_tol || use_drv_tol)) { + std::tie(dPMax, dSMax, dRsMax, dRvMax) = getMaxSolutionUpdate(ixCells); + relax_dsol_cnv = (!use_dp_tol || (dPMax > 0.0 && dPMax < this->param_.tolerance_max_dp_)) && + (!use_ds_tol || (dSMax > 0.0 && dSMax < this->param_.tolerance_max_ds_)) && + (!use_drs_tol || (dRsMax > 0.0 && dRsMax < this->param_.tolerance_max_drs_)) && + (!use_drv_tol || (dRvMax > 0.0 && dRvMax < this->param_.tolerance_max_drv_)); } // Determine if relaxed CNV tolerances should be used const bool use_relaxed_cnv = relax_final_iteration_cnv || relax_pv_fraction_cnv || relax_iter_cnv - || relax_dp_ds_cnv; + || relax_dsol_cnv; if ((use_relaxed_cnv || use_relaxed_mb) && this->terminal_output_) @@ -968,20 +989,21 @@ getReservoirConvergence(const double reportTime, else if (relax_pv_fraction_cnv) { message = "Sum of pore volumes below tolerance, "; } - else if (relax_dp_ds_cnv) { - if (use_dp_tol && !use_ds_tol) { - message = fmt::format("Max pressure change is below tolerance ({:.2e}<{:.2e}), ", - dPMax, this->param_.tolerance_max_dp_); + else if (relax_dsol_cnv) { + message = "Max. solution change below tolerance ( "; + if (use_dp_tol && dPMax < this->param_.tolerance_max_dp_ ) { + message += fmt::format("DP={:.2e}<{:.2e} ", dPMax, this->param_.tolerance_max_dp_); } - else if (!use_dp_tol && use_ds_tol) { - message = fmt::format("Max saturation change is below tolerance ({:.2e}<{:.2e}), ", - dSMax, this->param_.tolerance_max_ds_); + if (use_ds_tol && dSMax < this->param_.tolerance_max_ds_ ) { + message += fmt::format("DS={:.2e}<{:.2e} ", dSMax, this->param_.tolerance_max_ds_); } - else { - message = fmt::format("Both pressure and saturation changes are below tolerances ({:.2e}<{:.2e} " - "and {:.2e}<{:.2e}), ", - dPMax, this->param_.tolerance_max_dp_, dSMax, this->param_.tolerance_max_ds_); + if (use_drs_tol && dRsMax < this->param_.tolerance_max_drs_ ) { + message += fmt::format("DRS={:.2e}<{:.2e} ", dRsMax, this->param_.tolerance_max_drs_); } + if (use_drv_tol && dRvMax < this->param_.tolerance_max_drv_ ) { + message += fmt::format("DRV={:.2e}<{:.2e} ", dRvMax, this->param_.tolerance_max_drv_); + } + message += "), "; } message += "try to continue with relaxed tolerances:"; @@ -1077,6 +1099,8 @@ getReservoirConvergence(const double reportTime, msg += " DP "; msg += " DS "; + msg += " DRS "; + msg += " DRV "; OpmLog::debug(msg); } @@ -1097,11 +1121,14 @@ getReservoirConvergence(const double reportTime, if (iteration == 0) { ss << std::string(5, ' ') << "-" << std::string(5, ' '); ss << std::string(5, ' ') << "-" << std::string(5, ' '); + ss << std::string(5, ' ') << "-" << std::string(5, ' '); + ss << std::string(5, ' ') << "-" << std::string(5, ' '); } - else { ss << std::setw(11) << dPMax; ss << std::setw(11) << dSMax; + ss << std::setw(11) << dRsMax; + ss << std::setw(11) << dRvMax; } ss.precision(oprec); diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp index 0e5679a835..2be9f562cd 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp @@ -357,15 +357,19 @@ public: void updateTUNINGDP(const TuningDp& tuning_dp) { - // NOTE: If TUNINGDP is _not_ set, these are 0.0 + // NOTE: If TUNINGDP item is _not_ set it should be 0.0 modelParam_.tolerance_max_dp_ = tuning_dp.TRGDDP; modelParam_.tolerance_max_ds_ = tuning_dp.TRGDDS; + modelParam_.tolerance_max_drs_ = tuning_dp.TRGDDRS; + modelParam_.tolerance_max_drv_ = tuning_dp.TRGDDRV; // Terminal warnings if (terminalOutput_) { // TUNINGDP values - OpmLog::debug(fmt::format("TUNINGDP values: TRGDDP: {:.2e}, TRGDDS: {:.2e}", - tuning_dp.TRGDDP, tuning_dp.TRGDDS)); + OpmLog::debug( + fmt::format("TUNINGDP values: TRGDDP: {:.2e}, TRGDDS: {:.2e}, TRGDDRS: {:.2e}, TRGDDRV: {:.2e}", + tuning_dp.TRGDDP, tuning_dp.TRGDDS, tuning_dp.TRGDDRS, tuning_dp.TRGDDRV) + ); // Warnings unsupported items if (tuning_dp.TRGLCV_has_value) { From 283fbdd76d3477ec2ba6110ec6374796b19f7483 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Mon, 12 Jan 2026 11:17:31 +0100 Subject: [PATCH 05/11] Disregard CNV tolerances if TUNINGDP checks are met --- opm/simulators/flow/BlackoilModelParameters.cpp | 16 ++++++++-------- opm/simulators/flow/BlackoilModel_impl.hpp | 8 ++++++-- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelParameters.cpp b/opm/simulators/flow/BlackoilModelParameters.cpp index 9fdd1eb540..c6b9c94bc1 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -165,20 +165,20 @@ void BlackoilModelParameters::registerParameters() "after the iterations with the strict tolerance"); Parameters::Register> ("Tolerance for max pressure change during a Newton iteration. " - "A value greater than 0.0 allows for convergence according to " - "relaxed local convergence tolerance"); + "A value greater than 0.0 allows for convergence regardless " + "of residual tolerances. Use with care!"); Parameters::Register> ("Tolerance for max saturation change during a Newton iteration. " - "A value greater than 0.0 allows for convergence according to " - "relaxed local convergence tolerance."); + "A value greater than 0.0 allows for convergence regardless " + "of residual tolerances. Use with care!"); Parameters::Register> ("Tolerance for max RS change during a Newton iteration. " - "A value greater than 0.0 allows for convergence according to " - "relaxed local convergence tolerance."); + "A value greater than 0.0 allows for convergence regardless " + "of residual tolerances. Use with care!"); Parameters::Register> ("Tolerance for max RV change during a Newton iteration. " - "A value greater than 0.0 allows for convergence according to " - "relaxed local convergence tolerance."); + "A value greater than 0.0 allows for convergence regardless " + "of residual tolerances. Use with care!"); Parameters::Register> ("Well convergence tolerance"); Parameters::Register> diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index 90ebd07667..a8fc32170e 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -976,6 +976,10 @@ getReservoirConvergence(const double reportTime, || relax_iter_cnv || relax_dsol_cnv; + // Ensure that CNV convergence criteria is met when max. + // solution change tolerances have been fulfilled + Scalar tolerance_cnv_relaxed = relax_dsol_cnv ? 1e20 : param_.tolerance_cnv_relaxed_; + if ((use_relaxed_cnv || use_relaxed_mb) && this->terminal_output_) { @@ -1012,13 +1016,13 @@ getReservoirConvergence(const double reportTime, } if (use_relaxed_cnv) { - message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_); + message += fmt::format(" CNV: {:.1e}", tolerance_cnv_relaxed); } OpmLog::debug(message); } - const auto tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_; + const auto tol_cnv = use_relaxed_cnv ? tolerance_cnv_relaxed : param_.tolerance_cnv_; const auto tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_; const auto tol_cnv_energy = use_relaxed_cnv ? param_.tolerance_cnv_energy_relaxed_ : param_.tolerance_cnv_energy_; const auto tol_eb = use_relaxed_mb ? param_.tolerance_energy_balance_relaxed_ : param_.tolerance_energy_balance_; From 571397db967f05641edcd9178ef41b81d81d70bd Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Mon, 12 Jan 2026 12:54:24 +0100 Subject: [PATCH 06/11] Check max solution updates for all CNV>tol_cnv violations --- opm/simulators/flow/BlackoilModel_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index a8fc32170e..f5a469ffed 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -827,7 +827,7 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) ++cellCntPV[ix]; // For dP and dS check, we need cell indices of [1] violations - if ( ix == 1 && + if ( ix > 0 && (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0 ) ) { ixCells.push_back(cell_idx); From 16c85068964af9e52463839d529c11688e272e0a Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Tue, 13 Jan 2026 12:37:33 +0100 Subject: [PATCH 07/11] Removed unnecessary check for empty vector --- opm/simulators/flow/BlackoilModel_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index f5a469ffed..87e7adadaa 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -962,7 +962,7 @@ getReservoirConvergence(const double reportTime, bool use_drv_tol = iteration > 0 && this->param_.tolerance_max_drv_ > 0.0; bool relax_dsol_cnv = false; - if (!ixCells.empty() && (use_dp_tol || use_ds_tol || use_drs_tol || use_drv_tol)) { + if (use_dp_tol || use_ds_tol || use_drs_tol || use_drv_tol) { std::tie(dPMax, dSMax, dRsMax, dRvMax) = getMaxSolutionUpdate(ixCells); relax_dsol_cnv = (!use_dp_tol || (dPMax > 0.0 && dPMax < this->param_.tolerance_max_dp_)) && (!use_ds_tol || (dSMax > 0.0 && dSMax < this->param_.tolerance_max_ds_)) && From 6ecb37e19a5f51936f318293fcb188f52ee662e8 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Wed, 14 Jan 2026 10:09:32 +0100 Subject: [PATCH 08/11] Address review comments and warnings. Additionally, cleaned up debug output. --- opm/simulators/flow/BlackoilModel.hpp | 27 +++-- opm/simulators/flow/BlackoilModel_impl.hpp | 132 ++++++++++----------- 2 files changed, 81 insertions(+), 78 deletions(-) diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index 059c350792..4d7eb8857e 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -110,6 +110,19 @@ public: using ComponentName = ::Opm::ComponentName; + // Helper structs + struct CnvPvSplitData { + std::pair, std::vector> cnvPvSplit; + std::vector ixCells; + }; + + struct MaxSolutionUpdateData { + Scalar dPMax = 0.0; + Scalar dSMax = 0.0; + Scalar dRsMax = 0.0; + Scalar dRvMax = 0.0; + }; + // --------- Public methods --------- /// Construct the model. It will retain references to the @@ -184,11 +197,10 @@ public: /// Apply an update to the primary variables. void updateSolution(const BVector& dx); - /// Get solution update vector as a PrimaryVarible + /// Get solution update vector as a PrimaryVariable void prepareStoringSolutionUpdate(); void storeSolutionUpdate(const BVector& dx); - std::tuple - getMaxSolutionUpdate(const std::vector& ixCells); + MaxSolutionUpdateData getMaxSolutionUpdate(const std::vector& ixCells); /// Return true if output to cout is wanted. bool terminalOutputEnabled() const @@ -213,12 +225,9 @@ public: /// \brief Compute pore-volume/cell count split among "converged", /// "relaxed converged", "unconverged" cells based on CNV point - /// measures. - std::tuple< - std::pair, std::vector>, - std::vector - > - characteriseCnvPvSplit(const std::vector& B_avg, const double dt); + /// measures. Also returns list of cells where CNV is greater than + /// its strict tolerance + CnvPvSplitData characteriseCnvPvSplit(const std::vector& B_avg, const double dt); /// \brief Compute the number of Newtons required by each cell in order to /// satisfy the solution change convergence criteria at the last time step. diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index 87e7adadaa..43e4814211 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -566,9 +566,7 @@ prepareStoringSolutionUpdate() solUpd_[globalElemIdx] = simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; // Ensure each element is zero - for (auto& elem : solUpd_[globalElemIdx]) { - elem = 0.0; - } + std::fill(solUpd_[globalElemIdx].begin(), solUpd_[globalElemIdx].end(), 0.0); } } @@ -583,21 +581,16 @@ storeSolutionUpdate(const BVector& dx) // Get cell vectors unsigned globalElemIdx = elemMapper.index(elem); auto& value = solUpd_[globalElemIdx]; - auto& update = dx[globalElemIdx]; + const auto& update = dx[globalElemIdx]; assert(value.size() == update.size()); // Transfer update from dx to solution update container (SolutionVector type) - for (std::size_t i = 0; i < value.size(); ++i) { - value[i] = update[i]; - } + std::copy(update.begin(), update.end(), value.begin()); } } template -std::tuple::Scalar, - typename BlackoilModel::Scalar, - typename BlackoilModel::Scalar, - typename BlackoilModel::Scalar> +typename BlackoilModel::MaxSolutionUpdateData BlackoilModel:: getMaxSolutionUpdate(const std::vector& ixCells) { @@ -613,7 +606,7 @@ getMaxSolutionUpdate(const std::vector& ixCells) // Loop over solution update, get the correct variables and calculate max. for (const auto& ix : ixCells) { const auto& value = solUpd_[ix]; - for (std::size_t pvIdx = 0; pvIdx < value.size(); ++pvIdx) { + for (int pvIdx = 0; pvIdx < value.size(); ++pvIdx) { if (pvIdx == Indices::pressureSwitchIdx) { dPMax = std::max(dPMax, std::abs(value[pvIdx])); } @@ -760,10 +753,7 @@ localConvergenceData(std::vector& R_sum, } template -std::tuple< -std::pair, std::vector>, -std::vector -> +typename BlackoilModel::CnvPvSplitData BlackoilModel:: characteriseCnvPvSplit(const std::vector& B_avg, const double dt) { @@ -899,8 +889,9 @@ getReservoirConvergence(const double reportTime, numAquiferPvSumLocal, R_sum, maxCoeff, B_avg); - auto [cnvPvSplit, ixCells] = this->characteriseCnvPvSplit(B_avg, dt); - report.setCnvPoreVolSplit(cnvPvSplit, pvSum - numAquiferPvSum); + auto cnvSplitData = this->characteriseCnvPvSplit(B_avg, dt); + report.setCnvPoreVolSplit(cnvSplitData.cnvPvSplit, + pvSum - numAquiferPvSum); // For each iteration, we need to determine whether to use the // relaxed tolerances. To disable the usage of relaxed @@ -949,25 +940,20 @@ getReservoirConvergence(const double reportTime, // If tolerances for solution changes are met, we use the // relaxed cnv tolerance. Note that all tolerances > 0.0 // must be met to enable use of relaxed cnv tolerance. - Scalar dPMax = 0.0; - Scalar dSMax = 0.0; - Scalar dRsMax = 0.0; - Scalar dRvMax = 0.0; - bool use_dp_tol = iteration > 0 && - this->param_.tolerance_max_dp_ > 0.0; - bool use_ds_tol = iteration > 0 && - this->param_.tolerance_max_ds_ > 0.0; - bool use_drs_tol = iteration > 0 && - this->param_.tolerance_max_drs_ > 0.0; - bool use_drv_tol = iteration > 0 && - this->param_.tolerance_max_drv_ > 0.0; + MaxSolutionUpdateData maxSolUpd; + const bool use_dp_tol = this->param_.tolerance_max_dp_ > 0.0; + const bool use_ds_tol = this->param_.tolerance_max_ds_ > 0.0; + const bool use_drs_tol = this->param_.tolerance_max_drs_ > 0.0; + const bool use_drv_tol = this->param_.tolerance_max_drv_ > 0.0; + const bool use_dsol_tol = use_dp_tol || use_ds_tol || use_drs_tol || use_drv_tol; bool relax_dsol_cnv = false; - if (use_dp_tol || use_ds_tol || use_drs_tol || use_drv_tol) { - std::tie(dPMax, dSMax, dRsMax, dRvMax) = getMaxSolutionUpdate(ixCells); - relax_dsol_cnv = (!use_dp_tol || (dPMax > 0.0 && dPMax < this->param_.tolerance_max_dp_)) && - (!use_ds_tol || (dSMax > 0.0 && dSMax < this->param_.tolerance_max_ds_)) && - (!use_drs_tol || (dRsMax > 0.0 && dRsMax < this->param_.tolerance_max_drs_)) && - (!use_drv_tol || (dRvMax > 0.0 && dRvMax < this->param_.tolerance_max_drv_)); + if (iteration > 0 && use_dsol_tol) { + maxSolUpd = getMaxSolutionUpdate(cnvSplitData.ixCells); + relax_dsol_cnv = + (!use_dp_tol || (maxSolUpd.dPMax > 0.0 && maxSolUpd.dPMax < this->param_.tolerance_max_dp_)) && + (!use_ds_tol || (maxSolUpd.dSMax > 0.0 && maxSolUpd.dSMax < this->param_.tolerance_max_ds_)) && + (!use_drs_tol || (maxSolUpd.dRsMax > 0.0 && maxSolUpd.dRsMax < this->param_.tolerance_max_drs_)) && + (!use_drv_tol || (maxSolUpd.dRvMax > 0.0 && maxSolUpd.dRvMax < this->param_.tolerance_max_drv_)); } // Determine if relaxed CNV tolerances should be used @@ -984,30 +970,30 @@ getReservoirConvergence(const double reportTime, this->terminal_output_) { std::string message; - if (relax_final_iteration_mb || relax_final_iteration_cnv) { + if (relax_dsol_cnv) { + message = "Max. solution change below tolerance ( "; + if (use_dp_tol && maxSolUpd.dPMax < this->param_.tolerance_max_dp_ ) { + message += fmt::format("DP={:.2e}<{:.2e} ", maxSolUpd.dPMax, this->param_.tolerance_max_dp_); + } + if (use_ds_tol && maxSolUpd.dSMax < this->param_.tolerance_max_ds_ ) { + message += fmt::format("DS={:.2e}<{:.2e} ", maxSolUpd.dSMax, this->param_.tolerance_max_ds_); + } + if (use_drs_tol && maxSolUpd.dRsMax < this->param_.tolerance_max_drs_ ) { + message += fmt::format("DRS={:.2e}<{:.2e} ", maxSolUpd.dRsMax, this->param_.tolerance_max_drs_); + } + if (use_drv_tol && maxSolUpd.dRvMax < this->param_.tolerance_max_drv_ ) { + message += fmt::format("DRV={:.2e}<{:.2e} ", maxSolUpd.dRvMax, this->param_.tolerance_max_drv_); + } + message += "), "; + } + else if (relax_final_iteration_cnv || relax_final_iteration_mb) { message = "Number of newton iterations reached its maximum, "; } else if (relax_iter_cnv || relax_iter_mb) { message = "Number of newton iterations using strict tolerances have been reached, "; } else if (relax_pv_fraction_cnv) { - message = "Sum of pore volumes below tolerance, "; - } - else if (relax_dsol_cnv) { - message = "Max. solution change below tolerance ( "; - if (use_dp_tol && dPMax < this->param_.tolerance_max_dp_ ) { - message += fmt::format("DP={:.2e}<{:.2e} ", dPMax, this->param_.tolerance_max_dp_); - } - if (use_ds_tol && dSMax < this->param_.tolerance_max_ds_ ) { - message += fmt::format("DS={:.2e}<{:.2e} ", dSMax, this->param_.tolerance_max_ds_); - } - if (use_drs_tol && dRsMax < this->param_.tolerance_max_drs_ ) { - message += fmt::format("DRS={:.2e}<{:.2e} ", dRsMax, this->param_.tolerance_max_drs_); - } - if (use_drv_tol && dRvMax < this->param_.tolerance_max_drv_ ) { - message += fmt::format("DRV={:.2e}<{:.2e} ", dRvMax, this->param_.tolerance_max_drv_); - } - message += "), "; + message = "Fraction of pore volumes violating strict tolerances small, "; } message += "try to continue with relaxed tolerances:"; @@ -1101,10 +1087,12 @@ getReservoirConvergence(const double reportTime, msg += ") "; } - msg += " DP "; - msg += " DS "; - msg += " DRS "; - msg += " DRV "; + if (use_dsol_tol) { + msg += use_dp_tol ? " DP " : ""; + msg += use_ds_tol ? " DS " : ""; + msg += use_drs_tol ? " DRS " : ""; + msg += use_drv_tol ? " DRV " : ""; + } OpmLog::debug(msg); } @@ -1122,17 +1110,23 @@ getReservoirConvergence(const double reportTime, ss << std::setw(11) << CNV[compIdx]; } - if (iteration == 0) { - ss << std::string(5, ' ') << "-" << std::string(5, ' '); - ss << std::string(5, ' ') << "-" << std::string(5, ' '); - ss << std::string(5, ' ') << "-" << std::string(5, ' '); - ss << std::string(5, ' ') << "-" << std::string(5, ' '); - } - else { - ss << std::setw(11) << dPMax; - ss << std::setw(11) << dSMax; - ss << std::setw(11) << dRsMax; - ss << std::setw(11) << dRvMax; + if (use_dsol_tol) { + auto print_dsol = + [&] (bool use_tol, Scalar dsol) { + if (!use_tol) { + return; + } + if (iteration == 0 || dsol <= 0.0) { + ss << std::string(5, ' ') << "-" << std::string(5, ' '); + } + else { + ss << std::setw(11) << dsol; + } + }; + print_dsol(use_dp_tol, maxSolUpd.dPMax); + print_dsol(use_ds_tol, maxSolUpd.dSMax); + print_dsol(use_drs_tol, maxSolUpd.dRsMax); + print_dsol(use_drv_tol, maxSolUpd.dRvMax); } ss.precision(oprec); From c621f5254dc3239bb6264015fd0008b0c2e998d5 Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Thu, 22 Jan 2026 13:03:42 +0100 Subject: [PATCH 09/11] Removed debug message and fix tolerance parameter checks --- opm/simulators/flow/BlackoilModel_impl.hpp | 4 ++-- opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp | 6 ------ 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index 43e4814211..e06f20fa81 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -522,7 +522,7 @@ updateSolution(const BVector& dx) OPM_TIMEBLOCK(updateSolution); // Prepare to store solution update for convergence check if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 - || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_) { + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) { prepareStoringSolutionUpdate(); } @@ -544,7 +544,7 @@ updateSolution(const BVector& dx) // Store solution update if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0 - || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_) { + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) { storeSolutionUpdate(dx); } } diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp index 2be9f562cd..d7a5a82e93 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp @@ -365,12 +365,6 @@ public: // Terminal warnings if (terminalOutput_) { - // TUNINGDP values - OpmLog::debug( - fmt::format("TUNINGDP values: TRGDDP: {:.2e}, TRGDDS: {:.2e}, TRGDDRS: {:.2e}, TRGDDRV: {:.2e}", - tuning_dp.TRGDDP, tuning_dp.TRGDDS, tuning_dp.TRGDDRS, tuning_dp.TRGDDRV) - ); - // Warnings unsupported items if (tuning_dp.TRGLCV_has_value) { OpmLog::warning("TUNINGDP item 1 (TRGLCV) is not supported."); From c0e5851486bdcaa4165ccd1a803005bebbbfb15e Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Thu, 22 Jan 2026 13:05:37 +0100 Subject: [PATCH 10/11] Replace relaxed tol. messages with debug flags --- opm/simulators/flow/BlackoilModel.hpp | 7 +++ opm/simulators/flow/BlackoilModel_impl.hpp | 51 ++++------------------ 2 files changed, 16 insertions(+), 42 deletions(-) diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index 4d7eb8857e..276f7c44f7 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -123,6 +123,13 @@ public: Scalar dRvMax = 0.0; }; + // Output debug flags for which tolerances used + enum class DebugFlags { + STRICT = 0, + RELAXED = 1, + TUNINGDP = 2 + }; + // --------- Public methods --------- /// Construct the model. It will retain references to the diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index e06f20fa81..add68401c4 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -966,48 +966,6 @@ getReservoirConvergence(const double reportTime, // solution change tolerances have been fulfilled Scalar tolerance_cnv_relaxed = relax_dsol_cnv ? 1e20 : param_.tolerance_cnv_relaxed_; - if ((use_relaxed_cnv || use_relaxed_mb) && - this->terminal_output_) - { - std::string message; - if (relax_dsol_cnv) { - message = "Max. solution change below tolerance ( "; - if (use_dp_tol && maxSolUpd.dPMax < this->param_.tolerance_max_dp_ ) { - message += fmt::format("DP={:.2e}<{:.2e} ", maxSolUpd.dPMax, this->param_.tolerance_max_dp_); - } - if (use_ds_tol && maxSolUpd.dSMax < this->param_.tolerance_max_ds_ ) { - message += fmt::format("DS={:.2e}<{:.2e} ", maxSolUpd.dSMax, this->param_.tolerance_max_ds_); - } - if (use_drs_tol && maxSolUpd.dRsMax < this->param_.tolerance_max_drs_ ) { - message += fmt::format("DRS={:.2e}<{:.2e} ", maxSolUpd.dRsMax, this->param_.tolerance_max_drs_); - } - if (use_drv_tol && maxSolUpd.dRvMax < this->param_.tolerance_max_drv_ ) { - message += fmt::format("DRV={:.2e}<{:.2e} ", maxSolUpd.dRvMax, this->param_.tolerance_max_drv_); - } - message += "), "; - } - else if (relax_final_iteration_cnv || relax_final_iteration_mb) { - message = "Number of newton iterations reached its maximum, "; - } - else if (relax_iter_cnv || relax_iter_mb) { - message = "Number of newton iterations using strict tolerances have been reached, "; - } - else if (relax_pv_fraction_cnv) { - message = "Fraction of pore volumes violating strict tolerances small, "; - } - message += "try to continue with relaxed tolerances:"; - - if (use_relaxed_mb) { - message += fmt::format(" MB: {:.1e}", param_.tolerance_mb_relaxed_); - } - - if (use_relaxed_cnv) { - message += fmt::format(" CNV: {:.1e}", tolerance_cnv_relaxed); - } - - OpmLog::debug(message); - } - const auto tol_cnv = use_relaxed_cnv ? tolerance_cnv_relaxed : param_.tolerance_cnv_; const auto tol_mb = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_; const auto tol_cnv_energy = use_relaxed_cnv ? param_.tolerance_cnv_energy_relaxed_ : param_.tolerance_cnv_energy_; @@ -1094,6 +1052,9 @@ getReservoirConvergence(const double reportTime, msg += use_drv_tol ? " DRV " : ""; } + msg += " MBFLAG"; + msg += " CNVFLAG"; + OpmLog::debug(msg); } @@ -1129,6 +1090,12 @@ getReservoirConvergence(const double reportTime, print_dsol(use_drv_tol, maxSolUpd.dRvMax); } + int mb_flag = use_relaxed_mb ? static_cast(DebugFlags::RELAXED) : static_cast(DebugFlags::STRICT); + int cnv_flag = relax_dsol_cnv ? static_cast(DebugFlags::TUNINGDP) : + (use_relaxed_cnv ? static_cast(DebugFlags::RELAXED) : static_cast(DebugFlags::STRICT)); + ss << std::setw(7) << mb_flag; + ss << std::setw(8) << cnv_flag; + ss.precision(oprec); ss.flags(oflags); From 708af98b7ecad450edc09f783a5ea67527ffa58f Mon Sep 17 00:00:00 2001 From: Svenn Tveit Date: Thu, 22 Jan 2026 13:07:55 +0100 Subject: [PATCH 11/11] Only set relaxed CNV flag for PV fraction if > 0.0 --- opm/simulators/flow/BlackoilModel_impl.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index add68401c4..8b9713f301 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -933,8 +933,9 @@ getReservoirConvergence(const double reportTime, // [1]: tol < cnv <= relaxed // [2]: relaxed < cnv - return static_cast(cnvPvSplit[1] + cnvPvSplit[2]) < - this->param_.relaxed_max_pv_fraction_ * eligible; + Scalar cnvPvSum = static_cast(cnvPvSplit[1] + cnvPvSplit[2]); + return cnvPvSum < this->param_.relaxed_max_pv_fraction_ * eligible && + cnvPvSum > 0.0; }(); // If tolerances for solution changes are met, we use the