diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 801bc766364..b5008f76c69 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -532,6 +532,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_tpsa_preconditioner.cpp tests/test_tpsa_primaryvariables.cpp tests/test_vfpproperties.cpp + tests/test_SystemCprwPressureStage.cpp tests/test_WellMatrixMerger.cpp tests/test_WaterSatfuncConsistencyChecks.cpp tests/test_wellmodel.cpp @@ -707,6 +708,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 @@ -1142,6 +1147,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp opm/simulators/linalg/Preconditioner2InverseOperator.hpp opm/simulators/linalg/system/MultiComm.hpp + opm/simulators/linalg/system/SystemCprwPressureStage.hpp opm/simulators/linalg/system/SystemPreconditioner.hpp opm/simulators/linalg/system/SystemPreconditionerFactory.hpp opm/simulators/linalg/system/SystemTypes.hpp diff --git a/opm/simulators/linalg/FlowLinearSolverParameters.cpp b/opm/simulators/linalg/FlowLinearSolverParameters.cpp index 629ab7d4ef5..1e93160f58f 100644 --- a/opm/simulators/linalg/FlowLinearSolverParameters.cpp +++ b/opm/simulators/linalg/FlowLinearSolverParameters.cpp @@ -135,6 +135,7 @@ void FlowLinearSolverParameters::registerParameters() ("Scale linear system according to equation scale and primary variable types"); Parameters::Register ("Configuration of solver. Valid options are: cprw (default), system_cpr (CPU-only), " + "system_cprw (CPU-only, system_cpr with the wells in the pressure stage), " "ilu0, dilu, cpr (an alias for cprw), cpr_quasiimpes, " "cpr_trueimpes, cpr_trueimpesanalytic, amg or hybrid (experimental). " "Alternatively, you can request a configuration to be read from a " diff --git a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp index 40a0e57c386..3c6cbe20e41 100644 --- a/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp +++ b/opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp @@ -155,7 +155,7 @@ class ISTLSolverRuntimeOptionProxy : public AbstractISTLSolver(); - bool useSystemCpr = (linSolverConf == "system_cpr"); + bool useSystemCpr = (linSolverConf == "system_cpr") || (linSolverConf == "system_cprw"); if (!useSystemCpr && linSolverConf.size() > 5 && linSolverConf.ends_with(".json") && std::filesystem::exists(linSolverConf)) { diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 56a582e8d4d..265c5db2933 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -297,7 +297,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_ = cprDefaultMaxIter; } @@ -606,9 +608,10 @@ setupUMFPack([[maybe_unused]] const std::string& conf, const FlowLinearSolverPar PropertyTree -setupSystemCPR([[maybe_unused]] const std::string& conf, const FlowLinearSolverParameters& p) +setupSystemCPR(const std::string& conf, const FlowLinearSolverParameters& p) { using namespace std::string_literals; + const bool add_wells = (conf == "system_cprw"); PropertyTree prm; // Outer solver @@ -619,6 +622,45 @@ setupSystemCPR([[maybe_unused]] const std::string& conf, const FlowLinearSolverP // Top-level preconditioner: system_cpr prm.put("preconditioner.type", "system_cpr"s); + // How the well equations are contracted to the one coarse unknown each + // well carries in the CPRW pressure system. Only read when add_wells. + // cellavg - average of the reservoir weights over the well's + // perforated cells, on the conservation equations only. + // This is what cprw does (use_well_weights = false). + // cellblockavg - the same average taken per block row instead of per + // well. Identical to cellavg for a standard well, and + // worse for multisegment wells (Norne per-connection: + // 2604 against 2569). + // quasiimpes - inv(D)^T e_bhp, normalised. Catastrophic for + // multisegment wells; do not make it the default again + // without re-checking them. + // unit - the pressure row as-is; a debugging baseline. + prm.put("preconditioner.well_weight_type", "cellavg"s); + // How the well unknowns take part in the pressure-stage transfer: + // full - restrict the well residual, prolong the bhp correction + // no_prolongation - restrict, but discard the bhp correction + // classic - neither, i.e. the classic cprw formulation, so that + // the only remaining difference is numerics + // classic is the default: on full Norne with one segment per connection it + // measures best (2558, against 2569 for no_prolongation and 2576 for full), + // and on standard wells the three are within one iteration of each other. + // The margin is thin. It is also taken with an exact well solve, which + // nearly annihilates the well residual; restricting it may well pay once + // the well solve is inexact. + // Only read when add_wells. + prm.put("preconditioner.well_transfer", "classic"s); + // Give a pressure-controlled well a trivial coarse equation, matching + // StandardWellEquations::extractCPRPressureMatrix. Only read when add_wells. + prm.put("preconditioner.well_identity_on_pressure_control", "true"s); + // 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("preconditioner.well_coarse_diagonal", "contract_d"s); // --- Reservoir smoother --- prm.put("preconditioner.reservoir_smoother.maxiter", 1); @@ -640,7 +682,10 @@ setupSystemCPR([[maybe_unused]] const std::string& conf, const FlowLinearSolverP prm.put("preconditioner.reservoir_solver.preconditioner.type", "cpr"s); prm.put("preconditioner.reservoir_solver.preconditioner.relaxation", 1.0); prm.put("preconditioner.reservoir_solver.preconditioner.use_well_weights", "false"s); - prm.put("preconditioner.reservoir_solver.preconditioner.add_wells", "false"s); + // add_wells promotes the pressure stage from reservoir-only CPR to CPRW + // over the full (reservoir, well) system. + prm.put("preconditioner.reservoir_solver.preconditioner.add_wells", + add_wells ? "true"s : "false"s); prm.put("preconditioner.reservoir_solver.preconditioner.weight_type", "trueimpes"s); prm.put("preconditioner.reservoir_solver.preconditioner.pre_smooth", 0); prm.put("preconditioner.reservoir_solver.preconditioner.post_smooth", 0); @@ -685,6 +730,69 @@ void validateSystemCPRTree(const PropertyTree& prm) "In system_cpr configuration, the reservoir_solver must use the CPR preconditioner " "(preconditioner.reservoir_solver.preconditioner.type = 'cpr')."); } + const bool addWells = reservoir_solver->get("preconditioner.add_wells", false); + // 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. + const bool hasCoarseSolver + = reservoir_solver->get_child_optional("preconditioner.coarsesolver").has_value(); + if (addWells && !hasCoarseSolver) { + 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."); + } + // Fail here, at setup time, rather than inside + // SystemCprwPressureStage::wellTransferFromString/wellCoarseDiagonalFromString + // on the first linear solve of the run. Valid values mirror those two parsers. + if (addWells) { + const auto wellTransfer + = prm.get("preconditioner.well_transfer", std::string{"classic"}); + const bool wellTransferOk = (wellTransfer == "full") + || (wellTransfer == "no_prolongation") || (wellTransfer == "classic"); + if (!wellTransferOk) { + OPM_THROW(std::invalid_argument, + fmt::format("Unknown preconditioner.well_transfer '{}'. Valid " + "values are 'full', 'no_prolongation' and 'classic'.", + wellTransfer)); + } + const auto wellCoarseDiagonal + = prm.get("preconditioner.well_coarse_diagonal", std::string{"contract_d"}); + const bool wellCoarseDiagonalOk = (wellCoarseDiagonal == "auto") + || (wellCoarseDiagonal == "contract_d") || (wellCoarseDiagonal == "row_sum"); + if (!wellCoarseDiagonalOk) { + OPM_THROW(std::invalid_argument, + fmt::format("Unknown preconditioner.well_coarse_diagonal '{}'. Valid " + "values are 'auto', 'contract_d' and 'row_sum'.", + wellCoarseDiagonal)); + } + } + } + + // A Krylov well solver stops on a tolerance, so it performs a different + // number of inner iterations for each right-hand side. That makes the whole + // system preconditioner non-stationary, which Krylov methods with short + // recurrences (bicgstab, cg) and standard GMRES are not allowed to use: they + // assume a fixed preconditioning operator. The outer solver has to be a + // flexible one. + auto well_solver = prm.get_child_optional("preconditioner.well_solver"); + if (well_solver) { + const auto inner = well_solver->get("solver", "bicgstab"); + const bool inner_is_krylov = (inner == "bicgstab") || (inner == "gmres") + || (inner == "cg") || (inner == "flexgmres"); + const auto outer = prm.get("solver", "bicgstab"); + const bool outer_is_flexible = (outer == "flexgmres"); + if (inner_is_krylov && !outer_is_flexible) { + OPM_THROW(std::invalid_argument, + fmt::format("system_cpr is configured with an approximate (Krylov) well " + "solver, 'preconditioner.well_solver.solver' = '{}', which makes " + "the preconditioner vary between applications. The outer solver " + "must then be flexible: set 'solver' to 'flexgmres' (it is " + "currently '{}'), or use a stationary well solver such as " + "'umfpack' or 'preconditioner2inverseoperator'.", + inner, outer)); + } } } diff --git a/opm/simulators/linalg/system/ISTLSolverSystem.hpp b/opm/simulators/linalg/system/ISTLSolverSystem.hpp index 316440ee306..e23919fb98b 100644 --- a/opm/simulators/linalg/system/ISTLSolverSystem.hpp +++ b/opm/simulators/linalg/system/ISTLSolverSystem.hpp @@ -27,6 +27,20 @@ #include #include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + namespace Opm { @@ -96,6 +110,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 +145,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_; @@ -163,6 +193,24 @@ class ISTLSolverSystem : public ISTLSolver { OPM_TIMEBLOCK(flexibleSolverPrepare); + // Read before buildWellDofLayout(), which consults + // wellLayout_.identityOnPressureControl to decide whether to + // populate wellLayout_.pressureControlled: reading it afterwards + // would leave that flag one call stale. + const auto& prm = this->prm_[this->activeSolverNum_]; + wellWeightType_ = prm.get("preconditioner.well_weight_type", std::string{"cellavg"}); + if (wellWeightType_ != "unit" && wellWeightType_ != "cellavg" + && wellWeightType_ != "cellblockavg" && wellWeightType_ != "quasiimpes") { + OPM_THROW(std::invalid_argument, + "Unknown preconditioner.well_weight_type '" + wellWeightType_ + + "'. Valid values are 'unit', 'cellavg', " + "'cellblockavg' and 'quasiimpes'."); + } + // Give a pressure-controlled well a trivial coarse equation, as the + // classic CPRW does. Off keeps the contracted equation for every well. + wellLayout_.identityOnPressureControl + = prm.get("preconditioner.well_identity_on_pressure_control", false); + wellBMatrices_.clear(); wellCMatrices_.clear(); wellDMatrices_.clear(); @@ -171,6 +219,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_); @@ -188,8 +238,6 @@ class ISTLSolverSystem : public ISTLSolver #endif const bool needStructureRefresh = !sysInitialized_ || globalStructureChanged; - const auto& prm = this->prm_[this->activeSolverNum_]; - if (needStructureRefresh) { OPM_TIMEBLOCK(flexibleSolverCreate); merger.buildMatrices(mergedB_, mergedC_, mergedD_); @@ -197,6 +245,7 @@ class ISTLSolverSystem : public ISTLSolver sysMatrix_.B = &mergedB_; sysMatrix_.C = &mergedC_; sysMatrix_.D = &mergedD_; + sysMatrix_.wellLayout = &wellLayout_; cachedWellStructure_ = merger.buildStructure(); refreshSystemSolverForChangedWellStructure(prm); @@ -212,10 +261,150 @@ class ISTLSolverSystem : public ISTLSolver sysMatrix_.B = &mergedB_; sysMatrix_.C = &mergedC_; sysMatrix_.D = &mergedD_; + sysMatrix_.wellLayout = &wellLayout_; sysPrecond_->update(); } } + // Which merged well block rows belong to which well. The merged D matrix + // is the per-well D blocks concatenated, so this is a plain prefix sum + // over their dimensions: one block for a standard well, one per segment + // for a multisegment well. + void buildWellDofLayout() + { + auto& offsets = wellLayout_.wellBlockOffsets; + offsets.clear(); + offsets.reserve(wellDMatrices_.size() + 1); + offsets.push_back(0); + std::size_t total = 0; + for (const auto& d : wellDMatrices_) { + total += d.N(); + offsets.push_back(total); + } + + // Which wells are on pressure control. Asking this is the outer + // layer's job; below here it is just a flag per well. Index j is the + // position in this rank's BlackoilWellModel well container at this + // linear solve, the same loop addBCDMatrix uses to build the D blocks. + 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); + } + if (wellLayout_.pressureControlled.size() != wellDMatrices_.size()) { + OPM_THROW(std::logic_error, + "System CPRW: the well container and the extracted well " + "matrices disagree on the number of wells."); + } + } + } + + // Weights used to contract each well's equations down to the single scalar + // the CPRW pressure system carries for that well. Computed here rather + // than inside the preconditioner so that the linear-solver core never sees + // anything well-specific, and so that this can later be replaced by a + // value obtained from the well model without touching the core. + WellVector computeWellWeights(const ResVector& resWeights) const + { + const std::size_t numBlocks = mergedD_.N(); + const int q = wellLayout_.pressureDofIndex; + + WellVector weights(numBlocks); + std::optional cellavgWell; + 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 auto well = perWell ? wellLayout_.wellOfBlock(wb) : std::nullopt; + if (perWell && well.has_value() && well == cellavgWell) { + // A well's blocks are contiguous and share one average. + lambda = weights[wb - 1]; + continue; + } + cellavgWell = well; + const std::size_t first = perWell ? wellLayout_.firstBlock(*well) : wb; + const std::size_t last = perWell ? wellLayout_.endBlock(*well) : 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) { if (!sysInitialized_ || !sysPrecond_) { @@ -250,11 +439,15 @@ class ISTLSolverSystem : public ISTLSolver std::function()> resWeightCalc = this->getWeightsCalculator(resSolverPrm, this->getMatrix(), pressureIndex); + // The well part of the weights is filled here too: the CPRW pressure + // stage restricts the well rows with it, and re-reads it on every + // update, so it has to track the current merged D. std::function()> sysWeightCalc; if (resWeightCalc) { - sysWeightCalc = [resWeightCalc]() { + sysWeightCalc = [this, resWeightCalc]() { SystemVector w; w[_0] = resWeightCalc(); + w[_1] = this->computeWellWeights(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..b9a306bd9d0 --- /dev/null +++ b/opm/simulators/linalg/system/SystemCprwPressureStage.hpp @@ -0,0 +1,646 @@ +/* + 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 +#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; + // Sized from A alone; well rows are denser than that and land in the + // implicit-build overflow area, which is what overflowFraction is for. + // A.N() can be zero for a rank that locally owns no reservoir cells; + // guard it rather than feeding NaN to the implicit-build constructor. + const std::size_t averageEntriesPerRow = (A.N() == 0) + ? 0 + : 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 = wellLayout().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 = wellLayout().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 j} sum_i w1[wb][i] * D[wb][wb'][i][q] + // The merged D is block diagonal by well (WellMatrixMerger + // merges each well's own D block, never across wells), so + // every nonzero of D[wb] belongs to well j itself: accumulate + // straight into the diagonal. This 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 = wellLayout().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][wdof] += 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) { + if (layout.isPressureControlled(j)) { + // Identity coarse row (assembleCoarseMatrix): dp_j must be zero. + continue; + } + 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]; + } + } + } + + // Suffix identifying this rank, so that ranks sharing a working directory + // do not clobber each other's dump files. + std::string dumpRankSuffix() const + { + if constexpr (isParallel) { + return "_rank" + std::to_string(comm_->communicator().rank()); + } else { + return ""; + } + } + + // 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++) + + dumpRankSuffix() + ".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++) + + dumpRankSuffix() + ".mm"); + if (out) { + Dune::writeMatrixMarket(coarseRhs_, out); + } + } + + 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..07e8bcddaf7 100644 --- a/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp +++ b/opm/simulators/linalg/system/SystemPreconditionerFactory.cpp @@ -38,14 +38,8 @@ void addSystemCprSeq() [](const O& op, const P& prm, const std::function& sysWeightCalc, std::size_t pressureIndex) { - std::function()> resWeightCalc; - if (sysWeightCalc) { - resWeightCalc = [sysWeightCalc]() { - return sysWeightCalc()[Dune::Indices::_0]; - }; - } return std::make_shared>>( - op.getmat(), resWeightCalc, pressureIndex, prm); + op.getmat(), sysWeightCalc, pressureIndex, prm); }); } @@ -67,14 +61,8 @@ void addSystemCprParSeq() [](const O& op, const P& prm, const std::function& sysWeightCalc, std::size_t pressureIndex) { - std::function()> resWeightCalc; - if (sysWeightCalc) { - resWeightCalc = [sysWeightCalc]() { - return sysWeightCalc()[Dune::Indices::_0]; - }; - } return std::make_shared>>( - op.getmat(), resWeightCalc, pressureIndex, prm); + op.getmat(), sysWeightCalc, pressureIndex, prm); }); } @@ -91,15 +79,9 @@ void addSystemCprPar() const std::function& sysWeightCalc, std::size_t pressureIndex, const Opm::SystemComm& comm) { - std::function()> resWeightCalc; - if (sysWeightCalc) { - resWeightCalc = [sysWeightCalc]() { - return sysWeightCalc()[Dune::Indices::_0]; - }; - } const auto& resComm = comm[Dune::Indices::_0]; return std::make_shared, Opm::ParResComm>>( - op.getmat(), resWeightCalc, pressureIndex, prm, resComm); + op.getmat(), sysWeightCalc, pressureIndex, prm, resComm); }); } #endif diff --git a/opm/simulators/linalg/system/SystemTypes.hpp b/opm/simulators/linalg/system/SystemTypes.hpp index 093f1cd150d..a832fa76eca 100644 --- a/opm/simulators/linalg/system/SystemTypes.hpp +++ b/opm/simulators/linalg/system/SystemTypes.hpp @@ -26,6 +26,12 @@ #include #include +#include +#include +#include +#include +#include + namespace Opm { @@ -58,6 +64,87 @@ 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(); + } + + // Merged well block row -> well index, nullopt outside the layout. + std::optional wellOfBlock(const std::size_t blockRow) const + { + const auto it + = std::upper_bound(wellBlockOffsets.begin(), wellBlockOffsets.end(), blockRow); + if (it == wellBlockOffsets.begin() || it == wellBlockOffsets.end()) { + return std::nullopt; + } + return static_cast(std::distance(wellBlockOffsets.begin(), it) - 1); + } +}; + // -------------------------------------------------------------------------- // 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 +175,10 @@ class SystemMatrix const WRMatrix* B = nullptr; // (1,0) well–reservoir coupling const WWMatrix* D = nullptr; // (1,1) well + // Aggregation of the well block rows into wells. Only needed by the CPRW + // pressure stage; null when the well DOFs are not aggregated. + const WellDofLayout* wellLayout = nullptr; + // Sub-block access: S[_0][_0], S[_0][_1], S[_1][_0], S[_1][_1] inline SystemMatrixRow0 operator[](Dune::index_constant<0>) const; inline SystemMatrixRow1 operator[](Dune::index_constant<1>) const; diff --git a/tests/options_system_cprw_approx_wells.json b/tests/options_system_cprw_approx_wells.json new file mode 100644 index 00000000000..c96f0dcdb60 --- /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" +} 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..7b61377d6bc --- /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" + } + } + } +} diff --git a/tests/options_system_cprw_complete.json b/tests/options_system_cprw_complete.json new file mode 100644 index 00000000000..3d2926d13f1 --- /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" + } + } +} diff --git a/tests/options_system_cprw_missing_coarsesolver.json b/tests/options_system_cprw_missing_coarsesolver.json new file mode 100644 index 00000000000..1636d8fd732 --- /dev/null +++ b/tests/options_system_cprw_missing_coarsesolver.json @@ -0,0 +1,40 @@ +{ + "maxiter": "20", + "tol": "0.005", + "verbosity": "0", + "solver": "bicgstab", + "preconditioner": { + "type": "system_cpr", + "reservoir_smoother": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "paroverilu0", + "relaxation": "1" + } + }, + "reservoir_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "preconditioner2inverseoperator", + "preconditioner": { + "type": "cpr", + "relaxation": "1", + "use_well_weights": "false", + "add_wells": "true", + "weight_type": "trueimpes", + "pre_smooth": "0", + "post_smooth": "0" + } + }, + "well_solver": { + "maxiter": "1", + "tol": "0.005", + "verbosity": "0", + "solver": "umfpack" + } + } +} diff --git a/tests/test_SystemCprwPressureStage.cpp b/tests/test_SystemCprwPressureStage.cpp new file mode 100644 index 00000000000..fe55fc1bfac --- /dev/null +++ b/tests/test_SystemCprwPressureStage.cpp @@ -0,0 +1,548 @@ +/* + 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); + } + } +} + +// A pressure-controlled well gets an identity coarse row, so its coarse +// right-hand side must stay zero: otherwise the coarse solve returns the +// restricted residual as dp_j and full transfer prolongs it into the well. +BOOST_AUTO_TEST_CASE(PressureControlledWellHasZeroCoarseRhs) +{ + Fixture f; + f.layout.identityOnPressureControl = true; + f.layout.pressureControlled = {0, 1}; + + 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; + } + } + + Stage stage(f.S, Opm::PropertyTree(), pressureIndex, Opm::WellTransfer::Full); + stage.buildCoarseSystem(f.weights); + stage.moveToCoarseLevel(dRes, dWell, f.weights); + + BOOST_CHECK_NE(stage.coarseRhs()[numCells + 0][0], 0.0); + BOOST_CHECK_EQUAL(stage.coarseRhs()[numCells + 1][0], 0.0); + BOOST_CHECK_EQUAL(coarseEntry(stage, numCells + 1, numCells + 1), 1.0); +} + +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