Merge pull request #2990 from atgeirr/testing-precfactory-operators

Only add AMG preconditioners to factory if sensible.

CMakeLists_files.cmake

@@ -132,6 +132,8 @@ list (APPEND TEST_DATA_FILES
   tests/wells_no_perforation.data
   tests/matr33.txt
   tests/rhs3.txt
+  tests/matr33rep.txt
+  tests/rhs3rep.txt
   tests/options_flexiblesolver.json
   tests/options_flexiblesolver_simple.json
   )

opm/simulators/linalg/PreconditionerFactory.hpp

@@ -288,17 +288,25 @@ private:
             const double w = prm.get<double>("relaxation", 1.0);
             return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(comm, op.getmat(), n, w);
         });
-        doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
+
-            const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
+        // Only add AMG preconditioners to the factory if the operator
-            if (smoother == "ILU0" || smoother == "ParOverILU0") {
+        // is the overlapping schwarz operator. This could be extended
-                using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
+        // later, but at this point no other operators are compatible
-                auto crit = amgCriterion(prm);
+        // with the AMG hierarchy construction.
-                auto sargs = amgSmootherArgs<Smoother>(prm);
+        if constexpr (std::is_same_v<O, Dune::OverlappingSchwarzOperator<M, V, V, C>>) {
-                return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
+            doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
-            } else {
+                const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
-                OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
+                if (smoother == "ILU0" || smoother == "ParOverILU0") {
-            }
+                    using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
-        });
+                    auto crit = amgCriterion(prm);
+                    auto sargs = amgSmootherArgs<Smoother>(prm);
+                    return std::make_shared<Dune::Amg::AMGCPR<O, V, Smoother, C>>(op, crit, sargs, comm);
+                } else {
+                    OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << ".");
+                }
+            });
+        }
+
         doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
             assert(weightsCalculator);
             return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, comm);
@@ -355,74 +363,79 @@ private:
             const double w = prm.get<double>("relaxation", 1.0);
             return wrapPreconditioner<SeqSSOR<M, V, V>>(op.getmat(), n, w);
         });
-        doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
+
-            const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
+        // Only add AMG preconditioners to the factory if the operator
-            if (smoother == "ILU0" || smoother == "ParOverILU0") {
+        // is an actual matrix operator.
+        if constexpr (std::is_same_v<O, Dune::MatrixAdapter<M, V, V>>) {
+            doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
+                const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
+                if (smoother == "ILU0" || smoother == "ParOverILU0") {
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
-                using Smoother = SeqILU<M, V, V>;
+                    using Smoother = SeqILU<M, V, V>;
 #else
-                using Smoother = SeqILU0<M, V, V>;
+                    using Smoother = SeqILU0<M, V, V>;
 #endif
-                return makeAmgPreconditioner<Smoother>(op, prm);
+                    return makeAmgPreconditioner<Smoother>(op, prm);
-            } else if (smoother == "Jac") {
+                } else if (smoother == "Jac") {
-                using Smoother = SeqJac<M, V, V>;
+                    using Smoother = SeqJac<M, V, V>;
-                return makeAmgPreconditioner<Smoother>(op, prm);
+                    return makeAmgPreconditioner<Smoother>(op, prm);
-            } else if (smoother == "SOR") {
+                } else if (smoother == "SOR") {
-                using Smoother = SeqSOR<M, V, V>;
+                    using Smoother = SeqSOR<M, V, V>;
-                return makeAmgPreconditioner<Smoother>(op, prm);
+                    return makeAmgPreconditioner<Smoother>(op, prm);
-            } else if (smoother == "SSOR") {
+                } else if (smoother == "SSOR") {
-                using Smoother = SeqSSOR<M, V, V>;
+                    using Smoother = SeqSSOR<M, V, V>;
-                return makeAmgPreconditioner<Smoother>(op, prm);
+                    return makeAmgPreconditioner<Smoother>(op, prm);
-            } else if (smoother == "ILUn") {
+                } else if (smoother == "ILUn") {
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
-                using Smoother = SeqILU<M, V, V>;
+                    using Smoother = SeqILU<M, V, V>;
 #else
-                using Smoother = SeqILUn<M, V, V>;
+                            using Smoother = SeqILUn<M, V, V>;
 #endif
-                return makeAmgPreconditioner<Smoother>(op, prm);
+                    return makeAmgPreconditioner<Smoother>(op, prm);
-            } else {
+                } else {
-                OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
+                    OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << ".");
-            }
+                }
-        });
+            });
-	doAddCreator("kamg", [](const O& op, const P& prm, const std::function<Vector()>&) {
+            doAddCreator("kamg", [](const O& op, const P& prm, const std::function<Vector()>&) {
-            const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
+                const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
-            if (smoother == "ILU0" || smoother == "ParOverILU0") {
+                if (smoother == "ILU0" || smoother == "ParOverILU0") {
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
-                using Smoother = SeqILU<M, V, V>;
+                    using Smoother = SeqILU<M, V, V>;
 #else
-                using Smoother = SeqILU0<M, V, V>;
+                        using Smoother = SeqILU0<M, V, V>;
 #endif
-                return makeAmgPreconditioner<Smoother>(op, prm, true);
+                    return makeAmgPreconditioner<Smoother>(op, prm, true);
-            } else if (smoother == "Jac") {
+                } else if (smoother == "Jac") {
-                using Smoother = SeqJac<M, V, V>;
+                    using Smoother = SeqJac<M, V, V>;
-                return makeAmgPreconditioner<Smoother>(op, prm, true);
+                    return makeAmgPreconditioner<Smoother>(op, prm, true);
-            } else if (smoother == "SOR") {
+                } else if (smoother == "SOR") {
-                using Smoother = SeqSOR<M, V, V>;
+                    using Smoother = SeqSOR<M, V, V>;
-                return makeAmgPreconditioner<Smoother>(op, prm, true);
+                    return makeAmgPreconditioner<Smoother>(op, prm, true);
-	    // } else if (smoother == "GS") {
+                    // } else if (smoother == "GS") {
-            //     using Smoother = SeqGS<M, V, V>;
+                    //     using Smoother = SeqGS<M, V, V>;
-            //     return makeAmgPreconditioner<Smoother>(op, prm, true);
+                    //     return makeAmgPreconditioner<Smoother>(op, prm, true);
-            } else if (smoother == "SSOR") {
+                } else if (smoother == "SSOR") {
-                using Smoother = SeqSSOR<M, V, V>;
+                    using Smoother = SeqSSOR<M, V, V>;
-                return makeAmgPreconditioner<Smoother>(op, prm, true);
+                    return makeAmgPreconditioner<Smoother>(op, prm, true);
-            } else if (smoother == "ILUn") {
+                } else if (smoother == "ILUn") {
 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
-                using Smoother = SeqILU<M, V, V>;
+                    using Smoother = SeqILU<M, V, V>;
 #else
-                using Smoother = SeqILUn<M, V, V>;
+                        using Smoother = SeqILUn<M, V, V>;
 #endif
-                return makeAmgPreconditioner<Smoother>(op, prm, true);
+                    return makeAmgPreconditioner<Smoother>(op, prm, true);
-            } else {
+                } else {
-                OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<".");
+                    OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << ".");
-            }
+                }
-        });
+            });
-        doAddCreator("famg", [](const O& op, const P& prm, const std::function<Vector()>&) {
+            doAddCreator("famg", [](const O& op, const P& prm, const std::function<Vector()>&) {
-            auto crit = amgCriterion(prm);
+                auto crit = amgCriterion(prm);
-            Dune::Amg::Parameters parms;
+                Dune::Amg::Parameters parms;
-            parms.setNoPreSmoothSteps(1);
+                parms.setNoPreSmoothSteps(1);
-            parms.setNoPostSmoothSteps(1);
+                parms.setNoPostSmoothSteps(1);
-            return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
+                return wrapPreconditioner<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
-        });
+            });
+        }
         doAddCreator("cpr", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
             return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator);
         });

tests/matr33rep.txt

@@ -0,0 +1,84 @@
+%%MatrixMarket matrix coordinate real general
+% ISTL_STRUCT blocked 3 3
+9 9 81
+1 1 1
+2 1 -1
+3 1 -1
+4 1 -1
+5 1 -1
+6 1 -1
+7 1 -1
+8 1 -1
+9 1 -1
+1 2 -1
+2 2 1
+3 2 -1
+4 2 -1
+5 2 -1
+6 2 -1
+7 2 -1
+8 2 -1
+9 2 -1
+1 3 -1
+2 3 -1
+3 3 1
+4 3 -1
+5 3 -1
+6 3 -1
+7 3 -1
+8 3 -1
+9 3 -1
+1 4 -1
+2 4 -1
+3 4 -1
+4 4 1
+5 4 -1
+6 4 -1
+7 4 -1
+8 4 -1
+9 4 -1
+1 5 -1
+2 5 -1
+3 5 -1
+4 5 -1
+5 5 1
+6 5 -1
+7 5 -1
+8 5 -1
+9 5 -1
+1 6 -1
+2 6 -1
+3 6 -1
+4 6 -1
+5 6 -1
+6 6 1
+7 6 -1
+8 6 -1
+9 6 -1
+1 7 -1
+2 7 -1
+3 7 -1
+4 7 -1
+5 7 -1
+6 7 -1
+7 7 1
+8 7 -1
+9 7 -1
+1 8 -1
+2 8 -1
+3 8 -1
+4 8 -1
+5 8 -1
+6 8 -1
+7 8 -1
+8 8 1
+9 8 -1
+1 9 -1
+2 9 -1
+3 9 -1
+4 9 -1
+5 9 -1
+6 9 -1
+7 9 -1
+8 9 -1
+9 9 1

tests/rhs3rep.txt

@@ -0,0 +1,12 @@
+%%MatrixMarket matrix array real general
+% ISTL_STRUCT blocked 3 1
+9 1
+-1
+-1
+-1
+-3
+-3
+-3
+-3
+-3
+-3

tests/test_preconditionerfactory.cpp

@@ -232,6 +232,157 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner)
     test3(prm);
 }
 
+
+template<class Mat, class Vec>
+class RepeatingOperator : public Dune::AssembledLinearOperator<Mat, Vec, Vec>
+{
+public:
+    using matrix_type = Mat;
+    using domain_type = Vec;
+    using range_type = Vec;
+    using field_type = typename Vec::field_type;
+
+    Dune::SolverCategory::Category category() const override
+    {
+        return Dune::SolverCategory::sequential;
+    }
+
+    RepeatingOperator(const Mat& matrix, const int repeats)
+        : matrix_(matrix)
+        , repeats_(repeats)
+    {
+    }
+
+    // y = A*x;
+    virtual void apply(const Vec& x, Vec& y) const override
+    {
+        y = 0;
+        applyscaleadd(1.0, x, y);
+    }
+
+    // y += \alpha * A * x
+    virtual void applyscaleadd(field_type alpha, const Vec& x, Vec& y) const override
+    {
+        Vec temp1 = x;
+        Vec temp2 = x; // For size.
+        temp2 = 0.0;
+        for (int rr = 0; rr < repeats_; ++rr) {
+            // mv below means: temp2 = matrix_ * temp1;
+            matrix_.mv(temp1, temp2);
+            temp1 = temp2;
+        }
+        temp2 *= alpha;
+        y += temp2;
+    }
+
+    virtual const matrix_type& getmat() const override
+    {
+        return matrix_;
+    }
+
+protected:
+    const Mat& matrix_;
+    const int repeats_;
+};
+
+
+template <int bz>
+Dune::BlockVector<Dune::FieldVector<double, bz>>
+testPrecRepeating(const boost::property_tree::ptree& prm, const std::string& matrix_filename, const std::string& rhs_filename)
+{
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+    Matrix matrix;
+    {
+        std::ifstream mfile(matrix_filename);
+        if (!mfile) {
+            throw std::runtime_error("Could not read matrix file");
+        }
+        readMatrixMarket(matrix, mfile);
+    }
+    Vector rhs;
+    {
+        std::ifstream rhsfile(rhs_filename);
+        if (!rhsfile) {
+            throw std::runtime_error("Could not read rhs file");
+        }
+        readMatrixMarket(rhs, rhsfile);
+    }
+    using Operator = RepeatingOperator<Matrix, Vector>;
+    Operator op(matrix, 2);
+    using PrecFactory = Opm::PreconditionerFactory<Operator>;
+
+    // Add no-oppreconditioner to factory for block size 1.
+    PrecFactory::addCreator("nothing", [](const Operator&, const pt::ptree&, const std::function<Vector()>&) {
+        return Dune::wrapPreconditioner<NothingPreconditioner<Vector>>();
+    });
+
+    auto prec = PrecFactory::create(op, prm.get_child("preconditioner"));
+    Dune::BiCGSTABSolver<Vector> solver(op, *prec, prm.get<double>("tol"), prm.get<int>("maxiter"), prm.get<int>("verbosity"));
+    Vector x(rhs.size());
+    Dune::InverseOperatorResult res;
+    solver.apply(x, rhs, res);
+    return x;
+}
+
+void test1rep(const pt::ptree& prm)
+{
+    const int bz = 1;
+    auto sol = testPrecRepeating<bz>(prm, "matr33rep.txt", "rhs3rep.txt");
+    Dune::BlockVector<Dune::FieldVector<double, bz>> expected {0.285714285714286,
+                                                               0.285714285714286,
+                                                               0.285714285714286,
+                                                               -0.214285714285714,
+                                                               -0.214285714285714,
+                                                               -0.214285714285714,
+                                                               -0.214285714285714,
+                                                               -0.214285714285714,
+                                                               -0.214285714285714};
+    BOOST_REQUIRE_EQUAL(sol.size(), expected.size());
+    for (size_t i = 0; i < sol.size(); ++i) {
+        for (int row = 0; row < bz; ++row) {
+            BOOST_CHECK_CLOSE(sol[i][row], expected[i][row], 1e-3);
+        }
+    }
+}
+
+void test3rep(const pt::ptree& prm)
+{
+    const int bz = 3;
+    auto sol = testPrecRepeating<bz>(prm, "matr33rep.txt", "rhs3rep.txt");
+    Dune::BlockVector<Dune::FieldVector<double, bz>> expected {
+        {0.285714285714286, 0.285714285714286, 0.285714285714286},
+        {-0.214285714285714, -0.214285714285714, -0.214285714285714},
+        {-0.214285714285714, -0.214285714285714, -0.214285714285714}
+    };
+    BOOST_REQUIRE_EQUAL(sol.size(), expected.size());
+    for (size_t i = 0; i < sol.size(); ++i) {
+        for (int row = 0; row < bz; ++row) {
+            BOOST_CHECK_CLOSE(sol[i][row], expected[i][row], 1e-3);
+        }
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE(TestWithRepeatingOperator)
+{
+    pt::ptree prm;
+
+    // Read parameters.
+    {
+        std::ifstream file("options_flexiblesolver_simple.json");
+        pt::read_json(file, prm);
+    }
+
+    // Test with 1x1 block solvers.
+    test1rep(prm);
+
+    // Test with 3x3 block solvers.
+    test3rep(prm);
+}
+
+
+
 #else
 
 // Do nothing if we do not have at least Dune 2.6.