diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh index 9c3fe581765d50e38c2ac64cf5c11fa44b9a52c1..8ffdb1b165c404b29e39ffdd13d180e177d49900 100644 --- a/dumux/linear/seqsolverbackend.hh +++ b/dumux/linear/seqsolverbackend.hh @@ -26,6 +26,7 @@ #include <dune/istl/solvers.hh> #include <dumux/common/propertysystem.hh> +#include <dune/istl/superlu.hh> namespace Dumux { @@ -119,6 +120,46 @@ private: Dune::InverseOperatorResult result_; }; +template <class TypeTag> +class ILU0CGBackend +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; +public: + + ILU0CGBackend(Problem& problem) + {} + + template<class Matrix, class Vector> + void solve(const Matrix& A, Vector& x, const Vector& b) + { + int verbosity = GET_PROP_VALUE(TypeTag, PTAG(LSVerbosity)); + static const int maxIter = GET_PROP_VALUE(TypeTag, PTAG(LSMaxIterations)); + static const double residReduction = GET_PROP_VALUE(TypeTag, PTAG(LSResidualReduction)); + + Vector bTmp(b); + + typedef Dune::SeqILU0<Matrix, Vector, Vector> Preconditioner; + static const double relaxation = GET_PROP_VALUE(TypeTag, PTAG(PreconditionerRelaxation)); + Preconditioner precond(A, relaxation); + + typedef Dune::MatrixAdapter<Matrix, Vector, Vector> MatrixAdapter; + MatrixAdapter operatorA(A); + + typedef Dune::CGSolver<Vector> Solver; + Solver solver(operatorA, precond, residReduction, maxIter, verbosity); + + solver.apply(x, bTmp, result_); + } + + const Dune::InverseOperatorResult& result() const + { + return result_; + } + +private: + Dune::InverseOperatorResult result_; +}; + template <class TypeTag> class IterativePrecondSolverBackend { @@ -258,5 +299,143 @@ public: } }; +template <class TypeTag> +class ILUnCGBackend: public IterativePrecondSolverBackend<TypeTag> +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + typedef IterativePrecondSolverBackend<TypeTag> ParentType; +public: + + ILUnCGBackend(Problem& problem) + {} + + template<class Matrix, class Vector> + void solve(const Matrix& A, Vector& x, const Vector& b) + { + typedef Dune::SeqILUn<Matrix, Vector, Vector> Preconditioner; + typedef Dune::CGSolver<Vector> Solver; + + ParentType::template solve<Preconditioner, Solver>(A, x, b); + } +}; + +template <class TypeTag> +class SORCGBackend: public IterativePrecondSolverBackend<TypeTag> +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + typedef IterativePrecondSolverBackend<TypeTag> ParentType; +public: + + SORCGBackend(Problem& problem) + {} + + template<class Matrix, class Vector> + void solve(const Matrix& A, Vector& x, const Vector& b) + { + typedef Dune::SeqSOR<Matrix, Vector, Vector> Preconditioner; + typedef Dune::CGSolver<Vector> Solver; + + ParentType::template solve<Preconditioner, Solver>(A, x, b); + } +}; + +template <class TypeTag> +class SSORCGBackend: public IterativePrecondSolverBackend<TypeTag> +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + typedef IterativePrecondSolverBackend<TypeTag> ParentType; +public: + + SSORCGBackend(Problem& problem) + {} + + template<class Matrix, class Vector> + void solve(const Matrix& A, Vector& x, const Vector& b) + { + typedef Dune::SeqSSOR<Matrix, Vector, Vector> Preconditioner; + typedef Dune::CGSolver<Vector> Solver; + + ParentType::template solve<Preconditioner, Solver>(A, x, b); + } +}; + +template <class TypeTag> +class GSCGBackend: public IterativePrecondSolverBackend<TypeTag> +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + typedef IterativePrecondSolverBackend<TypeTag> ParentType; +public: + + GSCGBackend(Problem& problem) + {} + + template<class Matrix, class Vector> + void solve(const Matrix& A, Vector& x, const Vector& b) + { + typedef Dune::SeqGS<Matrix, Vector, Vector> Preconditioner; + typedef Dune::CGSolver<Vector> Solver; + + ParentType::template solve<Preconditioner, Solver>(A, x, b); + } +}; + +template <class TypeTag> +class JacCGBackend: public IterativePrecondSolverBackend<TypeTag> +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + typedef IterativePrecondSolverBackend<TypeTag> ParentType; +public: + + JacCGBackend(Problem& problem) + {} + + template<class Matrix, class Vector> + void solve(const Matrix& A, Vector& x, const Vector& b) + { + typedef Dune::SeqJac<Matrix, Vector, Vector> Preconditioner; + typedef Dune::CGSolver<Vector> Solver; + + ParentType::template solve<Preconditioner, Solver>(A, x, b); + } +}; + +#if HAVE_SUPERLU +template <class TypeTag> +class SuperLUBackend +{ + typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem; + +public: + + SuperLUBackend(const Problem& problem) + : problem_(problem) + {} + + template<class Matrix, class Vector> + int solve(const Matrix& A, Vector& x, const Vector& b) + { + int verbosity = GET_PROP_VALUE(TypeTag, PTAG(LSVerbosity)); + static const int maxIter = GET_PROP_VALUE(TypeTag, PTAG(LSMaxIterations)); + static const double residReduction = GET_PROP_VALUE(TypeTag, PTAG(LSResidualReduction)); + + Vector bTmp(b); + Dune::SuperLU<Matrix> solver(A); + + solver.apply(x, bTmp, result_); + + return result_.converged; + } + + const Dune::InverseOperatorResult& result() const + { + return result_; + } + +private: + Dune::InverseOperatorResult result_; + const Problem& problem_; +}; +#endif + } #endif