From d910d42f6e5596b1bf63c78d371f7ec721b47d42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 17 Dec 2020 14:51:05 +0100 Subject: [PATCH] Only add AMG preconditioners to factory if sensible. Also add test using a new operator class, that would not compile without the change. --- CMakeLists_files.cmake | 2 + .../linalg/PreconditionerFactory.hpp | 147 +++++++++-------- tests/matr33rep.txt | 84 ++++++++++ tests/rhs3rep.txt | 12 ++ tests/test_preconditionerfactory.cpp | 151 ++++++++++++++++++ 5 files changed, 329 insertions(+), 67 deletions(-) create mode 100644 tests/matr33rep.txt create mode 100644 tests/rhs3rep.txt diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 28c7ca555..70cc1c3df 100644 --- a/CMakeLists_files.cmake +++ b/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 ) diff --git a/opm/simulators/linalg/PreconditionerFactory.hpp b/opm/simulators/linalg/PreconditionerFactory.hpp index fc1647f6a..89fbf9ca5 100644 --- a/opm/simulators/linalg/PreconditionerFactory.hpp +++ b/opm/simulators/linalg/PreconditionerFactory.hpp @@ -288,17 +288,25 @@ private: const double w = prm.get("relaxation", 1.0); return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); }); - doAddCreator("amg", [](const O& op, const P& prm, const std::function&, const C& comm) { - const std::string smoother = prm.get("smoother", "ParOverILU0"); - if (smoother == "ILU0" || smoother == "ParOverILU0") { - using Smoother = Opm::ParallelOverlappingILU0; - auto crit = amgCriterion(prm); - auto sargs = amgSmootherArgs(prm); - return std::make_shared>(op, crit, sargs, comm); - } else { - OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<"."); - } - }); + + // Only add AMG preconditioners to the factory if the operator + // is the overlapping schwarz operator. This could be extended + // later, but at this point no other operators are compatible + // with the AMG hierarchy construction. + if constexpr (std::is_same_v>) { + doAddCreator("amg", [](const O& op, const P& prm, const std::function&, const C& comm) { + const std::string smoother = prm.get("smoother", "ParOverILU0"); + if (smoother == "ILU0" || smoother == "ParOverILU0") { + using Smoother = Opm::ParallelOverlappingILU0; + auto crit = amgCriterion(prm); + auto sargs = amgSmootherArgs(prm); + return std::make_shared>(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 weightsCalculator, const C& comm) { assert(weightsCalculator); return std::make_shared>(op, prm, weightsCalculator, comm); @@ -355,74 +363,79 @@ private: const double w = prm.get("relaxation", 1.0); return wrapPreconditioner>(op.getmat(), n, w); }); - doAddCreator("amg", [](const O& op, const P& prm, const std::function&) { - const std::string smoother = prm.get("smoother", "ParOverILU0"); - if (smoother == "ILU0" || smoother == "ParOverILU0") { + + // Only add AMG preconditioners to the factory if the operator + // is an actual matrix operator. + if constexpr (std::is_same_v>) { + doAddCreator("amg", [](const O& op, const P& prm, const std::function&) { + const std::string smoother = prm.get("smoother", "ParOverILU0"); + if (smoother == "ILU0" || smoother == "ParOverILU0") { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - using Smoother = SeqILU; + using Smoother = SeqILU; #else - using Smoother = SeqILU0; + using Smoother = SeqILU0; #endif - return makeAmgPreconditioner(op, prm); - } else if (smoother == "Jac") { - using Smoother = SeqJac; - return makeAmgPreconditioner(op, prm); - } else if (smoother == "SOR") { - using Smoother = SeqSOR; - return makeAmgPreconditioner(op, prm); - } else if (smoother == "SSOR") { - using Smoother = SeqSSOR; - return makeAmgPreconditioner(op, prm); - } else if (smoother == "ILUn") { + return makeAmgPreconditioner(op, prm); + } else if (smoother == "Jac") { + using Smoother = SeqJac; + return makeAmgPreconditioner(op, prm); + } else if (smoother == "SOR") { + using Smoother = SeqSOR; + return makeAmgPreconditioner(op, prm); + } else if (smoother == "SSOR") { + using Smoother = SeqSSOR; + return makeAmgPreconditioner(op, prm); + } else if (smoother == "ILUn") { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - using Smoother = SeqILU; + using Smoother = SeqILU; #else - using Smoother = SeqILUn; + using Smoother = SeqILUn; #endif - return makeAmgPreconditioner(op, prm); - } else { - OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<"."); - } - }); - doAddCreator("kamg", [](const O& op, const P& prm, const std::function&) { - const std::string smoother = prm.get("smoother", "ParOverILU0"); - if (smoother == "ILU0" || smoother == "ParOverILU0") { + return makeAmgPreconditioner(op, prm); + } else { + OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << "."); + } + }); + doAddCreator("kamg", [](const O& op, const P& prm, const std::function&) { + const std::string smoother = prm.get("smoother", "ParOverILU0"); + if (smoother == "ILU0" || smoother == "ParOverILU0") { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - using Smoother = SeqILU; + using Smoother = SeqILU; #else - using Smoother = SeqILU0; + using Smoother = SeqILU0; #endif - return makeAmgPreconditioner(op, prm, true); - } else if (smoother == "Jac") { - using Smoother = SeqJac; - return makeAmgPreconditioner(op, prm, true); - } else if (smoother == "SOR") { - using Smoother = SeqSOR; - return makeAmgPreconditioner(op, prm, true); - // } else if (smoother == "GS") { - // using Smoother = SeqGS; - // return makeAmgPreconditioner(op, prm, true); - } else if (smoother == "SSOR") { - using Smoother = SeqSSOR; - return makeAmgPreconditioner(op, prm, true); - } else if (smoother == "ILUn") { + return makeAmgPreconditioner(op, prm, true); + } else if (smoother == "Jac") { + using Smoother = SeqJac; + return makeAmgPreconditioner(op, prm, true); + } else if (smoother == "SOR") { + using Smoother = SeqSOR; + return makeAmgPreconditioner(op, prm, true); + // } else if (smoother == "GS") { + // using Smoother = SeqGS; + // return makeAmgPreconditioner(op, prm, true); + } else if (smoother == "SSOR") { + using Smoother = SeqSSOR; + return makeAmgPreconditioner(op, prm, true); + } else if (smoother == "ILUn") { #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7) - using Smoother = SeqILU; + using Smoother = SeqILU; #else - using Smoother = SeqILUn; + using Smoother = SeqILUn; #endif - return makeAmgPreconditioner(op, prm, true); - } else { - OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother <<"."); - } - }); - doAddCreator("famg", [](const O& op, const P& prm, const std::function&) { - auto crit = amgCriterion(prm); - Dune::Amg::Parameters parms; - parms.setNoPreSmoothSteps(1); - parms.setNoPostSmoothSteps(1); - return wrapPreconditioner>(op, crit, parms); - }); + return makeAmgPreconditioner(op, prm, true); + } else { + OPM_THROW(std::invalid_argument, "Properties: No smoother with name " << smoother << "."); + } + }); + doAddCreator("famg", [](const O& op, const P& prm, const std::function&) { + auto crit = amgCriterion(prm); + Dune::Amg::Parameters parms; + parms.setNoPreSmoothSteps(1); + parms.setNoPostSmoothSteps(1); + return wrapPreconditioner>(op, crit, parms); + }); + } doAddCreator("cpr", [](const O& op, const P& prm, const std::function& weightsCalculator) { return std::make_shared>(op, prm, weightsCalculator); }); diff --git a/tests/matr33rep.txt b/tests/matr33rep.txt new file mode 100644 index 000000000..1bffcb7c2 --- /dev/null +++ b/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 diff --git a/tests/rhs3rep.txt b/tests/rhs3rep.txt new file mode 100644 index 000000000..02b62f3e4 --- /dev/null +++ b/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 diff --git a/tests/test_preconditionerfactory.cpp b/tests/test_preconditionerfactory.cpp index ebd19b878..c6c0e18e1 100644 --- a/tests/test_preconditionerfactory.cpp +++ b/tests/test_preconditionerfactory.cpp @@ -232,6 +232,157 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner) test3(prm); } + +template +class RepeatingOperator : public Dune::AssembledLinearOperator +{ +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 +Dune::BlockVector> +testPrecRepeating(const boost::property_tree::ptree& prm, const std::string& matrix_filename, const std::string& rhs_filename) +{ + using Matrix = Dune::BCRSMatrix>; + using Vector = Dune::BlockVector>; + 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; + Operator op(matrix, 2); + using PrecFactory = Opm::PreconditionerFactory; + + // Add no-oppreconditioner to factory for block size 1. + PrecFactory::addCreator("nothing", [](const Operator&, const pt::ptree&, const std::function&) { + return Dune::wrapPreconditioner>(); + }); + + auto prec = PrecFactory::create(op, prm.get_child("preconditioner")); + Dune::BiCGSTABSolver solver(op, *prec, prm.get("tol"), prm.get("maxiter"), prm.get("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(prm, "matr33rep.txt", "rhs3rep.txt"); + Dune::BlockVector> 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(prm, "matr33rep.txt", "rhs3rep.txt"); + Dune::BlockVector> 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.