From 5c2d5c72f75351608c5ec68aedc249f2897eb51d Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 23 Nov 2009 18:15:25 +0000 Subject: [PATCH] Remove an implementation of a Stokes preconditioner that was too specific and never really worked. git-svn-id: https://svn.dealii.org/trunk@20154 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 10 + .../include/lac/trilinos_precondition_block.h | 374 +---------------- .../lac/source/trilinos_precondition_block.cc | 386 ------------------ 3 files changed, 12 insertions(+), 758 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 3e7848ee82..08dcf916c0 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -30,6 +30,16 @@ inconvenience this causes.

    +
  1. +

    + 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. +
    + (Martin Kronbichler 2009/11/23) +

    +
  2. Changed: The class MGSmootherRelaxation now instead of a diff --git a/deal.II/lac/include/lac/trilinos_precondition_block.h b/deal.II/lac/include/lac/trilinos_precondition_block.h index b3dc161707..65e6070455 100644 --- a/deal.II/lac/include/lac/trilinos_precondition_block.h +++ b/deal.II/lac/include/lac/trilinos_precondition_block.h @@ -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 > preconditioner; friend class SolverBlockBase; -}; - - /** - * Internal function that sets - * up an inverse matrix. - */ - inline - Thyra::ConstLinearOperator - 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 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) - * and D. Kay, D. Loghin, A. Wathen, A preconditioner for the - * steady-state Navier–Stokes equations, SIAM - * J. Sci. Comput. 24(1):237–256 (2002), 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 - * initialize 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 > &Av_constant_modes = std::vector > (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 > 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: 1e-9 - * 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 - * Av, that can be - * specified to be different - * from the (0,0) block in the - * system matrix, IC/SSOR on - * the pressure mass matrix - * Mp_matrix, and AMG - * on the pressure Laplace - * matrix - * Lp_matrix. 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 - * initialize - * 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 Mp_precondition_ssor; - - /** - * Pointer to preconditioner - * for the mass matrix. - */ - Teuchos::RCP Mp_precondition_ic; - - /** - * Pointer to preconditioner - * for the velocity-velocity - * matrix. - */ - Teuchos::RCP Av_precondition; - - /** - * Pointer to preconditioner - * for the pressure Laplace - * matrix. - */ - Teuchos::RCP Lp_precondition; }; } diff --git a/deal.II/lac/source/trilinos_precondition_block.cc b/deal.II/lac/source/trilinos_precondition_block.cc index 61fab8a232..ece2a4a611 100644 --- a/deal.II/lac/source/trilinos_precondition_block.cc +++ b/deal.II/lac/source/trilinos_precondition_block.cc @@ -19,16 +19,6 @@ # include # include -# include -# include - -# include -# include -# include -# include -# include -# include -# include DEAL_II_NAMESPACE_OPEN @@ -43,382 +33,6 @@ namespace TrilinosWrappers PreconditionBlockBase::~PreconditionBlockBase() {} - - - /** - * Internal function that sets - * up an inverse matrix. - */ - inline - Thyra::ConstLinearOperator - 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 > tmpM = - Thyra::epetraLinearOp(Teuchos::rcp(&input_M.trilinos_matrix(), false)); - - Teuchos::RCP > M - = Teuchos::rcp(new Thyra::DefaultLinearOpSource(tmpM), true); - - //Thyra::ConstLinearOperator opM = tmpM; - - // Create a Thyra version of - // the preconditioner. - Teuchos::RCP > tmpP = - Thyra::epetraLinearOp(Teuchos::rcp(const_cast(input_P), - false), - Thyra::NOTRANS, - Thyra::EPETRA_OP_APPLY_APPLY_INVERSE, - Thyra::EPETRA_OP_ADJOINT_SUPPORTED); - - //Thyra::ConstLinearOperator opPM = tmpP; - - Teuchos::RCP > PM = - Teuchos::rcp (new Thyra::DefaultPreconditioner(), true); - - { - Thyra::DefaultPreconditioner - *defaultPrec = &Teuchos::dyn_cast > - (*PM); - - (*defaultPrec).initializeUnspecified(tmpP); - } - - // Set up the solver. - Teuchos::RCP 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 > - aztecLowsFactory = Teuchos::rcp(new Thyra::AztecOOLinearOpWithSolveFactory(), true); - - aztecLowsFactory->setParameterList(aztecParams); - - Teuchos::RCP > 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 azM - // = new Meros::AztecSolveStrategy(*(aztecParams.get())); - //Thyra::LinearOperator Minverse - // = Thyra::inverse(*aztecLowsFactory, opM); - //return new Meros::InverseOperator(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 > Minverse = - Teuchos::rcp(new Thyra::DefaultInverseLinearOp(rcpAztec)); - - return Minverse; - } - - - - PreconditionStokes::AdditionalData:: - AdditionalData (const bool right_preconditioning, - const bool Av_is_symmetric, - const std::vector > &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 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 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 > - 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 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 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 (Mp_inverse.range()); - - // Create Fp Thyra operator. - Thyra::ConstLinearOperator Fp_op; - if (Fp_matrix.m() != system_matrix.block(1,1).m()) - { - Fp_op = Thyra::identity (Mp_inverse.range()); - } - else - { - Teuchos::RCP > - 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 Ivel - = Thyra::identity(Av_inverse.range()); - Thyra::ConstLinearOperator Ipress - = Thyra::identity(Mp_inverse.domain()); - - // Build zero operators. Need - // one that is pressure x - // velocity and one that is - // velocity x pressure - Thyra::ConstLinearOperator zero; - - // Build the composed Schur - // complement approximation - // inverse - Thyra::ConstLinearOperator 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 > tmpS01 = - Thyra::epetraLinearOp(Teuchos::rcp(&system_matrix.block(0,1).trilinos_matrix(), - false)); - Thyra::ConstLinearOperator S01 = tmpS01; - - Teuchos::RCP > tmpS10 = - Thyra::epetraLinearOp(Teuchos::rcp(&system_matrix.block(1,0).trilinos_matrix(), - false)); - Thyra::ConstLinearOperator 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 P1 - = block2x2( Av_inverse, zero, zero, Ipress ); - - Thyra::ConstLinearOperator 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 P3 - = block2x2( Ivel, zero, zero, (-1.0)*Xinv ); - - if (additional_data.right_preconditioning == true) - preconditioner = Teuchos::rcp (new - Thyra::ConstLinearOperator(P1 * P2 * P3)); - else - preconditioner = Teuchos::rcp (new - Thyra::ConstLinearOperator(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 -- 2.39.5