From db914f1c783a3f1595908114a00de0db04aa9306 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 12 Dec 2012 17:13:58 +0000
Subject: [PATCH] implicit branch: copy AMGBackend from devel. Copy and modify
the necessary PDELab backend files.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9821 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/linear/amgbackend.hh | 246 +++++
dumux/linear/istlvectorbackend.hh | 262 +++++
dumux/linear/novlpistlsolverbackend.hh | 1346 ++++++++++++++++++++++++
dumux/linear/ovlpistlsolverbackend.hh | 1000 ++++++++++++++++++
dumux/linear/seqistlsolverbackend.hh | 743 +++++++++++++
5 files changed, 3597 insertions(+)
create mode 100644 dumux/linear/amgbackend.hh
create mode 100644 dumux/linear/istlvectorbackend.hh
create mode 100644 dumux/linear/novlpistlsolverbackend.hh
create mode 100644 dumux/linear/ovlpistlsolverbackend.hh
create mode 100644 dumux/linear/seqistlsolverbackend.hh
diff --git a/dumux/linear/amgbackend.hh b/dumux/linear/amgbackend.hh
new file mode 100644
index 0000000000..aebb85536f
--- /dev/null
+++ b/dumux/linear/amgbackend.hh
@@ -0,0 +1,246 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Andreas Lauser *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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 Provides a linear solver for the stabilized BiCG method with
+ * an ILU-0 preconditioner.
+ */
+#ifndef DUMUX_AMGBACKEND_HH
+#define DUMUX_AMGBACKEND_HH
+
+#include <dumux/linear/boxlinearsolver.hh>
+#include <dumux/boxmodels/pdelab/pdelabadapter.hh>
+#include <dumux/boxmodels/pdelab/pdelabboxistlvectorbackend.hh>
+#include <dune/pdelab/backend/novlpistlsolverbackend.hh>
+#include <dune/pdelab/backend/ovlpistlsolverbackend.hh>
+#include <dune/pdelab/backend/seqistlsolverbackend.hh>
+
+namespace Dumux {
+
+namespace Properties
+{
+NEW_PROP_TAG(GridOperator);
+}
+
+template <class Matrix, class Vector>
+void scaleLinearSystem(Matrix& matrix, Vector& rhs)
+{
+ typename Matrix::RowIterator row = matrix.begin();
+ for(; row != matrix.end(); ++row)
+ {
+ typedef typename Matrix::size_type size_type;
+ size_type rowIdx = row.index();
+
+ typedef typename Matrix::block_type MatrixBlock;
+ MatrixBlock diagonal = matrix[rowIdx][rowIdx];
+ diagonal.invert();
+
+ typedef typename Vector::block_type VectorBlock;
+ const VectorBlock b = rhs[rowIdx];
+ diagonal.mv(b, rhs[rowIdx]);
+
+ typename Matrix::ColIterator col = row->begin();
+ for (; col != row->end(); ++col)
+ col->leftmultiply(diagonal);
+ }
+}
+
+template <class TypeTag>
+class AMGBackend
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridOperator) GridOperator;
+ typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_AMG_SSOR<GridOperator> Imp;
+public:
+
+ AMGBackend(const Problem& problem)
+ : problem_(problem)
+ {}
+
+ template<class Matrix, class Vector>
+ bool solve(Matrix& A, Vector& x, Vector& b)
+ {
+ imp_ = new Imp(problem_.gridFunctionSpace(),
+ GET_PROP_VALUE(TypeTag, LinearSolverMaxIterations),
+ GET_PROP_VALUE(TypeTag, LinearSolverVerbosity));
+
+ static const double residReduction = GET_PROP_VALUE(TypeTag, LinearSolverResidualReduction);
+ imp_->apply(A, x, b, residReduction);
+
+ result_.converged = imp_->result().converged;
+ result_.iterations = imp_->result().iterations;
+ result_.elapsed = imp_->result().elapsed;
+ result_.reduction = imp_->result().reduction;
+ result_.conv_rate = imp_->result().conv_rate;
+
+ delete imp_;
+
+ return result_.converged;
+ }
+
+ const Dune::InverseOperatorResult& result() const
+ {
+ return result_;
+ }
+
+private:
+ const Problem& problem_;
+ Imp *imp_;
+ Dune::InverseOperatorResult result_;
+};
+
+template <class TypeTag>
+class CCAMGBackend
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridOperator) GridOperator;
+ typedef Dune::PDELab::ISTLBackend_BCGS_AMG_SSOR<GridOperator> Imp;
+public:
+
+ CCAMGBackend(const Problem& problem)
+ : problem_(problem)
+ {}
+
+ template<class Matrix, class Vector>
+ bool solve(Matrix& A, Vector& x, Vector& b)
+ {
+ imp_ = new Imp(problem_.gridFunctionSpace(),
+ GET_PROP_VALUE(TypeTag, LinearSolverMaxIterations),
+ GET_PROP_VALUE(TypeTag, LinearSolverVerbosity));
+
+ static const double residReduction = GET_PROP_VALUE(TypeTag, LinearSolverResidualReduction);
+ imp_->apply(A, x, b, residReduction);
+
+ result_.converged = imp_->result().converged;
+ result_.iterations = imp_->result().iterations;
+ result_.elapsed = imp_->result().elapsed;
+ result_.reduction = imp_->result().reduction;
+ result_.conv_rate = imp_->result().conv_rate;
+
+ delete imp_;
+
+ return result_.converged;
+ }
+
+ const Dune::InverseOperatorResult& result() const
+ {
+ return result_;
+ }
+
+private:
+ const Problem& problem_;
+ Imp *imp_;
+ Dune::InverseOperatorResult result_;
+};
+
+template <class TypeTag>
+class SeqAMGBackend
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridOperator) GridOperator;
+ typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SSOR<GridOperator> Imp;
+public:
+
+ SeqAMGBackend(const Problem& problem)
+ : problem_(problem)
+ {}
+
+ template<class Matrix, class Vector>
+ bool solve(Matrix& A, Vector& x, Vector& b)
+ {
+ imp_ = new Imp(
+ GET_PROP_VALUE(TypeTag, LinearSolverMaxIterations),
+ GET_PROP_VALUE(TypeTag, LinearSolverVerbosity));
+
+ static const double residReduction = GET_PROP_VALUE(TypeTag, LinearSolverResidualReduction);
+ imp_->apply(A, x, b, residReduction);
+
+ result_.converged = imp_->result().converged;
+ result_.iterations = imp_->result().iterations;
+ result_.elapsed = imp_->result().elapsed;
+ result_.reduction = imp_->result().reduction;
+ result_.conv_rate = imp_->result().conv_rate;
+
+ delete imp_;
+
+ return result_.converged;
+ }
+
+ const Dune::InverseOperatorResult& result() const
+ {
+ return result_;
+ }
+
+private:
+ const Problem& problem_;
+ Imp *imp_;
+ Dune::InverseOperatorResult result_;
+};
+
+template <class TypeTag>
+class ScaledSeqAMGBackend
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridOperator) GridOperator;
+ typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_AMG_SSOR<GridOperator> Imp;
+public:
+
+ ScaledSeqAMGBackend(const Problem& problem)
+ : problem_(problem)
+ {}
+
+ template<class Matrix, class Vector>
+ bool solve(Matrix& A, Vector& x, Vector& b)
+ {
+ scaleLinearSystem(A, b);
+
+ imp_ = new Imp(
+ GET_PROP_VALUE(TypeTag, LinearSolverMaxIterations),
+ GET_PROP_VALUE(TypeTag, LinearSolverVerbosity));
+
+ static const double residReduction = GET_PROP_VALUE(TypeTag, LinearSolverResidualReduction);
+ imp_->apply(A, x, b, residReduction);
+
+ result_.converged = imp_->result().converged;
+ result_.iterations = imp_->result().iterations;
+ result_.elapsed = imp_->result().elapsed;
+ result_.reduction = imp_->result().reduction;
+ result_.conv_rate = imp_->result().conv_rate;
+
+ delete imp_;
+
+ return result_.converged;
+ }
+
+ const Dune::InverseOperatorResult& result() const
+ {
+ return result_;
+ }
+
+private:
+ const Problem& problem_;
+ Imp *imp_;
+ Dune::InverseOperatorResult result_;
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/linear/istlvectorbackend.hh b/dumux/linear/istlvectorbackend.hh
new file mode 100644
index 0000000000..cdbdde0705
--- /dev/null
+++ b/dumux/linear/istlvectorbackend.hh
@@ -0,0 +1,262 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+/*****************************************************************************
+ * This file is based on the file of DUNE-PDELab with the same name, *
+ * http://www.dune-project.org/pdelab/. *
+ * Some changes are made such that the backends work for the vectors and *
+ * matrices used in Dumux. *
+ *****************************************************************************/
+#ifndef DUNE_ISTLVECTORBACKEND_HH
+#define DUNE_ISTLVECTORBACKEND_HH
+
+#include<vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/pdelab/backend/backendselector.hh>
+#include <dune/pdelab/backend/istlmatrixbackend.hh>
+
+namespace Dumux {
+
+ template<int> class ISTLVectorBackend;
+
+ template<typename T, typename E, int BLOCKSIZE>
+ class ISTLBlockVectorContainer
+ {
+ public:
+ typedef E ElementType;
+ typedef Dune::BlockVector< Dune::FieldVector<E,BLOCKSIZE> > ContainerType;
+ typedef ContainerType BaseT;
+ typedef typename ContainerType::field_type field_type;
+ typedef typename ContainerType::iterator iterator;
+ typedef typename ContainerType::const_iterator const_iterator;
+ typedef typename ContainerType::block_type block_type;
+ typedef typename ContainerType::size_type size_type;
+ typedef ISTLVectorBackend<BLOCKSIZE> Backend;
+
+ ISTLBlockVectorContainer ()
+ {}
+ ISTLBlockVectorContainer (const T& t_) : container(t_.globalSize()/BLOCKSIZE)
+ {}
+ ISTLBlockVectorContainer (const T& t_, const E& e) : container(t_.globalSize()/BLOCKSIZE)
+ {
+ container=e;
+ }
+
+ size_type N() const
+ {
+ return container.N();
+ }
+
+
+ ISTLBlockVectorContainer& operator= (const E& e)
+ {
+ container=e;
+ return *this;
+ }
+
+ ISTLBlockVectorContainer& operator*= (const E& e)
+ {
+ container*=e;
+ return *this;
+ }
+
+
+ ISTLBlockVectorContainer& operator+= (const E& e)
+ {
+ container+=e;
+ return *this;
+ }
+
+ ISTLBlockVectorContainer& operator+= (const ISTLBlockVectorContainer& e)
+ {
+ container+=e;
+ return *this;
+ }
+
+ ISTLBlockVectorContainer& operator-= (const ISTLBlockVectorContainer& e)
+ {
+ container-=e;
+ return *this;
+ }
+
+ block_type& operator[](std::size_t i)
+ {
+ return container[i];
+ }
+
+ const block_type& operator[](std::size_t i) const
+ {
+ return container[i];
+ }
+
+ typename Dune::template FieldTraits<E>::real_type two_norm() const
+ {
+ return container.two_norm();
+ }
+
+ typename Dune::template FieldTraits<E>::real_type two_norm2() const
+ {
+ return container.two_norm2();
+ }
+
+ typename Dune::template FieldTraits<E>::real_type one_norm() const
+ {
+ return container.one_norm();
+ }
+
+ typename Dune::template FieldTraits<E>::real_type infinity_norm() const
+ {
+ return container.infinity_norm();
+ }
+
+ E operator*(const ISTLBlockVectorContainer& y) const
+ {
+ return container*y.base();
+ }
+
+ E dot(const ISTLBlockVectorContainer& y) const
+ {
+ return container.dot(y.base());
+ }
+
+ ISTLBlockVectorContainer& axpy(const E& a, const ISTLBlockVectorContainer& y)
+ {
+ container.axpy(a, y);
+ return *this;
+ }
+
+ // for debugging and AMG access
+ ContainerType& base ()
+ {
+ return container;
+ }
+
+ const ContainerType& base () const
+ {
+ return container;
+ }
+
+ operator ContainerType&()
+ {
+ return container;
+ }
+
+ operator const ContainerType&() const
+ {
+ return container;
+ }
+
+ iterator begin()
+ {
+ return container.begin();
+ }
+
+
+ const_iterator begin() const
+ {
+ return container.begin();
+ }
+
+ iterator end()
+ {
+ return container.end();
+ }
+
+
+ const_iterator end() const
+ {
+ return container.end();
+ }
+
+ size_t flatsize() const
+ {
+ return container.size()*BLOCKSIZE;
+ }
+
+ size_t size() const
+ {
+ return container.size();
+ }
+
+ void resize(size_t n)
+ {
+ container.resize(n);
+ }
+
+ template<typename X>
+ void std_copy_to (std::vector<X>& x) const
+ {
+ size_t n = flatsize();
+ x.resize(n);
+ for (size_t i=0; i<n; i++)
+ x[i] = container[i/BLOCKSIZE][i%BLOCKSIZE];
+ }
+
+ template<typename X>
+ void std_copy_from (const std::vector<X>& x)
+ {
+ //test if x has the same size as the container
+ assert (x.size() == flatsize());
+ for (size_t i=0; i<flatsize(); i++)
+ container[i/BLOCKSIZE][i%BLOCKSIZE] = x[i];
+ }
+
+ private:
+ Dune::BlockVector< Dune::FieldVector<E,BLOCKSIZE> > container;
+ };
+
+
+ //! ISTL backend for FunctionSpace
+ template<int BLOCKSIZE=1>
+ class ISTLVectorBackend
+ {
+ public:
+ enum{
+ //! \brief export the block size
+ BlockSize = BLOCKSIZE
+ };
+
+ //export Matrix Backend Type
+ typedef ISTLBCRSMatrixBackend<BLOCKSIZE,BLOCKSIZE> MatrixBackend;
+
+ //! container construction
+
+ // extract type of container element
+ template<class C>
+ struct Value
+ {
+ typedef typename C::field_type Type;
+ };
+
+ //! The size type
+ typedef typename Dune::BlockVector< Dune::FieldVector<float,BLOCKSIZE> >::size_type size_type;
+
+ // get const_reference to container element
+ // note: this method does not depend on T!
+ template<typename C>
+ static const typename C::field_type& access (const C& c, size_type i)
+ {
+ return c.base()[i/BLOCKSIZE][i%BLOCKSIZE];
+ }
+
+ // get non const_reference to container element
+ // note: this method does not depend on T!
+ template<typename C>
+ static typename C::field_type& access (C& c, size_type i)
+ {
+ return c.base()[i/BLOCKSIZE][i%BLOCKSIZE];
+ }
+ };
+
+ template<int BLOCKSIZE,typename T, typename E>
+ struct BackendVectorSelectorHelper<ISTLVectorBackend<BLOCKSIZE>, T, E>
+ {
+ typedef ISTLBlockVectorContainer<T,E,BLOCKSIZE> Type;
+ };
+
+
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/linear/novlpistlsolverbackend.hh b/dumux/linear/novlpistlsolverbackend.hh
new file mode 100644
index 0000000000..da5f1d4bdd
--- /dev/null
+++ b/dumux/linear/novlpistlsolverbackend.hh
@@ -0,0 +1,1346 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=8 sw=2 sts=2:
+/*****************************************************************************
+ * This file is based on the file of DUNE-PDELab with the same name, *
+ * http://www.dune-project.org/pdelab/. *
+ * Some changes are made such that the backends work for the vectors and *
+ * matrices used in Dumux. *
+ *****************************************************************************/
+#ifndef DUNE_NOVLPISTLSOLVERBACKEND_HH
+#define DUNE_NOVLPISTLSOLVERBACKEND_HH
+
+#include <cstddef>
+
+#include <dune/common/deprecated.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <dune/grid/common/gridenums.hh>
+
+#include <dune/istl/io.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/pinfo.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/solvercategory.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/superlu.hh>
+
+#include <dune/pdelab/backend/parallelistlhelper.hh>
+
+#include "istlvectorbackend.hh"
+#include "seqistlsolverbackend.hh"
+
+namespace Dumux {
+
+ //! \addtogroup Backend
+ //! \ingroup PDELab
+ //! \{
+
+ //========================================================
+ // Generic support for nonoverlapping grids
+ //========================================================
+
+ //! Operator for the non-overlapping parallel case
+ /**
+ * Calculate \f$y:=Ax\f$.
+ *
+ * \tparam GFS The GridFunctionSpace the vectors apply to.
+ * \tparam M Type of the matrix. Should be one of the ISTL matrix types.
+ * \tparam X Type of the vectors the matrix is applied to.
+ * \tparam Y Type of the result vectors.
+ */
+ template<class GFS, class M, class X, class Y>
+ class NonoverlappingOperator : public Dune::AssembledLinearOperator<M,X,Y>
+ {
+ public:
+ //! export type of matrix
+ typedef M matrix_type;
+ //! export type of vectors the matrix is applied to
+ typedef X domain_type;
+ //! export type of result vectors
+ typedef Y range_type;
+ //! export type of the entries for x
+ typedef typename X::field_type field_type;
+
+ //redefine the category, that is the only difference
+ enum {category=Dune::SolverCategory::nonoverlapping};
+
+ //! Construct a non-overlapping operator
+ /**
+ * \param gfs_ GridFunctionsSpace for the vectors.
+ * \param A Matrix for this operator. This should be the locally
+ * assembled matrix.
+ * \param helper_ Helper for parallel communication (not used).
+ *
+ * \note The constructed object stores references to all the objects
+ * given as parameters here. They should be valid for as long as
+ * the constructed object is used. They are not needed to
+ * destruct the constructed object.
+ *
+ * \deprecated The helper_ parameter is unused. Use the constructor
+ * without the helper_ parameter instead.
+ */
+ NonoverlappingOperator (const GFS& gfs_, const M& A,
+ const ParallelISTLHelper<GFS>& helper_)
+ DUNE_DEPRECATED
+ : gfs(gfs_), _A_(A)
+ {
+ }
+
+ //! Construct a non-overlapping operator
+ /**
+ * \param gfs_ GridFunctionsSpace for the vectors.
+ * \param A Matrix for this operator. This should be the locally
+ * assembled matrix.
+ *
+ * \note The constructed object stores references to all the objects
+ * given as parameters here. They should be valid for as long as
+ * the constructed object is used. They are not needed to
+ * destruct the constructed object.
+ */
+ NonoverlappingOperator (const GFS& gfs_, const M& A)
+ : gfs(gfs_), _A_(A)
+ { }
+
+ //! apply operator
+ /**
+ * Compute \f$y:=A(x)\f$ on this process, then make y consistent (sum up
+ * corresponding entries of y on the different processes and store the
+ * result back in y on each process).
+ */
+ virtual void apply (const X& x, Y& y) const
+ {
+ // apply local operator; now we have sum y_p = sequential y
+ _A_.mv(x,y);
+
+ // accumulate y on border
+ Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y);
+ if (gfs.gridView().comm().size()>1)
+ gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
+ }
+
+ //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
+ /**
+ * Compute \f$y:=\alpha A(x)\f$ on this process, then make y consistent
+ * (sum up corresponding entries of y on the different processes and
+ * store the result back in y on each process).
+ */
+ virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+ {
+ // apply local operator; now we have sum y_p = sequential y
+ _A_.usmv(alpha,x,y);
+
+ // accumulate y on border
+ Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y);
+ if (gfs.gridView().comm().size()>1)
+ gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
+ }
+
+ //! extract the matrix
+ virtual const M& getmat () const
+ {
+ return _A_;
+ }
+
+ private:
+ const GFS& gfs;
+ const M& _A_;
+ };
+
+ // parallel scalar product assuming no overlap
+ template<class GFS, class X>
+ class NonoverlappingScalarProduct : public Dune::ScalarProduct<X>
+ {
+ public:
+ //! export types
+ typedef X domain_type;
+ typedef typename X::ElementType field_type;
+
+ //! define the category
+ enum {category=Dune::SolverCategory::nonoverlapping};
+
+ /*! \brief Constructor needs to know the grid function space
+ */
+ NonoverlappingScalarProduct (const GFS& gfs_, const ParallelISTLHelper<GFS>& helper_)
+ : gfs(gfs_), helper(helper_)
+ {}
+
+ /*! \brief Dot product of two vectors.
+ It is assumed that the vectors are consistent on the interior+border
+ partition.
+ */
+ virtual field_type dot (const X& x, const X& y)
+ {
+ // do local scalar product on unique partition
+ field_type sum = 0;
+ for (typename X::size_type i=0; i<x.base().N(); ++i)
+ for (typename X::size_type j=0; j<x[i].N(); ++j)
+ sum += (x[i][j]*y[i][j])*helper.mask(i,j);
+
+ // do global communication
+ return gfs.gridView().comm().sum(sum);
+ }
+
+ /*! \brief Norm of a right-hand side vector.
+ The vector must be consistent on the interior+border partition
+ */
+ virtual double norm (const X& x)
+ {
+ return sqrt(static_cast<double>(this->dot(x,x)));
+ }
+
+ /*! \brief make additive vector consistent
+ */
+ void make_consistent (X& x) const
+ {
+ Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,x);
+ if (gfs.gridView().comm().size()>1)
+ gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
+ }
+
+ private:
+ const GFS& gfs;
+ const ParallelISTLHelper<GFS>& helper;
+ };
+
+ // parallel Richardson preconditioner
+ template<class GFS, class X, class Y>
+ class NonoverlappingRichardson : public Dune::Preconditioner<X,Y>
+ {
+ public:
+ //! \brief The domain type of the preconditioner.
+ typedef X domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::ElementType field_type;
+
+ // define the category
+ enum {
+ //! \brief The category the preconditioner is part of.
+ category=Dune::SolverCategory::nonoverlapping
+ };
+
+ //! \brief Constructor.
+ NonoverlappingRichardson (const GFS& gfs_, const ParallelISTLHelper<GFS>& helper_)
+ : gfs(gfs_), helper(helper_)
+ {
+ }
+
+ /*!
+ \brief Prepare the preconditioner.
+ */
+ virtual void pre (X& x, Y& b) {}
+
+ /*!
+ \brief Apply the precondioner.
+ */
+ virtual void apply (X& v, const Y& d)
+ {
+ v = d;
+ }
+
+ /*!
+ \brief Clean up.
+ */
+ virtual void post (X& x) {}
+
+ private:
+ const GFS& gfs;
+ const ParallelISTLHelper<GFS>& helper;
+ };
+
+ //! parallel non-overlapping Jacobi preconditioner
+ /**
+ * \tparam Diagonal Vector type used to store the diagonal of the matrix
+ * \tparam X Vector type used to store the result of applying the
+ * preconditioner.
+ * \tparam Y Vector type used to store the defect.
+ *
+ * The Jacobi preconditioner approximates the inverse of a matrix M by
+ * taking the diagonal diag(M) and inverting that. In the parallel case
+ * the matrix M is assumed to be inconsistent, so diagonal entries for
+ * dofs on the border are summed up over all relevant processes by this
+ * precoditioner before the inverse is computed.
+ */
+ template<class Diagonal, class X, class Y>
+ class NonoverlappingJacobi : public Dune::Preconditioner<X,Y>
+ {
+ typedef typename Diagonal::Backend DBackend;
+
+ std::size_t gfsSize;
+ Diagonal diagonal;
+
+ public:
+ //! The domain type of the operator.
+ /**
+ * The preconditioner is an inverse operator, so this is the output type
+ * of the preconditioner.
+ */
+ typedef X domain_type;
+ //! \brief The range type of the operator.
+ /**
+ * The preconditioner is an inverse operator, so this is the input type
+ * of the preconditioner.
+ */
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::ElementType field_type;
+
+ enum {
+ //! \brief The category the preconditioner is part of.
+ category=Dune::SolverCategory::nonoverlapping
+ };
+
+ //! \brief Constructor.
+ /**
+ * \param gfs The GridFunctionSpace the matrix and the vectors live on.
+ * \param m The matrix whose inverse the preconditioner should
+ * estimate. m is assumed to be inconsistent (i.e. rows for
+ * dofs on the border only contain the contribution of the
+ * local process).
+ *
+ * The preconditioner does not store any reference to the gfs or the
+ * matrix m. The diagonal of m is copied, since it has to be made
+ * consistent.
+ */
+ template<class GFS, class Matrix>
+ NonoverlappingJacobi(const GFS& gfs, const Matrix &m) :
+ gfsSize(gfs.size()), diagonal(gfs)
+ {
+ typedef typename Matrix::Backend MBackend;
+
+ for(std::size_t i = 0; i < gfsSize; ++i)
+ DBackend::access(diagonal, i) = MBackend::access(m, i, i);
+
+ AddDataHandle<GFS, Diagonal> addDH(gfs, diagonal);
+ gfs.gridView().communicate(addDH,
+ InteriorBorder_InteriorBorder_Interface,
+ ForwardCommunication);
+ }
+
+ //! Prepare the preconditioner.
+ virtual void pre (X& x, Y& b) {}
+
+ //! Apply the precondioner.
+ /*
+ * For this preconditioner, this method works with both consistent and
+ * inconsistent vectors: if d is consistent, v will be consistent, if d
+ * is inconsistent, v will be inconsistent.
+ */
+ virtual void apply (X& v, const Y& d)
+ {
+ typedef typename X::Backend XBackend;
+ typedef typename Y::Backend YBackend;
+
+ for(std::size_t i = 0; i < gfsSize; ++i)
+ XBackend::access(v, i) =
+ YBackend::access(d, i) / DBackend::access(diagonal, i);
+ }
+
+ //! Clean up.
+ virtual void post (X& x) {}
+ };
+
+ //! \addtogroup PDELab_novlpsolvers Nonoverlapping Solvers
+ //! \{
+
+ //! \brief Nonoverlapping parallel CG solver without preconditioner
+ template<class GFS>
+ class ISTLBackend_NOVLP_CG_NOPREC
+ {
+ typedef Dune::PDELab::ParallelISTLHelper<GFS> PHELPER;
+
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_NOVLP_CG_NOPREC (const GFS& gfs_,
+ unsigned maxiter_=5000,
+ int verbose_=1)
+ : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
+ {}
+
+ /*! \brief compute global norm of a vector
+
+ \param[in] v the given vector
+ */
+ template<class V>
+ typename V::ElementType norm (const V& v) const
+ {
+ V x(v); // make a copy because it has to be made consistent
+ typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+ psp.make_consistent(x);
+ return psp.norm(x);
+ }
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename V::ElementType reduction)
+ {
+ typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP;
+ POP pop(gfs,A);
+ typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+ typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH;
+ PRICH prich(gfs,phelper);
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ Dune::CGSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z,r,stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ /*! \brief Return access to result data */
+ const Dune::PDELab::LinearSolverResult<double>& result() const
+ {
+ return res;
+ }
+
+ private:
+ const GFS& gfs;
+ PHELPER phelper;
+ Dune::PDELab::LinearSolverResult<double> res;
+ unsigned maxiter;
+ int verbose;
+ };
+
+ //! \brief Nonoverlapping parallel CG solver with Jacobi preconditioner
+ template<class GFS>
+ class ISTLBackend_NOVLP_CG_Jacobi
+ {
+ typedef Dune::PDELab::ParallelISTLHelper<GFS> PHELPER;
+
+ const GFS& gfs;
+ PHELPER phelper;
+ LinearSolverResult<double> res;
+ unsigned maxiter;
+ int verbose;
+
+ public:
+ //! make a linear solver object
+ /**
+ * \param gfs_ A grid function space
+ * \param maxiter_ Maximum number of iterations to do.
+ * \param verbose_ Verbosity level, directly handed to the CGSolver.
+ */
+ explicit ISTLBackend_NOVLP_CG_Jacobi(const GFS& gfs_,
+ unsigned maxiter_ = 5000,
+ int verbose_ = 1) :
+ gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
+ {}
+
+ //! compute global norm of a vector
+ /**
+ * \param v The vector to compute the norm of. Should be an
+ * inconsistent vector (i.e. the entries corresponding a DoF on
+ * the border should only contain the summand of this process).
+ */
+ template<class V>
+ typename V::ElementType norm (const V& v) const
+ {
+ V x(v); // make a copy because it has to be made consistent
+ typedef NonoverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+ psp.make_consistent(x);
+ return psp.norm(x);
+ }
+
+ //! solve the given linear system
+ /**
+ * \param A The matrix to solve. Should be a matrix from one of
+ * PDELabs ISTL backends (only ISTLBCRSMatrixBackend at
+ * the moment).
+ * \param z The solution vector to be computed
+ * \param r Right hand side
+ * \param reduction to be achieved
+ *
+ * Solve the linear system A*z=r such that
+ * norm(A*z0-r)/norm(A*z-r) < reduction where z0 is the initial value of
+ * z.
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename V::ElementType reduction)
+ {
+ typedef NonoverlappingOperator<GFS,M,V,W> POP;
+ POP pop(gfs,A);
+ typedef NonoverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+
+ typedef typename M::ElementType MField;
+ typedef typename BackendVectorSelector<GFS,MField>::Type Diagonal;
+ typedef NonoverlappingJacobi<Diagonal,V,W> PPre;
+ PPre ppre(gfs,A);
+
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ CGSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb);
+ InverseOperatorResult stat;
+ solver.apply(z,r,stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ //! Return access to result data
+ const LinearSolverResult<double>& result() const
+ { return res; }
+ };
+
+ //! \brief Nonoverlapping parallel BiCGStab solver without preconditioner
+ template<class GFS>
+ class ISTLBackend_NOVLP_BCGS_NOPREC
+ {
+ typedef Dune::PDELab::ParallelISTLHelper<GFS> PHELPER;
+
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_NOVLP_BCGS_NOPREC (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1)
+ : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_)
+ {}
+
+ /*! \brief compute global norm of a vector
+
+ \param[in] v the given vector
+ */
+ template<class V>
+ typename V::ElementType norm (const V& v) const
+ {
+ V x(v); // make a copy because it has to be made consistent
+ typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+ psp.make_consistent(x);
+ return psp.norm(x);
+ }
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename V::ElementType reduction)
+ {
+ typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP;
+ POP pop(gfs,A);
+ typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+ typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH;
+ PRICH prich(gfs,phelper);
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ Dune::BiCGSTABSolver<V> solver(pop,psp,prich,reduction,maxiter,verb);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z,r,stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ /*! \brief Return access to result data */
+ const Dune::PDELab::LinearSolverResult<double>& result() const
+ {
+ return res;
+ }
+
+ private:
+ const GFS& gfs;
+ PHELPER phelper;
+ Dune::PDELab::LinearSolverResult<double> res;
+ unsigned maxiter;
+ int verbose;
+ };
+
+ //! Solver to be used for explicit time-steppers with (block-)diagonal mass matrix
+ template<typename GFS>
+ class ISTLBackend_NOVLP_ExplicitDiagonal
+ {
+ typedef Dune::PDELab::ParallelISTLHelper<GFS> PHELPER;
+
+ const GFS& gfs;
+ PHELPER phelper;
+ Dune::PDELab::LinearSolverResult<double> res;
+
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ GridFunctionSpace, used to identify DoFs for parallel
+ communication
+ */
+ explicit ISTLBackend_NOVLP_ExplicitDiagonal(const GFS& gfs_)
+ : gfs(gfs_), phelper(gfs)
+ {}
+
+ /*! \brief compute global norm of a vector
+
+ \param[in] v the given vector
+ */
+ template<class V>
+ typename V::ElementType norm (const V& v) const
+ {
+ typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
+ V x(v); // make a copy because it has to be made consistent
+ PSP psp(gfs,phelper);
+ psp.make_consistent(x);
+ return psp.norm(x);
+ }
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename W::ElementType reduction)
+ {
+ Dune::SeqJac<M,V,W> jac(A,1,1.0);
+ jac.pre(z,r);
+ jac.apply(z,r);
+ jac.post(z);
+ if (gfs.gridView().comm().size()>1)
+ {
+ Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,z);
+ gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
+ }
+ res.converged = true;
+ res.iterations = 1;
+ res.elapsed = 0.0;
+ res.reduction = reduction;
+ res.conv_rate = reduction; // pow(reduction,1.0/1)
+ }
+
+ /*! \brief Return access to result data */
+ const Dune::PDELab::LinearSolverResult<double>& result() const
+ {
+ return res;
+ }
+ };
+ //! \} Nonoverlapping Solvers
+
+ /**
+ * @brief Helper class for adding up matrix entries on border.
+ * @tparam GridOperator The grid operator to work on.
+ * @tparam MatrixType The MatrixType.
+ */
+ template<class GridOperator, class MatrixType>
+ class VertexExchanger
+ {
+ typedef MatrixType Matrix;
+ typedef typename GridOperator::Traits GridOperatorTraits;
+ typedef typename GridOperatorTraits::JacobianField Scalar;
+ typedef typename GridOperatorTraits::TrialGridFunctionSpace GFS;
+ typedef typename GFS::Traits::GridViewType GridView;
+ enum {dim = GridView::dimension};
+ typedef typename GridView::Traits::Grid Grid;
+ typedef typename Matrix::block_type BlockType;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ typedef typename Grid::Traits::GlobalIdSet IDS;
+ typedef typename IDS::IdType IdType;
+ typedef typename Matrix::RowIterator RowIterator;
+ typedef typename Matrix::ColIterator ColIterator;
+
+ public:
+ /*! \brief Constructor. Sets up the local to global relations.
+
+ \param[in] gridView The grid view to operate on.
+ */
+ VertexExchanger(const GridView& gridView)
+ : gridView_(gridView)
+ {
+ gid2Index_.clear();
+ index2GID_.clear();
+
+
+ VertexIterator vertexEndIt = gridView_.template end<dim>();
+ for (VertexIterator vertexIt = gridView_.template begin<dim>(); vertexIt != vertexEndIt; ++vertexIt)
+ {
+ if (vertexIt->partitionType() == BorderEntity)
+ {
+ int localIdx = gridView_.indexSet().index(*vertexIt);
+ IdType globalIdx = gridView_.grid().globalIdSet().id(*vertexIt);
+
+ std::pair<IdType,int> g2iPair(globalIdx, localIdx);
+ gid2Index_.insert(g2iPair);
+
+ std::pair<int,IdType> i2gPair(localIdx, globalIdx);
+ index2GID_.insert(i2gPair);
+
+ }
+ }
+ }
+
+ //! A DataHandle class to exchange matrix sparsity patterns
+ class MatPatternExchange
+ : public CommDataHandleIF<MatPatternExchange,IdType> {
+ typedef typename Matrix::RowIterator RowIterator;
+ typedef typename Matrix::ColIterator ColIterator;
+ public:
+ //! Export type of data for message buffer
+ typedef IdType DataType;
+
+ /** @brief Returns true if data for given valid codim should be communicated
+ */
+ bool contains (int dim, int codim) const
+ {
+ return (codim==dim);
+ }
+
+ /** @brief Returns true if size of data per entity of given dim and codim is a constant
+ */
+ bool fixedsize (int dim, int codim) const
+ {
+ return false;
+ }
+
+ /** @brief How many objects of type DataType have to be sent for a given entity
+ */
+ template<class EntityType>
+ size_t size (EntityType& e) const
+ {
+ int i = gridView_.indexSet().index(e);
+ int n = 0;
+ for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+ {
+ typename std::map<int,IdType>::const_iterator it = index2GID_.find(j.index());
+ if (it != index2GID_.end())
+ n++;
+ }
+
+ return n;
+ }
+
+ /** @brief Pack data from user to message buffer
+ */
+ template<class MessageBuffer, class EntityType>
+ void gather (MessageBuffer& buff, const EntityType& e) const
+ {
+ int i = gridView_.indexSet().index(e);
+ for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+ {
+ typename std::map<int,IdType>::const_iterator it=index2GID_.find(j.index());
+ if (it != index2GID_.end())
+ buff.write(it->second);
+ }
+
+ }
+
+ /** @brief Unpack data from message buffer to user
+ */
+ template<class MessageBuffer, class EntityType>
+ void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
+ {
+ int i = gridView_.indexSet().index(e);
+ for (size_t k = 0; k < n; k++)
+ {
+ IdType id;
+ buff.read(id);
+ // only add entries corresponding to border entities
+ typename std::map<IdType,int>::const_iterator it = gid2Index_.find(id);
+ if (it != gid2Index_.end()
+ && sparsity_[i].find(it->second) == sparsity_[i].end()
+ && helper_.ghost(it->second,0) != 1<<24)
+ sparsity_[i].insert(it->second);
+ }
+ }
+
+ /**
+ * @brief Get the communicated sparsity pattern
+ * @return the vector with the sparsity pattern
+ */
+ std::vector<std::set<int> >& sparsity ()
+ {
+ return sparsity_;
+ }
+
+ /** @brief Constructor
+ @param[in] gridView Grid view.
+ @param[in] g2i Global to local index map.
+ @param[in] i2g Local to global index map.
+ @param[in] A Matrix to operate on.
+ @param[in] helper parallel istl helper.
+ */
+ MatPatternExchange (const GridView& gridView,
+ const std::map<IdType,int>& g2i,
+ const std::map<int,IdType>& i2g, Matrix& A,
+ const ParallelISTLHelper<GFS>& helper)
+ : gridView_(gridView), gid2Index_(g2i), index2GID_(i2g),
+ sparsity_(A.N()), A_(A), helper_(helper)
+ {}
+
+ private:
+ const GridView& gridView_;
+ const std::map<IdType,int>& gid2Index_;
+ const std::map<int,IdType>& index2GID_;
+ std::vector<std::set<int> > sparsity_;
+ Matrix& A_;
+ const ParallelISTLHelper<GFS>& helper_;
+ };
+
+ //! Local matrix blocks associated with the global id set
+ struct MatEntry
+ {
+ IdType first;
+ BlockType second;
+ MatEntry (const IdType& f, const BlockType& s) : first(f),second(s) {}
+ MatEntry () {}
+ };
+
+ //! A DataHandle class to exchange matrix entries
+ class MatEntryExchange
+ : public CommDataHandleIF<MatEntryExchange,MatEntry> {
+ typedef typename Matrix::RowIterator RowIterator;
+ typedef typename Matrix::ColIterator ColIterator;
+ public:
+ //! Export type of data for message buffer
+ typedef MatEntry DataType;
+
+ /** @brief Returns true if data for given valid codim should be communicated
+ */
+ bool contains (int dim, int codim) const
+ {
+ return (codim==dim);
+ }
+
+ /** @brief Returns true if size of data per entity of given dim and codim is a constant
+ */
+ bool fixedsize (int dim, int codim) const
+ {
+ return false;
+ }
+
+ /** @brief How many objects of type DataType have to be sent for a given entity
+ */
+ template<class EntityType>
+ size_t size (EntityType& e) const
+ {
+ int i = gridView_.indexSet().index(e);
+ int n = 0;
+ for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+ {
+ typename std::map<int,IdType>::const_iterator it = index2GID_.find(j.index());
+ if (it != index2GID_.end())
+ n++;
+ }
+
+ return n;
+ }
+
+ /** @brief Pack data from user to message buffer
+ */
+ template<class MessageBuffer, class EntityType>
+ void gather (MessageBuffer& buff, const EntityType& e) const
+ {
+ int i = gridView_.indexSet().index(e);
+ for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
+ {
+ typename std::map<int,IdType>::const_iterator it=index2GID_.find(j.index());
+ if (it != index2GID_.end())
+ buff.write(MatEntry(it->second,*j));
+ }
+
+ }
+
+ /** @brief Unpack data from message buffer to user
+ */
+ template<class MessageBuffer, class EntityType>
+ void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
+ {
+ int i = gridView_.indexSet().index(e);
+ for (size_t k = 0; k < n; k++)
+ {
+ MatEntry m;
+ buff.read(m);
+ // only add entries corresponding to border entities
+ typename std::map<IdType,int>::const_iterator it = gid2Index_.find(m.first);
+ if (it != gid2Index_.end())
+ if (A_[i].find(it->second) != A_[i].end())
+ A_[i][it->second] += m.second;
+ }
+ }
+
+ /** @brief Constructor
+ @param[in] gridView Grid view.
+ @param[in] g2i Global to local index map.
+ @param[in] i2g Local to global index map.
+ @param[in] A Matrix to operate on.
+ */
+ MatEntryExchange (const GridView& gridView, const std::map<IdType,int>& g2i,
+ const std::map<int,IdType>& i2g,
+ Matrix& A)
+ : gridView_(gridView), gid2Index_(g2i), index2GID_(i2g), A_(A)
+ {}
+
+ private:
+ const GridView& gridView_;
+ const std::map<IdType,int>& gid2Index_;
+ const std::map<int,IdType>& index2GID_;
+ Matrix& A_;
+ };
+
+ /** @brief communicates values for the sparsity pattern of the new matrix.
+ @param A Matrix to operate on.
+ @param helper ParallelelISTLHelper.
+ */
+ void getextendedmatrix (Matrix& A,const ParallelISTLHelper<GFS>& helper)
+ {
+ if (gridView_.comm().size() > 1) {
+ Matrix tmp(A);
+ std::size_t nnz=0;
+ // get entries from other processes
+ MatPatternExchange datahandle(gridView_, gid2Index_, index2GID_, A, helper);
+ gridView_.communicate(datahandle, InteriorBorder_InteriorBorder_Interface, ForwardCommunication);
+ std::vector<std::set<int> >& sparsity = datahandle.sparsity();
+ // add own entries, count number of nonzeros
+ for (RowIterator i = A.begin(); i != A.end(); ++i){
+ for (ColIterator j = A[i.index()].begin(); j != A[i.index()].end(); ++j){
+ if (sparsity[i.index()].find(j.index()) == sparsity[i.index()].end())
+ sparsity[i.index()].insert(j.index());
+ }
+ nnz += sparsity[i.index()].size();
+ }
+ A.setSize(tmp.N(), tmp.N(), nnz);
+ A.setBuildMode(Matrix::row_wise);
+ typename Matrix::CreateIterator citer = A.createbegin();
+ typedef typename std::vector<std::set<int> >::const_iterator Iter;
+ for (Iter i = sparsity.begin(), end = sparsity.end(); i!=end; ++i, ++citer){
+ typedef typename std::set<int>::const_iterator SIter;
+ for (SIter si = i->begin(), send = i->end(); si!=send; ++si)
+ citer.insert(*si);
+ }
+ // set matrix old values
+ A = 0;
+ for (RowIterator i = tmp.begin(); i != tmp.end(); ++i)
+ for (ColIterator j = tmp[i.index()].begin(); j != tmp[i.index()].end(); ++j){
+ A[i.index()][j.index()] = tmp[i.index()][j.index()];
+ }
+ }
+ }
+
+ /** @brief Sums up the entries corresponding to border vertices.
+ @param A Matrix to operate on.
+ */
+ void sumEntries (Matrix& A)
+ {
+ if (gridView_.comm().size() > 1)
+ {
+ MatEntryExchange datahandle(gridView_, gid2Index_, index2GID_, A);
+ gridView_.communicate(datahandle, InteriorBorder_InteriorBorder_Interface, ForwardCommunication);
+ }
+ }
+
+ private:
+ const GridView& gridView_;
+ std::map<IdType,int> gid2Index_;
+ std::map<int,IdType> index2GID_;
+ };
+
+ template<class GO,
+ template<class,class,class,int> class Preconditioner,
+ template<class> class Solver>
+ class ISTLBackend_NOVLP_BASE_PREC
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+ typedef Dune::PDELab::ParallelISTLHelper<GFS> PHELPER;
+
+ public:
+ /*! \brief Constructor.
+
+ \param[in] gfs_ a grid function space
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] steps_ number of preconditioner steps to apply as inner iteration
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_NOVLP_BASE_PREC (const GFS& gfs_, unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1)
+ : gfs(gfs_), phelper(gfs_,verbose_),
+ maxiter(maxiter_), steps(steps_), verbose(verbose_)
+ {}
+
+ /*! \brief Compute global norm of a vector.
+
+ \param[in] v the given vector
+ */
+ template<class Vector>
+ typename Vector::ElementType norm (const Vector& v) const
+ {
+ Vector x(v); // make a copy because it has to be made consistent
+ typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,Vector> PSP;
+ PSP psp(gfs,phelper);
+ psp.make_consistent(x);
+ return psp.norm(x);
+ }
+
+ /*! \brief Solve the given linear system.
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename V::ElementType reduction)
+ {
+ typedef typename CommSelector<96,Dune::MPIHelper::isFake>::type Comm;
+ typedef typename M::BaseT MatrixType;
+ MatrixType& mat=A.base();
+ typedef typename BlockProcessor<GFS>::template AMGVectorTypeSelector<V>::Type VectorType;
+#if HAVE_MPI
+ Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
+ typedef VertexExchanger<GO,MatrixType> Exchanger;
+ Exchanger exchanger(gfs.gridView());
+ exchanger.getextendedmatrix(mat,phelper);
+ exchanger.sumEntries(mat);
+ phelper.createIndexSetAndProjectForAMG(mat, oocc);
+ typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
+ Smoother smoother(mat, steps, 1.0);
+ typedef Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> PSP;
+ PSP psp(oocc);
+ typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
+ Operator oop(mat,oocc);
+ typedef Dune::NonoverlappingBlockPreconditioner<Comm, Smoother> ParSmoother;
+ ParSmoother parsmoother(smoother, oocc);
+#else
+ Comm oocc(gfs.gridView().comm());
+ typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
+ ParSmoother parsmoother(mat, steps, 1.0);
+ typedef Dune::SeqScalarProduct<VectorType> PSP;
+ PSP psp;
+ typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
+ Operator oop(mat);
+#endif
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ Solver<VectorType> solver(oop,psp,parsmoother,reduction,maxiter,verb);
+ Dune::InverseOperatorResult stat;
+ //make r consistent
+ if (gfs.gridView().comm().size()>1){
+ Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r);
+ gfs.gridView().communicate(adddh,
+ Dune::InteriorBorder_InteriorBorder_Interface,
+ Dune::ForwardCommunication);
+ }
+
+ solver.apply(z,r,stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ /*! \brief Return access to result data. */
+ const Dune::PDELab::LinearSolverResult<double>& result() const
+ {
+ return res;
+ }
+
+ private:
+ const GFS& gfs;
+ PHELPER phelper;
+ Dune::PDELab::LinearSolverResult<double> res;
+ unsigned maxiter;
+ unsigned steps;
+ int verbose;
+ };
+
+ //! \addtogroup PDELab_novlpsolvers Nonoverlapping Solvers
+ //! \{
+
+ /**
+ * @brief Nonoverlapping parallel BiCGSTAB solver preconditioned by block SSOR.
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ *
+ * The solver uses a NonoverlappingBlockPreconditioner with underlying
+ * sequential SSOR preconditioner. The crucial step is to add up the matrix entries
+ * corresponding to the border vertices on each process. This is achieved by
+ * performing a VertexExchanger::sumEntries(Matrix&) before constructing the
+ * sequential SSOR.
+ */
+ template<class GO>
+ class ISTLBackend_NOVLP_BCGS_SSORk
+ : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] steps_ number of SSOR steps to apply as inner iteration
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_NOVLP_BCGS_SSORk (const GFS& gfs_, unsigned maxiter_=5000,
+ int steps_=5, int verbose_=1)
+ : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>(gfs_, maxiter_, steps_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Nonoverlapping parallel CG solver preconditioned by block SSOR.
+ */
+ template<class GO>
+ class ISTLBackend_NOVLP_CG_SSORk
+ : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] steps_ number of SSOR steps to apply as inner iteration
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_NOVLP_CG_SSORk (const GFS& gfs_, unsigned maxiter_=5000,
+ int steps_=5, int verbose_=1)
+ : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>(gfs_, maxiter_, steps_, verbose_)
+ {}
+ };
+ //! \} Nonoverlapping Solvers
+ //! \} group Backend
+
+ template<class GO,int s, template<class,class,class,int> class Preconditioner,
+ template<class> class Solver>
+ class ISTLBackend_AMG_NOVLP : public LinearResultStorage
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+ typedef typename Dune::PDELab::ParallelISTLHelper<GFS> PHELPER;
+ typedef typename GO::Traits::Jacobian M;
+ typedef typename M::BaseT MatrixType;
+ typedef typename GO::Traits::Domain V;
+ typedef typename BlockProcessor<GFS>::template AMGVectorTypeSelector<V>::Type VectorType;
+ typedef typename CommSelector<s,Dune::MPIHelper::isFake>::type Comm;
+#if HAVE_MPI
+ typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
+ typedef Dune::NonoverlappingBlockPreconditioner<Comm,Smoother> ParSmoother;
+ typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
+#else
+ typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
+ typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
+#endif
+ typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
+ typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
+ typedef Dune::Amg::Parameters Parameters;
+
+ public:
+ ISTLBackend_AMG_NOVLP(const GFS& gfs_, unsigned maxiter_=5000,
+ int verbose_=1, bool reuse_=false,
+ bool usesuperlu_=true)
+ : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_),
+ params(15,2000,1.2,1.6,Dune::Amg::noAccu/*Dune::Amg::atOnceAccu*/),
+ verbose(verbose_), reuse(reuse_), firstapply(true),
+ usesuperlu(usesuperlu_)
+ {
+ params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
+ params.setDebugLevel(verbose_);
+#if !HAVE_SUPERLU
+ if (gfs.gridView().comm().rank() == 0 && usesuperlu == true)
+ {
+ std::cout << "WARNING: You are using AMG without SuperLU!"
+ << " Please consider installing SuperLU,"
+ << " or set the usesuperlu flag to false"
+ << " to suppress this warning." << std::endl;
+ }
+#endif
+ }
+
+ /*! \brief set AMG parameters
+
+ \param[in] params_ a parameter object of Type Dune::Amg::Parameters
+ */
+ void setParameters(const Parameters& params_)
+ {
+ params = params_;
+ }
+
+ void setparams(Parameters params_)
+ {
+ params = params_;
+ }
+
+ /**
+ * @brief Get the parameters describing the behaviuour of AMG.
+ *
+ * The returned object can be adjusted to ones needs and then can be
+ * reset using setParameters.
+ * @return The object holding the parameters of AMG.
+ */
+ const Parameters& parameters() const
+ {
+ return params;
+ }
+
+ /*! \brief compute global norm of a vector
+
+ \param[in] v the given vector
+ */
+ typename V::ElementType norm (const V& v) const
+ {
+ V x(v); // make a copy because it has to be made consistent
+ typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+ psp.make_consistent(x);
+ return psp.norm(x);
+ }
+
+ void apply(MatrixType& A, V& z, V& r, typename V::ElementType reduction)
+ {
+ Timer watch;
+ MatrixType& mat=A;
+ typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
+ Dune::Amg::FirstDiagonal> > Criterion;
+#if HAVE_MPI
+ Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping);
+ typedef VertexExchanger<GO,MatrixType> Exchanger;
+ Exchanger exchanger(gfs.gridView());
+ exchanger.getextendedmatrix(mat,phelper);
+ exchanger.sumEntries(mat);
+ phelper.createIndexSetAndProjectForAMG(mat, oocc);
+ Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
+ Operator oop(mat, oocc);
+#else
+ Comm oocc(gfs.gridView().comm());
+ Operator oop(mat);
+ Dune::SeqScalarProduct<VectorType> sp;
+#endif
+ SmootherArgs smootherArgs;
+ smootherArgs.iterations = 1;
+ smootherArgs.relaxationFactor = 1;
+ //use noAccu or atOnceAccu
+ Criterion criterion(params);
+ stats.tprepare=watch.elapsed();
+ watch.reset();
+
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ //only construct a new AMG if the matrix changes
+ if (reuse==false || firstapply==true){
+ amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
+ firstapply = false;
+ stats.tsetup = watch.elapsed();
+ stats.levels = amg->maxlevels();
+ stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
+ }
+
+ Dune::InverseOperatorResult stat;
+ // make r consistent
+ if (gfs.gridView().comm().size()>1) {
+ Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r);
+ gfs.gridView().communicate(adddh,
+ Dune::InteriorBorder_InteriorBorder_Interface,
+ Dune::ForwardCommunication);
+ }
+ watch.reset();
+ Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb);
+ solver.apply(BlockProcessor<GFS>::getVector(z),BlockProcessor<GFS>::getVector(r),stat);
+ stats.tsolve= watch.elapsed();
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ /**
+ * @brief Get statistics of the AMG solver (no of levels, timings).
+ * @return statistis of the AMG solver.
+ */
+ const ISTLAMGStatistics& statistics() const
+ {
+ return stats;
+ }
+
+ private:
+ const GFS& gfs;
+ PHELPER phelper;
+ unsigned maxiter;
+ Parameters params;
+ int verbose;
+ bool reuse;
+ bool firstapply;
+ bool usesuperlu;
+ Dune::shared_ptr<AMG> amg;
+ ISTLAMGStatistics stats;
+ };
+
+ template<class GO, int s=96>
+ class ISTLBackend_NOVLP_CG_AMG_SSOR
+ : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+
+ public:
+ ISTLBackend_NOVLP_CG_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
+ int verbose_=1, bool reuse_=false,
+ bool usesuperlu_=true)
+ : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>
+ (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+
+ template<class GO, int s=96>
+ class ISTLBackend_NOVLP_BCGS_AMG_SSOR
+ : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+
+ public:
+ ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
+ int verbose_=1, bool reuse_=false,
+ bool usesuperlu_=true)
+ : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ template<class GO, int s=96>
+ class ISTLBackend_NOVLP_LS_AMG_SSOR
+ : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+
+ public:
+ ISTLBackend_NOVLP_LS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
+ int verbose_=1, bool reuse_=false,
+ bool usesuperlu_=true)
+ : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>
+ (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/linear/ovlpistlsolverbackend.hh b/dumux/linear/ovlpistlsolverbackend.hh
new file mode 100644
index 0000000000..6fd40b3b48
--- /dev/null
+++ b/dumux/linear/ovlpistlsolverbackend.hh
@@ -0,0 +1,1000 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=8 sw=2 sts=2:
+/*****************************************************************************
+ * This file is based on the file of DUNE-PDELab with the same name, *
+ * http://www.dune-project.org/pdelab/. *
+ * Some changes are made such that the backends work for the vectors and *
+ * matrices used in Dumux. *
+ *****************************************************************************/
+#ifndef DUNE_OVLPISTLSOLVERBACKEND_HH
+#define DUNE_OVLPISTLSOLVERBACKEND_HH
+
+#include <dune/common/deprecated.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/solvercategory.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/pinfo.hh>
+#include <dune/istl/io.hh>
+#include <dune/istl/superlu.hh>
+
+#include <dune/pdelab/backend/parallelistlhelper.hh>
+
+#include "istlvectorbackend.hh"
+#include "seqistlsolverbackend.hh"
+
+namespace Dumux {
+
+ //! \addtogroup Backend
+ //! \ingroup PDELab
+ //! \{
+
+ //========================================================
+ // Generic support for overlapping grids
+ // (need to be used with appropriate constraints)
+ //========================================================
+
+ // operator that resets result to zero at constrained DOFS
+ template<class CC, class M, class X, class Y>
+ class OverlappingOperator : public Dune::AssembledLinearOperator<M,X,Y>
+ {
+ public:
+ //! export types
+ typedef M matrix_type;
+ typedef X domain_type;
+ typedef Y range_type;
+ typedef typename X::field_type field_type;
+
+ //redefine the category, that is the only difference
+ enum {category=Dune::SolverCategory::overlapping};
+
+ OverlappingOperator (const CC& cc_, const M& A)
+ : cc(cc_), _A_(A)
+ {}
+
+ //! apply operator to x: \f$ y = A(x) \f$
+ virtual void apply (const X& x, Y& y) const
+ {
+ _A_.mv(x,y);
+ Dune::PDELab::set_constrained_dofs(cc,0.0,y);
+ }
+
+ //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
+ virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+ {
+ _A_.usmv(alpha,x,y);
+ Dune::PDELab::set_constrained_dofs(cc,0.0,y);
+ }
+
+ //! get matrix via *
+ virtual const M& getmat () const
+ {
+ return _A_;
+ }
+
+ private:
+ const CC& cc;
+ const M& _A_;
+ };
+
+ // new scalar product assuming at least overlap 1
+ // uses unique partitioning of nodes for parallelization
+ template<class GFS, class X>
+ class OverlappingScalarProduct : public Dune::ScalarProduct<X>
+ {
+ public:
+ //! export types
+ typedef X domain_type;
+ typedef typename X::ElementType field_type;
+
+ //! define the category
+ enum {category=Dune::SolverCategory::overlapping};
+
+ /*! \brief Constructor needs to know the grid function space
+ */
+ OverlappingScalarProduct (const GFS& gfs_, const ParallelISTLHelper<GFS>& helper_)
+ : gfs(gfs_), helper(helper_)
+ {}
+
+
+ /*! \brief Dot product of two vectors.
+ It is assumed that the vectors are consistent on the interior+border
+ partition.
+ */
+ virtual field_type dot (const X& x, const X& y)
+ {
+ // do local scalar product on unique partition
+ field_type sum = 0;
+ for (typename X::size_type i=0; i<x.base().N(); ++i)
+ for (typename X::size_type j=0; j<x.base()[i].N(); ++j)
+ sum += (x.base()[i][j]*y.base()[i][j])*helper.mask(i,j);
+
+ // do global communication
+ return gfs.gridView().comm().sum(sum);
+ }
+
+ /*! \brief Norm of a right-hand side vector.
+ The vector must be consistent on the interior+border partition
+ */
+ virtual double norm (const X& x)
+ {
+ return sqrt(static_cast<double>(this->dot(x,x)));
+ }
+
+ private:
+ const GFS& gfs;
+ const ParallelISTLHelper<GFS>& helper;
+ };
+
+ // wrapped sequential preconditioner
+ template<class CC, class GFS, class P>
+ class OverlappingWrappedPreconditioner
+ : public Dune::Preconditioner<typename Dune::PDELab::BackendVectorSelector<GFS,typename P::domain_type::field_type>::Type,
+ typename Dune::PDELab::BackendVectorSelector<GFS,typename P::range_type::field_type>::Type>
+ {
+ public:
+ //! \brief The domain type of the preconditioner.
+ typedef typename Dune::PDELab::BackendVectorSelector<GFS,typename P::domain_type::field_type>::Type
+ domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef typename Dune::PDELab::BackendVectorSelector<GFS,typename P::range_type::field_type>::Type
+ range_type;
+
+ // define the category
+ enum {
+ //! \brief The category the preconditioner is part of.
+ category=Dune::SolverCategory::overlapping
+ };
+
+ //! Constructor.
+ OverlappingWrappedPreconditioner (const GFS& gfs_, P& prec_, const CC& cc_,
+ const ParallelISTLHelper<GFS>& helper_)
+ : gfs(gfs_), prec(prec_), cc(cc_), helper(helper_)
+ {}
+
+ /*!
+ \brief Prepare the preconditioner.
+ */
+ virtual void pre (domain_type& x, range_type& b)
+ {
+ prec.pre(x,b);
+ }
+
+ /*!
+ \brief Apply the precondioner.
+ */
+ virtual void apply (domain_type& v, const range_type& d)
+ {
+ range_type dd(d);
+ set_constrained_dofs(cc,0.0,dd);
+ prec.apply(v,dd);
+ Dune::PDELab::AddDataHandle<GFS,domain_type,typename domain_type::field_type> adddh(gfs,v);
+ if (gfs.gridView().comm().size()>1)
+ gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
+ }
+
+ /*!
+ \brief Clean up.
+ */
+ virtual void post (domain_type& x)
+ {
+ prec.post(x);
+ }
+
+ private:
+ const GFS& gfs;
+ P& prec;
+ const CC& cc;
+ const ParallelISTLHelper<GFS>& helper;
+ };
+
+
+#if HAVE_SUPERLU
+ // exact subdomain solves with SuperLU as preconditioner
+ template<class GFS, class M, class X, class Y>
+ class SuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
+ {
+ typedef typename M::BaseT ISTLM;
+
+ public:
+ //! \brief The domain type of the preconditioner.
+ typedef X domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::ElementType field_type;
+
+
+ // define the category
+ enum {
+ //! \brief The category the preconditioner is part of.
+ category=Dune::SolverCategory::overlapping
+ };
+
+ /*! \brief Constructor.
+
+ Constructor gets all parameters to operate the prec.
+ \param gfs_ The grid function space.
+ \param A_ The matrix to operate on.
+ */
+ SuperLUSubdomainSolver (const GFS& gfs_, const M& A_)
+ : gfs(gfs_), A(A_), solver(A_,false) // this does the decomposition
+ {}
+
+ /*!
+ \brief Prepare the preconditioner.
+ */
+ virtual void pre (X& x, Y& b) {}
+
+ /*!
+ \brief Apply the precondioner.
+ */
+ virtual void apply (X& v, const Y& d)
+ {
+ Dune::InverseOperatorResult stat;
+ Y b(d); // need copy, since solver overwrites right hand side
+ solver.apply(v,b,stat);
+ Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,v);
+ if (gfs.gridview().comm().size()>1)
+ gfs.gridview().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
+ }
+
+ /*!
+ \brief Clean up.
+ */
+ virtual void post (X& x) {}
+
+ private:
+ const GFS& gfs;
+ const M& A;
+ Dune::SuperLU<ISTLM> solver;
+ };
+
+ // exact subdomain solves with SuperLU as preconditioner
+ template<class GFS, class M, class X, class Y>
+ class RestrictedSuperLUSubdomainSolver : public Dune::Preconditioner<X,Y>
+ {
+ typedef typename M::BaseT ISTLM;
+
+ public:
+ //! \brief The domain type of the preconditioner.
+ typedef X domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::ElementType field_type;
+
+
+ // define the category
+ enum {
+ //! \brief The category the preconditioner is part of.
+ category=Dune::SolverCategory::overlapping
+ };
+
+ /*! \brief Constructor.
+
+ Constructor gets all parameters to operate the prec.
+ \param gfs_ The grid function space.
+ \param A_ The matrix to operate on.
+ \param helper_ The parallel istl helper.
+ */
+ RestrictedSuperLUSubdomainSolver (const GFS& gfs_, const M& A_,
+ const ParallelISTLHelper<GFS>& helper_)
+ : gfs(gfs_), A(A_), solver(A_,false), helper(helper_) // this does the decomposition
+ {}
+
+ /*!
+ \brief Prepare the preconditioner.
+ */
+ virtual void pre (X& x, Y& b) {}
+
+ /*!
+ \brief Apply the precondioner.
+ */
+ virtual void apply (X& v, const Y& d)
+ {
+ Dune::InverseOperatorResult stat;
+ Y b(d); // need copy, since solver overwrites right hand side
+ solver.apply(v,b,stat);
+ helper.mask(v);
+ Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,v);
+ if (gfs.gridView().comm().size()>1)
+ gfs.gridView().communicate(adddh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
+ }
+
+ /*!
+ \brief Clean up.
+ */
+ virtual void post (X& x) {}
+
+ private:
+ const GFS& gfs;
+ const M& A;
+ Dune::SuperLU<ISTLM> solver;
+ const ParallelISTLHelper<GFS>& helper;
+ };
+#endif
+
+ template<typename GFS>
+ class OVLPScalarProductImplementation
+ {
+ public:
+ OVLPScalarProductImplementation(const GFS& gfs_)
+ : gfs(gfs_), helper(gfs_)
+ {}
+
+ /*! \brief Dot product of two vectors.
+ It is assumed that the vectors are consistent on the interior+border
+ partition.
+ */
+ template<typename X>
+ typename X::ElementType dot (const X& x, const X& y) const
+ {
+ // do local scalar product on unique partition
+ typename X::ElementType sum = 0;
+ for (typename X::size_type i=0; i<x.base().N(); ++i)
+ for (typename X::size_type j=0; j<x[i].N(); ++j)
+ sum += (x[i][j]*y[i][j])*helper.mask(i,j);
+
+ // do global communication
+ return gfs.gridView().comm().sum(sum);
+ }
+
+ /*! \brief Norm of a right-hand side vector.
+ The vector must be consistent on the interior+border partition
+ */
+ template<typename X>
+ typename X::ElementType norm (const X& x) const
+ {
+ return sqrt(static_cast<double>(this->dot(x,x)));
+ }
+
+ const ParallelISTLHelper<GFS>& parallelHelper()
+ {
+ return helper;
+ }
+
+ private:
+ const GFS& gfs;
+ ParallelISTLHelper<GFS> helper;
+ };
+
+
+ template<typename GFS, typename X>
+ class OVLPScalarProduct
+ : public ScalarProduct<X>
+ {
+ public:
+ enum {category=Dune::SolverCategory::overlapping};
+ OVLPScalarProduct(const OVLPScalarProductImplementation<GFS>& implementation_)
+ : implementation(implementation_)
+ {}
+ virtual typename X::ElementType dot(const X& x, const X& y)
+ {
+ return implementation.dot(x,y);
+ }
+
+ virtual typename X::ElementType norm (const X& x)
+ {
+ return sqrt(static_cast<double>(this->dot(x,x)));
+ }
+
+ private:
+ const OVLPScalarProductImplementation<GFS>& implementation;
+ };
+
+ template<class GFS, class C,
+ template<class,class,class,int> class Preconditioner,
+ template<class> class Solver>
+ class ISTLBackend_OVLP_Base
+ : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] c_ a constraints object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] steps_ number of SSOR steps to apply as inner iteration
+ \param[in] verbose_ print messages if true
+ */
+ ISTLBackend_OVLP_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
+ int steps_=5, int verbose_=1)
+ : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), steps(steps_), verbose(verbose_)
+ {}
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename V::ElementType reduction)
+ {
+ typedef Dune::PDELab::OverlappingOperator<C,M,V,W> POP;
+ POP pop(c,A);
+ typedef OVLPScalarProduct<GFS,V> PSP;
+ PSP psp(*this);
+ typedef Preconditioner<M,V,W,1> SeqPrec;
+ SeqPrec seqprec(A,steps,1.0);
+ typedef Dune::PDELab::OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
+ WPREC wprec(gfs,seqprec,c,this->parallelHelper());
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z,r,stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+ private:
+ const GFS& gfs;
+ const C& c;
+ unsigned maxiter;
+ int steps;
+ int verbose;
+ };
+
+ // Base class for ILU0 as preconditioner
+ template<class GFS, class C,
+ template<class> class Solver>
+ class ISTLBackend_OVLP_ILU0_Base
+ : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] c_ a constraints object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ ISTLBackend_OVLP_ILU0_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000, int verbose_=1)
+ : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
+ {}
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename V::ElementType reduction)
+ {
+ typedef Dune::PDELab::OverlappingOperator<C,M,V,W> POP;
+ POP pop(c,A);
+ typedef OVLPScalarProduct<GFS,V> PSP;
+ PSP psp(*this);
+ typedef SeqILU0<M,V,W,1> SeqPrec;
+ SeqPrec seqprec(A,1.0);
+ typedef Dune::PDELab::OverlappingWrappedPreconditioner<C,GFS,SeqPrec> WPREC;
+ WPREC wprec(gfs,seqprec,c,this->parallelHelper());
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ Solver<V> solver(pop,psp,wprec,reduction,maxiter,verb);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z,r,stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+ private:
+ const GFS& gfs;
+ const C& c;
+ unsigned maxiter;
+ int steps;
+ int verbose;
+ };
+
+ //! \addtogroup PDELab_ovlpsolvers Overlapping Solvers
+ //! \{
+
+ /**
+ * @brief Overlapping parallel BiCGStab solver with SSOR preconditioner
+ * @tparam GFS The Type of the GridFunctionSpace.
+ * @tparam CC The Type of the Constraints Container.
+ */
+ template<class GFS, class CC>
+ class ISTLBackend_OVLP_BCGS_SSORk
+ : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs a grid function space
+ \param[in] cc a constraints container object
+ \param[in] maxiter maximum number of iterations to do
+ \param[in] steps number of SSOR steps to apply as inner iteration
+ \param[in] verbose print messages if true
+ */
+ ISTLBackend_OVLP_BCGS_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
+ int steps=5, int verbose=1)
+ : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::BiCGSTABSolver>(gfs, cc, maxiter, steps, verbose)
+ {}
+ };
+ /**
+ * @brief Overlapping parallel BiCGStab solver with ILU0 preconditioner
+ * @tparam GFS The Type of the GridFunctionSpace.
+ * @tparam CC The Type of the Constraints Container.
+ */
+ template<class GFS, class CC>
+ class ISTLBackend_OVLP_BCGS_ILU0
+ : public ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs a grid function space
+ \param[in] cc a constraints container object
+ \param[in] maxiter maximum number of iterations to do
+ \param[in] verbose print messages if true
+ */
+ ISTLBackend_OVLP_BCGS_ILU0 (const GFS& gfs, const CC& cc, unsigned maxiter=5000, int verbose=1)
+ : ISTLBackend_OVLP_ILU0_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs, cc, maxiter, verbose)
+ {}
+ };
+ /**
+ * @brief Overlapping parallel CGS solver with SSOR preconditioner
+ * @tparam GFS The Type of the GridFunctionSpace.
+ * @tparam CC The Type of the Constraints Container.
+ */
+ template<class GFS, class CC>
+ class ISTLBackend_OVLP_CG_SSORk
+ : public ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs a grid function space
+ \param[in] cc a constraints container object
+ \param[in] maxiter maximum number of iterations to do
+ \param[in] steps number of SSOR steps to apply as inner iteration
+ \param[in] verbose print messages if true
+ */
+ ISTLBackend_OVLP_CG_SSORk (const GFS& gfs, const CC& cc, unsigned maxiter=5000,
+ int steps=5, int verbose=1)
+ : ISTLBackend_OVLP_Base<GFS,CC,Dune::SeqSSOR, Dune::CGSolver>(gfs, cc, maxiter, steps, verbose)
+ {}
+ };
+
+ //! \} Solver
+
+ template<class GFS, class C, template<typename> class Solver>
+ class ISTLBackend_OVLP_SuperLU_Base
+ : public OVLPScalarProductImplementation<GFS>, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] c_ a constraints object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ ISTLBackend_OVLP_SuperLU_Base (const GFS& gfs_, const C& c_, unsigned maxiter_=5000,
+ int verbose_=1)
+ : OVLPScalarProductImplementation<GFS>(gfs_), gfs(gfs_), c(c_), maxiter(maxiter_), verbose(verbose_)
+ {}
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename V::ElementType reduction)
+ {
+ typedef Dune::PDELab::OverlappingOperator<C,M,V,W> POP;
+ POP pop(c,A);
+ typedef OVLPScalarProduct<GFS,V> PSP;
+ PSP psp(*this);
+#if HAVE_SUPERLU
+ typedef Dune::PDELab::SuperLUSubdomainSolver<GFS,M,V,W> PREC;
+ PREC prec(gfs,A);
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ Solver<V> solver(pop,psp,prec,reduction,maxiter,verbose);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z,r,stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+#else
+ std::cout << "No superLU support, please install and configure it." << std::endl;
+#endif
+ }
+
+ private:
+ const GFS& gfs;
+ const C& c;
+ unsigned maxiter;
+ int verbose;
+ };
+
+ //! \addtogroup PDELab_ovlpsolvers Overlapping Solvers
+ //! \{
+ /**
+ * @brief Overlapping parallel BiCGStab solver with SuperLU preconditioner
+ * @tparam GFS The Type of the GridFunctionSpace.
+ * @tparam CC The Type of the Constraints Container.
+ */
+ template<class GFS, class CC>
+ class ISTLBackend_OVLP_BCGS_SuperLU
+ : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>
+ {
+ public:
+
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] cc_ a constraints container object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ ISTLBackend_OVLP_BCGS_SuperLU (const GFS& gfs_, const CC& cc_, unsigned maxiter_=5000,
+ int verbose_=1)
+ : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::BiCGSTABSolver>(gfs_,cc_,maxiter_,verbose_)
+ {}
+ };
+
+ /**
+ * @brief Overlapping parallel CG solver with SuperLU preconditioner
+ * @tparam GFS The Type of the GridFunctionSpace.
+ * @tparam CC The Type of the Constraints Container.
+ */
+ template<class GFS, class CC>
+ class ISTLBackend_OVLP_CG_SuperLU
+ : public ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>
+ {
+ public:
+
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ \param[in] cc_ a constraints object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ ISTLBackend_OVLP_CG_SuperLU (const GFS& gfs_, const CC& cc_,
+ unsigned maxiter_=5000,
+ int verbose_=1)
+ : ISTLBackend_OVLP_SuperLU_Base<GFS,CC,Dune::CGSolver>(gfs_,cc_,maxiter_,verbose_)
+ {}
+ };
+
+
+ /** @brief Solver to be used for explicit time-steppers with (block-)diagonal mass matrix
+ * @tparam GFS The Type of the GridFunctionSpace.
+ */
+ template<class GFS>
+ class ISTLBackend_OVLP_ExplicitDiagonal
+ : public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] gfs_ a grid function space
+ */
+ explicit ISTLBackend_OVLP_ExplicitDiagonal (const GFS& gfs_)
+ : gfs(gfs_)
+ {}
+
+ explicit ISTLBackend_OVLP_ExplicitDiagonal (const ISTLBackend_OVLP_ExplicitDiagonal& other_)
+ : gfs(other_.gfs)
+ {}
+
+ /*! \brief compute global norm of a vector
+
+ \param[in] v the given vector
+ */
+ template<class V>
+ typename V::ElementType norm(const V& v) const
+ {
+ dune_static_assert
+ (AlwaysFalse<V>::value,
+ "ISTLBackend_OVLP_ExplicitDiagonal::norm() should not be "
+ "neccessary, so we skipped the implementation. If you have a "
+ "scenario where you need it, please implement it or report back to "
+ "us.");
+ }
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename W::ElementType reduction)
+ {
+ Dune::SeqJac<M,V,W> jac(A,1,1.0);
+ jac.pre(z,r);
+ jac.apply(z,r);
+ jac.post(z);
+ if (gfs.gridView().comm().size()>1)
+ {
+ Dune::PDELab::CopyDataHandle<GFS,V> copydh(gfs,z);
+ gfs.gridView().communicate(copydh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
+ }
+ res.converged = true;
+ res.iterations = 1;
+ res.elapsed = 0.0;
+ res.reduction = reduction;
+ res.conv_rate = reduction; // pow(reduction,1.0/1)
+ }
+
+ private:
+ const GFS& gfs;
+ };
+ //! \} Overlapping Solvers
+
+ /**
+ * @brief traits class for using AMG with the old grid operator space
+ * instead of the new grid operator
+ * @tparam MatrixFT The field type of the matrix.
+ * @tparam V The type of the BackendVectorSelector
+ * @tparam GOS The type of the old grid operator space.
+ */
+ template<class MatrixFT, class V, class GOS>
+ class fakeGOTraits
+ {
+ public:
+ struct Traits{
+ typedef V Domain;
+ typedef MatrixFT JacobianField;
+ typedef typename V::ElementType DomainField;
+ typedef typename GOS::Traits GOSTraits;
+ typedef typename GOSTraits::TrialGridFunctionSpace TrialGridFunctionSpace;
+ typedef typename GOS::template MatrixContainer<MatrixFT>::Type Jacobian;
+ };
+ };
+
+ template<class GO, int s, template<class,class,class,int> class Preconditioner,
+ template<class> class Solver>
+ class ISTLBackend_AMG : public LinearResultStorage
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+ typedef typename Dune::PDELab::ParallelISTLHelper<GFS> PHELPER;
+ typedef typename GO::Traits::Jacobian M;
+ typedef typename M::BaseT MatrixType;
+ typedef typename GO::Traits::Domain V;
+ typedef typename BlockProcessor<GFS>::template AMGVectorTypeSelector<V>::Type VectorType;
+ typedef typename CommSelector<s,Dune::MPIHelper::isFake>::type Comm;
+#if HAVE_MPI
+ typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
+ typedef Dune::BlockPreconditioner<VectorType,VectorType,Comm,Smoother> ParSmoother;
+ typedef Dune::OverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator;
+#else
+ typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother;
+ typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
+#endif
+ typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
+ typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG;
+
+ public:
+
+ /**
+ * @brief Parameters object to customize matrix hierachy building.
+ */
+ typedef Dune::Amg::Parameters Parameters;
+
+ ISTLBackend_AMG(const GFS& gfs_, unsigned maxiter_=5000,
+ int verbose_=1, bool reuse_=false,
+ bool usesuperlu_=true)
+ : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu),
+ verbose(verbose_), reuse(reuse_), firstapply(true),
+ usesuperlu(usesuperlu_)
+ {
+ params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
+ params.setDebugLevel(verbose_);
+#if !HAVE_SUPERLU
+ if (gfs.gridView().comm().rank() == 0 && usesuperlu == true)
+ {
+ std::cout << "WARNING: You are using AMG without SuperLU!"
+ << " Please consider installing SuperLU,"
+ << " or set the usesuperlu flag to false"
+ << " to suppress this warning." << std::endl;
+ }
+#endif
+ }
+
+ /*! \brief set AMG parameters
+
+ \param[in] params_ a parameter object of Type Dune::Amg::Parameters
+ */
+ void setParameters(const Parameters& params_)
+ {
+ params = params_;
+ }
+
+ void setparams(Parameters params_)
+ {
+ params = params_;
+ }
+
+ /**
+ * @brief Get the parameters describing the behaviuour of AMG.
+ *
+ * The returned object can be adjusted to ones needs and then can be
+ * reset using setParameters.
+ * @return The object holding the parameters of AMG.
+ */
+ const Parameters& parameters() const
+ {
+ return params;
+ }
+
+ /*! \brief compute global norm of a vector
+
+ \param[in] v the given vector
+ */
+ typename V::ElementType norm (const V& v) const
+ {
+ typedef Dune::PDELab::OverlappingScalarProduct<GFS,V> PSP;
+ PSP psp(gfs,phelper);
+ return psp.norm(v);
+ }
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template <class Vector>
+ void apply(MatrixType& A, Vector& z, Vector& r, typename V::ElementType reduction)
+ {
+ Timer watch;
+ Comm oocc(gfs.gridView().comm());
+ MatrixType& mat=A;
+ typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
+ Dune::Amg::FirstDiagonal> > Criterion;
+#if HAVE_MPI
+ phelper.createIndexSetAndProjectForAMG(mat, oocc);
+ Operator oop(mat, oocc);
+ Dune::OverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc);
+#else
+ Operator oop(mat);
+ Dune::SeqScalarProduct<VectorType> sp;
+#endif
+ SmootherArgs smootherArgs;
+ smootherArgs.iterations = 1;
+ smootherArgs.relaxationFactor = 1;
+ Criterion criterion(params);
+ stats.tprepare=watch.elapsed();
+ watch.reset();
+
+ int verb=0;
+ if (gfs.gridView().comm().rank()==0) verb=verbose;
+ //only construct a new AMG if the matrix changes
+ if (reuse==false || firstapply==true){
+ amg.reset(new AMG(oop, criterion, smootherArgs, oocc));
+ firstapply = false;
+ stats.tsetup = watch.elapsed();
+ stats.levels = amg->maxlevels();
+ stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
+ }
+ watch.reset();
+ Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb);
+ Dune::InverseOperatorResult stat;
+
+// solver.apply(BlockProcessor<GFS>::getVector(z),BlockProcessor<GFS>::getVector(r),stat);
+ solver.apply(z,r,stat);
+ stats.tsolve= watch.elapsed();
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ /**
+ * @brief Get statistics of the AMG solver (no of levels, timings).
+ * @return statistis of the AMG solver.
+ */
+ const ISTLAMGStatistics& statistics() const
+ {
+ return stats;
+ }
+
+ private:
+ const GFS& gfs;
+ PHELPER phelper;
+ unsigned maxiter;
+ Parameters params;
+ int verbose;
+ bool reuse;
+ bool firstapply;
+ bool usesuperlu;
+ Dune::shared_ptr<AMG> amg;
+ ISTLAMGStatistics stats;
+ };
+
+ //! \addtogroup PDELab_ovlpsolvers Overlapping Solvers
+ //! \{
+
+ /**
+ * @brief Overlapping parallel conjugate gradient solver preconditioned with AMG smoothed by SSOR
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ * @tparam s The bits to use for the global index.
+ */
+ template<class GO, int s=96>
+ class ISTLBackend_CG_AMG_SSOR
+ : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+ public:
+ /**
+ * @brief Constructor
+ * @param gfs_ The grid function space used.
+ * @param maxiter_ The maximum number of iterations allowed.
+ * @param verbose_ The verbosity level to use.
+ * @param reuse_ Set true, if the Matrix to be used is always identical
+ * (AMG aggregation is then only performed once).
+ * @param usesuperlu_ Set false, to suppress the no SuperLU warning
+ */
+ ISTLBackend_CG_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
+ int verbose_=1, bool reuse_=false,
+ bool usesuperlu_=true)
+ : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::CGSolver>
+ (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ /**
+ * @brief Overlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR.
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ * @tparam s The bits to use for the globale index.
+ */
+ template<class GO, int s=96>
+ class ISTLBackend_BCGS_AMG_SSOR
+ : public ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+ public:
+ /**
+ * @brief Constructor
+ * @param gfs_ The grid function space used.
+ * @param maxiter_ The maximum number of iterations allowed.
+ * @param verbose_ The verbosity level to use.
+ * @param reuse_ Set true, if the Matrix to be used is always identical
+ * (AMG aggregation is then only performed once).
+ * @param usesuperlu_ Set false, to suppress the no SuperLU warning
+ */
+ ISTLBackend_BCGS_AMG_SSOR(const GFS& gfs_, unsigned maxiter_=5000,
+ int verbose_=1, bool reuse_=false,
+ bool usesuperlu_=true)
+ : ISTLBackend_AMG<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ (gfs_, maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ //! \} Overlapping Solvers
+
+ //! \} group Backend
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/linear/seqistlsolverbackend.hh b/dumux/linear/seqistlsolverbackend.hh
new file mode 100644
index 0000000000..3aa5e06254
--- /dev/null
+++ b/dumux/linear/seqistlsolverbackend.hh
@@ -0,0 +1,743 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=8 sw=2 sts=2:
+/*****************************************************************************
+ * This file is based on the file of DUNE-PDELab with the same name, *
+ * http://www.dune-project.org/pdelab/. *
+ * Some changes are made such that the backends work for the vectors and *
+ * matrices used in Dumux. *
+ *****************************************************************************/
+#ifndef DUNE_SEQISTLSOLVERBACKEND_HH
+#define DUNE_SEQISTLSOLVERBACKEND_HH
+
+#include <dune/common/deprecated.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/solvercategory.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/scalarproducts.hh>
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/pinfo.hh>
+#include <dune/istl/io.hh>
+#include <dune/istl/superlu.hh>
+
+#include <dune/pdelab/constraints/constraints.hh>
+#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
+#include <dune/pdelab/backend/solver.hh>
+#include <dune/pdelab/backend/parallelistlhelper.hh>
+
+#include "istlvectorbackend.hh"
+
+namespace Dumux {
+
+ //! \addtogroup Backend
+ //! \ingroup PDELab
+ //! \{
+
+ template<typename X, typename Y, typename GOS>
+ class OnTheFlyOperator : public Dune::LinearOperator<X,Y>
+ {
+ public:
+ typedef X domain_type;
+ typedef Y range_type;
+ typedef typename X::field_type field_type;
+
+ enum {category=Dune::SolverCategory::sequential};
+
+ OnTheFlyOperator (GOS& gos_)
+ : gos(gos_)
+ {}
+
+ virtual void apply (const X& x, Y& y) const
+ {
+ y = 0.0;
+ gos.jacobian_apply(x,y);
+ }
+
+ virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+ {
+ Y temp(y);
+ temp = 0.0;
+ gos.jacobian_apply(x,temp);
+ y.axpy(alpha,temp);
+ }
+
+ private:
+ GOS& gos;
+ };
+
+ //==============================================================================
+ // Here we add some standard linear solvers conforming to the linear solver
+ // interface required to solve linear and nonlinear problems.
+ //==============================================================================
+
+ template<template<class,class,class,int> class Preconditioner,
+ template<class> class Solver>
+ class ISTLBackend_SEQ_Base
+ : public SequentialNorm, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_Base(unsigned maxiter_=5000, int verbose_=1)
+ : maxiter(maxiter_), verbose(verbose_)
+ {}
+
+
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename W::ElementType reduction)
+ {
+ Dune::MatrixAdapter<typename M::BaseT,
+ typename V::BaseT,
+ typename W::BaseT> opa(A);
+ Preconditioner<typename M::BaseT,
+ typename V::BaseT,
+ typename W::BaseT,1> ssor(A, 3, 1.0);
+ Solver<typename V::BaseT> solver(opa, ssor, reduction, maxiter, verbose);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z, r, stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ private:
+ unsigned maxiter;
+ int verbose;
+ };
+
+ template<template<typename> class Solver>
+ class ISTLBackend_SEQ_ILU0
+ : public SequentialNorm, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_ILU0 (unsigned maxiter_=5000, int verbose_=1)
+ : maxiter(maxiter_), verbose(verbose_)
+ {}
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction)
+ {
+ Dune::MatrixAdapter<typename M::BaseT,
+ typename V::BaseT,
+ typename W::BaseT> opa(A);
+ Dune::SeqILU0<typename M::BaseT,
+ typename V::BaseT,
+ typename W::BaseT> ilu0(A, 1.0);
+ Solver<typename V::BaseT> solver(opa, ilu0, reduction, maxiter, verbose);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z, r, stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+ private:
+ unsigned maxiter;
+ int verbose;
+ };
+
+ template<template<typename> class Solver>
+ class ISTLBackend_SEQ_ILUn
+ : public SequentialNorm, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+ \param[in] n The number of levels to be used.
+ \param[in] w The relaxation factor.
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ ISTLBackend_SEQ_ILUn (int n, double w, unsigned maxiter_=5000, int verbose_=1)
+ : n_(n), w_(w), maxiter(maxiter_), verbose(verbose_)
+ {}
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename W::ElementType reduction)
+ {
+ Dune::MatrixAdapter<M,V,W> opa(A);
+ Dune::SeqILUn<typename M::BaseT,V,W> ilu0(A.base(), n_, w_);
+ Solver<V> solver(opa, ilu0, reduction, maxiter, verbose);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z, r, stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+ private:
+ int n_;
+ double w_;
+
+ unsigned maxiter;
+ int verbose;
+ };
+
+ //! \addtogroup PDELab_seqsolvers Sequential Solvers
+ //! \{
+
+ /**
+ * @brief Backend for sequential loop solver with Jacobi preconditioner.
+ */
+ class ISTLBackend_SEQ_LOOP_Jac
+ : public ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::LoopSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_LOOP_Jac (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::LoopSolver>(maxiter_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Backend for sequential BiCGSTAB solver with Jacobi preconditioner.
+ */
+ class ISTLBackend_SEQ_BCGS_Jac
+ : public ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::BiCGSTABSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_BCGS_Jac (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::BiCGSTABSolver>(maxiter_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Backend for sequential BiCGSTAB solver with SSOR preconditioner.
+ */
+ class ISTLBackend_SEQ_BCGS_SSOR
+ : public ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_BCGS_SSOR (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::BiCGSTABSolver>(maxiter_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Backend for sequential BiCGSTAB solver with ILU0 preconditioner.
+ */
+ class ISTLBackend_SEQ_BCGS_ILU0
+ : public ISTLBackend_SEQ_ILU0<Dune::BiCGSTABSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_BCGS_ILU0 (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_ILU0<Dune::BiCGSTABSolver>(maxiter_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Backend for sequential conjugate gradient solver with ILU0 preconditioner.
+ */
+ class ISTLBackend_SEQ_CG_ILU0
+ : public ISTLBackend_SEQ_ILU0<Dune::CGSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_CG_ILU0 (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_ILU0<Dune::CGSolver>(maxiter_, verbose_)
+ {}
+ };
+
+ //! \brief Sequential BiCGStab solver with ILU0 preconditioner
+ class ISTLBackend_SEQ_BCGS_ILUn
+ : public ISTLBackend_SEQ_ILUn<Dune::BiCGSTABSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+
+ \param[in] n_ The number of levels to be used.
+ \param[in] w_ The relaxation factor.
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_BCGS_ILUn (int n_, double w_=1.0, unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_ILUn<Dune::BiCGSTABSolver>(n_, w_, maxiter_, verbose_)
+ {}
+ };
+
+ //! \brief Sequential congute gradient solver with ILU0 preconditioner
+ class ISTLBackend_SEQ_CG_ILUn
+ : public ISTLBackend_SEQ_ILUn<Dune::CGSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+
+ \param[in] n_ The number of levels to be used.
+ \param[in] w_ The relaxation factor.
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_CG_ILUn (int n_, double w_=1.0, unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_ILUn<Dune::CGSolver>(n_, w_, maxiter_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Backend for sequential conjugate gradient solver with SSOR preconditioner.
+ */
+ class ISTLBackend_SEQ_CG_SSOR
+ : public ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::CGSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_CG_SSOR (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::CGSolver>(maxiter_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Backend using a MINRes solver preconditioned by SSOR.
+ */
+ class ISTLBackend_SEQ_MINRES_SSOR
+ : public ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::MINRESSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_MINRES_SSOR (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_Base<Dune::SeqSSOR, Dune::MINRESSolver>(maxiter_, verbose_)
+ {}
+ };
+
+ /**
+ * @brief Backend for conjugate gradient solver with Jacobi preconditioner.
+ */
+ class ISTLBackend_SEQ_CG_Jac
+ : public ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::CGSolver>
+ {
+ public:
+ /*! \brief make a linear solver object
+ \param[in] maxiter_ maximum number of iterations to do
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_CG_Jac (unsigned maxiter_=5000, int verbose_=1)
+ : ISTLBackend_SEQ_Base<Dune::SeqJac, Dune::CGSolver>(maxiter_, verbose_)
+ {}
+ };
+
+#if HAVE_SUPERLU
+ /**
+ * @brief Solver backend using SuperLU as a direct solver.
+ */
+ class ISTLBackend_SEQ_SuperLU
+ : public SequentialNorm, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+
+ \param[in] verbose_ print messages if true
+ */
+ explicit ISTLBackend_SEQ_SuperLU (int verbose_=1)
+ : verbose(verbose_)
+ {}
+
+
+ /*! \brief make a linear solver object
+
+ \param[in] maxiter Maximum number of allowed steps (ignored)
+ \param[in] verbose_ print messages if true
+ */
+ ISTLBackend_SEQ_SuperLU (int maxiter, int verbose_)
+ : verbose(verbose_)
+ {}
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename W::ElementType reduction)
+ {
+ typedef typename M::BaseT ISTLM;
+ Dune::SuperLU<ISTLM> solver(A, verbose);
+ Dune::InverseOperatorResult stat;
+ solver.apply(z, r, stat);
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+ private:
+ int verbose;
+ };
+#endif // HAVE_SUPERLU
+
+ //! Solver to be used for explicit time-steppers with (block-)diagonal mass matrix
+ class ISTLBackend_SEQ_ExplicitDiagonal
+ : public SequentialNorm, public LinearResultStorage
+ {
+ public:
+ /*! \brief make a linear solver object
+ */
+ ISTLBackend_SEQ_ExplicitDiagonal ()
+ {}
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ template<class M, class V, class W>
+ void apply(M& A, V& z, W& r, typename W::ElementType reduction)
+ {
+ Dune::SeqJac<typename M::BaseT,
+ typename V::BaseT,
+ typename W::BaseT> jac(A,1,1.0);
+ jac.pre(z,r);
+ jac.apply(z,r);
+ jac.post(z);
+ res.converged = true;
+ res.iterations = 1;
+ res.elapsed = 0.0;
+ res.reduction = reduction;
+ res.conv_rate = reduction; // pow(reduction,1.0/1)
+ }
+ };
+
+ //! \} Sequential Solvers
+
+ /**
+ * @brief Class providing some statistics of the AMG solver.
+ *
+ */
+ struct ISTLAMGStatistics
+ {
+ /**
+ * @brief The needed for computing the parallel information and
+ * for adapting the linear system.
+ */
+ double tprepare;
+ /** @brief the number of levels in the AMG hierarchy. */
+ int levels;
+ /** @brief The time spent in solving the system (without building the hierarchy. */
+ double tsolve;
+ /** @brief The time needed for building the AMG hierarchy (coarsening). */
+ double tsetup;
+ /** @brief The number of iterations performed until convergence was reached. */
+ int iterations;
+ /** @brief True if a direct solver was used on the coarset level. */
+ bool directCoarseLevelSolver;
+ };
+
+ template<class GO, template<class,class,class,int> class Preconditioner, template<class> class Solver>
+ class ISTLBackend_SEQ_AMG : public LinearResultStorage
+ {
+ typedef typename GO::Traits::TrialGridFunctionSpace GFS;
+ typedef typename GO::Traits::Jacobian M;
+ typedef typename M::BaseT MatrixType;
+ typedef typename GO::Traits::Domain V;
+ typedef typename BlockProcessor<GFS>::template AMGVectorTypeSelector<V>::Type VectorType;
+ typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother;
+ typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator;
+ typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+ typedef Dune::Amg::AMG<Operator,VectorType,Smoother> AMG;
+ typedef Dune::Amg::Parameters Parameters;
+
+ public:
+ ISTLBackend_SEQ_AMG(unsigned maxiter_=5000, int verbose_=1,
+ bool reuse_=false, bool usesuperlu_=true)
+ : maxiter(maxiter_), params(15,2000), verbose(verbose_),
+ reuse(reuse_), firstapply(true), usesuperlu(usesuperlu_)
+ {
+ params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
+ params.setDebugLevel(verbose_);
+#if !HAVE_SUPERLU
+ if (usesuperlu == true)
+ {
+ std::cout << "WARNING: You are using AMG without SuperLU!"
+ << " Please consider installing SuperLU,"
+ << " or set the usesuperlu flag to false"
+ << " to suppress this warning." << std::endl;
+ }
+#endif
+ }
+
+ /*! \brief set AMG parameters
+
+ \param[in] params_ a parameter object of Type Dune::Amg::Parameters
+ */
+ void setparams(Parameters params_)
+ {
+ params = params_;
+ }
+
+ /*! \brief compute global norm of a vector
+
+ \param[in] v the given vector
+ */
+ typename V::ElementType norm (const V& v) const
+ {
+ return v.base().two_norm();
+ }
+
+ /*! \brief solve the given linear system
+
+ \param[in] A the given matrix
+ \param[out] z the solution vector to be computed
+ \param[in] r right hand side
+ \param[in] reduction to be achieved
+ */
+ void apply(MatrixType& A, V& z, V& r, typename V::ElementType reduction)
+ {
+ Timer watch;
+ MatrixType& mat=A;
+ typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType,
+ Dune::Amg::FirstDiagonal> > Criterion;
+ SmootherArgs smootherArgs;
+ smootherArgs.iterations = 1;
+ smootherArgs.relaxationFactor = 1;
+
+ Criterion criterion(params);
+ Operator oop(mat);
+ //only construct a new AMG if the matrix changes
+ if (reuse==false || firstapply==true){
+ amg.reset(new AMG(oop, criterion, smootherArgs));
+ firstapply = false;
+ stats.tsetup = watch.elapsed();
+ stats.levels = amg->maxlevels();
+ stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver();
+ }
+ watch.reset();
+ Dune::InverseOperatorResult stat;
+
+ Solver<VectorType> solver(oop,*amg,reduction,maxiter,verbose);
+ solver.apply(BlockProcessor<GFS>::getVector(z),BlockProcessor<GFS>::getVector(r),stat);
+ stats.tsolve= watch.elapsed();
+ res.converged = stat.converged;
+ res.iterations = stat.iterations;
+ res.elapsed = stat.elapsed;
+ res.reduction = stat.reduction;
+ res.conv_rate = stat.conv_rate;
+ }
+
+
+ /**
+ * @brief Get statistics of the AMG solver (no of levels, timings).
+ * @return statistis of the AMG solver.
+ */
+ const ISTLAMGStatistics& statistics() const
+ {
+ return stats;
+ }
+
+ private:
+ unsigned maxiter;
+ Parameters params;
+ int verbose;
+ bool reuse;
+ bool firstapply;
+ bool usesuperlu;
+ Dune::shared_ptr<AMG> amg;
+ ISTLAMGStatistics stats;
+ };
+
+ //! \addtogroup PDELab_seqsolvers Sequential Solvers
+ //! \{
+
+ /**
+ * @brief Sequential conjugate gradient solver preconditioned with AMG smoothed by SSOR
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ */
+ template<class GO>
+ class ISTLBackend_SEQ_CG_AMG_SSOR
+ : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::CGSolver>
+ {
+
+ public:
+ /**
+ * @brief Constructor
+ * @param maxiter_ The maximum number of iterations allowed.
+ * @param verbose_ The verbosity level to use.
+ * @param reuse_ Set true, if the Matrix to be used is always identical
+ * (AMG aggregation is then only performed once).
+ * @param usesuperlu_ Set false, to suppress the no SuperLU warning
+ */
+ ISTLBackend_SEQ_CG_AMG_SSOR(unsigned maxiter_=5000, int verbose_=1,
+ bool reuse_=false, bool usesuperlu_=true)
+ : ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::CGSolver>
+ (maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ /**
+ * @brief Sequential BiCGStab solver preconditioned with AMG smoothed by SSOR
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ */
+ template<class GO>
+ class ISTLBackend_SEQ_BCGS_AMG_SSOR
+ : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ {
+
+ public:
+ /**
+ * @brief Constructor
+ * @param maxiter_ The maximum number of iterations allowed.
+ * @param verbose_ The verbosity level to use.
+ * @param reuse_ Set true, if the Matrix to be used is always identical
+ * (AMG aggregation is then only performed once).
+ * @param usesuperlu_ Set false, to suppress the no SuperLU warning
+ */
+ ISTLBackend_SEQ_BCGS_AMG_SSOR(unsigned maxiter_=5000, int verbose_=1,
+ bool reuse_=false, bool usesuperlu_=true)
+ : ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::BiCGSTABSolver>
+ (maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ /**
+ * @brief Sequential BiCGSTAB solver preconditioned with AMG smoothed by SOR
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ */
+ template<class GO>
+ class ISTLBackend_SEQ_BCGS_AMG_SOR
+ : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::BiCGSTABSolver>
+ {
+
+ public:
+ /**
+ * @brief Constructor
+ * @param maxiter_ The maximum number of iterations allowed.
+ * @param verbose_ The verbosity level to use.
+ * @param reuse_ Set true, if the Matrix to be used is always identical
+ * (AMG aggregation is then only performed once).
+ * @param usesuperlu_ Set false, to suppress the no SuperLU warning
+ */
+ ISTLBackend_SEQ_BCGS_AMG_SOR(unsigned maxiter_=5000, int verbose_=1,
+ bool reuse_=false, bool usesuperlu_=true)
+ : ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::BiCGSTABSolver>
+ (maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ /**
+ * @brief Sequential Loop solver preconditioned with AMG smoothed by SSOR
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ */
+ template<class GO>
+ class ISTLBackend_SEQ_LS_AMG_SSOR
+ : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::LoopSolver>
+ {
+
+ public:
+ /**
+ * @brief Constructor
+ * @param maxiter_ The maximum number of iterations allowed.
+ * @param verbose_ The verbosity level to use.
+ * @param reuse_ Set true, if the Matrix to be used is always identical
+ * (AMG aggregation is then only performed once).
+ * @param usesuperlu_ Set false, to suppress the no SuperLU warning
+ */
+ ISTLBackend_SEQ_LS_AMG_SSOR(unsigned maxiter_=5000, int verbose_=1,
+ bool reuse_=false, bool usesuperlu_=true)
+ : ISTLBackend_SEQ_AMG<GO, Dune::SeqSSOR, Dune::LoopSolver>
+ (maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ /**
+ * @brief Sequential Loop solver preconditioned with AMG smoothed by SOR
+ * @tparam GO The type of the grid operator
+ * (or the fakeGOTraits class for the old grid operator space).
+ */
+ template<class GO>
+ class ISTLBackend_SEQ_LS_AMG_SOR
+ : public ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::LoopSolver>
+ {
+
+ public:
+ /**
+ * @brief Constructor
+ * @param maxiter_ The maximum number of iterations allowed.
+ * @param verbose_ The verbosity level to use.
+ * @param reuse_ Set true, if the Matrix to be used is always identical
+ * (AMG aggregation is then only performed once).
+ * @param usesuperlu_ Set false, to suppress the no SuperLU warning
+ */
+ ISTLBackend_SEQ_LS_AMG_SOR(unsigned maxiter_=5000, int verbose_=1,
+ bool reuse_=false, bool usesuperlu_=true)
+ : ISTLBackend_SEQ_AMG<GO, Dune::SeqSOR, Dune::LoopSolver>
+ (maxiter_, verbose_, reuse_, usesuperlu_)
+ {}
+ };
+
+ //! \} group Sequential Solvers
+ //! \} group Backend
+
+} // namespace Dumux
+
+#endif
--
GitLab