mirror of
https://github.com/OPM/opm-simulators.git
synced 2024-11-26 03:00:17 -06:00
Make update method of preconditioners parameter-less again.
Previously, it got passed the weights only needed for CPR. Additionally those were passed with the parameter tree to the update method and constructor. Now the CPR constructor gets a function to use for recalculating the weights and the property is not changed. Unfortunately this means that the preconditioner creators of the factory get another parameter.
This commit is contained in:
Markus Blatt
parent
95a1e1ca0e
commit
dcb316f442
@ -46,16 +46,27 @@ public:
|
||||
using VectorType = VectorTypeT;
|
||||
|
||||
/// Create a sequential solver.
|
||||
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix)
|
||||
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix,
|
||||
const std::function<VectorTypeT()>& weightsCalculator = std::function<VectorTypeT()>())
|
||||
{
|
||||
init(prm, matrix, Dune::Amg::SequentialInformation());
|
||||
init(prm, matrix, weightsCalculator, Dune::Amg::SequentialInformation());
|
||||
}
|
||||
|
||||
/// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
|
||||
template <class Comm>
|
||||
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix, const Comm& comm)
|
||||
FlexibleSolver(const boost::property_tree::ptree& prm,
|
||||
const typename std::enable_if<std::is_function<Comm>::value,MatrixType>::type& matrix,
|
||||
const Comm& comm)
|
||||
{
|
||||
init(prm, matrix, comm);
|
||||
init(prm, matrix, std::function<VectorTypeT()>(), comm);
|
||||
}
|
||||
|
||||
/// Create a parallel solver (if Comm is e.g. OwnerOverlapCommunication).
|
||||
template <class Comm>
|
||||
FlexibleSolver(const boost::property_tree::ptree& prm, const MatrixType& matrix,
|
||||
const std::function<VectorTypeT()>& weightsCalculator, const Comm& comm)
|
||||
{
|
||||
init(prm, matrix, weightsCalculator, comm);
|
||||
}
|
||||
|
||||
virtual void apply(VectorType& x, VectorType& rhs, Dune::InverseOperatorResult& res) override
|
||||
@ -89,24 +100,28 @@ private:
|
||||
|
||||
// Machinery for making sequential or parallel operators/preconditioners/scalar products.
|
||||
template <class Comm>
|
||||
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm, const Comm& comm)
|
||||
void initOpPrecSp(const MatrixType& matrix, const boost::property_tree::ptree& prm,
|
||||
const std::function<VectorTypeT()> weightsCalculator, const Comm& comm)
|
||||
{
|
||||
// Parallel case.
|
||||
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
|
||||
auto linop = std::make_shared<ParOperatorType>(matrix, comm);
|
||||
linearoperator_ = linop;
|
||||
preconditioner_
|
||||
= Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, prm.get_child("preconditioner"), comm);
|
||||
= Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, prm.get_child("preconditioner"),
|
||||
weightsCalculator, comm);
|
||||
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, linearoperator_->category());
|
||||
}
|
||||
|
||||
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 std::function<VectorTypeT()> weightsCalculator, const Dune::Amg::SequentialInformation&)
|
||||
{
|
||||
// Sequential case.
|
||||
using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
|
||||
auto linop = std::make_shared<SeqOperatorType>(matrix);
|
||||
linearoperator_ = linop;
|
||||
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, prm.get_child("preconditioner"));
|
||||
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, prm.get_child("preconditioner"),
|
||||
weightsCalculator);
|
||||
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
|
||||
}
|
||||
|
||||
@ -155,9 +170,10 @@ private:
|
||||
// Main initialization routine.
|
||||
// Call with Comm == Dune::Amg::SequentialInformation to get a serial solver.
|
||||
template <class Comm>
|
||||
void init(const boost::property_tree::ptree& prm, const MatrixType& matrix, const Comm& comm)
|
||||
void init(const boost::property_tree::ptree& prm, const MatrixType& matrix,
|
||||
const std::function<VectorTypeT()> weightsCalculator, const Comm& comm)
|
||||
{
|
||||
initOpPrecSp(matrix, prm, comm);
|
||||
initOpPrecSp(matrix, prm, weightsCalculator, comm);
|
||||
initSolver(prm);
|
||||
}
|
||||
|
||||
|
||||
@ -151,7 +151,7 @@ public:
|
||||
// Never recreate solver.
|
||||
}
|
||||
|
||||
VectorType weights;
|
||||
std::function<VectorType()> weightsCalculator;
|
||||
|
||||
if( prm_.get<std::string>("preconditioner.type") == "cpr" ||
|
||||
prm_.get<std::string>("preconditioner.type") == "cprt"
|
||||
@ -163,20 +163,21 @@ public:
|
||||
}
|
||||
|
||||
if(prm_.get<std::string>("preconditioner.weight_type") == "quasiimpes") {
|
||||
if( not( recreate_solver || !solver_) ){
|
||||
// weighs will be created as default in the solver
|
||||
weights = Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
|
||||
mat.istlMatrix(),
|
||||
prm_.get<int>("preconditioner.pressure_var_index"), transpose);
|
||||
}
|
||||
// weighs will be created as default in the solver
|
||||
weightsCalculator =
|
||||
[&mat, this, transpose](){
|
||||
return Opm::Amg::getQuasiImpesWeights<MatrixType,
|
||||
VectorType>(
|
||||
mat.istlMatrix(),
|
||||
this->prm_.get<int>("preconditioner.pressure_var_index"),
|
||||
transpose);
|
||||
};
|
||||
|
||||
}else if(prm_.get<std::string>("preconditioner.weight_type") == "trueimpes" ){
|
||||
weights =
|
||||
this->getTrueImpesWeights(b, prm_.get<int>("preconditioner.pressure_var_index"));
|
||||
if( recreate_solver || !solver_){
|
||||
// need weights for the constructor
|
||||
prm_.put("preconditioner.weights",weights);
|
||||
}
|
||||
weightsCalculator =
|
||||
[this, &b](){
|
||||
return this->getTrueImpesWeights(b, this->prm_.get<int>("preconditioner.pressure_var_index"));
|
||||
};
|
||||
}else{
|
||||
throw std::runtime_error("no such weights implemented for cpr");
|
||||
}
|
||||
@ -185,14 +186,14 @@ public:
|
||||
if (recreate_solver || !solver_) {
|
||||
if (isParallel()) {
|
||||
#if HAVE_MPI
|
||||
solver_.reset(new SolverType(prm_, mat.istlMatrix(), *comm_));
|
||||
solver_.reset(new SolverType(prm_, mat.istlMatrix(), weightsCalculator, *comm_));
|
||||
#endif
|
||||
} else {
|
||||
solver_.reset(new SolverType(prm_, mat.istlMatrix()));
|
||||
solver_.reset(new SolverType(prm_, mat.istlMatrix(), weightsCalculator));
|
||||
}
|
||||
rhs_ = b;
|
||||
} else {
|
||||
solver_->preconditioner().update(weights);
|
||||
solver_->preconditioner().update();
|
||||
rhs_ = b;
|
||||
}
|
||||
}
|
||||
|
||||
@ -63,9 +63,9 @@ public:
|
||||
}
|
||||
|
||||
// The update() function does nothing for a wrapped preconditioner.
|
||||
virtual void update(const X& w) override
|
||||
virtual void update() override
|
||||
{
|
||||
orig_precond_.update(w);
|
||||
orig_precond_.update();
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
@ -84,17 +84,13 @@ public:
|
||||
using MatrixType = typename OperatorType::matrix_type;
|
||||
using PrecFactory = Opm::PreconditionerFactory<OperatorType, Communication>;
|
||||
|
||||
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm)
|
||||
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm,
|
||||
const std::function<VectorType()> weightsCalculator)
|
||||
: linear_operator_(linearoperator)
|
||||
, finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother")))
|
||||
, comm_(nullptr)
|
||||
, weights_(
|
||||
(prm.get<std::string>("weight_type") != "quasiimpes") ?
|
||||
prm.get<VectorType>("weights") :
|
||||
Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(linearoperator.getmat(),
|
||||
prm.get<int>("pressure_var_index"),
|
||||
transpose)
|
||||
)
|
||||
, weightsCalculator_(weightsCalculator)
|
||||
, weights_(weightsCalculator())
|
||||
, levelTransferPolicy_(dummy_comm_, weights_, prm.get<int>("pressure_var_index"))
|
||||
, coarseSolverPolicy_(prm.get_child("coarsesolver"))
|
||||
, twolevel_method_(linearoperator,
|
||||
@ -114,17 +110,13 @@ public:
|
||||
}
|
||||
}
|
||||
|
||||
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm, const Communication& comm)
|
||||
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm,
|
||||
const std::function<VectorType()> weightsCalculator, const Communication& comm)
|
||||
: linear_operator_(linearoperator)
|
||||
, finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother"), comm))
|
||||
, comm_(&comm)
|
||||
, weights_(
|
||||
(prm.get<std::string>("weight_type") != "quasiimpes") ?
|
||||
prm.get<VectorType>("weights") :
|
||||
Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(linearoperator.getmat(),
|
||||
prm.get<int>("pressure_var_index"),
|
||||
transpose)
|
||||
)
|
||||
, weightsCalculator_(weightsCalculator)
|
||||
, weights_(weightsCalculator())
|
||||
, levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index"))
|
||||
, coarseSolverPolicy_(prm.get_child("coarsesolver"))
|
||||
, twolevel_method_(linearoperator,
|
||||
@ -161,9 +153,9 @@ public:
|
||||
twolevel_method_.post(x);
|
||||
}
|
||||
|
||||
virtual void update(const VectorType& weights) override
|
||||
virtual void update() override
|
||||
{
|
||||
weights_ = weights;
|
||||
weights_ = weightsCalculator_();
|
||||
updateImpl(comm_);
|
||||
}
|
||||
|
||||
@ -207,6 +199,7 @@ private:
|
||||
const OperatorType& linear_operator_;
|
||||
std::shared_ptr<Dune::Preconditioner<VectorType, VectorType>> finesmoother_;
|
||||
const Communication* comm_;
|
||||
std::function<VectorType()> weightsCalculator_;
|
||||
VectorType weights_;
|
||||
LevelTransferPolicy levelTransferPolicy_;
|
||||
CoarseSolverPolicy coarseSolverPolicy_;
|
||||
|
||||
@ -921,7 +921,7 @@ public:
|
||||
DUNE_UNUSED_PARAMETER(x);
|
||||
}
|
||||
|
||||
virtual void update(const Range& = Range()) override
|
||||
virtual void update() override
|
||||
{
|
||||
// (For older DUNE versions the communicator might be
|
||||
// invalid if redistribution in AMG happened on the coarset level.
|
||||
|
||||
@ -57,16 +57,30 @@ public:
|
||||
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&)>;
|
||||
using Creator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&, const std::function<Vector()>&)>;
|
||||
using ParCreator = std::function<PrecPtr(const Operator&, const boost::property_tree::ptree&, const std::function<Vector()>&, 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.
|
||||
/// \param weightsCalculator Calculator for weights used in CPR.
|
||||
/// \return (smart) pointer to the created preconditioner.
|
||||
static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm)
|
||||
static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm,
|
||||
const std::function<Vector()>& weightsCalculator = std::function<Vector()>())
|
||||
{
|
||||
return instance().doCreate(op, prm);
|
||||
return instance().doCreate(op, prm, weightsCalculator);
|
||||
}
|
||||
|
||||
/// 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).
|
||||
/// \param weightsCalculator Calculator for weights used in CPR.
|
||||
/// \return (smart) pointer to the created preconditioner.
|
||||
static PrecPtr create(const Operator& op, const boost::property_tree::ptree& prm,
|
||||
const std::function<Vector()>& weightsCalculator, const Comm& comm)
|
||||
{
|
||||
return instance().doCreate(op, prm, weightsCalculator, comm);
|
||||
}
|
||||
|
||||
/// Create a new parallel preconditioner and return a pointer to it.
|
||||
@ -76,9 +90,8 @@ public:
|
||||
/// \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);
|
||||
return instance().doCreate(op, prm, std::function<Vector()>(), 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
|
||||
@ -154,44 +167,45 @@ private:
|
||||
using V = Vector;
|
||||
using P = boost::property_tree::ptree;
|
||||
using C = Comm;
|
||||
doAddCreator("ILU0", [](const O& op, const P& prm, const C& comm) {
|
||||
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
|
||||
op.getmat(), comm, 0, w, Opm::MILU_VARIANT::ILU);
|
||||
});
|
||||
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const C& comm) {
|
||||
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&, 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 std::make_shared<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) {
|
||||
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
|
||||
const int n = prm.get<int>("ilulevel");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V, C>>(
|
||||
op.getmat(), comm, n, w, Opm::MILU_VARIANT::ILU);
|
||||
});
|
||||
doAddCreator("Jac", [](const O& op, const P& prm, const C& comm) {
|
||||
doAddCreator("Jac", [](const O& op, const P& prm, const std::function<Vector()>&,
|
||||
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) {
|
||||
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&, 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) {
|
||||
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&, 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) {
|
||||
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&, 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) {
|
||||
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
|
||||
const std::string smoother = prm.get<std::string>("smoother");
|
||||
if (smoother == "ILU0") {
|
||||
using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
|
||||
@ -204,11 +218,13 @@ private:
|
||||
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("cpr", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
|
||||
assert(weightsCalculator);
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false, Comm>>(op, prm, weightsCalculator, comm);
|
||||
});
|
||||
doAddCreator("cprt", [](const O& op, const P& prm, const C& comm) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, comm);
|
||||
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()> weightsCalculator, const C& comm) {
|
||||
assert(weightsCalculator);
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true, Comm>>(op, prm, weightsCalculator, comm);
|
||||
});
|
||||
}
|
||||
|
||||
@ -221,43 +237,43 @@ private:
|
||||
using M = Matrix;
|
||||
using V = Vector;
|
||||
using P = boost::property_tree::ptree;
|
||||
doAddCreator("ILU0", [](const O& op, const P& prm) {
|
||||
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
|
||||
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
|
||||
});
|
||||
doAddCreator("ParOverILU0", [](const O& op, const P& prm) {
|
||||
doAddCreator("ParOverILU0", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
|
||||
op.getmat(), 0, w, Opm::MILU_VARIANT::ILU);
|
||||
});
|
||||
doAddCreator("ILUn", [](const O& op, const P& prm) {
|
||||
doAddCreator("ILUn", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const int n = prm.get<int>("ilulevel");
|
||||
const double w = prm.get<double>("relaxation");
|
||||
return std::make_shared<Opm::ParallelOverlappingILU0<M, V, V>>(
|
||||
op.getmat(), n, w, Opm::MILU_VARIANT::ILU);
|
||||
});
|
||||
doAddCreator("Jac", [](const O& op, const P& prm) {
|
||||
doAddCreator("Jac", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
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) {
|
||||
doAddCreator("GS", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
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) {
|
||||
doAddCreator("SOR", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
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) {
|
||||
doAddCreator("SSOR", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
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) {
|
||||
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
const std::string smoother = prm.get<std::string>("smoother");
|
||||
if (smoother == "ILU0") {
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
|
||||
@ -288,18 +304,18 @@ private:
|
||||
throw std::runtime_error(msg);
|
||||
}
|
||||
});
|
||||
doAddCreator("famg", [](const O& op, const P& prm) {
|
||||
doAddCreator("famg", [](const O& op, const P& prm, const std::function<Vector()>&) {
|
||||
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("cpr", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, false>>(op, prm, weightsCalculator);
|
||||
});
|
||||
doAddCreator("cprt", [](const O& op, const P& prm) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(op, prm);
|
||||
doAddCreator("cprt", [](const O& op, const P& prm, const std::function<Vector()>& weightsCalculator) {
|
||||
return std::make_shared<OwningTwoLevelPreconditioner<O, V, true>>(op, prm, weightsCalculator);
|
||||
});
|
||||
}
|
||||
|
||||
@ -320,7 +336,8 @@ private:
|
||||
}
|
||||
|
||||
// Actually creates the product object.
|
||||
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm)
|
||||
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm,
|
||||
const std::function<Vector()> weightsCalculator)
|
||||
{
|
||||
const std::string& type = prm.get<std::string>("type");
|
||||
auto it = creators_.find(type);
|
||||
@ -333,10 +350,11 @@ private:
|
||||
msg << std::endl;
|
||||
throw std::runtime_error(msg.str());
|
||||
}
|
||||
return it->second(op, prm);
|
||||
return it->second(op, prm, weightsCalculator);
|
||||
}
|
||||
|
||||
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm, const Comm& comm)
|
||||
PrecPtr doCreate(const Operator& op, const boost::property_tree::ptree& prm,
|
||||
const std::function<Vector()> weightsCalculator, const Comm& comm)
|
||||
{
|
||||
const std::string& type = prm.get<std::string>("type");
|
||||
auto it = parallel_creators_.find(type);
|
||||
@ -350,7 +368,7 @@ private:
|
||||
msg << std::endl;
|
||||
throw std::runtime_error(msg.str());
|
||||
}
|
||||
return it->second(op, prm, comm);
|
||||
return it->second(op, prm, weightsCalculator, comm);
|
||||
}
|
||||
|
||||
// Actually adds the creator.
|
||||
|
||||
@ -31,7 +31,7 @@ template <class X, class Y>
|
||||
class PreconditionerWithUpdate : public Preconditioner<X, Y>
|
||||
{
|
||||
public:
|
||||
virtual void update(const X& w) = 0;
|
||||
virtual void update() = 0;
|
||||
};
|
||||
|
||||
template <class OriginalPreconditioner>
|
||||
@ -69,7 +69,7 @@ public:
|
||||
}
|
||||
|
||||
// The update() function does nothing for a wrapped preconditioner.
|
||||
virtual void update(const X& /*w*/) override
|
||||
virtual void update() override
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
@ -61,9 +61,9 @@ namespace Amg
|
||||
: linsolver_()
|
||||
{
|
||||
if (op.category() == Dune::SolverCategory::overlapping) {
|
||||
linsolver_.reset(new Solver(prm, op.getmat(), comm));
|
||||
linsolver_.reset(new Solver(prm, op.getmat(), std::function<X()>(), comm));
|
||||
} else {
|
||||
linsolver_.reset(new Solver(prm, op.getmat()));
|
||||
linsolver_.reset(new Solver(prm, op.getmat(), std::function<X()>()));
|
||||
}
|
||||
}
|
||||
|
||||
@ -84,8 +84,7 @@ namespace Amg
|
||||
|
||||
void updatePreconditioner()
|
||||
{
|
||||
X w;
|
||||
linsolver_->preconditioner().update(w);
|
||||
linsolver_->preconditioner().update();
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
@ -239,9 +239,8 @@ namespace Dune
|
||||
/**
|
||||
* @brief Update the coarse solver and the hierarchies.
|
||||
*/
|
||||
virtual void update(const X& w);
|
||||
|
||||
virtual void update();
|
||||
|
||||
/**
|
||||
* @brief Check whether the coarse solver used is a direct solver.
|
||||
* @return True if the coarse level solver is a direct solver.
|
||||
@ -469,11 +468,6 @@ namespace Dune
|
||||
update();
|
||||
}
|
||||
|
||||
template<class M, class X, class S, class PI, class A>
|
||||
void AMGCPR<M,X,S,PI,A>::update(const X& /*w*/)
|
||||
{
|
||||
update();
|
||||
}
|
||||
template<class M, class X, class S, class PI, class A>
|
||||
void AMGCPR<M,X,S,PI,A>::update()
|
||||
{
|
||||
|
||||
@ -58,7 +58,19 @@ testSolver(const boost::property_tree::ptree& prm, const std::string& matrix_fil
|
||||
}
|
||||
readMatrixMarket(rhs, rhsfile);
|
||||
}
|
||||
Dune::FlexibleSolver<Matrix, Vector> solver(prm, matrix);
|
||||
bool transpose = false;
|
||||
|
||||
if(prm.get<std::string>("preconditioner.type") == "cprt"){
|
||||
transpose = true;
|
||||
}
|
||||
auto wc = [&matrix, &prm, transpose]()
|
||||
{
|
||||
return Opm::Amg::getQuasiImpesWeights<Matrix,
|
||||
Vector>(matrix,
|
||||
prm.get<int>("preconditioner.pressure_var_index"),
|
||||
transpose);
|
||||
};
|
||||
Dune::FlexibleSolver<Matrix, Vector> solver(prm, matrix, wc);
|
||||
Vector x(rhs.size());
|
||||
Dune::InverseOperatorResult res;
|
||||
solver.apply(x, rhs, res);
|
||||
|
||||
@ -93,7 +93,19 @@ testPrec(const boost::property_tree::ptree& prm, const std::string& matrix_filen
|
||||
using Operator = Dune::MatrixAdapter<Matrix, Vector, Vector>;
|
||||
Operator op(matrix);
|
||||
using PrecFactory = Opm::PreconditionerFactory<Operator>;
|
||||
auto prec = PrecFactory::create(op, prm.get_child("preconditioner"));
|
||||
bool transpose = false;
|
||||
|
||||
if(prm.get<std::string>("preconditioner.type") == "cprt"){
|
||||
transpose = true;
|
||||
}
|
||||
auto wc = [&matrix, &prm, transpose]()
|
||||
{
|
||||
return Opm::Amg::getQuasiImpesWeights<Matrix,
|
||||
Vector>(matrix,
|
||||
prm.get<int>("preconditioner.pressure_var_index"),
|
||||
transpose);
|
||||
};
|
||||
auto prec = PrecFactory::create(op, prm.get_child("preconditioner"), wc);
|
||||
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;
|
||||
@ -194,7 +206,7 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner)
|
||||
|
||||
|
||||
// Add preconditioner to factory for block size 1.
|
||||
PF<1>::addCreator("nothing", [](const O<1>&, const pt::ptree&) {
|
||||
PF<1>::addCreator("nothing", [](const O<1>&, const pt::ptree&, const std::function<V<1>()>&) {
|
||||
return Dune::wrapPreconditioner<NothingPreconditioner<V<1>>>();
|
||||
});
|
||||
|
||||
@ -209,7 +221,7 @@ BOOST_AUTO_TEST_CASE(TestAddingPreconditioner)
|
||||
}
|
||||
|
||||
// Add preconditioner to factory for block size 3.
|
||||
PF<3>::addCreator("nothing", [](const O<3>&, const pt::ptree&) {
|
||||
PF<3>::addCreator("nothing", [](const O<3>&, const pt::ptree&, const std::function<V<3>()>&) {
|
||||
return Dune::wrapPreconditioner<NothingPreconditioner<V<3>>>();
|
||||
});
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user