From: Denis Davydov Date: Thu, 17 Nov 2016 09:39:55 +0000 (+0100) Subject: add MGInterfaceOperator and use it in a few tests X-Git-Tag: v8.5.0-rc1~377^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a94688bafba07de676a6daa993fae34b43a8c918;p=dealii.git add MGInterfaceOperator and use it in a few tests --- diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 87282f1ba8..668e408a46 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -245,6 +245,86 @@ namespace MatrixFreeOperators + /** + * Auxiliary class to provide interface vmult/Tvmult methods required in + * adaptive geometric multgrids. @p OperatorType class should be derived + * from MatrixFreeOperators::Base class. + * + * The adaptive multigrid realization in deal.II implements an approach + * called local smoothing. This means that the smoothing on the finest level + * only covers the local part of the mesh defined by the fixed (finest) grid + * level and ignores parts of the computational domain where the terminal + * cells are coarser than this level. As the method progresses to coarser + * levels, more and more of the global mesh will be covered. At some coarser + * level, the whole mesh will be covered. Since all level matrices in the + * multigrid method cover a single level in the mesh, no hanging nodes + * appear on the level matrices. At the interface between multigrid levels, + * homogeneous Dirichlet boundary conditions are set while smoothing. When + * the residual is transferred to the next coarser level, however, the + * coupling over the multigrid interface needs to be taken into account. + * This is done by the so-called interface (or edge) matrices that compute + * the part of the residual that is missed by the level matrix with + * homogeneous Dirichlet conditions. We refer to the @ref mg_paper + * "Multigrid paper by Janssen and Kanschat" for more details. + * + * For the implementation of those interface matrices, most infrastructure + * is already in place and provided by MatrixFreeOperators::Base through the + * two multiplication routines vmult_interface_down() and + * vmult_interface_up(). The only thing MGInterfaceOperator does is + * wrapping those operations and make them accessible via + * @p vmult() and @p Tvmult interface as expected by the multigrid routines + * (that were originally written for matrices, hence expecting those names). + * Note that the vmult_interface_down is used during the restriction phase of + * the multigrid V-cycle, whereas vmult_interface_up is used during the + * prolongation phase. + * + * @author Martin Kronbichler, 2016 + */ + template + class MGInterfaceOperator : public Subscriptor + { + public: + /** + * Number typedef. + */ + typedef typename OperatorType::value_type value_type; + + /** + * Default constructor. + */ + MGInterfaceOperator(); + + /** + * Clear the pointer to the OperatorType object. + */ + void clear(); + + /** + * Initialize this class with an operator @p operator_in. + */ + void initialize (const OperatorType &operator_in); + + /** + * vmult operator, see class description for more info. + */ + void vmult (LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &src) const; + + /** + * Tvmult operator, see class description for more info. + */ + void Tvmult (LinearAlgebra::distributed::Vector &dst, + const LinearAlgebra::distributed::Vector &src) const; + + private: + /** + * Const pointer to the operator class. + */ + SmartPointer mf_base_operator; + }; + + + /** * This class implements the operation of the action of the inverse of a * mass matrix on an element for the special case of an evaluation object @@ -965,6 +1045,62 @@ namespace MatrixFreeOperators + //------------------------- MGInterfaceOperator ------------------------------ + + template + MGInterfaceOperator::MGInterfaceOperator () + : + Subscriptor(), + mf_base_operator(NULL) + { + } + + + + template + void + MGInterfaceOperator::clear () + { + mf_base_operator = NULL; + } + + + + template + void + MGInterfaceOperator::initialize (const OperatorType &operator_in) + { + mf_base_operator = &operator_in; + } + + + + template + void + MGInterfaceOperator::vmult (LinearAlgebra::distributed::Vector::value_type> &dst, + const LinearAlgebra::distributed::Vector::value_type> &src) const + { + Assert(mf_base_operator != NULL, + ExcNotInitialized()); + + mf_base_operator->vmult_interface_down(dst, src); + } + + + + template + void + MGInterfaceOperator::Tvmult (LinearAlgebra::distributed::Vector::value_type> &dst, + const LinearAlgebra::distributed::Vector::value_type> &src) const + { + Assert(mf_base_operator != NULL, + ExcNotInitialized()); + + mf_base_operator->vmult_interface_up(dst, src); + } + + + //-----------------------------MassOperator---------------------------------- template diff --git a/tests/matrix_free/parallel_multigrid_adaptive_02.cc b/tests/matrix_free/parallel_multigrid_adaptive_02.cc index 1e99b8af18..ccbf07a91d 100644 --- a/tests/matrix_free/parallel_multigrid_adaptive_02.cc +++ b/tests/matrix_free/parallel_multigrid_adaptive_02.cc @@ -49,34 +49,6 @@ std::ofstream logfile("output"); using namespace dealii::MatrixFreeOperators; - -template -class MGInterfaceMatrix : public Subscriptor -{ -public: - void initialize (const LAPLACEOPERATOR &laplace) - { - this->laplace = &laplace; - } - - void vmult (LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src) const - { - laplace->vmult_interface_down(dst, src); - } - - void Tvmult (LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src) const - { - laplace->vmult_interface_up(dst, src); - } - -private: - SmartPointer laplace; -}; - - - template class MGTransferMF : public MGTransferMatrixFree { @@ -239,7 +211,7 @@ void do_test (const DoFHandler &dof) level); mg_matrices[level].compute_diagonal(); } - MGLevelObject > mg_interface_matrices; + MGLevelObject > mg_interface_matrices; mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1); for (unsigned int level=0; level -class MGInterfaceMatrix : public Subscriptor -{ -public: - void initialize (const LAPLACEOPERATOR &laplace) - { - this->laplace = &laplace; - } - - void vmult (LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src) const - { - laplace->vmult_interface_down(dst, src); - } - - void Tvmult (LinearAlgebra::distributed::Vector &dst, - const LinearAlgebra::distributed::Vector &src) const - { - laplace->vmult_interface_up(dst, src); - } - -private: - SmartPointer laplace; -}; - - - template class MGTransferMF : public MGTransferMatrixFree { @@ -238,7 +211,7 @@ void do_test (const DoFHandler &dof) level); mg_matrices[level].compute_diagonal(); } - MGLevelObject > mg_interface_matrices; + MGLevelObject > mg_interface_matrices; mg_interface_matrices.resize(0, dof.get_triangulation().n_global_levels()-1); for (unsigned int level=0; level