diff --git a/opm/simulators/flow/ExtraConvergenceOutputThread.cpp b/opm/simulators/flow/ExtraConvergenceOutputThread.cpp index e877a160956..fc99e695a73 100644 --- a/opm/simulators/flow/ExtraConvergenceOutputThread.cpp +++ b/opm/simulators/flow/ExtraConvergenceOutputThread.cpp @@ -71,6 +71,23 @@ namespace { }; } + // The trailing diagnostic columns, each with the longest token it can emit. + // Sizing them by their headings alone is not enough -- "GROUP+NETWORK" is + // wider than "Gate". "WellStatus" appends free-form text after its token, + // which is why it must stay last. + auto diagnosticHeaders() noexcept + { + using namespace std::literals::string_view_literals; + + return std::array { + std::pair { "CnvRelax"sv, "FINALIT"sv }, // to_string(CnvRelaxSource) + std::pair { "CnvTolUsed"sv, ""sv }, // a number, never wider + std::pair { "OscSource"sv, "WELLCTRL"sv }, // to_string(OscillationSource) + std::pair { "Gate"sv, "GROUP+NETWORK"sv }, + std::pair { "WellStatus"sv, "CONV"sv }, + }; + } + template int maxHeaderSize(const HeaderSequence& headers) { @@ -126,17 +143,27 @@ namespace { } } - const auto& metrics = firstRequest.reports.front().reservoirConvergence(); - const auto headerSize = maxColHeaderSize(minColSize, getPhaseName, metrics); + const auto& metrics = firstRequest.reports.front().reservoirConvergence(); + + auto headerSize = maxColHeaderSize(minColSize, getPhaseName, metrics); + for (const auto& [heading, longestValue] : diagnosticHeaders()) { + headerSize = std::max({ headerSize, + static_cast(heading.size() + 1), + static_cast(longestValue.size() + 1) }); + } for (const auto& metric : metrics) { os << std::right << std::setw(headerSize) << formatMetricColumn(getPhaseName, metric); } + for (const auto& [heading, _] : diagnosticHeaders()) { + os << std::right << std::setw(headerSize) << heading; + } + // Note: Newline character intentionally placed in separate output // request to not influence right-justification of column header. - os << std::right << std::setw(headerSize) << "WellStatus" << '\n'; + os << '\n'; return { minColSize, headerSize }; } @@ -224,6 +251,46 @@ namespace { } } + void writeCnvRelaxation(std::ostream& os, + const int colSize, + const Opm::ConvergenceReport& report) + { + // Which condition granted the relaxed CNV tolerance, and the tolerance actually + // applied. Without this the accepted state's quality is invisible: the same run + // may accept some steps at tolerance-cnv and others at its relaxed twin. + os << std::right << std::setw(colSize) << to_string(report.cnvRelaxSource()); + + const auto tol = report.cnvToleranceApplied(); + if (tol > 0.0) { + os << std::right << std::setw(colSize) << tol; + } + else { + os << std::right << std::setw(colSize) << '-'; + } + } + + void writeOscillationSource(std::ostream& os, + const int colSize, + const Opm::ConvergenceReport& report) + { + os << std::right << std::setw(colSize) << to_string(report.oscillationSource()); + } + + void writeConvergenceGate(std::ostream& os, + const int colSize, + const Opm::ConvergenceReport& report) + { + // Which secondary gate, if any, kept this iteration from counting as + // converged. These gates fire when the residual metrics themselves + // are already satisfied, so they identify iterations spent purely on + // well/group control changes or network balancing. + const auto gate = report.wellGroupTargetsViolated() + ? (report.networkNeedsMoreBalancing() ? "GROUP+NETWORK" : "GROUP") + : (report.networkNeedsMoreBalancing() ? "NETWORK" : "NONE"); + + os << std::right << std::setw(colSize) << gate; + } + void writeConvergenceRequest(std::ostream& os, const Opm::ConvergenceOutputThread::ConvertToTimeUnits& convertTime, const int firstColSize, @@ -243,6 +310,9 @@ namespace { writePenaltyCount(os, firstColSize, report); writeReservoirConvergence(os, colSize, report); + writeCnvRelaxation(os, colSize, report); + writeOscillationSource(os, colSize, report); + writeConvergenceGate(os, colSize, report); writeWellConvergence(os, colSize, report); os << '\n'; diff --git a/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp b/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp index e39e61b4df0..a315b0603e7 100644 --- a/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp +++ b/opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp @@ -173,6 +173,11 @@ initialLinearization(SimulatorReportSingle& report, auto convrep = getConvergence(timer, maxIter, residual_norms); report.converged = convrep.converged() && this->simulator_.problem().iterationContext().iteration() >= minIter; + if (report.converged && + convrep.cnvRelaxSource() != ConvergenceReport::CnvRelaxSource::None) + { + ++report.relaxed_cnv_acceptances; + } ConvergenceReport::Severity severity = convrep.severityOfWorstFailure(); this->convergence_reports_.back().report.push_back(std::move(convrep)); @@ -290,8 +295,21 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer, if (isOscillate) { this->current_relaxation_ -= nonlinear_solver.relaxIncrement(); this->current_relaxation_ = std::max(this->current_relaxation_, nonlinear_solver.relaxMax()); + // The detector reads reservoir residual history and the response damps the + // reservoir update, but the oscillation may originate in the well/group + // control layer. Record what was unsatisfied here so the two can be told + // apart. + auto source = ConvergenceReport::OscillationSource::NotDetected; + if (!this->convergence_reports_.empty() && + !this->convergence_reports_.back().report.empty()) + { + auto& convrep = this->convergence_reports_.back().report.back(); + source = classifyOscillationSource(convrep); + convrep.setOscillationSource(source); + } if (this->terminalOutputEnabled()) { - OpmLog::info(" Oscillating behavior detected: Relaxation set to " + OpmLog::info(" Oscillating behavior detected (" + to_string(source) + + "): Relaxation set to " + std::to_string(this->current_relaxation_)); } } @@ -816,6 +834,18 @@ getReservoirConvergence(const double reportTime, const auto tol_cnv = use_relaxed_cnv ? tolerance_cnv_relaxed : this->param_.tolerance_cnv_; const auto tol_mb = use_relaxed_mb ? this->param_.tolerance_mb_relaxed_ : this->param_.tolerance_mb_; + + // Record which condition granted the relaxation, so the accepted state's quality is + // visible in the output. Order mirrors the code's own precedence. + { + using RS = ConvergenceReport::CnvRelaxSource; + const auto source = relax_dsol_cnv ? RS::SolChange + : relax_pv_fraction_cnv ? RS::PvFraction + : relax_final_iteration_cnv ? RS::FinalIter + : relax_iter_cnv ? RS::IterCount + : RS::None; + report.setCnvRelaxation(source, static_cast(tol_cnv)); + } const auto tol_cnv_energy = use_relaxed_cnv ? this->param_.tolerance_cnv_energy_relaxed_ : this->param_.tolerance_cnv_energy_; const auto tol_eb = use_relaxed_mb ? this->param_.tolerance_energy_balance_relaxed_ : this->param_.tolerance_energy_balance_; diff --git a/opm/simulators/timestepping/ConvergenceReport.cpp b/opm/simulators/timestepping/ConvergenceReport.cpp index 9bf81495b79..9f3b583a33d 100644 --- a/opm/simulators/timestepping/ConvergenceReport.cpp +++ b/opm/simulators/timestepping/ConvergenceReport.cpp @@ -19,6 +19,7 @@ #include +#include #include #include #include @@ -107,4 +108,57 @@ namespace Opm pc.nonConverged, pc.distanceDecay, pc.largeWellResiduals, pc.total()); } + std::string to_string(const ConvergenceReport::CnvRelaxSource s) + { + using S = ConvergenceReport::CnvRelaxSource; + switch (s) { + case S::None: return "NONE"; + case S::PvFraction: return "PVFRAC"; + case S::SolChange: return "DSOL"; + case S::FinalIter: return "FINALIT"; + case S::IterCount: return "ITER"; + } + return "?"; + } + + std::string to_string(const ConvergenceReport::OscillationSource s) + { + using S = ConvergenceReport::OscillationSource; + switch (s) { + case S::NotDetected: return "NONE"; + case S::Reservoir: return "RES"; + case S::WellControl: return "WELLCTRL"; + case S::Mixed: return "MIXED"; + } + return "?"; + } + + ConvergenceReport::OscillationSource + classifyOscillationSource(const ConvergenceReport& report) + { + using WT = ConvergenceReport::WellFailure::Type; + + const auto& wellFailures = report.wellFailures(); + const bool controlUnsatisfied = + std::any_of(wellFailures.begin(), wellFailures.end(), + [](const ConvergenceReport::WellFailure& wf) + { + return (wf.type() == WT::ControlBHP) + || (wf.type() == WT::ControlTHP) + || (wf.type() == WT::ControlRate); + }) + || report.wellGroupTargetsViolated() + || report.networkNeedsMoreBalancing(); + + const bool reservoirUnsatisfied = report.reservoirFailed(); + + if (controlUnsatisfied && reservoirUnsatisfied) { + return ConvergenceReport::OscillationSource::Mixed; + } + if (controlUnsatisfied) { + return ConvergenceReport::OscillationSource::WellControl; + } + return ConvergenceReport::OscillationSource::Reservoir; + } + } // namespace Opm diff --git a/opm/simulators/timestepping/ConvergenceReport.hpp b/opm/simulators/timestepping/ConvergenceReport.hpp index dedd9e940d5..bd456a3889e 100644 --- a/opm/simulators/timestepping/ConvergenceReport.hpp +++ b/opm/simulators/timestepping/ConvergenceReport.hpp @@ -240,6 +240,29 @@ namespace Opm std::string well_name_ {}; }; + /// Which condition, if any, granted the relaxed CNV tolerance this iteration. + /// Ordered by precedence in the code so a single value can be reported even + /// though several conditions may hold at once. + enum struct CnvRelaxSource { + None, //!< strict tolerance applied + PvFraction, //!< non-strict pore volume below relaxed_max_pv_fraction + SolChange, //!< TUNINGDP-style solution-change targets met (disables CNV) + FinalIter, //!< min_strict_cnv_iter < 0 and iteration == NewtonMaxIterations + IterCount, //!< min_strict_cnv_iter >= 0 and iteration >= it + }; + + /// What was unsatisfied at the iteration where the Newton solver detected an + /// oscillation. An association, not a proof of cause. Note the base rate when + /// reading this: a non-final iteration almost always has unsatisfied reservoir + /// residuals, so Reservoir is the common label and the discriminating signal is + /// how often the control layer *also* participates (Mixed + WellControl). + enum struct OscillationSource { + NotDetected, + Reservoir, //!< reservoir residuals only + WellControl, //!< well control equations or a group/network gate only + Mixed, //!< both + }; + // ----------- Mutating member functions ----------- ConvergenceReport() @@ -305,6 +328,21 @@ namespace Opm this->eligiblePoreVolume_ = eligiblePoreVolume; } + //! \brief Record that the Newton solver damped this iteration for oscillation, + //! together with what was unsatisfied at that point. + void setOscillationSource(const OscillationSource source) + { + this->oscillationSource_ = source; + } + + //! \brief Record which condition granted the relaxed CNV tolerance and the + //! tolerance actually applied, so an accepted step's quality is visible. + void setCnvRelaxation(const CnvRelaxSource source, const double toleranceApplied) + { + this->cnvRelaxSource_ = source; + this->cnvToleranceApplied_ = toleranceApplied; + } + ConvergenceReport& operator+=(const ConvergenceReport& other) { reportTime_ = std::max(reportTime_, other.reportTime_); @@ -331,6 +369,13 @@ namespace Opm this->eligiblePoreVolume_ = other.eligiblePoreVolume_; } + // Same reasoning as the pore-volume split: only the reservoir report sets + // these, so take 'other's values whenever it carries them. + if (other.cnvToleranceApplied_ > 0.0) { + this->cnvRelaxSource_ = other.cnvRelaxSource_; + this->cnvToleranceApplied_ = other.cnvToleranceApplied_; + } + return *this; } @@ -351,6 +396,24 @@ namespace Opm return this->cnvPvSplit_; } + //! \brief What was unsatisfied when oscillation damping was triggered, if it was. + OscillationSource oscillationSource() const + { + return this->oscillationSource_; + } + + //! \brief Which condition granted the relaxed CNV tolerance this iteration. + CnvRelaxSource cnvRelaxSource() const + { + return this->cnvRelaxSource_; + } + + //! \brief The CNV tolerance actually applied (0 if not recorded). + double cnvToleranceApplied() const + { + return this->cnvToleranceApplied_; + } + bool converged() const { return (status_ == AllGood) @@ -358,6 +421,21 @@ namespace Opm && !network_needs_more_balancing_force_another_newton_iteration_; } + //! \brief Whether well/group control targets were violated (a control + //! or target changed this iteration), forcing another Newton iteration + //! even when all residual metrics are converged. + bool wellGroupTargetsViolated() const + { + return wellGroupTargetsViolated_; + } + + //! \brief Whether the network balance forced another Newton iteration + //! even when all residual metrics are converged. + bool networkNeedsMoreBalancing() const + { + return network_needs_more_balancing_force_another_newton_iteration_; + } + bool reservoirFailed() const { return status_ & ReservoirFailed; @@ -437,6 +515,9 @@ namespace Opm serializer(this->network_needs_more_balancing_force_another_newton_iteration_); serializer(this->cnvPvSplit_); serializer(this->eligiblePoreVolume_); + serializer(this->cnvRelaxSource_); + serializer(this->cnvToleranceApplied_); + serializer(this->oscillationSource_); serializer(this->penaltyCard_); } @@ -454,6 +535,9 @@ namespace Opm bool network_needs_more_balancing_force_another_newton_iteration_; CnvPvSplit cnvPvSplit_{}; double eligiblePoreVolume_{}; + CnvRelaxSource cnvRelaxSource_{CnvRelaxSource::None}; + double cnvToleranceApplied_{}; + OscillationSource oscillationSource_{OscillationSource::NotDetected}; PenaltyCard penaltyCard_; }; @@ -474,6 +558,13 @@ namespace Opm std::string to_string(const ConvergenceReport::PenaltyCard& pc); + std::string to_string(const ConvergenceReport::CnvRelaxSource s); + + std::string to_string(const ConvergenceReport::OscillationSource s); + + /// Classify what was unsatisfied in a report, for oscillation-source reporting. + ConvergenceReport::OscillationSource classifyOscillationSource(const ConvergenceReport& report); + } // namespace Opm diff --git a/opm/simulators/timestepping/SimulatorReport.cpp b/opm/simulators/timestepping/SimulatorReport.cpp index 4dce35e4fd9..223fe0b3ba7 100644 --- a/opm/simulators/timestepping/SimulatorReport.cpp +++ b/opm/simulators/timestepping/SimulatorReport.cpp @@ -34,9 +34,9 @@ namespace Opm { return SimulatorReportSingle{1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, - 13, 14, 15, 16, 17, 18, - true, false, false, 19, 20.0, 21.0, - 22, 23, 24, 25, 26, 27, 28, 29}; + 13, 14, 15, 16, 17, 18, 19, + true, false, false, 20, 21.0, 22.0, + 23, 24, 25, 26, 27, 28, 29, 30}; } bool SimulatorReportSingle::operator==(const SimulatorReportSingle& rhs) const @@ -56,6 +56,7 @@ namespace Opm this->total_well_iterations == rhs.total_well_iterations && this->total_linearizations == rhs.total_linearizations && this->total_newton_iterations == rhs.total_newton_iterations && + this->relaxed_cnv_acceptances == rhs.relaxed_cnv_acceptances && this->total_linear_iterations == rhs.total_linear_iterations && this->min_linear_iterations == rhs.min_linear_iterations && this->max_linear_iterations == rhs.max_linear_iterations && @@ -92,6 +93,7 @@ namespace Opm total_well_iterations += sr.total_well_iterations; total_linearizations += sr.total_linearizations; total_newton_iterations += sr.total_newton_iterations; + relaxed_cnv_acceptances += sr.relaxed_cnv_acceptances; total_linear_iterations += sr.total_linear_iterations; if (sr.total_linear_iterations > 0) { min_linear_iterations = std::min(min_linear_iterations, sr.total_linear_iterations); @@ -237,6 +239,13 @@ namespace Opm 100.0*failureReport->total_linear_iterations/noZero(n)); } os << std::endl; + + // A step accepted only under the relaxed CNV tolerance carries a local + // volumetric error up to tolerance-cnv-relaxed rather than tolerance-cnv. + if (relaxed_cnv_acceptances > 0) { + os << fmt::format("Relaxed CNV acceptances: {:7}", relaxed_cnv_acceptances) + << std::endl; + } } diff --git a/opm/simulators/timestepping/SimulatorReport.hpp b/opm/simulators/timestepping/SimulatorReport.hpp index 5d473f2292a..cd56a1a3f6e 100644 --- a/opm/simulators/timestepping/SimulatorReport.hpp +++ b/opm/simulators/timestepping/SimulatorReport.hpp @@ -51,6 +51,8 @@ namespace Opm unsigned int total_linear_iterations = 0; unsigned int min_linear_iterations = std::numeric_limits::max(); unsigned int max_linear_iterations = 0; + // Accepted substeps whose CNV criterion was met only under a relaxed tolerance. + unsigned int relaxed_cnv_acceptances = 0; bool converged = false; bool time_step_rejected = false; @@ -101,6 +103,7 @@ namespace Opm serializer(total_linear_iterations); serializer(min_linear_iterations); serializer(max_linear_iterations); + serializer(relaxed_cnv_acceptances); serializer(converged); serializer(time_step_rejected); serializer(well_group_control_changed); diff --git a/tests/test_tuning_TRGMBE.cpp b/tests/test_tuning_TRGMBE.cpp index 90abe554784..8f5b1cbedd1 100644 --- a/tests/test_tuning_TRGMBE.cpp +++ b/tests/test_tuning_TRGMBE.cpp @@ -83,7 +83,7 @@ struct ColumnData while (iss >> colname) { column_names.push_back(colname); raw_columns.emplace_back(colname); - columns[colname] = &(raw_columns.back()); + columns[colname] = raw_columns.size() - 1; } const int num_columns = column_names.size(); @@ -104,14 +104,16 @@ struct ColumnData } // Get data vectors of different types - std::vector get_dvector(const std::string& colname) const { return columns.at(colname)->dvalues(); } - std::vector get_ivector(const std::string& colname) const { return columns.at(colname)->ivalues(); } + std::vector get_dvector(const std::string& colname) const { return raw_columns[columns.at(colname)].dvalues(); } + std::vector get_ivector(const std::string& colname) const { return raw_columns[columns.at(colname)].ivalues(); } // Default is to return double values - std::vector operator[](const std::string& colname) const { return columns.at(colname)->dvalues(); } + std::vector operator[](const std::string& colname) const { return raw_columns[columns.at(colname)].dvalues(); } std::vector column_names; std::vector raw_columns; - std::map columns; + // Index, not pointer: raw_columns reallocates as soon as the file has more + // columns than the reserve, which dangles every pointer stored before it. + std::map columns; }; BOOST_AUTO_TEST_CASE(CheckMassBalanceWithinTRGMBE)