From baa0ff6af7381b6ebf2e4c704b73e7fae2755d72 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Sun, 26 Jan 2020 09:48:04 +0100
Subject: [PATCH] Added a generic linear backend using the solver factories
from istl.
This requires a version newer than 2.7.0. The backend translates
the configuration parameters from Dumux CamelCase to the dune-istl
names. If a Dumux parameter is not found we fall back to using
the dune-istl parameter name.
We fall back to the parallel AMG backend if the DUNE version is
not sufficient.
---
dumux/common/loggingparametertree.hh | 51 ++++++
dumux/common/parameters.hh | 7 +
dumux/linear/CMakeLists.txt | 1 +
dumux/linear/genericistlbackend.hh | 261 +++++++++++++++++++++++++++
4 files changed, 320 insertions(+)
create mode 100644 dumux/linear/genericistlbackend.hh
diff --git a/dumux/common/loggingparametertree.hh b/dumux/common/loggingparametertree.hh
index 8117172a29..0dc766d7ef 100644
--- a/dumux/common/loggingparametertree.hh
+++ b/dumux/common/loggingparametertree.hh
@@ -66,6 +66,16 @@ public:
bool hasKey(const std::string& key) const
{ return params_.hasKey(key); }
+ /** \brief test for key (even in default parameters
+ *
+ * Tests whether given key exists.
+ *
+ * \param key key name
+ * \return true if key exists in structure, otherwise false
+ */
+ bool hasKeyOrDefaultKey(const std::string& key) const
+ { return params_.hasKey(key) || defaultParams_.hasKey(key); }
+
/** \brief test for key in group
*
* Tests whether given key exists in a group.
@@ -101,6 +111,47 @@ public:
return false;
}
+ /** \brief test for subgroup in group
+ *
+ * Tests whether given sub group exists in a group.
+ * Given a group this function starts to look from the back
+ * for dots. In G1.G2.G3 the function first looks if the key
+ * "G3.Key" exists, then "G2.Key", ...
+ *
+ * \param key key name
+ * \param groupPrefix the group prefix name
+ * \return a vector of fully qualified groups ordered by decresing relevance
+ */
+ std::vector<std::string> getSubGroups(const std::string& groupName,
+ std::string groupPrefix) const
+ {
+ std::vector<std::string> groupNames;
+ /*
+ if (groupPrefix == "" &&
+ (params_.hasSub(groupName) || defaultParams_.hasSub(groupName))
+ {
+ groupNames.push_back(groupName);
+ return groupNames;
+ }
+ */
+ auto compoundGroup = groupPrefix.empty()? groupName : groupPrefix + "." + groupName;
+ auto dot = groupPrefix.rfind(".");
+
+ while (dot != std::string::npos)
+ {
+ groupPrefix = groupPrefix.substr(0, dot);
+ compoundGroup = groupPrefix + "." + groupName;
+ if (params_.hasSub(compoundGroup) || defaultParams_.hasSub(compoundGroup))
+ groupNames.push_back(compoundGroup);
+ dot = groupPrefix.rfind(".");
+ }
+
+ if (params_.hasSub(groupName) || defaultParams_.hasSub(groupName))
+ groupNames.push_back(groupName);
+
+ return groupNames;
+ }
+
/** \brief print the hierarchical parameter tree to stream
*
* \param stream the output stream to print to
diff --git a/dumux/common/parameters.hh b/dumux/common/parameters.hh
index b01b9aa1ad..4ae00d23f8 100644
--- a/dumux/common/parameters.hh
+++ b/dumux/common/parameters.hh
@@ -454,6 +454,13 @@ inline bool hasParam(const std::string& param)
inline bool hasParamInGroup(const std::string& paramGroup, const std::string& param)
{ return Parameters::getTree().hasKeyInGroup(param, paramGroup); }
+/*!
+ * \ingroup Common
+ * \brief Get a list of sub groups from the parameter tree sorted by relevance
+ * \return A vector of fully qualified subGroup names sorted by descending relevance.
+ */
+inline std::vector<std::string> getSubGroups(const std::string& paramGroup, const std::string& subGroupName)
+{ return Parameters::getTree().getSubGroups(subGroupName, paramGroup); }
} // namespace Dumux
#endif
diff --git a/dumux/linear/CMakeLists.txt b/dumux/linear/CMakeLists.txt
index 6cdfe68ee2..14e0505b49 100644
--- a/dumux/linear/CMakeLists.txt
+++ b/dumux/linear/CMakeLists.txt
@@ -2,6 +2,7 @@ install(FILES
amgbackend.hh
amgparallelhelpers.hh
amgtraits.hh
+genericistlbackend.hh
linearsolvertraits.hh
linearsolveracceptsmultitypematrix.hh
matrixconverter.hh
diff --git a/dumux/linear/genericistlbackend.hh b/dumux/linear/genericistlbackend.hh
new file mode 100644
index 0000000000..75f10fc35b
--- /dev/null
+++ b/dumux/linear/genericistlbackend.hh
@@ -0,0 +1,261 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program 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. *
+ * *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Linear
+ * \brief Provides a generic linear solver based on the ISTL that chooses the
+ * solver and preconditioner at runtime. Needs dune version 2.7.1 or
+ * higher
+ */
+
+#ifndef DUMUX_GENERIC_ISTLBACKEND_HH
+#define DUMUX_GENERIC_ISTLBACKEND_HH
+
+#include <dumux/linear/solver.hh>
+#include <dumux/linear/amgtraits.hh>
+#include <dumux/linear/amgparallelhelpers.hh>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+// preconditioners
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/paamg/amg.hh>
+
+// solvers
+#include <dune/istl/solvers.hh>
+
+#if DUNE_VERSION_NEWER_REV(DUNE_ISTL,2,7,1)
+
+namespace Dumux
+{
+/*!
+ * \ingroup Linear
+ * \brief A linear solver using the ISTL solver factory
+ * allowing choosing the solver and preconditioner
+ * at runtime.
+ */
+template <class Matrix, class Vector, class GridGeometry>
+class GenericIstlBackend : public LinearSolver
+{
+ using GridView = typename GridGeometry::GridView;
+ using AMGTraits = AmgTraits<Matrix, Vector, GridGeometry>;
+ using Grid = typename GridView::Grid;
+ using LinearOperator = typename AMGTraits::LinearOperator;
+ using ScalarProduct = typename AMGTraits::ScalarProduct;
+ using VType = typename AMGTraits::VType;
+ using Comm = typename AMGTraits::Comm;
+ using BCRSMat = typename AMGTraits::LinearOperator::matrix_type;
+ using DofMapper = typename AMGTraits::DofMapper;
+
+public:
+ /*!
+ * \brief Construct the backend for the sequential case only
+ *
+ * \param paramGroup the parameter group for parameter lookup
+ */
+ GenericIstlBackend(const std::string& paramGroup = "")
+ : firstCall_(true)
+ {
+ if (Dune::MPIHelper::getCollectiveCommunication().size() > 1)
+ DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
+ convertParameterTree(paramGroup);
+ }
+ /*!
+ * \brief Construct the backend for parallel or sequential runs
+ *
+ * \param gridView the grid view on which we are performing the multi-grid
+ * \param dofMapper an index mapper for dof entities
+ * \param paramGroup the parameter group for parameter lookup
+ */
+ GenericIstlBackend(const GridView& gridView,
+ const DofMapper& dofMapper,
+ const std::string& paramGroup = "")
+ : phelper_(std::make_shared<ParallelISTLHelper<GridView, AMGTraits>>(gridView, dofMapper))
+ , firstCall_(true)
+ {
+ convertParameterTree(paramGroup);
+ }
+ /*!
+ * \brief Solve a linear system.
+ *
+ * \param A the matrix
+ * \param x the seeked solution vector, containing the initial solution upon entry
+ * \param b the right hand side vector
+ */
+ bool solve(Matrix& A, Vector& x, Vector& b)
+ {
+ int rank = 0;
+ std::shared_ptr<Comm> comm;
+ std::shared_ptr<LinearOperator> fop;
+ std::shared_ptr<ScalarProduct> sp; // not used.
+
+#if HAVE_MPI
+ if constexpr (AMGTraits::isParallel)
+ prepareLinearAlgebraParallel<AMGTraits>(A, b, comm, fop, sp, *phelper_, firstCall_);
+ else
+ prepareLinearAlgebraSequential<AMGTraits>(A, comm, fop, sp);
+#else
+ prepareLinearAlgebraSequential<AMGTraits>(A, comm, fop, sp);
+#endif
+
+ if (firstCall_)
+ {
+ Dune::initSolverFactories<typename AMGTraits::LinearOperator>();
+ }
+ try{
+ std::shared_ptr<Dune::InverseOperator<Vector, Vector>> solver = getSolverFromFactory(fop, params_);
+ Dune::InverseOperatorResult res;
+ solver->apply(x,b,result_);
+ }catch(Dune::Exception& e){
+ std::cerr << "Could not create solver" << std::endl;
+ std::cerr << e.what() << std::endl;
+ throw e;
+ }
+ firstCall_ = false;
+ return result_.converged;
+ }
+private:
+
+ void convertParameterTree(const std::string& paramGroup="")
+ {
+ const auto& loggingTree = Parameters::getTree();
+ auto matchingGroups = loggingTree.getSubGroups("LinearSolver",
+ paramGroup);
+
+ bool doThrow{};
+
+ for (const auto& transPair : istl2DumuxSolverParams)
+ {
+ for (const auto fullGroup : matchingGroups)
+ {
+ auto istlName = fullGroup+"."+transPair[0];
+ auto dumuxName = fullGroup+"."+transPair[1];
+ if(loggingTree.hasKeyOrDefaultKey(dumuxName))
+ {
+ if(loggingTree.hasKeyOrDefaultKey(istlName))
+ {
+ std::cerr << "Found equivalent keys " << istlName
+ << " " << dumuxName << std::endl
+ << "Please use only one (e.g. " << dumuxName
+ << ")." << std::endl;
+ doThrow = true;
+ }
+ params_[transPair[0]] = loggingTree.get<std::string>(dumuxName);
+ break;
+ }
+ else if (loggingTree.hasKey(istlName))
+ {
+ params_[transPair[0]] = loggingTree.get<std::string>(istlName);
+ break;
+ }
+ }
+ }
+
+ for (const auto& transPair : istl2DumuxPreconditionerParams)
+ {
+ for (const auto fullGroup : matchingGroups)
+ {
+ auto istlName = fullGroup+".preconditioner."+transPair[0];
+ auto dumuxName = fullGroup+"."+transPair[1];
+ if(loggingTree.hasKey(dumuxName))
+ {
+ if(loggingTree.hasKeyOrDefaultKey(istlName))
+ {
+ std::cerr << "Found equivalent keys " << istlName
+ << " " << dumuxName << std::endl
+ << "Please use only one (e.g. " << dumuxName
+ << ")." << std::endl;
+ doThrow = true;
+ }
+ params_["preconditioner."+transPair[0]] = loggingTree.get<std::string>(dumuxName);
+ break;
+ }
+ else if (loggingTree.hasKeyOrDefaultKey(istlName))
+ {
+ params_["preconditioner."+transPair[0]] = loggingTree.get<std::string>(istlName);
+ break;
+ }
+ }
+ }
+ params_.report();
+ if (!params_.hasKey("type"))
+ // prevent throw in solve
+ DUNE_THROW(Dune::InvalidStateException, "Solverfactory needs a specified \"type\" key to select the solver");
+
+ if (doThrow)
+ DUNE_THROW(Dune::InvalidStateException, "Ambiguous parameters used for linear solver");
+ }
+
+ static std::vector<std::array<std::string,2> > istl2DumuxSolverParams;
+ static std::vector<std::array<std::string,2> > istl2DumuxPreconditionerParams;
+ std::shared_ptr<ParallelISTLHelper<GridView, AMGTraits>> phelper_;
+ bool firstCall_;
+ Dune::InverseOperatorResult result_;
+ Dune::ParameterTree params_;
+};
+
+template<class Matrix, class Vector, class Geometry>
+std::vector<std::array<std::string,2> > GenericIstlBackend<Matrix, Vector, Geometry>::istl2DumuxSolverParams =
+ {
+ {"verbose", "Verbosity"}, {"maxit", "MaxIterations"},
+ {"reduction", "ResidualReduction"}, {"type", "Type"},
+ {"restart", "Restart"}, // cycles before restarting
+ // maximum number of vectors to store for orthogonalization
+ {"mmax", "MaxOrthogonalizationVectors"}
+ };
+
+template<class Matrix, class Vector, class Geometry>
+std::vector<std::array<std::string,2> > GenericIstlBackend<Matrix, Vector, Geometry>::istl2DumuxPreconditionerParams =
+ {
+ {"verbosity", "PreconditionerVerbosity"}, {"type", "PreconditionerType"},
+ {"iterations", "PreconditionerIterations"}, {"relaxation", "PreconditionerRelaxation"},
+ {"n", "ILUOrder"}, {"resort", "ILUResort"},
+ {"smootherRelaxation", "AmgSmootherRelaxation"},
+ {"smootherIterations", "AmgSmootherIterations"},
+ {"maxLevel", "AmgMaxLevel"}, {"coarsenTarget", "AmgCoarsenTarget"},
+ {"minCoarseningRate", "MinCoarseningRate"},
+ {"prolongationDampingFactor", "AmgProlongationDampingFactor"},
+ {"alpha", "AmgAlpha"}, {"beta", "AmgBeta"},
+ {"additive", "AmgAdditive"}, {"gamma", "AmgGamma"},
+ {"preSteps", "AmgPreSmoothingSteps"}, {"postSteps", "AmgPostSmoothingSteps"},
+ {"criterionSymmetric", "AmgCriterionSymmetric"}, {"strengthMeasure", "AmgStrengthMeasure"},
+ {"diagonalRowIndex", "AmgDiagonalRowIndex"},
+ {"defaultAggregationSizeMode", "DefaultAggregationSizeMode"},
+ {"defaultAggregationDimension", "defaultAggregationDimension"},
+ {"maxAggregateDistance", "MaxAggregateDistance"},
+ {"minAggregateSize", "MinAggregateSize"},
+ {"maxAggregateSize", "MaxAggregateSize"}
+ };
+
+template<class TypeTag>
+using IstlSolverBackend = GenericIstlBackend<
+ GetPropType<TypeTag, Properties::JacobianMatrix>,
+ GetPropType<TypeTag, Properties::SolutionVector>,
+ GetPropType<TypeTag, Properties::GridGeometry>>;
+} // end namespace Dumux
+#else
+#warn "Generic ISTL backend needs dune-istl > 2.7.0!"
+#warn "Ignoring configuration parameters and falling back to AMG backend."
+#include "amgbackend.hh"
+template<class TypeTag>
+using IstlSolverBackend = AMGBackend<TypeTag>;
+#endif // DUNE version check
+#endif // DUMUX_GENERIC_ISTLBACKEND_HH
--
GitLab