// $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
*
* @ingroup TrilinosWrappers
* @ingroup Preconditioners
- * @author Martin Kronbichler, 2008
+ * @author Martin Kronbichler, 2008, 2009
*/
class PreconditionBlockBase : public Subscriptor
{
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–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–1367 (1994)</em>
- * and <em> D. Kay, D. Loghin, A. Wathen, A preconditioner for the
- * steady-state Navier–Stokes equations, SIAM
- * J. Sci. Comput. 24(1):237–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–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;
};
}
# 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
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