From c237ab491fbd4df7ca985ae3dd928b53cfafd521 Mon Sep 17 00:00:00 2001 From: David Landa Marban Date: Wed, 16 Jul 2025 21:19:04 +0200 Subject: [PATCH] Support for CONV_NEW when CONV in RPTRST --- opm/simulators/flow/BlackoilModel.hpp | 8 ++++ opm/simulators/flow/BlackoilModel_impl.hpp | 50 ++++++++++++++++++++++ opm/simulators/flow/RSTConv.cpp | 38 ++++++++++++++++ opm/simulators/flow/RSTConv.hpp | 16 ++++++- tests/test_rstconv.cpp | 11 +++++ 5 files changed, 122 insertions(+), 1 deletion(-) diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index d65994e45b..25289fced8 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -213,6 +213,14 @@ public: std::pair, std::vector> 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. + void convergencePerCell(const std::vector& B_avg, + const double tol_cnv, + const double dt, + const double tol_cnv_energy, + const int iteration); + void updateTUNING(const Tuning& tuning); ConvergenceReport diff --git a/opm/simulators/flow/BlackoilModel_impl.hpp b/opm/simulators/flow/BlackoilModel_impl.hpp index d1b80278a1..139c991b35 100644 --- a/opm/simulators/flow/BlackoilModel_impl.hpp +++ b/opm/simulators/flow/BlackoilModel_impl.hpp @@ -894,6 +894,9 @@ getReservoirConvergence(const double reportTime, } } + // Compute the Newton convergence per cell. + this->convergencePerCell(B_avg, dt, tol_cnv, tol_cnv_energy, iteration); + // Output of residuals. if (this->terminal_output_) { // Only rank 0 does print to std::cout @@ -936,6 +939,53 @@ getReservoirConvergence(const double reportTime, return report; } +template +void +BlackoilModel:: +convergencePerCell(const std::vector& B_avg, + const double dt, + const double tol_cnv, + const double tol_cnv_energy, + const int iteration) +{ + auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv(); + if (!rst_conv.hasConv()) { + return; + } + + if (iteration == 0) { + rst_conv.prepareConv(); + } + + const auto& residual = simulator_.model().linearizer().residual(); + const auto& gridView = this->simulator().gridView(); + const IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid()); + ElementContext elemCtx(this->simulator()); + std::vector convNewt(residual.size(), 0); + OPM_BEGIN_PARALLEL_TRY_CATCH(); + unsigned idx = 0; + const int numComp = B_avg.size(); + for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { + elemCtx.updatePrimaryStencil(elem); + + const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto pvValue = simulator_.problem().referencePorosity(cell_idx, /*timeIdx=*/0) * + simulator_.model().dofTotalVolume(cell_idx); + for (int compIdx = 0; compIdx < numComp; ++compIdx) { + const auto tol = (has_energy_ && compIdx == contiEnergyEqIdx) ? tol_cnv_energy : tol_cnv; + const Scalar cnv = std::abs(B_avg[compIdx] * residual[cell_idx][compIdx]) * dt / pvValue; + if (std::isnan(cnv) || cnv > maxResidualAllowed() || cnv < 0.0 || cnv > tol) { + convNewt[idx] = 1; + break; + } + } + ++idx; + } + OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::convergencePerCell() failed: ", + this->grid_.comm()); + rst_conv.updateNewton(convNewt); +} + template ConvergenceReport BlackoilModel:: diff --git a/opm/simulators/flow/RSTConv.cpp b/opm/simulators/flow/RSTConv.cpp index 45f5d9ec21..10f25321f6 100644 --- a/opm/simulators/flow/RSTConv.cpp +++ b/opm/simulators/flow/RSTConv.cpp @@ -42,6 +42,7 @@ void RSTConv::init(const std::size_t numCells, if (kw == rst_config.keywords.end()) { N_ = 0; cnv_X_.clear(); + conv_new_.clear(); return; } @@ -54,6 +55,7 @@ void RSTConv::init(const std::size_t numCells, cnv_X_[i].resize(numCells); } } + conv_new_.resize(numCells, 1); } void RSTConv::outputRestart(data::Solution& sol) @@ -69,6 +71,42 @@ void RSTConv::outputRestart(data::Solution& sol) } ++i; }); + sol.insert("CONV_NEW", std::move(conv_new_), data::TargetType::RESTART_SOLUTION); + } +} + +bool RSTConv::hasConv() const +{ + return !conv_new_.empty(); +} + +void RSTConv::prepareConv() +{ + std::fill(conv_new_.begin(), conv_new_.end(), 1); +} + +void RSTConv::updateNewton(const std::vector& convNewt) +{ + const std::size_t numGloCells = conv_new_.size(); + + if (comm_.size() == 1) { + for (std::size_t n = 0; n < numGloCells; ++n) { + conv_new_[n] += convNewt[n]; + } + return; + } + + std::vector values(numGloCells, 0); + const std::size_t numLocCells = convNewt.size(); + + for (std::size_t i = 0; i < numLocCells; ++i) { + values[globalCell_(i)] = convNewt[i]; + } + + comm_.sum(values.data(), numGloCells); + + for (std::size_t n = 0; n < numGloCells; ++n) { + conv_new_[n] += values[n]; } } diff --git a/opm/simulators/flow/RSTConv.hpp b/opm/simulators/flow/RSTConv.hpp index 4e188b6701..a169d7f3c1 100644 --- a/opm/simulators/flow/RSTConv.hpp +++ b/opm/simulators/flow/RSTConv.hpp @@ -60,10 +60,23 @@ public: template void update(const ResidualVector& resid); - //! \brief Obtain a const-ref to the accumulated data. + //! \brief Adds the CONV_NEW output. + void updateNewton(const std::vector& convNewt); + + //! \brief Obtain a const-ref to the accumulated data (only used for unit testing). const std::vector>& getData() const { return cnv_X_; } + //! \brief Obtain a const-ref to the CONV_NEW (only used for unit testing). + const std::vector& getConvNew() const + { return conv_new_; } + + //! \brief Check if CONV is required. + bool hasConv() const; + + //! \brief CONV_NEW is restarted at iteration 0. + void prepareConv(); + private: //! \brief Gathers and accumulates results to the CONV arrays. //! \param lIdx Vector of local indices (N first is used) @@ -77,6 +90,7 @@ private: LocalToGlobalCellFunc globalCell_; //!< Global cell indices Parallel::Communication comm_; //!< Communicator std::vector> cnv_X_{}; //!< Counts of worst cells for RPTRST CONV + std::vector conv_new_{}; //!< Number of Newtons by each cell std::array compIdx_{}; //!< Component indices int N_ = 0; //!< Number of cells to consider }; diff --git a/tests/test_rstconv.cpp b/tests/test_rstconv.cpp index e8028818da..9f09d3c5ac 100644 --- a/tests/test_rstconv.cpp +++ b/tests/test_rstconv.cpp @@ -176,6 +176,17 @@ BOOST_AUTO_TEST_CASE(RstConvTest) } } } + + std::vector conv_new(10, 0); + cnv.updateNewton(conv_new); + std::vector conv_new_ite1(10, 1); + cnv.updateNewton(conv_new_ite1); + if (cc.rank() == 0) { + BOOST_CHECK_EQUAL(cnv.getConvNew().size(), 10*cc.size()); + for (int i = 0; i < cnv.getConvNew().size(); ++i) { + BOOST_CHECK_EQUAL(cnv.getConvNew()[i], 2); + } + } #if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 67 } #endif