]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added interfaces to Trilinos block solvers and a block preconditioner (for Stokes...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 13 Oct 2008 09:20:41 +0000 (09:20 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 13 Oct 2008 09:20:41 +0000 (09:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@17188 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_precondition.h
deal.II/lac/include/lac/trilinos_precondition_block.h [new file with mode: 0644]
deal.II/lac/include/lac/trilinos_solver.h
deal.II/lac/include/lac/trilinos_solver_block.h [new file with mode: 0644]
deal.II/lac/source/trilinos_precondition.cc
deal.II/lac/source/trilinos_precondition_block.cc [new file with mode: 0755]
deal.II/lac/source/trilinos_solver.cc
deal.II/lac/source/trilinos_solver_block.cc [new file with mode: 0644]

index 54db649b4d172ae4c82d867912d4bfdada364b8c..ede32aa72376995c2b7af842e7d3c3b52233efee 100755 (executable)
@@ -53,6 +53,7 @@ namespace TrilinosWrappers
 
   class VectorBase;
   class SparseMatrix;
+  class BlockSparseMatrix;
   class SolverBase;
 
 /**
@@ -124,9 +125,9 @@ namespace TrilinosWrappers
                      << ". Check preconditioner and matrix setup.");
 
       friend class SolverBase;
-      friend class SolverSaddlePoint;
+      friend class PreconditionStokes;
 
-    protected:
+      //protected:
                                       /**
                                        * This is a pointer to the
                                        * preconditioner object that
diff --git a/deal.II/lac/include/lac/trilinos_precondition_block.h b/deal.II/lac/include/lac/trilinos_precondition_block.h
new file mode 100644 (file)
index 0000000..25aff11
--- /dev/null
@@ -0,0 +1,448 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_precondition_block_h
+#define __deal2__trilinos_precondition_block_h
+
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/trilinos_sparse_matrix.h>
+#include <lac/trilinos_block_sparse_matrix.h>
+#include <lac/trilinos_precondition.h>
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Teuchos_RefCountPtr.hpp>
+#include <Thyra_LinearOperatorDecl.hpp>
+#include <Epetra_Operator.h>
+
+// some forward declarations
+class Ifpack_Preconditioner;
+
+
+DEAL_II_NAMESPACE_OPEN
+
+/*! @addtogroup TrilinosWrappers
+ *@{
+ */
+
+namespace TrilinosWrappers
+{
+
+/**
+ * This implements an interface class of preconditioners for block
+ * matrices based on Trilinos, which are intended to be used with the
+ * solver classes SolverBlockCG and SolverBlockGMRES. This is based on
+ * the Trilinos Thyra abstract interface of linear operators. This
+ * class does not implement any method, and derived classes will need
+ * to do that job.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionBlockBase : public Subscriptor
+  {
+    public:
+                                      /**
+                                       * Empty constructor.
+                                       */
+      PreconditionBlockBase();
+
+                                      /**
+                                       * Destructor.
+                                       */
+      ~PreconditionBlockBase();
+
+    protected:
+                                      /**
+                                       * Preconditioner object.
+                                       */
+      Teuchos::RCP<Thyra::ConstLinearOperator<double> > preconditioner;
+
+    friend class SolverBlockBase;
+};
+
+                                      /**
+                                       * Internal function that sets
+                                       * up an inverse matrix.
+                                       */
+  inline
+  Thyra::ConstLinearOperator<double>
+  inverse_matrix (const SparseMatrix    &M,
+                 const Epetra_Operator *P,
+                 const bool             is_symmetric,
+                 const unsigned int     n_iterations,
+                 const double           solve_tolerance,
+                 const bool             output_details);
+
+
+                                       // Forward declarations.
+
+/**  
+ * This class implements a black box preconditioner for saddle points
+ * systems arising from the Stokes or Navier&ndash;Stokes equations as
+ * specified by the papers <em>D. Silvester, A. Wathen, Fast iterative
+ * solution of stabilised Stokes systems part II. Using general block
+ * preconditioners, SIAM J. Numer. Anal. 31:1352&ndash1367 (1994)</em>
+ * and <em> D. Kay, D. Loghin, A. Wathen, A preconditioner for the
+ * steady-state Navier&ndash;Stokes equations, SIAM
+ * J. Sci. Comput. 24(1):237&ndash;256 (2002)</em>, respectively.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Preconditioners
+ * @author Martin Kronbichler, 2008
+ */
+  class PreconditionStokes : public PreconditionBlockBase
+  {
+    public:
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional data to the
+                                        * solver.
+                                        */
+      struct AdditionalData
+      {
+                                       /**
+                                       * Constructor. This sets the
+                                       * additional data to its
+                                       * default values. For lazy
+                                       * initializations, the user
+                                       * will just set the parameters
+                                       * in here and let the
+                                       * <tt>initialize</tt> function
+                                       * set up good preconditioners,
+                                       * where some good values have
+                                       * to be provided here. The
+                                       * default ones mean that we
+                                       * assume the velocity-velocity
+                                       * matrix to be symmetric (only
+                                       * true for Stokes problems!),
+                                       * we use a relative inner
+                                       * tolerance of 1e-9 for the
+                                       * inversion of matrices, and
+                                       * take in IC preconditioner
+                                       * for the inversion of the
+                                       * pressure mass
+                                       * matrix. Moreover, there is
+                                       * no inner CG/GMRES solve
+                                       * performed on the
+                                       * velocity-velocity block by
+                                       * default. Otherwise, we do a
+                                       * solve with the specified
+                                       * solver and the number of
+                                       * iterations and
+                                       * tolerance. Note that too
+                                       * many iterations on Av can
+                                       * slow down the overall
+                                       * procedure considerably.
+                                       */
+       AdditionalData (const bool          right_preconditioning = false,
+                       const bool          Av_is_symmetric = true,
+                       const std::vector<std::vector<bool> > &Av_null_space = std::vector<std::vector<bool> > (1),
+                       const double        inner_solve_tolerance = 1e-9,
+                       const bool          use_ssor_on_mass_matrix = false,
+                       const bool          velocity_uses_higher_order_elements = false,
+                       const bool          pressure_uses_higher_order_elements = false,
+                       const bool          Av_do_solve   = false,
+                       const unsigned int  Av_n_iters    = 1,
+                       const double        Av_tolerance  = 1e-3,
+                       const bool          ouput_details = false);
+
+                                      /**
+                                       * This flag specifies whether
+                                       * the preconditioner should be
+                                       * build as a left
+                                       * preconditioner or right
+                                       * preconditioner. Note that
+                                       * this setting must be the
+                                       * same as the one used in the
+                                       * solver call.
+                                       */
+       bool right_preconditioning;
+
+                                      /**
+                                       * This flag determines whether
+                                       * the velocity-velocity matrix
+                                       * is symmetric (and positive
+                                       * definite) or not. This
+                                       * choice will influence the
+                                       * AMG preconditioner that is
+                                       * built for that system as
+                                       * well as the inner solve on
+                                       * Av (if that flag is
+                                       * enabled).
+                                       */
+       bool Av_is_symmetric;
+
+                                      /**
+                                       * This determines the near
+                                       * null space for the operator
+                                       * Av, the vector-valued
+                                       * velocity-velocity space. In
+                                       * order to get good
+                                       * performance, this vector has
+                                       * to be specified using the
+                                       * DoFTools::extract_constant_modes()
+                                       * function.
+                                       */
+       std::vector<std::vector<bool> > Av_null_space;
+
+                                      /**
+                                       * Tolerance level that the
+                                       * inner solvers on the
+                                       * pressure mass matrix (and
+                                       * the pressure Laplace matrix
+                                       * in case of a Navier-Stokes
+                                       * problem) should be solve
+                                       * to. It is essential that
+                                       * this tolerance is in the
+                                       * same order of magnitude than
+                                       * the one chosen for the outer
+                                       * GMRES solver (or little
+                                       * less). Otherwise unexpected
+                                       * variations in the number of
+                                       * iterations can occur from
+                                       * solve to solve. 
+                                       *
+                                       * Default value: <tt>1e-9</tt>
+                                       * residual relative to initial
+                                       * residual.
+                                       */
+       double inner_solve_tolerance;
+
+                                      /**
+                                       * Determines whether the inner
+                                       * solve for the pressure mass
+                                       * matrix should use an IC
+                                       * preconditioner (incomplete
+                                       * Cholesky factorization) or
+                                       * an SSOR preconditioner. The
+                                       * former tends to be faster in
+                                       * most situations, whereas the
+                                       * latter uses almost no
+                                       * additional memory besides
+                                       * the matrix storage.
+                                       */
+       bool use_ssor_on_mass_matrix;
+
+                                      /**
+                                       * Determines whether the
+                                       * underlying velocity
+                                       * discretization uses linear
+                                       * elements or higher order
+                                       * elements, which has
+                                       * consequences for the set up
+                                       * of the internal setup of the
+                                       * ML AMG preconditioner.
+                                       */
+       bool velocity_uses_higher_order_elements;
+
+                                      /**
+                                       * Determines whether the
+                                       * underlying pressure
+                                       * discretization uses linear
+                                       * elements or higher order
+                                       * elements, which has
+                                       * consequences for the set up
+                                       * of the internal setup of the
+                                       * ML AMG preconditioner.
+                                       */
+       bool pressure_uses_higher_order_elements;
+
+                                       /**
+                                       * Flag to determine whether we
+                                       * should do an inner solve on
+                                       * the velocity-velocity
+                                       * matrix. By default, this
+                                       * option is not on, since it
+                                       * is not necessary for a good
+                                       * performance of the outer
+                                       * solver. However, it can be
+                                       * beneficial to perform some
+                                       * iterations on this block
+                                       * only and not on the whole
+                                       * block, so that the number of
+                                       * outer iterations is
+                                       * reduced. See the discussion
+                                       * in the @step_22 step-22
+                                       * tutorial program.
+                                       */
+       unsigned int Av_do_solve;
+         
+                                      /**
+                                       * The number of iterations in
+                                       * the inner solve.
+                                       */
+       unsigned int Av_n_iters;
+         
+                                      /**
+                                       * The residual tolerance that
+                                       * is going to be used for the
+                                       * inner solve on the
+                                       * velocity-velocity matrix.
+                                       */
+       double Av_tolerance;
+
+                                      /**
+                                       * This defines whether
+                                       * internal details about the
+                                       * solve process should be
+                                       * written to screen. This can
+                                       * be a lot of information.
+                                       */
+       bool output_details;
+      };
+
+                                      /**
+                                       * Lazy setup of the block
+                                       * preconditioner for the
+                                       * (Navier-) Stokes
+                                       * system. This function takes
+                                       * the matrices given here and
+                                       * firs calculates good
+                                       * preconditioners, i.e.,
+                                       * algebraic multigrid
+                                       * preconditioners for the
+                                       * velocity-velocity matrix
+                                       * <tt>Av</tt>, that can be
+                                       * specified to be different
+                                       * from the (0,0) block in the
+                                       * system matrix, IC/SSOR on
+                                       * the pressure mass matrix
+                                       * <tt>Mp_matrix</tt>, and AMG
+                                       * on the pressure Laplace
+                                       * matrix
+                                       * <tt>Lp_matrix</tt>. Next,
+                                       * these preconditioners are
+                                       * used to generate a block
+                                       * operator.
+                                       */
+      void initialize (const BlockSparseMatrix  &system_matrix,
+                      const SparseMatrix       &Mp_matrix,
+                      const AdditionalData     &additional_data,
+                      const SparseMatrix       &Av_matrix = SparseMatrix(),
+                      const SparseMatrix       &Lp_matrix = SparseMatrix(),
+                      const SparseMatrix       &Fp_matrix = SparseMatrix());
+
+                                      /**
+                                       * Set up the block
+                                       * preconditioner for the
+                                       * (Navier-) Stokes system. In
+                                       * contrast to the other
+                                       * <tt>initialize</tt>
+                                       * function, this one expects
+                                       * the user to specify
+                                       * preconditioner objects for
+                                       * the individual matrices,
+                                       * that will then be used for
+                                       * the inner inversion of the
+                                       * matrices when building the
+                                       * Schur complement
+                                       * preconditioner. If no
+                                       * preconditioner is specified,
+                                       * the inner solvers will use
+                                       * an identity preconditioner
+                                       * object.
+                                       */
+      void initialize (const BlockSparseMatrix  &system_matrix,
+                      const SparseMatrix       &Mp_matrix,
+                      const AdditionalData     &additional_data,
+                      const PreconditionBase   &Mp_preconditioner = PreconditionBase(),
+                      const PreconditionBase   &Av_preconditioner = PreconditionBase(),
+                      const SparseMatrix       &Lp_matrix = SparseMatrix(),
+                      const PreconditionBase   &Lp_preconditioner = PreconditionBase(),
+                      const SparseMatrix       &Fp_matrix = SparseMatrix());
+
+                                      /**
+                                       * This function can be used
+                                       * for a faster recalculation
+                                       * of the preconditioner
+                                       * construction when the system
+                                       * matrix entries underlying
+                                       * the preconditioner have
+                                       * changed but not the sparsity
+                                       * pattern (this means that the
+                                       * grid has remained the same
+                                       * and the matrix structure has
+                                       * not been
+                                       * reinitialized). Moreover, it
+                                       * is assumed that the PDE
+                                       * parameters have not changed
+                                       * drastically. This function
+                                       * is only needed for the lazy
+                                       * variant of the
+                                       * preconditioner. In the other
+                                       * case, the preconditioner can
+                                       * be modified outside this
+                                       * function, which will be
+                                       * recongnized in this
+                                       * preconditioner as well when
+                                       * doing a solve, without the
+                                       * need to restructure anything
+                                       * here.
+                                       */
+      void reinit_lazy ();
+
+                                       /**
+                                       * Exception.
+                                       */
+      DeclException1 (ExcNonMatchingMaps,
+                     std::string,
+                     << "The sparse matrix the preconditioner is based on "
+                     << "uses a map that is not compatible to the one in vector"
+                     << arg1
+                     << ". Check preconditioner and matrix setup.");
+
+    private:
+
+                                      /**
+                                       * Pointer to preconditioner
+                                       * for the mass matrix.
+                                       */
+      Teuchos::RCP<PreconditionSSOR> Mp_precondition_ssor;
+
+                                      /**
+                                       * Pointer to preconditioner
+                                       * for the mass matrix.
+                                       */
+      Teuchos::RCP<PreconditionIC>   Mp_precondition_ic;
+
+                                      /**
+                                       * Pointer to preconditioner
+                                       * for the velocity-velocity
+                                       * matrix.
+                                       */
+      Teuchos::RCP<PreconditionAMG> Av_precondition;
+
+                                      /**
+                                       * Pointer to preconditioner
+                                       * for the pressure Laplace
+                                       * matrix.
+                                       */
+      Teuchos::RCP<PreconditionAMG> Lp_precondition;
+  };
+}
+
+/*@}*/
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*--------------------------   trilinos_precondition_block.h   ---------------------------*/
+
+#endif
+/*--------------------------   trilinos_precondition_block.h   ---------------------------*/
index fbddcef679cfbcd123919478b5fc01f204cc44f6..f7a01c17f16b398d0f0b24ab5cecd2636aa01821 100644 (file)
@@ -47,6 +47,12 @@ namespace TrilinosWrappers
  * "http://trilinos.sandia.gov/packages/aztecoo/AztecOOUserGuide.pdf">AztecOO
  * user guide</a>.
  *
+ * This solver class can also be used as a standalone class, where the
+ * respective Krylov method is set via the flag
+ * <tt>solver_name</tt>. This can be done at runtime (e.g., when
+ * parsing the solver from a ParameterList) and is similar to the
+ * deal.II class SolverSelector.
+ *
  * @ingroup TrilinosWrappers
  * @author Martin Kronbichler, 2008
  */
@@ -54,6 +60,22 @@ namespace TrilinosWrappers
   {
     public:
 
+                                      /**
+                                       * Enumeration object that is
+                                       * set in the constructor of
+                                       * the derived classes and
+                                       * tells Trilinos which solver
+                                       * to use. This option can also
+                                       * be set in the user program,
+                                       * so one might use this base
+                                       * class instead of one of the
+                                       * specialized derived classes
+                                       * when the solver should be
+                                       * set at runtime. Currently
+                                       * enabled options are:
+                                       */
+      enum SolverName {cg, cgs, gmres, bicgstab, tfqmr} solver_name;
+
                                        /**
                                         * Standardized data struct to
                                         * pipe additional data to the
@@ -77,11 +99,20 @@ namespace TrilinosWrappers
                                        /**
                                         * Constructor. Takes the
                                         * solver control object and
-                                        * the MPI communicator over
-                                        * which parallel computations
-                                        * are to happen.
+                                        * creates the solver.
                                         */
       SolverBase (SolverControl  &cn);
+
+                                       /**
+                                        * Second constructor. This
+                                        * constructor takes an enum
+                                        * object that specifies the
+                                        * solver name and sets the
+                                        * appropriate Krylov
+                                        * method.
+                                        */
+      SolverBase (const enum SolverName  solver_name,
+                 SolverControl         &cn);
       
                                        /**
                                         * Destructor.
@@ -135,18 +166,6 @@ namespace TrilinosWrappers
                                         */
       SolverControl &solver_control;
 
-                                      /**
-                                       * String object that is set in
-                                       * the constructor of the
-                                       * derived classes and tells
-                                       * Trilinos which solver to
-                                       * use. Currently enabled
-                                       * options are <tt>"cg", "cgs",
-                                       * "gmres", "bicgstab"</tt>,
-                                       * and <tt>"tfqmr"</tt>.
-                                       */
-      enum SolverName {cg, cgs, gmres, bicgstab, tfqmr} solver_name;
-
     private:
 
                                        /**
diff --git a/deal.II/lac/include/lac/trilinos_solver_block.h b/deal.II/lac/include/lac/trilinos_solver_block.h
new file mode 100644 (file)
index 0000000..5b13a87
--- /dev/null
@@ -0,0 +1,366 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_block_solver_h
+#define __deal2__trilinos_block_solver_h
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/solver_control.h>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace TrilinosWrappers
+{
+                                   // forward declarations
+  class BlockSparseMatrix;
+  class BlockVector;
+  namespace MPI
+  {
+    class BlockVector;
+  }
+  class SparseMatrix;
+  class PreconditionBase;
+  class PreconditionBlockBase;
+  
+  
+/**
+ * Base class for solver classes using the Trilinos solvers on block
+ * matrices based on the Trilinos abstract solver interface Thyra.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+  class SolverBlockBase
+  {
+    public:
+
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional data to the
+                                        * solver.
+                                        */
+      struct AdditionalData
+      {
+                                      /**
+                                       * Restart parameter in case
+                                       * the derived class is
+                                       * GMRES.
+                                       * 
+                                       * TODO: Find a better way for
+                                       * doing this.
+                                       */
+       AdditionalData (const unsigned int gmres_restart_parameter = 30,
+                       const bool         right_preconditioning   = false,
+                       const bool         output_details          = false);
+
+                                      /**
+                                       * The number of vectors that
+                                       * should be used in the GMRES
+                                       * algorithm before the solver
+                                       * is restarted.
+                                       */
+       unsigned int gmres_restart_parameter;
+
+                                      /**
+                                       * Flag that turns right
+                                       * preconditioning on.
+                                       */
+       bool right_preconditioning;
+
+                                      /**
+                                       * A flag to determine whether
+                                       * solver details should be
+                                       * written to screen.
+                                       */
+       bool output_details;
+      };
+
+                                       /**
+                                        * Constructor. Takes the
+                                        * solver control object and
+                                        * the MPI communicator over
+                                        * which parallel computations
+                                        * are to happen.
+                                        */
+      SolverBlockBase (SolverControl  &cn);
+      
+                                       /**
+                                        * Destructor.
+                                        */
+      virtual ~SolverBlockBase ();
+
+                                       /**
+                                        * Solve the linear system
+                                        * <tt>Ax=b</tt>. Depending on
+                                        * the information provided by
+                                        * derived classes and the
+                                        * object passed as a
+                                        * preconditioner, one of the
+                                        * linear solvers and
+                                        * preconditioners of Trilinos
+                                        * is chosen. This interface is
+                                        * designed for Trilinos
+                                        * running in parallel.
+                                        */
+      void
+      solve (const BlockSparseMatrix     &A,
+             MPI::BlockVector            &x,
+             const MPI::BlockVector      &b,
+             const PreconditionBlockBase &preconditioner);
+
+                                       /**
+                                        * Solve the linear system
+                                        * <tt>Ax=b</tt>. Depending on
+                                        * the information provided by
+                                        * derived classes and the
+                                        * object passed as a
+                                        * preconditioner, one of the
+                                        * linear solvers and
+                                        * preconditioners of Trilinos
+                                        * is chosen. This interface is
+                                        * designed for Trilinos
+                                        * running in serial.
+                                        */
+      void
+      solve (const BlockSparseMatrix     &A,
+             BlockVector                 &x,
+             const BlockVector           &b,
+             const PreconditionBlockBase &preconditioner);
+
+                                       /**
+                                        * Access to object that controls
+                                        * convergence.
+                                        */
+      SolverControl & control() const;
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException1 (ExcTrilinosError,
+                      int,
+                      << "An error with error number " << arg1
+                      << " occured while calling a Trilinos function");
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException2 (ExcOverlappingMaps,
+                      std::string, std::string,
+                      << "The Epetra map in the " << arg1 << " " << arg2
+                      << " is overlapping between the processors, which"
+                     << " is forbidden when doing a solve.");
+
+                                       /**
+                                       * Exception.
+                                       */
+      DeclException1 (ExcNonMatchingMaps,
+                     std::string,
+                     << "The Epetra map in the vector "
+                     << arg1
+                     << " does not match the one used in the matrix. "
+                     << "Check vector and matrix setup.");
+
+    protected:
+
+                                       /**
+                                        * Reference to the object that
+                                        * controls convergence of the
+                                        * iterative solver. In fact,
+                                        * for these Trilinos wrappers,
+                                        * Trilinos does so itself, but
+                                        * we copy the data from this
+                                        * object before starting the
+                                        * solution process, and copy
+                                        * the data back into it
+                                        * afterwards.
+                                        */
+      SolverControl &solver_control;
+
+                                      /**
+                                       * String object that is set in
+                                       * the constructor of the
+                                       * derived classes and tells
+                                       * Trilinos which solver to
+                                       * use. Currently enabled
+                                       * options are <tt>"cg",
+                                       * "gmres"</tt>.
+                                       */
+      enum SolverBlockName {cg, gmres} solver_name;
+
+    protected:
+
+                                       /**
+                                        * Store a copy of the flags for this
+                                        * particular solver.
+                                        */
+      AdditionalData additional_data;
+
+  };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos CG
+ * solver for block matrices and vectors.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+  class SolverBlockCG : public SolverBlockBase
+  {
+    public:
+
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional data to the
+                                        * solver.
+                                        */
+      struct AdditionalData
+      {
+                                      /**
+                                       * Restart parameter in case
+                                       * the derived class is
+                                       * GMRES.
+                                       * 
+                                       * TODO: Find a better way for
+                                       * doing this.
+                                       */
+       AdditionalData (const bool         right_preconditioning   = false,
+                       const bool         output_details          = false);
+
+                                      /**
+                                       * Flag that turns right
+                                       * preconditioning on.
+                                       */
+       bool right_preconditioning;
+
+                                      /**
+                                       * A flag to determine whether
+                                       * solver details should be
+                                       * written to screen.
+                                       */
+       bool output_details;
+      };
+
+                                       /**
+                                        * Constructor. In contrast to
+                                        * deal.II's own solvers, there is no
+                                        * need to give a vector memory
+                                        * object.
+                                       *
+                                        * The last argument takes a structure
+                                        * with additional, solver dependent
+                                        * flags for tuning.
+                                        */
+      SolverBlockCG (SolverControl        &cn,
+                    const AdditionalData &data = AdditionalData());
+
+    protected:
+
+                                       /**
+                                        * Store a copy of the flags for this
+                                        * particular solver.
+                                        */
+      const AdditionalData data;
+  };
+
+
+
+/**
+ * An implementation of the solver interface using the Trilinos GMRES
+ * solver for block matrices and vectors.
+ *
+ * @author Martin Kronbichler, 2008
+ */
+  class SolverBlockGMRES : public SolverBlockBase
+  {
+    public:
+
+                                       /**
+                                        * Standardized data struct to
+                                        * pipe additional data to the
+                                        * solver.
+                                        */
+      struct AdditionalData
+      {
+                                      /**
+                                       * Restart parameter in case
+                                       * the derived class is
+                                       * GMRES.
+                                       * 
+                                       * TODO: Find a better way for
+                                       * doing this.
+                                       */
+       AdditionalData (const unsigned int gmres_restart_parameter = 30,
+                       const bool         right_preconditioning   = false,
+                       const bool         output_details          = false);
+
+                                      /**
+                                       * The number of vectors that
+                                       * should be used in the GMRES
+                                       * algorithm before the solver
+                                       * is restarted.
+                                       */
+       unsigned int gmres_restart_parameter;
+
+                                      /**
+                                       * Flag that turns right
+                                       * preconditioning on.
+                                       */
+       bool right_preconditioning;
+
+                                      /**
+                                       * A flag to determine whether
+                                       * solver details should be
+                                       * written to screen.
+                                       */
+       bool output_details;
+      };
+
+                                       /**
+                                        * Constructor. In contrast to
+                                        * deal.II's own solvers, there is no
+                                        * need to give a vector memory
+                                        * object.
+                                        *
+                                        * The last argument takes a structure
+                                        * with additional, solver dependent
+                                        * flags for tuning.
+                                        */
+      SolverBlockGMRES (SolverControl        &cn,
+                       const AdditionalData &data = AdditionalData());
+
+    protected:
+
+                                       /**
+                                        * Store a copy of the flags for this
+                                        * particular solver.
+                                        */
+      const AdditionalData data;
+  };
+
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*----------------------   trilinos_solver_block.h     -----------------------*/
+
+#endif
+/*----------------------   trilinos_solver_block.h     -----------------------*/
index b09cae729504a988fbf2761cfb49de6582fd478f..acf6fcfb677c399ce69abb0d138e5312a0cbbc73 100755 (executable)
@@ -224,8 +224,8 @@ namespace TrilinosWrappers
     ifpack.release();
 
     ifpack = Teuchos::rcp (Ifpack().Create ("point relaxation",
-                                                   &*matrix.matrix, 
-                                                   additional_data.overlap));
+                                           &*matrix.matrix, 
+                                           additional_data.overlap));
 
     Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
                                       "preconditioner"));
@@ -281,7 +281,7 @@ namespace TrilinosWrappers
                                            additional_data.overlap));
 
     Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
-                                              "preconditioner"));
+                                      "preconditioner"));
 
     int ierr;
 
@@ -332,7 +332,7 @@ namespace TrilinosWrappers
                                            additional_data.overlap));
 
     Assert (&*ifpack != 0, ExcMessage ("Trilinos could not create this "
-                                              "preconditioner"));
+                                      "preconditioner"));
 
     int ierr;
 
diff --git a/deal.II/lac/source/trilinos_precondition_block.cc b/deal.II/lac/source/trilinos_precondition_block.cc
new file mode 100755 (executable)
index 0000000..433a891
--- /dev/null
@@ -0,0 +1,423 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/trilinos_precondition_block.h>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Thyra_EpetraLinearOp.hpp>
+#include <Thyra_EpetraThyraWrappers.hpp>
+
+#include <Thyra_VectorSpaceImpl.hpp>
+#include <Thyra_LinearOperatorDecl.hpp>
+#include <Thyra_DefaultBlockedLinearOpDecl.hpp>
+#include <Thyra_DefaultBlockedLinearOp.hpp>
+#include <Thyra_LinearOpWithSolveFactoryHelpers.hpp>
+#include <Thyra_InverseLinearOperator.hpp>
+#include <Thyra_DefaultInverseLinearOp.hpp>
+#include <Thyra_AztecOOLinearOpWithSolveFactory.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+
+  PreconditionBlockBase::PreconditionBlockBase()
+  {}
+
+
+
+  PreconditionBlockBase::~PreconditionBlockBase()
+  {}
+
+
+
+                                      /**
+                                       * Internal function that sets
+                                       * up an inverse matrix.
+                                       */
+  inline
+  Thyra::ConstLinearOperator<double>
+  inverse_matrix (const SparseMatrix    &input_M,
+                 const Epetra_Operator *input_P,
+                 const bool             is_symmetric,
+                 const unsigned int     n_iterations,
+                 const double           solve_tolerance,
+                 const bool             output_details)
+  {
+
+                                       // Create discrete
+                                       // Epetra operator
+                                       // (Thyra wrapper around the
+                                       // matrix).
+    Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpM = 
+      Teuchos::rcp(new Thyra::EpetraLinearOp(Teuchos::rcp(&*input_M.matrix, 
+                                                         false)));
+
+    Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > M 
+      = Teuchos::rcp(new Thyra::DefaultLinearOpSource<double>(tmpM), true);
+
+    //Thyra::ConstLinearOperator<double> opM = tmpM;
+
+                                       // Create a Thyra version of
+                                       // the preconditioner.
+    Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpP = 
+      Teuchos::rcp(new Thyra::EpetraLinearOp(Teuchos::rcp(const_cast<Epetra_Operator*>(input_P),
+                                                         false), 
+                                            Thyra::NOTRANS,
+                                            Thyra::EPETRA_OP_APPLY_APPLY_INVERSE,
+                                            Thyra::EPETRA_OP_ADJOINT_SUPPORTED));
+       
+    //Thyra::ConstLinearOperator<double> opPM = tmpP;
+
+    Teuchos::RCP<Thyra::PreconditionerBase<double> > PM = 
+      Teuchos::rcp (new Thyra::DefaultPreconditioner<double>(), true);
+
+    {
+      Thyra::DefaultPreconditioner<double>
+       *defaultPrec = &Teuchos::dyn_cast<Thyra::DefaultPreconditioner<double> >
+                      (*PM);
+
+      (*defaultPrec).initializeRight(tmpP);
+    }
+
+                                       // Set up the solver.
+    Teuchos::RCP<Teuchos::ParameterList> aztecParams 
+      = Teuchos::rcp(new Teuchos::ParameterList("aztecOOFSolverFactory"), true);
+
+    aztecParams->sublist("Forward Solve").set("Max Iterations", 
+                                             (int)n_iterations);
+    aztecParams->sublist("Forward Solve").set("Tolerance", 
+                                             solve_tolerance);
+
+                                       // aztecOO solver settings
+    if (is_symmetric == true)
+      {
+       aztecParams->sublist("Forward Solve")
+         .sublist("AztecOO Settings").set("Aztec Solver", "CG");
+      }
+    else
+      {
+       aztecParams->sublist("Forward Solve")
+         .sublist("AztecOO Settings").set("Aztec Solver", "GMRES");
+       aztecParams->sublist("Forward Solve")
+         .sublist("AztecOO Settings").set("Size of Krylov Subspace", 
+                                          (int)n_iterations);
+      }
+
+    aztecParams->sublist("Forward Solve")
+      .sublist("AztecOO Settings").set("Aztec Preconditioner", "none");
+    aztecParams->sublist("Forward Solve")
+      .sublist("AztecOO Settings").set("Convergence Test", "r0");
+
+                                       // This sets if solver details
+                                       // (the residual in each
+                                       // iterations) should be
+                                       // written to screen.
+    if (output_details == true)
+      {
+       aztecParams->sublist("Forward Solve")
+         .sublist("AztecOO Settings").set("Output Frequency", 1);
+      }
+    else
+      {
+       aztecParams->sublist("VerboseObject").set("Verbosity Level", "none");
+      }
+
+    Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
+      aztecLowsFactory = Teuchos::rcp(new Thyra::AztecOOLinearOpWithSolveFactory(), true);
+
+    aztecLowsFactory->setParameterList(aztecParams);
+
+    Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> > rcpAztec 
+      = aztecLowsFactory->createOp();
+
+                                       // Does not work like
+                                       // this. The preconditioner
+                                       // will only appear in the
+                                       // operator, but won't be set
+                                       // correctly to the right hand
+                                       // side vector, which makes
+                                       // the results wrong. So try
+                                       // something else...
+    //Meros::LinearSolveStrategy<double> azM 
+    //  = new Meros::AztecSolveStrategy(*(aztecParams.get()));
+    //Thyra::LinearOperator<double> Minverse 
+    //  = Thyra::inverse(*aztecLowsFactory, opM);
+    //return new Meros::InverseOperator<double>(opM * opPM, azM);
+
+                                       // Attach the preconditioner
+                                       // to the solver.
+    aztecLowsFactory->initializePreconditionedOp(M, PM, &*rcpAztec, 
+                                                Thyra::SUPPORT_SOLVE_FORWARD_ONLY);
+                                       // Create an inverse matrix
+                                       // object and return.
+    Teuchos::RCP<const Thyra::LinearOpBase<double> > Minverse = 
+      Teuchos::rcp(new Thyra::DefaultInverseLinearOp<double>(rcpAztec));
+
+    return Minverse;
+  }
+
+
+
+  PreconditionStokes::AdditionalData::
+  AdditionalData (const bool          right_preconditioning,
+                 const bool          Av_is_symmetric,
+                 const std::vector<std::vector<bool> > &Av_null_space,
+                 const double        inner_solve_tolerance,
+                 const bool          use_ssor_on_mass_matrix,
+                 const bool          velocity_uses_higher_order_elements,
+                 const bool          pressure_uses_higher_order_elements,
+                 const bool          Av_do_solve,
+                 const unsigned int  Av_n_iters,
+                 const double        Av_tolerance,
+                 const bool          output_details)
+                  :
+                  right_preconditioning (right_preconditioning),
+                  Av_is_symmetric (Av_is_symmetric),
+                 Av_null_space (Av_null_space),
+                 inner_solve_tolerance (inner_solve_tolerance),
+                 use_ssor_on_mass_matrix (use_ssor_on_mass_matrix),
+                 velocity_uses_higher_order_elements (velocity_uses_higher_order_elements),
+                 pressure_uses_higher_order_elements (pressure_uses_higher_order_elements),
+                 Av_do_solve (Av_do_solve),
+                 Av_n_iters (Av_n_iters),
+                 Av_tolerance (Av_tolerance),
+                 output_details (output_details)
+  {}
+
+
+
+  void
+  PreconditionStokes::initialize(const BlockSparseMatrix  &system_matrix,
+                                const SparseMatrix       &Mp_matrix,
+                                const AdditionalData     &additional_data,
+                                const SparseMatrix       &Av_matrix,
+                                const SparseMatrix       &Lp_matrix,
+                                const SparseMatrix       &Fp_matrix)
+  {
+    AssertThrow (system_matrix.n_block_rows() == 2,
+                ExcDimensionMismatch(system_matrix.n_block_rows(),
+                                     2));
+    AssertThrow (system_matrix.n_block_cols() == 2,
+                ExcDimensionMismatch(system_matrix.n_block_cols(),
+                                     2));
+
+                                       // Create Mp preconditioner.
+    Teuchos::RCP<PreconditionBase> Mp_precondition;
+    {
+      if (additional_data.use_ssor_on_mass_matrix == true)
+       {
+         Mp_precondition_ssor = Teuchos::rcp (new PreconditionSSOR(), true);
+         Mp_precondition_ssor->initialize (Mp_matrix,
+                                        PreconditionSSOR::AdditionalData (1.2,0,0));
+         Mp_precondition = Teuchos::rcp(&*Mp_precondition_ssor, false);
+       }
+      else
+       {
+         Mp_precondition_ic = Teuchos::rcp (new PreconditionIC(), true);
+         Mp_precondition_ic->initialize (Mp_matrix);
+         Mp_precondition = Teuchos::rcp(&*Mp_precondition_ic, false);
+       }
+    }
+
+                                       // Create Av AMG
+                                       // preconditioner.
+    {
+      Av_precondition = Teuchos::rcp (new PreconditionAMG(), true);
+      PreconditionAMG::AdditionalData av_amg_data;
+      av_amg_data.elliptic = additional_data.Av_is_symmetric;
+      av_amg_data.higher_order_elements = additional_data.velocity_uses_higher_order_elements;
+      av_amg_data.null_space = additional_data.Av_null_space;
+
+      if (Av_matrix.m() == system_matrix.block(0,0).m())
+       Av_precondition->initialize (Av_matrix, av_amg_data);
+      else
+       Av_precondition->initialize (system_matrix.block(0,0), av_amg_data);
+    }
+
+                                       // Create Lp AMG
+                                       // preconditioner.
+    if (Lp_matrix.m() == system_matrix.block(1,1).m())
+      {
+       Lp_precondition = Teuchos::rcp (new PreconditionAMG(), true);
+       PreconditionAMG::AdditionalData lp_amg_data;
+       lp_amg_data.elliptic = true;
+       lp_amg_data.higher_order_elements = additional_data.pressure_uses_higher_order_elements;
+
+       Lp_precondition->initialize (Lp_matrix, lp_amg_data);
+      }
+
+    initialize (system_matrix, Mp_matrix, additional_data, *Mp_precondition,
+               *Av_precondition, Lp_matrix, *Lp_precondition, Fp_matrix);
+
+  }
+
+  void 
+  PreconditionStokes::initialize (const BlockSparseMatrix  &system_matrix,
+                                 const SparseMatrix       &Mp_matrix,
+                                 const AdditionalData     &additional_data,
+                                 const PreconditionBase   &Mp_preconditioner,
+                                 const PreconditionBase   &Av_preconditioner,
+                                 const SparseMatrix       &Lp_matrix,
+                                 const PreconditionBase   &Lp_preconditioner,
+                                 const SparseMatrix       &Fp_matrix)
+  {
+    AssertThrow (system_matrix.n_block_rows() == 2,
+                ExcDimensionMismatch(system_matrix.n_block_rows(),
+                                     2));
+    AssertThrow (system_matrix.n_block_cols() == 2,
+                ExcDimensionMismatch(system_matrix.n_block_cols(),
+                                     2));
+
+                                       // Create Av inverse matrix
+    Thyra::ConstLinearOperator<double> Av_inverse;
+    if (additional_data.Av_do_solve)
+      Av_inverse = inverse_matrix(system_matrix.block(0,0),
+                                 &*Av_preconditioner.preconditioner,
+                                 additional_data.Av_is_symmetric,
+                                 additional_data.Av_n_iters,
+                                 additional_data.Av_tolerance,
+                                 additional_data.output_details);
+    else
+      {
+       Teuchos::RCP<const Thyra::LinearOpBase<double> >
+         tmpAvp = Thyra::epetraLinearOp(Teuchos::rcp(&*Av_preconditioner.preconditioner,
+                                                     false), 
+                                        Thyra::NOTRANS,
+                                        Thyra::EPETRA_OP_APPLY_APPLY_INVERSE,
+                                        Thyra::EPETRA_OP_ADJOINT_SUPPORTED);
+       Av_inverse = tmpAvp;
+      }
+
+
+                                       // Create Mp inverse matrix.
+    Thyra::ConstLinearOperator<double> Mp_inverse 
+      = inverse_matrix(Mp_matrix, &*Mp_preconditioner.preconditioner, true, 
+                      Mp_matrix.m(),
+                      additional_data.inner_solve_tolerance,
+                      additional_data.output_details);
+
+                                       // Create Lp preconditioner
+                                       // and inverse matrix in case
+                                       // there is a matrix given.
+    Thyra::ConstLinearOperator<double> Lp_inverse;
+
+    if (Lp_matrix.m() == system_matrix.block(1,1).m())
+      {
+                                       // Create Lp inverse matrix.
+       Lp_inverse = inverse_matrix(Lp_matrix, &*Lp_preconditioner.preconditioner,
+                                   true, Lp_matrix.m(), 
+                                   additional_data.inner_solve_tolerance,
+                                   additional_data.output_details);
+      }
+    else
+      Lp_inverse = Thyra::identity<double> (Mp_inverse.range());
+
+                                       // Create Fp Thyra operator.
+    Thyra::ConstLinearOperator<double> Fp_op;
+    if (Fp_matrix.m() != system_matrix.block(1,1).m())
+      {
+       Fp_op = Thyra::identity<double> (Mp_inverse.range());
+      }
+    else
+      {
+       Teuchos::RCP<const Thyra::LinearOpBase<double> >
+         tmpFp = Thyra::epetraLinearOp(Teuchos::rcp(&*Fp_matrix.matrix,false));
+       Fp_op = tmpFp;
+      }
+
+                                       // Build the PCD block
+                                       // preconditioner factory.
+                                       // Build identity matrices on
+                                       // the velocity and pressure
+                                       // spaces
+    Thyra::ConstLinearOperator<double> Ivel 
+      = Thyra::identity<double>(Av_inverse.range());
+    Thyra::ConstLinearOperator<double> Ipress 
+      = Thyra::identity<double>(Mp_inverse.domain());
+
+                                       // Build zero operators. Need
+                                       // one that is pressure x
+                                       // velocity and one that is
+                                       // velocity x pressure
+    Thyra::ConstLinearOperator<double> zero;
+    
+                                       // Build the composed Schur
+                                       // complement approximation
+                                       // inverse
+    Thyra::ConstLinearOperator<double> Xinv;
+    if (additional_data.right_preconditioning == true)
+      Xinv = Mp_inverse * Fp_op * Lp_inverse;
+    else
+      Xinv = Lp_inverse * Fp_op * Mp_inverse;
+
+                                       // Create discrete pressure
+                                       // gradient and divergence
+                                       // operators in Thyra format.
+    Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpS01 = 
+      Thyra::epetraLinearOp(Teuchos::rcp(&*system_matrix.block(0,1).matrix, 
+                                        false));
+    Thyra::ConstLinearOperator<double> S01 = tmpS01;
+
+    Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpS10 = 
+      Thyra::epetraLinearOp(Teuchos::rcp(&*system_matrix.block(1,0).matrix, 
+                                        false));
+    Thyra::ConstLinearOperator<double> S10 = tmpS10;
+
+                                       // Build the 3 block operators
+                                       // for the
+                                       // preconditioner. Here we
+                                       // have to set different
+                                       // matrices depending on
+                                       // whether we do right or left
+                                       // preconditioning.
+    Thyra::ConstLinearOperator<double> P1 
+      = block2x2( Av_inverse, zero, zero, Ipress );
+
+    Thyra::ConstLinearOperator<double> P2;
+    if (additional_data.right_preconditioning == true)
+      P2 = block2x2( Ivel, (-1.0)*S01, zero, Ipress  );
+    else
+      P2 = block2x2( Ivel, zero, (-1.0)*S10, Ipress );
+
+    Thyra::ConstLinearOperator<double> P3 
+      = block2x2( Ivel,   zero, zero, (-1.0)*Xinv );
+
+    if (additional_data.right_preconditioning == true)    
+      preconditioner = Teuchos::rcp (new 
+                      Thyra::ConstLinearOperator<double>(P1 * P2 * P3));
+    else
+      preconditioner = Teuchos::rcp (new 
+                      Thyra::ConstLinearOperator<double>(P3 * P2 * P1));
+
+  }
+
+
+
+  void
+  PreconditionStokes::reinit_lazy ()
+  {
+    Assert (&*Av_precondition != 0,
+           ExcMessage ("Can only be done when the outer preconditioner owns"
+                       " the AMG object!"));
+
+    Av_precondition->reinit();
+  }
+
+} /* end of namespace TrilinosWrappers */
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
index 57ba88ac89cfbd6fc17a3b61de9c23e85d765a87..4ca0f6ed8c88cfb81c16f5cf1e89dd2f6ce0872c 100644 (file)
@@ -37,8 +37,17 @@ namespace TrilinosWrappers
   
   SolverBase::SolverBase (SolverControl  &cn)
                   :
-                  solver_control (cn),
-                 solver_name    (gmres)
+                 solver_name    (gmres),
+                  solver_control (cn)
+  {}
+
+  
+
+  SolverBase::SolverBase (const enum SolverBase::SolverName  solver_name,
+                         SolverControl                     &cn)
+                  :
+                 solver_name    (solver_name),
+                  solver_control (cn)
   {}
 
   
diff --git a/deal.II/lac/source/trilinos_solver_block.cc b/deal.II/lac/source/trilinos_solver_block.cc
new file mode 100644 (file)
index 0000000..122ef6e
--- /dev/null
@@ -0,0 +1,420 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/trilinos_solver_block.h>
+#include <lac/trilinos_block_sparse_matrix.h>
+#include <lac/trilinos_block_vector.h>
+#include <lac/trilinos_precondition_block.h>
+
+#include <base/utilities.h>
+
+#include <cmath>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Teuchos_RCP.hpp>
+#include <Teuchos_ParameterList.hpp>
+
+#include <Epetra_Vector.h>
+
+#include <Thyra_VectorDecl.hpp>
+#include <Thyra_VectorImpl.hpp>
+#include <Thyra_VectorSpaceImpl.hpp>
+#include <Thyra_MultiVectorBase.hpp>
+#include <Thyra_MultiVectorDefaultBase.hpp>
+#include <Thyra_LinearOperatorDecl.hpp>
+#include <Thyra_LinearOperatorImpl.hpp>
+#include <Thyra_DefaultBlockedLinearOpDecl.hpp>
+#include <Thyra_DefaultBlockedLinearOp.hpp>
+#include <Thyra_EpetraLinearOp.hpp>
+#include <Thyra_EpetraThyraWrappers.hpp>
+#include <Thyra_LinearOpWithSolveFactoryHelpers.hpp>
+#include <Thyra_AztecOOLinearOpWithSolveFactory.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+
+  SolverBlockBase::AdditionalData::
+  AdditionalData (const unsigned int gmres_restart_parameter,
+                 const bool         right_preconditioning,
+                 const bool         output_details)
+                  :
+                  gmres_restart_parameter (gmres_restart_parameter),
+                 right_preconditioning   (right_preconditioning),
+                 output_details          (output_details)
+  {}
+
+  
+  
+  SolverBlockBase::SolverBlockBase (SolverControl  &cn)
+                  :
+                  solver_control (cn),
+                 solver_name    (gmres)
+  {}
+
+  
+
+  SolverBlockBase::~SolverBlockBase ()
+  {}
+
+  
+
+  SolverControl &
+  SolverBlockBase::control() const
+  {
+    return solver_control;
+  }
+
+
+
+  void
+  SolverBlockBase::solve (const BlockSparseMatrix     &A,
+                         BlockVector                 &x,
+                         const BlockVector           &b,
+                         const PreconditionBlockBase &preconditioner)
+  {
+
+    Assert (x.n_blocks() == b.n_blocks(),
+           ExcDimensionMismatch(x.n_blocks(), b.n_blocks()));
+
+                                       // Just copy everything over
+                                       // to a distributed vector...
+    MPI::BlockVector tmpx (x.n_blocks());
+    MPI::BlockVector tmpb (b.n_blocks());
+    for (unsigned int i=0; i<x.n_blocks(); ++i)
+      {
+       tmpx.block(i).reinit (A.block(i,i).matrix->DomainMap(), false);
+       tmpx.block(i) = x.block(i);
+       tmpb.block(i).reinit (A.block(i,i).matrix->RangeMap(), false);
+       tmpb.block(i) = b.block(i);
+      }
+    tmpx.collect_sizes();
+    tmpb.collect_sizes();
+
+    solve (A, tmpx, tmpb, preconditioner);
+    x = tmpx;
+  }
+  
+    
+
+  void
+  SolverBlockBase::solve (const BlockSparseMatrix     &input_A,
+                         MPI::BlockVector            &input_x,
+                         const MPI::BlockVector      &input_b,
+                         const PreconditionBlockBase &preconditioner)
+  {
+
+                                       // Get block information about
+                                       // the preconditioner.
+    Thyra::ConstLinearOperator<double> P = *preconditioner.preconditioner;
+    const unsigned int n_rows = P.numBlockRows();
+    const unsigned int n_cols = P.numBlockCols();
+
+                                       // Check the dimensions of the
+                                       // block vectors and matrix
+                                       // and check if the parallel
+                                       // distribution of the
+                                       // matrices and vectors is the
+                                       // same.
+    AssertThrow (n_rows == n_cols,
+                ExcDimensionMismatch (n_rows, n_cols));
+
+    AssertThrow (input_x.n_blocks() == input_A.n_block_cols(),
+                ExcDimensionMismatch(input_x.n_blocks(),
+                                     input_A.n_block_cols()));
+    AssertThrow (input_b.n_blocks() == input_A.n_block_rows(),
+                ExcDimensionMismatch(input_b.n_blocks(),
+                                     input_A.n_block_rows()));
+    AssertThrow (input_A.n_block_rows() == n_rows,
+                ExcDimensionMismatch(input_A.n_block_rows(), n_rows));
+    AssertThrow (input_A.n_block_cols() == n_cols,
+                ExcDimensionMismatch(input_A.n_block_cols(), n_cols));
+
+    for (unsigned int i=0; i<n_rows; ++i)
+      {
+       AssertThrow (input_x.block(i).vector->Map().UniqueGIDs() == true,
+                    ExcOverlappingMaps("vector", "x"));
+       AssertThrow (input_b.block(i).vector->Map().UniqueGIDs() == true,
+                    ExcOverlappingMaps("vector", "b"));
+       for (unsigned int j=0; j<n_cols; ++j)
+         {
+           {
+             std::ostringstream i_str;
+             i_str << i;
+             std::ostringstream error_component;
+             error_component << "x.block(" << i_str.str() << ")";
+             AssertThrow (input_x.block(j).vector->Map().SameAs(
+                            input_A.block(i,j).matrix->DomainMap()) == true,
+                          ExcNonMatchingMaps (error_component.str()));
+           }
+           {
+             std::ostringstream i_str;
+             i_str << j;
+             std::ostringstream error_component;
+             error_component << "b.block(" << i_str.str() << ")";
+             AssertThrow (input_b.block(i).vector->Map().SameAs(
+                            input_A.block(i,j).matrix->RangeMap()) == true,
+                          ExcNonMatchingMaps (error_component.str()));
+           }
+         }
+      }
+
+
+                                       // Wrap Epetra vectors into
+                                       // Thyra vectors in order to
+                                       // be able to work with the
+                                       // abstract Thyra/Stratimikos
+                                       // interfaces to the AztecOO
+                                       // solver.
+                                       // 
+                                       // As a first step, create a
+                                       // copy of the vector spaces.
+    std::vector<Teuchos::RCP<const Thyra::VectorSpaceBase<double> > > 
+      epetra_vector_spaces;
+    for (unsigned int i=0; i<n_rows; ++i)
+      {
+       Teuchos::RCP<const Thyra::VectorSpaceBase<double> > tmp_space
+         = Thyra::create_VectorSpace(
+                   Teuchos::rcp(&input_A.block(i,i).matrix->DomainMap(), 
+                                false));
+
+       epetra_vector_spaces.push_back(tmp_space);
+      }
+
+                                       // Convert the block matrix to
+                                       // a Thyra linear operator
+                                       // that acts on blocks.
+    Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<double> >
+      tmpA = Teuchos::rcp(new Thyra::DefaultBlockedLinearOp<double>());
+    tmpA->beginBlockFill(n_rows,n_cols);
+    for (unsigned int i=0; i<n_rows; ++i)
+      for (unsigned int j=0; j<n_cols; ++j)
+       {
+         Teuchos::RCP<const Thyra::LinearOpBase<double> > 
+           A_ij = Thyra::epetraLinearOp(Teuchos::rcp(&*input_A.block(i,j).matrix, 
+                                                     false));
+         tmpA->setBlock(i, j, A_ij);
+       }
+    tmpA->endBlockFill();
+    tmpA->setObjectLabel("Blocked system matrix");
+    Thyra::RCP<const Thyra::LinearOpBase<double> > baseA = tmpA;
+
+    Thyra::ConstLinearOperator<double> A = baseA;
+
+                                       // Convert block vectors for
+                                       // right hand side and
+                                       // solution to Thyra objects.
+    Thyra::Vector<double> rhs = A.range().createMember();
+    Thyra::Vector<double> sol = A.domain().createMember();
+    for (unsigned int i=0; i<n_rows; ++i)
+      {
+       Epetra_Vector *block_ptr = (*input_b.block(i).vector)(0);
+       Teuchos::RCP<Thyra::VectorBase<double> > tmp_rhs_i
+         = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), 
+                                epetra_vector_spaces[i]);
+       const Thyra::Vector<double> rhs_i = tmp_rhs_i;
+       rhs.setBlock(i, rhs_i);
+
+       block_ptr = (*input_x.block(i).vector)(0);
+       Teuchos::RCP<Thyra::VectorBase<double> > tmp_sol_i
+         = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), 
+                                epetra_vector_spaces[i]);
+       Thyra::Vector<double> sol_i = tmp_sol_i;
+       sol.setBlock(i, sol_i);
+      }
+
+    // ------- Now build a solver factory for the block problem ------- 
+
+                                       // Set up a parameter list for
+                                       // the AztecOO solver.
+    Teuchos::RCP<Teuchos::ParameterList> aztecBlockParams 
+      = rcp(new Teuchos::ParameterList("AztecOO block solver parameters"));
+      
+                                       // Set the number of
+                                       // iterations and the solver
+                                       // tolerance from the solver
+                                       // control object.
+    aztecBlockParams->sublist("Forward Solve").set("Max Iterations", 
+                                                  (int)solver_control.max_steps());
+    aztecBlockParams->sublist("Forward Solve").set("Tolerance",
+                                                   solver_control.tolerance());
+    aztecBlockParams->sublist("Forward Solve")
+      .sublist("AztecOO Settings").set("Convergence Test", "no scaling");
+
+                                       // The settings about the
+                                       // actual solver that is going
+                                       // to be used - let this be
+                                       // determined by the
+                                       // <tt>solver_name</tt> set in
+                                       // a derived class.
+    switch (solver_name)
+      {
+        case cg:
+         aztecBlockParams->sublist("Forward Solve")
+           .sublist("AztecOO Settings").set("Aztec Solver", "CG");
+         break;
+        case gmres:
+         aztecBlockParams->sublist("Forward Solve")
+           .sublist("AztecOO Settings").set("Aztec Solver", "GMRES");
+         aztecBlockParams->sublist("Forward Solve")
+           .sublist("AztecOO Settings").set("Size of Krylov Subspace", 
+                                            (int)additional_data.gmres_restart_parameter);
+         break;
+        default:
+         ExcNotImplemented();
+      }
+
+                                       // We use an externally
+                                       // provided preconditioner, so
+                                       // do not set anything here.
+    aztecBlockParams->sublist("Forward Solve")
+      .sublist("AztecOO Settings").set("Aztec Preconditioner", "none");
+
+    if (additional_data.output_details == true)
+      aztecBlockParams->sublist("Forward Solve")
+       .sublist("AztecOO Settings").set("Output Frequency", 1);
+    else
+      aztecBlockParams->sublist("VerboseObject").set("Verbosity Level", 
+                                                    "none");
+
+                                       // Create solve factory from
+                                       // parameter list.
+    Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
+      aztecBlockLowsFactory = Teuchos::rcp
+        (new Thyra::AztecOOLinearOpWithSolveFactory());
+
+    aztecBlockLowsFactory->setParameterList(aztecBlockParams);
+
+    Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> > rcpBlockAztec 
+      = aztecBlockLowsFactory->createOp();
+
+
+    // ---- Set up the preconditioned inverse object and do the solve ----
+
+                                       // Extract base operators from
+                                       // block matrix A and block
+                                       // preconditioner P.
+    Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > A_ptr 
+      = Teuchos::rcp(new Thyra::DefaultLinearOpSource<double>(A.constPtr()), 
+                    true);
+    Teuchos::RCP<Thyra::PreconditionerBase<double> > P_ptr = 
+      Teuchos::rcp (new Thyra::DefaultPreconditioner<double>(), true);
+    {
+      Thyra::DefaultPreconditioner<double>
+       *defaultPrec = &Teuchos::dyn_cast<Thyra::DefaultPreconditioner<double> >
+                     (*P_ptr);
+      if (additional_data.right_preconditioning == true)
+       (*defaultPrec).initializeRight(P.constPtr());
+      else
+       (*defaultPrec).initializeUnspecified(P.constPtr());
+    }
+
+                                       // Attach the matrix and the
+                                       // preconditioner to the
+                                       // solver object.
+    aztecBlockLowsFactory->initializePreconditionedOp(A_ptr, P_ptr, 
+                                                     &*rcpBlockAztec, 
+                                                     Thyra::SUPPORT_SOLVE_FORWARD_ONLY);
+
+                                       // Now do the solve.
+    Thyra::SolveStatus<double>
+      solveStatus = Thyra::solve( *rcpBlockAztec, Thyra::NOTRANS, 
+                                 *rhs.constPtr().get(), &*sol.ptr().get() );
+
+                                       // Extract the number of
+                                       // needed iterations from the
+                                       // solver message.
+    unsigned int n_iterations;
+    {
+      std::string output_text = solveStatus.message;
+      unsigned int position = output_text.find ("using ");
+      const std::pair<int,unsigned int> tmp 
+       = Utilities::get_integer_at_position (output_text, position+6);
+      n_iterations = tmp.first;
+    }
+
+                                       // Do a dummy check of
+                                       // convergence to tell the
+                                       // SolverControl object the
+                                       // number of iterations. We do
+                                       // not impose the achieved
+                                       // tolerance here since that
+                                       // would require additional
+                                       // work (computing a
+                                       // residual).
+    solver_control.check (n_iterations, 0);
+
+  }
+  
+    
+
+
+  
+/* ------------------------ SolverCG -------------------------- */
+
+  SolverBlockCG::AdditionalData::
+  AdditionalData (const bool         right_preconditioning,
+                 const bool         output_details)
+                  :
+                  right_preconditioning   (right_preconditioning),
+                 output_details          (output_details)
+  {}
+
+  
+  
+  SolverBlockCG::SolverBlockCG (SolverControl        &cn,
+                               const AdditionalData &data)
+                  :
+                  SolverBlockBase (cn),
+                  data (data)
+  {
+    solver_name = cg;
+    additional_data.right_preconditioning = data.right_preconditioning;
+    additional_data.output_details = data.output_details;
+  }
+  
+
+/* ---------------------- SolverGMRES ------------------------- */
+
+  SolverBlockGMRES::AdditionalData::
+  AdditionalData (const unsigned int gmres_restart_parameter,
+                 const bool         right_preconditioning,
+                 const bool         output_details)
+                  :
+                  gmres_restart_parameter (gmres_restart_parameter),
+                 right_preconditioning   (right_preconditioning),
+                 output_details          (output_details)
+  {}
+
+  
+
+  SolverBlockGMRES::SolverBlockGMRES (SolverControl        &cn,
+                                     const AdditionalData &data)
+                  :
+                  SolverBlockBase (cn),
+                  data (data)
+  {
+    solver_name = gmres;
+    additional_data.gmres_restart_parameter = data.gmres_restart_parameter;
+    additional_data.right_preconditioning = data.right_preconditioning;
+    additional_data.output_details = data.output_details;
+  }
+
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_PETSC

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.