diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh index 1ada67218fff0a76cb3444fd2d8699c5ee673532..e14ca263852ce64196667d2619bf512a8db57fe7 100644 --- a/dumux/linear/seqsolverbackend.hh +++ b/dumux/linear/seqsolverbackend.hh @@ -33,6 +33,7 @@ #include <dumux/common/parameters.hh> #include <dumux/common/basicproperties.hh> #include <dumux/linear/linearsolverproperties.hh> +#include <dumux/linear/solver.hh> namespace Dumux { @@ -61,24 +62,17 @@ class IterativePreconditionedSolverImpl { public: - template<class Preconditioner, class Solver, class Matrix, class Vector> - static bool solve(const Matrix& A, Vector& x, const Vector& b, + template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector> + static bool solve(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b, const std::string& modelParamGroup = "") { - // get some runtime configurable parameters - const int verbosity = getParamFromGroup<int>(modelParamGroup, "LinearSolver.Verbosity"); - const int maxIter = getParamFromGroup<int>(modelParamGroup, "LinearSolver.MaxIterations"); - const double residReduction = getParamFromGroup<double>(modelParamGroup, "LinearSolver.ResidualReduction"); - const double relaxation = getParamFromGroup<double>(modelParamGroup, "LinearSolver.PreconditionerRelaxation"); - const int precondIter = getParamFromGroup<int>(modelParamGroup, "LinearSolver.PreconditionerIterations"); - - Preconditioner precond(A, precondIter, relaxation); + Preconditioner precond(A, s.precondIter(), s.relaxation()); // make a linear operator from a matrix using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>; MatrixAdapter linearOperator(A); - Solver solver(linearOperator, precond, residReduction, maxIter, verbosity); + Solver solver(linearOperator, precond, s.residReduction(), s.maxIter(), s.verbosity()); Vector bTmp(b); @@ -88,25 +82,20 @@ public: return result.converged; } - template<class Preconditioner, class Solver, class Matrix, class Vector> - static bool solveWithGMRes(const Matrix& A, Vector& x, const Vector& b, + template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector> + static bool solveWithGMRes(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b, const std::string& modelParamGroup = "") { - // get some runtime configurable parameters - const int verbosity = getParamFromGroup<int>(modelParamGroup, "LinearSolver.Verbosity"); - const int maxIter = getParamFromGroup<int>(modelParamGroup, "LinearSolver.MaxIterations"); - const double residReduction = getParamFromGroup<double>(modelParamGroup, "LinearSolver.ResidualReduction"); - const double relaxation = getParamFromGroup<double>(modelParamGroup, "LinearSolver.PreconditionerRelaxation"); - const int precondIter = getParamFromGroup<int>(modelParamGroup, "LinearSolver.PreconditionerIterations"); + // get the restart threshold const int restartGMRes = getParamFromGroup<double>(modelParamGroup, "LinearSolver.GMResRestart"); - Preconditioner precond(A, precondIter, relaxation); + Preconditioner precond(A, s.precondIter(), s.relaxation()); // make a linear operator from a matrix using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>; MatrixAdapter linearOperator(A); - Solver solver(linearOperator, precond, residReduction, restartGMRes, maxIter, verbosity); + Solver solver(linearOperator, precond, s.residReduction(), restartGMRes, s.maxIter(), s.verbosity()); Vector bTmp(b); @@ -116,22 +105,16 @@ public: return result.converged; } - template<class Preconditioner, class Solver, class Matrix, class Vector> - static bool solveWithILU0Prec(const Matrix& A, Vector& x, const Vector& b, + template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector> + static bool solveWithILU0Prec(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b, const std::string& modelParamGroup = "") { - // get some runtime configurable parameters - const int verbosity = getParamFromGroup<int>(modelParamGroup, "LinearSolver.Verbosity"); - const int maxIter = getParamFromGroup<int>(modelParamGroup, "LinearSolver.MaxIterations"); - const double residReduction = getParamFromGroup<double>(modelParamGroup, "LinearSolver.ResidualReduction"); - const double relaxation = getParamFromGroup<double>(modelParamGroup, "LinearSolver.PreconditionerRelaxation"); - - Preconditioner precond(A, relaxation); + Preconditioner precond(A, s.relaxation()); using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>; MatrixAdapter operatorA(A); - Solver solver(operatorA, precond, residReduction, maxIter, verbosity); + Solver solver(operatorA, precond, s.residReduction(), s.maxIter(), s.verbosity()); Vector bTmp(b); @@ -142,23 +125,19 @@ public: } // solve with RestartedGMRes (needs restartGMRes as additional argument) - template<class Preconditioner, class Solver, class Matrix, class Vector> - static bool solveWithILU0PrecGMRes(const Matrix& A, Vector& x, const Vector& b, + template<class Preconditioner, class Solver, class SolverInterface, class Matrix, class Vector> + static bool solveWithILU0PrecGMRes(const SolverInterface& s, const Matrix& A, Vector& x, const Vector& b, const std::string& modelParamGroup = "") { - // get some runtime configurable parameters - const int verbosity = getParamFromGroup<int>(modelParamGroup, "LinearSolver.Verbosity"); - const int maxIter = getParamFromGroup<int>(modelParamGroup, "LinearSolver.MaxIterations"); - const double residReduction = getParamFromGroup<double>(modelParamGroup, "LinearSolver.ResidualReduction"); - const double relaxation = getParamFromGroup<double>(modelParamGroup, "LinearSolver.PreconditionerRelaxation"); + // get the restart threshold const int restartGMRes = getParamFromGroup<int>(modelParamGroup, "LinearSolver.GMResRestart"); - Preconditioner precond(A, relaxation); + Preconditioner precond(A, s.relaxation()); using MatrixAdapter = Dune::MatrixAdapter<Matrix, Vector, Vector>; MatrixAdapter operatorA(A); - Solver solver(operatorA, precond, residReduction, restartGMRes, maxIter, verbosity); + Solver solver(operatorA, precond, s.residReduction(), restartGMRes, s.maxIter(), s.verbosity()); Vector bTmp(b); @@ -187,7 +166,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class ILUnBiCGSTABBackend +class ILUnBiCGSTABBackend : public LinearSolver<TypeTag> { public: @@ -199,7 +178,7 @@ public: using Preconditioner = Dune::SeqILUn<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::BiCGSTABSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -226,7 +205,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class SORBiCGSTABBackend +class SORBiCGSTABBackend : public LinearSolver<TypeTag> { public: @@ -238,7 +217,7 @@ public: using Preconditioner = Dune::SeqSOR<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::BiCGSTABSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -265,7 +244,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class SSORBiCGSTABBackend +class SSORBiCGSTABBackend : public LinearSolver<TypeTag> { public: @@ -277,7 +256,7 @@ public: using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::BiCGSTABSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -304,7 +283,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class GSBiCGSTABBackend +class GSBiCGSTABBackend : public LinearSolver<TypeTag> { public: @@ -316,7 +295,7 @@ public: using Preconditioner = Dune::SeqGS<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::BiCGSTABSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -342,7 +321,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class JacBiCGSTABBackend +class JacBiCGSTABBackend : public LinearSolver<TypeTag> { public: @@ -354,7 +333,7 @@ public: using Preconditioner = Dune::SeqJac<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::BiCGSTABSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -380,7 +359,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class ILUnCGBackend +class ILUnCGBackend : public LinearSolver<TypeTag> { public: @@ -392,7 +371,7 @@ public: using Preconditioner = Dune::SeqILUn<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::CGSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -418,7 +397,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class SORCGBackend +class SORCGBackend : public LinearSolver<TypeTag> { public: @@ -430,7 +409,7 @@ public: using Preconditioner = Dune::SeqSOR<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::CGSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -456,7 +435,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class SSORCGBackend +class SSORCGBackend : public LinearSolver<TypeTag> { public: @@ -468,7 +447,7 @@ public: using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::CGSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -494,7 +473,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class GSCGBackend +class GSCGBackend : public LinearSolver<TypeTag> { public: @@ -506,7 +485,7 @@ public: using Preconditioner = Dune::SeqGS<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::CGSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -531,7 +510,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class JacCGBackend +class JacCGBackend : public LinearSolver<TypeTag> { public: @@ -543,7 +522,7 @@ public: using Preconditioner = Dune::SeqJac<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::CGSolver<Vector>; - return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solve<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -570,7 +549,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class SSORRestartedGMResBackend +class SSORRestartedGMResBackend : public LinearSolver<TypeTag> { public: @@ -582,7 +561,7 @@ public: using Preconditioner = Dune::SeqSSOR<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::RestartedGMResSolver<Vector>; - return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -608,7 +587,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class ILU0BiCGSTABBackend +class ILU0BiCGSTABBackend : public LinearSolver<TypeTag> { public: @@ -620,7 +599,7 @@ public: using Preconditioner = Dune::SeqILU0<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::BiCGSTABSolver<Vector>; - return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -645,7 +624,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class ILU0CGBackend +class ILU0CGBackend : public LinearSolver<TypeTag> { public: @@ -657,7 +636,7 @@ public: using Preconditioner = Dune::SeqILU0<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::CGSolver<Vector>; - return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solveWithILU0Prec<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -683,7 +662,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class ILU0RestartedGMResBackend +class ILU0RestartedGMResBackend : public LinearSolver<TypeTag> { public: @@ -695,7 +674,7 @@ public: using Preconditioner = Dune::SeqILU0<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::RestartedGMResSolver<Vector>; - return IterativePreconditionedSolverImpl::template solveWithILU0PrecGMRes<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solveWithILU0PrecGMRes<Preconditioner, Solver>(*this, A, x, b, paramGroup); } std::string name() const @@ -722,7 +701,7 @@ public: * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press. */ template <class TypeTag> -class ILUnRestartedGMResBackend +class ILUnRestartedGMResBackend : public LinearSolver<TypeTag> { public: @@ -734,7 +713,7 @@ public: using Preconditioner = Dune::SeqILUn<Matrix, Vector, Vector, blockLevel>; using Solver = Dune::RestartedGMResSolver<Vector>; - return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(A, x, b, paramGroup); + return IterativePreconditionedSolverImpl::template solveWithGMRes<Preconditioner, Solver>(*this, A, x, b, paramGroup); } }; @@ -748,7 +727,7 @@ public: * http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ */ template <class TypeTag> -class SuperLUBackend +class SuperLUBackend : public LinearSolver<TypeTag> { public: @@ -761,8 +740,7 @@ public: using ISTLMatrix = typename Dune::BCRSMatrix<MatrixBlock>; static_assert(std::is_same<Matrix, ISTLMatrix>::value, "SuperLU only works with BCRS matrices!"); - const int verbosity = getParamFromGroup<int>(paramGroup, "LinearSolver.Verbosity"); - Dune::SuperLU<ISTLMatrix> solver(A, verbosity > 0); + Dune::SuperLU<ISTLMatrix> solver(A, this->verbosity() > 0); Vector bTmp(b); solver.apply(x, bTmp, result_); @@ -810,7 +788,7 @@ private: * http://faculty.cse.tamu.edu/davis/suitesparse.html */ template <class TypeTag> -class UMFPackBackend +class UMFPackBackend : public LinearSolver<TypeTag> { public: @@ -823,8 +801,7 @@ public: using ISTLMatrix = typename Dune::BCRSMatrix<MatrixBlock>; static_assert(std::is_same<Matrix, ISTLMatrix>::value, "UMFPack only works with BCRS matrices!"); - const int verbosity = getParamFromGroup<int>(paramGroup, "LinearSolver.Verbosity"); - Dune::UMFPack<ISTLMatrix> solver(A, verbosity > 0); + Dune::UMFPack<ISTLMatrix> solver(A, this->verbosity() > 0); Vector bTmp(b); solver.apply(x, bTmp, result_); diff --git a/dumux/linear/solver.hh b/dumux/linear/solver.hh new file mode 100644 index 0000000000000000000000000000000000000000..1c74fb10b7771d99bd4e63e75fac5b07f751e876 --- /dev/null +++ b/dumux/linear/solver.hh @@ -0,0 +1,114 @@ +// -*- 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 2 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 + * \brief Dumux sequential solver backends + */ +#ifndef DUMUX_LINEAR_SOLVER_HH +#define DUMUX_LINEAR_SOLVER_HH + +#include <dune/common/exceptions.hh> +#include <dumux/common/parameters.hh> +#include <dumux/common/basicproperties.hh> +#include <dumux/linear/linearsolverproperties.hh> + +namespace Dumux +{ + +/*! + * \ingroup Linear + * \brief Base class for linear solvers + */ +template <class TypeTag> +class LinearSolver +{ + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); +public: + //! default constructor sets some parameters + LinearSolver() + { + static const std::string modelParamGroup = GET_PROP_VALUE(TypeTag, ModelParameterGroup); + verbosity_ = getParamFromGroup<int>(modelParamGroup, "LinearSolver.Verbosity"); + maxIter_ = getParamFromGroup<int>(modelParamGroup, "LinearSolver.MaxIterations"); + residReduction_ = getParamFromGroup<Scalar>(modelParamGroup, "LinearSolver.ResidualReduction"); + relaxation_ = getParamFromGroup<Scalar>(modelParamGroup, "LinearSolver.PreconditionerRelaxation"); + precondIter_ = getParamFromGroup<int>(modelParamGroup, "LinearSolver.PreconditionerIterations"); + } + + template<class Matrix, class Vector> + bool solve(const Matrix& A, Vector& x, const Vector& b) + { + DUNE_THROW(Dune::NotImplemented, "Linear solver doesn't implement a solve method!"); + } + + //! the name of the linear solver + std::string name() const + { return "unknown solver"; } + + //! the verbosity level + int verbosity() const + { return verbosity_; } + + //! set the verbosity level + void setVerbosity(int v) + { verbosity_ = v; } + + //! the maximum number of linear solver iterations + int maxIter() const + { return maxIter_; } + + //! set the maximum number of linear solver iterations + void setMaxIter(int i) + { maxIter_ = i; } + + //! the linear solver residual reduction + Scalar residReduction() const + { return residReduction_; } + + //! set the linear solver residual reduction + void setResidualReduction(Scalar r) + { residReduction_ = r; } + + //! the linear solver relaxation factor + Scalar relaxation() const + { return relaxation_; } + + //! set the linear solver relaxation factor + void setRelaxation(Scalar r) + { relaxation_ = r; } + + //! the number of preconditioner iterations + int precondIter() const + { return precondIter_; } + + //! set the number of preconditioner iterations + void setPrecondIter(int i) + { precondIter_ = i; } + +private: + int verbosity_; + int maxIter_; + Scalar residReduction_; + Scalar relaxation_; + int precondIter_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh index 0cb7e4a77ba93a92323a168e77596b9c112e6b68..b0be643d2d6aebfd7a9d3f0d9b23c30b8f72ce0c 100644 --- a/dumux/nonlinear/newtonmethod.hh +++ b/dumux/nonlinear/newtonmethod.hh @@ -55,6 +55,7 @@ NEW_PROP_TAG(JacobianMatrix); template <class TypeTag, class NewtonController, class JacobianAssembler, class LinearSolver> class NewtonMethod { + using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector); using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix); @@ -70,6 +71,11 @@ public: { // set the linear system (matrix & residual) in the assembler assembler_->setLinearSystem(matrix_, residual_); + + // set a different default for the linear solver residual reduction + // within the Newton the linear solver doesn't need to solve too exact + static const std::string modelParamGroup = GET_PROP_VALUE(TypeTag, ModelParameterGroup); + linearSolver_->setResidualReduction(getParamFromGroup<Scalar>(modelParamGroup, "LinearSolver.ResidualReduction", 1e-6)); } /*!