From c21fd6f26d45269d96ccf4ae2ab836400fd4caf0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 29 May 2019 16:21:34 +0200 Subject: [PATCH] Add and use PreconditionerFactory class. This replaces the makePreconditoner() function. The main advantage is that it is extensible, making it easy to for example add new preconditioners to the factory at runtime. It supports both parallel and serial preconditioners. --- CMakeLists_files.cmake | 4 +- opm/simulators/linalg/FlexibleSolver.hpp | 12 +- .../linalg/OwningTwoLevelPreconditioner.hpp | 31 +- .../linalg/PreconditionerFactory.hpp | 362 ++++++++++++++++++ .../linalg/PreconditionerWithUpdate.hpp | 6 +- opm/simulators/linalg/makePreconditioner.hpp | 333 ---------------- tests/options_flexiblesolver_simple.json | 9 + tests/test_preconditionerfactory.cpp | 215 +++++++++++ 8 files changed, 614 insertions(+), 358 deletions(-) create mode 100644 opm/simulators/linalg/PreconditionerFactory.hpp delete mode 100644 opm/simulators/linalg/makePreconditioner.hpp create mode 100644 tests/options_flexiblesolver_simple.json create mode 100644 tests/test_preconditionerfactory.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index a1a54fb81..87de3f91f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -57,6 +57,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_blackoil_amg.cpp tests/test_convergencereport.cpp tests/test_flexiblesolver.cpp + tests/test_preconditionerfactory.cpp tests/test_graphcoloring.cpp tests/test_vfpproperties.cpp tests/test_milu.cpp @@ -115,6 +116,7 @@ list (APPEND TEST_DATA_FILES tests/matr33.txt tests/rhs3.txt tests/options_flexiblesolver.json + tests/options_flexiblesolver_simple.json ) @@ -176,8 +178,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/ParallelIstlInformation.hpp opm/simulators/linalg/PressureSolverPolicy.hpp opm/simulators/linalg/PressureTransferPolicy.hpp + opm/simulators/linalg/PreconditionerFactory.hpp opm/simulators/linalg/PreconditionerWithUpdate.hpp - opm/simulators/linalg/makePreconditioner.hpp opm/simulators/linalg/setupPropertyTree.hpp opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp diff --git a/opm/simulators/linalg/FlexibleSolver.hpp b/opm/simulators/linalg/FlexibleSolver.hpp index f334745da..335a265c0 100644 --- a/opm/simulators/linalg/FlexibleSolver.hpp +++ b/opm/simulators/linalg/FlexibleSolver.hpp @@ -21,7 +21,7 @@ #ifndef OPM_FLEXIBLE_SOLVER_HEADER_INCLUDED #define OPM_FLEXIBLE_SOLVER_HEADER_INCLUDED -#include +#include #include #include @@ -83,7 +83,6 @@ public: } private: - using AbstractOperatorType = Dune::AssembledLinearOperator; using AbstractPrecondType = Dune::PreconditionerWithUpdate; using AbstractScalarProductType = Dune::ScalarProduct; @@ -97,17 +96,20 @@ private: using ParOperatorType = Dune::OverlappingSchwarzOperator; auto linop = std::make_shared(matrix, comm); linearoperator_ = linop; - preconditioner_ = Dune::makePreconditioner(*linop, prm.get_child("preconditioner"), comm); + preconditioner_ + = Dune::PreconditionerFactory::create(*linop, prm.get_child("preconditioner"), comm); scalarproduct_ = Dune::createScalarProduct(comm, linearoperator_->category()); } template <> - void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, const Dune::Amg::SequentialInformation&) + void initOpPrecSp(const MatrixType& matrix, + const boost::property_tree::ptree& prm, + const Dune::Amg::SequentialInformation&) { // Sequential case. using SeqOperatorType = Dune::MatrixAdapter; auto linop = std::make_shared(matrix); linearoperator_ = linop; - preconditioner_ = Dune::makePreconditioner(*linop, prm.get_child("preconditioner")); + preconditioner_ = Dune::PreconditionerFactory::create(*linop, prm.get_child("preconditioner")); scalarproduct_ = std::make_shared>(); } diff --git a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp index b092201fe..2de80f592 100644 --- a/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp +++ b/opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp @@ -38,22 +38,22 @@ namespace Dune { -// Circular dependency between makePreconditioner() [which can make an OwningTwoLevelPreconditioner] -// and OwningTwoLevelPreconditioner [which uses makePreconditioner() to choose the fine-level smoother] +// Circular dependency between PreconditionerFactory [which can make an OwningTwoLevelPreconditioner] +// and OwningTwoLevelPreconditioner [which uses PreconditionerFactory to choose the fine-level smoother] // must be broken, accomplished by forward-declaration here. -template -std::shared_ptr> -makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm); - -template -std::shared_ptr> -makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm); +template +class PreconditionerFactory; // Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system. template class FlexibleSolver; + +/// A version of the two-level preconditioner that is: +/// - Self-contained, because it owns its policy components. +/// - Flexible, because it uses the runtime-flexible solver +/// and preconditioner factory. template ; OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm) : linear_operator_(linearoperator) - , finesmoother_(makePreconditioner(linearoperator, prm.get_child("finesmoother"))) + , finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"))) , comm_(nullptr) , weights_(Opm::Amg::getQuasiImpesWeights( linearoperator.getmat(), prm.get("pressure_var_index"), transpose)) @@ -91,8 +92,7 @@ public: OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm, const Communication& comm) : linear_operator_(linearoperator) - , finesmoother_(makePreconditioner( - linearoperator, prm.get_child("finesmoother"), comm)) + , finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"), comm)) , comm_(&comm) , weights_(Opm::Amg::getQuasiImpesWeights( linearoperator.getmat(), prm.get("pressure_var_index"), transpose)) @@ -158,13 +158,12 @@ private: using TwoLevelMethod = Dune::Amg::TwoLevelMethodCpr>; - // Handling parallel vs serial instantiation of makePreconditioner(). + // Handling parallel vs serial instantiation of preconditioner factory. template void updateImpl() { // Parallel case. - finesmoother_ = makePreconditioner( - linear_operator_, prm_.get_child("finesmoother"), *comm_); + finesmoother_ = PrecFactory::create(linear_operator_, prm_.get_child("finesmoother"), *comm_); twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_); } @@ -172,7 +171,7 @@ private: void updateImpl() { // Serial case. - finesmoother_ = makePreconditioner(linear_operator_, prm_.get_child("finesmoother")); + finesmoother_ = PrecFactory::create(linear_operator_, prm_.get_child("finesmoother")); twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_); } diff --git a/opm/simulators/linalg/PreconditionerFactory.hpp b/opm/simulators/linalg/PreconditionerFactory.hpp new file mode 100644 index 000000000..e73c217d3 --- /dev/null +++ b/opm/simulators/linalg/PreconditionerFactory.hpp @@ -0,0 +1,362 @@ + +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2019 SINTEF Digital, Mathematics and Cybernetics. + + 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_PRECONDITIONERFACTORY_HEADER +#define OPM_PRECONDITIONERFACTORY_HEADER + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +namespace Dune +{ + +/// This is an object factory for creating preconditioners. The +/// user need only interact with the factory through the static +/// methods addStandardPreconditioners() and create(). In addition +/// a user can call the addCreator() static method to add further +/// preconditioners. +template +class PreconditionerFactory +{ +public: + /// Linear algebra types. + using Matrix = typename Operator::matrix_type; + using Vector = typename Operator::domain_type; // Assuming symmetry: that domain and range types are the same. + + /// The type of pointer returned by create(). + using PrecPtr = std::shared_ptr>; + + /// The type of creator functions passed to addCreator(). + using Creator = std::function; + using ParCreator = std::function; + + /// Create a new serial preconditioner and return a pointer to it. + /// \param op operator to be preconditioned. + /// \param prm parameters for the preconditioner, in particular its type. + /// \return (smart) pointer to the created preconditioner. + static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm) + { + return instance().doCreate(op, prm); + } + + /// Create a new parallel preconditioner and return a pointer to it. + /// \param op operator to be preconditioned. + /// \param prm parameters for the preconditioner, in particular its type. + /// \param comm communication object (typically OwnerOverlapCopyCommunication). + /// \return (smart) pointer to the created preconditioner. + static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm, const Comm& comm) + { + return instance().doCreate(op, prm, comm); + } + + /// Add a creator for a serial preconditioner to the PreconditionerFactory. + /// After the call, the user may obtain a preconditioner by + /// calling create() with the given type string as a parameter + /// contained in the property_tree. + /// \param type the type string we want the PreconditionerFactory to + /// associate with the preconditioner. + /// \param creator a function or lambda creating a preconditioner. + static void addCreator(const std::string& type, Creator creator) + { + instance().doAddCreator(type, creator); + } + + /// Add a creator for a parallel preconditioner to the PreconditionerFactory. + /// After the call, the user may obtain a preconditioner by + /// calling create() with the given type string as a parameter + /// contained in the property_tree. + /// \param type the type string we want the PreconditionerFactory to + /// associate with the preconditioner. + /// \param creator a function or lambda creating a preconditioner. + static void addCreator(const std::string& type, ParCreator creator) + { + instance().doAddCreator(type, creator); + } + +private: + using CriterionBase + = Dune::Amg::AggregationCriterion>; + using Criterion = Dune::Amg::CoarsenCriterion; + + // Helpers for creation of AMG preconditioner. + static Criterion amgCriterion(const boost::property_tree::ptree& prm) + { + Criterion criterion(15, prm.get("coarsenTarget")); + criterion.setDefaultValuesIsotropic(2); + criterion.setAlpha(prm.get("alpha")); + criterion.setBeta(prm.get("beta")); + criterion.setMaxLevel(prm.get("maxlevel")); + criterion.setSkipIsolated(false); + criterion.setDebugLevel(prm.get("verbosity")); + return criterion; + } + + template + static auto amgSmootherArgs(const boost::property_tree::ptree& prm) + { + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + SmootherArgs smootherArgs; + smootherArgs.iterations = prm.get("iterations"); + // smootherArgs.overlap=SmootherArgs::vertex; + // smootherArgs.overlap=SmootherArgs::none; + // smootherArgs.overlap=SmootherArgs::aggregate; + smootherArgs.relaxationFactor = prm.get("relaxation"); + return smootherArgs; + } + + template + static PrecPtr makeAmgPreconditioner(const Operator& op, const boost::property_tree::ptree& prm) + { + auto crit = amgCriterion(prm); + auto sargs = amgSmootherArgs(prm); + return std::make_shared>(op, crit, sargs); + } + + // Add a useful default set of preconditioners to the factory. + // This is the default template, used for parallel preconditioners. + // (Serial specialization below). + template + void addStandardPreconditioners() + { + using namespace Dune; + using O = Operator; + using M = Matrix; + using V = Vector; + using P = boost::property_tree::ptree; + using C = Comm; + doAddCreator("ILU0", [](const O& op, const P& prm, const C& comm) { + const double w = prm.get("relaxation"); + return wrapBlockPreconditioner>>(comm, op.getmat(), w); + }); + doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) { + const double w = prm.get("relaxation"); + // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. + return wrapPreconditioner>( + op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU); + }); + doAddCreator("ILUn", [](const O& op, const P& prm, const C& comm) { + const int n = prm.get("ilulevel"); + const double w = prm.get("relaxation"); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + doAddCreator("GS", [](const O& op, const P& prm, const C& comm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + doAddCreator("SOR", [](const O& op, const P& prm, const C& comm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + doAddCreator("SSOR", [](const O& op, const P& prm, const C& comm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapBlockPreconditioner>>(comm, op.getmat(), n, w); + }); + doAddCreator("amg", [](const O& op, const P& prm, const C& comm) { + const std::string smoother = prm.get("smoother"); + if (smoother == "ILU0") { + using Smoother = Opm::ParallelOverlappingILU0; + auto crit = amgCriterion(prm); + auto sargs = amgSmootherArgs(prm); + return std::make_shared>(op, crit, sargs, comm); + } else { + std::string msg("No such smoother: "); + msg += smoother; + throw std::runtime_error(msg); + } + }); + doAddCreator("cpr", [](const O& op, const P& prm, const C& comm) { + return std::make_shared>(op, prm, comm); + }); + doAddCreator("cprt", [](const O& op, const P& prm, const C& comm) { + return std::make_shared>(op, prm, comm); + }); + } + + // Add a useful default set of preconditioners to the factory. + // This is the specialization for the serial case. + template <> + void addStandardPreconditioners() + { + using namespace Dune; + using O = Operator; + using M = Matrix; + using V = Vector; + using P = boost::property_tree::ptree; + doAddCreator("ILU0", [](const O& op, const P& prm) { + const double w = prm.get("relaxation"); + return wrapPreconditioner>(op.getmat(), w); + }); + doAddCreator("ParOverILU0", [](const O& op, const P& prm) { + const double w = prm.get("relaxation"); + return wrapPreconditioner>(op.getmat(), 0, w, Opm::MILU_VARIANT::ILU); + }); + doAddCreator("ILUn", [](const O& op, const P& prm) { + const int n = prm.get("ilulevel"); + const double w = prm.get("relaxation"); + return wrapPreconditioner>(op.getmat(), n, w); + }); + doAddCreator("Jac", [](const O& op, const P& prm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapPreconditioner>(op.getmat(), n, w); + }); + doAddCreator("GS", [](const O& op, const P& prm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapPreconditioner>(op.getmat(), n, w); + }); + doAddCreator("SOR", [](const O& op, const P& prm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapPreconditioner>(op.getmat(), n, w); + }); + doAddCreator("SSOR", [](const O& op, const P& prm) { + const int n = prm.get("repeats"); + const double w = prm.get("relaxation"); + return wrapPreconditioner>(op.getmat(), n, w); + }); + doAddCreator("amg", [](const O& op, const P& prm) { + const std::string smoother = prm.get("smoother"); + if (smoother == "ILU0") { + using Smoother = SeqILU0; + 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") { + using Smoother = SeqILUn; + return makeAmgPreconditioner(op, prm); + } else { + std::string msg("No such smoother: "); + msg += smoother; + throw std::runtime_error(msg); + } + }); + doAddCreator("famg", [](const O& op, const P& prm) { + 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) { + return std::make_shared>(op, prm); + }); + doAddCreator("cprt", [](const O& op, const P& prm) { + return std::make_shared>(op, prm); + }); + } + + + // The method that implements the singleton pattern, + // using the Meyers singleton technique. + static PreconditionerFactory& instance() + { + static PreconditionerFactory singleton; + return singleton; + } + + // Private constructor, to keep users from creating a PreconditionerFactory. + PreconditionerFactory() + { + addStandardPreconditioners(); + } + + // Actually creates the product object. + PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm) + { + const std::string& type = prm.get("type"); + auto it = creators_.find(type); + if (it == creators_.end()) { + std::ostringstream msg; + msg << "Preconditioner type " << type << " is not registered in the factory. Available types are: "; + for (const auto& prec : creators_) { + msg << prec.first << ' '; + } + msg << std::endl; + throw std::runtime_error(msg.str()); + } + return it->second(op, prm); + } + + PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm, const Comm& comm) + { + const std::string& type = prm.get("type"); + auto it = parallel_creators_.find(type); + if (it == parallel_creators_.end()) { + std::ostringstream msg; + msg << "Parallel preconditioner type " << type + << " is not registered in the factory. Available types are: "; + for (const auto& prec : parallel_creators_) { + msg << prec.first << ' '; + } + msg << std::endl; + throw std::runtime_error(msg.str()); + } + return it->second(op, prm, comm); + } + + // Actually adds the creator. + void doAddCreator(const std::string& type, Creator c) + { + creators_[type] = c; + } + + // Actually adds the creator. + void doAddCreator(const std::string& type, ParCreator c) + { + parallel_creators_[type] = c; + } + + // This map contains the whole factory, i.e. all the Creators. + std::map creators_; + std::map parallel_creators_; +}; + +} // namespace Dune + +#endif // OPM_PRECONDITIONERFACTORY_HEADER diff --git a/opm/simulators/linalg/PreconditionerWithUpdate.hpp b/opm/simulators/linalg/PreconditionerWithUpdate.hpp index 4f63c69ea..ba9111dfd 100644 --- a/opm/simulators/linalg/PreconditionerWithUpdate.hpp +++ b/opm/simulators/linalg/PreconditionerWithUpdate.hpp @@ -48,17 +48,17 @@ public: using X = typename OriginalPreconditioner::domain_type; using Y = typename OriginalPreconditioner::range_type; - virtual void pre (X& x, Y& b) override + virtual void pre(X& x, Y& b) override { orig_precond_.pre(x, b); } - virtual void apply (X& v, const Y& d) override + virtual void apply(X& v, const Y& d) override { orig_precond_.apply(v, d); } - virtual void post (X& x) override + virtual void post(X& x) override { orig_precond_.post(x); } diff --git a/opm/simulators/linalg/makePreconditioner.hpp b/opm/simulators/linalg/makePreconditioner.hpp deleted file mode 100644 index c3f535dd2..000000000 --- a/opm/simulators/linalg/makePreconditioner.hpp +++ /dev/null @@ -1,333 +0,0 @@ -/* - Copyright 2019 SINTEF Digital, Mathematics and Cybernetics. - - 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_MAKEPRECONDITIONER_HEADER_INCLUDED -#define OPM_MAKEPRECONDITIONER_HEADER_INCLUDED - -#include -#include -#include -#include -#include - -#include -#include -#include - -#include - -namespace Dune -{ - -template -std::shared_ptr> -makeSeqPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm) -{ - using MatrixType = typename OperatorType::matrix_type; - const MatrixType& matrix = linearoperator.getmat(); - const double w = prm.get("relaxation"); // Used by all the preconditioners. - const std::string& precond(prm.get("type")); - if (precond == "ILU0") { - return wrapPreconditioner>(matrix, w); - } else if (precond == "ParOverILU0") { - return wrapPreconditioner>( - matrix, 0, w, Opm::MILU_VARIANT::ILU); - } else if (precond == "ILUn") { - const int n = prm.get("ilulevel"); - return wrapPreconditioner>(matrix, n, w); - } - const int n = prm.get("repeats"); - if (precond == "Jac") { - return wrapPreconditioner>(matrix, n, w); - } else if (precond == "GS") { - return wrapPreconditioner>(matrix, n, w); - } else if (precond == "SOR") { - return wrapPreconditioner>(matrix, n, w); - } else if (precond == "SSOR") { - return wrapPreconditioner>(matrix, n, w); - } else { - std::string msg("No such preconditioner : "); - msg += precond; - throw std::runtime_error(msg); - } -} - -template -std::shared_ptr> -makeParPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm) -{ - using MatrixType = typename OperatorType::matrix_type; - const MatrixType& matrix = linearoperator.getmat(); - const double w = prm.get("relaxation"); // Used by all the preconditioners. - const std::string& precond(prm.get("type")); - if (precond == "ILU0") { - return wrapBlockPreconditioner>>( - comm, matrix, w); - } else if (precond == "ParOverILU0") { - // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. - return wrapPreconditioner< - DummyUpdatePreconditioner>>( - matrix, comm, 0, w, Opm::MILU_VARIANT::ILU); - } else if (precond == "ParOverILUn") { - // Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner. - const int n = prm.get("ilulevel"); - return wrapPreconditioner< - DummyUpdatePreconditioner>>( - matrix, comm, n, w, Opm::MILU_VARIANT::ILU); - } else if (precond == "ILUn") { - const int n = prm.get("ilulevel"); - return wrapBlockPreconditioner>>( - comm, matrix, n, w); - } - const int n = prm.get("repeats"); - if (precond == "Jac") { - return wrapBlockPreconditioner>>( - comm, matrix, n, w); - } else if (precond == "GS") { - return wrapBlockPreconditioner>>( - comm, matrix, n, w); - } else if (precond == "SOR") { - return wrapBlockPreconditioner>>( - comm, matrix, n, w); - } else if (precond == "SSOR") { - return wrapBlockPreconditioner>>( - comm, matrix, n, w); - } else { - std::string msg("No such preconditioner : "); - msg += precond; - throw std::runtime_error(msg); - } -} - -template -std::shared_ptr> -makeAmgPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm) -{ - using MatrixType = typename OperatorType::matrix_type; - using CriterionBase - = Dune::Amg::AggregationCriterion>; - using Criterion = Dune::Amg::CoarsenCriterion; - int coarsenTarget = prm.get("coarsenTarget"); - int ml = prm.get("maxlevel"); - Criterion criterion(15, coarsenTarget); - criterion.setDefaultValuesIsotropic(2); - criterion.setAlpha(prm.get("alpha")); - criterion.setBeta(prm.get("beta")); - criterion.setMaxLevel(ml); - criterion.setSkipIsolated(false); - criterion.setDebugLevel(prm.get("verbosity")); - if (prm.get("type") == "famg") { - Dune::Amg::Parameters parms; - parms.setNoPreSmoothSteps(1); - parms.setNoPostSmoothSteps(1); - return wrapPreconditioner>(linearoperator, criterion, parms); - } else { - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - SmootherArgs smootherArgs; - smootherArgs.iterations = prm.get("iterations"); - // smootherArgs.overlap=SmootherArgs::vertex; - // smootherArgs.overlap=SmootherArgs::none; - // smootherArgs.overlap=SmootherArgs::aggregate; - smootherArgs.relaxationFactor = prm.get("relaxation"); - return std::make_shared>( - linearoperator, criterion, smootherArgs); - } - return std::shared_ptr>(nullptr); -} - -template -std::shared_ptr> -makeParAmgPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm) -{ - using MatrixType = typename OperatorType::matrix_type; - using CriterionBase - = Dune::Amg::AggregationCriterion>; - using Criterion = Dune::Amg::CoarsenCriterion; - int coarsenTarget = prm.get("coarsenTarget"); - int ml = prm.get("maxlevel"); - Criterion criterion(15, coarsenTarget); - criterion.setDefaultValuesIsotropic(2); - criterion.setAlpha(prm.get("alpha")); - criterion.setBeta(prm.get("beta")); - criterion.setMaxLevel(ml); - criterion.setSkipIsolated(false); - criterion.setDebugLevel(prm.get("verbosity")); - if (prm.get("type") == "famg") { - throw std::runtime_error("The FastAMG preconditioner cannot be used in parallel."); - } else { - typedef typename Dune::Amg::SmootherTraits::Arguments SmootherArgs; - SmootherArgs smootherArgs; - smootherArgs.iterations = prm.get("iterations"); - // smootherArgs.overlap=SmootherArgs::vertex; - // smootherArgs.overlap=SmootherArgs::none; - // smootherArgs.overlap=SmootherArgs::aggregate; - smootherArgs.relaxationFactor = prm.get("relaxation"); - return std::make_shared>( - linearoperator, criterion, smootherArgs, comm); - } -} - - - - - -template -std::shared_ptr> -makeAmgPreconditioners(const OperatorType& linearoperator, const boost::property_tree::ptree& prm) -{ - using MatrixType = typename OperatorType::matrix_type; - if (prm.get("type") == "famg") { - // smoother type should not be used - return makeAmgPreconditioner, OperatorType, VectorType>( - linearoperator, prm); - } - - std::string precond = prm.get("smoother"); - if (precond == "ILU0") { - return makeAmgPreconditioner, OperatorType, VectorType>( - linearoperator, prm); - } else if (precond == "Jac") { - return makeAmgPreconditioner, OperatorType, VectorType>( - linearoperator, prm); - // } else if (precond == "GS") { - // return makeAmgPreconditioner< - // Dune::SeqGS, MatrixType, VectorType>(linearoperator, prm); - } else if (precond == "SOR") { - return makeAmgPreconditioner, OperatorType, VectorType>( - linearoperator, prm); - } else if (precond == "SSOR") { - return makeAmgPreconditioner, OperatorType, VectorType>( - linearoperator, prm); - } else if (precond == "ILUn") { - return makeAmgPreconditioner, OperatorType, VectorType>( - linearoperator, prm); - } else { - std::string msg("No such sequential preconditioner : "); - msg += precond; - throw std::runtime_error(msg); - } -} - - -template -std::shared_ptr> -makeParAmgPreconditioners(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm) -{ - if (prm.get("type") == "famg") { - throw std::runtime_error("The FastAMG preconditioner cannot be used in parallel."); - } - - using MatrixType = typename OperatorType::matrix_type; - std::string precond = prm.get("smoother"); - if (precond == "ILU0") { - using SmootherType = Opm::ParallelOverlappingILU0; - return makeParAmgPreconditioner(linearoperator, prm, comm); - /* - return makeParAmgPreconditioner, MatrixType, VectorType>( - linearoperator, prm, comm); - } else if (precond == "Jac") { - return makeParAmgPreconditioner, MatrixType, VectorType>( - linearoperator, prm, comm); - // } else if (precond == "GS") { - // return makeParAmgPreconditioner< - // Dune::SeqGS, MatrixType, VectorType>(linearoperator, prm, comm); - } else if (precond == "SOR") { - return makeParAmgPreconditioner, MatrixType, VectorType>( - linearoperator, prm, comm); - } else if (precond == "SSOR") { - return makeParAmgPreconditioner, MatrixType, VectorType>( - linearoperator, prm, comm); - } else if (precond == "ILUn") { - return makeParAmgPreconditioner, MatrixType, VectorType>( - linearoperator, prm, comm); - */ - } else { - std::string msg("No such parallel preconditioner : "); - msg += precond; - throw std::runtime_error(msg); - } -} - - - -template -std::shared_ptr> -makeTwoLevelPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm) -{ - if (prm.get("type") == "cpr") { - return std::make_shared>(linearoperator, prm); - } else if (prm.get("type") == "cprt") { - return std::make_shared>(linearoperator, prm); - } else { - std::string msg("Wrong cpr Should not happen"); - throw std::runtime_error(msg); - } -} - -template -std::shared_ptr> -makeParTwoLevelPreconditioner(const OperatorType& linearoperator, - const boost::property_tree::ptree& prm, - const Comm& comm) -{ - if (prm.get("type") == "cpr") { - return std::make_shared>( - linearoperator, prm, comm); - } else if (prm.get("type") == "cprt") { - return std::make_shared>( - linearoperator, prm, comm); - } else { - std::string msg("Wrong cpr Should not happen"); - throw std::runtime_error(msg); - } -} - -template -std::shared_ptr> -makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm) -{ - if ((prm.get("type") == "famg") or (prm.get("type") == "amg")) { - return makeAmgPreconditioners(linearoperator, prm); - } else if ((prm.get("type") == "cpr") - or (prm.get("type") == "cprt")) { - return makeTwoLevelPreconditioner(linearoperator, prm); - } else { - return makeSeqPreconditioner(linearoperator, prm); - } -} - -template -std::shared_ptr> -makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm) -{ - if ((prm.get("type") == "famg") or (prm.get("type") == "amg")) { - return makeParAmgPreconditioners(linearoperator, prm, comm); - } else if ((prm.get("type") == "cpr") - or (prm.get("type") == "cprt")) { - return makeParTwoLevelPreconditioner(linearoperator, prm, comm); - } else { - return makeParPreconditioner(linearoperator, prm, comm); - } -} - - -} // namespace Dune - - -#endif // OPM_MAKEPRECONDITIONER_HEADER_INCLUDED diff --git a/tests/options_flexiblesolver_simple.json b/tests/options_flexiblesolver_simple.json new file mode 100644 index 000000000..059007d9f --- /dev/null +++ b/tests/options_flexiblesolver_simple.json @@ -0,0 +1,9 @@ +{ + "tol": "1e-12", + "maxiter": "200", + "verbosity": "0", + "solver": "bicgstab", + "preconditioner": { + "type": "nothing" + } +} diff --git a/tests/test_preconditionerfactory.cpp b/tests/test_preconditionerfactory.cpp new file mode 100644 index 000000000..5545dc328 --- /dev/null +++ b/tests/test_preconditionerfactory.cpp @@ -0,0 +1,215 @@ +/* + Copyright 2019 SINTEF Digital, Mathematics and Cybernetics. + + 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_PreconditionerFactory +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + + +template +class NothingPreconditioner : public Dune::Preconditioner +{ +public: + virtual void pre(X&, X&) override + { + } + + virtual void apply(X& v, const X& d) override + { + v = d; + } + + virtual void post(X&) override + { + } + + virtual Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::sequential; + } +}; + + +template +Dune::BlockVector> +testPrec(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 = Dune::MatrixAdapter; + Operator op(matrix); + using PrecFactory = Dune::PreconditionerFactory; + 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; +} + +namespace pt = boost::property_tree; + +void test1(const pt::ptree& prm) +{ + const int bz = 1; + auto sol = testPrec(prm, "matr33.txt", "rhs3.txt"); + Dune::BlockVector> expected {-1.62493, + -1.76435e-06, + 1.86991e-10, + -458.542, + 2.28308e-06, + -2.45341e-07, + -1.48005, + -5.02264e-07, + -1.049e-05}; + 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 test3(const pt::ptree& prm) +{ + const int bz = 3; + auto sol = testPrec(prm, "matr33.txt", "rhs3.txt"); + Dune::BlockVector> expected {{-1.62493, -1.76435e-06, 1.86991e-10}, + {-458.542, 2.28308e-06, -2.45341e-07}, + {-1.48005, -5.02264e-07, -1.049e-05}}; + 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(TestDefaultPreconditionerFactory) +{ + pt::ptree prm; + + // Read parameters. + { + std::ifstream file("options_flexiblesolver.json"); + pt::read_json(file, prm); + } + + // Test with 1x1 block solvers. + test1(prm); + + // Test with 3x3 block solvers. + test3(prm); +} + + +template +using M = Dune::BCRSMatrix>; +template +using V = Dune::BlockVector>; +template +using O = Dune::MatrixAdapter, V, V>; +template +using PF = Dune::PreconditionerFactory>; + + +BOOST_AUTO_TEST_CASE(TestAddingPreconditioner) +{ + namespace pt = boost::property_tree; + pt::ptree prm; + + // Read parameters. + { + std::ifstream file("options_flexiblesolver_simple.json"); // Requests "nothing" for preconditioner type. + pt::read_json(file, prm); + } + + // Test with 1x1 block solvers. + { + const int bz = 1; + BOOST_CHECK_THROW(testPrec(prm, "matr33.txt", "rhs3.txt"), std::runtime_error); + } + + // Test with 3x3 block solvers. + { + const int bz = 3; + BOOST_CHECK_THROW(testPrec(prm, "matr33.txt", "rhs3.txt"), std::runtime_error); + } + + + // Add preconditioner to factory for block size 1. + PF<1>::addCreator("nothing", [](const O<1>&, const pt::ptree&) { + return Dune::wrapPreconditioner>>(); + }); + + + // Test with 1x1 block solvers. + test1(prm); + + // Test with 3x3 block solvers. + { + const int bz = 3; + BOOST_CHECK_THROW(testPrec(prm, "matr33.txt", "rhs3.txt"), std::runtime_error); + } + + // Add preconditioner to factory for block size 3. + PF<3>::addCreator("nothing", [](const O<3>&, const pt::ptree&) { + return Dune::wrapPreconditioner>>(); + }); + + // Test with 1x1 block solvers. + test1(prm); + + // Test with 3x3 block solvers. + test3(prm); +}