System CPRW: THP well treatment and trivial-row scaling

CMakeLists_files.cmake

@@ -515,6 +515,9 @@ 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_MswCprWellDiagonal.cpp
+  tests/test_SystemCprwPressureStage.cpp
   tests/test_WellMatrixMerger.cpp
   tests/test_WaterSatfuncConsistencyChecks.cpp
   tests/test_wellmodel.cpp
@@ -690,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
@@ -1122,6 +1129,10 @@ 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/SystemPreconditioner.hpp
   opm/simulators/linalg/system/SystemPreconditionerFactory.hpp
   opm/simulators/linalg/system/SystemTypes.hpp

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<Parameters::LinearSolver>
         ("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 "

opm/simulators/linalg/ISTLSolver.hpp

@@ -56,7 +56,9 @@
 
 #include <any>
 #include <cstddef>
+#include <fstream>
 #include <functional>
+#include <iomanip>
 #include <memory>
 #include <set>
 #include <sstream>
@@ -374,6 +376,37 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
             if (isParallel() && type != "paroverilu0") {
                 detail::makeOverlapRowsInvalid(getMatrix(), overlapRows_);
             }
+
+            dumpPrimaryVariables();
+        }
+
+        // Diagnostic (verbosity > 10): the primary variables and their meanings at
+        // every linearization, so two solver variants can be compared cell by cell.
+        void dumpPrimaryVariables() const
+        {
+            if (prm_[activeSolverNum_].template get<int>("verbosity", 0) <= 10) {
+                return;
+            }
+            static int counter = 0;
+            std::ofstream out("primvars_" + std::to_string(counter++) + ".txt");
+            if (!out) {
+                return;
+            }
+            const auto& sol = simulator_.model().solution(/*timeIdx=*/0);
+            out << std::setprecision(16);
+            for (std::size_t i = 0; i < sol.size(); ++i) {
+                const auto& pv = sol[i];
+                out << i;
+                for (unsigned j = 0; j < std::decay_t<decltype(pv)>::dimension; ++j) {
+                    out << ' ' << pv[j];
+                }
+                if constexpr (requires { pv.primaryVarsMeaningGas(); }) {
+                    out << " gas=" << static_cast<int>(pv.primaryVarsMeaningGas())
+                        << " wat=" << static_cast<int>(pv.primaryVarsMeaningWater())
+                        << " pre=" << static_cast<int>(pv.primaryVarsMeaningPressure());
+                }
+                out << '\n';
+            }
         }
 
         void prepare(const SparseMatrixAdapter& M, Vector& b) override
@@ -588,18 +621,30 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
                                                      const Matrix& matrix,
                                                      std::size_t pressIndex) const
         {
-            std::function<Vector()> weightsCalculator;
-
             using namespace std::string_literals;
 
             auto preconditionerType = prm.get("preconditioner.type"s, "cpr"s);
             // We use lower case as the internal canonical representation of solver names
             std::ranges::transform(preconditionerType, preconditionerType.begin(), ::tolower);
-            if (preconditionerType == "cpr" || preconditionerType == "cprt"
-                || preconditionerType == "cprw" || preconditionerType == "cprwt") {
-                const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
+            if (preconditionerType != "cpr" && preconditionerType != "cprt"
+                && preconditionerType != "cprw" && preconditionerType != "cprwt") {
+                return {};
+            }
+            const bool transpose = preconditionerType == "cprt" || preconditionerType == "cprwt";
+            return makeWeightsCalculator(prm.get("preconditioner.weight_type"s, "quasiimpes"s),
+                                         matrix, pressIndex, transpose);
+        }
+
+        // The same, for a caller that knows the weighting directly rather than
+        // through a CPR preconditioner sub-tree.
+        std::function<Vector()> makeWeightsCalculator(const std::string& weightsType,
+                                                      const Matrix& matrix,
+                                                      std::size_t pressIndex,
+                                                      const bool transpose = false) const
+        {
+            std::function<Vector()> weightsCalculator;
+            {
                 const bool enableThreadParallel = this->parameters_[0].cpr_weights_thread_parallel_;
-                const auto weightsType = prm.get("preconditioner.weight_type"s, "quasiimpes"s);
                 if (weightsType == "quasiimpes") {
                     // weights will be created as default in the solver
                     // assignment p = pressureIndex prevent compiler warning about

opm/simulators/linalg/ISTLSolverRuntimeOptionProxy.hpp

@@ -155,13 +155,15 @@ class ISTLSolverRuntimeOptionProxy : public AbstractISTLSolver<GetPropType<TypeT
     void createSolver(const Simulator& simulator, Args&&... args)
     {
         auto linSolverConf = Parameters::Get<Parameters::LinearSolver>();
-        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<std::string>("preconditioner.type", "") == "system_cpr");
+                const auto type = prm.get<std::string>("preconditioner.type", "");
+                useSystemCpr = (type == "system_cpr") || (type == "general_system_cpr");
             } catch (...) {}
         }
         if (useSystemCpr) {

opm/simulators/linalg/PressureBhpTransferPolicy.hpp

@@ -20,9 +20,16 @@
 
 #pragma once
 
+#include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 
 #include <opm/simulators/linalg/matrixblock.hh>
+#include <opm/simulators/linalg/MatrixMarketSpecializations.hpp>
+
+#include <dune/istl/matrixmarket.hh>
+
+#include <fstream>
+#include <string>
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/twolevelmethodcpr.hh>
 
@@ -199,14 +206,46 @@ namespace Opm
         if (prm_.get<bool>("add_wells")) {
             OPM_TIMEBLOCK(cprwAddWellEquation);
             assert(transpose == false); // not implemented
-            bool use_well_weights = prm_.get<bool>("use_well_weights");
-            fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_, use_well_weights);
+            const bool use_well_weights = prm_.get<bool>("use_well_weights");
+            // "contract_d" takes the coarse well diagonal from lambda' D(:,bhp)
+            // for every well; "auto" (the default) leaves the multisegment path
+            // on its historical row-sum diagonal. Standard wells contract D
+            // either way, so this only changes multisegment wells.
+            const auto diagonal = prm_.get<std::string>("well_coarse_diagonal", "auto");
+            if (diagonal != "auto" && diagonal != "contract_d") {
+                OPM_THROW(std::invalid_argument,
+                          "Unknown well_coarse_diagonal '" + diagonal +
+                          "'. Valid values are 'auto' and 'contract_d'.");
+            }
+            const bool contract_d_diagonal = (diagonal == "contract_d");
+            fineOperator.addWellPressureEquations(*coarseLevelMatrix_, weights_,
+                                                  use_well_weights, contract_d_diagonal);
 #ifndef NDEBUG
             std::advance(rowCoarse, fineOperator.getNumberOfExtraEquations());
             assert(rowCoarse == coarseLevelMatrix_->end());
 #endif
 
         }
+        dumpCoarseMatrix();
+    }
+
+    // verbosity above 10 writes the coarse system out, so that it can be
+    // compared entry by entry against another CPRW implementation.
+    void dumpCoarseMatrix() const
+    {
+        if (prm_.get<int>("verbosity", 0) <= 10) {
+            return;
+        }
+        static int counter = 0;
+        const int id = counter++;
+        std::ofstream out("cprw_coarse_" + std::to_string(id) + ".mm");
+        if (out) {
+            Dune::writeMatrixMarket(*coarseLevelMatrix_, out);
+        }
+        std::ofstream wout("cprw_weights_" + std::to_string(id) + ".mm");
+        if (wout) {
+            Dune::writeMatrixMarket(weights_, wout);
+        }
     }
 
     void moveToCoarseLevel(const typename ParentType::FineRangeType& fine) override

opm/simulators/linalg/PropertyTree.cpp

@@ -114,6 +114,36 @@ std::vector<std::string> PropertyTree::get_child_keys() const
     return keys;
 }
 
+void PropertyTree::put_child_list(const std::string& key,
+                                  const std::vector<PropertyTree>& 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<std::vector<PropertyTree>>
+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<PropertyTree> items;
+    items.reserve(subTree->size());
+    for (const auto& item : *subTree) {
+        items.emplace_back(PropertyTree(item.second));
+    }
+
+    return items;
+}
+
 template <typename T>
 std::optional<std::vector<T>>
 PropertyTree::get_child_items_as_vector(const std::string& child) const

opm/simulators/linalg/PropertyTree.hpp

@@ -137,6 +137,31 @@ class PropertyTree
     /// \return Vector of strings containing the names of all immediate children
     std::vector<std::string> 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<std::vector<PropertyTree>>
+    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<PropertyTree>& items);
+
     /// Retrieve node items as linearised vector.
     ///
     /// Assumes that the node's child is an array type of homongeneous

opm/simulators/linalg/WellOperators.hpp

@@ -60,7 +60,8 @@ class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
     using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
     virtual void addWellPressureEquations(PressureMatrix& jacobian,
                                           const X& weights,
-                                          const bool use_well_weights) const = 0;
+                                          const bool use_well_weights,
+                                          const bool contract_d_diagonal) const = 0;
     virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian) const = 0;
     virtual int getNumberOfExtraEquations() const = 0;
 };
@@ -120,10 +121,11 @@ class WellModelAsLinearOperator : public LinearOperatorExtra<X, Y>
 
     void addWellPressureEquations(PressureMatrix& jacobian,
                                   const X& weights,
-                                  const bool use_well_weights) const override
+                                  const bool use_well_weights,
+                                  const bool contract_d_diagonal) const override
     {
         OPM_TIMEBLOCK(addWellPressureEquations);
-        wellMod_.addWellPressureEquations(jacobian, weights, use_well_weights);
+        wellMod_.addWellPressureEquations(jacobian, weights, use_well_weights, contract_d_diagonal);
     }
 
     void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
@@ -197,12 +199,14 @@ class DomainWellModelAsLinearOperator : public WellModelAsLinearOperator<WellMod
 
     void addWellPressureEquations(PressureMatrix& jacobian,
                                   const X& weights,
-                                  const bool use_well_weights) const override
+                                  const bool use_well_weights,
+                                  const bool contract_d_diagonal) const override
     {
         OPM_TIMEBLOCK(addWellPressureEquations);
         this->wellMod_.addWellPressureEquationsDomain(jacobian,
                                                       weights,
                                                       use_well_weights,
+                                                      contract_d_diagonal,
                                                       domainIndex_);
     }
 
@@ -264,10 +268,11 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
 
     void addWellPressureEquations(PressureMatrix& jacobian,
                                   const X& weights,
-                                  const bool use_well_weights) const
+                                  const bool use_well_weights,
+                                  const bool contract_d_diagonal) const
     {
         OPM_TIMEBLOCK(addWellPressureEquations);
-        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
+        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights, contract_d_diagonal);
     }
 
     void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
@@ -359,10 +364,11 @@ class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X
 
     void addWellPressureEquations(PressureMatrix& jacobian,
                                   const X& weights,
-                                  const bool use_well_weights) const
+                                  const bool use_well_weights,
+                                  const bool contract_d_diagonal) const
     {
         OPM_TIMEBLOCK(addWellPressureEquations);
-        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
+        wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights, contract_d_diagonal);
     }
 
     void addWellPressureEquationsStruct(PressureMatrix& jacobian) const