From b445db2ab8c55dfa9bbba656dda458109a737bbb Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 14:24:41 +0200 Subject: [PATCH 01/22] Add a CPRW pressure stage to the system solver The system solver (--linear-solver=system_cpr) solves the coupled (reservoir, well) system [A C; B D] explicitly, but its pressure stage was still reservoir-only: stage 1 of SystemPreconditioner ran a plain CPR on A with add_wells = false, so the well unknowns never entered the coarse system. That is exactly the coupling CPRW exists to supply. Classic CPRW cannot be reused here. PressureBhpTransferPolicy calls back into the live well model from inside the preconditioner (addWellPressureEquations -> per-well extractCPRPressureMatrix, which needs WellState just to ask isPressureControlled). That is what ties WellModelAsLinearOperator to TypeTag and what keeps CPRW out of NLDD and the GPU path. Instead the coarse system is assembled from the B/C/D blocks the outer layer already extracts, once, via addBCDMatrix. Two pieces of plain data are added to carry what is still missing: - WellDofLayout: a prefix sum over the per-well D dimensions, so the merged well block rows can be mapped back to wells (one row per standard well, one per segment for a multisegment well), plus the index of the pressure-like unknown inside a well block. - the well half of the weights. The weights calculator was already std::function; it now fills w[_1] as well as w[_0]. ISTLSolverSystem computes those weights (quasi-IMPES from the diagonal D block by default, 'unit' for debugging) and can later obtain them from the well model without the core changing. Both are produced in ISTLSolverSystem. Below that point the pressure stage reads nothing but sparse matrices, weights and integers -- no part of the well model is visible. The preconditioner factory previously discarded w[_1]; it now passes the whole SystemVector through and SystemPreconditioner derives the reservoir-only calculator for its sub-solvers. SystemCprwPressureStage builds the (Nres + nWells) scalar system as R*S*P with R = blockdiag(w0^T ; sum over the well's block rows of w1^T) and P = blockdiag(e_p ; e_q on the well's top block row), reusing Details::CoarseOperatorType and extendCommunicatorWithWells from PressureBhpTransferPolicy. One coarse unknown per well, as in classic CPRW; one per segment is left as a separate question. Two details differ from the classic policy: the well residual is really restricted, and the coarse well correction is prolonged back rather than discarded. Both are selectable through preconditioner.well_transfer, whose 'classic' value reproduces the classic formulation so that the two differ only in numerics. On SPE1 that lands on the classic cprw iteration count exactly. Selected with --linear-solver=system_cprw, or by setting preconditioner.reservoir_solver.preconditioner.add_wells in a JSON configuration. Linear iterations, serial: deck cprw system_cpr system_cprw SPE1 443 571 429 SPE9_CP 439 682 470 BASE2_MSW_HFA 30 82 52 and on SPE9_CP_SHORT the count stays flat under decomposition (110/107/105/106 for 1/2/4/8 ranks) where system_cpr drifts up (177/178/179/180). Also reject an approximate (Krylov) well solver combined with a non-flexible outer solver: such a well solve stops on a tolerance, so the preconditioner varies between applications and bicgstab or plain gmres are no longer valid. flexgmres is required. Co-Authored-By: Claude Opus 5 (cherry picked from commit 6dd1454876a9e19f3db046675fb418af56d98a0a) --- CMakeLists_files.cmake | 6 + .../linalg/FlowLinearSolverParameters.cpp | 1 + .../linalg/ISTLSolverRuntimeOptionProxy.hpp | 2 +- opm/simulators/linalg/setupPropertyTree.cpp | 58 ++- .../linalg/system/ISTLSolverSystem.hpp | 98 +++- .../linalg/system/SystemCprwPressureStage.hpp | 470 +++++++++++++++++ .../linalg/system/SystemPreconditioner.hpp | 108 +++- .../system/SystemPreconditionerFactory.cpp | 24 +- opm/simulators/linalg/system/SystemTypes.hpp | 59 +++ tests/options_system_cprw_approx_wells.json | 77 +++ ...ns_system_cprw_approx_wells_bad_outer.json | 76 +++ tests/options_system_cprw_complete.json | 72 +++ ...ions_system_cprw_missing_coarsesolver.json | 40 ++ tests/test_SystemCprwPressureStage.cpp | 486 ++++++++++++++++++ tests/test_setuppropertytree.cpp | 52 ++ 15 files changed, 1590 insertions(+), 39 deletions(-) create mode 100644 opm/simulators/linalg/system/SystemCprwPressureStage.hpp create mode 100644 tests/options_system_cprw_approx_wells.json create mode 100644 tests/options_system_cprw_approx_wells_bad_outer.json create mode 100644 tests/options_system_cprw_complete.json create mode 100644 tests/options_system_cprw_missing_coarsesolver.json create mode 100644 tests/test_SystemCprwPressureStage.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f24ae5d6a6a..1531b7adb68 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -516,6 +516,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_tpsa_localresidual.cpp tests/test_tpsa_primaryvariables.cpp tests/test_vfpproperties.cpp + tests/test_SystemCprwPressureStage.cpp tests/test_WellMatrixMerger.cpp tests/test_WaterSatfuncConsistencyChecks.cpp tests/test_wellmodel.cpp @@ -691,6 +692,10 @@ list (APPEND TEST_DATA_FILES tests/options_system_cpr_missing_smoother.json tests/options_system_cpr_missing_well.json tests/options_system_cpr_res_precond_not_cpr.json + tests/options_system_cprw_approx_wells.json + tests/options_system_cprw_approx_wells_bad_outer.json + tests/options_system_cprw_complete.json + tests/options_system_cprw_missing_coarsesolver.json tests/GCONSUMP.DATA tests/GCONSUMP_COMPLEX.DATA tests/GROUP_HIGHER_CONSTRAINTS.DATA @@ -1123,6 +1128,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp opm/simulators/linalg/Preconditioner2InverseOperator.hpp opm/simulators/linalg/system/MultiComm.hpp + opm/simulators/linalg/system/SystemCprwPressureStage.hpp opm/simulators/linalg/system/SystemPreconditioner.hpp opm/simulators/linalg/system/SystemPreconditionerFactory.hpp opm/simulators/linalg/system/SystemTypes.hpp diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.cpp b/opm/simulators/linalg/FlowLinearSolverParameters.cpp index 629ab7d4ef5..1e93160f58f 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.cpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.cpp @@ -135,6 +135,7 @@ void FlowLinearSolverParameters::registerParameters() ("Scale linear system according to equation scale and primary variable types"); Parameters::Register ("Configuration of solver. Valid options are: cprw (default), system_cpr (CPU-only), " + "system_cprw (CPU-only, system_cpr with the wells in the pressure stage), " "ilu0, dilu, cpr (an alias for cprw), cpr_quasiimpes, " "cpr_trueimpes, cpr_trueimpesanalytic, amg or hybrid (experimental). " "Alternatively, you can request a configuration to be read from a " diff --git a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp index 40a0e57c386..3c6cbe20e41 100644 --- a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp +++ b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp @@ -155,7 +155,7 @@ class ISTLSolverRuntimeOptionProxy : public AbstractISTLSolver(); - bool useSystemCpr = (linSolverConf == "system_cpr"); + bool useSystemCpr = (linSolverConf == "system_cpr") || (linSolverConf == "system_cprw"); if (!useSystemCpr && linSolverConf.size() > 5 && linSolverConf.ends_with(".json") && std::filesystem::exists(linSolverConf)) { diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index b406e32ec30..becfe8135ea 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -232,7 +232,9 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters } // System CPR configuration (coupled reservoir-well system solver). - if (conf == "system_cpr") { + // system_cprw differs only in that its pressure stage carries the well + // unknowns, exactly as cprw does relative to cpr. + if ((conf == "system_cpr") || (conf == "system_cprw")) { if (!linearSolverMaxIterSet) { p.linear_solver_maxiter_ = 20; } @@ -514,9 +516,10 @@ setupUMFPack([[maybe_unused]] const std::string& conf, const FlowLinearSolverPar PropertyTree -setupSystemCPR([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p) +setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) { using namespace std::string_literals; + const bool add_wells = (conf == "system_cprw"); PropertyTree prm; // Outer solver @@ -527,6 +530,16 @@ setupSystemCPR([[maybe_unused]] const std::string& conf, const FlowLinearSolverP // Top-level preconditioner: system_cpr prm.put("preconditioner.type", "system_cpr"s); + // How the well equations are contracted to the one coarse unknown each + // well carries in the CPRW pressure system. Only read when add_wells. + prm.put("preconditioner.well_weight_type", "quasiimpes"s); + // How the well unknowns take part in the pressure-stage transfer: + // full - restrict the well residual, prolong the bhp correction + // no_prolongation - restrict, but discard the bhp correction + // classic - neither, i.e. the classic cprw formulation, so that + // the only remaining difference is numerics + // Only read when add_wells. + prm.put("preconditioner.well_transfer", "full"s); // --- Reservoir smoother --- prm.put("preconditioner.reservoir_smoother.maxiter", 1); @@ -548,7 +561,10 @@ setupSystemCPR([[maybe_unused]] const std::string& conf, const FlowLinearSolverP prm.put("preconditioner.reservoir_solver.preconditioner.type", "cpr"s); prm.put("preconditioner.reservoir_solver.preconditioner.relaxation", 1.0); prm.put("preconditioner.reservoir_solver.preconditioner.use_well_weights", "false"s); - prm.put("preconditioner.reservoir_solver.preconditioner.add_wells", "false"s); + // add_wells promotes the pressure stage from reservoir-only CPR to CPRW + // over the full (reservoir, well) system. + prm.put("preconditioner.reservoir_solver.preconditioner.add_wells", + add_wells ? "true"s : "false"s); prm.put("preconditioner.reservoir_solver.preconditioner.weight_type", "trueimpes"s); prm.put("preconditioner.reservoir_solver.preconditioner.pre_smooth", 0); prm.put("preconditioner.reservoir_solver.preconditioner.post_smooth", 0); @@ -593,6 +609,42 @@ void validateSystemCPRTree(const PropertyTree& prm) "In system_cpr configuration, the reservoir_solver must use the CPR preconditioner " "(preconditioner.reservoir_solver.preconditioner.type = 'cpr')."); } + // With add_wells the pressure stage is assembled and solved directly by + // the system preconditioner, which takes its solver settings from the + // coarsesolver sub-tree rather than from the reservoir_solver wrapper. + if (reservoir_solver->get("preconditioner.add_wells", false) + && !reservoir_solver->get_child_optional("preconditioner.coarsesolver").has_value()) { + OPM_THROW(std::invalid_argument, + "In system_cpr configuration with " + "preconditioner.reservoir_solver.preconditioner.add_wells = true, the " + "'preconditioner.reservoir_solver.preconditioner.coarsesolver' sub-tree is " + "required: it configures the solver for the CPRW pressure system."); + } + } + + // A Krylov well solver stops on a tolerance, so it performs a different + // number of inner iterations for each right-hand side. That makes the whole + // system preconditioner non-stationary, which Krylov methods with short + // recurrences (bicgstab, cg) and standard GMRES are not allowed to use: they + // assume a fixed preconditioning operator. The outer solver has to be a + // flexible one. + auto well_solver = prm.get_child_optional("preconditioner.well_solver"); + if (well_solver) { + const auto inner = well_solver->get("solver", "bicgstab"); + const bool inner_is_krylov + = (inner == "bicgstab") || (inner == "gmres") || (inner == "cg") || (inner == "flexgmres"); + const auto outer = prm.get("solver", "bicgstab"); + const bool outer_is_flexible = (outer == "flexgmres"); + if (inner_is_krylov && !outer_is_flexible) { + OPM_THROW(std::invalid_argument, + fmt::format("system_cpr is configured with an approximate (Krylov) well " + "solver, 'preconditioner.well_solver.solver' = '{}', which makes " + "the preconditioner vary between applications. The outer solver " + "must then be flexible: set 'solver' to 'flexgmres' (it is " + "currently '{}'), or use a stationary well solver such as " + "'umfpack' or 'preconditioner2inverseoperator'.", + inner, outer)); + } } } diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 316440ee306..5ed463df168 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -27,6 +27,15 @@ #include #include +#include +#include + +#include +#include +#include +#include +#include + namespace Opm { @@ -122,6 +131,13 @@ class ISTLSolverSystem : public ISTLSolver bool sysInitialized_ = false; WellMatrixStructure cachedWellStructure_; + // Aggregation of merged well block rows into wells, and the well weights + // used by the CPRW pressure stage. Both are produced here, in the outer + // layer, from data already extracted from the well model; the + // preconditioner consumes them as plain numbers. + WellDofLayout wellLayout_; + std::string wellWeightType_ = "quasiimpes"; + // Current per-well B/C/D blocks for the explicit 2x2 system matrix. std::vector> wellBMatrices_; std::vector> wellCMatrices_; @@ -171,6 +187,8 @@ class ISTLSolverSystem : public ISTLSolver this->simulator_.problem().wellModel().addBCDMatrix( wellBMatrices_, wellCMatrices_, wellDMatrices_, wellCells_); + buildWellDofLayout(); + const Opm::WellMatrixMerger merger( Parent::matrix_->N(), wellBMatrices_, wellCMatrices_, wellDMatrices_, wellCells_); @@ -189,6 +207,7 @@ class ISTLSolverSystem : public ISTLSolver const bool needStructureRefresh = !sysInitialized_ || globalStructureChanged; const auto& prm = this->prm_[this->activeSolverNum_]; + wellWeightType_ = prm.get("preconditioner.well_weight_type", std::string{"quasiimpes"}); if (needStructureRefresh) { OPM_TIMEBLOCK(flexibleSolverCreate); @@ -197,6 +216,7 @@ class ISTLSolverSystem : public ISTLSolver sysMatrix_.B = &mergedB_; sysMatrix_.C = &mergedC_; sysMatrix_.D = &mergedD_; + sysMatrix_.wellLayout = &wellLayout_; cachedWellStructure_ = merger.buildStructure(); refreshSystemSolverForChangedWellStructure(prm); @@ -212,10 +232,82 @@ class ISTLSolverSystem : public ISTLSolver sysMatrix_.B = &mergedB_; sysMatrix_.C = &mergedC_; sysMatrix_.D = &mergedD_; + sysMatrix_.wellLayout = &wellLayout_; sysPrecond_->update(); } } + // Which merged well block rows belong to which well. The merged D matrix + // is the per-well D blocks concatenated, so this is a plain prefix sum + // over their dimensions: one block for a standard well, one per segment + // for a multisegment well. + void buildWellDofLayout() + { + auto& offsets = wellLayout_.wellBlockOffsets; + offsets.clear(); + offsets.reserve(wellDMatrices_.size() + 1); + offsets.push_back(0); + std::size_t total = 0; + for (const auto& d : wellDMatrices_) { + total += d.N(); + offsets.push_back(total); + } + } + + // Weights used to contract each well's equations down to the single scalar + // the CPRW pressure system carries for that well. Computed here rather + // than inside the preconditioner so that the linear-solver core never sees + // anything well-specific, and so that this can later be replaced by a + // value obtained from the well model without touching the core. + WellVector computeWellWeights() const + { + const std::size_t numBlocks = mergedD_.N(); + const int q = wellLayout_.pressureDofIndex; + + WellVector weights(numBlocks); + for (std::size_t wb = 0; wb < numBlocks; ++wb) { + auto& lambda = weights[wb]; + lambda = 0.0; + + if (wellWeightType_ == "unit") { + // Pick the pressure row of the well equations as-is. + lambda[q] = 1.0; + continue; + } + + // Quasi-IMPES well weights: lambda = D_ii^-T e_q, scaled to unit + // max norm. This is the analogue of the use_well_weights=true + // branch of StandardWellEquations::extractCPRPressureMatrix, and + // it needs no knowledge of the well's control mode. + Dune::FieldVector rhs(0.0); + rhs[q] = 1.0; + bool ok = false; + if (mergedD_.exists(wb, wb)) { + try { + const auto dt = mergedD_[wb][wb].transposed(); + dt.solve(lambda, rhs); + Scalar absMax = 0.0; + for (int i = 0; i < numWellDofs; ++i) { + absMax = std::max(absMax, std::abs(lambda[i])); + } + if (absMax > 0.0 && std::isfinite(absMax)) { + lambda /= absMax; + ok = true; + } + } catch (const Dune::FMatrixError&) { + ok = false; + } + } + if (!ok) { + // Singular or degenerate well block: fall back to the plain + // pressure row rather than poisoning the coarse system. + lambda = 0.0; + lambda[q] = 1.0; + } + } + return weights; + } + void refreshSystemSolverForChangedWellStructure(const Opm::PropertyTree& prm) { if (!sysInitialized_ || !sysPrecond_) { @@ -250,11 +342,15 @@ class ISTLSolverSystem : public ISTLSolver std::function()> resWeightCalc = this->getWeightsCalculator(resSolverPrm, this->getMatrix(), pressureIndex); + // The well part of the weights is filled here too: the CPRW pressure + // stage restricts the well rows with it, and re-reads it on every + // update, so it has to track the current merged D. std::function()> sysWeightCalc; if (resWeightCalc) { - sysWeightCalc = [resWeightCalc]() { + sysWeightCalc = [this, resWeightCalc]() { SystemVector w; w[_0] = resWeightCalc(); + w[_1] = this->computeWellWeights(); return w; }; } diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp new file mode 100644 index 00000000000..5c5e453d779 --- /dev/null +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -0,0 +1,470 @@ +/* + Copyright Equinor ASA 2026 + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_SYSTEMCPRWPRESSURESTAGE_HEADER_INCLUDED +#define OPM_SYSTEMCPRWPRESSURESTAGE_HEADER_INCLUDED + +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + +// -------------------------------------------------------------------------- +// CPRW pressure stage for the coupled (reservoir, well) system. +// +// Builds and solves the scalar pressure system of dimension +// +// Nres + nWells +// +// obtained by contracting the full system matrix +// +// S = [ A C ] +// [ B D ] +// +// with a restriction R and a prolongation P: +// +// R = blockdiag( w0^T ; sum_{wb in well j} w1[wb]^T ) +// P = blockdiag( e_p ; e_q placed on the top block row of well j ) +// +// where w0/w1 are the reservoir/well weights supplied from the outer layer, +// p is the reservoir pressure variable and q is the well pressure variable +// (bhp, or top segment pressure for a multisegment well). Every well +// contributes exactly one coarse unknown, as in the classic CPRW. +// +// The point of doing it this way is that the assembly reads nothing but the +// four sparse blocks, the weights and the WellDofLayout. No part of the well +// model is visible here. +// +// How the well part of the fine vectors takes part in the transfer. The +// coarse matrix is the same in every case; only the vector transfers differ. +// +// The classic PressureBhpTransferPolicy has no well unknowns on the fine level +// at all (the wells are Schur-eliminated in the operator), so it can neither +// restrict a well residual nor apply a coarse bhp correction. Selecting +// Classic here reproduces that, which makes the system solver and the classic +// cprw differ only in numerics rather than in formulation. +enum class WellTransfer +{ + Full, // restrict the well residual and prolong the bhp correction + NoProlongation, // restrict the well residual, discard the bhp correction + Classic, // neither -- as in PressureBhpTransferPolicy +}; + +inline WellTransfer wellTransferFromString(const std::string& name) +{ + if (name == "full") { + return WellTransfer::Full; + } + if (name == "no_prolongation") { + return WellTransfer::NoProlongation; + } + if (name == "classic") { + return WellTransfer::Classic; + } + OPM_THROW(std::invalid_argument, + "Unknown well_transfer '" + name + + "'. Valid values are 'full', 'no_prolongation' and 'classic'."); +} + +// -------------------------------------------------------------------------- +template +class SystemCprwPressureStage +{ +public: + static constexpr bool isParallel = !std::is_same_v; + + using CoarseOperator = Details::CoarseOperatorType; + using CoarseMatrix = typename CoarseOperator::matrix_type; + using CoarseVector = Details::PressureVectorType; + using CoarseSolver = Dune::FlexibleSolver; + + SystemCprwPressureStage(const SystemMatrix& S, + const PropertyTree& coarseSolverPrm, + const int pressureIndex, + const WellTransfer wellTransfer = WellTransfer::Full, + const Comm* comm = nullptr) + : S_(S) + , prm_(coarseSolverPrm) + , pressureIndex_(pressureIndex) + , wellTransfer_(wellTransfer) + , comm_(comm) + { + } + + // Whether a coarse correction reaches the well unknowns at all. The + // caller can skip the C and D defect updates when it does not. + bool prolongatesWellPressure() const + { + return wellTransfer_ == WellTransfer::Full; + } + + // (Re)create the coarse sparsity pattern, communication and entries. + // Separate from the coarse solver so that the assembly can be exercised on + // its own. + void buildCoarseSystem(const SystemVector& weights) + { + OPM_TIMEBLOCK(systemCprwBuildCoarseSystem); + const auto& layout = wellLayout(); + const std::size_t numRes = S_.A->N(); + const std::size_t numWells = layout.numWells(); + + buildCoarsePattern(numRes, numWells); + buildCoarseCommunication(numWells); + assembleCoarseMatrix(weights); + + coarseRhs_.resize(coarseMatrix_->N()); + coarseSol_.resize(coarseMatrix_->M()); + } + + // (Re)create the coarse system and the solver acting on it. Must be + // called whenever the well structure changes. + void buildStructure(const SystemVector& weights) + { + OPM_TIMEBLOCK(systemCprwBuildStructure); + buildCoarseSystem(weights); + + using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; + OperatorArgs oargs(coarseMatrix_, *coarseComm_); + coarseOperator_ = Dune::Amg::ConstructionTraits::construct(oargs); + + std::function noWeights; + if constexpr (isParallel) { + coarseSolver_ = std::make_unique(*coarseOperator_, *coarseComm_, + prm_, noWeights, /*pressureIndex=*/1); + } else { + coarseSolver_ = std::make_unique(*coarseOperator_, prm_, + noWeights, /*pressureIndex=*/1); + } + } + + // Recompute the coarse entries from the current system matrix and weights. + void update(const SystemVector& weights) + { + OPM_TIMEBLOCK(systemCprwUpdate); + assembleCoarseMatrix(weights); + coarseSolver_->preconditioner().update(); + } + + // One pressure-stage application: restrict, coarse solve, prolong. + void apply(const ResVector& dRes, + const WellVector& dWell, + const SystemVector& weights, + ResVector& vRes, + WellVector& vWell) + { + OPM_TIMEBLOCK(systemCprwApply); + moveToCoarseLevel(dRes, dWell, weights); + + coarseSol_ = 0.0; + Dune::InverseOperatorResult result; + coarseSolver_->apply(coarseSol_, coarseRhs_, result); + + moveToFineLevel(vRes, vWell); + } + + // Restriction: coarseRhs = R * (dRes, dWell). Unlike the classic CPRW + // transfer policy, the well residual really is carried to the coarse + // level rather than dropped. + void moveToCoarseLevel(const ResVector& dRes, + const WellVector& dWell, + const SystemVector& weights) + { + using namespace Dune::Indices; + const auto& layout = wellLayout(); + const auto& w0 = weights[_0]; + const auto& w1 = weights[_1]; + const std::size_t numRes = dRes.size(); + + coarseRhs_ = 0.0; + + for (std::size_t c = 0; c < numRes; ++c) { + const auto& bw = w0[c]; + Scalar el = 0.0; + for (std::size_t i = 0; i < bw.size(); ++i) { + el += dRes[c][i] * bw[i]; + } + coarseRhs_[c] = el; + } + + if (wellTransfer_ == WellTransfer::Classic) { + // The classic policy leaves the well rows of the coarse right-hand + // side at zero; keep them zero so that the two formulations agree. + return; + } + + for (std::size_t j = 0; j < layout.numWells(); ++j) { + Scalar el = 0.0; + for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { + const auto& lw = w1[wb]; + for (std::size_t i = 0; i < lw.size(); ++i) { + el += lw[i] * dWell[wb][i]; + } + } + coarseRhs_[numRes + j] = el; + } + } + + // Prolongation: (vRes, vWell) = P * coarseSol. The coarse well + // correction lands on the pressure unknown of each well's top block -- + // the classic CPRW throws it away instead. + void moveToFineLevel(ResVector& vRes, WellVector& vWell) const + { + const auto& layout = wellLayout(); + const std::size_t numRes = vRes.size(); + const int q = layout.pressureDofIndex; + + vRes = 0.0; + for (std::size_t c = 0; c < numRes; ++c) { + vRes[c][pressureIndex_] = coarseSol_[c][0]; + } + + vWell = 0.0; + if (!prolongatesWellPressure()) { + // The coarse bhp correction is computed but discarded; it acts only + // through its influence on the reservoir pressure, as in the + // classic policy. Stages 2 and 3 correct the well unknowns. + return; + } + for (std::size_t j = 0; j < layout.numWells(); ++j) { + vWell[layout.firstBlock(j)][q] = coarseSol_[numRes + j][0]; + } + } + + const CoarseMatrix& coarseMatrix() const + { + return *coarseMatrix_; + } + + const CoarseVector& coarseRhs() const + { + return coarseRhs_; + } + + CoarseVector& coarseSolution() + { + return coarseSol_; + } + +private: + const WellDofLayout& wellLayout() const + { + if (S_.wellLayout == nullptr) { + OPM_THROW(std::logic_error, + "SystemCprwPressureStage requires a WellDofLayout on the system matrix. " + "It is filled by ISTLSolverSystem; a null layout means the CPRW pressure " + "stage was constructed outside that path."); + } + return *S_.wellLayout; + } + + // Pattern: the reservoir block keeps A's pattern; each well j adds one row + // and one column, coupled to exactly the cells its B/C rows touch, plus a + // diagonal. The merged D is block diagonal by well, so the well-well part + // of the coarse system is diagonal. + void buildCoarsePattern(const std::size_t numRes, const std::size_t numWells) + { + const auto& A = *S_.A; + const auto& B = *S_.B; + const auto& C = *S_.C; + const auto& layout = wellLayout(); + + const std::size_t dim = numRes + numWells; + const std::size_t averageEntriesPerRow + = static_cast(std::ceil(static_cast(A.nonzeroes()) / A.N())); + const double overflowFraction = 1.2; + coarseMatrix_ = std::make_shared(dim, dim, + averageEntriesPerRow, + overflowFraction, + CoarseMatrix::implicit); + + // Reservoir-reservoir: A's pattern. + for (auto row = A.begin(), rowEnd = A.end(); row != rowEnd; ++row) { + for (auto col = row->begin(), colEnd = row->end(); col != colEnd; ++col) { + coarseMatrix_->entry(row.index(), col.index()) = 0.0; + } + } + + for (std::size_t j = 0; j < numWells; ++j) { + const std::size_t wdof = numRes + j; + coarseMatrix_->entry(wdof, wdof) = 0.0; + + // Well row: the cells reached by any block row of this well. + for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { + for (auto col = B[wb].begin(), colEnd = B[wb].end(); col != colEnd; ++col) { + coarseMatrix_->entry(wdof, col.index()) = 0.0; + } + } + } + + // Well column: the cells whose C row references this well's top block. + for (std::size_t c = 0; c < numRes; ++c) { + for (auto col = C[c].begin(), colEnd = C[c].end(); col != colEnd; ++col) { + const auto j = wellOfBlock(col.index()); + if (j.has_value()) { + coarseMatrix_->entry(c, numRes + *j) = 0.0; + } + } + } + + coarseMatrix_->compress(); + } + + void buildCoarseCommunication([[maybe_unused]] const std::size_t numWells) + { + if constexpr (isParallel) { + coarseComm_ = std::make_shared(comm_->communicator(), comm_->category(), false); + // Well DOFs are rank local and owned, appended after the reservoir + // DOFs -- the same convention the classic CPRW coarse system uses. + extendCommunicatorWithWells(*comm_, coarseComm_, static_cast(numWells)); + } else { + coarseComm_ = std::make_shared(); + } + } + + void assembleCoarseMatrix(const SystemVector& weights) + { + using namespace Dune::Indices; + OPM_TIMEBLOCK(systemCprwAssemble); + + const auto& A = *S_.A; + const auto& B = *S_.B; + const auto& C = *S_.C; + const auto& D = *S_.D; + const auto& layout = wellLayout(); + const auto& w0 = weights[_0]; + const auto& w1 = weights[_1]; + + const std::size_t numRes = A.N(); + const int p = pressureIndex_; + const int q = layout.pressureDofIndex; + + *coarseMatrix_ = 0.0; + + // Reservoir rows, reservoir columns: sum_i w0[c][i] * A[c][c'][i][p] + for (auto row = A.begin(), rowEnd = A.end(); row != rowEnd; ++row) { + const auto& bw = w0[row.index()]; + for (auto col = row->begin(), colEnd = row->end(); col != colEnd; ++col) { + Scalar el = 0.0; + for (std::size_t i = 0; i < bw.size(); ++i) { + el += (*col)[i][p] * bw[i]; + } + (*coarseMatrix_)[row.index()][col.index()] = el; + } + } + + // Reservoir rows, well columns: sum_i w0[c][i] * C[c][b(j)][i][q]. + // Only the top block of each well carries a coarse unknown, because + // that is the only place the prolongation writes. + for (std::size_t c = 0; c < numRes; ++c) { + const auto& bw = w0[c]; + for (auto col = C[c].begin(), colEnd = C[c].end(); col != colEnd; ++col) { + const auto j = wellOfBlock(col.index()); + if (!j.has_value() || layout.firstBlock(*j) != col.index()) { + continue; + } + Scalar el = 0.0; + for (std::size_t i = 0; i < bw.size(); ++i) { + el += (*col)[i][q] * bw[i]; + } + (*coarseMatrix_)[c][numRes + *j] += el; + } + } + + // Well rows. + for (std::size_t j = 0; j < layout.numWells(); ++j) { + const std::size_t wdof = numRes + j; + for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { + const auto& lw = w1[wb]; + + // Well row, reservoir columns: + // sum_{wb in j} sum_i w1[wb][i] * B[wb][c][i][p] + for (auto col = B[wb].begin(), colEnd = B[wb].end(); col != colEnd; ++col) { + Scalar el = 0.0; + for (std::size_t i = 0; i < lw.size(); ++i) { + el += lw[i] * (*col)[i][p]; + } + (*coarseMatrix_)[wdof][col.index()] += el; + } + + // Well row, well columns: + // sum_{wb in j} sum_i w1[wb][i] * D[wb][b(k)][i][q] + for (auto col = D[wb].begin(), colEnd = D[wb].end(); col != colEnd; ++col) { + const auto k = wellOfBlock(col.index()); + if (!k.has_value() || layout.firstBlock(*k) != col.index()) { + continue; + } + Scalar el = 0.0; + for (std::size_t i = 0; i < lw.size(); ++i) { + el += lw[i] * (*col)[i][q]; + } + (*coarseMatrix_)[wdof][numRes + *k] += el; + } + } + } + } + + // Merged well block row -> well index. Linear scan is fine: the offsets + // are sorted and this is only used while walking sparse rows of the well + // blocks, which are short. + std::optional wellOfBlock(const std::size_t blockRow) const + { + const auto& offsets = wellLayout().wellBlockOffsets; + const auto it = std::upper_bound(offsets.begin(), offsets.end(), blockRow); + if (it == offsets.begin() || it == offsets.end()) { + return std::nullopt; + } + return static_cast(std::distance(offsets.begin(), it) - 1); + } + + const SystemMatrix& S_; + PropertyTree prm_; + int pressureIndex_ = 0; + WellTransfer wellTransfer_ = WellTransfer::Full; + const Comm* comm_ = nullptr; + + std::shared_ptr coarseComm_; + std::shared_ptr coarseMatrix_; + std::shared_ptr coarseOperator_; + std::unique_ptr coarseSolver_; + + CoarseVector coarseRhs_; + CoarseVector coarseSol_; +}; + +} // namespace Opm + +#endif // OPM_SYSTEMCPRWPRESSURESTAGE_HEADER_INCLUDED diff --git a/opm/simulators/linalg/system/SystemPreconditioner.hpp b/opm/simulators/linalg/system/SystemPreconditioner.hpp index 8a8d66acf92..79e5651062a 100644 --- a/opm/simulators/linalg/system/SystemPreconditioner.hpp +++ b/opm/simulators/linalg/system/SystemPreconditioner.hpp @@ -20,14 +20,21 @@ #define OPM_SYSTEMPRECONDITIONER_HEADER_INCLUDED #include +#include #include #include #include #include +#include + #include #include +#include +#include +#include + namespace Opm { // Reservoir operator/comm types used as template arguments. @@ -65,7 +72,7 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdate& S, - const std::function()>& weightsCalculator, + const std::function()>& weightsCalculator, int pressureIndex, const Opm::PropertyTree& prm) requires (!isParallel) @@ -78,7 +85,7 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdate& S, - const std::function()>& weightsCalculator, + const std::function()>& weightsCalculator, int pressureIndex, const Opm::PropertyTree& prm, const ResComm& resComm) @@ -109,14 +116,26 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdatepreconditioner().update(); + if (cprwStage_) { + weights_ = weightsCalculator_(); + cprwStage_->update(weights_); + } else { + resSolver_->preconditioner().update(); + } resSmoother_->preconditioner().update(); wellSolver_->preconditioner().update(); } void updateForChangedWellStructure() { - resSolver_->preconditioner().update(); + if (cprwStage_) { + // The coarse system carries one unknown per well, so a changed + // well structure changes its dimension and pattern. + weights_ = weightsCalculator_(); + cprwStage_->buildStructure(weights_); + } else { + resSolver_->preconditioner().update(); + } resSmoother_->preconditioner().update(); initWellSolver(); resizeWellWorkVectors(); @@ -149,8 +168,27 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdateapply(tmp_resRes_, wRes_, weights_, dresSol_, dwSol_); + resSol_ += dresSol_; + // resRes_ -= A * dresSol_ + A.mmv(dresSol_, resRes_); + // wRes_ -= B * dresSol_ + B.mmv(dresSol_, wRes_); + if (cprwStage_->prolongatesWellPressure()) { + wSol_ += dwSol_; + // resRes_ -= C * dwSol_ ; wRes_ -= D * dwSol_ + C.mmv(dwSol_, resRes_); + D.mmv(dwSol_, wRes_); + } + } else { Dune::InverseOperatorResult res_result; dresSol_ = 0.0; tmp_resRes_ = resRes_; @@ -212,6 +250,13 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdate resSmoother_; std::unique_ptr wellSolver_; + // Non-null when the pressure stage includes the well unknowns (CPRW). + // Then resSolver_ is not built and stage 1 goes through cprwStage_. + using CprwStage = SystemCprwPressureStage; + std::unique_ptr cprwStage_; + std::function()> weightsCalculator_; + SystemVector weights_; + WellVector wSol_; ResVector resSol_; ResVector dresSol_; @@ -237,24 +282,61 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdate()>& weightsCalculator) + const std::function()>& weightsCalculator) { auto resprm = prm.get_child("reservoir_solver"); auto resprmsmoother = prm.get_child("reservoir_smoother"); wellprm_ = prm.get_child("well_solver"); + // add_wells is the same switch the classic CPR/CPRW pair uses: it + // promotes the pressure stage from reservoir-only CPR to CPRW over the + // full (reservoir, well) system. + const bool addWells = resprm.get("preconditioner.add_wells", false); + + // The weights arrive from the outer layer for the whole system; the + // reservoir-only sub-solvers want just their own part of them. + std::function()> resWeightCalc; + if (weightsCalculator) { + resWeightCalc = [weightsCalculator]() { + return weightsCalculator()[_0]; + }; + } + if constexpr (isParallel) { rop_ = std::make_unique(S_[_0][_0], *resComm_); - resSolver_ = std::make_unique( - *rop_, *resComm_, resprm, weightsCalculator, pressureIndex_); resSmoother_ = std::make_unique( - *rop_, *resComm_, resprmsmoother, weightsCalculator, pressureIndex_); + *rop_, *resComm_, resprmsmoother, resWeightCalc, pressureIndex_); } else { rop_ = std::make_unique(S_[_0][_0]); - resSolver_ = std::make_unique( - *rop_, resprm, weightsCalculator, pressureIndex_); resSmoother_ = std::make_unique( - *rop_, resprmsmoother, weightsCalculator, pressureIndex_); + *rop_, resprmsmoother, resWeightCalc, pressureIndex_); + } + + if (addWells) { + if (!weightsCalculator) { + OPM_THROW(std::invalid_argument, + "The CPRW pressure stage (add_wells) needs a weights calculator, but " + "none was configured. Set reservoir_solver.preconditioner.weight_type."); + } + weightsCalculator_ = weightsCalculator; + weights_ = weightsCalculator_(); + + auto coarseprm = resprm.get_child_optional("preconditioner.coarsesolver") + ? resprm.get_child("preconditioner.coarsesolver") + : PropertyTree(); + const auto wellTransfer = wellTransferFromString( + prm.get("well_transfer", std::string{"full"})); + cprwStage_ = std::make_unique(S_, coarseprm, pressureIndex_, + wellTransfer, resComm_); + cprwStage_->buildStructure(weights_); + } else { + if constexpr (isParallel) { + resSolver_ = std::make_unique( + *rop_, *resComm_, resprm, resWeightCalc, pressureIndex_); + } else { + resSolver_ = std::make_unique( + *rop_, resprm, resWeightCalc, pressureIndex_); + } } initWellSolver(); diff --git a/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp b/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp index 7bd59137f99..07e8bcddaf7 100644 --- a/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp +++ b/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp @@ -38,14 +38,8 @@ void addSystemCprSeq() [](const O& op, const P& prm, const std::function& sysWeightCalc, std::size_t pressureIndex) { - std::function()> resWeightCalc; - if (sysWeightCalc) { - resWeightCalc = [sysWeightCalc]() { - return sysWeightCalc()[Dune::Indices::_0]; - }; - } return std::make_shared>>( - op.getmat(), resWeightCalc, pressureIndex, prm); + op.getmat(), sysWeightCalc, pressureIndex, prm); }); } @@ -67,14 +61,8 @@ void addSystemCprParSeq() [](const O& op, const P& prm, const std::function& sysWeightCalc, std::size_t pressureIndex) { - std::function()> resWeightCalc; - if (sysWeightCalc) { - resWeightCalc = [sysWeightCalc]() { - return sysWeightCalc()[Dune::Indices::_0]; - }; - } return std::make_shared>>( - op.getmat(), resWeightCalc, pressureIndex, prm); + op.getmat(), sysWeightCalc, pressureIndex, prm); }); } @@ -91,15 +79,9 @@ void addSystemCprPar() const std::function& sysWeightCalc, std::size_t pressureIndex, const Opm::SystemComm& comm) { - std::function()> resWeightCalc; - if (sysWeightCalc) { - resWeightCalc = [sysWeightCalc]() { - return sysWeightCalc()[Dune::Indices::_0]; - }; - } const auto& resComm = comm[Dune::Indices::_0]; return std::make_shared, Opm::ParResComm>>( - op.getmat(), resWeightCalc, pressureIndex, prm, resComm); + op.getmat(), sysWeightCalc, pressureIndex, prm, resComm); }); } #endif diff --git a/opm/simulators/linalg/system/SystemTypes.hpp b/opm/simulators/linalg/system/SystemTypes.hpp index 093f1cd150d..7feeb63fa48 100644 --- a/opm/simulators/linalg/system/SystemTypes.hpp +++ b/opm/simulators/linalg/system/SystemTypes.hpp @@ -26,6 +26,9 @@ #include #include +#include +#include + namespace Opm { @@ -58,6 +61,58 @@ using WellVector = Dune::BlockVector>; template using SystemVector = Dune::MultiTypeBlockVector, WellVector>; +// -------------------------------------------------------------------------- +// WellDofLayout: which block rows of the merged well matrices belong to which +// well, plus the position of the pressure-like unknown inside a well block. +// +// The merged D matrix is block diagonal by well (WellMatrixMerger simply +// concatenates the per-well blocks), with one block row per standard well and +// one per segment of a multisegment well. Everything the preconditioner needs +// in order to aggregate well DOFs back to wells is therefore a prefix sum over +// the per-well D dimensions, which the outer layer already has. +// +// This is deliberately plain data: it is filled by ISTLSolverSystem from the +// matrices it already extracted, so that nothing below that point has to know +// anything about the well model. +// -------------------------------------------------------------------------- +struct WellDofLayout +{ + // Size numWells()+1, prefix sum of the per-well D_j.N(). + std::vector wellBlockOffsets; + + // Index of the pressure-like unknown (bhp for a standard well, segment + // pressure for a multisegment well) inside a well block. numWellDofs-1 is + // correct for the only configuration ISTLSolverSystem supports + // (Indices::numEq == 3, no energy): StandardWellPrimaryVariables::Bhp is + // numStaticWellEq - numWellControlEq == 3 and + // MultisegmentWellPrimaryVariables::SPres is + // has_wfrac + has_gfrac + 1 + enable_energy == 3. Carried as data so that + // generalising later is a change in the outer layer only. + int pressureDofIndex = numWellDofs - 1; + + std::size_t numWells() const + { + return wellBlockOffsets.empty() ? 0 : wellBlockOffsets.size() - 1; + } + + // First (top) block row of well j. For a multisegment well this is the + // top segment, whose pressure plays the role of the bhp. + std::size_t firstBlock(const std::size_t j) const + { + return wellBlockOffsets[j]; + } + + std::size_t endBlock(const std::size_t j) const + { + return wellBlockOffsets[j + 1]; + } + + std::size_t totalWellBlocks() const + { + return wellBlockOffsets.empty() ? 0 : wellBlockOffsets.back(); + } +}; + // -------------------------------------------------------------------------- // SystemMatrix: a lightweight read-only view over a 2×2 block-matrix // structure. All four sub-blocks are stored as const pointers; the actual @@ -88,6 +143,10 @@ class SystemMatrix const WRMatrix* B = nullptr; // (1,0) well–reservoir coupling const WWMatrix* D = nullptr; // (1,1) well + // Aggregation of the well block rows into wells. Only needed by the CPRW + // pressure stage; null when the well DOFs are not aggregated. + const WellDofLayout* wellLayout = nullptr; + // Sub-block access: S[_0][_0], S[_0][_1], S[_1][_0], S[_1][_1] inline SystemMatrixRow0 operator[](Dune::index_constant<0>) const; inline SystemMatrixRow1 operator[](Dune::index_constant<1>) const; diff --git a/tests/options_system_cprw_approx_wells.json b/tests/options_system_cprw_approx_wells.json new file mode 100644 index 00000000000..2ab9229acba --- /dev/null +++ b/tests/options_system_cprw_approx_wells.json @@ -0,0 +1,77 @@ +{ + "maxiter": "20", + "tol": "0.005", + "verbosity": "0", + "solver": "flexgmres", + "preconditioner": { + "type": "system_cpr", + "well_weight_type": "quasiimpes", + "reservoir_smoother": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "paroverilu0", + "relaxation": "1" + } + }, + "reservoir_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "cpr", + "relaxation": "1", + "use_well_weights": "false", + "add_wells": "true", + "weight_type": "trueimpes", + "pre_smooth": "0", + "post_smooth": "0", + "verbosity": "0", + "finesmoother": { + "type": "jac", + "relaxation": "1" + }, + "coarsesolver": { + "maxiter": "1", + "tol": "0.1", + "solver": "loopsolver", + "verbosity": "0", + "preconditioner": { + "type": "amg", + "alpha": "0.333333333333", + "relaxation": "1", + "iterations": "1", + "coarsenTarget": "1200", + "pre_smooth": "1", + "post_smooth": "1", + "beta": "0", + "smoother": "ilu0", + "verbosity": "0", + "maxlevel": "15", + "skip_isolated": "0", + "accumulate": "1", + "prolongationdamping": "1", + "maxdistance": "2", + "maxconnectivity": "15", + "maxaggsize": "6", + "minaggsize": "4" + } + } + } + }, + "well_solver": { + "maxiter": "5", + "tol": "0.01", + "verbosity": "0", + "solver": "bicgstab", + "preconditioner": { + "type": "ilu0", + "relaxation": "1" + } + } + }, + "restart": "20" +} \ No newline at end of file diff --git a/tests/options_system_cprw_approx_wells_bad_outer.json b/tests/options_system_cprw_approx_wells_bad_outer.json new file mode 100644 index 00000000000..94693733f67 --- /dev/null +++ b/tests/options_system_cprw_approx_wells_bad_outer.json @@ -0,0 +1,76 @@ +{ + "maxiter": "20", + "tol": "0.005", + "verbosity": "0", + "solver": "bicgstab", + "preconditioner": { + "type": "system_cpr", + "well_weight_type": "quasiimpes", + "reservoir_smoother": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "paroverilu0", + "relaxation": "1" + } + }, + "reservoir_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "cpr", + "relaxation": "1", + "use_well_weights": "false", + "add_wells": "true", + "weight_type": "trueimpes", + "pre_smooth": "0", + "post_smooth": "0", + "verbosity": "0", + "finesmoother": { + "type": "jac", + "relaxation": "1" + }, + "coarsesolver": { + "maxiter": "1", + "tol": "0.1", + "solver": "loopsolver", + "verbosity": "0", + "preconditioner": { + "type": "amg", + "alpha": "0.333333333333", + "relaxation": "1", + "iterations": "1", + "coarsenTarget": "1200", + "pre_smooth": "1", + "post_smooth": "1", + "beta": "0", + "smoother": "ilu0", + "verbosity": "0", + "maxlevel": "15", + "skip_isolated": "0", + "accumulate": "1", + "prolongationdamping": "1", + "maxdistance": "2", + "maxconnectivity": "15", + "maxaggsize": "6", + "minaggsize": "4" + } + } + } + }, + "well_solver": { + "maxiter": "5", + "tol": "0.01", + "verbosity": "0", + "solver": "bicgstab", + "preconditioner": { + "type": "ilu0", + "relaxation": "1" + } + } + } +} \ No newline at end of file diff --git a/tests/options_system_cprw_complete.json b/tests/options_system_cprw_complete.json new file mode 100644 index 00000000000..cee3c6b1468 --- /dev/null +++ b/tests/options_system_cprw_complete.json @@ -0,0 +1,72 @@ +{ + "maxiter": "20", + "tol": "0.005", + "verbosity": "0", + "solver": "bicgstab", + "preconditioner": { + "type": "system_cpr", + "well_weight_type": "quasiimpes", + "reservoir_smoother": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "paroverilu0", + "relaxation": "1" + } + }, + "reservoir_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "cpr", + "relaxation": "1", + "use_well_weights": "false", + "add_wells": "true", + "weight_type": "trueimpes", + "pre_smooth": "0", + "post_smooth": "0", + "verbosity": "0", + "finesmoother": { + "type": "jac", + "relaxation": "1" + }, + "coarsesolver": { + "maxiter": "1", + "tol": "0.1", + "solver": "loopsolver", + "verbosity": "0", + "preconditioner": { + "type": "amg", + "alpha": "0.333333333333", + "relaxation": "1", + "iterations": "1", + "coarsenTarget": "1200", + "pre_smooth": "1", + "post_smooth": "1", + "beta": "0", + "smoother": "ilu0", + "verbosity": "0", + "maxlevel": "15", + "skip_isolated": "0", + "accumulate": "1", + "prolongationdamping": "1", + "maxdistance": "2", + "maxconnectivity": "15", + "maxaggsize": "6", + "minaggsize": "4" + } + } + } + }, + "well_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "umfpack" + } + } +} diff --git a/tests/options_system_cprw_missing_coarsesolver.json b/tests/options_system_cprw_missing_coarsesolver.json new file mode 100644 index 00000000000..1636d8fd732 --- /dev/null +++ b/tests/options_system_cprw_missing_coarsesolver.json @@ -0,0 +1,40 @@ +{ + "maxiter": "20", + "tol": "0.005", + "verbosity": "0", + "solver": "bicgstab", + "preconditioner": { + "type": "system_cpr", + "reservoir_smoother": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "paroverilu0", + "relaxation": "1" + } + }, + "reservoir_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "cpr", + "relaxation": "1", + "use_well_weights": "false", + "add_wells": "true", + "weight_type": "trueimpes", + "pre_smooth": "0", + "post_smooth": "0" + } + }, + "well_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "umfpack" + } + } +} diff --git a/tests/test_SystemCprwPressureStage.cpp b/tests/test_SystemCprwPressureStage.cpp new file mode 100644 index 00000000000..42d61093baa --- /dev/null +++ b/tests/test_SystemCprwPressureStage.cpp @@ -0,0 +1,486 @@ +/* + Copyright Equinor ASA 2026 + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include +#define BOOST_TEST_MODULE OPM_test_SystemCprwPressureStage +#include + +#include +#include + +#include +#include +#include + +namespace { + +using Scalar = double; +using RRMatrix = Opm::RRMatrix; +using RWMatrix = Opm::RWMatrix; +using WRMatrix = Opm::WRMatrix; +using WWMatrix = Opm::WWMatrix; +using SystemVector = Opm::SystemVector; +using Stage = Opm::SystemCprwPressureStage; + +constexpr int numRes = Opm::numResDofs; // 3 +constexpr int numWell = Opm::numWellDofs; // 4 +constexpr int pressureIndex = 0; +constexpr int wellPressureIndex = numWell - 1; + +// The test fixture below models 4 reservoir cells and 2 wells: +// well 0 - a standard well, 1 block row, perforating cells 0 and 1 +// well 1 - a multisegment well, 3 block rows (segments), perforating +// cells 2 and 3 from different segments +// so that the per-well aggregation over block rows is actually exercised. +constexpr std::size_t numCells = 4; +constexpr std::size_t numWellBlocks = 4; // 1 + 3 + +struct BlockSpec +{ + std::size_t column; + Scalar base; +}; + +using MatrixPattern = std::vector>; + +// Distinct, non-symmetric values so that a transposed or mis-indexed +// contraction cannot accidentally produce the right answer. +template +Block makeBlock(const Scalar base) +{ + Block block; + for (int row = 0; row < Block::rows; ++row) { + for (int col = 0; col < Block::cols; ++col) { + block[row][col] = base + 3.0 * row + 7.0 * col + 0.25 * row * col; + } + } + return block; +} + +template +Matrix buildMatrix(const std::size_t rows, const std::size_t cols, const MatrixPattern& pattern) +{ + std::size_t nonzeroes = 0; + for (const auto& row : pattern) { + nonzeroes += row.size(); + } + + Matrix matrix(rows, cols, nonzeroes, Matrix::row_wise); + for (auto row = matrix.createbegin(); row != matrix.createend(); ++row) { + for (const auto& e : pattern[row.index()]) { + row.insert(e.column); + } + } + for (std::size_t row = 0; row < rows; ++row) { + for (const auto& e : pattern[row]) { + matrix[row][e.column] = makeBlock(e.base); + } + } + return matrix; +} + +struct Fixture +{ + RRMatrix A; + RWMatrix C; + WRMatrix B; + WWMatrix D; + Opm::WellDofLayout layout; + Opm::SystemMatrix S; + SystemVector weights; + + Fixture() + { + // A: tridiagonal over the 4 cells. + A = buildMatrix(numCells, numCells, + {{{0, 1.0}, {1, 2.0}}, + {{0, 3.0}, {1, 4.0}, {2, 5.0}}, + {{1, 6.0}, {2, 7.0}, {3, 8.0}}, + {{2, 9.0}, {3, 10.0}}}); + + // C: cell -> well block. Cells 0,1 see well 0 (block 0); cells 2,3 + // see well 1 through segments 1 and 2 (blocks 2 and 3). Only the top + // block of each well carries a coarse unknown, so the entries on + // blocks 2 and 3 must be ignored by the assembly. + C = buildMatrix(numCells, numWellBlocks, + {{{0, 11.0}}, + {{0, 12.0}}, + {{1, 13.0}, {2, 14.0}}, + {{1, 15.0}, {3, 16.0}}}); + + // B: well block -> cell. + B = buildMatrix(numWellBlocks, numCells, + {{{0, 17.0}, {1, 18.0}}, + {{2, 19.0}}, + {{2, 20.0}}, + {{3, 21.0}}}); + + // D: block diagonal by well. Well 0 is a single block; well 1 is a + // 3x3 segment coupling with blocks 1..3. + D = buildMatrix(numWellBlocks, numWellBlocks, + {{{0, 22.0}}, + {{1, 23.0}, {2, 24.0}}, + {{1, 25.0}, {2, 26.0}, {3, 27.0}}, + {{2, 28.0}, {3, 29.0}}}); + + layout.wellBlockOffsets = {0, 1, 4}; + layout.pressureDofIndex = wellPressureIndex; + + S.A = &A; + S.C = &C; + S.B = &B; + S.D = &D; + S.wellLayout = &layout; + + weights[Dune::Indices::_0].resize(numCells); + for (std::size_t c = 0; c < numCells; ++c) { + for (int i = 0; i < numRes; ++i) { + weights[Dune::Indices::_0][c][i] = 0.5 + 0.1 * c + 0.3 * i; + } + } + weights[Dune::Indices::_1].resize(numWellBlocks); + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + for (int i = 0; i < numWell; ++i) { + weights[Dune::Indices::_1][wb][i] = 0.2 + 0.7 * wb - 0.15 * i; + } + } + } +}; + +// Independent brute-force reference: densify S, build R and P explicitly and +// form R*S*P. Deliberately written without reusing any production helper. +std::vector> referenceCoarseMatrix(const Fixture& f) +{ + const std::size_t nWells = f.layout.numWells(); + const std::size_t fineDim = numCells * numRes + numWellBlocks * numWell; + const std::size_t coarseDim = numCells + nWells; + + // Dense fine system. + std::vector> S(fineDim, std::vector(fineDim, 0.0)); + const auto resOff = [](const std::size_t c, const int i) { return c * numRes + i; }; + const auto wellOff = [](const std::size_t wb, const int i) { + return numCells * numRes + wb * numWell + i; + }; + + for (auto row = f.A.begin(); row != f.A.end(); ++row) { + for (auto col = row->begin(); col != row->end(); ++col) { + for (int i = 0; i < numRes; ++i) { + for (int j = 0; j < numRes; ++j) { + S[resOff(row.index(), i)][resOff(col.index(), j)] = (*col)[i][j]; + } + } + } + } + for (auto row = f.C.begin(); row != f.C.end(); ++row) { + for (auto col = row->begin(); col != row->end(); ++col) { + for (int i = 0; i < numRes; ++i) { + for (int j = 0; j < numWell; ++j) { + S[resOff(row.index(), i)][wellOff(col.index(), j)] = (*col)[i][j]; + } + } + } + } + for (auto row = f.B.begin(); row != f.B.end(); ++row) { + for (auto col = row->begin(); col != row->end(); ++col) { + for (int i = 0; i < numWell; ++i) { + for (int j = 0; j < numRes; ++j) { + S[wellOff(row.index(), i)][resOff(col.index(), j)] = (*col)[i][j]; + } + } + } + } + for (auto row = f.D.begin(); row != f.D.end(); ++row) { + for (auto col = row->begin(); col != row->end(); ++col) { + for (int i = 0; i < numWell; ++i) { + for (int j = 0; j < numWell; ++j) { + S[wellOff(row.index(), i)][wellOff(col.index(), j)] = (*col)[i][j]; + } + } + } + } + + // R (coarseDim x fineDim) and P (fineDim x coarseDim). + std::vector> R(coarseDim, std::vector(fineDim, 0.0)); + std::vector> P(fineDim, std::vector(coarseDim, 0.0)); + + for (std::size_t c = 0; c < numCells; ++c) { + for (int i = 0; i < numRes; ++i) { + R[c][resOff(c, i)] = f.weights[Dune::Indices::_0][c][i]; + } + P[resOff(c, pressureIndex)][c] = 1.0; + } + for (std::size_t j = 0; j < nWells; ++j) { + for (std::size_t wb = f.layout.firstBlock(j); wb < f.layout.endBlock(j); ++wb) { + for (int i = 0; i < numWell; ++i) { + R[numCells + j][wellOff(wb, i)] = f.weights[Dune::Indices::_1][wb][i]; + } + } + P[wellOff(f.layout.firstBlock(j), wellPressureIndex)][numCells + j] = 1.0; + } + + std::vector> coarse(coarseDim, std::vector(coarseDim, 0.0)); + for (std::size_t r = 0; r < coarseDim; ++r) { + for (std::size_t c = 0; c < coarseDim; ++c) { + Scalar sum = 0.0; + for (std::size_t k = 0; k < fineDim; ++k) { + if (R[r][k] == 0.0) { + continue; + } + for (std::size_t l = 0; l < fineDim; ++l) { + if (P[l][c] != 0.0) { + sum += R[r][k] * S[k][l] * P[l][c]; + } + } + } + coarse[r][c] = sum; + } + } + return coarse; +} + +Scalar coarseEntry(const Stage& stage, const std::size_t row, const std::size_t col) +{ + const auto& m = stage.coarseMatrix(); + if (!m.exists(row, col)) { + return 0.0; + } + return m[row][col][0][0]; +} + +} // anonymous namespace + +// The coarse system must be exactly R*S*P. This fails if the C or B +// contraction is dropped, if the wrong block row is taken as a well's coarse +// unknown, or if the reservoir/well pressure index is wrong. +BOOST_AUTO_TEST_CASE(CoarseMatrixEqualsRestrictedSystem) +{ + const Fixture f; + Stage stage(f.S, Opm::PropertyTree(), pressureIndex); + stage.buildCoarseSystem(f.weights); + + const auto expected = referenceCoarseMatrix(f); + const std::size_t coarseDim = numCells + f.layout.numWells(); + + BOOST_REQUIRE_EQUAL(stage.coarseMatrix().N(), coarseDim); + BOOST_REQUIRE_EQUAL(stage.coarseMatrix().M(), coarseDim); + + for (std::size_t r = 0; r < coarseDim; ++r) { + for (std::size_t c = 0; c < coarseDim; ++c) { + BOOST_CHECK_CLOSE(coarseEntry(stage, r, c), expected[r][c], 1e-10); + } + } +} + +// The well coupling must actually be present: every well row/column pair that +// the fixture perforates has to be non-zero. Guards against a coarse system +// that is merely the reservoir block padded with an identity. +BOOST_AUTO_TEST_CASE(WellCouplingIsPresent) +{ + const Fixture f; + Stage stage(f.S, Opm::PropertyTree(), pressureIndex); + stage.buildCoarseSystem(f.weights); + + // Well 0 perforates cells 0 and 1, well 1 cells 2 and 3. + const std::vector> perfs = {{0, 1}, {2, 3}}; + for (std::size_t j = 0; j < perfs.size(); ++j) { + const std::size_t wdof = numCells + j; + BOOST_CHECK_NE(coarseEntry(stage, wdof, wdof), 0.0); + for (const auto c : perfs[j]) { + BOOST_CHECK_NE(coarseEntry(stage, wdof, c), 0.0); + BOOST_CHECK_NE(coarseEntry(stage, c, wdof), 0.0); + } + } + + // Cells belonging to one well must not couple to the other well. + BOOST_CHECK_EQUAL(coarseEntry(stage, numCells + 0, 2), 0.0); + BOOST_CHECK_EQUAL(coarseEntry(stage, numCells + 1, 0), 0.0); +} + +// The reservoir block of the coarse system must be the plain CPR contraction, +// i.e. adding the wells must not perturb the reservoir rows. +BOOST_AUTO_TEST_CASE(ReservoirBlockIsPlainCprContraction) +{ + const Fixture f; + Stage stage(f.S, Opm::PropertyTree(), pressureIndex); + stage.buildCoarseSystem(f.weights); + + for (auto row = f.A.begin(); row != f.A.end(); ++row) { + const auto& bw = f.weights[Dune::Indices::_0][row.index()]; + for (auto col = row->begin(); col != row->end(); ++col) { + Scalar expected = 0.0; + for (int i = 0; i < numRes; ++i) { + expected += (*col)[i][pressureIndex] * bw[i]; + } + BOOST_CHECK_CLOSE(coarseEntry(stage, row.index(), col.index()), expected, 1e-10); + } + } +} + +// Restriction and prolongation must be transposes of each other in the sense +// that R*P is the identity when the weights select exactly the unknowns the +// prolongation writes. +BOOST_AUTO_TEST_CASE(RestrictOfProlongIsIdentityForSelectingWeights) +{ + Fixture f; + // Weights that pick out precisely the prolonged components. + for (std::size_t c = 0; c < numCells; ++c) { + f.weights[Dune::Indices::_0][c] = 0.0; + f.weights[Dune::Indices::_0][c][pressureIndex] = 1.0; + } + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + f.weights[Dune::Indices::_1][wb] = 0.0; + } + for (std::size_t j = 0; j < f.layout.numWells(); ++j) { + f.weights[Dune::Indices::_1][f.layout.firstBlock(j)][wellPressureIndex] = 1.0; + } + + Stage stage(f.S, Opm::PropertyTree(), pressureIndex); + stage.buildCoarseSystem(f.weights); + + const std::size_t coarseDim = numCells + f.layout.numWells(); + auto& coarseSol = stage.coarseSolution(); + coarseSol.resize(coarseDim); + for (std::size_t i = 0; i < coarseDim; ++i) { + coarseSol[i] = 1.0 + 2.0 * i; + } + + Opm::ResVector vRes(numCells); + Opm::WellVector vWell(numWellBlocks); + stage.moveToFineLevel(vRes, vWell); + stage.moveToCoarseLevel(vRes, vWell, f.weights); + + for (std::size_t i = 0; i < coarseDim; ++i) { + BOOST_CHECK_CLOSE(stage.coarseRhs()[i][0], 1.0 + 2.0 * i, 1e-10); + } +} + +// well_transfer only changes the vector transfers -- the coarse matrix is the +// same in every mode. +BOOST_AUTO_TEST_CASE(WellTransferModeDoesNotChangeCoarseMatrix) +{ + const Fixture f; + const std::size_t coarseDim = numCells + f.layout.numWells(); + + Stage full(f.S, Opm::PropertyTree(), pressureIndex, Opm::WellTransfer::Full); + Stage classic(f.S, Opm::PropertyTree(), pressureIndex, Opm::WellTransfer::Classic); + full.buildCoarseSystem(f.weights); + classic.buildCoarseSystem(f.weights); + + for (std::size_t r = 0; r < coarseDim; ++r) { + for (std::size_t c = 0; c < coarseDim; ++c) { + BOOST_CHECK_EQUAL(coarseEntry(full, r, c), coarseEntry(classic, r, c)); + } + } +} + +// Classic mode must leave the well rows of the coarse right-hand side at zero +// and must not write any well correction, matching PressureBhpTransferPolicy. +BOOST_AUTO_TEST_CASE(ClassicTransferDropsWellResidualAndCorrection) +{ + const Fixture f; + const std::size_t coarseDim = numCells + f.layout.numWells(); + + Opm::ResVector dRes(numCells); + for (std::size_t c = 0; c < numCells; ++c) { + for (int i = 0; i < numRes; ++i) { + dRes[c][i] = 1.0 + c + i; + } + } + Opm::WellVector dWell(numWellBlocks); + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + for (int i = 0; i < numWell; ++i) { + dWell[wb][i] = 2.0 + wb - i; + } + } + + for (const auto mode : {Opm::WellTransfer::Full, + Opm::WellTransfer::NoProlongation, + Opm::WellTransfer::Classic}) { + Stage stage(f.S, Opm::PropertyTree(), pressureIndex, mode); + stage.buildCoarseSystem(f.weights); + stage.moveToCoarseLevel(dRes, dWell, f.weights); + + const bool restricts = (mode != Opm::WellTransfer::Classic); + for (std::size_t j = 0; j < f.layout.numWells(); ++j) { + if (restricts) { + BOOST_CHECK_NE(stage.coarseRhs()[numCells + j][0], 0.0); + } else { + BOOST_CHECK_EQUAL(stage.coarseRhs()[numCells + j][0], 0.0); + } + } + // The reservoir rows are restricted identically in every mode. + for (std::size_t c = 0; c < numCells; ++c) { + Scalar expected = 0.0; + for (int i = 0; i < numRes; ++i) { + expected += dRes[c][i] * f.weights[Dune::Indices::_0][c][i]; + } + BOOST_CHECK_CLOSE(stage.coarseRhs()[c][0], expected, 1e-10); + } + + auto& coarseSol = stage.coarseSolution(); + coarseSol.resize(coarseDim); + for (std::size_t i = 0; i < coarseDim; ++i) { + coarseSol[i] = 1.0 + i; + } + Opm::ResVector vRes(numCells); + Opm::WellVector vWell(numWellBlocks); + stage.moveToFineLevel(vRes, vWell); + + const bool prolongs = (mode == Opm::WellTransfer::Full); + BOOST_CHECK_EQUAL(stage.prolongatesWellPressure(), prolongs); + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + for (int i = 0; i < numWell; ++i) { + if (!prolongs) { + BOOST_CHECK_EQUAL(vWell[wb][i], 0.0); + } + } + } + if (prolongs) { + for (std::size_t j = 0; j < f.layout.numWells(); ++j) { + BOOST_CHECK_NE(vWell[f.layout.firstBlock(j)][wellPressureIndex], 0.0); + } + } + // The reservoir prolongation is the same in every mode. + for (std::size_t c = 0; c < numCells; ++c) { + BOOST_CHECK_CLOSE(vRes[c][pressureIndex], 1.0 + c, 1e-10); + } + } +} + +BOOST_AUTO_TEST_CASE(WellTransferFromStringRejectsUnknownValues) +{ + BOOST_CHECK(Opm::wellTransferFromString("full") == Opm::WellTransfer::Full); + BOOST_CHECK(Opm::wellTransferFromString("no_prolongation") == Opm::WellTransfer::NoProlongation); + BOOST_CHECK(Opm::wellTransferFromString("classic") == Opm::WellTransfer::Classic); + BOOST_CHECK_THROW(Opm::wellTransferFromString("nonsense"), std::invalid_argument); +} + +// A layout with no wells must reproduce a reservoir-only coarse system. +BOOST_AUTO_TEST_CASE(NoWellsGivesReservoirOnlyCoarseSystem) +{ + Fixture f; + Opm::WellDofLayout emptyLayout; + emptyLayout.wellBlockOffsets = {0}; + f.S.wellLayout = &emptyLayout; + + Stage stage(f.S, Opm::PropertyTree(), pressureIndex); + stage.buildCoarseSystem(f.weights); + + BOOST_CHECK_EQUAL(stage.coarseMatrix().N(), numCells); + BOOST_CHECK_EQUAL(stage.coarseMatrix().M(), numCells); +} diff --git a/tests/test_setuppropertytree.cpp b/tests/test_setuppropertytree.cpp index df040972a04..ee4996d840c 100644 --- a/tests/test_setuppropertytree.cpp +++ b/tests/test_setuppropertytree.cpp @@ -27,10 +27,12 @@ #include +#include #include #include #include +#include BOOST_AUTO_TEST_SUITE(SystemCPR) @@ -87,4 +89,54 @@ BOOST_AUTO_TEST_CASE(MatrixAddWellContributionsIncompatible) BOOST_CHECK_NO_THROW(Opm::checkSystemCPRMatrixAddWell(false)); } +// With add_wells the pressure stage is assembled and solved by the system +// preconditioner itself, taking its solver settings from the coarsesolver +// sub-tree. Without that sub-tree there is nothing to solve the CPRW pressure +// system with, so it must be rejected at setup time rather than falling over +// inside SystemCprwPressureStage::buildStructure. +BOOST_AUTO_TEST_CASE(JSONAddWellsRequiresCoarseSolver) +{ + Opm::PropertyTree prm("options_system_cprw_missing_coarsesolver.json"); + BOOST_CHECK_THROW(Opm::validateSystemCPRTree(prm), std::invalid_argument); + + Opm::PropertyTree complete("options_system_cprw_complete.json"); + BOOST_CHECK_NO_THROW(Opm::validateSystemCPRTree(complete)); +} + +// An approximate (Krylov) well solver stops on a tolerance and therefore does +// a different number of inner iterations per right-hand side, so the system +// preconditioner is no longer a fixed operator. Only a flexible outer solver +// may be combined with it; bicgstab must be rejected. +BOOST_AUTO_TEST_CASE(ApproximateWellSolverRequiresFlexibleOuterSolver) +{ + Opm::PropertyTree ok("options_system_cprw_approx_wells.json"); + BOOST_CHECK_NO_THROW(Opm::validateSystemCPRTree(ok)); + + Opm::PropertyTree bad("options_system_cprw_approx_wells_bad_outer.json"); + BOOST_CHECK_THROW(Opm::validateSystemCPRTree(bad), std::invalid_argument); + + // The stationary well solvers stay valid with the default outer solver. + Opm::PropertyTree exact("options_system_cprw_complete.json"); + BOOST_CHECK_EQUAL(exact.get("preconditioner.well_solver.solver"), "umfpack"); + BOOST_CHECK_NO_THROW(Opm::validateSystemCPRTree(exact)); +} + +// system_cprw must produce the same tree as system_cpr apart from add_wells. +BOOST_AUTO_TEST_CASE(SystemCPRWEnablesAddWells) +{ + const Opm::FlowLinearSolverParameters p; + const auto cpr = Opm::setupSystemCPR("system_cpr", p); + const auto cprw = Opm::setupSystemCPR("system_cprw", p); + + const std::string key = "preconditioner.reservoir_solver.preconditioner.add_wells"; + BOOST_CHECK_EQUAL(cpr.get(key), false); + BOOST_CHECK_EQUAL(cprw.get(key), true); + + // Same coarse solver in both, since that is what solves the pressure + // system in either case. + const std::string coarse + = "preconditioner.reservoir_solver.preconditioner.coarsesolver.preconditioner.type"; + BOOST_CHECK_EQUAL(cpr.get(coarse), cprw.get(coarse)); +} + BOOST_AUTO_TEST_SUITE_END() // SystemCPR From 4b5bd0e3f929b248af97133b7829d96d3512e4e1 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 14:32:27 +0200 Subject: [PATCH 02/22] Add cell-average well weights and regularise the coarse well diagonal Two pieces of #7209 that were missing here. The classic use_well_weights = false weighting is now available as well_weight_type = cellavg: average the reservoir weights over the cells a well block perforates and apply them to the conservation equations only, with weight zero on the control equation, falling back to unit weights when a well has no perforations on this rank. It is not a strictly worse or better choice than the quasi-IMPES default, so both are kept -- on SPE9_CP it needs 428 linear iterations against 470 for quasi-IMPES (and 439 for the classic cprw), while on SPE1 quasi-IMPES wins with 429 against 455. A well whose contraction cancels exactly leaves a zero on the coarse diagonal and makes the pressure system singular. Regularise that row to a unit diagonal instead of handing AMG a singular system. Co-Authored-By: Claude Opus 5 (cherry picked from commit bf21d5df04cd37b03126cc52b55d95d44fc79922) --- .../linalg/system/ISTLSolverSystem.hpp | 32 +++++++++++++++++-- .../linalg/system/SystemCprwPressureStage.hpp | 10 ++++++ tests/test_SystemCprwPressureStage.cpp | 26 +++++++++++++++ 3 files changed, 66 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 5ed463df168..7043ebebed4 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -259,7 +259,7 @@ class ISTLSolverSystem : public ISTLSolver // than inside the preconditioner so that the linear-solver core never sees // anything well-specific, and so that this can later be replaced by a // value obtained from the well model without touching the core. - WellVector computeWellWeights() const + WellVector computeWellWeights(const ResVector& resWeights) const { const std::size_t numBlocks = mergedD_.N(); const int q = wellLayout_.pressureDofIndex; @@ -275,6 +275,34 @@ class ISTLSolverSystem : public ISTLSolver continue; } + if (wellWeightType_ == "cellavg") { + // The classic CPRW default (use_well_weights = false): average + // the reservoir weights over the cells this block row + // perforates, and use them on the conservation equations only. + // The control equation gets weight zero. + int nperf = 0; + for (auto col = mergedB_[wb].begin(), end = mergedB_[wb].end(); col != end; ++col) { + const auto& cw = resWeights[col.index()]; + for (int i = 0; i < numResDofs; ++i) { + lambda[i] += cw[i]; + } + ++nperf; + } + if (nperf > 0) { + for (int i = 0; i < numResDofs; ++i) { + lambda[i] /= nperf; + } + } else { + // No perforations of this well on this rank; regularise + // rather than leaving an empty row. + for (int i = 0; i < numResDofs; ++i) { + lambda[i] = 1.0; + } + } + lambda[q] = 0.0; + continue; + } + // Quasi-IMPES well weights: lambda = D_ii^-T e_q, scaled to unit // max norm. This is the analogue of the use_well_weights=true // branch of StandardWellEquations::extractCPRPressureMatrix, and @@ -350,7 +378,7 @@ class ISTLSolverSystem : public ISTLSolver sysWeightCalc = [this, resWeightCalc]() { SystemVector w; w[_0] = resWeightCalc(); - w[_1] = this->computeWellWeights(); + w[_1] = this->computeWellWeights(w[_0]); return w; }; } diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp index 5c5e453d779..7121ba1fa5d 100644 --- a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -434,6 +434,16 @@ class SystemCprwPressureStage (*coarseMatrix_)[wdof][numRes + *k] += el; } } + + // A well whose contraction cancels exactly would leave a zero on + // the coarse diagonal and make the pressure system singular. That + // can happen for a well with no local perforations, or when the + // weights annihilate the pressure column. Regularise to a unit row + // rather than handing a singular system to AMG. + auto& diag = (*coarseMatrix_)[wdof][wdof][0][0]; + if (!(std::abs(diag) > 0.0)) { + diag = 1.0; + } } } diff --git a/tests/test_SystemCprwPressureStage.cpp b/tests/test_SystemCprwPressureStage.cpp index 42d61093baa..ab2fe358d33 100644 --- a/tests/test_SystemCprwPressureStage.cpp +++ b/tests/test_SystemCprwPressureStage.cpp @@ -470,6 +470,32 @@ BOOST_AUTO_TEST_CASE(WellTransferFromStringRejectsUnknownValues) BOOST_CHECK_THROW(Opm::wellTransferFromString("nonsense"), std::invalid_argument); } +// A well whose weights annihilate its pressure column would leave a zero on +// the coarse diagonal and hand AMG a singular system. It must be regularised +// to a unit diagonal instead. +BOOST_AUTO_TEST_CASE(ZeroCoarseWellDiagonalIsRegularised) +{ + Fixture f; + // Zero weights on every well block: the whole well row, including the + // diagonal, contracts to zero. + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + f.weights[Dune::Indices::_1][wb] = 0.0; + } + + Stage stage(f.S, Opm::PropertyTree(), pressureIndex); + stage.buildCoarseSystem(f.weights); + + for (std::size_t j = 0; j < f.layout.numWells(); ++j) { + const std::size_t wdof = numCells + j; + BOOST_CHECK_EQUAL(coarseEntry(stage, wdof, wdof), 1.0); + // The off-diagonal well row really is zero -- it is only the diagonal + // that gets regularised. + for (std::size_t c = 0; c < numCells; ++c) { + BOOST_CHECK_EQUAL(coarseEntry(stage, wdof, c), 0.0); + } + } +} + // A layout with no wells must reproduce a reservoir-only coarse system. BOOST_AUTO_TEST_CASE(NoWellsGivesReservoirOnlyCoarseSystem) { From e71c16a618a84b1953e8330d1632dc6c92d3d5e7 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 14:44:52 +0200 Subject: [PATCH 03/22] Default the CPRW pressure stage to the standard solver's weighting The pressure stage now defaults to well_weight_type = cellavg, so together with the trueimpes reservoir weights it already used, its weighting matches what cprw does by default (use_well_weights = false). Quasi-IMPES well weights remain available and are better on some cases, but the sensible default is the one the standard solver uses. A test pins both halves against setupCPRW so they cannot drift apart. Co-Authored-By: Claude Opus 5 (cherry picked from commit dce700140e17139f4bd0245819136d4de5adbf6d) --- opm/simulators/linalg/setupPropertyTree.cpp | 9 ++++++++- tests/options_system_cprw_approx_wells.json | 2 +- ...options_system_cprw_approx_wells_bad_outer.json | 2 +- tests/options_system_cprw_complete.json | 4 ++-- tests/test_setuppropertytree.cpp | 14 ++++++++++++++ 5 files changed, 26 insertions(+), 5 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index becfe8135ea..4f3943ded9c 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -532,7 +532,14 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) prm.put("preconditioner.type", "system_cpr"s); // How the well equations are contracted to the one coarse unknown each // well carries in the CPRW pressure system. Only read when add_wells. - prm.put("preconditioner.well_weight_type", "quasiimpes"s); + // cellavg - average of the reservoir weights over the well's + // perforated cells, on the conservation equations only. + // This is what cprw does (use_well_weights = false) and, + // together with the trueimpes reservoir weights below, + // makes the pressure stage match the standard solver. + // quasiimpes - inv(D)^T e_bhp, normalised. + // unit - the pressure row as-is; a debugging baseline. + prm.put("preconditioner.well_weight_type", "cellavg"s); // How the well unknowns take part in the pressure-stage transfer: // full - restrict the well residual, prolong the bhp correction // no_prolongation - restrict, but discard the bhp correction diff --git a/tests/options_system_cprw_approx_wells.json b/tests/options_system_cprw_approx_wells.json index 2ab9229acba..b885cc0c344 100644 --- a/tests/options_system_cprw_approx_wells.json +++ b/tests/options_system_cprw_approx_wells.json @@ -5,7 +5,7 @@ "solver": "flexgmres", "preconditioner": { "type": "system_cpr", - "well_weight_type": "quasiimpes", + "well_weight_type": "cellavg", "reservoir_smoother": { "maxiter": "1", "tol": "0.005", diff --git a/tests/options_system_cprw_approx_wells_bad_outer.json b/tests/options_system_cprw_approx_wells_bad_outer.json index 94693733f67..670707ee467 100644 --- a/tests/options_system_cprw_approx_wells_bad_outer.json +++ b/tests/options_system_cprw_approx_wells_bad_outer.json @@ -5,7 +5,7 @@ "solver": "bicgstab", "preconditioner": { "type": "system_cpr", - "well_weight_type": "quasiimpes", + "well_weight_type": "cellavg", "reservoir_smoother": { "maxiter": "1", "tol": "0.005", diff --git a/tests/options_system_cprw_complete.json b/tests/options_system_cprw_complete.json index cee3c6b1468..1cbad0c6c8c 100644 --- a/tests/options_system_cprw_complete.json +++ b/tests/options_system_cprw_complete.json @@ -5,7 +5,7 @@ "solver": "bicgstab", "preconditioner": { "type": "system_cpr", - "well_weight_type": "quasiimpes", + "well_weight_type": "cellavg", "reservoir_smoother": { "maxiter": "1", "tol": "0.005", @@ -69,4 +69,4 @@ "solver": "umfpack" } } -} +} \ No newline at end of file diff --git a/tests/test_setuppropertytree.cpp b/tests/test_setuppropertytree.cpp index ee4996d840c..eeef5ad1b44 100644 --- a/tests/test_setuppropertytree.cpp +++ b/tests/test_setuppropertytree.cpp @@ -137,6 +137,20 @@ BOOST_AUTO_TEST_CASE(SystemCPRWEnablesAddWells) const std::string coarse = "preconditioner.reservoir_solver.preconditioner.coarsesolver.preconditioner.type"; BOOST_CHECK_EQUAL(cpr.get(coarse), cprw.get(coarse)); + + // The pressure stage must default to the same weighting the standard cprw + // solver uses: trueimpes on the reservoir equations and the perforated-cell + // average on the well equations (cprw's use_well_weights = false). + BOOST_CHECK_EQUAL( + cprw.get("preconditioner.reservoir_solver.preconditioner.weight_type"), + "trueimpes"); + BOOST_CHECK_EQUAL(cprw.get("preconditioner.well_weight_type"), "cellavg"); + + const auto classic = Opm::setupCPRW("cprw", p); + BOOST_CHECK_EQUAL(classic.get("preconditioner.weight_type"), + cprw.get( + "preconditioner.reservoir_solver.preconditioner.weight_type")); + BOOST_CHECK_EQUAL(classic.get("preconditioner.use_well_weights"), false); } BOOST_AUTO_TEST_SUITE_END() // SystemCPR From 04ee7631fdb8f535e361919d254a1dbc7398fce3 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 22:41:53 +0200 Subject: [PATCH 04/22] Share one implementation for the pressure-stage transfers The restriction and prolongation grow caller-supplied-output forms, so that a caller which already owns the target vectors does not need a copy, and the two existing forms delegate to them instead of repeating the loops. The coarse matrix, communication and rhs also get handles, so the coarse level can be driven from outside the stage. No functional change: the transfers compute what they computed before, and the stage is still the only thing that assembles the coarse system. (cherry picked from commit 66d642a94132e12abff4c68dc1b51fb7e6461aa2, with the composable preconditioner parts and the Dune transfer policy left out -- they belong with the general system preconditioner, which is not part of this PR.) Co-Authored-By: Claude Opus 5 --- .../linalg/system/SystemCprwPressureStage.hpp | 163 ++++++++++++------ 1 file changed, 109 insertions(+), 54 deletions(-) diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp index 7121ba1fa5d..4618f75bc3f 100644 --- a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -200,42 +200,10 @@ class SystemCprwPressureStage // transfer policy, the well residual really is carried to the coarse // level rather than dropped. void moveToCoarseLevel(const ResVector& dRes, - const WellVector& dWell, - const SystemVector& weights) + const WellVector& dWell, + const SystemVector& weights) { - using namespace Dune::Indices; - const auto& layout = wellLayout(); - const auto& w0 = weights[_0]; - const auto& w1 = weights[_1]; - const std::size_t numRes = dRes.size(); - - coarseRhs_ = 0.0; - - for (std::size_t c = 0; c < numRes; ++c) { - const auto& bw = w0[c]; - Scalar el = 0.0; - for (std::size_t i = 0; i < bw.size(); ++i) { - el += dRes[c][i] * bw[i]; - } - coarseRhs_[c] = el; - } - - if (wellTransfer_ == WellTransfer::Classic) { - // The classic policy leaves the well rows of the coarse right-hand - // side at zero; keep them zero so that the two formulations agree. - return; - } - - for (std::size_t j = 0; j < layout.numWells(); ++j) { - Scalar el = 0.0; - for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { - const auto& lw = w1[wb]; - for (std::size_t i = 0; i < lw.size(); ++i) { - el += lw[i] * dWell[wb][i]; - } - } - coarseRhs_[numRes + j] = el; - } + restrictInto(dRes, dWell, weights, coarseRhs_); } // Prolongation: (vRes, vWell) = P * coarseSol. The coarse well @@ -243,25 +211,7 @@ class SystemCprwPressureStage // the classic CPRW throws it away instead. void moveToFineLevel(ResVector& vRes, WellVector& vWell) const { - const auto& layout = wellLayout(); - const std::size_t numRes = vRes.size(); - const int q = layout.pressureDofIndex; - - vRes = 0.0; - for (std::size_t c = 0; c < numRes; ++c) { - vRes[c][pressureIndex_] = coarseSol_[c][0]; - } - - vWell = 0.0; - if (!prolongatesWellPressure()) { - // The coarse bhp correction is computed but discarded; it acts only - // through its influence on the reservoir pressure, as in the - // classic policy. Stages 2 and 3 correct the well unknowns. - return; - } - for (std::size_t j = 0; j < layout.numWells(); ++j) { - vWell[layout.firstBlock(j)][q] = coarseSol_[numRes + j][0]; - } + prolongFrom(coarseSol_, vRes, vWell); } const CoarseMatrix& coarseMatrix() const @@ -269,6 +219,40 @@ class SystemCprwPressureStage return *coarseMatrix_; } + // Handles needed when the coarse level is driven from outside, e.g. by a + // Dune two-level transfer policy. + const std::shared_ptr& coarseMatrixPtr() const + { + return coarseMatrix_; + } + + const Comm& coarseCommunication() const + { + return *coarseComm_; + } + + void assembleCoarseEntries(const SystemVector& weights) + { + assembleCoarseMatrix(weights); + } + + // Transfer forms writing into a caller-supplied coarse vector, so that a + // transfer policy can use its own storage without an extra copy. + void moveToCoarseLevel(const ResVector& dRes, + const WellVector& dWell, + const SystemVector& weights, + CoarseVector& out) const + { + restrictInto(dRes, dWell, weights, out); + } + + void moveToFineLevel(const CoarseVector& in, + ResVector& vRes, + WellVector& vWell) const + { + prolongFrom(in, vRes, vWell); + } + const CoarseVector& coarseRhs() const { return coarseRhs_; @@ -447,6 +431,77 @@ class SystemCprwPressureStage } } + // Restriction: out = R * (dRes, dWell). Unlike the classic CPRW transfer + // policy the well residual really is carried to the coarse level, unless + // the classic transfer was asked for. + void restrictInto(const ResVector& dRes, + const WellVector& dWell, + const SystemVector& weights, + CoarseVector& out) const + { + using namespace Dune::Indices; + const auto& layout = wellLayout(); + const auto& w0 = weights[_0]; + const auto& w1 = weights[_1]; + const std::size_t numRes = dRes.size(); + + out = 0.0; + + for (std::size_t c = 0; c < numRes; ++c) { + const auto& bw = w0[c]; + Scalar el = 0.0; + for (std::size_t i = 0; i < bw.size(); ++i) { + el += dRes[c][i] * bw[i]; + } + out[c] = el; + } + + if (wellTransfer_ == WellTransfer::Classic) { + // The classic policy leaves the well rows of the coarse right-hand + // side at zero; keep them zero so that the two formulations agree. + return; + } + + for (std::size_t j = 0; j < layout.numWells(); ++j) { + Scalar el = 0.0; + for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { + const auto& lw = w1[wb]; + for (std::size_t i = 0; i < lw.size(); ++i) { + el += lw[i] * dWell[wb][i]; + } + } + out[numRes + j] = el; + } + } + + // Prolongation: (vRes, vWell) = P * in. The coarse well correction lands + // on the pressure unknown of each well's top block; the classic policy + // throws it away instead. + void prolongFrom(const CoarseVector& in, + ResVector& vRes, + WellVector& vWell) const + { + const auto& layout = wellLayout(); + const std::size_t numRes = vRes.size(); + const int q = layout.pressureDofIndex; + + vRes = 0.0; + for (std::size_t c = 0; c < numRes; ++c) { + vRes[c][pressureIndex_] = in[c][0]; + } + + vWell = 0.0; + if (!prolongatesWellPressure()) { + // The coarse bhp correction is computed but discarded; it acts only + // through its influence on the reservoir pressure, as in the + // classic policy. A following well solve corrects the wells. + return; + } + for (std::size_t j = 0; j < layout.numWells(); ++j) { + vWell[layout.firstBlock(j)][q] = in[numRes + j][0]; + } + } + // Merged well block row -> well index. Linear scan is fine: the offsets // are sorted and this is only used while walking sparse rows of the well // blocks, which are short. From 9e66699f7b4084f80df843f7c982da0224c5bf15 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 16:19:03 +0200 Subject: [PATCH 05/22] Sum the coarse well column over all of a well's blocks For a multisegment well the coarse column was taken from the well's top block alone, while MultisegmentWellEquations::extractCPRPressureMatrix accumulates over every segment row. With one segment per connection that discarded most of the well: on Norne with --convert-to-multisegment-well=per-connection the coupled pressure stage needed 4939 linear iterations against 2716 for the classic cprw, i.e. it lost its whole advantage, while standard wells were unaffected because one block per well makes the two conventions coincide. Sum the C and D column contractions over all of a well's block rows. That is the Galerkin column for a prolongation spreading a well's coarse unknown over all of its segment pressures by a constant, so the prolongation now spreads to match and the coarse operator stays R*S*P. The same case now takes 2494 linear iterations in 152.0 s against 2716 in 156.4 s for cprw. Default well_transfer becomes no_prolongation: the restriction is kept, the segment pressure correction is discarded and the trailing well solve corrects the wells. Prolonging one coarse value onto segment pressures as well as rates and compositions is not well defined, and it measures worse -- 2685 with the prolongation against 2494 without, while dropping the restriction too costs ~3 % (2565). Co-Authored-By: Claude Opus 5 (cherry picked from commit 496fb6d09fb9b7df7c3eb007cb5b5377891bd006) --- opm/simulators/linalg/setupPropertyTree.cpp | 2 +- .../linalg/system/SystemCprwPressureStage.hpp | 35 ++++++++++++++----- tests/test_SystemCprwPressureStage.cpp | 6 +++- 3 files changed, 33 insertions(+), 10 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 4f3943ded9c..46b743439ba 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -546,7 +546,7 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) // classic - neither, i.e. the classic cprw formulation, so that // the only remaining difference is numerics // Only read when add_wells. - prm.put("preconditioner.well_transfer", "full"s); + prm.put("preconditioner.well_transfer", "no_prolongation"s); // --- Reservoir smoother --- prm.put("preconditioner.reservoir_smoother.maxiter", 1); diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp index 4618f75bc3f..de4e5b4136f 100644 --- a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -314,7 +314,7 @@ class SystemCprwPressureStage } } - // Well column: the cells whose C row references this well's top block. + // Well column: every cell whose C row references any block of the well. for (std::size_t c = 0; c < numRes; ++c) { for (auto col = C[c].begin(), colEnd = C[c].end(); col != colEnd; ++col) { const auto j = wellOfBlock(col.index()); @@ -370,14 +370,22 @@ class SystemCprwPressureStage } } - // Reservoir rows, well columns: sum_i w0[c][i] * C[c][b(j)][i][q]. - // Only the top block of each well carries a coarse unknown, because - // that is the only place the prolongation writes. + // Reservoir rows, well columns: sum over every block row of well j of + // sum_i w0[c][i] * C[c][wb][i][q]. + // + // Summing over all of a well's blocks is the Galerkin column for a + // prolongation that spreads a well's coarse unknown over all of its + // segment pressures, which is what MultisegmentWellEquations:: + // extractCPRPressureMatrix does (it accumulates over every segment + // row). Taking the top block alone instead loses every segment but + // the first: on Norne with one segment per connection that is most of + // the well, and it is what made the coarse system far weaker than the + // classic cprw one for multisegment wells. for (std::size_t c = 0; c < numRes; ++c) { const auto& bw = w0[c]; for (auto col = C[c].begin(), colEnd = C[c].end(); col != colEnd; ++col) { const auto j = wellOfBlock(col.index()); - if (!j.has_value() || layout.firstBlock(*j) != col.index()) { + if (!j.has_value()) { continue; } Scalar el = 0.0; @@ -405,10 +413,15 @@ class SystemCprwPressureStage } // Well row, well columns: - // sum_{wb in j} sum_i w1[wb][i] * D[wb][b(k)][i][q] + // sum_{wb in j} sum_{wb' in k} sum_i w1[wb][i] * D[wb][wb'][i][q] + // Summed over all of well k's blocks, to match the column + // convention above. Because the merged D is block diagonal by + // well this only ever contributes to k == j, but it now picks + // up the full segment-to-segment coupling rather than just the + // top segment's column. for (auto col = D[wb].begin(), colEnd = D[wb].end(); col != colEnd; ++col) { const auto k = wellOfBlock(col.index()); - if (!k.has_value() || layout.firstBlock(*k) != col.index()) { + if (!k.has_value()) { continue; } Scalar el = 0.0; @@ -497,8 +510,14 @@ class SystemCprwPressureStage // classic policy. A following well solve corrects the wells. return; } + // Spread the well's coarse value over all of its segment pressures by + // a constant. This is the P the coarse matrix is assembled for; only + // the segment pressures are set, the other well unknowns (rates, + // compositions) are left to the well solve that follows. for (std::size_t j = 0; j < layout.numWells(); ++j) { - vWell[layout.firstBlock(j)][q] = in[numRes + j][0]; + for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { + vWell[wb][q] = in[numRes + j][0]; + } } } diff --git a/tests/test_SystemCprwPressureStage.cpp b/tests/test_SystemCprwPressureStage.cpp index ab2fe358d33..4930a48abb8 100644 --- a/tests/test_SystemCprwPressureStage.cpp +++ b/tests/test_SystemCprwPressureStage.cpp @@ -230,7 +230,11 @@ std::vector> referenceCoarseMatrix(const Fixture& f) R[numCells + j][wellOff(wb, i)] = f.weights[Dune::Indices::_1][wb][i]; } } - P[wellOff(f.layout.firstBlock(j), wellPressureIndex)][numCells + j] = 1.0; + // The prolongation spreads a well's coarse value over all of its + // segment pressures by a constant. + for (std::size_t wb = f.layout.firstBlock(j); wb < f.layout.endBlock(j); ++wb) { + P[wellOff(wb, wellPressureIndex)][numCells + j] = 1.0; + } } std::vector> coarse(coarseDim, std::vector(coarseDim, 0.0)); From add7a8871b5be9b3391d01961b4f86eaf622c9c0 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 16:50:03 +0200 Subject: [PATCH 06/22] Make cellavg the classic per-well weighting, add cellblockavg MultisegmentWellEquations::extractCPRPressureMatrix computes one weight vector per well, averaged over every perforation of every segment, and applies it to all of the well's rows. cellavg now does the same. The previous per-block-row averaging is kept as cellblockavg. On Norne with all wells converted to multisegment the two are close -- 2547 against 2565 linear iterations -- and both beat the classic cprw's 2716, so this is not why the system solver and cprw disagree: both still diverge from cprw at the sixth linear solve rather than in the roundoff regime. The remaining difference is the coarse diagonal, which classic builds as the negated row sum without touching D while this contracts D. Co-Authored-By: Claude Opus 5 (cherry picked from commit cfde17bac23f5a93e9e3c8b2a425f31f06e6f05b) --- .../linalg/system/ISTLSolverSystem.hpp | 41 ++++++++++++++----- 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 7043ebebed4..565d3afd010 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -31,6 +31,7 @@ #include #include +#include #include #include #include @@ -259,6 +260,15 @@ class ISTLSolverSystem : public ISTLSolver // than inside the preconditioner so that the linear-solver core never sees // anything well-specific, and so that this can later be replaced by a // value obtained from the well model without touching the core. + // Merged well block row -> well index. + std::size_t wellOfBlock(const std::size_t blockRow) const + { + const auto& off = wellLayout_.wellBlockOffsets; + const auto it = std::upper_bound(off.begin(), off.end(), blockRow); + assert(it != off.begin() && it != off.end()); + return static_cast(std::distance(off.begin(), it) - 1); + } + WellVector computeWellWeights(const ResVector& resWeights) const { const std::size_t numBlocks = mergedD_.N(); @@ -275,18 +285,29 @@ class ISTLSolverSystem : public ISTLSolver continue; } - if (wellWeightType_ == "cellavg") { - // The classic CPRW default (use_well_weights = false): average - // the reservoir weights over the cells this block row - // perforates, and use them on the conservation equations only. - // The control equation gets weight zero. + if (wellWeightType_ == "cellavg" || wellWeightType_ == "cellblockavg") { + // The classic CPRW weighting (use_well_weights = false): + // average the reservoir weights over perforated cells and use + // them on the conservation equations only, weight zero on the + // control equation. + // + // "cellavg" averages over every perforation of the whole well + // and gives every block of that well the same weights, which is + // what MultisegmentWellEquations::extractCPRPressureMatrix + // does. "cellblockavg" averages per block row instead, which + // is a finer but non-classic variant. + const bool perWell = (wellWeightType_ == "cellavg"); + const std::size_t first = perWell ? wellLayout_.firstBlock(wellOfBlock(wb)) : wb; + const std::size_t last = perWell ? wellLayout_.endBlock(wellOfBlock(wb)) : wb + 1; int nperf = 0; - for (auto col = mergedB_[wb].begin(), end = mergedB_[wb].end(); col != end; ++col) { - const auto& cw = resWeights[col.index()]; - for (int i = 0; i < numResDofs; ++i) { - lambda[i] += cw[i]; + for (std::size_t b = first; b < last; ++b) { + for (auto col = mergedB_[b].begin(), end = mergedB_[b].end(); col != end; ++col) { + const auto& cw = resWeights[col.index()]; + for (int i = 0; i < numResDofs; ++i) { + lambda[i] += cw[i]; + } + ++nperf; } - ++nperf; } if (nperf > 0) { for (int i = 0; i < numResDofs; ++i) { From d57761e82318d202389d52b60f83011fd5f45d97 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 19:17:08 +0200 Subject: [PATCH 07/22] Optionally give pressure-controlled wells a trivial coarse equation StandardWellEquations::extractCPRPressureMatrix gives a pressure-controlled well a unit diagonal and skips its B and C contributions, so its coarse equation is dp = 0 rather than a contracted well equation. The system CPRW stage contracted the matrix regardless of control mode. Add preconditioner.well_identity_on_pressure_control, on by default, doing the same. Which wells are pressure controlled is decided in the outer layer and handed down as one flag per well on the WellDofLayout, so the preconditioner still sees nothing but plain data. On SPE1 this takes the system solver from 459 to 440 linear iterations against 443 for the classic cprw. It does not, however, explain why the two diverge: the per-iteration residual norms agree to six digits for the first eight and then differ by ~5e-4, at a solve where no well is pressure controlled, and that divergence is unchanged by this option, by fixing the AMG setup with --cpr-reuse-setup=3, and by using ILU0 rather than UMFPack for the well solve. Co-Authored-By: Claude Opus 5 (cherry picked from commit 35dfe3b942db2c31fda61888b9c795299325c6ca) --- opm/simulators/linalg/setupPropertyTree.cpp | 3 +++ .../linalg/system/ISTLSolverSystem.hpp | 20 ++++++++++++++++++- .../linalg/system/SystemCprwPressureStage.hpp | 9 ++++++++- opm/simulators/linalg/system/SystemTypes.hpp | 18 +++++++++++++++++ 4 files changed, 48 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 46b743439ba..4d1b0a74aea 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -547,6 +547,9 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) // the only remaining difference is numerics // Only read when add_wells. prm.put("preconditioner.well_transfer", "no_prolongation"s); + // Give a pressure-controlled well a trivial coarse equation, matching + // StandardWellEquations::extractCPRPressureMatrix. Only read when add_wells. + prm.put("preconditioner.well_identity_on_pressure_control", "true"s); // --- Reservoir smoother --- prm.put("preconditioner.reservoir_smoother.maxiter", 1); diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 565d3afd010..34d68636751 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -208,7 +208,11 @@ class ISTLSolverSystem : public ISTLSolver const bool needStructureRefresh = !sysInitialized_ || globalStructureChanged; const auto& prm = this->prm_[this->activeSolverNum_]; - wellWeightType_ = prm.get("preconditioner.well_weight_type", std::string{"quasiimpes"}); + wellWeightType_ = prm.get("preconditioner.well_weight_type", std::string{"cellavg"}); + // Give a pressure-controlled well a trivial coarse equation, as the + // classic CPRW does. Off keeps the contracted equation for every well. + wellLayout_.identityOnPressureControl + = prm.get("preconditioner.well_identity_on_pressure_control", false); if (needStructureRefresh) { OPM_TIMEBLOCK(flexibleSolverCreate); @@ -253,6 +257,20 @@ class ISTLSolverSystem : public ISTLSolver total += d.N(); offsets.push_back(total); } + + // Which wells are on pressure control. Asking this is the outer + // layer's job; below here it is just a flag per well. The order + // matches addBCDMatrix, which walks the same well container. + wellLayout_.pressureControlled.clear(); + if (wellLayout_.identityOnPressureControl) { + const auto& wellModel = this->simulator_.problem().wellModel(); + const auto& wellState = wellModel.wellState(); + wellLayout_.pressureControlled.reserve(wellDMatrices_.size()); + for (const auto& well : wellModel) { + wellLayout_.pressureControlled.push_back( + well->isPressureControlled(wellState) ? 1 : 0); + } + } } // Weights used to contract each well's equations down to the single scalar diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp index de4e5b4136f..d2ff36cf4bb 100644 --- a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -385,7 +385,7 @@ class SystemCprwPressureStage const auto& bw = w0[c]; for (auto col = C[c].begin(), colEnd = C[c].end(); col != colEnd; ++col) { const auto j = wellOfBlock(col.index()); - if (!j.has_value()) { + if (!j.has_value() || layout.isPressureControlled(*j)) { continue; } Scalar el = 0.0; @@ -399,6 +399,13 @@ class SystemCprwPressureStage // Well rows. for (std::size_t j = 0; j < layout.numWells(); ++j) { const std::size_t wdof = numRes + j; + if (layout.isPressureControlled(j)) { + // A pressure-controlled well has a trivial coarse equation: + // its bhp is prescribed, so the coarse system carries dp = 0 + // rather than a contracted well equation. + (*coarseMatrix_)[wdof][wdof] = 1.0; + continue; + } for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { const auto& lw = w1[wb]; diff --git a/opm/simulators/linalg/system/SystemTypes.hpp b/opm/simulators/linalg/system/SystemTypes.hpp index 7feeb63fa48..f8cfe665cd6 100644 --- a/opm/simulators/linalg/system/SystemTypes.hpp +++ b/opm/simulators/linalg/system/SystemTypes.hpp @@ -90,6 +90,24 @@ struct WellDofLayout // generalising later is a change in the outer layer only. int pressureDofIndex = numWellDofs - 1; + // Per well: is it currently on pressure (bhp/thp) control? Filled in the + // outer layer, which is the only place that can ask. When + // identityOnPressureControl is set, such a well gets a trivial coarse + // equation instead of a contracted one -- what the classic CPRW does in + // StandardWellEquations::extractCPRPressureMatrix, where a + // pressure-controlled well is given a unit diagonal and its B and C + // contributions are skipped. Empty means "nothing is pressure + // controlled", so the flag is safe to leave unset. + std::vector pressureControlled; + bool identityOnPressureControl = false; + + bool isPressureControlled(const std::size_t j) const + { + return identityOnPressureControl + && j < pressureControlled.size() + && pressureControlled[j] != 0; + } + std::size_t numWells() const { return wellBlockOffsets.empty() ? 0 : wellBlockOffsets.size() - 1; From 6cf3c318de49c562ee92507ba4664a04109a7428 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 19:32:59 +0200 Subject: [PATCH 08/22] Add the classic coarse diagonal convention for multisegment wells Classic CPRW uses two conventions: StandardWellEquations contracts D for the coarse diagonal, while MultisegmentWellEquations sets it to minus the sum of the well row's reservoir entries and never reads D. The system stage always contracted D, which made the multisegment coarse matrix differ from the very first preconditioner application. Add preconditioner.well_coarse_diagonal with auto (contract D for single-block wells, row sum for multisegment ones, i.e. classic), contract_d and row_sum. Also dump the coarse matrix and right-hand side when verbosity exceeds 10, mirroring the classic path, so the two coarse systems can be diffed directly. On SPE1 with every well converted to multisegment this brings the system solver to 436 linear iterations, exactly the classic cprw count, with the per-iteration residual norms agreeing for ten lines and first differing at line 11 -- the same point at which the well-free cpr pair diverges. The multisegment formulation difference is therefore gone and only roundoff remains. Default stays contract_d pending a decision on which is wanted generally. Co-Authored-By: Claude Opus 5 (cherry picked from commit 15a984a44ed62c2f56962e6bf3b6e9d60f38d34e) --- opm/simulators/linalg/setupPropertyTree.cpp | 6 ++ .../linalg/system/SystemCprwPressureStage.hpp | 78 ++++++++++++++++++- .../linalg/system/SystemPreconditioner.hpp | 5 +- 3 files changed, 87 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 4d1b0a74aea..520900b2231 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -550,6 +550,12 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) // Give a pressure-controlled well a trivial coarse equation, matching // StandardWellEquations::extractCPRPressureMatrix. Only read when add_wells. prm.put("preconditioner.well_identity_on_pressure_control", "true"s); + // How a well's coarse diagonal is formed: + // auto - contract D for single-block wells, minus the row sum for + // multisegment ones, i.e. what classic cprw does + // contract_d - always contract D + // row_sum - always minus the row sum + prm.put("preconditioner.well_coarse_diagonal", "contract_d"s); // --- Reservoir smoother --- prm.put("preconditioner.reservoir_smoother.maxiter", 1); diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp index d2ff36cf4bb..0fe81e5c8dd 100644 --- a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -30,7 +30,12 @@ #include +#include + +#include + #include +#include #include #include #include @@ -84,6 +89,27 @@ enum class WellTransfer Classic, // neither -- as in PressureBhpTransferPolicy }; +// How the coarse diagonal of a well equation is formed. Classic CPRW uses +// two different conventions: StandardWellEquations contracts D, while +// MultisegmentWellEquations sets the diagonal to minus the sum of the well +// row's reservoir entries and never reads D at all. +enum class WellCoarseDiagonal +{ + Auto, // contract D for single-block wells, row sum for multi-block: as classic + ContractD, // always contract D + RowSum, // always minus the row sum +}; + +inline WellCoarseDiagonal wellCoarseDiagonalFromString(const std::string& name) +{ + if (name == "auto") { return WellCoarseDiagonal::Auto; } + if (name == "contract_d") { return WellCoarseDiagonal::ContractD; } + if (name == "row_sum") { return WellCoarseDiagonal::RowSum; } + OPM_THROW(std::invalid_argument, + "Unknown well_coarse_diagonal '" + name + + "'. Valid values are 'auto', 'contract_d' and 'row_sum'."); +} + inline WellTransfer wellTransferFromString(const std::string& name) { if (name == "full") { @@ -116,12 +142,16 @@ class SystemCprwPressureStage const PropertyTree& coarseSolverPrm, const int pressureIndex, const WellTransfer wellTransfer = WellTransfer::Full, - const Comm* comm = nullptr) + const Comm* comm = nullptr, + const WellCoarseDiagonal diagonal = WellCoarseDiagonal::ContractD, + const int verbosity = 0) : S_(S) , prm_(coarseSolverPrm) , pressureIndex_(pressureIndex) , wellTransfer_(wellTransfer) , comm_(comm) + , diagonal_(diagonal) + , verbosity_(verbosity) { } @@ -148,6 +178,7 @@ class SystemCprwPressureStage coarseRhs_.resize(coarseMatrix_->N()); coarseSol_.resize(coarseMatrix_->M()); + dumpCoarseMatrix(); } // (Re)create the coarse system and the solver acting on it. Must be @@ -189,6 +220,7 @@ class SystemCprwPressureStage OPM_TIMEBLOCK(systemCprwApply); moveToCoarseLevel(dRes, dWell, weights); + dumpCoarseRhs(); coarseSol_ = 0.0; Dune::InverseOperatorResult result; coarseSolver_->apply(coarseSol_, coarseRhs_, result); @@ -406,6 +438,12 @@ class SystemCprwPressureStage (*coarseMatrix_)[wdof][wdof] = 1.0; continue; } + const bool rowSumDiag + = (diagonal_ == WellCoarseDiagonal::RowSum) + || (diagonal_ == WellCoarseDiagonal::Auto + && layout.endBlock(j) - layout.firstBlock(j) > 1); + Scalar rowSum = 0.0; + for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { const auto& lw = w1[wb]; @@ -417,6 +455,13 @@ class SystemCprwPressureStage el += lw[i] * (*col)[i][p]; } (*coarseMatrix_)[wdof][col.index()] += el; + rowSum += el; + } + + if (rowSumDiag) { + // The classic multisegment convention takes the diagonal + // from the row sum and never reads D. + continue; } // Well row, well columns: @@ -444,6 +489,9 @@ class SystemCprwPressureStage // can happen for a well with no local perforations, or when the // weights annihilate the pressure column. Regularise to a unit row // rather than handing a singular system to AMG. + if (rowSumDiag) { + (*coarseMatrix_)[wdof][wdof] = -rowSum; + } auto& diag = (*coarseMatrix_)[wdof][wdof][0][0]; if (!(std::abs(diag) > 0.0)) { diag = 1.0; @@ -528,6 +576,32 @@ class SystemCprwPressureStage } } + // Same convention as the classic path: verbosity above 10 writes the + // coarse system out so the two can be compared entry by entry. + void dumpCoarseMatrix() const + { + if (verbosity_ <= 10) { + return; + } + static int counter = 0; + std::ofstream out("system_cprw_coarse_" + std::to_string(counter++) + ".mm"); + if (out) { + Dune::writeMatrixMarket(*coarseMatrix_, out); + } + } + + void dumpCoarseRhs() const + { + if (verbosity_ <= 10) { + return; + } + static int counter = 0; + std::ofstream out("system_cprw_rhs_" + std::to_string(counter++) + ".mm"); + if (out) { + Dune::writeMatrixMarket(coarseRhs_, out); + } + } + // Merged well block row -> well index. Linear scan is fine: the offsets // are sorted and this is only used while walking sparse rows of the well // blocks, which are short. @@ -546,6 +620,8 @@ class SystemCprwPressureStage int pressureIndex_ = 0; WellTransfer wellTransfer_ = WellTransfer::Full; const Comm* comm_ = nullptr; + WellCoarseDiagonal diagonal_ = WellCoarseDiagonal::ContractD; + int verbosity_ = 0; std::shared_ptr coarseComm_; std::shared_ptr coarseMatrix_; diff --git a/opm/simulators/linalg/system/SystemPreconditioner.hpp b/opm/simulators/linalg/system/SystemPreconditioner.hpp index 79e5651062a..2b4f83fc12c 100644 --- a/opm/simulators/linalg/system/SystemPreconditioner.hpp +++ b/opm/simulators/linalg/system/SystemPreconditioner.hpp @@ -326,8 +326,11 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdate(S_, coarseprm, pressureIndex_, - wellTransfer, resComm_); + wellTransfer, resComm_, diagonal, + prm.get("verbosity", 0)); cprwStage_->buildStructure(weights_); } else { if constexpr (isParallel) { From 3e8ac72677b10d3fea39d704a0f7450bf23fb1aa Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 4 Aug 2026 22:29:30 +0200 Subject: [PATCH 09/22] Dump the fine system on the system-solver path too ISTLSolverSystem::solve() overrides ISTLSolver::solve() and dropped the verbosity > 10 writeSystem block along with it, so --linear-solver=system_cpr silently wrote no reports/ dump. Restore it. Co-Authored-By: Claude Opus 5 (cherry picked from commit 397825c4111f6536eaf5897d74643b8688818afd) --- opm/simulators/linalg/system/ISTLSolverSystem.hpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 34d68636751..ddac29efe59 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -106,6 +106,15 @@ class ISTLSolverSystem : public ISTLSolver OPM_TIMEBLOCK(istlSolverSolve); ++this->solveCount_; + // Same fine-system dump as ISTLSolver::solve(), which this overrides, + // so the reservoir matrix and rhs can be diffed against the classic path. + if (this->prm_[this->activeSolverNum_].get("verbosity", 0) > 10) { + Helper::writeSystem(this->simulator_, + this->getMatrix(), + *Parent::rhs_, + this->comm_.get()); + } + const std::size_t numRes = Parent::matrix_->N(); const std::size_t numWell = cachedWellStructure_.totalWellBlocks; From ace10319fc5e2ee633fb009ed90738a1a2439dbd Mon Sep 17 00:00:00 2001 From: hnil Date: Wed, 5 Aug 2026 01:38:03 +0200 Subject: [PATCH 10/22] Default the CPRW pressure stage to the classic well transfer A settings sweep on full Norne, one factor at a time off the shipped defaults, serial, at the deck's own outer solver settings (maxiter 20, tol 5e-3): standard wells one segment per connection cprw 2262 2716 system_cpr 4303 4839 system_cprw, shipped defaults 2247 2569 well_transfer = classic 2248 2558 well_transfer = full 2246 2576 well_coarse_diagonal = row_sum 2242 2580 well_weight_type = cellblockavg 2247 2604 well_weight_type = quasiimpes 2245 5728 well_identity_on_pressure_control=0 2301 2574 Only well_transfer is worth changing, and only for multisegment wells: 2558 against 2569. Everything else is either flat (standard wells span 2242-2248) or worse. The margin is thin and comes from a case with an exact well solve, which nearly annihilates the well residual there is to restrict -- the other two modes stay, and are the ones to try if the well solve is ever made inexact. cellblockavg and quasiimpes are the two that must not become the default: 2604 and 5728 against 2569. Co-Authored-By: Claude Opus 5 (cherry picked from commit bd08be1a537281f612f57bca27fd046fb51514be) --- opm/simulators/linalg/setupPropertyTree.cpp | 29 +++++++++++++++------ 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 520900b2231..db459156b4b 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -532,21 +532,31 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) prm.put("preconditioner.type", "system_cpr"s); // How the well equations are contracted to the one coarse unknown each // well carries in the CPRW pressure system. Only read when add_wells. - // cellavg - average of the reservoir weights over the well's - // perforated cells, on the conservation equations only. - // This is what cprw does (use_well_weights = false) and, - // together with the trueimpes reservoir weights below, - // makes the pressure stage match the standard solver. - // quasiimpes - inv(D)^T e_bhp, normalised. - // unit - the pressure row as-is; a debugging baseline. + // cellavg - average of the reservoir weights over the well's + // perforated cells, on the conservation equations only. + // This is what cprw does (use_well_weights = false). + // cellblockavg - the same average taken per block row instead of per + // well. Identical to cellavg for a standard well, and + // worse for multisegment wells (Norne per-connection: + // 2604 against 2569). + // quasiimpes - inv(D)^T e_bhp, normalised. Catastrophic for + // multisegment wells; do not make it the default again + // without re-checking them. + // unit - the pressure row as-is; a debugging baseline. prm.put("preconditioner.well_weight_type", "cellavg"s); // How the well unknowns take part in the pressure-stage transfer: // full - restrict the well residual, prolong the bhp correction // no_prolongation - restrict, but discard the bhp correction // classic - neither, i.e. the classic cprw formulation, so that // the only remaining difference is numerics + // classic is the default: on full Norne with one segment per connection it + // measures best (2558, against 2569 for no_prolongation and 2576 for full), + // and on standard wells the three are within one iteration of each other. + // The margin is thin. It is also taken with an exact well solve, which + // nearly annihilates the well residual; restricting it may well pay once + // the well solve is inexact. // Only read when add_wells. - prm.put("preconditioner.well_transfer", "no_prolongation"s); + prm.put("preconditioner.well_transfer", "classic"s); // Give a pressure-controlled well a trivial coarse equation, matching // StandardWellEquations::extractCPRPressureMatrix. Only read when add_wells. prm.put("preconditioner.well_identity_on_pressure_control", "true"s); @@ -555,6 +565,9 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) // multisegment ones, i.e. what classic cprw does // contract_d - always contract D // row_sum - always minus the row sum + // contract_d is the default. On full Norne with one segment per connection + // it is worth ~0.4% over row_sum (2569 against 2580), and it is what makes + // the classic cprw path 2646 rather than 2716 on the same case. prm.put("preconditioner.well_coarse_diagonal", "contract_d"s); // --- Reservoir smoother --- From a61a260045b020715b43aff03113ffa046e9a280 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 14:09:26 +0200 Subject: [PATCH 11/22] Add composable sweep steps for the system preconditioner Each step reads the residual, corrects one block, and subtracts that correction from both halves of the residual, so a step is appendable to any other and composing them is a multiplicative sweep. SystemSweepPreconditioner runs an ordered list of them as one Dune preconditioner, carrying its defect internally -- which is what a two-level method needs, since Dune hands a smoother a single (correction, defect) pair. The residual is maintained incrementally rather than recomputed as d - S v. The two are equal in exact arithmetic but not bit for bit, and the incremental form is the one the existing three-stage SystemPreconditioner uses. Co-Authored-By: Claude Opus 5 (cherry picked from commit 51d443f04a4bff27d172005d594ab7bf05fbe680) --- CMakeLists_files.cmake | 1 + .../system/SystemPreconditionerParts.hpp | 416 ++++++++++++++++++ 2 files changed, 417 insertions(+) create mode 100644 opm/simulators/linalg/system/SystemPreconditionerParts.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 1531b7adb68..2e69f3f322d 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1129,6 +1129,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/Preconditioner2InverseOperator.hpp opm/simulators/linalg/system/MultiComm.hpp opm/simulators/linalg/system/SystemCprwPressureStage.hpp + opm/simulators/linalg/system/SystemPreconditionerParts.hpp opm/simulators/linalg/system/SystemPreconditioner.hpp opm/simulators/linalg/system/SystemPreconditionerFactory.hpp opm/simulators/linalg/system/SystemTypes.hpp diff --git a/opm/simulators/linalg/system/SystemPreconditionerParts.hpp b/opm/simulators/linalg/system/SystemPreconditionerParts.hpp new file mode 100644 index 00000000000..c3325b684a3 --- /dev/null +++ b/opm/simulators/linalg/system/SystemPreconditionerParts.hpp @@ -0,0 +1,416 @@ +/* + Copyright Equinor ASA 2026 + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_SYSTEMPRECONDITIONERPARTS_HEADER_INCLUDED +#define OPM_SYSTEMPRECONDITIONERPARTS_HEADER_INCLUDED + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include + +namespace Opm +{ + +// -------------------------------------------------------------------------- +// Building blocks for a preconditioner on the coupled system +// +// S = [ A C ] +// [ B D ] +// +// A sweep carries two things: the correction accumulated so far, and the +// residual that correction leaves. Each step reads the residual, produces a +// correction for one block, and subtracts that correction's effect from both +// halves of the residual. +// +// The residual is maintained incrementally rather than recomputed as d - S v. +// The two are equal in exact arithmetic but not bit for bit, and the point of +// this decomposition is that composing the steps reproduces the hand-written +// 3-stage SystemPreconditioner exactly -- see GeneralSystemPreconditioner. +// +// Composition is therefore multiplicative: every step sees what the step +// before it left behind. "Reservoir smoother followed by a well solve" is a +// block Gauss-Seidel sweep, not a block Jacobi one. +// -------------------------------------------------------------------------- + +template +struct SystemSweepState +{ + // Correction accumulated so far. + SystemVector* v = nullptr; + // Residual left by that correction, maintained incrementally. + SystemVector* res = nullptr; +}; + +template +class SystemSweepStep +{ +public: + virtual ~SystemSweepStep() = default; + + virtual void apply(SystemSweepState& state) = 0; + + // Refresh against changed matrix values, pattern unchanged. + virtual void update() = 0; + + // Rebuild against a changed well structure, which changes D's dimension. + // A changed well structure comes with changed matrix values, so a step + // that owns nothing sized by the wells still has to refresh itself; only + // steps that must be rebuilt outright override this. + virtual void rebuildForChangedWellStructure() + { + update(); + } +}; + +// -------------------------------------------------------------------------- +// Solve the well block: corr = D^-1 res[_1], v[_1] += corr, +// res[_0] -= C corr, res[_1] -= D corr. +// +// Appending this to anything -- a coarse pressure correction, a reservoir +// smoother -- cleans up the well equations that step disturbed. +// -------------------------------------------------------------------------- +template +class SystemWellStep : public SystemSweepStep +{ +public: + using WellOperator = Dune::MatrixAdapter, WellVector, WellVector>; + using WellSolver = Dune::FlexibleSolver; + + static constexpr auto _0 = Dune::Indices::_0; + static constexpr auto _1 = Dune::Indices::_1; + + SystemWellStep(const SystemMatrix& S, const PropertyTree& prm) + : S_(S) + , prm_(prm) + { + rebuildForChangedWellStructure(); + } + + void rebuildForChangedWellStructure() override + { + wop_ = std::make_unique(*S_.D); + std::function()> noWeights; + solver_ = std::make_unique(*wop_, prm_, noWeights, dummyPressureIndex); + rhs_.resize(S_.D->N()); + corr_.resize(S_.D->N()); + } + + void update() override + { + solver_->preconditioner().update(); + } + + void apply(SystemSweepState& state) override + { + auto& v = *state.v; + auto& res = *state.res; + + rhs_ = res[_1]; + corr_ = 0.0; + Dune::InverseOperatorResult result; + solver_->apply(corr_, rhs_, result); + + v[_1] += corr_; + S_.C->mmv(corr_, res[_0]); + S_.D->mmv(corr_, res[_1]); + } + +private: + // The well solver never uses a pressure index; pass something that would + // fail loudly if it ever did. + static constexpr std::size_t dummyPressureIndex = static_cast(-1); + + const SystemMatrix& S_; + PropertyTree prm_; + std::unique_ptr wop_; + std::unique_ptr solver_; + WellVector rhs_; + WellVector corr_; +}; + +// -------------------------------------------------------------------------- +// Apply a reservoir-only solver to the reservoir block: +// corr = R^-1 res[_0], v[_0] += corr, +// res[_0] -= A corr, res[_1] -= B corr. +// +// Lifts anything acting on a ResVector (ILU0, a CPR solve, ...) into a step on +// the coupled system. It leaves the well unknowns untouched; pair it with +// SystemWellStep for a full sweep. +// +// In parallel the residual handed to the solver is made consistent first; the +// correction is not, and the caller synchronises the accumulated reservoir +// correction once at the end instead. +// -------------------------------------------------------------------------- +template +class SystemReservoirStep : public SystemSweepStep +{ +public: + static constexpr bool isParallel = !std::is_same_v; + + using ResSolver = Dune::FlexibleSolver; + + static constexpr auto _0 = Dune::Indices::_0; + static constexpr auto _1 = Dune::Indices::_1; + + SystemReservoirStep(const SystemMatrix& S, + const PropertyTree& prm, + const std::function()>& weightsCalculator, + const int pressureIndex, + const ResComm* comm = nullptr) + : S_(S) + { + if constexpr (isParallel) { + comm_ = comm; + rop_ = std::make_unique(*S_.A, *comm_); + solver_ = std::make_unique(*rop_, *comm_, prm, + weightsCalculator, pressureIndex); + } else { + rop_ = std::make_unique(*S_.A); + solver_ = std::make_unique(*rop_, prm, + weightsCalculator, pressureIndex); + } + rhs_.resize(S_.A->N()); + corr_.resize(S_.A->N()); + } + + void update() override + { + solver_->preconditioner().update(); + } + + void apply(SystemSweepState& state) override + { + auto& v = *state.v; + auto& res = *state.res; + + rhs_ = res[_0]; + if constexpr (isParallel) { + comm_->copyOwnerToAll(rhs_, rhs_); + } + + corr_ = 0.0; + Dune::InverseOperatorResult result; + solver_->apply(corr_, rhs_, result); + + v[_0] += corr_; + S_.A->mmv(corr_, res[_0]); + S_.B->mmv(corr_, res[_1]); + } + +private: + const SystemMatrix& S_; + const ResComm* comm_ = nullptr; + std::unique_ptr rop_; + std::unique_ptr solver_; + ResVector rhs_; + ResVector corr_; +}; + +// -------------------------------------------------------------------------- +// The CPRW pressure stage as a sweep step: restrict the coupled residual to +// the scalar pressure system, solve it, prolong back. +// +// Whether the coarse correction reaches the well unknowns at all is the +// stage's own business (well_transfer); when it does not, there is nothing to +// subtract from the residual through C and D. +// -------------------------------------------------------------------------- +template +class SystemCprwStep : public SystemSweepStep +{ +public: + static constexpr bool isParallel = !std::is_same_v; + + using Stage = SystemCprwPressureStage; + + static constexpr auto _0 = Dune::Indices::_0; + static constexpr auto _1 = Dune::Indices::_1; + + SystemCprwStep(const SystemMatrix& S, + const PropertyTree& coarseSolverPrm, + const std::function()>& weightsCalculator, + const int pressureIndex, + const WellTransfer wellTransfer, + const WellCoarseDiagonal diagonal, + const int verbosity, + const ResComm* comm = nullptr) + : S_(S) + , weightsCalculator_(weightsCalculator) + { + if (!weightsCalculator_) { + OPM_THROW(std::invalid_argument, + "The CPRW pressure stage (add_wells) needs a weights calculator, but " + "none was configured. Set the reservoir solver's weight_type."); + } + if constexpr (isParallel) { + comm_ = comm; + } + stage_ = std::make_unique(S_, coarseSolverPrm, pressureIndex, + wellTransfer, comm_, diagonal, verbosity); + weights_ = weightsCalculator_(); + stage_->buildStructure(weights_); + rhs_.resize(S_.A->N()); + corrRes_.resize(S_.A->N()); + corrWell_.resize(S_.D->N()); + } + + void update() override + { + weights_ = weightsCalculator_(); + stage_->update(weights_); + } + + void rebuildForChangedWellStructure() override + { + // The coarse system carries one unknown per well, so a changed well + // structure changes its dimension and pattern. + weights_ = weightsCalculator_(); + stage_->buildStructure(weights_); + corrWell_.resize(S_.D->N()); + } + + void apply(SystemSweepState& state) override + { + auto& v = *state.v; + auto& res = *state.res; + + rhs_ = res[_0]; + if constexpr (isParallel) { + comm_->copyOwnerToAll(rhs_, rhs_); + } + + corrRes_ = 0.0; + corrWell_ = 0.0; + stage_->apply(rhs_, res[_1], weights_, corrRes_, corrWell_); + + v[_0] += corrRes_; + S_.A->mmv(corrRes_, res[_0]); + S_.B->mmv(corrRes_, res[_1]); + + if (stage_->prolongatesWellPressure()) { + v[_1] += corrWell_; + S_.C->mmv(corrWell_, res[_0]); + S_.D->mmv(corrWell_, res[_1]); + } + } + +private: + const SystemMatrix& S_; + const ResComm* comm_ = nullptr; + std::function()> weightsCalculator_; + SystemVector weights_; + std::unique_ptr stage_; + ResVector rhs_; + ResVector corrRes_; + WellVector corrWell_; +}; + + +// -------------------------------------------------------------------------- +// An ordered list of steps applied as one multiplicative sweep, presented as a +// Dune preconditioner so it can serve as the fine smoother of a two-level +// method. +// +// Dune hands a smoother a single (correction, defect) pair, so a sweep of more +// than one block has to carry its own defect internally. That is what this +// does: it starts from the defect it is given, and each step reduces it. +// -------------------------------------------------------------------------- +template +class SystemSweepPreconditioner + : public Dune::PreconditionerWithUpdate, SystemVector> +{ +public: + static constexpr auto _0 = Dune::Indices::_0; + static constexpr auto _1 = Dune::Indices::_1; + + explicit SystemSweepPreconditioner(std::vector>> steps, + const bool parallel = false) + : steps_(std::move(steps)) + , parallel_(parallel) + { + } + + void pre(SystemVector&, SystemVector&) override {} + void post(SystemVector&) override {} + + Dune::SolverCategory::Category category() const override + { + return parallel_ ? Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential; + } + + bool hasPerfectUpdate() const override + { + return true; + } + + void update() override + { + for (const auto& step : steps_) { + step->update(); + } + } + + void rebuildForChangedWellStructure() + { + for (const auto& step : steps_) { + step->rebuildForChangedWellStructure(); + } + } + + bool empty() const + { + return steps_.empty(); + } + + void apply(SystemVector& v, const SystemVector& d) override + { + v[_0].resize(d[_0].size()); + v[_1].resize(d[_1].size()); + v[_0] = 0.0; + v[_1] = 0.0; + res_[_0] = d[_0]; + res_[_1] = d[_1]; + + SystemSweepState state{&v, &res_}; + for (const auto& step : steps_) { + step->apply(state); + } + } + +private: + std::vector>> steps_; + bool parallel_ = false; + SystemVector res_; +}; + +} // namespace Opm + +#endif // OPM_SYSTEMPRECONDITIONERPARTS_HEADER_INCLUDED From 35b9e005576250aa94337e30b894ffa49cb4d255 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 14:09:44 +0200 Subject: [PATCH 12/22] TwoLevelMethodCpr: only resize the work vectors when the type allows it pre() resizes u_ and rhs_, which a MultiTypeBlockVector cannot do, and pre() is a virtual override so it is instantiated whether or not it is called. That kept the two-level method from being used on a block system at all. Guarded with if constexpr. apply() assigns to both vectors, which sizes them, so nothing is lost for a fine level that cannot resize. Co-Authored-By: Claude Opus 5 (cherry picked from commit a4325932d9eef70d48bfe0f2649d4d4231d5a237) --- opm/simulators/linalg/twolevelmethodcpr.hh | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/opm/simulators/linalg/twolevelmethodcpr.hh b/opm/simulators/linalg/twolevelmethodcpr.hh index 972e68c272f..4620461aa19 100644 --- a/opm/simulators/linalg/twolevelmethodcpr.hh +++ b/opm/simulators/linalg/twolevelmethodcpr.hh @@ -486,11 +486,15 @@ public: void pre(FineDomainType& x, FineRangeType& b) { - if (x.dim() != u_.dim()) { - u_.resize(x.dim()); - } - if (b.dim() != rhs_.dim()) { - rhs_.resize(b.dim()); + // A MultiTypeBlockVector has no resize(); apply() assigns to u_ and rhs_, + // which sizes them, so there is nothing to do for such a fine level. + if constexpr (requires { u_.resize(x.dim()); rhs_.resize(b.dim()); }) { + if (x.dim() != u_.dim()) { + u_.resize(x.dim()); + } + if (b.dim() != rhs_.dim()) { + rhs_.resize(b.dim()); + } } smoother_->pre(x,b); } From c2a58997420c64d17241984a362ab51dad2e2b05 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 14:09:44 +0200 Subject: [PATCH 13/22] Add the general system preconditioner, a two-level method on the system SystemPressureBhpTransferPolicy presents the CPRW pressure system, of dimension nCells + nWells, as a Dune level transfer policy; PressureSolverPolicy builds its solver; TwoLevelMethodCpr drives them with preSteps = 0 and postSteps = 1. Everything after the coarse correction -- the coarse solver's own well solve and then the smoother's parts -- is composed into one SystemSweepPreconditioner as the fine smoother, so the sequence applied is coarse -> well -> reservoir smoother -> well which is what the fixed three-stage SystemPreconditioner does. The parts are named by the property tree, coarse_solver { reservoir_solver, well_solver } smoother { reservoir_smoother, well_solver } and every sub-tree is optional, so parts can be left out -- which the fixed sequence cannot express. Without add_wells there is no system-wide coarse space: the reservoir solve is an ordinary block step in the sweep and no two-level method is built. updateForChangedWellStructure rebuilds everything rather than refreshing in place. That is not cosmetic: refreshing a CPR reservoir solve keeps its existing hierarchy while rebuilding re-aggregates it, and the two are different preconditioners. Rebuilding is the conservative choice for now. Registered as preconditioner type general_system_cpr on the same three operator/communication combinations as system_cpr. Co-Authored-By: Claude Opus 5 (cherry picked from commit 4ca0973031ed420f77314ff35ca79e3274a39868) --- CMakeLists_files.cmake | 2 + .../system/GeneralSystemPreconditioner.hpp | 321 ++++++++++++++++++ .../linalg/system/ISTLSolverSystem.hpp | 20 +- .../system/SystemPreconditionerFactory.cpp | 27 ++ .../SystemPressureBhpTransferPolicy.hpp | 134 ++++++++ 5 files changed, 503 insertions(+), 1 deletion(-) create mode 100644 opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp create mode 100644 opm/simulators/linalg/system/SystemPressureBhpTransferPolicy.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2e69f3f322d..12b669ebe70 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -1129,6 +1129,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/Preconditioner2InverseOperator.hpp opm/simulators/linalg/system/MultiComm.hpp opm/simulators/linalg/system/SystemCprwPressureStage.hpp + opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp + opm/simulators/linalg/system/SystemPressureBhpTransferPolicy.hpp opm/simulators/linalg/system/SystemPreconditionerParts.hpp opm/simulators/linalg/system/SystemPreconditioner.hpp opm/simulators/linalg/system/SystemPreconditionerFactory.hpp diff --git a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp new file mode 100644 index 00000000000..6807829e66d --- /dev/null +++ b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp @@ -0,0 +1,321 @@ +/* + Copyright Equinor ASA 2026 + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_GENERALSYSTEMPRECONDITIONER_HEADER_INCLUDED +#define OPM_GENERALSYSTEMPRECONDITIONER_HEADER_INCLUDED + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include + +namespace Opm +{ + +// -------------------------------------------------------------------------- +// A preconditioner for the coupled (reservoir, well) system assembled from +// named parts: +// +// coarse_solver { reservoir_solver, well_solver } +// smoother { reservoir_smoother, well_solver } +// +// When the coarse solver is the CPRW pressure stage (add_wells), this is a +// genuine two-level method on the system and is run by Dune's +// TwoLevelMethodCpr: SystemPressureBhpTransferPolicy is the transfer, the +// scalar pressure system of dimension nCells + nWells is the coarse level, and +// everything after the coarse correction is the fine smoother, with +// preSteps = 0 and postSteps = 1. +// +// Dune hands a smoother one (correction, defect) pair, so the parts that +// follow the coarse correction -- the coarse solver's own well solve, then the +// smoother's parts -- are composed into a single SystemSweepPreconditioner +// that carries its defect internally. The sequence applied is therefore +// +// coarse -> well -> reservoir smoother -> well +// +// which is what SystemPreconditioner does, and the two agree bit for bit. +// That works because TwoLevelMethodCpr's bookkeeping is the same one the +// sweep uses: postsmooth() reduces the defect by applyscaleadd(-1, lhs, rhs), +// which on SystemMatrix accumulates A,C into the reservoir half and B,D into +// the well half -- the same operations in the same order as the hand-written +// version's four mmv calls. +// +// Without add_wells there is no system-wide coarse space: the reservoir solve +// is an ordinary block solve, so all the parts go into one sweep and no +// two-level method is built. +// -------------------------------------------------------------------------- +template +class GeneralSystemPreconditioner + : public Dune::PreconditionerWithUpdate, SystemVector> +{ +public: + static constexpr bool isParallel = !std::is_same_v; + + static constexpr auto _0 = Dune::Indices::_0; + static constexpr auto _1 = Dune::Indices::_1; + + using ReservoirStep = SystemReservoirStep; + using WellStep = SystemWellStep; + using Sweep = SystemSweepPreconditioner; + + // The fine level of the two-level method is the whole coupled system. + using FineOperator = Dune::MatrixAdapter, + SystemVector, + SystemVector>; + using TransferPolicy = SystemPressureBhpTransferPolicy; + using CoarseOperator = typename TransferPolicy::CoarseOperator; + using CoarseSolverPolicy = Dune::Amg::PressureSolverPolicy, + TransferPolicy>; + using TwoLevel = Dune::Amg::TwoLevelMethodCpr< + FineOperator, + CoarseSolverPolicy, + Dune::PreconditionerWithUpdate, SystemVector>>; + + GeneralSystemPreconditioner(const SystemMatrix& S, + const std::function()>& weightsCalculator, + int pressureIndex, + const PropertyTree& prm) + requires(!isParallel) + : S_(S) + , pressureIndex_(pressureIndex) + , prm_(prm) + , weightsCalculator_(weightsCalculator) + { + build(); + } + + GeneralSystemPreconditioner(const SystemMatrix& S, + const std::function()>& weightsCalculator, + int pressureIndex, + const PropertyTree& prm, + const ResComm& resComm) + requires(isParallel) + : S_(S) + , resComm_(&resComm) + , pressureIndex_(pressureIndex) + , prm_(prm) + , weightsCalculator_(weightsCalculator) + { + build(); + } + + void pre(SystemVector&, SystemVector&) override + { + // Not forwarded to the two-level method: its pre() resizes u_ and rhs_, + // and MultiTypeBlockVector has no resize. apply() assigns to both, which + // sizes them, so there is nothing for pre() to do here. + } + + void post(SystemVector&) override + { + } + + Dune::SolverCategory::Category category() const override + { + if constexpr (isParallel) { + return Dune::SolverCategory::overlapping; + } else { + return Dune::SolverCategory::sequential; + } + } + + bool hasPerfectUpdate() const override + { + return true; + } + + void update() override + { + if (twoLevel_) { + weights_ = weightsCalculator_(); + twoLevel_->updatePreconditioner(sweep_, *coarsePolicy_); + } else { + sweep_->update(); + } + } + + void updateForChangedWellStructure() + { + // A changed well structure changes D's dimension and the dimension of + // the coarse system, so everything is built again from the property + // tree. That is deliberately blunt: refreshing the reservoir solves in + // place instead would keep a CPR hierarchy rather than re-aggregate it, + // and the two are different preconditioners -- on Norne that choice + // alone moves general_system_cpr between 4303 and 4253 linear + // iterations. Rebuilding is the safe default until it is measured. + build(); + } + + void apply(SystemVector& v, const SystemVector& d) override + { + v[_0].resize(d[_0].size()); + v[_1].resize(d[_1].size()); + v[_0] = 0.0; + v[_1] = 0.0; + + if (twoLevel_) { + // TwoLevelMethodCpr accumulates into v, so it has to start at zero. + twoLevel_->apply(v, d); + } else { + sweep_->apply(v, d); + } + + // The steps leave the reservoir correction with stale overlap entries; + // make it consistent once, as the fixed three-stage version does. + if constexpr (isParallel) { + resComm_->copyOwnerToAll(v[_0], v[_0]); + } + } + +private: + const SystemMatrix& S_; + const ResComm* resComm_ = nullptr; + int pressureIndex_ = 0; + PropertyTree prm_; + std::function()> weightsCalculator_; + SystemVector weights_; + + std::shared_ptr sweep_; + std::unique_ptr fineOp_; + std::unique_ptr transfer_; + std::unique_ptr coarsePolicy_; + std::unique_ptr twoLevel_; + // Kept so the coarse level can be built again when the well count changes. + std::optional coarseResPrm_; + + void build() + { + twoLevel_.reset(); + coarsePolicy_.reset(); + transfer_.reset(); + fineOp_.reset(); + sweep_.reset(); + coarseResPrm_.reset(); + + // The weights arrive from the outer layer for the whole system; the + // reservoir-only sub-solvers want just their own part of them. + std::function()> resWeightCalc; + if (weightsCalculator_) { + const auto calc = weightsCalculator_; + resWeightCalc = [calc]() { return calc()[_0]; }; + } + + const auto coarse = prm_.get_child_optional("coarse_solver"); + const auto smoother = prm_.get_child_optional("smoother"); + if (!coarse && !smoother) { + OPM_THROW(std::invalid_argument, + "general_system_cpr needs at least one of the sub-trees " + "'coarse_solver' and 'smoother'."); + } + + std::vector>> steps; + const auto addReservoir = [&](const PropertyTree& p) { + steps.push_back(std::make_unique( + S_, p, resWeightCalc, pressureIndex_, resComm_)); + }; + const auto addWell = [&](const PropertyTree& p) { + steps.push_back(std::make_unique(S_, p)); + }; + + bool twoLevelCoarse = false; + if (coarse) { + if (const auto res = coarse->get_child_optional("reservoir_solver")) { + // add_wells is the same switch the classic CPR/CPRW pair uses: + // it promotes the pressure stage from reservoir-only CPR to + // CPRW over the full (reservoir, well) system, and only then + // is there a coarse space for the whole system. + if (res->get("preconditioner.add_wells", false)) { + coarseResPrm_ = *res; + buildCoarse(); + twoLevelCoarse = true; + } else { + addReservoir(*res); + } + } + if (const auto well = coarse->get_child_optional("well_solver")) { + addWell(*well); + } + } + if (smoother) { + if (const auto res = smoother->get_child_optional("reservoir_smoother")) { + addReservoir(*res); + } + if (const auto well = smoother->get_child_optional("well_solver")) { + addWell(*well); + } + } + + if (steps.empty() && !twoLevelCoarse) { + OPM_THROW(std::invalid_argument, + "general_system_cpr was configured with no parts at all."); + } + + sweep_ = std::make_shared(std::move(steps), isParallel); + + if (twoLevelCoarse) { + twoLevel_ = std::make_unique(*fineOp_, sweep_, *transfer_, + *coarsePolicy_, + /*preSteps=*/0, /*postSteps=*/1); + } + } + + void buildCoarse() + { + const auto& resPrm = *coarseResPrm_; + const auto coarsePrm = resPrm.get_child_optional("preconditioner.coarsesolver") + ? resPrm.get_child("preconditioner.coarsesolver") + : PropertyTree(); + const auto wellTransfer = wellTransferFromString( + prm_.get("well_transfer", std::string{"full"})); + const auto diagonal = wellCoarseDiagonalFromString( + prm_.get("well_coarse_diagonal", std::string{"contract_d"})); + + if (!weightsCalculator_) { + OPM_THROW(std::invalid_argument, + "The CPRW pressure stage (add_wells) needs a weights calculator, but " + "none was configured. Set the reservoir solver's weight_type."); + } + weights_ = weightsCalculator_(); + + fineOp_ = std::make_unique(S_); + transfer_ = std::make_unique(S_, weights_, coarsePrm, pressureIndex_, + wellTransfer, diagonal, + prm_.get("verbosity", 0), resComm_); + coarsePolicy_ = std::make_unique(coarsePrm); + } +}; + +} // namespace Opm + +#endif // OPM_GENERALSYSTEMPRECONDITIONER_HEADER_INCLUDED diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index ddac29efe59..1ec4b40d55e 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -20,6 +20,7 @@ #define OPM_ISTLSOLVERSYSTEM_HEADER_INCLUDED #include +#include #include #include #include @@ -179,8 +180,10 @@ class ISTLSolverSystem : public ISTLSolver using SysSolverType = Dune::InverseOperator, SystemVector>; using SysPrecondType = Dune::PreconditionerWithUpdate, SystemVector>; using SeqSysPrecondType = SystemPreconditioner>; + using SeqGeneralSysPrecondType = GeneralSystemPreconditioner>; #if HAVE_MPI using ParSysPrecondType = SystemPreconditioner, ParResComm>; + using ParGeneralSysPrecondType = GeneralSystemPreconditioner, ParResComm>; #endif SysSolverType* sysSolver_ = nullptr; SysPrecondType* sysPrecond_ = nullptr; @@ -395,6 +398,8 @@ class ISTLSolverSystem : public ISTLSolver if (this->comm_->communicator().size() > 1) { if (auto* precond = dynamic_cast(sysPrecond_)) { precond->updateForChangedWellStructure(); + } else if (auto* general = dynamic_cast(sysPrecond_)) { + general->updateForChangedWellStructure(); } else { // Rebuild the parallel solver if the parallel preconditioner cannot be updated in-place. createSystemSolver(prm); @@ -405,16 +410,29 @@ class ISTLSolverSystem : public ISTLSolver if (auto* precond = dynamic_cast(sysPrecond_)) { precond->updateForChangedWellStructure(); + } else if (auto* general = dynamic_cast(sysPrecond_)) { + general->updateForChangedWellStructure(); } else { // Rebuild the solver if the sequential preconditioner cannot be updated in-place createSystemSolver(prm); } } + // Where the reservoir sub-solver sits depends on the preconditioner: the + // fixed three-stage one keeps it at the top, the general one nests it under + // the coarse solver. Only the weights are read from it here. + Opm::PropertyTree reservoirSolverTree(const Opm::PropertyTree& prm) const + { + if (auto general = prm.get_child_optional("preconditioner.coarse_solver.reservoir_solver")) { + return *general; + } + return prm.get_child("preconditioner.reservoir_solver"); + } + void createSystemSolver(const Opm::PropertyTree& prm) { // Derive weights from the reservoir sub-block config (which uses CPR internally) - auto resSolverPrm = prm.get_child("preconditioner.reservoir_solver"); + auto resSolverPrm = reservoirSolverTree(prm); std::function()> resWeightCalc = this->getWeightsCalculator(resSolverPrm, this->getMatrix(), pressureIndex); diff --git a/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp b/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp index 07e8bcddaf7..94c68834b5c 100644 --- a/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp +++ b/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -41,6 +42,14 @@ void addSystemCprSeq() return std::make_shared>>( op.getmat(), sysWeightCalc, pressureIndex, prm); }); + + F::addCreator("general_system_cpr", + [](const O& op, const P& prm, + const std::function& sysWeightCalc, + std::size_t pressureIndex) { + return std::make_shared>>( + op.getmat(), sysWeightCalc, pressureIndex, prm); + }); } #if HAVE_MPI @@ -64,6 +73,14 @@ void addSystemCprParSeq() return std::make_shared>>( op.getmat(), sysWeightCalc, pressureIndex, prm); }); + + F::addCreator("general_system_cpr", + [](const O& op, const P& prm, + const std::function& sysWeightCalc, + std::size_t pressureIndex) { + return std::make_shared>>( + op.getmat(), sysWeightCalc, pressureIndex, prm); + }); } template @@ -83,6 +100,16 @@ void addSystemCprPar() return std::make_shared, Opm::ParResComm>>( op.getmat(), sysWeightCalc, pressureIndex, prm, resComm); }); + + F::addCreator("general_system_cpr", + [](const O& op, const P& prm, + const std::function& sysWeightCalc, + std::size_t pressureIndex, + const Opm::SystemComm& comm) { + const auto& resComm = comm[Dune::Indices::_0]; + return std::make_shared, Opm::ParResComm>>( + op.getmat(), sysWeightCalc, pressureIndex, prm, resComm); + }); } #endif diff --git a/opm/simulators/linalg/system/SystemPressureBhpTransferPolicy.hpp b/opm/simulators/linalg/system/SystemPressureBhpTransferPolicy.hpp new file mode 100644 index 00000000000..410ccd4da11 --- /dev/null +++ b/opm/simulators/linalg/system/SystemPressureBhpTransferPolicy.hpp @@ -0,0 +1,134 @@ +/* + Copyright Equinor ASA 2026 + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#ifndef OPM_SYSTEMPRESSUREBHPTRANSFERPOLICY_HEADER_INCLUDED +#define OPM_SYSTEMPRESSUREBHPTRANSFERPOLICY_HEADER_INCLUDED + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +namespace Opm +{ + +// -------------------------------------------------------------------------- +// The CPRW pressure system of the coupled (reservoir, well) system, presented +// as a Dune level transfer policy so that it can drive TwoLevelMethodCpr. +// +// The assembly and the transfers themselves live in SystemCprwPressureStage; +// this only adapts them to the interface the two-level machinery expects. The +// fine level is the whole system (a SystemVector), the coarse level is the +// scalar pressure system of dimension nCells + nWells. +// +// This is the counterpart of PressureBhpTransferPolicy for the system solver, +// with the essential difference that it needs no well model: the coarse system +// comes from the B/C/D blocks and the weights alone. +// -------------------------------------------------------------------------- +template +class SystemPressureBhpTransferPolicy + : public Dune::Amg::LevelTransferPolicyCpr> +{ +public: + using CoarseOperator = Details::CoarseOperatorType; + using ParentType = Dune::Amg::LevelTransferPolicyCpr; + using ParallelInformation = Comm; + using Stage = SystemCprwPressureStage; + + static constexpr auto _0 = Dune::Indices::_0; + static constexpr auto _1 = Dune::Indices::_1; + + SystemPressureBhpTransferPolicy(const SystemMatrix& S, + const SystemVector& weights, + const PropertyTree& coarseSolverPrm, + const int pressureIndex, + const WellTransfer wellTransfer, + const WellCoarseDiagonal diagonal, + const int verbosity, + const Comm* comm = nullptr) + : weights_(&weights) + , stage_(std::make_shared(S, coarseSolverPrm, pressureIndex, + wellTransfer, comm, diagonal, verbosity)) + { + } + + void createCoarseLevelSystem(const FineOperator&) override + { + stage_->buildCoarseSystem(*weights_); + const auto& coarse = stage_->coarseMatrixPtr(); + this->lhs_.resize(coarse->M()); + this->rhs_.resize(coarse->N()); + using OperatorArgs = typename Dune::Amg::ConstructionTraits::Arguments; + OperatorArgs oargs(coarse, stage_->coarseCommunication()); + this->operator_ = Dune::Amg::ConstructionTraits::construct(oargs); + } + + void calculateCoarseEntries(const FineOperator&) override + { + stage_->assembleCoarseEntries(*weights_); + } + + void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override + { + stage_->moveToCoarseLevel(fine[_0], fine[_1], *weights_, this->rhs_); + this->lhs_ = 0; + } + + void moveToFineLevel(typename ParentType::FineDomainType& fine) override + { + // The well half comes back zero when the stage does not prolong to it, + // so the two-level defect update subtracts C*0 and D*0 -- wasted + // arithmetic, but exactly no change, which is what lets this reproduce + // the hand-written sequence bit for bit. + stage_->moveToFineLevel(this->lhs_, fine[_0], fine[_1]); + } + + SystemPressureBhpTransferPolicy* clone() const override + { + return new SystemPressureBhpTransferPolicy(*this); + } + + const Comm& getCoarseLevelCommunication() const + { + return stage_->coarseCommunication(); + } + + Stage& stage() + { + return *stage_; + } + +private: + // Points at the weights the owning preconditioner refreshes, so an update + // is seen here without copying. + const SystemVector* weights_; + // Shared so that clone() is a shallow copy, as it is for + // PressureBhpTransferPolicy. + std::shared_ptr stage_; +}; + +} // namespace Opm + +#endif // OPM_SYSTEMPRESSUREBHPTRANSFERPOLICY_HEADER_INCLUDED From daba299805914e6f9d93df5e4341e7af3924711e Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 14:10:03 +0200 Subject: [PATCH 14/22] Add --linear-solver=general_system_cpr and general_system_cprw setupGeneralSystemCPR writes the same sub-solvers setupSystemCPR does, just under coarse_solver/ and smoother/ rather than side by side, so the two configurations are the same algorithm with the same settings. The sub-solvers are now written by shared helpers, so they cannot drift apart. The JSON validation for the general layout only checks the composition, since which parts are present is the point of it. Co-Authored-By: Claude Opus 5 (cherry picked from commit a5792b9c0e5dca2f9b703cbccf4e3c22c0aa14c2) --- .../linalg/FlowLinearSolverParameters.cpp | 1 + .../linalg/ISTLSolverRuntimeOptionProxy.hpp | 6 +- opm/simulators/linalg/setupPropertyTree.cpp | 225 +++++++++++++----- opm/simulators/linalg/setupPropertyTree.hpp | 2 + 4 files changed, 176 insertions(+), 58 deletions(-) diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.cpp b/opm/simulators/linalg/FlowLinearSolverParameters.cpp index 1e93160f58f..c5bf38702fd 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.cpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.cpp @@ -136,6 +136,7 @@ void FlowLinearSolverParameters::registerParameters() Parameters::Register ("Configuration of solver. Valid options are: cprw (default), system_cpr (CPU-only), " "system_cprw (CPU-only, system_cpr with the wells in the pressure stage), " + "general_system_cpr / general_system_cprw (the same composed from named parts), " "ilu0, dilu, cpr (an alias for cprw), cpr_quasiimpes, " "cpr_trueimpes, cpr_trueimpesanalytic, amg or hybrid (experimental). " "Alternatively, you can request a configuration to be read from a " diff --git a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp index 3c6cbe20e41..8cb975d48ab 100644 --- a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp +++ b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp @@ -155,13 +155,15 @@ class ISTLSolverRuntimeOptionProxy : public AbstractISTLSolver(); - bool useSystemCpr = (linSolverConf == "system_cpr") || (linSolverConf == "system_cprw"); + bool useSystemCpr = (linSolverConf == "system_cpr") || (linSolverConf == "system_cprw") + || (linSolverConf == "general_system_cpr") || (linSolverConf == "general_system_cprw"); if (!useSystemCpr && linSolverConf.size() > 5 && linSolverConf.ends_with(".json") && std::filesystem::exists(linSolverConf)) { try { PropertyTree prm(linSolverConf); - useSystemCpr = (prm.get("preconditioner.type", "") == "system_cpr"); + const auto type = prm.get("preconditioner.type", ""); + useSystemCpr = (type == "system_cpr") || (type == "general_system_cpr"); } catch (...) {} } if (useSystemCpr) { diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index db459156b4b..e05ddea7bc0 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -243,6 +243,18 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters } return setupSystemCPR(conf, p); } + + // The same algorithm composed from named parts rather than a fixed + // three-stage sequence; identical results with these settings. + if ((conf == "general_system_cpr") || (conf == "general_system_cprw")) { + if (!linearSolverMaxIterSet) { + p.linear_solver_maxiter_ = 20; + } + if (!linearSolverReductionSet) { + p.linear_solver_reduction_ = 0.005; + } + return setupGeneralSystemCPR(conf, p); + } } if (conf == "amg") { @@ -299,7 +311,8 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters else { OPM_THROW(std::invalid_argument, conf + " is not a valid setting for --linear-solver-configuration." - " Please use ilu0, dilu, isai, cpr, cprw, cpr_trueimpes, cpr_quasiimpes, cpr_trueimpesanalytic, or system_cpr"); + " Please use ilu0, dilu, isai, cpr, cprw, cpr_trueimpes, cpr_quasiimpes, cpr_trueimpesanalytic," + " system_cpr, system_cprw, general_system_cpr or general_system_cprw"); } @@ -515,21 +528,79 @@ setupUMFPack([[maybe_unused]] const std::string& conf, const FlowLinearSolverPar } -PropertyTree -setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) +namespace { + +// The sub-solvers of the system preconditioner, each written at a caller +// chosen key. Both the fixed three-stage layout and the general one use +// these, so the trees they produce are identical and the two preconditioners +// can be compared setting for setting. +void setupSystemReservoirSmoother(PropertyTree& prm, + const std::string& at, + const FlowLinearSolverParameters& p) { using namespace std::string_literals; - const bool add_wells = (conf == "system_cprw"); - PropertyTree prm; + prm.put(at + ".maxiter", 1); + prm.put(at + ".tol", p.linear_solver_reduction_); + prm.put(at + ".verbosity", 0); + // Apply the configured smoother once without loopsolver's extra defect + // updates and convergence bookkeeping. + prm.put(at + ".solver", "preconditioner2inverseoperator"s); + prm.put(at + ".preconditioner.type", "paroverilu0"s); + prm.put(at + ".preconditioner.relaxation", 1.0); +} - // Outer solver - prm.put("maxiter", p.linear_solver_maxiter_); - prm.put("tol", p.linear_solver_reduction_); - prm.put("verbosity", p.linear_solver_verbosity_); - prm.put("solver", getSolverString(p)); +void setupSystemReservoirSolver(PropertyTree& prm, + const std::string& at, + const FlowLinearSolverParameters& p, + const bool add_wells) +{ + using namespace std::string_literals; + prm.put(at + ".maxiter", 1); + prm.put(at + ".tol", p.linear_solver_reduction_); + prm.put(at + ".verbosity", 0); + // Apply the configured reservoir preconditioner once without loopsolver's + // extra defect updates and convergence bookkeeping. + prm.put(at + ".solver", "preconditioner2inverseoperator"s); + prm.put(at + ".preconditioner.type", "cpr"s); + prm.put(at + ".preconditioner.relaxation", 1.0); + prm.put(at + ".preconditioner.use_well_weights", "false"s); + // add_wells promotes the pressure stage from reservoir-only CPR to CPRW + // over the full (reservoir, well) system. + prm.put(at + ".preconditioner.add_wells", add_wells ? "true"s : "false"s); + prm.put(at + ".preconditioner.weight_type", "trueimpes"s); + prm.put(at + ".preconditioner.pre_smooth", 0); + prm.put(at + ".preconditioner.post_smooth", 0); + // Set unused finesmoother to jac to avoid spending time setuping an ILU smoother that won't be used. + prm.put(at + ".preconditioner.finesmoother.type", "jac"s); + prm.put(at + ".preconditioner.finesmoother.relaxation", 1.0); + prm.put(at + ".preconditioner.verbosity", 0); + prm.put(at + ".preconditioner.coarsesolver.maxiter", 1); + prm.put(at + ".preconditioner.coarsesolver.tol", 1e-1); + prm.put(at + ".preconditioner.coarsesolver.solver", "loopsolver"s); + prm.put(at + ".preconditioner.coarsesolver.verbosity", 0); + prm.put(at + ".preconditioner.coarsesolver.preconditioner.type", "amg"s); + setupDuneAMG(prm, at + ".preconditioner.coarsesolver.preconditioner."); +} - // Top-level preconditioner: system_cpr - prm.put("preconditioner.type", "system_cpr"s); +void setupSystemWellSolver(PropertyTree& prm, + const std::string& at, + const FlowLinearSolverParameters& p) +{ + using namespace std::string_literals; + prm.put(at + ".maxiter", 1); + prm.put(at + ".tol", p.linear_solver_reduction_); + prm.put(at + ".verbosity", 0); + prm.put(at + ".solver", "umfpack"s); +} + +} // anonymous namespace + +namespace { + +// The knobs the CPRW pressure stage reads, shared by both layouts. +void setupSystemCPRWellOptions(PropertyTree& prm) +{ + using namespace std::string_literals; // How the well equations are contracted to the one coarse unknown each // well carries in the CPRW pressure system. Only read when add_wells. // cellavg - average of the reservoir weights over the well's @@ -570,56 +641,98 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) // the classic cprw path 2646 rather than 2716 on the same case. prm.put("preconditioner.well_coarse_diagonal", "contract_d"s); - // --- Reservoir smoother --- - prm.put("preconditioner.reservoir_smoother.maxiter", 1); - prm.put("preconditioner.reservoir_smoother.tol", p.linear_solver_reduction_); - prm.put("preconditioner.reservoir_smoother.verbosity", 0); - // Apply the configured smoother once without loopsolver's extra defect - // updates and convergence bookkeeping. - prm.put("preconditioner.reservoir_smoother.solver", "preconditioner2inverseoperator"s); - prm.put("preconditioner.reservoir_smoother.preconditioner.type", "paroverilu0"s); - prm.put("preconditioner.reservoir_smoother.preconditioner.relaxation", 1.0); - - // --- Reservoir solver (CPR with AMG coarse solver) --- - prm.put("preconditioner.reservoir_solver.maxiter", 1); - prm.put("preconditioner.reservoir_solver.tol", p.linear_solver_reduction_); - prm.put("preconditioner.reservoir_solver.verbosity", 0); - // Apply the configured reservoir preconditioner once without loopsolver's - // extra defect updates and convergence bookkeeping. - prm.put("preconditioner.reservoir_solver.solver", "preconditioner2inverseoperator"s); - prm.put("preconditioner.reservoir_solver.preconditioner.type", "cpr"s); - prm.put("preconditioner.reservoir_solver.preconditioner.relaxation", 1.0); - prm.put("preconditioner.reservoir_solver.preconditioner.use_well_weights", "false"s); - // add_wells promotes the pressure stage from reservoir-only CPR to CPRW - // over the full (reservoir, well) system. - prm.put("preconditioner.reservoir_solver.preconditioner.add_wells", - add_wells ? "true"s : "false"s); - prm.put("preconditioner.reservoir_solver.preconditioner.weight_type", "trueimpes"s); - prm.put("preconditioner.reservoir_solver.preconditioner.pre_smooth", 0); - prm.put("preconditioner.reservoir_solver.preconditioner.post_smooth", 0); - // Set unused finesmoother to jac to avoid spending time setuping an ILU smoother that won't be used. - prm.put("preconditioner.reservoir_solver.preconditioner.finesmoother.type", "jac"s); - prm.put("preconditioner.reservoir_solver.preconditioner.finesmoother.relaxation", 1.0); - prm.put("preconditioner.reservoir_solver.preconditioner.verbosity", 0); - prm.put("preconditioner.reservoir_solver.preconditioner.coarsesolver.maxiter", 1); - prm.put("preconditioner.reservoir_solver.preconditioner.coarsesolver.tol", 1e-1); - prm.put("preconditioner.reservoir_solver.preconditioner.coarsesolver.solver", "loopsolver"s); - prm.put("preconditioner.reservoir_solver.preconditioner.coarsesolver.verbosity", 0); - prm.put("preconditioner.reservoir_solver.preconditioner.coarsesolver.preconditioner.type", "amg"s); - setupDuneAMG(prm, "preconditioner.reservoir_solver.preconditioner.coarsesolver.preconditioner."); - - // --- Well solver --- - prm.put("preconditioner.well_solver.maxiter", 1); - prm.put("preconditioner.well_solver.tol", p.linear_solver_reduction_); - prm.put("preconditioner.well_solver.verbosity", 0); - prm.put("preconditioner.well_solver.solver", "umfpack"s); +} + +} // anonymous namespace + +PropertyTree +setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) +{ + using namespace std::string_literals; + const bool add_wells = (conf == "system_cprw"); + PropertyTree prm; + + // Outer solver + prm.put("maxiter", p.linear_solver_maxiter_); + prm.put("tol", p.linear_solver_reduction_); + prm.put("verbosity", p.linear_solver_verbosity_); + prm.put("solver", getSolverString(p)); + + // Top-level preconditioner: system_cpr + prm.put("preconditioner.type", "system_cpr"s); + setupSystemCPRWellOptions(prm); + + setupSystemReservoirSmoother(prm, "preconditioner.reservoir_smoother"s, p); + setupSystemReservoirSolver(prm, "preconditioner.reservoir_solver"s, p, add_wells); + setupSystemWellSolver(prm, "preconditioner.well_solver"s, p); + return prm; +} + +// The same algorithm as setupSystemCPR, expressed as the parts the general +// preconditioner composes: +// +// coarse_solver { reservoir_solver, well_solver } +// smoother { reservoir_smoother, well_solver } +// +// applied in that order. The sub-trees are byte for byte the ones +// setupSystemCPR writes, and with this layout the general preconditioner +// reproduces the fixed three-stage one exactly. It exists so that the parts +// can be left out or repeated, which the fixed sequence cannot express. +PropertyTree +setupGeneralSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) +{ + using namespace std::string_literals; + const bool add_wells = (conf == "general_system_cprw"); + PropertyTree prm; + + prm.put("maxiter", p.linear_solver_maxiter_); + prm.put("tol", p.linear_solver_reduction_); + prm.put("verbosity", p.linear_solver_verbosity_); + prm.put("solver", getSolverString(p)); + + prm.put("preconditioner.type", "general_system_cpr"s); + setupSystemCPRWellOptions(prm); + + // Each stage gets the well solver the fixed layout shares between its + // stages, so the sequence is coarse -> well -> smoother -> well. + setupSystemReservoirSolver(prm, "preconditioner.coarse_solver.reservoir_solver"s, p, add_wells); + setupSystemWellSolver(prm, "preconditioner.coarse_solver.well_solver"s, p); + setupSystemReservoirSmoother(prm, "preconditioner.smoother.reservoir_smoother"s, p); + setupSystemWellSolver(prm, "preconditioner.smoother.well_solver"s, p); return prm; } +void validateGeneralSystemCPRTree(const PropertyTree& prm) +{ + // Only the composition is checked here: which parts exist is the point of + // this layout, so nearly everything is optional. The parts themselves are + // validated by whatever builds them. + if (!prm.get_child_optional("preconditioner.coarse_solver").has_value() + && !prm.get_child_optional("preconditioner.smoother").has_value()) { + OPM_THROW(std::invalid_argument, + "general_system_cpr JSON configuration needs at least one of the " + "'preconditioner.coarse_solver' and 'preconditioner.smoother' sub-trees."); + } + const auto coarse = prm.get_child_optional("preconditioner.coarse_solver.reservoir_solver"); + if (coarse && coarse->get("preconditioner.add_wells", false) + && !coarse->get_child_optional("preconditioner.coarsesolver").has_value()) { + OPM_THROW(std::invalid_argument, + "In general_system_cpr configuration with " + "preconditioner.coarse_solver.reservoir_solver.preconditioner.add_wells = true, " + "the '...reservoir_solver.preconditioner.coarsesolver' sub-tree is required: " + "it configures the solver for the CPRW pressure system."); + } +} + void validateSystemCPRTree(const PropertyTree& prm) { - if (prm.get("preconditioner.type", std::string{}) != "system_cpr") { + const auto type = prm.get("preconditioner.type", std::string{}); + if (type == "general_system_cpr") { + validateGeneralSystemCPRTree(prm); + return; + } + if (type != "system_cpr") { return; } for (const char* sub : {"reservoir_solver", "reservoir_smoother", "well_solver"}) { diff --git a/opm/simulators/linalg/setupPropertyTree.hpp b/opm/simulators/linalg/setupPropertyTree.hpp index 3d5da679706..3a15c72ce98 100644 --- a/opm/simulators/linalg/setupPropertyTree.hpp +++ b/opm/simulators/linalg/setupPropertyTree.hpp @@ -38,6 +38,8 @@ PropertyTree setupCPRW(const std::string& conf, const FlowLinearSolverParameters PropertyTree setupCPR(const std::string& conf, const FlowLinearSolverParameters& p); PropertyTree setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p); +PropertyTree setupGeneralSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p); + // Validates that a PropertyTree with preconditioner.type == "system_cpr" contains // all three required sub-trees. Throws std::invalid_argument if any are missing. // No-op when the preconditioner type is not system_cpr. From e767cf34c0a17c2a4e0a7e31e78641cefcaaf191 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 14:10:03 +0200 Subject: [PATCH 15/22] Assert the general preconditioner reproduces the fixed one test_GeneralSystemPreconditioner applies both to the same right-hand side on a system with a standard and a multisegment well. For well_transfer classic and no_prolongation, and for the reservoir-only pressure stage, it compares with BOOST_CHECK_EQUAL rather than a tolerance: TwoLevelMethodCpr reduces the defect with applyscaleadd(-1, lhs, rhs), and SystemMatrix::usmv accumulates A,C into the reservoir half and B,D into the well half -- the same operations in the same order as the fixed version's four mmv calls, so the two agree exactly. well_transfer = full is the exception and gets a tolerance. There the well half is a sum of three non-zero terms and the two associate it differently, ((0 + dw) + c1) + c2 against dw + ((0 + c1) + c2), because Dune adds the coarse correction to the update and then adds the smoother's separately accumulated correction. The reservoir half has two terms and stays exact; the other transfer modes stay exact because the coarse well correction is zero there. Dropping either well solve changes the result, which is what makes the equality meaningful. Measured on decks, general against fixed, linear iterations: system_cprw system_cpr SPE1CASE1 440 = 440 571 = 571 SPE9_CP 468 = 468 682 = 682 Norne 2248 = 2248 4303 vs 4253 Norne msw 2558 = 2558 4839 vs 4972 with all 2831 physical summary vectors identical on Norne for the cprw pairs (TCPU is the only one that differs) and the SPE files identical byte for byte. The system_cpr column differs on Norne alone, the only deck here that opens and closes wells: that path has no coarse space, so its reservoir CPR solve is rebuilt rather than refreshed. Refreshing in place reproduces 4303 and 4839. Co-Authored-By: Claude Opus 5 (cherry picked from commit 8e873e0936522790731680b7db21d9fe5a2509d6) --- CMakeLists_files.cmake | 2 + tests/test_GeneralSystemPreconditioner.cpp | 396 +++++++++++++++++++++ 2 files changed, 398 insertions(+) create mode 100644 tests/test_GeneralSystemPreconditioner.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 12b669ebe70..ca0d2145e5c 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -516,6 +516,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_tpsa_localresidual.cpp tests/test_tpsa_primaryvariables.cpp tests/test_vfpproperties.cpp + tests/test_GeneralSystemPreconditioner.cpp tests/test_SystemCprwPressureStage.cpp tests/test_WellMatrixMerger.cpp tests/test_WaterSatfuncConsistencyChecks.cpp @@ -1130,6 +1131,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/system/MultiComm.hpp opm/simulators/linalg/system/SystemCprwPressureStage.hpp opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp + opm/simulators/linalg/system/SystemPreconditionerParts.hpp opm/simulators/linalg/system/SystemPressureBhpTransferPolicy.hpp opm/simulators/linalg/system/SystemPreconditionerParts.hpp opm/simulators/linalg/system/SystemPreconditioner.hpp diff --git a/tests/test_GeneralSystemPreconditioner.cpp b/tests/test_GeneralSystemPreconditioner.cpp new file mode 100644 index 00000000000..baf8b461143 --- /dev/null +++ b/tests/test_GeneralSystemPreconditioner.cpp @@ -0,0 +1,396 @@ +/* + Copyright Equinor ASA 2026 + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include +#define BOOST_TEST_MODULE OPM_test_GeneralSystemPreconditioner +#include + +#include +#include +#include + +#include + +#include +#include +#include + +namespace { + +using Scalar = double; +using RRMatrix = Opm::RRMatrix; +using RWMatrix = Opm::RWMatrix; +using WRMatrix = Opm::WRMatrix; +using WWMatrix = Opm::WWMatrix; +using SystemVector = Opm::SystemVector; + +using Fixed = Opm::SystemPreconditioner>; +using General = Opm::GeneralSystemPreconditioner>; + +constexpr int numRes = Opm::numResDofs; +constexpr int numWell = Opm::numWellDofs; +constexpr int pressureIndex = 0; + +// 4 reservoir cells and 2 wells: well 0 a standard well (1 block row, +// perforating cells 0 and 1), well 1 a multisegment well (3 block rows), so +// the per-well aggregation is exercised. +constexpr std::size_t numCells = 4; +constexpr std::size_t numWellBlocks = 4; + +struct BlockSpec +{ + std::size_t column; + Scalar base; +}; + +using MatrixPattern = std::vector>; + +// Diagonally dominant so that every sub-solve is well posed, and +// non-symmetric so a transposed application cannot pass by accident. +template +Block makeBlock(const Scalar base, const bool diagonal) +{ + Block block; + for (int row = 0; row < Block::rows; ++row) { + for (int col = 0; col < Block::cols; ++col) { + block[row][col] = 0.05 * (base + 3.0 * row + 7.0 * col); + if (diagonal && row == col) { + block[row][col] += 40.0 + base; + } + } + } + return block; +} + +template +Matrix buildMatrix(const std::size_t rows, const std::size_t cols, + const MatrixPattern& pattern, const bool diagonal) +{ + std::size_t nonzeroes = 0; + for (const auto& row : pattern) { + nonzeroes += row.size(); + } + Matrix matrix(rows, cols, nonzeroes, Matrix::row_wise); + for (auto row = matrix.createbegin(); row != matrix.createend(); ++row) { + for (const auto& e : pattern[row.index()]) { + row.insert(e.column); + } + } + for (std::size_t row = 0; row < rows; ++row) { + for (const auto& e : pattern[row]) { + matrix[row][e.column] = + makeBlock(e.base, diagonal && e.column == row); + } + } + return matrix; +} + +struct Fixture +{ + RRMatrix A; + RWMatrix C; + WRMatrix B; + WWMatrix D; + Opm::WellDofLayout layout; + Opm::SystemMatrix S; + SystemVector weights; + + Fixture() + { + A = buildMatrix(numCells, numCells, + {{{0, 1.0}, {1, 2.0}}, + {{0, 3.0}, {1, 4.0}, {2, 5.0}}, + {{1, 6.0}, {2, 7.0}, {3, 8.0}}, + {{2, 9.0}, {3, 10.0}}}, + /*diagonal=*/true); + + C = buildMatrix(numCells, numWellBlocks, + {{{0, 11.0}}, + {{0, 12.0}}, + {{1, 13.0}, {2, 14.0}}, + {{1, 15.0}, {3, 16.0}}}, + false); + + B = buildMatrix(numWellBlocks, numCells, + {{{0, 17.0}, {1, 18.0}}, + {{2, 19.0}}, + {{2, 20.0}}, + {{3, 21.0}}}, + false); + + D = buildMatrix(numWellBlocks, numWellBlocks, + {{{0, 22.0}}, + {{1, 23.0}, {2, 24.0}}, + {{1, 25.0}, {2, 26.0}, {3, 27.0}}, + {{2, 28.0}, {3, 29.0}}}, + /*diagonal=*/true); + + layout.wellBlockOffsets = {0, 1, 4}; + layout.pressureDofIndex = numWell - 1; + + S.A = &A; + S.C = &C; + S.B = &B; + S.D = &D; + S.wellLayout = &layout; + + weights[Dune::Indices::_0].resize(numCells); + for (std::size_t c = 0; c < numCells; ++c) { + for (int i = 0; i < numRes; ++i) { + weights[Dune::Indices::_0][c][i] = 0.5 + 0.1 * c + 0.3 * i; + } + } + weights[Dune::Indices::_1].resize(numWellBlocks); + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + for (int i = 0; i < numWell; ++i) { + weights[Dune::Indices::_1][wb][i] = 0.2 + 0.7 * wb - 0.15 * i; + } + } + } + + std::function weightsCalculator() const + { + const auto w = weights; + return [w]() { return w; }; + } +}; + +// One "apply once" sub-solve: run the configured preconditioner a single time +// without a Krylov method around it. +void putOnceSolver(Opm::PropertyTree& prm, const std::string& at, const std::string& type) +{ + prm.put(at + ".maxiter", 1); + prm.put(at + ".tol", 1e-2); + prm.put(at + ".verbosity", 0); + prm.put(at + ".solver", std::string{"preconditioner2inverseoperator"}); + prm.put(at + ".preconditioner.type", type); + prm.put(at + ".preconditioner.relaxation", 1.0); +} + +void putReservoirSolver(Opm::PropertyTree& prm, const std::string& at, const bool addWells) +{ + putOnceSolver(prm, at, "cpr"); + prm.put(at + ".preconditioner.use_well_weights", std::string{"false"}); + prm.put(at + ".preconditioner.add_wells", addWells ? std::string{"true"} : std::string{"false"}); + prm.put(at + ".preconditioner.weight_type", std::string{"trueimpes"}); + prm.put(at + ".preconditioner.pre_smooth", 0); + prm.put(at + ".preconditioner.post_smooth", 0); + prm.put(at + ".preconditioner.finesmoother.type", std::string{"jac"}); + prm.put(at + ".preconditioner.finesmoother.relaxation", 1.0); + prm.put(at + ".preconditioner.verbosity", 0); + // The coarse pressure system is tiny here, so a single ILU0 apply is both + // adequate and free of any optional dependency. + prm.put(at + ".preconditioner.coarsesolver.maxiter", 1); + prm.put(at + ".preconditioner.coarsesolver.tol", 1e-1); + prm.put(at + ".preconditioner.coarsesolver.solver", std::string{"preconditioner2inverseoperator"}); + prm.put(at + ".preconditioner.coarsesolver.verbosity", 0); + prm.put(at + ".preconditioner.coarsesolver.preconditioner.type", std::string{"ilu0"}); + prm.put(at + ".preconditioner.coarsesolver.preconditioner.relaxation", 1.0); +} + +void putWellOptions(Opm::PropertyTree& prm, const std::string& transfer) +{ + prm.put("well_transfer", transfer); + prm.put("well_coarse_diagonal", std::string{"contract_d"}); + prm.put("well_identity_on_pressure_control", std::string{"false"}); + prm.put("verbosity", 0); +} + +Opm::PropertyTree fixedTree(const bool addWells, const std::string& transfer) +{ + Opm::PropertyTree prm; + prm.put("type", std::string{"system_cpr"}); + putWellOptions(prm, transfer); + putReservoirSolver(prm, "reservoir_solver", addWells); + putOnceSolver(prm, "reservoir_smoother", "ilu0"); + putOnceSolver(prm, "well_solver", "ilu0"); + return prm; +} + +Opm::PropertyTree generalTree(const bool addWells, const std::string& transfer, + const bool coarseWellSolve = true, + const bool smootherWellSolve = true) +{ + Opm::PropertyTree prm; + prm.put("type", std::string{"general_system_cpr"}); + putWellOptions(prm, transfer); + putReservoirSolver(prm, "coarse_solver.reservoir_solver", addWells); + if (coarseWellSolve) { + putOnceSolver(prm, "coarse_solver.well_solver", "ilu0"); + } + putOnceSolver(prm, "smoother.reservoir_smoother", "ilu0"); + if (smootherWellSolve) { + putOnceSolver(prm, "smoother.well_solver", "ilu0"); + } + return prm; +} + +SystemVector makeRhs() +{ + SystemVector d; + d[Dune::Indices::_0].resize(numCells); + d[Dune::Indices::_1].resize(numWellBlocks); + for (std::size_t c = 0; c < numCells; ++c) { + for (int i = 0; i < numRes; ++i) { + d[Dune::Indices::_0][c][i] = 1.0 + 0.37 * c - 0.11 * i; + } + } + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + for (int i = 0; i < numWell; ++i) { + d[Dune::Indices::_1][wb][i] = 2.0 - 0.23 * wb + 0.41 * i; + } + } + return d; +} + +SystemVector zeroLike(const SystemVector& d) +{ + SystemVector v; + v[Dune::Indices::_0].resize(d[Dune::Indices::_0].size()); + v[Dune::Indices::_1].resize(d[Dune::Indices::_1].size()); + v[Dune::Indices::_0] = 0.0; + v[Dune::Indices::_1] = 0.0; + return v; +} + +// Bit for bit, not merely close: the point of the decomposition is that the +// composed sweep performs the same floating point operations in the same order. +void checkIdentical(const SystemVector& a, const SystemVector& b) +{ + BOOST_REQUIRE_EQUAL(a[Dune::Indices::_0].size(), b[Dune::Indices::_0].size()); + BOOST_REQUIRE_EQUAL(a[Dune::Indices::_1].size(), b[Dune::Indices::_1].size()); + for (std::size_t c = 0; c < a[Dune::Indices::_0].size(); ++c) { + for (int i = 0; i < numRes; ++i) { + BOOST_CHECK_EQUAL(a[Dune::Indices::_0][c][i], b[Dune::Indices::_0][c][i]); + } + } + for (std::size_t wb = 0; wb < a[Dune::Indices::_1].size(); ++wb) { + for (int i = 0; i < numWell; ++i) { + BOOST_CHECK_EQUAL(a[Dune::Indices::_1][wb][i], b[Dune::Indices::_1][wb][i]); + } + } +} + +bool differs(const SystemVector& a, const SystemVector& b) +{ + for (std::size_t c = 0; c < a[Dune::Indices::_0].size(); ++c) { + for (int i = 0; i < numRes; ++i) { + if (a[Dune::Indices::_0][c][i] != b[Dune::Indices::_0][c][i]) { + return true; + } + } + } + for (std::size_t wb = 0; wb < a[Dune::Indices::_1].size(); ++wb) { + for (int i = 0; i < numWell; ++i) { + if (a[Dune::Indices::_1][wb][i] != b[Dune::Indices::_1][wb][i]) { + return true; + } + } + } + return false; +} + +// Both preconditioners applied once to the same right-hand side. +std::pair +applyBoth(const Fixture& f, const Opm::PropertyTree& fixedPrm, const Opm::PropertyTree& generalPrm) +{ + const auto d = makeRhs(); + + Fixed fixed(f.S, f.weightsCalculator(), pressureIndex, fixedPrm); + General general(f.S, f.weightsCalculator(), pressureIndex, generalPrm); + + auto vFixed = zeroLike(d); + auto vGeneral = zeroLike(d); + fixed.apply(vFixed, d); + general.apply(vGeneral, d); + return {vFixed, vGeneral}; +} + +} // anonymous namespace + +BOOST_AUTO_TEST_CASE(MatchesFixedThreeStageWithCprwPressureStage) +{ + const Fixture f; + for (const auto* transfer : {"classic", "no_prolongation"}) { + const auto [vFixed, vGeneral] = + applyBoth(f, fixedTree(true, transfer), generalTree(true, transfer)); + BOOST_TEST_CONTEXT("well_transfer = " << transfer) { + checkIdentical(vFixed, vGeneral); + } + } +} + +// well_transfer = full is the one mode where the coarse correction reaches the +// well unknowns, so the well half of the result is a sum of three terms rather +// than two: the coarse correction and the two well solves. The fixed sequence +// accumulates them as ((0 + dw) + c1) + c2, while the two-level method adds the +// coarse correction to the update and then adds the smoother's own accumulated +// correction, dw + ((0 + c1) + c2). Same terms, same order, different +// parenthesisation -- so the two agree to a rounding step and no closer. The +// reservoir half has only two terms and stays exact, and so do the other +// transfer modes, where dw is exactly zero. +BOOST_AUTO_TEST_CASE(MatchesFixedThreeStageToRoundoffWithProlongedWellPressure) +{ + const Fixture f; + const auto [vFixed, vGeneral] = applyBoth(f, fixedTree(true, "full"), generalTree(true, "full")); + + for (std::size_t c = 0; c < numCells; ++c) { + for (int i = 0; i < numRes; ++i) { + BOOST_CHECK_EQUAL(vFixed[Dune::Indices::_0][c][i], vGeneral[Dune::Indices::_0][c][i]); + } + } + for (std::size_t wb = 0; wb < numWellBlocks; ++wb) { + for (int i = 0; i < numWell; ++i) { + BOOST_CHECK_CLOSE(vFixed[Dune::Indices::_1][wb][i], + vGeneral[Dune::Indices::_1][wb][i], 1e-12); + } + } +} + +BOOST_AUTO_TEST_CASE(MatchesFixedThreeStageWithReservoirOnlyPressureStage) +{ + const Fixture f; + const auto [vFixed, vGeneral] = + applyBoth(f, fixedTree(false, "classic"), generalTree(false, "classic")); + checkIdentical(vFixed, vGeneral); +} + +BOOST_AUTO_TEST_CASE(DroppingAWellSolveChangesTheResult) +{ + // Guards the equivalence above: if the composed parts were somehow being + // ignored, leaving one out would not change anything. + const Fixture f; + const auto both = applyBoth(f, fixedTree(true, "classic"), generalTree(true, "classic")); + const auto noCoarseWell = + applyBoth(f, fixedTree(true, "classic"), generalTree(true, "classic", false, true)); + const auto noSmootherWell = + applyBoth(f, fixedTree(true, "classic"), generalTree(true, "classic", true, false)); + + BOOST_CHECK(differs(both.second, noCoarseWell.second)); + BOOST_CHECK(differs(both.second, noSmootherWell.second)); +} + +BOOST_AUTO_TEST_CASE(RejectsAnEmptyComposition) +{ + const Fixture f; + Opm::PropertyTree prm; + prm.put("type", std::string{"general_system_cpr"}); + putWellOptions(prm, "classic"); + BOOST_CHECK_THROW(General(f.S, f.weightsCalculator(), pressureIndex, prm), + std::invalid_argument); +} From 3acf44f379c49b10aab34949e1a4e08e73c58059 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 15:21:28 +0200 Subject: [PATCH 16/22] Skip the residual updates a sweep never reads again A step subtracted its correction from both halves of the residual whether or not anything read them afterwards. For a reservoir step that is a full A*corr over every cell, and in the default composition -- coarse, well, reservoir smoother, well -- nothing reads the reservoir half after the smoother. The hand-written three-stage preconditioner leaves the same updates out by hand. Which halves are still read is decided once, walking the step list backwards: a step needs the reservoir half only if a later step is a reservoir step, the well half only if a later step is a well step, and the last step needs neither. Skipping an update nothing reads cannot change a result, and the equality against the fixed preconditioner still holds. On full Norne it takes the general variant's overhead from 3.8% to about 1.7%, which is within the run-to-run spread. Co-Authored-By: Claude Opus 5 (cherry picked from commit e969719fd6309425d0f5ab4a828b9f23cd1a6b5c) --- .../system/SystemPreconditionerParts.hpp | 117 ++++++++++++++++-- 1 file changed, 109 insertions(+), 8 deletions(-) diff --git a/opm/simulators/linalg/system/SystemPreconditionerParts.hpp b/opm/simulators/linalg/system/SystemPreconditionerParts.hpp index c3325b684a3..908aca65f46 100644 --- a/opm/simulators/linalg/system/SystemPreconditionerParts.hpp +++ b/opm/simulators/linalg/system/SystemPreconditionerParts.hpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -61,6 +62,43 @@ namespace Opm // block Gauss-Seidel sweep, not a block Jacobi one. // -------------------------------------------------------------------------- +// What a well change did to the structure the solver was built against. The +// distinction is between a change that stays inside what the initial build +// anticipated and one that does not: only the latter forces anything sized by +// the wells to be created again. +enum class WellStructureChange +{ + // Same wells, same patterns, new numbers. + Values, + // Same dimensions, different sparsity -- a connection opened or closed + // inside a well that was already there. + Pattern, + // A well or a segment appeared or vanished, so D and the coarse system + // change size. Nothing sized by the wells survives this. + Dimension, +}; + +// What to do about it. +enum class WellStructureUpdate +{ + // Build every part again from the property tree. Blunt, and the default: + // it cannot be wrong, only wasteful. + Rebuild, + // Create only the parts the wells size -- the well solves and the coarse + // system -- and refresh the rest in place. This is what the fixed + // three-stage SystemPreconditioner does, so it reproduces it exactly. + Refresh, +}; + +inline WellStructureUpdate wellStructureUpdateFromString(const std::string& name) +{ + if (name == "rebuild") { return WellStructureUpdate::Rebuild; } + if (name == "refresh") { return WellStructureUpdate::Refresh; } + OPM_THROW(std::invalid_argument, + "Unknown well_structure_update '" + name + + "'. Valid values are 'rebuild' and 'refresh'."); +} + template struct SystemSweepState { @@ -78,6 +116,17 @@ class SystemSweepStep virtual void apply(SystemSweepState& state) = 0; + // Which halves of the residual a later step still reads. A step only + // subtracts its correction from those: updating a half nothing reads again + // cannot change any result, and the matvec that does it is not free -- for + // a reservoir step it is a full A*corr over every cell. The hand-written + // three-stage preconditioner leaves the same updates out by hand. + void setResidualNeeded(const bool res0, const bool res1) + { + needRes0_ = res0; + needRes1_ = res1; + } + // Refresh against changed matrix values, pattern unchanged. virtual void update() = 0; @@ -89,6 +138,14 @@ class SystemSweepStep { update(); } + + // Which block a step corrects, so a sweep can work out what its residual + // updates are still good for. + virtual bool correctsReservoir() const = 0; + +protected: + bool needRes0_ = true; + bool needRes1_ = true; }; // -------------------------------------------------------------------------- @@ -140,8 +197,17 @@ class SystemWellStep : public SystemSweepStep solver_->apply(corr_, rhs_, result); v[_1] += corr_; - S_.C->mmv(corr_, res[_0]); - S_.D->mmv(corr_, res[_1]); + if (this->needRes0_) { + S_.C->mmv(corr_, res[_0]); + } + if (this->needRes1_) { + S_.D->mmv(corr_, res[_1]); + } + } + + bool correctsReservoir() const override + { + return false; } private: @@ -222,8 +288,17 @@ class SystemReservoirStep : public SystemSweepStep solver_->apply(corr_, rhs_, result); v[_0] += corr_; - S_.A->mmv(corr_, res[_0]); - S_.B->mmv(corr_, res[_1]); + if (this->needRes0_) { + S_.A->mmv(corr_, res[_0]); + } + if (this->needRes1_) { + S_.B->mmv(corr_, res[_1]); + } + } + + bool correctsReservoir() const override + { + return true; } private: @@ -312,16 +387,29 @@ class SystemCprwStep : public SystemSweepStep stage_->apply(rhs_, res[_1], weights_, corrRes_, corrWell_); v[_0] += corrRes_; - S_.A->mmv(corrRes_, res[_0]); - S_.B->mmv(corrRes_, res[_1]); + if (this->needRes0_) { + S_.A->mmv(corrRes_, res[_0]); + } + if (this->needRes1_) { + S_.B->mmv(corrRes_, res[_1]); + } if (stage_->prolongatesWellPressure()) { v[_1] += corrWell_; - S_.C->mmv(corrWell_, res[_0]); - S_.D->mmv(corrWell_, res[_1]); + if (this->needRes0_) { + S_.C->mmv(corrWell_, res[_0]); + } + if (this->needRes1_) { + S_.D->mmv(corrWell_, res[_1]); + } } } + bool correctsReservoir() const override + { + return true; + } + private: const SystemMatrix& S_; const ResComm* comm_ = nullptr; @@ -356,6 +444,19 @@ class SystemSweepPreconditioner : steps_(std::move(steps)) , parallel_(parallel) { + // Walking backwards, a step needs to update the half of the residual + // that some later step reads: reservoir steps read the reservoir half, + // well steps the well half. Nothing reads either after the last step. + bool laterReservoir = false; + bool laterWell = false; + for (auto step = steps_.rbegin(); step != steps_.rend(); ++step) { + (*step)->setResidualNeeded(laterReservoir, laterWell); + if ((*step)->correctsReservoir()) { + laterReservoir = true; + } else { + laterWell = true; + } + } } void pre(SystemVector&, SystemVector&) override {} From 22bd29a000925a4b6b08cb4e69297b9a02c89a41 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 15:21:28 +0200 Subject: [PATCH 17/22] Make the response to a changed well structure configurable Two things were conflated. What the change was: Values same wells, same patterns, new numbers Pattern same dimensions, different sparsity -- a connection opened or closed inside a well that was already there Dimension a well or segment appeared or vanished, so D and the coarse system change size and what to do about it, now preconditioner.well_structure_update: rebuild build every part again from the property tree (default) refresh create only what the wells size -- the well solves and the coarse system -- and refresh the rest in place These are different preconditioners, not two routes to one: a CPR reservoir solve keeps its hierarchy when refreshed and re-aggregates it when rebuilt. On full Norne, general_system_cpr gives 4253 linear iterations rebuilding and 4303 refreshing, and 4303 is what the fixed three-stage system_cpr gives. With refresh both variants reproduce it exactly -- all 2831 physical summary vectors identical, TCPU aside -- at 143s against 142.5s for cprw and 172s against 175.5s for cpr. WellMatrixStructure::hasSameDimensions separates Pattern from Dimension and ISTLSolverSystem passes the verdict down. Both still force the coarse system to be built again, since a changed sparsity invalidates its pattern as surely as a changed dimension does; the distinction is carried so that a cheaper Pattern path can be added without re-deriving it. Co-Authored-By: Claude Opus 5 (cherry picked from commit ac602cce5c7eee63e8517346904fa12236a5146a) --- opm/simulators/linalg/setupPropertyTree.cpp | 10 +++++ .../system/GeneralSystemPreconditioner.hpp | 42 +++++++++++++++---- .../linalg/system/ISTLSolverSystem.hpp | 20 ++++++--- .../linalg/system/WellMatrixMerger.hpp | 10 +++++ 4 files changed, 68 insertions(+), 14 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index e05ddea7bc0..133673516be 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -692,6 +692,16 @@ setupGeneralSystemCPR(const std::string& conf, const FlowLinearSolverParameters& prm.put("preconditioner.type", "general_system_cpr"s); setupSystemCPRWellOptions(prm); + // What to do when a well change introduces something the initial build of + // the structure did not have: + // rebuild - build every part again from this tree + // refresh - build only the parts the wells size (the well solves and the + // coarse system) and refresh the rest in place, which is what + // the fixed three-stage system_cpr does + // These are different preconditioners, not two routes to the same one: a + // CPR reservoir solve keeps its hierarchy when refreshed and re-aggregates + // it when rebuilt. + prm.put("preconditioner.well_structure_update", "rebuild"s); // Each stage gets the well solver the fixed layout shares between its // stages, so the sequence is coarse -> well -> smoother -> well. diff --git a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp index 6807829e66d..247842c80ad 100644 --- a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp +++ b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp @@ -165,16 +165,36 @@ class GeneralSystemPreconditioner } } - void updateForChangedWellStructure() + // A well change that the initial structure did not anticipate. What to do + // about it is preconditioner.well_structure_update; see WellStructureUpdate. + // + // This is not a cosmetic choice. Refreshing a CPR reservoir solve keeps its + // existing hierarchy while rebuilding re-aggregates it, and the two are + // different preconditioners: on Norne it moves general_system_cpr between + // 4303 iterations (refresh, which is what the fixed three-stage version + // does) and 4253 (rebuild). + void updateForChangedWellStructure(const WellStructureChange change + = WellStructureChange::Dimension) { - // A changed well structure changes D's dimension and the dimension of - // the coarse system, so everything is built again from the property - // tree. That is deliberately blunt: refreshing the reservoir solves in - // place instead would keep a CPR hierarchy rather than re-aggregate it, - // and the two are different preconditioners -- on Norne that choice - // alone moves general_system_cpr between 4303 and 4253 linear - // iterations. Rebuilding is the safe default until it is measured. - build(); + if (update_ == WellStructureUpdate::Rebuild) { + build(); + return; + } + + // Refresh: the well solves are created again because D changed size, + // the reservoir solves only refreshed. + sweep_->rebuildForChangedWellStructure(); + + // The coarse system carries one unknown per well and one row per cell, + // so both a changed dimension and a changed sparsity invalidate it -- + // only Values leaves it alone, and that never reaches this function. + if (coarseResPrm_ && change != WellStructureChange::Values) { + weights_ = weightsCalculator_(); + buildCoarse(); + twoLevel_ = std::make_unique(*fineOp_, sweep_, *transfer_, + *coarsePolicy_, + /*preSteps=*/0, /*postSteps=*/1); + } } void apply(SystemVector& v, const SystemVector& d) override @@ -213,6 +233,7 @@ class GeneralSystemPreconditioner std::unique_ptr twoLevel_; // Kept so the coarse level can be built again when the well count changes. std::optional coarseResPrm_; + WellStructureUpdate update_ = WellStructureUpdate::Rebuild; void build() { @@ -223,6 +244,9 @@ class GeneralSystemPreconditioner sweep_.reset(); coarseResPrm_.reset(); + update_ = wellStructureUpdateFromString( + prm_.get("well_structure_update", std::string{"rebuild"})); + // The weights arrive from the outer layer for the whole system; the // reservoir-only sub-solvers want just their own part of them. std::function()> resWeightCalc; diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 1ec4b40d55e..69842a2750d 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -234,9 +234,18 @@ class ISTLSolverSystem : public ISTLSolver sysMatrix_.C = &mergedC_; sysMatrix_.D = &mergedD_; sysMatrix_.wellLayout = &wellLayout_; - cachedWellStructure_ = merger.buildStructure(); - refreshSystemSolverForChangedWellStructure(prm); + const auto newStructure = merger.buildStructure(); + // A connection opening or closing inside an existing well keeps + // every dimension; a well or segment appearing or vanishing does + // not, and only the latter introduces unknowns the initial build + // never saw. + const auto change = (sysInitialized_ && newStructure.hasSameDimensions(cachedWellStructure_)) + ? WellStructureChange::Pattern + : WellStructureChange::Dimension; + cachedWellStructure_ = newStructure; + + refreshSystemSolverForChangedWellStructure(prm, change); sysInitialized_ = true; } else { OPM_TIMEBLOCK(flexibleSolverUpdate); @@ -387,7 +396,8 @@ class ISTLSolverSystem : public ISTLSolver return weights; } - void refreshSystemSolverForChangedWellStructure(const Opm::PropertyTree& prm) + void refreshSystemSolverForChangedWellStructure(const Opm::PropertyTree& prm, + const WellStructureChange change) { if (!sysInitialized_ || !sysPrecond_) { createSystemSolver(prm); @@ -399,7 +409,7 @@ class ISTLSolverSystem : public ISTLSolver if (auto* precond = dynamic_cast(sysPrecond_)) { precond->updateForChangedWellStructure(); } else if (auto* general = dynamic_cast(sysPrecond_)) { - general->updateForChangedWellStructure(); + general->updateForChangedWellStructure(change); } else { // Rebuild the parallel solver if the parallel preconditioner cannot be updated in-place. createSystemSolver(prm); @@ -411,7 +421,7 @@ class ISTLSolverSystem : public ISTLSolver if (auto* precond = dynamic_cast(sysPrecond_)) { precond->updateForChangedWellStructure(); } else if (auto* general = dynamic_cast(sysPrecond_)) { - general->updateForChangedWellStructure(); + general->updateForChangedWellStructure(change); } else { // Rebuild the solver if the sequential preconditioner cannot be updated in-place createSystemSolver(prm); diff --git a/opm/simulators/linalg/system/WellMatrixMerger.hpp b/opm/simulators/linalg/system/WellMatrixMerger.hpp index 77b270cb9ee..55039fb5eae 100644 --- a/opm/simulators/linalg/system/WellMatrixMerger.hpp +++ b/opm/simulators/linalg/system/WellMatrixMerger.hpp @@ -82,6 +82,16 @@ struct WellMatrixStructure { return !(*this == other); } + + // Whether everything sized by the wells still has the same size. A change + // that keeps these can be absorbed by rebuilding patterns; a change that + // does not introduces unknowns the initial build never saw. + bool hasSameDimensions(const WellMatrixStructure& other) const + { + return numResDofs == other.numResDofs + && totalWellBlocks == other.totalWellBlocks + && wellCells.size() == other.wellCells.size(); + } }; template From 81d30720abdf4a92d85e749580ff1b1dd3961af8 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 16:05:10 +0200 Subject: [PATCH 18/22] PropertyTree: read and write arrays of sub trees A JSON array parses to a node whose children all carry an empty key, so they cannot be reached by name, and get_child_items_as_vector reads scalars rather than sub trees. get_child_list and put_child_list handle the object-array case, which a configuration describing an ordered list of steps needs. Co-Authored-By: Claude Opus 5 (cherry picked from commit 87bcc1e009cfef3ad5020b9a5458864cfb56d4fa) --- opm/simulators/linalg/PropertyTree.cpp | 30 ++++++++++++++++++++++++++ opm/simulators/linalg/PropertyTree.hpp | 25 +++++++++++++++++++++ 2 files changed, 55 insertions(+) diff --git a/opm/simulators/linalg/PropertyTree.cpp b/opm/simulators/linalg/PropertyTree.cpp index 4c1ab1159dd..7a317f690ba 100644 --- a/opm/simulators/linalg/PropertyTree.cpp +++ b/opm/simulators/linalg/PropertyTree.cpp @@ -114,6 +114,36 @@ std::vector PropertyTree::get_child_keys() const return keys; } +void PropertyTree::put_child_list(const std::string& key, + const std::vector& items) +{ + // An array is a node whose children all have an empty key. + boost::property_tree::ptree array; + for (const auto& item : items) { + array.push_back(std::make_pair(std::string{}, *item.tree_)); + } + tree_->put_child(key, array); +} + +std::optional> +PropertyTree::get_child_list(const std::string& child) const +{ + auto subTree = this->tree_->get_child_optional(child); + if (! subTree) { + return std::nullopt; + } + + // A JSON array parses to a node whose children all have an empty key, so + // they cannot be reached by name; walk them in order instead. + std::vector items; + items.reserve(subTree->size()); + for (const auto& item : *subTree) { + items.emplace_back(PropertyTree(item.second)); + } + + return items; +} + template std::optional> PropertyTree::get_child_items_as_vector(const std::string& child) const diff --git a/opm/simulators/linalg/PropertyTree.hpp b/opm/simulators/linalg/PropertyTree.hpp index 9fb3fdbfea3..a6b9673ee99 100644 --- a/opm/simulators/linalg/PropertyTree.hpp +++ b/opm/simulators/linalg/PropertyTree.hpp @@ -137,6 +137,31 @@ class PropertyTree /// \return Vector of strings containing the names of all immediate children std::vector get_child_keys() const; + /// Retrieve a node's children as sub trees, in order. + /// + /// For a JSON array of objects. Such an array parses to a node whose + /// children all have an empty key, so they cannot be reached by name and + /// get_child_items_as_vector() cannot read them either -- that one reads + /// scalars. + /// + /// \param[in] child Property key. Expected to be in hierarchical + /// notation for subtrees--i.e., using periods ('.') to separate + /// hierarchy levels. + /// + /// \return The child sub trees in the order they appear. Nullopt if no + /// node named by \p child exists. + std::optional> + get_child_list(const std::string& child) const; + + /// Store a sequence of sub trees as a JSON array. + /// + /// \param[in] key Property key. Expected to be in hierarchical + /// notation for subtrees--i.e., using periods ('.') to separate + /// hierarchy levels. + /// + /// \param[in] items Sub trees, stored in the order given. + void put_child_list(const std::string& key, const std::vector& items); + /// Retrieve node items as linearised vector. /// /// Assumes that the node's child is an array type of homongeneous From 87863646cc463810d7ddd491cead47ea9dafa9d3 Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 6 Aug 2026 16:05:10 +0200 Subject: [PATCH 19/22] Give general_system_cpr a JSON layout shaped like cpr's The composition was fixed in four named slots, so the one thing "general" ought to mean could not be said in JSON. The smoother is now an ordered list: "preconditioner": { "type": "general_system_cpr", "weight_type": "trueimpes", "pre_smooth": 0, "post_smooth": 1, "well_structure_update": "rebuild", "coarsesolver": { "type": "cprw_pressure", "well_weight_type": ..., "well_transfer": ..., "well_coarse_diagonal": ..., "well_identity_on_pressure_control": ..., "maxiter": 1, "tol": 0.1, "solver": "loopsolver", "preconditioner": { "type": "amg", ... } }, "finesmoother": { "steps": [ { "block": "well", ... }, { "block": "reservoir", ... }, { "block": "well", ... } ] } } so steps can be reordered, repeated or left out. The names follow cpr -- coarsesolver, finesmoother, weight_type, pre_smooth, post_smooth in the same places -- and the sub-trees are ordinary solver specs, as they were. Three things this removes. add_wells nested inside a reservoir solver's preconditioner was deciding whether a two-level method got built at all; that is now coarsesolver.type, and leaving coarsesolver out means no coarse space. The coarse pressure solver was buried at coarse_solver.reservoir_solver.preconditioner.coarsesolver, a coarsesolver inside a reservoir_solver that was itself the coarse solver; it is now the coarsesolver node directly. And the keys describing the coarse space sat at the top of the preconditioner while being read by ISTLSolverSystem rather than the preconditioner; they now sit beside the coarse solver they belong to, and both readers look in one place. SPE1CASE1 440/571 and SPE9_CP 468/682, matching system_cprw and system_cpr as before. Co-Authored-By: Claude Opus 5 (cherry picked from commit e79a9e92acaff09f12f7366f3c534d78fde39f46) --- opm/simulators/linalg/ISTLSolver.hpp | 24 ++- opm/simulators/linalg/setupPropertyTree.cpp | 163 ++++++++++++------ .../system/GeneralSystemPreconditioner.hpp | 107 ++++++------ .../linalg/system/ISTLSolverSystem.hpp | 36 ++-- tests/test_GeneralSystemPreconditioner.cpp | 88 +++++++--- 5 files changed, 274 insertions(+), 144 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index 109d52b90cb..bcdf4ddad5e 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -596,18 +596,30 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, const Matrix& matrix, std::size_t pressIndex) const { - std::function weightsCalculator; - using namespace std::string_literals; auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s); // We use lower case as the internal canonical representation of solver names std::ranges::transform(preconditionerType, preconditionerType.begin(), ::tolower); - if (preconditionerType == "cpr" || preconditionerType == "cprt" - || preconditionerType == "cprw" || preconditionerType == "cprwt") { - const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt"; + if (preconditionerType != "cpr" && preconditionerType != "cprt" + && preconditionerType != "cprw" && preconditionerType != "cprwt") { + return {}; + } + const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt"; + return makeWeightsCalculator(prm.get("preconditioner.weight_type"s, "quasiimpes"s), + matrix, pressIndex, transpose); + } + + // The same, for a caller that knows the weighting directly rather than + // through a CPR preconditioner sub-tree. + std::function makeWeightsCalculator(const std::string& weightsType, + const Matrix& matrix, + std::size_t pressIndex, + const bool transpose = false) const + { + std::function weightsCalculator; + { const bool enableThreadParallel = this->parameters_[0].cpr_weights_thread_parallel_; - const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s); if (weightsType == "quasiimpes") { // weights will be created as default in the solver // assignment p = pressureIndex prevent compiler warning about diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 133673516be..707ad8a3d40 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -530,6 +530,12 @@ setupUMFPack([[maybe_unused]] const std::string& conf, const FlowLinearSolverPar namespace { +// Join a node key with a leaf name, allowing the node to be the tree root. +std::string at_(const std::string& at, const std::string& leaf) +{ + return at.empty() ? leaf : at + "." + leaf; +} + // The sub-solvers of the system preconditioner, each written at a caller // chosen key. Both the fixed three-stage layout and the general one use // these, so the trees they produce are identical and the two preconditioners @@ -539,14 +545,14 @@ void setupSystemReservoirSmoother(PropertyTree& prm, const FlowLinearSolverParameters& p) { using namespace std::string_literals; - prm.put(at + ".maxiter", 1); - prm.put(at + ".tol", p.linear_solver_reduction_); - prm.put(at + ".verbosity", 0); + prm.put(at_(at, "maxiter"), 1); + prm.put(at_(at, "tol"), p.linear_solver_reduction_); + prm.put(at_(at, "verbosity"), 0); // Apply the configured smoother once without loopsolver's extra defect // updates and convergence bookkeeping. - prm.put(at + ".solver", "preconditioner2inverseoperator"s); - prm.put(at + ".preconditioner.type", "paroverilu0"s); - prm.put(at + ".preconditioner.relaxation", 1.0); + prm.put(at_(at, "solver"), "preconditioner2inverseoperator"s); + prm.put(at_(at, "preconditioner.type"), "paroverilu0"s); + prm.put(at_(at, "preconditioner.relaxation"), 1.0); } void setupSystemReservoirSolver(PropertyTree& prm, @@ -555,31 +561,31 @@ void setupSystemReservoirSolver(PropertyTree& prm, const bool add_wells) { using namespace std::string_literals; - prm.put(at + ".maxiter", 1); - prm.put(at + ".tol", p.linear_solver_reduction_); - prm.put(at + ".verbosity", 0); + prm.put(at_(at, "maxiter"), 1); + prm.put(at_(at, "tol"), p.linear_solver_reduction_); + prm.put(at_(at, "verbosity"), 0); // Apply the configured reservoir preconditioner once without loopsolver's // extra defect updates and convergence bookkeeping. - prm.put(at + ".solver", "preconditioner2inverseoperator"s); - prm.put(at + ".preconditioner.type", "cpr"s); - prm.put(at + ".preconditioner.relaxation", 1.0); - prm.put(at + ".preconditioner.use_well_weights", "false"s); + prm.put(at_(at, "solver"), "preconditioner2inverseoperator"s); + prm.put(at_(at, "preconditioner.type"), "cpr"s); + prm.put(at_(at, "preconditioner.relaxation"), 1.0); + prm.put(at_(at, "preconditioner.use_well_weights"), "false"s); // add_wells promotes the pressure stage from reservoir-only CPR to CPRW // over the full (reservoir, well) system. - prm.put(at + ".preconditioner.add_wells", add_wells ? "true"s : "false"s); - prm.put(at + ".preconditioner.weight_type", "trueimpes"s); - prm.put(at + ".preconditioner.pre_smooth", 0); - prm.put(at + ".preconditioner.post_smooth", 0); + prm.put(at_(at, "preconditioner.add_wells"), add_wells ? "true"s : "false"s); + prm.put(at_(at, "preconditioner.weight_type"), "trueimpes"s); + prm.put(at_(at, "preconditioner.pre_smooth"), 0); + prm.put(at_(at, "preconditioner.post_smooth"), 0); // Set unused finesmoother to jac to avoid spending time setuping an ILU smoother that won't be used. - prm.put(at + ".preconditioner.finesmoother.type", "jac"s); - prm.put(at + ".preconditioner.finesmoother.relaxation", 1.0); - prm.put(at + ".preconditioner.verbosity", 0); - prm.put(at + ".preconditioner.coarsesolver.maxiter", 1); - prm.put(at + ".preconditioner.coarsesolver.tol", 1e-1); - prm.put(at + ".preconditioner.coarsesolver.solver", "loopsolver"s); - prm.put(at + ".preconditioner.coarsesolver.verbosity", 0); - prm.put(at + ".preconditioner.coarsesolver.preconditioner.type", "amg"s); - setupDuneAMG(prm, at + ".preconditioner.coarsesolver.preconditioner."); + prm.put(at_(at, "preconditioner.finesmoother.type"), "jac"s); + prm.put(at_(at, "preconditioner.finesmoother.relaxation"), 1.0); + prm.put(at_(at, "preconditioner.verbosity"), 0); + prm.put(at_(at, "preconditioner.coarsesolver.maxiter"), 1); + prm.put(at_(at, "preconditioner.coarsesolver.tol"), 1e-1); + prm.put(at_(at, "preconditioner.coarsesolver.solver"), "loopsolver"s); + prm.put(at_(at, "preconditioner.coarsesolver.verbosity"), 0); + prm.put(at_(at, "preconditioner.coarsesolver.preconditioner.type"), "amg"s); + setupDuneAMG(prm, at_(at, "preconditioner.coarsesolver.preconditioner.")); } void setupSystemWellSolver(PropertyTree& prm, @@ -587,10 +593,10 @@ void setupSystemWellSolver(PropertyTree& prm, const FlowLinearSolverParameters& p) { using namespace std::string_literals; - prm.put(at + ".maxiter", 1); - prm.put(at + ".tol", p.linear_solver_reduction_); - prm.put(at + ".verbosity", 0); - prm.put(at + ".solver", "umfpack"s); + prm.put(at_(at, "maxiter"), 1); + prm.put(at_(at, "tol"), p.linear_solver_reduction_); + prm.put(at_(at, "verbosity"), 0); + prm.put(at_(at, "solver"), "umfpack"s); } } // anonymous namespace @@ -598,7 +604,7 @@ void setupSystemWellSolver(PropertyTree& prm, namespace { // The knobs the CPRW pressure stage reads, shared by both layouts. -void setupSystemCPRWellOptions(PropertyTree& prm) +void setupSystemCPRWellOptions(PropertyTree& prm, const std::string& at = "preconditioner.") { using namespace std::string_literals; // How the well equations are contracted to the one coarse unknown each @@ -614,7 +620,7 @@ void setupSystemCPRWellOptions(PropertyTree& prm) // multisegment wells; do not make it the default again // without re-checking them. // unit - the pressure row as-is; a debugging baseline. - prm.put("preconditioner.well_weight_type", "cellavg"s); + prm.put(at + "well_weight_type", "cellavg"s); // How the well unknowns take part in the pressure-stage transfer: // full - restrict the well residual, prolong the bhp correction // no_prolongation - restrict, but discard the bhp correction @@ -627,10 +633,10 @@ void setupSystemCPRWellOptions(PropertyTree& prm) // nearly annihilates the well residual; restricting it may well pay once // the well solve is inexact. // Only read when add_wells. - prm.put("preconditioner.well_transfer", "classic"s); + prm.put(at + "well_transfer", "classic"s); // Give a pressure-controlled well a trivial coarse equation, matching // StandardWellEquations::extractCPRPressureMatrix. Only read when add_wells. - prm.put("preconditioner.well_identity_on_pressure_control", "true"s); + prm.put(at + "well_identity_on_pressure_control", "true"s); // How a well's coarse diagonal is formed: // auto - contract D for single-block wells, minus the row sum for // multisegment ones, i.e. what classic cprw does @@ -639,7 +645,7 @@ void setupSystemCPRWellOptions(PropertyTree& prm) // contract_d is the default. On full Norne with one segment per connection // it is worth ~0.4% over row_sum (2569 against 2580), and it is what makes // the classic cprw path 2646 rather than 2716 on the same case. - prm.put("preconditioner.well_coarse_diagonal", "contract_d"s); + prm.put(at + "well_coarse_diagonal", "contract_d"s); } @@ -691,7 +697,11 @@ setupGeneralSystemCPR(const std::string& conf, const FlowLinearSolverParameters& prm.put("solver", getSolverString(p)); prm.put("preconditioner.type", "general_system_cpr"s); - setupSystemCPRWellOptions(prm); + prm.put("preconditioner.verbosity", 0); + // Sweeps before and after the coarse correction, as for cpr. 0/1 is the + // composition the fixed three-stage system_cpr uses. + prm.put("preconditioner.pre_smooth", 0); + prm.put("preconditioner.post_smooth", 1); // What to do when a well change introduces something the initial build of // the structure did not have: // rebuild - build every part again from this tree @@ -702,36 +712,83 @@ setupGeneralSystemCPR(const std::string& conf, const FlowLinearSolverParameters& // CPR reservoir solve keeps its hierarchy when refreshed and re-aggregates // it when rebuilt. prm.put("preconditioner.well_structure_update", "rebuild"s); + // The reservoir weighting, at the top of the preconditioner as it is for + // cpr: every part that contracts the reservoir equations uses it. + prm.put("preconditioner.weight_type", "trueimpes"s); - // Each stage gets the well solver the fixed layout shares between its - // stages, so the sequence is coarse -> well -> smoother -> well. - setupSystemReservoirSolver(prm, "preconditioner.coarse_solver.reservoir_solver"s, p, add_wells); - setupSystemWellSolver(prm, "preconditioner.coarse_solver.well_solver"s, p); - setupSystemReservoirSmoother(prm, "preconditioner.smoother.reservoir_smoother"s, p); - setupSystemWellSolver(prm, "preconditioner.smoother.well_solver"s, p); + if (add_wells) { + // The coarse space, and the solver for it. As for cpr, this node is + // itself the solver spec; the keys beside it describe the space. + prm.put("preconditioner.coarsesolver.type", "cprw_pressure"s); + setupSystemCPRWellOptions(prm, "preconditioner.coarsesolver."s); + prm.put("preconditioner.coarsesolver.maxiter", 1); + prm.put("preconditioner.coarsesolver.tol", 1e-1); + prm.put("preconditioner.coarsesolver.solver", "loopsolver"s); + prm.put("preconditioner.coarsesolver.verbosity", 0); + prm.put("preconditioner.coarsesolver.preconditioner.type", "amg"s); + setupDuneAMG(prm, "preconditioner.coarsesolver.preconditioner."s); + } + + // The sweep: a reservoir smoother followed by a well solve, which together + // with the coarse correction is what system_cpr does. Unlike cpr's single + // finesmoother this is a list, because a coupled system has more than one + // block to visit and the order matters. + PropertyTree res; + setupSystemReservoirSmoother(res, ""s, p); + res.put("block", "reservoir"s); + PropertyTree well; + setupSystemWellSolver(well, ""s, p); + well.put("block", "well"s); + + if (add_wells) { + // coarse -> well -> reservoir smoother -> well, the sequence the fixed + // three-stage system_cprw applies. The leading well solve cleans up + // after the coarse correction. + prm.put_child_list("preconditioner.finesmoother.steps", {well, res, well}); + } else { + // No coarse space: the reservoir CPR solve is an ordinary block step + // and leads the sweep instead. + PropertyTree resSolver; + setupSystemReservoirSolver(resSolver, ""s, p, /*add_wells=*/false); + resSolver.put("block", "reservoir"s); + prm.put_child_list("preconditioner.finesmoother.steps", + {resSolver, well, res, well}); + } return prm; } - void validateGeneralSystemCPRTree(const PropertyTree& prm) { // Only the composition is checked here: which parts exist is the point of // this layout, so nearly everything is optional. The parts themselves are // validated by whatever builds them. - if (!prm.get_child_optional("preconditioner.coarse_solver").has_value() - && !prm.get_child_optional("preconditioner.smoother").has_value()) { + const auto coarse = prm.get_child_optional("preconditioner.coarsesolver"); + const auto smoother = prm.get_child_optional("preconditioner.finesmoother"); + if (!coarse.has_value() && !smoother.has_value()) { OPM_THROW(std::invalid_argument, "general_system_cpr JSON configuration needs at least one of the " - "'preconditioner.coarse_solver' and 'preconditioner.smoother' sub-trees."); + "'preconditioner.coarsesolver' and 'preconditioner.finesmoother' sub-trees."); } - const auto coarse = prm.get_child_optional("preconditioner.coarse_solver.reservoir_solver"); - if (coarse && coarse->get("preconditioner.add_wells", false) - && !coarse->get_child_optional("preconditioner.coarsesolver").has_value()) { + if (coarse && !coarse->get_child_optional("preconditioner").has_value()) { OPM_THROW(std::invalid_argument, - "In general_system_cpr configuration with " - "preconditioner.coarse_solver.reservoir_solver.preconditioner.add_wells = true, " - "the '...reservoir_solver.preconditioner.coarsesolver' sub-tree is required: " - "it configures the solver for the CPRW pressure system."); + "general_system_cpr's 'preconditioner.coarsesolver' is the solver for the " + "CPRW pressure system and needs its own 'preconditioner' sub-tree."); + } + if (smoother) { + const auto steps = smoother->get_child_list("steps"); + if (!steps.has_value() || steps->empty()) { + OPM_THROW(std::invalid_argument, + "general_system_cpr's 'preconditioner.finesmoother' needs a non-empty " + "'steps' array, each entry naming a \"block\" of 'reservoir' or 'well'."); + } + for (const auto& step : *steps) { + const auto block = step.get("block", std::string{}); + if (block != "reservoir" && block != "well") { + OPM_THROW(std::invalid_argument, + "Every 'preconditioner.finesmoother.steps' entry needs \"block\" set " + "to 'reservoir' or 'well'; got '" + block + "'."); + } + } } } diff --git a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp index 247842c80ad..9d78ccaa842 100644 --- a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp +++ b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp @@ -191,9 +191,7 @@ class GeneralSystemPreconditioner if (coarseResPrm_ && change != WellStructureChange::Values) { weights_ = weightsCalculator_(); buildCoarse(); - twoLevel_ = std::make_unique(*fineOp_, sweep_, *transfer_, - *coarsePolicy_, - /*preSteps=*/0, /*postSteps=*/1); + makeTwoLevel(); } } @@ -234,6 +232,8 @@ class GeneralSystemPreconditioner // Kept so the coarse level can be built again when the well count changes. std::optional coarseResPrm_; WellStructureUpdate update_ = WellStructureUpdate::Rebuild; + std::size_t preSteps_ = 0; + std::size_t postSteps_ = 1; void build() { @@ -255,80 +255,91 @@ class GeneralSystemPreconditioner resWeightCalc = [calc]() { return calc()[_0]; }; } - const auto coarse = prm_.get_child_optional("coarse_solver"); - const auto smoother = prm_.get_child_optional("smoother"); + // Same shape as cpr: an optional coarsesolver, and a finesmoother. + // Unlike cpr the smoother is a list, since a sweep over a coupled + // system has more than one block to visit and the order matters. + const auto coarse = prm_.get_child_optional("coarsesolver"); + const auto smoother = prm_.get_child_optional("finesmoother"); if (!coarse && !smoother) { OPM_THROW(std::invalid_argument, "general_system_cpr needs at least one of the sub-trees " - "'coarse_solver' and 'smoother'."); + "'coarsesolver' and 'finesmoother'."); } - std::vector>> steps; - const auto addReservoir = [&](const PropertyTree& p) { - steps.push_back(std::make_unique( - S_, p, resWeightCalc, pressureIndex_, resComm_)); - }; - const auto addWell = [&](const PropertyTree& p) { - steps.push_back(std::make_unique(S_, p)); - }; - - bool twoLevelCoarse = false; if (coarse) { - if (const auto res = coarse->get_child_optional("reservoir_solver")) { - // add_wells is the same switch the classic CPR/CPRW pair uses: - // it promotes the pressure stage from reservoir-only CPR to - // CPRW over the full (reservoir, well) system, and only then - // is there a coarse space for the whole system. - if (res->get("preconditioner.add_wells", false)) { - coarseResPrm_ = *res; - buildCoarse(); - twoLevelCoarse = true; - } else { - addReservoir(*res); - } - } - if (const auto well = coarse->get_child_optional("well_solver")) { - addWell(*well); + const auto type = coarse->get("type", std::string{"cprw_pressure"}); + if (type != "cprw_pressure") { + OPM_THROW(std::invalid_argument, + "Unknown coarsesolver type '" + type + + "'. The only coarse space on the coupled system is " + "'cprw_pressure'; leave coarsesolver out for none."); } + coarseResPrm_ = *coarse; + buildCoarse(); } + + std::vector>> steps; if (smoother) { - if (const auto res = smoother->get_child_optional("reservoir_smoother")) { - addReservoir(*res); + const auto list = smoother->get_child_list("steps"); + if (!list || list->empty()) { + OPM_THROW(std::invalid_argument, + "general_system_cpr's finesmoother needs a non-empty 'steps' array."); } - if (const auto well = smoother->get_child_optional("well_solver")) { - addWell(*well); + for (const auto& step : *list) { + const auto block = step.get("block", std::string{}); + if (block == "reservoir") { + steps.push_back(std::make_unique( + S_, step, resWeightCalc, pressureIndex_, resComm_)); + } else if (block == "well") { + steps.push_back(std::make_unique(S_, step)); + } else { + OPM_THROW(std::invalid_argument, + "A finesmoother step needs \"block\" set to 'reservoir' or 'well', got '" + + block + "'."); + } } } - if (steps.empty() && !twoLevelCoarse) { + if (steps.empty() && !coarseResPrm_) { OPM_THROW(std::invalid_argument, "general_system_cpr was configured with no parts at all."); } sweep_ = std::make_shared(std::move(steps), isParallel); - if (twoLevelCoarse) { - twoLevel_ = std::make_unique(*fineOp_, sweep_, *transfer_, - *coarsePolicy_, - /*preSteps=*/0, /*postSteps=*/1); + // How many sweeps run before and after the coarse correction. The + // default 0/1 is the composition the fixed three-stage preconditioner + // uses: coarse first, then one sweep. A pre-sweep sees the original + // defect and feeds a coarse correction built from what it leaves. + preSteps_ = prm_.get("pre_smooth", 0); + postSteps_ = prm_.get("post_smooth", 1); + + if (coarseResPrm_) { + makeTwoLevel(); } } + void makeTwoLevel() + { + twoLevel_ = std::make_unique(*fineOp_, sweep_, *transfer_, *coarsePolicy_, + preSteps_, postSteps_); + } + void buildCoarse() { - const auto& resPrm = *coarseResPrm_; - const auto coarsePrm = resPrm.get_child_optional("preconditioner.coarsesolver") - ? resPrm.get_child("preconditioner.coarsesolver") - : PropertyTree(); + // The coarsesolver node is itself the solver spec for the coarse + // pressure system, exactly as it is for cpr; the extra keys beside it + // describe the coarse space and are ignored by FlexibleSolver. + const auto& coarsePrm = *coarseResPrm_; const auto wellTransfer = wellTransferFromString( - prm_.get("well_transfer", std::string{"full"})); + coarsePrm.get("well_transfer", std::string{"full"})); const auto diagonal = wellCoarseDiagonalFromString( - prm_.get("well_coarse_diagonal", std::string{"contract_d"})); + coarsePrm.get("well_coarse_diagonal", std::string{"contract_d"})); if (!weightsCalculator_) { OPM_THROW(std::invalid_argument, - "The CPRW pressure stage (add_wells) needs a weights calculator, but " - "none was configured. Set the reservoir solver's weight_type."); + "The CPRW pressure stage needs a weights calculator, but none was " + "configured. Set coarsesolver.weight_type."); } weights_ = weightsCalculator_(); diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 69842a2750d..932317095b1 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -220,11 +220,14 @@ class ISTLSolverSystem : public ISTLSolver const bool needStructureRefresh = !sysInitialized_ || globalStructureChanged; const auto& prm = this->prm_[this->activeSolverNum_]; - wellWeightType_ = prm.get("preconditioner.well_weight_type", std::string{"cellavg"}); + // The keys describing the coarse space live beside the coarse solver + // for general_system_cpr and at the top for system_cpr. + const auto wellOpts = coarseSpaceTree(prm); + wellWeightType_ = wellOpts.get("well_weight_type", std::string{"cellavg"}); // Give a pressure-controlled well a trivial coarse equation, as the // classic CPRW does. Off keeps the contracted equation for every well. wellLayout_.identityOnPressureControl - = prm.get("preconditioner.well_identity_on_pressure_control", false); + = wellOpts.get("well_identity_on_pressure_control", false); if (needStructureRefresh) { OPM_TIMEBLOCK(flexibleSolverCreate); @@ -428,23 +431,32 @@ class ISTLSolverSystem : public ISTLSolver } } - // Where the reservoir sub-solver sits depends on the preconditioner: the - // fixed three-stage one keeps it at the top, the general one nests it under - // the coarse solver. Only the weights are read from it here. - Opm::PropertyTree reservoirSolverTree(const Opm::PropertyTree& prm) const + // The keys describing the coarse space and its weighting. general_system_cpr + // keeps them beside the coarse solver; system_cpr keeps them at the top of + // the preconditioner and its weighting inside the reservoir solver. + Opm::PropertyTree coarseSpaceTree(const Opm::PropertyTree& prm) const { - if (auto general = prm.get_child_optional("preconditioner.coarse_solver.reservoir_solver")) { + if (auto general = prm.get_child_optional("preconditioner.coarsesolver")) { return *general; } - return prm.get_child("preconditioner.reservoir_solver"); + return prm.get_child("preconditioner"); } void createSystemSolver(const Opm::PropertyTree& prm) { - // Derive weights from the reservoir sub-block config (which uses CPR internally) - auto resSolverPrm = reservoirSolverTree(prm); - std::function()> resWeightCalc - = this->getWeightsCalculator(resSolverPrm, this->getMatrix(), pressureIndex); + std::function()> resWeightCalc; + if (prm.get("preconditioner.type", std::string{}) == "general_system_cpr") { + // The general layout names the weighting directly, at the top of + // the preconditioner as cpr does, rather than through a nested CPR + // preconditioner sub-tree. + resWeightCalc = this->makeWeightsCalculator( + prm.get("preconditioner.weight_type", std::string{"trueimpes"}), + this->getMatrix(), pressureIndex); + } else { + // Derive weights from the reservoir sub-block config (which uses CPR internally) + resWeightCalc = this->getWeightsCalculator( + prm.get_child("preconditioner.reservoir_solver"), this->getMatrix(), pressureIndex); + } // The well part of the weights is filled here too: the CPRW pressure // stage restricts the well rows with it, and re-reads it on every diff --git a/tests/test_GeneralSystemPreconditioner.cpp b/tests/test_GeneralSystemPreconditioner.cpp index baf8b461143..89140ee9d96 100644 --- a/tests/test_GeneralSystemPreconditioner.cpp +++ b/tests/test_GeneralSystemPreconditioner.cpp @@ -172,35 +172,40 @@ struct Fixture // One "apply once" sub-solve: run the configured preconditioner a single time // without a Krylov method around it. +std::string at_(const std::string& at, const std::string& leaf) +{ + return at.empty() ? leaf : at + "." + leaf; +} + void putOnceSolver(Opm::PropertyTree& prm, const std::string& at, const std::string& type) { - prm.put(at + ".maxiter", 1); - prm.put(at + ".tol", 1e-2); - prm.put(at + ".verbosity", 0); - prm.put(at + ".solver", std::string{"preconditioner2inverseoperator"}); - prm.put(at + ".preconditioner.type", type); - prm.put(at + ".preconditioner.relaxation", 1.0); + prm.put(at_(at, "maxiter"), 1); + prm.put(at_(at, "tol"), 1e-2); + prm.put(at_(at, "verbosity"), 0); + prm.put(at_(at, "solver"), std::string{"preconditioner2inverseoperator"}); + prm.put(at_(at, "preconditioner.type"), type); + prm.put(at_(at, "preconditioner.relaxation"), 1.0); } void putReservoirSolver(Opm::PropertyTree& prm, const std::string& at, const bool addWells) { putOnceSolver(prm, at, "cpr"); - prm.put(at + ".preconditioner.use_well_weights", std::string{"false"}); - prm.put(at + ".preconditioner.add_wells", addWells ? std::string{"true"} : std::string{"false"}); - prm.put(at + ".preconditioner.weight_type", std::string{"trueimpes"}); - prm.put(at + ".preconditioner.pre_smooth", 0); - prm.put(at + ".preconditioner.post_smooth", 0); - prm.put(at + ".preconditioner.finesmoother.type", std::string{"jac"}); - prm.put(at + ".preconditioner.finesmoother.relaxation", 1.0); - prm.put(at + ".preconditioner.verbosity", 0); + prm.put(at_(at, "preconditioner.use_well_weights"), std::string{"false"}); + prm.put(at_(at, "preconditioner.add_wells"), addWells ? std::string{"true"} : std::string{"false"}); + prm.put(at_(at, "preconditioner.weight_type"), std::string{"trueimpes"}); + prm.put(at_(at, "preconditioner.pre_smooth"), 0); + prm.put(at_(at, "preconditioner.post_smooth"), 0); + prm.put(at_(at, "preconditioner.finesmoother.type"), std::string{"jac"}); + prm.put(at_(at, "preconditioner.finesmoother.relaxation"), 1.0); + prm.put(at_(at, "preconditioner.verbosity"), 0); // The coarse pressure system is tiny here, so a single ILU0 apply is both // adequate and free of any optional dependency. - prm.put(at + ".preconditioner.coarsesolver.maxiter", 1); - prm.put(at + ".preconditioner.coarsesolver.tol", 1e-1); - prm.put(at + ".preconditioner.coarsesolver.solver", std::string{"preconditioner2inverseoperator"}); - prm.put(at + ".preconditioner.coarsesolver.verbosity", 0); - prm.put(at + ".preconditioner.coarsesolver.preconditioner.type", std::string{"ilu0"}); - prm.put(at + ".preconditioner.coarsesolver.preconditioner.relaxation", 1.0); + prm.put(at_(at, "preconditioner.coarsesolver.maxiter"), 1); + prm.put(at_(at, "preconditioner.coarsesolver.tol"), 1e-1); + prm.put(at_(at, "preconditioner.coarsesolver.solver"), std::string{"preconditioner2inverseoperator"}); + prm.put(at_(at, "preconditioner.coarsesolver.verbosity"), 0); + prm.put(at_(at, "preconditioner.coarsesolver.preconditioner.type"), std::string{"ilu0"}); + prm.put(at_(at, "preconditioner.coarsesolver.preconditioner.relaxation"), 1.0); } void putWellOptions(Opm::PropertyTree& prm, const std::string& transfer) @@ -211,6 +216,8 @@ void putWellOptions(Opm::PropertyTree& prm, const std::string& transfer) prm.put("verbosity", 0); } +// putOnceSolver at the root of its own tree, for a steps entry. + Opm::PropertyTree fixedTree(const bool addWells, const std::string& transfer) { Opm::PropertyTree prm; @@ -222,21 +229,52 @@ Opm::PropertyTree fixedTree(const bool addWells, const std::string& transfer) return prm; } +// One entry of finesmoother.steps: a solver spec plus the block it corrects. +Opm::PropertyTree step(const std::string& block, const std::string& type) +{ + Opm::PropertyTree s; + putOnceSolver(s, "", type); + s.put("block", block); + return s; +} + Opm::PropertyTree generalTree(const bool addWells, const std::string& transfer, const bool coarseWellSolve = true, const bool smootherWellSolve = true) { Opm::PropertyTree prm; prm.put("type", std::string{"general_system_cpr"}); - putWellOptions(prm, transfer); - putReservoirSolver(prm, "coarse_solver.reservoir_solver", addWells); + prm.put("verbosity", 0); + + std::vector steps; + if (addWells) { + // The coarse space and its solver, with the transfer knobs beside it. + prm.put("coarsesolver.type", std::string{"cprw_pressure"}); + prm.put("coarsesolver.weight_type", std::string{"trueimpes"}); + prm.put("coarsesolver.well_transfer", transfer); + prm.put("coarsesolver.well_coarse_diagonal", std::string{"contract_d"}); + prm.put("coarsesolver.well_identity_on_pressure_control", std::string{"false"}); + prm.put("coarsesolver.maxiter", 1); + prm.put("coarsesolver.tol", 1e-1); + prm.put("coarsesolver.solver", std::string{"preconditioner2inverseoperator"}); + prm.put("coarsesolver.verbosity", 0); + prm.put("coarsesolver.preconditioner.type", std::string{"ilu0"}); + prm.put("coarsesolver.preconditioner.relaxation", 1.0); + } else { + // No coarse space: the reservoir CPR solve leads the sweep instead. + auto res = step("reservoir", "cpr"); + putReservoirSolver(res, "", false); + res.put("block", std::string{"reservoir"}); + steps.push_back(res); + } if (coarseWellSolve) { - putOnceSolver(prm, "coarse_solver.well_solver", "ilu0"); + steps.push_back(step("well", "ilu0")); } - putOnceSolver(prm, "smoother.reservoir_smoother", "ilu0"); + steps.push_back(step("reservoir", "ilu0")); if (smootherWellSolve) { - putOnceSolver(prm, "smoother.well_solver", "ilu0"); + steps.push_back(step("well", "ilu0")); } + prm.put_child_list("finesmoother.steps", steps); return prm; } From 84c601ba19d640ba2642aac5788ca065fe98f0bc Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 11 Aug 2026 10:34:39 +0200 Subject: [PATCH 20/22] Address review: match the shipped well_transfer default in both paths Co-Authored-By: Claude Opus 5 --- opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp | 3 ++- opm/simulators/linalg/system/SystemPreconditioner.hpp | 4 +++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp index 9d78ccaa842..df2b6a12f69 100644 --- a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp +++ b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp @@ -332,7 +332,8 @@ class GeneralSystemPreconditioner // describe the coarse space and are ignored by FlexibleSolver. const auto& coarsePrm = *coarseResPrm_; const auto wellTransfer = wellTransferFromString( - coarsePrm.get("well_transfer", std::string{"full"})); + // Same default as setupPropertyTree ships for general_system_cpr. + coarsePrm.get("well_transfer", std::string{"classic"})); const auto diagonal = wellCoarseDiagonalFromString( coarsePrm.get("well_coarse_diagonal", std::string{"contract_d"})); diff --git a/opm/simulators/linalg/system/SystemPreconditioner.hpp b/opm/simulators/linalg/system/SystemPreconditioner.hpp index 2b4f83fc12c..05305b0dc8a 100644 --- a/opm/simulators/linalg/system/SystemPreconditioner.hpp +++ b/opm/simulators/linalg/system/SystemPreconditioner.hpp @@ -325,7 +325,9 @@ class SystemPreconditioner : public Dune::PreconditionerWithUpdate(S_, coarseprm, pressureIndex_, From 8fb69d8eb6fad20b95fada2c8bc1360628dc9fc4 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 11 Aug 2026 11:22:56 +0200 Subject: [PATCH 21/22] Keep ISTLSolver untouched, and let the coarse-matrix dump fire The general layout names the weight type at the top of the preconditioner rather than in a nested CPR sub-tree. Rather than give ISTLSolver a second entry point for that, the derived solver hands the base class the two keys its existing public getWeightsCalculator() reads, so the shared cpr/cprw path is byte-identical to master. Also pass the solver verbosity down: the pressure stage dumps its coarse matrix above 10, but the sub-tree carried a hardcoded 0. Co-Authored-By: Claude Opus 5 --- opm/simulators/linalg/ISTLSolver.hpp | 24 +++++-------------- opm/simulators/linalg/setupPropertyTree.cpp | 5 +++- .../linalg/system/ISTLSolverSystem.hpp | 15 +++++++----- 3 files changed, 19 insertions(+), 25 deletions(-) diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index bcdf4ddad5e..109d52b90cb 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -596,30 +596,18 @@ std::unique_ptr blockJacobiAdjacency(const Grid& grid, const Matrix& matrix, std::size_t pressIndex) const { + std::function weightsCalculator; + using namespace std::string_literals; auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s); // We use lower case as the internal canonical representation of solver names std::ranges::transform(preconditionerType, preconditionerType.begin(), ::tolower); - if (preconditionerType != "cpr" && preconditionerType != "cprt" - && preconditionerType != "cprw" && preconditionerType != "cprwt") { - return {}; - } - const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt"; - return makeWeightsCalculator(prm.get("preconditioner.weight_type"s, "quasiimpes"s), - matrix, pressIndex, transpose); - } - - // The same, for a caller that knows the weighting directly rather than - // through a CPR preconditioner sub-tree. - std::function makeWeightsCalculator(const std::string& weightsType, - const Matrix& matrix, - std::size_t pressIndex, - const bool transpose = false) const - { - std::function weightsCalculator; - { + if (preconditionerType == "cpr" || preconditionerType == "cprt" + || preconditionerType == "cprw" || preconditionerType == "cprwt") { + const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt"; const bool enableThreadParallel = this->parameters_[0].cpr_weights_thread_parallel_; + const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s); if (weightsType == "quasiimpes") { // weights will be created as default in the solver // assignment p = pressureIndex prevent compiler warning about diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 707ad8a3d40..dfe1f08ad85 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -666,6 +666,7 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) // Top-level preconditioner: system_cpr prm.put("preconditioner.type", "system_cpr"s); + prm.put("preconditioner.verbosity", p.linear_solver_verbosity_); setupSystemCPRWellOptions(prm); setupSystemReservoirSmoother(prm, "preconditioner.reservoir_smoother"s, p); @@ -697,7 +698,9 @@ setupGeneralSystemCPR(const std::string& conf, const FlowLinearSolverParameters& prm.put("solver", getSolverString(p)); prm.put("preconditioner.type", "general_system_cpr"s); - prm.put("preconditioner.verbosity", 0); + // The pressure stage dumps its coarse matrix above 10, same convention as + // the solver-level dump; hardcoding 0 here made that unreachable. + prm.put("preconditioner.verbosity", p.linear_solver_verbosity_); // Sweeps before and after the coarse correction, as for cpr. 0/1 is the // composition the fixed three-stage system_cpr uses. prm.put("preconditioner.pre_smooth", 0); diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 932317095b1..2f96672a831 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -446,12 +446,15 @@ class ISTLSolverSystem : public ISTLSolver { std::function()> resWeightCalc; if (prm.get("preconditioner.type", std::string{}) == "general_system_cpr") { - // The general layout names the weighting directly, at the top of - // the preconditioner as cpr does, rather than through a nested CPR - // preconditioner sub-tree. - resWeightCalc = this->makeWeightsCalculator( - prm.get("preconditioner.weight_type", std::string{"trueimpes"}), - this->getMatrix(), pressureIndex); + // The general layout names the weighting directly at the top of the + // preconditioner rather than through a nested CPR sub-tree, so hand + // the base class the two keys it reads instead of teaching it a + // second entry point. + Opm::PropertyTree cprPrm; + cprPrm.put("preconditioner.type", std::string{"cpr"}); + cprPrm.put("preconditioner.weight_type", + prm.get("preconditioner.weight_type", std::string{"trueimpes"})); + resWeightCalc = this->getWeightsCalculator(cprPrm, this->getMatrix(), pressureIndex); } else { // Derive weights from the reservoir sub-block config (which uses CPR internally) resWeightCalc = this->getWeightsCalculator( From ef7584a7a45acdb49d4257b026b6826473017e8a Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 11 Aug 2026 11:30:32 +0200 Subject: [PATCH 22/22] Say that the coarse dump is JSON-only, on purpose --linear-solver-verbosity drives the Krylov output and the full-system dump; the coarse dump is a developer aid reached by setting preconditioner.verbosity in a JSON config. Reverts the config-side change from the previous commit. Co-Authored-By: Claude Opus 5 --- opm/simulators/linalg/setupPropertyTree.cpp | 5 +---- opm/simulators/linalg/system/SystemCprwPressureStage.hpp | 6 ++++-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index dfe1f08ad85..707ad8a3d40 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -666,7 +666,6 @@ setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) // Top-level preconditioner: system_cpr prm.put("preconditioner.type", "system_cpr"s); - prm.put("preconditioner.verbosity", p.linear_solver_verbosity_); setupSystemCPRWellOptions(prm); setupSystemReservoirSmoother(prm, "preconditioner.reservoir_smoother"s, p); @@ -698,9 +697,7 @@ setupGeneralSystemCPR(const std::string& conf, const FlowLinearSolverParameters& prm.put("solver", getSolverString(p)); prm.put("preconditioner.type", "general_system_cpr"s); - // The pressure stage dumps its coarse matrix above 10, same convention as - // the solver-level dump; hardcoding 0 here made that unreachable. - prm.put("preconditioner.verbosity", p.linear_solver_verbosity_); + prm.put("preconditioner.verbosity", 0); // Sweeps before and after the coarse correction, as for cpr. 0/1 is the // composition the fixed three-stage system_cpr uses. prm.put("preconditioner.pre_smooth", 0); diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp index 0fe81e5c8dd..4f6fa3f1345 100644 --- a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -576,8 +576,10 @@ class SystemCprwPressureStage } } - // Same convention as the classic path: verbosity above 10 writes the - // coarse system out so the two can be compared entry by entry. + // Developer aid: verbosity above 10 writes the coarse system out so it can + // be compared entry by entry with the classic path. This reads the + // preconditioner sub-tree, which only a JSON configuration sets, so + // --linear-solver-verbosity does not reach here and is not meant to. void dumpCoarseMatrix() const { if (verbosity_ <= 10) {