Add sensible defaults for flexible solver properties...

to prevent throwing. The get methods will throw if called without a
default value. This quite unfortunate and not very user friendly, as
there are many properties and the throwing will happen during the
linear solve and result in time step chopping.

This commit should prevent such throws and allow users to provide
jsdon files omitting some options.
This commit is contained in:
Markus Blatt 2020-03-25 21:05:01 +01:00
parent fce9f5e57c
commit ff0d54d4ea
4 changed files with 87 additions and 74 deletions

View File

@ -105,10 +105,12 @@ private:
{
// Parallel case.
using ParOperatorType = Dune::OverlappingSchwarzOperator<MatrixType, VectorType, VectorType, Comm>;
using pt = const boost::property_tree::ptree;
auto linop = std::make_shared<ParOperatorType>(matrix, comm);
linearoperator_ = linop;
auto child = prm.get_child_optional("preconditioner");
preconditioner_
= Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, prm.get_child("preconditioner"),
= Opm::PreconditionerFactory<ParOperatorType, Comm>::create(*linop, child? *child : pt(),
weightsCalculator, comm);
scalarproduct_ = Dune::createScalarProduct<VectorType, Comm>(comm, linearoperator_->category());
}
@ -118,19 +120,21 @@ private:
{
// Sequential case.
using SeqOperatorType = Dune::MatrixAdapter<MatrixType, VectorType, VectorType>;
using pt = const boost::property_tree::ptree;
auto linop = std::make_shared<SeqOperatorType>(matrix);
linearoperator_ = linop;
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, prm.get_child("preconditioner"),
auto child = prm.get_child_optional("preconditioner");
preconditioner_ = Opm::PreconditionerFactory<SeqOperatorType>::create(*linop, child? *child : pt(),
weightsCalculator);
scalarproduct_ = std::make_shared<Dune::SeqScalarProduct<VectorType>>();
}
void initSolver(const boost::property_tree::ptree& prm)
{
const double tol = prm.get<double>("tol");
const int maxiter = prm.get<int>("maxiter");
const int verbosity = prm.get<int>("verbosity");
const std::string solver_type = prm.get<std::string>("solver");
const double tol = prm.get<double>("tol", 1e-2);
const int maxiter = prm.get<int>("maxiter", 200);
const int verbosity = prm.get<int>("verbosity", 0);
const std::string solver_type = prm.get<std::string>("solver", "bicgstab");
if (solver_type == "bicgstab") {
linsolver_.reset(new Dune::BiCGSTABSolver<VectorType>(*linearoperator_,
*scalarproduct_,
@ -146,7 +150,7 @@ private:
maxiter, // maximum number of iterations
verbosity));
} else if (solver_type == "gmres") {
int restart = prm.get<int>("restart");
int restart = prm.get<int>("restart", 15);
linsolver_.reset(new Dune::RestartedGMResSolver<VectorType>(*linearoperator_,
*scalarproduct_,
*preconditioner_,

View File

@ -153,30 +153,33 @@ public:
std::function<VectorType()> weightsCalculator;
if( prm_.get<std::string>("preconditioner.type") == "cpr" ||
prm_.get<std::string>("preconditioner.type") == "cprt"
auto preconditionerType = prm_.get("preconditioner.type", "cpr");
if( preconditionerType == "cpr" ||
preconditionerType == "cprt"
)
{
bool transpose = false;
if(prm_.get<std::string>("preconditioner.type") == "cprt"){
if(preconditionerType == "cprt"){
transpose = true;
}
if(prm_.get<std::string>("preconditioner.weight_type") == "quasiimpes") {
auto weightsType = prm_.get("preconditioner.weight_type", "quasiimpes");
auto pressureIndex = this->prm_.get("preconditioner.pressure_var_index", 1);
if(weightsType == "quasiimpes") {
// weighs will be created as default in the solver
weightsCalculator =
[&mat, this, transpose](){
[&mat, this, transpose, pressureIndex](){
return Opm::Amg::getQuasiImpesWeights<MatrixType,
VectorType>(
mat.istlMatrix(),
this->prm_.get<int>("preconditioner.pressure_var_index"),
pressureIndex,
transpose);
};
}else if(prm_.get<std::string>("preconditioner.weight_type") == "trueimpes" ){
}else if(weightsType == "trueimpes" ){
weightsCalculator =
[this, &b](){
return this->getTrueImpesWeights(b, this->prm_.get<int>("preconditioner.pressure_var_index"));
[this, &b, pressureIndex](){
return this->getTrueImpesWeights(b, pressureIndex);
};
}else{
throw std::runtime_error("no such weights implemented for cpr");

View File

@ -87,22 +87,24 @@ public:
OwningTwoLevelPreconditioner(const OperatorType& linearoperator, const pt& prm,
const std::function<VectorType()> weightsCalculator)
: linear_operator_(linearoperator)
, finesmoother_(PrecFactory::create(linearoperator, prm.get_child("finesmoother")))
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt()))
, comm_(nullptr)
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(dummy_comm_, weights_, prm.get<int>("pressure_var_index"))
, coarseSolverPolicy_(prm.get_child("coarsesolver"))
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
coarseSolverPolicy_,
prm.get<int>("pre_smooth"),
prm.get<int>("post_smooth"))
prm.get<int>("pre_smooth", transpose? 1 : 0),
prm.get<int>("post_smooth", transpose? 0 : 1))
, prm_(prm)
{
if (prm.get<int>("verbosity") > 10) {
std::ofstream outfile(prm.get<std::string>("weights_filename"));
if (prm.get<int>("verbosity", 0) > 10) {
std::ofstream outfile(prm.get<std::string>("weights_filename", "impes_weights.txt"));
if (!outfile) {
throw std::runtime_error("Could not write weights");
}
@ -113,24 +115,26 @@ public:
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))
, finesmoother_(PrecFactory::create(linearoperator,
prm.get_child_optional("finesmoother")?
prm.get_child("finesmoother"): pt(), comm))
, comm_(&comm)
, weightsCalculator_(weightsCalculator)
, weights_(weightsCalculator())
, levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index"))
, coarseSolverPolicy_(prm.get_child("coarsesolver"))
, levelTransferPolicy_(*comm_, weights_, prm.get<int>("pressure_var_index", 1))
, coarseSolverPolicy_(prm.get_child_optional("coarsesolver")? prm.get_child("coarsesolver") : pt())
, twolevel_method_(linearoperator,
finesmoother_,
levelTransferPolicy_,
coarseSolverPolicy_,
transpose ? 1 : 0,
transpose ? 0 : 1)
prm.get<int>("pre_smooth", transpose? 1 : 0),
prm.get<int>("post_smooth", transpose? 0 : 1))
, prm_(prm)
{
//Opm::Amg::getQuasiImpesWeights<MatrixType, VectorType>(
// linearoperator.getmat(), prm.get<int>("pressure_var_index"), transpose))
if (prm.get<int>("verbosity") > 10) {
std::ofstream outfile(prm.get<std::string>("weights_filename"));
if (prm.get<int>("verbosity", 0) > 10) {
std::ofstream outfile(prm.get<std::string>("weights_filename", "impes_weights.txt"));
if (!outfile) {
throw std::runtime_error("Could not write weights");
}
@ -185,14 +189,16 @@ private:
void updateImpl(const Comm*)
{
// Parallel case.
finesmoother_ = PrecFactory::create(linear_operator_, prm_.get_child("finesmoother"), *comm_);
auto child = prm_.get_child_optional("finesmoother");
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt(), *comm_);
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
}
void updateImpl(const Dune::Amg::SequentialInformation*)
{
// Serial case.
finesmoother_ = PrecFactory::create(linear_operator_, prm_.get_child("finesmoother"));
auto child = prm_.get_child_optional("finesmoother");
finesmoother_ = PrecFactory::create(linear_operator_, child ? *child : pt());
twolevel_method_.updatePreconditioner(finesmoother_, coarseSolverPolicy_);
}

View File

@ -125,15 +125,15 @@ private:
// Helpers for creation of AMG preconditioner.
static Criterion amgCriterion(const boost::property_tree::ptree& prm)
{
Criterion criterion(15, prm.get<int>("coarsenTarget"));
Criterion criterion(15, prm.get<int>("coarsenTarget", 1200));
criterion.setDefaultValuesIsotropic(2);
criterion.setAlpha(prm.get<double>("alpha"));
criterion.setBeta(prm.get<double>("beta"));
criterion.setMaxLevel(prm.get<int>("maxlevel"));
criterion.setSkipIsolated(prm.get<bool>("skip_isolated"));
criterion.setNoPreSmoothSteps(prm.get<int>("pre_smooth"));
criterion.setNoPostSmoothSteps(prm.get<int>("post_smooth"));
criterion.setDebugLevel(prm.get<int>("verbosity"));
criterion.setAlpha(prm.get<double>("alpha", 0.33));
criterion.setBeta(prm.get<double>("beta", 1e-5));
criterion.setMaxLevel(prm.get<int>("maxlevel", 15));
criterion.setSkipIsolated(prm.get<bool>("skip_isolated", false));
criterion.setNoPreSmoothSteps(prm.get<int>("pre_smooth", 1));
criterion.setNoPostSmoothSteps(prm.get<int>("post_smooth", 1));
criterion.setDebugLevel(prm.get<int>("verbosity", 0));
return criterion;
}
@ -142,11 +142,11 @@ private:
{
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
SmootherArgs smootherArgs;
smootherArgs.iterations = prm.get<int>("iterations");
smootherArgs.iterations = prm.get<int>("iterations", 1);
// smootherArgs.overlap=SmootherArgs::vertex;
// smootherArgs.overlap=SmootherArgs::none;
// smootherArgs.overlap=SmootherArgs::aggregate;
smootherArgs.relaxationFactor = prm.get<double>("relaxation");
smootherArgs.relaxationFactor = prm.get<double>("relaxation", 0.9);
return smootherArgs;
}
@ -161,8 +161,8 @@ private:
Dune::Amg::KAMG< Operator, Vector, Smoother>
>
>(op, crit, sargs,
prm.get<size_t>("max_krylov"),
prm.get<double>("min_reduction") );
prm.get<size_t>("max_krylov", 1),
prm.get<double>("min_reduction", 1e-1) );
}else{
return std::make_shared<Dune::Amg::AMGCPR<Operator, Vector, Smoother>>(op, crit, sargs);
}
@ -181,45 +181,45 @@ private:
using P = boost::property_tree::ptree;
using C = Comm;
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const double w = prm.get<double>("relaxation");
const double w = prm.get<double>("relaxation", 1.0);
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 std::function<Vector()>&, const C& comm) {
const double w = prm.get<double>("relaxation");
const double w = prm.get<double>("relaxation", 1.0);
// 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 std::function<Vector()>&, const C& comm) {
const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation");
const int n = prm.get<int>("ilulevel", 0);
const double w = prm.get<double>("relaxation", 1.0);
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 std::function<Vector()>&,
const C& comm) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqJac<M, V, V>>>(comm, op.getmat(), n, w);
});
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");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqGS<M, V, V>>>(comm, op.getmat(), n, w);
});
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");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSOR<M, V, V>>>(comm, op.getmat(), n, w);
});
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");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapBlockPreconditioner<DummyUpdatePreconditioner<SeqSSOR<M, V, V>>>(comm, op.getmat(), n, w);
});
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&, const C& comm) {
const std::string smoother = prm.get<std::string>("smoother");
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
using Smoother = Opm::ParallelOverlappingILU0<M, V, V, C>;
auto crit = amgCriterion(prm);
@ -251,44 +251,44 @@ private:
using V = Vector;
using P = boost::property_tree::ptree;
doAddCreator("ILU0", [](const O& op, const P& prm, const std::function<Vector()>&) {
const double w = prm.get<double>("relaxation");
const double w = prm.get<double>("relaxation", 1.0);
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, const std::function<Vector()>&) {
const double w = prm.get<double>("relaxation");
const double w = prm.get<double>("relaxation", 1.0);
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, const std::function<Vector()>&) {
const int n = prm.get<int>("ilulevel");
const double w = prm.get<double>("relaxation");
const int n = prm.get<int>("ilulevel", 0);
const double w = prm.get<double>("relaxation", 1.0);
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, const std::function<Vector()>&) {
const int n = prm.get<int>("repeats");
const double w = prm.get<double>("relaxation");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqJac<M, V, V>>(op.getmat(), n, w);
});
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");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqGS<M, V, V>>(op.getmat(), n, w);
});
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");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqSOR<M, V, V>>(op.getmat(), n, w);
});
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");
const int n = prm.get<int>("repeats", 1);
const double w = prm.get<double>("relaxation", 1.0);
return wrapPreconditioner<SeqSSOR<M, V, V>>(op.getmat(), n, w);
});
doAddCreator("amg", [](const O& op, const P& prm, const std::function<Vector()>&) {
const std::string smoother = prm.get<std::string>("smoother");
if (smoother == "ILU0") {
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
using Smoother = SeqILU<M, V, V>;
#else
@ -318,8 +318,8 @@ private:
}
});
doAddCreator("kamg", [](const O& op, const P& prm, const std::function<Vector()>&) {
const std::string smoother = prm.get<std::string>("smoother");
if (smoother == "ILU0") {
const std::string smoother = prm.get<std::string>("smoother", "ParOverILU0");
if (smoother == "ILU0" || smoother == "ParOverILU0") {
using Smoother = SeqILU0<M, V, V>;
return makeAmgPreconditioner<Smoother>(op, prm, true);
} else if (smoother == "Jac") {
@ -378,7 +378,7 @@ private:
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");
const std::string& type = prm.get<std::string>("type", "ParOverILU0");
auto it = creators_.find(type);
if (it == creators_.end()) {
std::ostringstream msg;
@ -395,7 +395,7 @@ private:
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");
const std::string& type = prm.get<std::string>("type", "ParOverILU0");
auto it = parallel_creators_.find(type);
if (it == parallel_creators_.end()) {
std::ostringstream msg;