]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove an implementation of a Stokes preconditioner that was too specific and never...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 23 Nov 2009 18:15:25 +0000 (18:15 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 23 Nov 2009 18:15:25 +0000 (18:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@20154 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/lac/include/lac/trilinos_precondition_block.h
deal.II/lac/source/trilinos_precondition_block.cc

index 3e7848ee826e18474bc763b2158326810bb32dde..08dcf916c04b67f7197ec61f35b33fde47c2dd3a 100644 (file)
@@ -30,6 +30,16 @@ inconvenience this causes.
 </p>
 
 <ol>
+  <li>
+  <p>
+  Removed: The class TrilinosWrappers::PreconditionStokes has been deleted
+  because it did not work properly, and it is too specific to be part of
+  the general library. A preconditioner that has the same properties is 
+  explained in the @ref step_31 step-31 tutorial program.
+  <br>
+  (Martin Kronbichler 2009/11/23)
+  </p>
+
   <li>
   <p>
   Changed: The class MGSmootherRelaxation now instead of a
index b3dc161707fbb773478e760dcd17e1d72fde0441..65e60704554b6ef846b8cd76317281faa59d176d 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2008 by the deal.II authors
+//    Copyright (C) 2008, 2009 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -51,7 +51,7 @@ namespace TrilinosWrappers
  *
  * @ingroup TrilinosWrappers
  * @ingroup Preconditioners
- * @author Martin Kronbichler, 2008
+ * @author Martin Kronbichler, 2008, 2009
  */
   class PreconditionBlockBase : public Subscriptor
   {
@@ -73,376 +73,6 @@ namespace TrilinosWrappers
       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&ndash;1367 (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.
- *
- * The preconditioner is based an approximation to the Schur
- * complement of the block matrix. The Schur complement $S=B
- * A_{\mathbf u}^{-1} B^T$ is approximated by a mass matrix $M_p$ on
- * the pressure space in the case of the Stokes equations, and as a
- * product $S^{-1} = L_p^{-1} F_p M_p^{-1}$ with pressure Laplace
- * matrix $L_p$, pressure convection-diffusion operator $F_p$
- * (corresponding to the sum of time derivative, convection and
- * diffusion), and pressure mass matrix $M_p$ in the case of the
- * Navier&ndash;Stokes equations.
- *
- * @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_constant_modes = 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_constant_modes;
-
-                                      /**
-                                       * 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 @ref 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;
   };
 }
 
index 61fab8a23204d5dd0da2b653bde23485963274aa..ece2a4a6111467e73fd0ab8a0af377157efd7e7e 100644 (file)
 #  include <lac/trilinos_sparse_matrix.h>
 #  include <lac/trilinos_block_sparse_matrix.h>
 
-#  include <Thyra_EpetraLinearOp.hpp>
-#  include <Thyra_EpetraThyraWrappers.hpp>
-
-#  include <Thyra_VectorSpaceImpl.hpp>
-#  include <Thyra_LinearOperatorDecl.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
 
@@ -43,382 +33,6 @@ namespace TrilinosWrappers
   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 = 
-      Thyra::epetraLinearOp(Teuchos::rcp(&input_M.trilinos_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 = 
-      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).initializeUnspecified(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_constant_modes,
-                 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_constant_modes (Av_constant_modes),
-                 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.constant_modes = additional_data.Av_constant_modes;
-
-      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.trilinos_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).trilinos_matrix(), 
-                                        false));
-    Thyra::ConstLinearOperator<double> S01 = tmpS01;
-
-    Teuchos::RCP<const Thyra::LinearOpBase<double> > tmpS10 = 
-      Thyra::epetraLinearOp(Teuchos::rcp(&system_matrix.block(1,0).trilinos_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

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.