diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f24ae5d6a6a..ca0d2145e5c 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -516,6 +516,8 @@ 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 tests/test_wellmodel.cpp @@ -691,6 +693,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 +1129,11 @@ 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/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 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..c5bf38702fd 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.cpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.cpp @@ -135,6 +135,8 @@ 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), " + "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 40a0e57c386..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"); + 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/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 diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index b406e32ec30..707ad8a3d40 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; } @@ -241,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") { @@ -297,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"); } @@ -513,10 +528,134 @@ 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 +// can be compared setting for setting. +void setupSystemReservoirSmoother(PropertyTree& prm, + const std::string& at, + const FlowLinearSolverParameters& p) +{ + using namespace std::string_literals; + 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_(at, "solver"), "preconditioner2inverseoperator"s); + prm.put(at_(at, "preconditioner.type"), "paroverilu0"s); + prm.put(at_(at, "preconditioner.relaxation"), 1.0); +} + +void setupSystemReservoirSolver(PropertyTree& prm, + const std::string& at, + const FlowLinearSolverParameters& p, + const bool add_wells) +{ + using namespace std::string_literals; + 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_(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_(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_(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, + const std::string& at, + const FlowLinearSolverParameters& p) +{ + using namespace std::string_literals; + 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 + +namespace { + +// The knobs the CPRW pressure stage reads, shared by both layouts. +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 + // 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). + // 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(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 + // 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(at + "well_transfer", "classic"s); + // Give a pressure-controlled well a trivial coarse equation, matching + // StandardWellEquations::extractCPRPressureMatrix. Only read when add_wells. + 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 + // 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(at + "well_coarse_diagonal", "contract_d"s); + +} + +} // anonymous namespace + 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,54 +666,140 @@ setupSystemCPR([[maybe_unused]] const std::string& conf, const FlowLinearSolverP // Top-level preconditioner: system_cpr prm.put("preconditioner.type", "system_cpr"s); + setupSystemCPRWellOptions(prm); - // --- 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); - prm.put("preconditioner.reservoir_solver.preconditioner.add_wells", "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); + 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); + 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 + // 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); + // 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); + + 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. + 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.coarsesolver' and 'preconditioner.finesmoother' sub-trees."); + } + if (coarse && !coarse->get_child_optional("preconditioner").has_value()) { + OPM_THROW(std::invalid_argument, + "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 + "'."); + } + } + } +} 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"}) { @@ -593,6 +818,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/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. diff --git a/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp new file mode 100644 index 00000000000..df2b6a12f69 --- /dev/null +++ b/opm/simulators/linalg/system/GeneralSystemPreconditioner.hpp @@ -0,0 +1,357 @@ +/* + 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(); + } + } + + // 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) + { + 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(); + makeTwoLevel(); + } + } + + 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_; + WellStructureUpdate update_ = WellStructureUpdate::Rebuild; + std::size_t preSteps_ = 0; + std::size_t postSteps_ = 1; + + void build() + { + twoLevel_.reset(); + coarsePolicy_.reset(); + transfer_.reset(); + fineOp_.reset(); + 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; + if (weightsCalculator_) { + const auto calc = weightsCalculator_; + resWeightCalc = [calc]() { return calc()[_0]; }; + } + + // 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 " + "'coarsesolver' and 'finesmoother'."); + } + + if (coarse) { + 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) { + 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."); + } + 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() && !coarseResPrm_) { + OPM_THROW(std::invalid_argument, + "general_system_cpr was configured with no parts at all."); + } + + sweep_ = std::make_shared(std::move(steps), isParallel); + + // 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() + { + // 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( + // 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"})); + + if (!weightsCalculator_) { + OPM_THROW(std::invalid_argument, + "The CPRW pressure stage needs a weights calculator, but none was " + "configured. Set coarsesolver.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 316440ee306..2f96672a831 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 @@ -27,6 +28,16 @@ #include #include +#include +#include + +#include +#include +#include +#include +#include +#include + namespace Opm { @@ -96,6 +107,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; @@ -122,6 +142,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_; @@ -153,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; @@ -171,6 +200,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 +220,14 @@ class ISTLSolverSystem : public ISTLSolver const bool needStructureRefresh = !sysInitialized_ || globalStructureChanged; const auto& prm = this->prm_[this->activeSolverNum_]; + // 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 + = wellOpts.get("well_identity_on_pressure_control", false); if (needStructureRefresh) { OPM_TIMEBLOCK(flexibleSolverCreate); @@ -197,9 +236,19 @@ class ISTLSolverSystem : public ISTLSolver sysMatrix_.B = &mergedB_; sysMatrix_.C = &mergedC_; sysMatrix_.D = &mergedD_; - cachedWellStructure_ = merger.buildStructure(); - - refreshSystemSolverForChangedWellStructure(prm); + sysMatrix_.wellLayout = &wellLayout_; + + 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); @@ -212,11 +261,146 @@ class ISTLSolverSystem : public ISTLSolver sysMatrix_.B = &mergedB_; sysMatrix_.C = &mergedC_; sysMatrix_.D = &mergedD_; + sysMatrix_.wellLayout = &wellLayout_; sysPrecond_->update(); } } - void refreshSystemSolverForChangedWellStructure(const Opm::PropertyTree& prm) + // 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); + } + + // 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 + // 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. + // 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(); + 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; + } + + 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 (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; + } + } + 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 + // 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, + const WellStructureChange change) { if (!sysInitialized_ || !sysPrecond_) { createSystemSolver(prm); @@ -227,6 +411,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(change); } else { // Rebuild the parallel solver if the parallel preconditioner cannot be updated in-place. createSystemSolver(prm); @@ -237,24 +423,53 @@ class ISTLSolverSystem : public ISTLSolver if (auto* precond = dynamic_cast(sysPrecond_)) { precond->updateForChangedWellStructure(); + } else if (auto* general = dynamic_cast(sysPrecond_)) { + general->updateForChangedWellStructure(change); } else { // Rebuild the solver if the sequential preconditioner cannot be updated in-place createSystemSolver(prm); } } + // 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.coarsesolver")) { + return *general; + } + 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 = prm.get_child("preconditioner.reservoir_solver"); - 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 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( + 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 + // 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(w[_0]); return w; }; } diff --git a/opm/simulators/linalg/system/SystemCprwPressureStage.hpp b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp new file mode 100644 index 00000000000..4f6fa3f1345 --- /dev/null +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -0,0 +1,639 @@ +/* + 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 +#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 +}; + +// 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") { + 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, + const WellCoarseDiagonal diagonal = WellCoarseDiagonal::ContractD, + const int verbosity = 0) + : S_(S) + , prm_(coarseSolverPrm) + , pressureIndex_(pressureIndex) + , wellTransfer_(wellTransfer) + , comm_(comm) + , diagonal_(diagonal) + , verbosity_(verbosity) + { + } + + // 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()); + dumpCoarseMatrix(); + } + + // (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); + + dumpCoarseRhs(); + 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) + { + restrictInto(dRes, dWell, weights, coarseRhs_); + } + + // 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 + { + prolongFrom(coarseSol_, vRes, vWell); + } + + const CoarseMatrix& coarseMatrix() const + { + 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_; + } + + 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: 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()); + 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 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.isPressureControlled(*j)) { + 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; + 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; + } + 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]; + + // 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; + rowSum += el; + } + + if (rowSumDiag) { + // The classic multisegment convention takes the diagonal + // from the row sum and never reads D. + continue; + } + + // Well row, well columns: + // 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()) { + 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; + } + } + + // 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. + if (rowSumDiag) { + (*coarseMatrix_)[wdof][wdof] = -rowSum; + } + auto& diag = (*coarseMatrix_)[wdof][wdof][0][0]; + if (!(std::abs(diag) > 0.0)) { + diag = 1.0; + } + } + } + + // 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; + } + // 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) { + for (std::size_t wb = layout.firstBlock(j); wb < layout.endBlock(j); ++wb) { + vWell[wb][q] = in[numRes + j][0]; + } + } + } + + // 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) { + 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. + 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; + WellCoarseDiagonal diagonal_ = WellCoarseDiagonal::ContractD; + int verbosity_ = 0; + + 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..05305b0dc8a 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,66 @@ 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( + // Same default as setupPropertyTree ships, so a JSON that omits + // the key gets the same preconditioner as the built-in setup. + prm.get("well_transfer", std::string{"classic"})); + const auto diagonal = wellCoarseDiagonalFromString( + prm.get("well_coarse_diagonal", std::string{"contract_d"})); + cprwStage_ = std::make_unique(S_, coarseprm, pressureIndex_, + wellTransfer, resComm_, diagonal, + prm.get("verbosity", 0)); + 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..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 @@ -38,14 +39,16 @@ 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); + }); + + 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); }); } @@ -67,14 +70,16 @@ 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); + }); + + 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); }); } @@ -91,15 +96,19 @@ 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); + }); + + 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/SystemPreconditionerParts.hpp b/opm/simulators/linalg/system/SystemPreconditionerParts.hpp new file mode 100644 index 00000000000..908aca65f46 --- /dev/null +++ b/opm/simulators/linalg/system/SystemPreconditionerParts.hpp @@ -0,0 +1,517 @@ +/* + 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 +#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. +// -------------------------------------------------------------------------- + +// 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 +{ + // 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; + + // 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; + + // 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(); + } + + // 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; +}; + +// -------------------------------------------------------------------------- +// 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_; + 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: + // 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_; + 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: + 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_; + if (this->needRes0_) { + S_.A->mmv(corrRes_, res[_0]); + } + if (this->needRes1_) { + S_.B->mmv(corrRes_, res[_1]); + } + + if (stage_->prolongatesWellPressure()) { + v[_1] += corrWell_; + 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; + 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) + { + // 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 {} + 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 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 diff --git a/opm/simulators/linalg/system/SystemTypes.hpp b/opm/simulators/linalg/system/SystemTypes.hpp index 093f1cd150d..f8cfe665cd6 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,76 @@ 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; + + // 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; + } + + // 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 +161,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/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 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); } diff --git a/tests/options_system_cprw_approx_wells.json b/tests/options_system_cprw_approx_wells.json new file mode 100644 index 00000000000..b885cc0c344 --- /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": "cellavg", + "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..670707ee467 --- /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": "cellavg", + "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..1cbad0c6c8c --- /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": "cellavg", + "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" + } + } +} \ No newline at end of file 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_GeneralSystemPreconditioner.cpp b/tests/test_GeneralSystemPreconditioner.cpp new file mode 100644 index 00000000000..89140ee9d96 --- /dev/null +++ b/tests/test_GeneralSystemPreconditioner.cpp @@ -0,0 +1,434 @@ +/* + 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. +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_(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_(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_(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) +{ + 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); +} + +// 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; + 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; +} + +// 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"}); + 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) { + steps.push_back(step("well", "ilu0")); + } + steps.push_back(step("reservoir", "ilu0")); + if (smootherWellSolve) { + steps.push_back(step("well", "ilu0")); + } + prm.put_child_list("finesmoother.steps", steps); + 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); +} diff --git a/tests/test_SystemCprwPressureStage.cpp b/tests/test_SystemCprwPressureStage.cpp new file mode 100644 index 00000000000..4930a48abb8 --- /dev/null +++ b/tests/test_SystemCprwPressureStage.cpp @@ -0,0 +1,516 @@ +/* + 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]; + } + } + // 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)); + 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 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) +{ + 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..eeef5ad1b44 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,68 @@ 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)); + + // 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