]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
experiments
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 13 Aug 2002 16:02:37 +0000 (16:02 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 13 Aug 2002 16:02:37 +0000 (16:02 +0000)
git-svn-id: https://svn.dealii.org/trunk@6324 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/fas.h [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/multigrid/fas.h b/deal.II/deal.II/include/multigrid/fas.h
new file mode 100644 (file)
index 0000000..d7e999a
--- /dev/null
@@ -0,0 +1,524 @@
+//----------------------------  multigrid.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  multigrid.h  ---------------------------
+#ifndef __deal2__fas_h
+#define __deal2__fas_h
+
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <lac/sparse_matrix.h>
+#include <lac/vector.h>
+#include <multigrid/mg_base.h>
+
+#include <vector>
+
+
+/**
+ * Implementation of the multigrid method.
+ *
+ * The function actually performing a multi-level cycle,
+ * @p{level_mgstep}, as well as the function @p{vcycle}, calling it,
+ * require several helper classes handed over as template parameters.
+ * These classes have to meet the following requirements:
+ *
+ * @p{MATRIX} is a matrix as specified for LAC-solvers, that is, it
+ * has a function
+ * \begin{verbatim}
+ *   void MATRIX::vmult(VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * performing the matrix vector product @p{x=Ab}.
+ *
+ * @p{SMOOTH} is a class with a function
+ * \begin{verbatim}
+ *   void SMOOTH::smooth (unsigned int level, VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * modifying the vector @p{x} such that the residual @{b-Ax}
+ * on level @p{l} is smoothened. Refer to @ref{MGSmootherRelaxation}
+ * for an example.
+ *
+ * @p{COARSE} has an
+ * \begin{verbatim}
+ *   void COARSE::operator() (VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * returning the solution to the system @p{Ax=b} on the coarsest level
+ * in @p{x}.
+ * 
+ * @author Guido Kanschat, 1999, 2001
+ */
+template <class VECTOR>
+class FAS : public Subscriptor
+{
+  public:
+  typedef VECTOR vector_type;
+  typedef const VECTOR const_vector_type;
+  
+                                    /**
+                                     * Constructor. The
+                                     * @p{MGDoFHandler} is used to
+                                     * determine the highest possible
+                                     * level. @p{transfer} is an
+                                     * object performing prolongation
+                                     * and restriction.
+                                     *
+                                     * The V-cycle will start on
+                                     * level @p{maxlevel} and goes
+                                     * down to level @p{minlevel},
+                                     * where the coarse grid solver
+                                     * will be used.
+                                     *
+                                     * This function already
+                                     * initializes the vectors which
+                                     * will be used later on in the
+                                     * course of the
+                                     * computations. You should
+                                     * therefore create objects of
+                                     * this type as late as possible.
+                                     */
+  template <int dim>
+    FAS(const MGDoFHandler<dim>& mg_dof_handler,
+             const unsigned int            minlevel = 0,
+             const unsigned int            maxlevel = 1000000);
+    
+
+                                    /**
+                                     * Execute one step of the
+                                     * v-cycle algorithm.  This
+                                     * function assumes, that the
+                                     * vector @p{defect} is properly
+                                     * filled with the residual in
+                                     * the outer defect correction
+                                     * scheme (usually performed by
+                                     * @ref{PreconditionMG}). After
+                                     * execution of @p{vcycle()}, the
+                                     * result is in the vector
+                                     * @p{solution}. See
+                                     * @p{copy_*_mg} in class
+                                     * @p{MGTools} if you want to use
+                                     * these vectors yourself.
+                                     *
+                                     * The actual work for this
+                                     * function is done in
+                                     * @p{level_mgstep}.
+                                     */
+    template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+    void vcycle(const MGLevelObject<MATRIX>& matrix,
+               const TRANSFER&              transfer,
+               const SMOOTHER&              pre_smooth,
+               const SMOOTHER&              post_smooth,
+               const COARSE&                cgs);
+
+                                    /**
+                                     * Print a level vector using
+                                     * @ref{DoFHandler}.
+                                     */
+/*      virtual void print_vector (const unsigned int level, */
+/*                            const VECTOR& v, */
+/*                            const char* name) const; */
+    
+
+  private:
+    
+                                    /**
+                                     * The actual v-cycle multigrid method.
+                                     * @p{level} is the level the
+                                     * function should work on. It
+                                     * will usually be called for the
+                                     * highest level from outside,
+                                     * but will then call itself
+                                     * recursively for @p{level-1},
+                                     * unless we are on @p{minlevel}
+                                     * where instead of recursing
+                                     * deeper, the coarse grid solver
+                                     * is used to solve the matrix of
+                                     * this level.
+                                     */
+    template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+    void level_mgstep (const unsigned int           level,
+                      const MGLevelObject<MATRIX>& matrix,
+                      const TRANSFER&              transfer,
+                      const SMOOTHER&              pre_smooth,
+                      const SMOOTHER&              post_smooth,
+                      const COARSE&                cgs);
+                                    /**
+                                     * Level for coarse grid solution.
+                                     */
+    unsigned int minlevel;
+
+                                    /**
+                                     * Highest level of cells.
+                                     */
+    unsigned int maxlevel;
+    
+                                    /**
+                                     * Auxiliary vector. Contains the
+                                     * defect to be corrected before
+                                     * the multigrid step.
+                                     */
+    MGLevelObject<VECTOR> defect;
+
+                                    /**
+                                     * Auxiliary vector. Contains the
+                                     * updated solution after the
+                                     * multigrid step.
+                                     */
+    MGLevelObject<VECTOR> solution;
+    
+                                    /**
+                                     * Auxiliary vector.
+                                     */
+    MGLevelObject<VECTOR> t;    
+
+
+  private:
+    
+
+                                    /**
+                                     * Exception.
+                                     */
+    DeclException2(ExcSwitchedLevels, int, int,
+                  << "minlevel and maxlevel switched, should be: "
+                  << arg1 << "<=" << arg2);
+    template<int dim, class MATRIX, class VECTOR2, class TRANSFER, class SMOOTHER, class COARSE> friend class PreconditionFAS;
+};
+
+
+/**
+ * Multi-level preconditioner.
+ * Here, we collect all information needed for multi-level preconditioning
+ * and provide the standard interface for LAC iterative methods.
+ *
+ * The template parameter class @p{MG} is required to inherit @p{MGBase}.
+ * Furthermore, it needs functions @p{void copy_to_mg(const VECTOR&)}
+ * to store @p{src} in the right hand side of the multi-level method and
+ * @p{void copy_from_mg(VECTOR&)} to store the result of the v-cycle in @p{dst}.
+ *
+ * @author Guido Kanschat, 1999
+ */
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+class PreconditionFAS : public Subscriptor
+{
+  public:
+                                    /**
+                                     * Constructor.
+                                     * Arguments are the multigrid object,
+                                     * pre-smoother, post-smoother and
+                                     * coarse grid solver.
+                                     */
+    PreconditionFAS(const MGDoFHandler<dim>&     mg_dof,
+                  FAS<VECTOR>&           mg,
+                  const MGLevelObject<MATRIX>& matrix,
+                  const TRANSFER&              transfer,
+                  const SMOOTHER&              pre,
+                  const SMOOTHER&              post,
+                  const COARSE&                coarse);
+    
+                                    /**
+                                     * Preconditioning operator.
+                                     * Calls the @p{vcycle} function
+                                     * of the @p{MG} object passed to
+                                     * the constructor.
+                                     *
+                                     * This is the operator used by
+                                     * LAC iterative solvers.
+                                     */
+    void vmult (VECTOR       &dst,
+               const VECTOR &src) const;
+    
+                                    /**
+                                     * Preconditioning operator.
+                                     * Calls the @p{vcycle} function
+                                     * of the @p{MG} object passed to
+                                     * the constructor.
+                                     */
+  void vmult_add (VECTOR       &dst,
+                 const VECTOR &src) const;
+    
+                                    /**
+                                     * Tranposed preconditioning operator.
+                                     *
+                                     * Not implemented, but the
+                                     * definition may be needed.
+                                     */
+    void Tvmult (VECTOR       &dst,
+               const VECTOR &src) const;
+    
+                                    /**
+                                     * Tranposed preconditioning operator.
+                                     *
+                                     * Not implemented, but the
+                                     * definition may be needed.
+                                     */
+    void Tvmult_add (VECTOR       &dst,
+                    const VECTOR &src) const;
+    
+  private:
+                                    /**
+                                     * Associated @p{MGDoFHandler}.
+                                     */
+    SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
+
+                                    /**
+                                     * The multigrid object.
+                                     */
+    SmartPointer<FAS<VECTOR> > multigrid;
+    
+                                    /**
+                                     * Matrices for each level. The
+                                     * matrices are prepared by the
+                                     * constructor of @p{FAS} and
+                                     * can be accessed for
+                                     * assembling.
+                                     */
+    SmartPointer<const MGLevelObject<MATRIX> > level_matrices;
+
+                                  /**
+                                   * Object for grid tranfer.
+                                   */
+    SmartPointer<const TRANSFER> transfer;
+  
+                                    /**
+                                     * The pre-smoothing object.
+                                     */
+    SmartPointer<const SMOOTHER> pre;
+
+                                    /**
+                                     * The post-smoothing object.
+                                     */
+    SmartPointer<const SMOOTHER> post;
+    
+                                    /**
+                                     * The coarse grid solver.
+                                     */
+    SmartPointer<const COARSE> coarse;
+};
+
+
+
+
+/* --------------------------- inline functions --------------------- */
+
+
+template <class VECTOR>
+template <int dim>
+FAS<VECTOR>::FAS (const MGDoFHandler<dim>& mg_dof_handler,
+                             const unsigned int min_level,
+                             const unsigned int max_level)
+               :
+               minlevel(min_level),
+                                    maxlevel(std::min(mg_dof_handler.get_tria().n_levels()-1,
+                                                      max_level)),
+                                    defect(minlevel,maxlevel),
+                                    solution(minlevel,maxlevel),
+                                    t(minlevel,maxlevel)
+{
+};
+
+
+template <class VECTOR>
+template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+void
+FAS<VECTOR>::level_mgstep(const unsigned int        level,
+                               const MGLevelObject<MATRIX>& matrix,
+                               const TRANSFER&              transfer,
+                               const SMOOTHER     &pre_smooth,
+                               const SMOOTHER     &post_smooth,
+                               const COARSE &coarse_grid_solver)
+{
+#ifdef MG_DEBUG
+  char *name = new char[100];
+  sprintf(name, "MG%d-defect",level);
+  print_vector(level, defect[level], name);
+#endif
+
+  solution[level] = 0.;
+  
+  if (level == minlevel)
+    {
+      coarse_grid_solver(level, solution[level], defect[level]);
+#ifdef MG_DEBUG
+      sprintf(name, "MG%d-solution",level);
+      print_vector(level, solution[level], name);
+#endif
+      return;
+    }
+
+                          // smoothing of the residual by modifying s
+  pre_smooth.smooth(level, solution[level], defect[level]);
+
+#ifdef MG_DEBUG
+  sprintf(name, "MG%d-pre",level);
+  print_vector(level, solution[level], name);
+#endif
+  
+                                  // t = -A*solution[level]
+  matrix[level].vmult(t[level], solution[level]);
+  t[level].scale(-1);
+  
+                                  // make t rhs of lower level
+                                  // The non-refined parts of the
+                                  // coarse-level defect already contain
+                                  // the global defect, the refined parts
+                                  // its restriction.
+  for (unsigned int l = level;l>minlevel;--l)
+    {
+      t[l-1] = 0.;
+      transfer.restrict_and_add (l, t[l-1], t[l]);
+      defect[l-1] += t[l-1];
+    }
+
+                                  // add additional DG contribution
+//  edge_vmult(level, defect[level-1], defect[level]);
+  
+                                  // do recursion
+  level_mgstep(level-1, matrix, transfer, pre_smooth, post_smooth, coarse_grid_solver);
+
+                                  // reset size of the auxiliary
+                                  // vector, since it has been
+                                  // resized in the recursive call to
+                                  // level_mgstep directly above
+  t[level] = 0.;
+
+                                  // do coarse grid correction
+
+  transfer.prolongate(level, t[level], solution[level-1]);
+
+#ifdef MG_DEBUG
+  sprintf(name, "MG%d-cgc",level);
+  print_vector(level, t, name);
+#endif
+
+  solution[level] += t[level];
+  
+                                  // post-smoothing
+  post_smooth.smooth(level, solution[level], defect[level]);
+
+#ifdef MG_DEBUG
+  sprintf(name, "MG%d-post",level);
+  print_vector(level, solution[level], name);
+
+  delete[] name;
+#endif
+}
+
+
+template <class VECTOR>
+template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+void
+FAS<VECTOR>::vcycle(const MGLevelObject<MATRIX>& matrix,
+                         const TRANSFER&              transfer,
+                         const SMOOTHER     &pre_smooth,
+                         const SMOOTHER     &post_smooth,
+                         const COARSE &coarse_grid_solver)
+{
+                                  // The defect vector has been
+                                  // initialized by copy_to_mg. Now
+                                  // adjust the other vectors.
+  solution.resize(minlevel, maxlevel);
+  t.resize(minlevel, maxlevel);
+  
+  for (unsigned int level=minlevel; level<=maxlevel;++level)
+    {
+      solution[level].reinit(defect[level]);
+      t[level].reinit(defect[level]);
+    }
+
+  level_mgstep (maxlevel, matrix, transfer, pre_smooth, post_smooth, coarse_grid_solver);
+//  abort ();
+}
+
+
+/* ------------------------------------------------------------------- */
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+PreconditionFAS<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>
+::PreconditionFAS(const MGDoFHandler<dim>&      mg_dof_handler,
+                FAS<VECTOR>& mg,
+                const MGLevelObject<MATRIX>& matrix,
+                const TRANSFER&              transfer,
+                const SMOOTHER&              pre,
+                const SMOOTHER&              post,
+                const COARSE&                coarse)
+               :
+  mg_dof_handler(&mg_dof_handler),
+  multigrid(&mg),
+  level_matrices(&matrix),
+  transfer(&transfer),
+  pre(&pre),
+  post(&post),
+  coarse(&coarse)
+{}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionFAS<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::vmult (
+  VECTOR& dst,
+  const VECTOR& src) const
+{
+  transfer->copy_to_mg(*mg_dof_handler,
+                      multigrid->defect,
+                      src);
+  
+  multigrid->vcycle(*level_matrices,
+                   *transfer,
+                   *pre, *post, *coarse);
+
+  transfer->copy_from_mg(*mg_dof_handler,
+                        dst,
+                        multigrid->solution);
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionFAS<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::vmult_add (
+  VECTOR& dst,
+  const VECTOR& src) const
+{
+  transfer->copy_to_mg(*mg_dof_handler,
+                      multigrid->defect,
+                      src);
+  
+  multigrid->vcycle(*level_matrices,
+                   *transfer,
+                   *pre, *post, *coarse);
+
+  transfer->copy_from_mg_add(*mg_dof_handler,
+                            dst,
+                            multigrid->solution);
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionFAS<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::Tvmult (
+  VECTOR&,
+  const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionFAS<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::Tvmult_add (
+  VECTOR&,
+  const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+#endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.