mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-21 08:54:08 -06:00
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.
This commit is contained in:
Atgeirr Flø Rasmussen
parent
d44b974dc8
commit
c21fd6f26d
@ -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
|
||||
|
||||
@ -21,7 +21,7 @@
|
||||
#ifndef OPM_FLEXIBLE_SOLVER_HEADER_INCLUDED
|
||||
#define OPM_FLEXIBLE_SOLVER_HEADER_INCLUDED
|
||||
|
||||
#include <opm/simulators/linalg/makePreconditioner.hpp>
|
||||
#include <opm/simulators/linalg/PreconditionerFactory.hpp>
|
||||
|
||||
#include <dune/common/fmatrix.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
@ -83,7 +83,6 @@ public:
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
using AbstractOperatorType = Dune::AssembledLinearOperator<MatrixType, VectorType, VectorType>;
|
||||
using AbstractPrecondType = Dune::PreconditionerWithUpdate<VectorType, VectorType>;
|
||||
using AbstractScalarProductType = Dune::ScalarProduct<VectorType>;
|
||||
@ -97,17 +96,20 @@ private:
|
||||
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
|
||||
auto linop = std::make_shared<ParOperatorType>(matrix, comm);
|
||||
linearoperator_ = linop;
|
||||
preconditioner_ = Dune::makePreconditioner<ParOperatorType, VectorType, Comm>(*linop, prm.get_child("preconditioner"), comm);
|
||||
preconditioner_
|
||||
= Dune::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, prm.get_child("preconditioner"), comm);
|
||||
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, linearoperator_->category());
|
||||
}
|
||||
template <>
|
||||
void initOpPrecSp<Dune::Amg::SequentialInformation>(const MatrixType& matrix, const boost::property_tree::ptree& prm, const Dune::Amg::SequentialInformation&)
|
||||
void initOpPrecSp<Dune::Amg::SequentialInformation>(const MatrixType& matrix,
|
||||
const boost::property_tree::ptree& prm,
|
||||
const Dune::Amg::SequentialInformation&)
|
||||
{
|
||||
// Sequential case.
|
||||
using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
|
||||
auto linop = std::make_shared<SeqOperatorType>(matrix);
|
||||
linearoperator_ = linop;
|
||||
preconditioner_ = Dune::makePreconditioner<SeqOperatorType, VectorType>(*linop, prm.get_child("preconditioner"));
|
||||
preconditioner_ = Dune::PreconditionerFactory<SeqOperatorType>::create(*linop, prm.get_child("preconditioner"));
|
||||
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
|
||||
}
|
||||
|
||||
|
||||
@ -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 <class OperatorType, class VectorType>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm);
|
||||
|
||||
template <class OperatorType, class VectorType, class Comm>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm);
|
||||
template <class Operator, class Comm = Dune::Amg::SequentialInformation>
|
||||
class PreconditionerFactory;
|
||||
|
||||
|
||||
// Must forward-declare FlexibleSolver as we want to use it as solver for the pressure system.
|
||||
template <class MatrixTypeT, class VectorTypeT>
|
||||
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 <class OperatorType,
|
||||
class VectorType,
|
||||
bool transpose = false,
|
||||
@ -63,10 +63,11 @@ class OwningTwoLevelPreconditioner : public Dune::PreconditionerWithUpdate<Vecto
|
||||
public:
|
||||
using pt = boost::property_tree::ptree;
|
||||
using MatrixType = typename OperatorType::matrix_type;
|
||||
using PrecFactory = PreconditionerFactory<OperatorType, Communication>;
|
||||
|
||||
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm)
|
||||
: linear_operator_(linearoperator)
|
||||
, finesmoother_(makePreconditioner<OperatorType, VectorType>(linearoperator, prm.get_child("finesmoother")))
|
||||
, finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother")))
|
||||
, comm_(nullptr)
|
||||
, weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
|
||||
linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
|
||||
@ -91,8 +92,7 @@ public:
|
||||
|
||||
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm, const Communication& comm)
|
||||
: linear_operator_(linearoperator)
|
||||
, finesmoother_(makePreconditioner<OperatorType, VectorType, Communication>(
|
||||
linearoperator, prm.get_child("finesmoother"), comm))
|
||||
, finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"), comm))
|
||||
, comm_(&comm)
|
||||
, weights_(Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
|
||||
linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
|
||||
@ -158,13 +158,12 @@ private:
|
||||
using TwoLevelMethod
|
||||
= Dune::Amg::TwoLevelMethodCpr<OperatorType, CoarseSolverPolicy, Dune::Preconditioner<VectorType, VectorType>>;
|
||||
|
||||
// Handling parallel vs serial instantiation of makePreconditioner().
|
||||
// Handling parallel vs serial instantiation of preconditioner factory.
|
||||
template <class Comm>
|
||||
void updateImpl()
|
||||
{
|
||||
// Parallel case.
|
||||
finesmoother_ = makePreconditioner<OperatorType, VectorType, Communication>(
|
||||
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<Dune::Amg::SequentialInformation>()
|
||||
{
|
||||
// Serial case.
|
||||
finesmoother_ = makePreconditioner<OperatorType, VectorType>(linear_operator_, prm_.get_child("finesmoother"));
|
||||
finesmoother_ = PrecFactory::create(linear_operator_, prm_.get_child("finesmoother"));
|
||||
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
|
||||
}
|
||||
|
||||
|
||||
362
opm/simulators/linalg/PreconditionerFactory.hpp
Normal file
362
opm/simulators/linalg/PreconditionerFactory.hpp
Normal file
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_PRECONDITIONERFACTORY_HEADER
|
||||
#define OPM_PRECONDITIONERFACTORY_HEADER
|
||||
|
||||
#include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
|
||||
#include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
|
||||
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
|
||||
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
|
||||
#include <opm/simulators/linalg/amgcpr.hh>
|
||||
|
||||
#include <dune/istl/paamg/amg.hh>
|
||||
#include <dune/istl/paamg/fastamg.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
|
||||
#include <boost/property_tree/ptree.hpp>
|
||||
|
||||
#include <map>
|
||||
#include <memory>
|
||||
|
||||
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 Operator, class Comm>
|
||||
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<Dune::PreconditionerWithUpdate<Vector, Vector>>;
|
||||
|
||||
/// The type of creator functions passed to addCreator().
|
||||
using Creator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&)>;
|
||||
using ParCreator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&, const Comm&)>;
|
||||
|
||||
/// 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<Dune::Amg::SymmetricMatrixDependency<Matrix, Dune::Amg::FirstDiagonal>>;
|
||||
using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
|
||||
|
||||
// Helpers for creation of AMG preconditioner.
|
||||
static Criterion amgCriterion(const boost::property_tree::ptree& prm)
|
||||
{
|
||||
Criterion criterion(15, prm.get<int>("coarsenTarget"));
|
||||
criterion.setDefaultValuesIsotropic(2);
|
||||
criterion.setAlpha(prm.get<double>("alpha"));
|
||||
criterion.setBeta(prm.get<double>("beta"));
|
||||
criterion.setMaxLevel(prm.get<int>("maxlevel"));
|
||||
criterion.setSkipIsolated(false);
|
||||
criterion.setDebugLevel(prm.get<int>("verbosity"));
|
||||
return criterion;
|
||||
}
|
||||
|
||||
template <typename Smoother>
|
||||
static auto amgSmootherArgs(const boost::property_tree::ptree& prm)
|
||||
{
|
||||
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
|
||||
SmootherArgs smootherArgs;
|
||||
smootherArgs.iterations = prm.get<int>("iterations");
|
||||
// smootherArgs.overlap=SmootherArgs::vertex;
|
||||
// smootherArgs.overlap=SmootherArgs::none;
|
||||
// smootherArgs.overlap=SmootherArgs::aggregate;
|
||||
smootherArgs.relaxationFactor = prm.get<double>("relaxation");
|
||||
return smootherArgs;
|
||||
}
|
||||
|
||||
template <class Smoother>
|
||||
static PrecPtr makeAmgPreconditioner(const Operator& op, const boost::property_tree::ptree& prm)
|
||||
{
|
||||
auto crit = amgCriterion(prm);
|
||||
auto sargs = amgSmootherArgs<Smoother>(prm);
|
||||
return std::make_shared<Dune::Amg::AMGCPR<Operator, Vector, Smoother>>(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 <class CommArg>
|
||||
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<double>("relaxation");
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILU0<M, V, V>>>(comm, op.getmat(), w);
|
||||
});
|
||||
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) {
|
||||
const double w = prm.get<double>("relaxation");
|
||||
// Already a parallel preconditioner. Need to pass comm, but no need to wrap it in a BlockPreconditioner.
|
||||
return wrapPreconditioner<Opm::ParallelOverlappingILU0<M, V, V, C>>(
|
||||
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<int>("ilulevel");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqILUn<M, V, V>>>(comm, op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqJac<M, V, V>>>(comm, op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("GS", [](const O& op, const P& prm, const C& comm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqGS<M, V, V>>>(comm, op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("SOR", [](const O& op, const P& prm, const C& comm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSOR<M, V, V>>>(comm, op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("SSOR", [](const O& op, const P& prm, const C& comm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(comm, op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("amg", [](const O& op, const P& prm, const C& comm) {
|
||||
const std::string smoother = prm.get<std::string>("smoother");
|
||||
if (smoother == "ILU0") {
|
||||
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 {
|
||||
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<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, comm);
|
||||
});
|
||||
doAddCreator("cprt", [](const O& op, const P& prm, const C& comm) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, comm);
|
||||
});
|
||||
}
|
||||
|
||||
// Add a useful default set of preconditioners to the factory.
|
||||
// This is the specialization for the serial case.
|
||||
template <>
|
||||
void addStandardPreconditioners<Dune::Amg::SequentialInformation>()
|
||||
{
|
||||
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<double>("relaxation");
|
||||
return wrapPreconditioner<SeqILU0<M, V, V>>(op.getmat(), w);
|
||||
});
|
||||
doAddCreator("ParOverILU0", [](const O& op, const P& prm) {
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapPreconditioner<Opm::ParallelOverlappingILU0<M, V, V>>(op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
|
||||
});
|
||||
doAddCreator("ILUn", [](const O& op, const P& prm) {
|
||||
const int n = prm.get<int>("ilulevel");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapPreconditioner<SeqILUn<M, V, V>>(op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("Jac", [](const O& op, const P& prm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapPreconditioner<SeqJac<M, V, V>>(op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("GS", [](const O& op, const P& prm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapPreconditioner<SeqGS<M, V, V>>(op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("SOR", [](const O& op, const P& prm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapPreconditioner<SeqSOR<M, V, V>>(op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("SSOR", [](const O& op, const P& prm) {
|
||||
const int n = prm.get<int>("repeats");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return wrapPreconditioner<SeqSSOR<M, V, V>>(op.getmat(), n, w);
|
||||
});
|
||||
doAddCreator("amg", [](const O& op, const P& prm) {
|
||||
const std::string smoother = prm.get<std::string>("smoother");
|
||||
if (smoother == "ILU0") {
|
||||
using Smoother = SeqILU0<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "Jac") {
|
||||
using Smoother = SeqJac<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "SOR") {
|
||||
using Smoother = SeqSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "SSOR") {
|
||||
using Smoother = SeqSSOR<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(op, prm);
|
||||
} else if (smoother == "ILUn") {
|
||||
using Smoother = SeqILUn<M, V, V>;
|
||||
return makeAmgPreconditioner<Smoother>(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<Dune::Amg::FastAMG<O, V>>(op, crit, parms);
|
||||
});
|
||||
doAddCreator("cpr", [](const O& op, const P& prm) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm);
|
||||
});
|
||||
doAddCreator("cprt", [](const O& op, const P& prm) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(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<Comm>();
|
||||
}
|
||||
|
||||
// Actually creates the product object.
|
||||
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm)
|
||||
{
|
||||
const std::string& type = prm.get<std::string>("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<std::string>("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<std::string, Creator> creators_;
|
||||
std::map<std::string, ParCreator> parallel_creators_;
|
||||
};
|
||||
|
||||
} // namespace Dune
|
||||
|
||||
#endif // OPM_PRECONDITIONERFACTORY_HEADER
|
||||
@ -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);
|
||||
}
|
||||
|
||||
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_MAKEPRECONDITIONER_HEADER_INCLUDED
|
||||
#define OPM_MAKEPRECONDITIONER_HEADER_INCLUDED
|
||||
|
||||
#include <opm/simulators/linalg/OwningBlockPreconditioner.hpp>
|
||||
#include <opm/simulators/linalg/OwningTwoLevelPreconditioner.hpp>
|
||||
#include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
|
||||
#include <opm/simulators/linalg/PreconditionerWithUpdate.hpp>
|
||||
#include <opm/simulators/linalg/amgcpr.hh>
|
||||
|
||||
#include <dune/istl/paamg/amg.hh>
|
||||
#include <dune/istl/paamg/fastamg.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
|
||||
#include <boost/property_tree/ptree.hpp>
|
||||
|
||||
namespace Dune
|
||||
{
|
||||
|
||||
template <class OperatorType, class VectorType>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
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<double>("relaxation"); // Used by all the preconditioners.
|
||||
const std::string& precond(prm.get<std::string>("type"));
|
||||
if (precond == "ILU0") {
|
||||
return wrapPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>>(matrix, w);
|
||||
} else if (precond == "ParOverILU0") {
|
||||
return wrapPreconditioner<Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType>>(
|
||||
matrix, 0, w, Opm::MILU_VARIANT::ILU);
|
||||
} else if (precond == "ILUn") {
|
||||
const int n = prm.get<int>("ilulevel");
|
||||
return wrapPreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>>(matrix, n, w);
|
||||
}
|
||||
const int n = prm.get<int>("repeats");
|
||||
if (precond == "Jac") {
|
||||
return wrapPreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>>(matrix, n, w);
|
||||
} else if (precond == "GS") {
|
||||
return wrapPreconditioner<Dune::SeqGS<MatrixType, VectorType, VectorType>>(matrix, n, w);
|
||||
} else if (precond == "SOR") {
|
||||
return wrapPreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>>(matrix, n, w);
|
||||
} else if (precond == "SSOR") {
|
||||
return wrapPreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>>(matrix, n, w);
|
||||
} else {
|
||||
std::string msg("No such preconditioner : ");
|
||||
msg += precond;
|
||||
throw std::runtime_error(msg);
|
||||
}
|
||||
}
|
||||
|
||||
template <class OperatorType, class VectorType, class Comm>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
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<double>("relaxation"); // Used by all the preconditioners.
|
||||
const std::string& precond(prm.get<std::string>("type"));
|
||||
if (precond == "ILU0") {
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>>>(
|
||||
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<Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType, Comm>>>(
|
||||
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<int>("ilulevel");
|
||||
return wrapPreconditioner<
|
||||
DummyUpdatePreconditioner<Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType, Comm>>>(
|
||||
matrix, comm, n, w, Opm::MILU_VARIANT::ILU);
|
||||
} else if (precond == "ILUn") {
|
||||
const int n = prm.get<int>("ilulevel");
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>>>(
|
||||
comm, matrix, n, w);
|
||||
}
|
||||
const int n = prm.get<int>("repeats");
|
||||
if (precond == "Jac") {
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>>>(
|
||||
comm, matrix, n, w);
|
||||
} else if (precond == "GS") {
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqGS<MatrixType, VectorType, VectorType>>>(
|
||||
comm, matrix, n, w);
|
||||
} else if (precond == "SOR") {
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>>>(
|
||||
comm, matrix, n, w);
|
||||
} else if (precond == "SSOR") {
|
||||
return wrapBlockPreconditioner<DummyUpdatePreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>>>(
|
||||
comm, matrix, n, w);
|
||||
} else {
|
||||
std::string msg("No such preconditioner : ");
|
||||
msg += precond;
|
||||
throw std::runtime_error(msg);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Smoother, class OperatorType, class VectorType>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makeAmgPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm)
|
||||
{
|
||||
using MatrixType = typename OperatorType::matrix_type;
|
||||
using CriterionBase
|
||||
= Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<MatrixType, Dune::Amg::FirstDiagonal>>;
|
||||
using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
|
||||
int coarsenTarget = prm.get<int>("coarsenTarget");
|
||||
int ml = prm.get<int>("maxlevel");
|
||||
Criterion criterion(15, coarsenTarget);
|
||||
criterion.setDefaultValuesIsotropic(2);
|
||||
criterion.setAlpha(prm.get<double>("alpha"));
|
||||
criterion.setBeta(prm.get<double>("beta"));
|
||||
criterion.setMaxLevel(ml);
|
||||
criterion.setSkipIsolated(false);
|
||||
criterion.setDebugLevel(prm.get<int>("verbosity"));
|
||||
if (prm.get<std::string>("type") == "famg") {
|
||||
Dune::Amg::Parameters parms;
|
||||
parms.setNoPreSmoothSteps(1);
|
||||
parms.setNoPostSmoothSteps(1);
|
||||
return wrapPreconditioner<Dune::Amg::FastAMG<OperatorType, VectorType>>(linearoperator, criterion, parms);
|
||||
} else {
|
||||
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
|
||||
SmootherArgs smootherArgs;
|
||||
smootherArgs.iterations = prm.get<int>("iterations");
|
||||
// smootherArgs.overlap=SmootherArgs::vertex;
|
||||
// smootherArgs.overlap=SmootherArgs::none;
|
||||
// smootherArgs.overlap=SmootherArgs::aggregate;
|
||||
smootherArgs.relaxationFactor = prm.get<double>("relaxation");
|
||||
return std::make_shared<Dune::Amg::AMGCPR<OperatorType, VectorType, Smoother>>(
|
||||
linearoperator, criterion, smootherArgs);
|
||||
}
|
||||
return std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>(nullptr);
|
||||
}
|
||||
|
||||
template <class Smoother, class OperatorType, class VectorType, class Comm>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makeParAmgPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm)
|
||||
{
|
||||
using MatrixType = typename OperatorType::matrix_type;
|
||||
using CriterionBase
|
||||
= Dune::Amg::AggregationCriterion<Dune::Amg::SymmetricMatrixDependency<MatrixType, Dune::Amg::FirstDiagonal>>;
|
||||
using Criterion = Dune::Amg::CoarsenCriterion<CriterionBase>;
|
||||
int coarsenTarget = prm.get<int>("coarsenTarget");
|
||||
int ml = prm.get<int>("maxlevel");
|
||||
Criterion criterion(15, coarsenTarget);
|
||||
criterion.setDefaultValuesIsotropic(2);
|
||||
criterion.setAlpha(prm.get<double>("alpha"));
|
||||
criterion.setBeta(prm.get<double>("beta"));
|
||||
criterion.setMaxLevel(ml);
|
||||
criterion.setSkipIsolated(false);
|
||||
criterion.setDebugLevel(prm.get<int>("verbosity"));
|
||||
if (prm.get<std::string>("type") == "famg") {
|
||||
throw std::runtime_error("The FastAMG preconditioner cannot be used in parallel.");
|
||||
} else {
|
||||
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
|
||||
SmootherArgs smootherArgs;
|
||||
smootherArgs.iterations = prm.get<int>("iterations");
|
||||
// smootherArgs.overlap=SmootherArgs::vertex;
|
||||
// smootherArgs.overlap=SmootherArgs::none;
|
||||
// smootherArgs.overlap=SmootherArgs::aggregate;
|
||||
smootherArgs.relaxationFactor = prm.get<double>("relaxation");
|
||||
return std::make_shared<Dune::Amg::AMGCPR<OperatorType, VectorType, Smoother, Comm>>(
|
||||
linearoperator, criterion, smootherArgs, comm);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template <class OperatorType, class VectorType>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makeAmgPreconditioners(const OperatorType& linearoperator, const boost::property_tree::ptree& prm)
|
||||
{
|
||||
using MatrixType = typename OperatorType::matrix_type;
|
||||
if (prm.get<std::string>("type") == "famg") {
|
||||
// smoother type should not be used
|
||||
return makeAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
|
||||
linearoperator, prm);
|
||||
}
|
||||
|
||||
std::string precond = prm.get<std::string>("smoother");
|
||||
if (precond == "ILU0") {
|
||||
return makeAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
|
||||
linearoperator, prm);
|
||||
} else if (precond == "Jac") {
|
||||
return makeAmgPreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
|
||||
linearoperator, prm);
|
||||
// } else if (precond == "GS") {
|
||||
// return makeAmgPreconditioner<
|
||||
// Dune::SeqGS<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(linearoperator, prm);
|
||||
} else if (precond == "SOR") {
|
||||
return makeAmgPreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
|
||||
linearoperator, prm);
|
||||
} else if (precond == "SSOR") {
|
||||
return makeAmgPreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
|
||||
linearoperator, prm);
|
||||
} else if (precond == "ILUn") {
|
||||
return makeAmgPreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>, OperatorType, VectorType>(
|
||||
linearoperator, prm);
|
||||
} else {
|
||||
std::string msg("No such sequential preconditioner : ");
|
||||
msg += precond;
|
||||
throw std::runtime_error(msg);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template <class OperatorType, class VectorType, class Comm>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makeParAmgPreconditioners(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm)
|
||||
{
|
||||
if (prm.get<std::string>("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<std::string>("smoother");
|
||||
if (precond == "ILU0") {
|
||||
using SmootherType = Opm::ParallelOverlappingILU0<MatrixType, VectorType, VectorType, Comm>;
|
||||
return makeParAmgPreconditioner<SmootherType, OperatorType, VectorType, Comm>(linearoperator, prm, comm);
|
||||
/*
|
||||
return makeParAmgPreconditioner<Dune::SeqILU0<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
|
||||
linearoperator, prm, comm);
|
||||
} else if (precond == "Jac") {
|
||||
return makeParAmgPreconditioner<Dune::SeqJac<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
|
||||
linearoperator, prm, comm);
|
||||
// } else if (precond == "GS") {
|
||||
// return makeParAmgPreconditioner<
|
||||
// Dune::SeqGS<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(linearoperator, prm, comm);
|
||||
} else if (precond == "SOR") {
|
||||
return makeParAmgPreconditioner<Dune::SeqSOR<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
|
||||
linearoperator, prm, comm);
|
||||
} else if (precond == "SSOR") {
|
||||
return makeParAmgPreconditioner<Dune::SeqSSOR<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
|
||||
linearoperator, prm, comm);
|
||||
} else if (precond == "ILUn") {
|
||||
return makeParAmgPreconditioner<Dune::SeqILUn<MatrixType, VectorType, VectorType>, MatrixType, VectorType>(
|
||||
linearoperator, prm, comm);
|
||||
*/
|
||||
} else {
|
||||
std::string msg("No such parallel preconditioner : ");
|
||||
msg += precond;
|
||||
throw std::runtime_error(msg);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class OperatorType, class VectorType>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makeTwoLevelPreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm)
|
||||
{
|
||||
if (prm.get<std::string>("type") == "cpr") {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, false>>(linearoperator, prm);
|
||||
} else if (prm.get<std::string>("type") == "cprt") {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, true>>(linearoperator, prm);
|
||||
} else {
|
||||
std::string msg("Wrong cpr Should not happen");
|
||||
throw std::runtime_error(msg);
|
||||
}
|
||||
}
|
||||
|
||||
template <class OperatorType, class VectorType, class Comm>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makeParTwoLevelPreconditioner(const OperatorType& linearoperator,
|
||||
const boost::property_tree::ptree& prm,
|
||||
const Comm& comm)
|
||||
{
|
||||
if (prm.get<std::string>("type") == "cpr") {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, false, Comm>>(
|
||||
linearoperator, prm, comm);
|
||||
} else if (prm.get<std::string>("type") == "cprt") {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<OperatorType, VectorType, true, Comm>>(
|
||||
linearoperator, prm, comm);
|
||||
} else {
|
||||
std::string msg("Wrong cpr Should not happen");
|
||||
throw std::runtime_error(msg);
|
||||
}
|
||||
}
|
||||
|
||||
template <class OperatorType, class VectorType>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm)
|
||||
{
|
||||
if ((prm.get<std::string>("type") == "famg") or (prm.get<std::string>("type") == "amg")) {
|
||||
return makeAmgPreconditioners<OperatorType, VectorType>(linearoperator, prm);
|
||||
} else if ((prm.get<std::string>("type") == "cpr")
|
||||
or (prm.get<std::string>("type") == "cprt")) {
|
||||
return makeTwoLevelPreconditioner<OperatorType, VectorType>(linearoperator, prm);
|
||||
} else {
|
||||
return makeSeqPreconditioner<OperatorType, VectorType>(linearoperator, prm);
|
||||
}
|
||||
}
|
||||
|
||||
template <class OperatorType, class VectorType, class Comm>
|
||||
std::shared_ptr<Dune::PreconditionerWithUpdate<VectorType, VectorType>>
|
||||
makePreconditioner(const OperatorType& linearoperator, const boost::property_tree::ptree& prm, const Comm& comm)
|
||||
{
|
||||
if ((prm.get<std::string>("type") == "famg") or (prm.get<std::string>("type") == "amg")) {
|
||||
return makeParAmgPreconditioners<OperatorType, VectorType, Comm>(linearoperator, prm, comm);
|
||||
} else if ((prm.get<std::string>("type") == "cpr")
|
||||
or (prm.get<std::string>("type") == "cprt")) {
|
||||
return makeParTwoLevelPreconditioner<OperatorType, VectorType, Comm>(linearoperator, prm, comm);
|
||||
} else {
|
||||
return makeParPreconditioner<OperatorType, VectorType, Comm>(linearoperator, prm, comm);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} // namespace Dune
|
||||
|
||||
|
||||
#endif // OPM_MAKEPRECONDITIONER_HEADER_INCLUDED
|
||||
9
tests/options_flexiblesolver_simple.json
Normal file
9
tests/options_flexiblesolver_simple.json
Normal file
@ -0,0 +1,9 @@
|
||||
{
|
||||
"tol": "1e-12",
|
||||
"maxiter": "200",
|
||||
"verbosity": "0",
|
||||
"solver": "bicgstab",
|
||||
"preconditioner": {
|
||||
"type": "nothing"
|
||||
}
|
||||
}
|
||||
215
tests/test_preconditionerfactory.cpp
Normal file
215
tests/test_preconditionerfactory.cpp
Normal file
@ -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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <config.h>
|
||||
|
||||
#define BOOST_TEST_MODULE OPM_test_PreconditionerFactory
|
||||
#include <boost/test/unit_test.hpp>
|
||||
|
||||
#include <opm/simulators/linalg/PreconditionerFactory.hpp>
|
||||
#include <opm/simulators/linalg/FlexibleSolver.hpp>
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/istl/bvector.hh>
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/matrixmarket.hh>
|
||||
#include <dune/istl/solvers.hh>
|
||||
|
||||
#include <boost/property_tree/json_parser.hpp>
|
||||
#include <boost/property_tree/ptree.hpp>
|
||||
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
|
||||
|
||||
template <class X>
|
||||
class NothingPreconditioner : public Dune::Preconditioner<X, X>
|
||||
{
|
||||
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 <int bz>
|
||||
Dune::BlockVector<Dune::FieldVector<double, bz>>
|
||||
testPrec(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 = Dune::MatrixAdapter<Matrix, Vector, Vector>;
|
||||
Operator op(matrix);
|
||||
using PrecFactory = Dune::PreconditionerFactory<Operator>;
|
||||
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;
|
||||
}
|
||||
|
||||
namespace pt = boost::property_tree;
|
||||
|
||||
void test1(const pt::ptree& prm)
|
||||
{
|
||||
const int bz = 1;
|
||||
auto sol = testPrec<bz>(prm, "matr33.txt", "rhs3.txt");
|
||||
Dune::BlockVector<Dune::FieldVector<double, bz>> 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<bz>(prm, "matr33.txt", "rhs3.txt");
|
||||
Dune::BlockVector<Dune::FieldVector<double, bz>> 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 <int bz>
|
||||
using M = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
|
||||
template <int bz>
|
||||
using V = Dune::BlockVector<Dune::FieldVector<double, bz>>;
|
||||
template <int bz>
|
||||
using O = Dune::MatrixAdapter<M<bz>, V<bz>, V<bz>>;
|
||||
template <int bz>
|
||||
using PF = Dune::PreconditionerFactory<O<bz>>;
|
||||
|
||||
|
||||
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<bz>(prm, "matr33.txt", "rhs3.txt"), std::runtime_error);
|
||||
}
|
||||
|
||||
// Test with 3x3 block solvers.
|
||||
{
|
||||
const int bz = 3;
|
||||
BOOST_CHECK_THROW(testPrec<bz>(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<NothingPreconditioner<V<1>>>();
|
||||
});
|
||||
|
||||
|
||||
// Test with 1x1 block solvers.
|
||||
test1(prm);
|
||||
|
||||
// Test with 3x3 block solvers.
|
||||
{
|
||||
const int bz = 3;
|
||||
BOOST_CHECK_THROW(testPrec<bz>(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<NothingPreconditioner<V<3>>>();
|
||||
});
|
||||
|
||||
// Test with 1x1 block solvers.
|
||||
test1(prm);
|
||||
|
||||
// Test with 3x3 block solvers.
|
||||
test3(prm);
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user