diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index c7bbc5038d..af0b830757 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; @@ -108,6 +110,26 @@ 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; + }; + + // 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 @@ -177,6 +199,11 @@ public: /// Apply an update to the primary variables. void updateSolution(const BVector& dx); + /// Get solution update vector as a PrimaryVariable + void prepareStoringSolutionUpdate(); + void storeSolutionUpdate(const BVector& dx); + MaxSolutionUpdateData getMaxSolutionUpdate(const std::vector& ixCells); + /// Return true if output to cout is wanted. bool terminalOutputEnabled() const { return terminal_output_; } @@ -200,9 +227,9 @@ public: /// \brief Compute pore-volume/cell count split among "converged", /// "relaxed converged", "unconverged" cells based on CNV point - /// measures. - std::pair, 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. @@ -213,6 +240,7 @@ public: const int iteration); void updateTUNING(const Tuning& tuning); + void updateTUNINGDP(const TuningDp& tuning_dp); ConvergenceReport getReservoirConvergence(const double reportTime, @@ -339,6 +367,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..c6b9c94bc1 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -49,6 +49,10 @@ 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_max_drs_ = Parameters::Get>(); + tolerance_max_drv_ = Parameters::Get>(); tolerance_wells_ = Parameters::Get>(); tolerance_well_control_ = Parameters::Get>(); max_welleq_iter_ = Parameters::Get(); @@ -159,6 +163,22 @@ 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 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 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 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 regardless " + "of residual tolerances. Use with care!"); Parameters::Register> ("Well convergence tolerance"); Parameters::Register> diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index 55ed751d5e..9ac306abbc 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,18 @@ 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 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; }; @@ -212,6 +224,14 @@ 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_; + /// 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 59b6402fdb..fb810cd687 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -507,6 +507,12 @@ 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 + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) { + prepareStoringSolutionUpdate(); + } + auto& newtonMethod = simulator_.model().newtonMethod(); SolutionVector& solution = simulator_.model().solution(/*timeIdx=*/0); @@ -522,6 +528,103 @@ 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 + || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 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 + std::fill(solUpd_[globalElemIdx].begin(), solUpd_[globalElemIdx].end(), 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]; + const auto& update = dx[globalElemIdx]; + assert(value.size() == update.size()); + + // Transfer update from dx to solution update container (SolutionVector type) + std::copy(update.begin(), update.end(), value.begin()); + } +} + +template +typename BlackoilModel::MaxSolutionUpdateData +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; + Scalar dRsMax = 0.0; + Scalar dRvMax = 0.0; + + // Loop over solution update, get the correct variables and calculate max. + for (const auto& ix : ixCells) { + const auto& value = solUpd_[ix]; + for (int pvIdx = 0; pvIdx < value.size(); ++pvIdx) { + if (pvIdx == Indices::pressureSwitchIdx) { + dPMax = std::max(dPMax, std::abs(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, 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, dRsMax, dRvMax }; } template @@ -637,7 +740,7 @@ localConvergenceData(std::vector& R_sum, } template -std::pair, std::vector> +typename BlackoilModel::CnvPvSplitData BlackoilModel:: characteriseCnvPvSplit(const std::vector& B_avg, const double dt) { @@ -677,6 +780,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 @@ -697,6 +802,13 @@ characteriseCnvPvSplit(const std::vector& B_avg, const double dt) splitPV[ix] += static_cast(pvValue); ++cellCntPV[ix]; + + // For dP and dS check, we need cell indices of [1] violations + 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); + } } OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::characteriseCnvPvSplit() failed: ", @@ -705,7 +817,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 @@ -721,6 +833,18 @@ updateTUNING(const Tuning& tuning) this->param_.newton_min_iter_ = tuning.NEWTMN; } +template +void +BlackoilModel:: +updateTUNINGDP(const TuningDp& tuning_dp) +{ + // 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 ConvergenceReport BlackoilModel:: @@ -752,7 +876,8 @@ getReservoirConvergence(const double reportTime, numAquiferPvSumLocal, R_sum, maxCoeff, B_avg); - report.setCnvPoreVolSplit(this->characteriseCnvPvSplit(B_avg, dt), + auto cnvSplitData = this->characteriseCnvPvSplit(B_avg, dt); + report.setCnvPoreVolSplit(cnvSplitData.cnvPvSplit, pvSum - numAquiferPvSum); // For each iteration, we need to determine whether to use the @@ -765,9 +890,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 @@ -792,33 +920,41 @@ 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; }(); - const bool use_relaxed_cnv = relax_final_iteration_cnv - || relax_pv_fraction_cnv - || relax_iter_cnv; - - if ((relax_final_iteration_mb || relax_final_iteration_cnv) && - this->terminal_output_) - { - std::string message = - "Number of newton iterations reached its maximum " - "try to continue with relaxed tolerances:"; - - if (relax_final_iteration_mb) { - message += fmt::format(" MB: {:.1e}", param_.tolerance_mb_relaxed_); - } - - if (relax_final_iteration_cnv) { - message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_); - } - - OpmLog::debug(message); + // 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. + 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 (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_)); } - const auto tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_; + // 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_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_; + + 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_; @@ -897,6 +1033,16 @@ getReservoirConvergence(const double reportTime, msg += ") "; } + if (use_dsol_tol) { + msg += use_dp_tol ? " DP " : ""; + msg += use_ds_tol ? " DS " : ""; + msg += use_drs_tol ? " DRS " : ""; + msg += use_drv_tol ? " DRV " : ""; + } + + msg += " MBFLAG"; + msg += " CNVFLAG"; + OpmLog::debug(msg); } @@ -913,6 +1059,31 @@ getReservoirConvergence(const double reportTime, ss << std::setw(11) << CNV[compIdx]; } + 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); + } + + 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); diff --git a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp index 6f802c2f25..369f9bb844 100644 --- a/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp @@ -355,6 +355,26 @@ public: } } + void updateTUNINGDP(const TuningDp& tuning_dp) + { + // 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_) { + // 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 +469,19 @@ public: result = max_next_tstep > 0; } + 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); + + // 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}},