From 3c63a7aa6d7ee780a7b5d595416221f2a5a64a7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 8 Dec 2022 17:16:34 +0100 Subject: [PATCH] Capture Timestep's Non-Linear Convergence History This enables outputting a formatted record of the limiting MB and CNV quantities as time and non-linear iterations progress. This, in turn, is intended for diagnostic and analysis purposes and will not be output unless specifically requested. In particular, add a new type, ConvergenceReport::ReservoirConvergenceMetric which captures the convergence metric type (MB or CNV) along with the associate phase and numerical value of the convergence metric. We add a vector of these convergence metric objects as a new data member of the ConvergenceReport. Finally, foreshadowing the intended use case, also store the report time in the ConvergenceReport object. --- opm/simulators/flow/BlackoilModelEbos.hpp | 10 ++-- .../timestepping/ConvergenceReport.hpp | 45 +++++++++++++++- .../timestepping/gatherConvergenceReport.cpp | 52 +++++++++++++++++-- 3 files changed, 99 insertions(+), 8 deletions(-) diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index ed1fe2883..8e40d1321 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -888,7 +888,8 @@ namespace Opm { return grid_.comm().sum(errorPV); } - ConvergenceReport getReservoirConvergence(const double dt, + ConvergenceReport getReservoirConvergence(const double reportTime, + const double dt, const int iteration, std::vector& B_avg, std::vector& residual_norms) @@ -927,7 +928,7 @@ namespace Opm { } // Create convergence report. - ConvergenceReport report; + ConvergenceReport report{reportTime}; using CR = ConvergenceReport; for (int compIdx = 0; compIdx < numComp; ++compIdx) { double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] }; @@ -953,6 +954,7 @@ namespace Opm { } else if (res[ii] > tol[ii]) { report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx}); } + report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]); } } @@ -1003,7 +1005,9 @@ namespace Opm { { // Get convergence reports for reservoir and wells. std::vector B_avg(numEq, 0.0); - auto report = getReservoirConvergence(timer.currentStepLength(), iteration, B_avg, residual_norms); + auto report = getReservoirConvergence(timer.simulationTimeElapsed(), + timer.currentStepLength(), + iteration, B_avg, residual_norms); report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged()); return report; diff --git a/opm/simulators/timestepping/ConvergenceReport.hpp b/opm/simulators/timestepping/ConvergenceReport.hpp index 978b020f2..76562cd84 100644 --- a/opm/simulators/timestepping/ConvergenceReport.hpp +++ b/opm/simulators/timestepping/ConvergenceReport.hpp @@ -21,9 +21,11 @@ #ifndef OPM_CONVERGENCEREPORT_HEADER_INCLUDED #define OPM_CONVERGENCEREPORT_HEADER_INCLUDED +#include #include #include #include +#include #include namespace Opm @@ -61,6 +63,21 @@ namespace Opm Severity severity_; int phase_; }; + class ReservoirConvergenceMetric + { + public: + ReservoirConvergenceMetric(ReservoirFailure::Type t, int phase, double value) + : type_(t), phase_(phase), value_(value) + { + } + ReservoirFailure::Type type() const { return type_; } + int phase() const { return phase_; } + double value() const { return value_; } + private: + ReservoirFailure::Type type_; + int phase_; + double value_; + }; class WellFailure { public: @@ -83,7 +100,13 @@ namespace Opm // ----------- Mutating member functions ----------- ConvergenceReport() - : status_{AllGood} + : ConvergenceReport{0.0} + { + } + + explicit ConvergenceReport(const double reportTime) + : reportTime_{reportTime} + , status_{AllGood} , res_failures_{} , well_failures_{} , wellGroupTargetsViolated_(false) @@ -110,6 +133,12 @@ namespace Opm well_failures_.push_back(wf); } + template + void setReservoirConvergenceMetric(Args&&... args) + { + this->res_convergence_.emplace_back(std::forward(args)...); + } + void setWellGroupTargetsViolated(const bool wellGroupTargetsViolated) { wellGroupTargetsViolated_ = wellGroupTargetsViolated; @@ -117,9 +146,11 @@ namespace Opm ConvergenceReport& operator+=(const ConvergenceReport& other) { + reportTime_ = std::max(reportTime_, other.reportTime_); status_ = static_cast(status_ | other.status_); res_failures_.insert(res_failures_.end(), other.res_failures_.begin(), other.res_failures_.end()); well_failures_.insert(well_failures_.end(), other.well_failures_.begin(), other.well_failures_.end()); + res_convergence_.insert(res_convergence_.end(), other.res_convergence_.begin(), other.res_convergence_.end()); assert(reservoirFailed() != res_failures_.empty()); assert(wellFailed() != well_failures_.empty()); wellGroupTargetsViolated_ = (wellGroupTargetsViolated_ || other.wellGroupTargetsViolated_); @@ -128,6 +159,11 @@ namespace Opm // ----------- Const member functions (queries) ----------- + double reportTime() const + { + return reportTime_; + } + bool converged() const { return (status_ == AllGood) && !wellGroupTargetsViolated_; @@ -148,6 +184,11 @@ namespace Opm return res_failures_; } + const std::vector& reservoirConvergence() const + { + return res_convergence_; + } + const std::vector& wellFailures() const { return well_failures_; @@ -172,9 +213,11 @@ namespace Opm private: // ----------- Member variables ----------- + double reportTime_; Status status_; std::vector res_failures_; std::vector well_failures_; + std::vector res_convergence_; bool wellGroupTargetsViolated_; }; diff --git a/opm/simulators/timestepping/gatherConvergenceReport.cpp b/opm/simulators/timestepping/gatherConvergenceReport.cpp index f781ac239..514a76f6f 100644 --- a/opm/simulators/timestepping/gatherConvergenceReport.cpp +++ b/opm/simulators/timestepping/gatherConvergenceReport.cpp @@ -1,6 +1,6 @@ /* + Copyright 2018, 2022 Equinor ASA. Copyright 2018 SINTEF Digital, Mathematics and Cybernetics. - Copyright 2018 Equinor. This file is part of the Open Porous Media project (OPM). @@ -43,6 +43,18 @@ namespace MPI_Pack(&phase, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator); } + void packReservoirConvergenceMetric(const ConvergenceReport::ReservoirConvergenceMetric& m, + std::vector& buf, + int& offset, MPI_Comm mpi_communicator) + { + int type = static_cast(m.type()); + int phase = m.phase(); + double value = m.value(); + MPI_Pack(&type, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator); + MPI_Pack(&phase, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator); + MPI_Pack(&value, 1, MPI_DOUBLE, buf.data(), buf.size(), &offset, mpi_communicator); + } + void packWellFailure(const ConvergenceReport::WellFailure& f, std::vector& buf, int& offset, MPI_Comm mpi_communicator) @@ -65,12 +77,21 @@ namespace // Pack the data. // Status will not be packed, it is possible to deduce from the other data. // Reservoir failures. + double reportTime = local_report.reportTime(); + MPI_Pack(&reportTime, 1, MPI_DOUBLE, buf.data(), buf.size(), &offset, mpi_communicator); const auto rf = local_report.reservoirFailures(); int num_rf = rf.size(); MPI_Pack(&num_rf, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator); for (const auto& f : rf) { packReservoirFailure(f, buf, offset, mpi_communicator); } + // Reservoir convergence metrics. + const auto rm = local_report.reservoirConvergence(); + int num_rm = rm.size(); + MPI_Pack(&num_rm, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator); + for (const auto& m : rm) { + packReservoirConvergenceMetric(m, buf, offset, mpi_communicator); + } // Well failures. const auto wf = local_report.wellFailures(); int num_wf = wf.size(); @@ -84,13 +105,16 @@ namespace { int int_pack_size = 0; MPI_Pack_size(1, MPI_INT, mpi_communicator, &int_pack_size); + int double_pack_size = 0; + MPI_Pack_size(1, MPI_DOUBLE, mpi_communicator, &double_pack_size); const int num_rf = local_report.reservoirFailures().size(); + const int num_rm = local_report.reservoirConvergence().size(); const int num_wf = local_report.wellFailures().size(); int wellnames_length = 0; for (const auto& f : local_report.wellFailures()) { wellnames_length += (f.wellName().size() + 1); } - return (2 + 3*num_rf + 4*num_wf) * int_pack_size + wellnames_length; + return (3 + 3*num_rf + 2*num_rm + 4*num_wf)*int_pack_size + (1 + 1*num_rm)*double_pack_size + wellnames_length; } ConvergenceReport::ReservoirFailure unpackReservoirFailure(const std::vector& recv_buffer, int& offset, MPI_Comm mpi_communicator) @@ -107,6 +131,19 @@ namespace phase); } + ConvergenceReport::ReservoirConvergenceMetric + unpackReservoirConvergenceMetric(const std::vector& recv_buffer, int& offset, MPI_Comm mpi_communicator) + { + int type = -1; + int phase = -1; + double value = -1.0; + auto* data = const_cast(recv_buffer.data()); + MPI_Unpack(data, recv_buffer.size(), &offset, &type, 1, MPI_INT, mpi_communicator); + MPI_Unpack(data, recv_buffer.size(), &offset, &phase, 1, MPI_INT, mpi_communicator); + MPI_Unpack(data, recv_buffer.size(), &offset, &value, 1, MPI_DOUBLE, mpi_communicator); + return { static_cast(type), phase, value }; + } + ConvergenceReport::WellFailure unpackWellFailure(const std::vector& recv_buffer, int& offset, MPI_Comm mpi_communicator) { int type = -1; @@ -129,14 +166,21 @@ namespace ConvergenceReport unpackSingleConvergenceReport(const std::vector& recv_buffer, int& offset, MPI_Comm mpi_communicator) { - ConvergenceReport cr; - int num_rf = -1; auto* data = const_cast(recv_buffer.data()); + double reportTime{0.0}; + MPI_Unpack(data, recv_buffer.size(), &offset, &reportTime, 1, MPI_DOUBLE, mpi_communicator); + ConvergenceReport cr{reportTime}; + int num_rf = -1; MPI_Unpack(data, recv_buffer.size(), &offset, &num_rf, 1, MPI_INT, mpi_communicator); for (int rf = 0; rf < num_rf; ++rf) { ConvergenceReport::ReservoirFailure f = unpackReservoirFailure(recv_buffer, offset, mpi_communicator); cr.setReservoirFailed(f); } + int num_rm = -1; + MPI_Unpack(data, recv_buffer.size(), &offset, &num_rm, 1, MPI_INT, mpi_communicator); + for (int rm = 0; rm < num_rm; ++rm) { + cr.setReservoirConvergenceMetric(unpackReservoirConvergenceMetric(recv_buffer, offset, mpi_communicator)); + } int num_wf = -1; MPI_Unpack(data, recv_buffer.size(), &offset, &num_wf, 1, MPI_INT, mpi_communicator); for (int wf = 0; wf < num_wf; ++wf) {