From d16d88cf884915376c0477dc0d837e5240dcffa2 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 13 Oct 2008 09:20:41 +0000 Subject: [PATCH] Added interfaces to Trilinos block solvers and a block preconditioner (for Stokes systems). git-svn-id: https://svn.dealii.org/trunk@17188 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_precondition.h | 5 +- .../include/lac/trilinos_precondition_block.h | 448 ++++++++++++++++++ deal.II/lac/include/lac/trilinos_solver.h | 49 +- .../lac/include/lac/trilinos_solver_block.h | 366 ++++++++++++++ deal.II/lac/source/trilinos_precondition.cc | 8 +- .../lac/source/trilinos_precondition_block.cc | 423 +++++++++++++++++ deal.II/lac/source/trilinos_solver.cc | 13 +- deal.II/lac/source/trilinos_solver_block.cc | 420 ++++++++++++++++ 8 files changed, 1709 insertions(+), 23 deletions(-) create mode 100644 deal.II/lac/include/lac/trilinos_precondition_block.h create mode 100644 deal.II/lac/include/lac/trilinos_solver_block.h create mode 100755 deal.II/lac/source/trilinos_precondition_block.cc create mode 100644 deal.II/lac/source/trilinos_solver_block.cc diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index 54db649b4d..ede32aa723 100755 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -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 index 0000000000..25aff1116f --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_precondition_block.h @@ -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 +#include +#include +#include +#include + + +#ifdef DEAL_II_USE_TRILINOS + +#include +#include +#include + +// 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 > 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&ndash1367 (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. + * + * @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_null_space = 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_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: 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 @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; + }; +} + +/*@}*/ + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS + +/*-------------------------- trilinos_precondition_block.h ---------------------------*/ + +#endif +/*-------------------------- trilinos_precondition_block.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/trilinos_solver.h b/deal.II/lac/include/lac/trilinos_solver.h index fbddcef679..f7a01c17f1 100644 --- a/deal.II/lac/include/lac/trilinos_solver.h +++ b/deal.II/lac/include/lac/trilinos_solver.h @@ -47,6 +47,12 @@ namespace TrilinosWrappers * "http://trilinos.sandia.gov/packages/aztecoo/AztecOOUserGuide.pdf">AztecOO * user guide. * + * This solver class can also be used as a standalone class, where the + * respective Krylov method is set via the flag + * solver_name. 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 "cg", "cgs", - * "gmres", "bicgstab", - * and "tfqmr". - */ - 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 index 0000000000..5b13a870e3 --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_solver_block.h @@ -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 +#include +#include + +#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 + * Ax=b. 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 + * Ax=b. 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 "cg", + * "gmres". + */ + 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 -----------------------*/ diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index b09cae7295..acf6fcfb67 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -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 index 0000000000..433a891a78 --- /dev/null +++ b/deal.II/lac/source/trilinos_precondition_block.cc @@ -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 + +#ifdef DEAL_II_USE_TRILINOS + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + + PreconditionBlockBase::PreconditionBlockBase() + {} + + + + 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 = + Teuchos::rcp(new Thyra::EpetraLinearOp(Teuchos::rcp(&*input_M.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 = + Teuchos::rcp(new 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).initializeRight(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_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 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 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.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).matrix, + false)); + Thyra::ConstLinearOperator S01 = tmpS01; + + Teuchos::RCP > tmpS10 = + Thyra::epetraLinearOp(Teuchos::rcp(&*system_matrix.block(1,0).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 + +#endif // DEAL_II_USE_TRILINOS diff --git a/deal.II/lac/source/trilinos_solver.cc b/deal.II/lac/source/trilinos_solver.cc index 57ba88ac89..4ca0f6ed8c 100644 --- a/deal.II/lac/source/trilinos_solver.cc +++ b/deal.II/lac/source/trilinos_solver.cc @@ -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 index 0000000000..122ef6ea75 --- /dev/null +++ b/deal.II/lac/source/trilinos_solver_block.cc @@ -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 +#include +#include +#include + +#include + +#include + +#ifdef DEAL_II_USE_TRILINOS + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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; iDomainMap(), 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 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; iMap().UniqueGIDs() == true, + ExcOverlappingMaps("vector", "x")); + AssertThrow (input_b.block(i).vector->Map().UniqueGIDs() == true, + ExcOverlappingMaps("vector", "b")); + for (unsigned int j=0; jMap().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 > > + epetra_vector_spaces; + for (unsigned int i=0; i > 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 > + tmpA = Teuchos::rcp(new Thyra::DefaultBlockedLinearOp()); + tmpA->beginBlockFill(n_rows,n_cols); + for (unsigned int i=0; i > + 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 > baseA = tmpA; + + Thyra::ConstLinearOperator A = baseA; + + // Convert block vectors for + // right hand side and + // solution to Thyra objects. + Thyra::Vector rhs = A.range().createMember(); + Thyra::Vector sol = A.domain().createMember(); + for (unsigned int i=0; i > tmp_rhs_i + = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), + epetra_vector_spaces[i]); + const Thyra::Vector rhs_i = tmp_rhs_i; + rhs.setBlock(i, rhs_i); + + block_ptr = (*input_x.block(i).vector)(0); + Teuchos::RCP > tmp_sol_i + = Thyra::create_Vector(Teuchos::rcp(block_ptr, false), + epetra_vector_spaces[i]); + Thyra::Vector 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 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 + // solver_name 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 > + aztecBlockLowsFactory = Teuchos::rcp + (new Thyra::AztecOOLinearOpWithSolveFactory()); + + aztecBlockLowsFactory->setParameterList(aztecBlockParams); + + Teuchos::RCP > 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 > A_ptr + = Teuchos::rcp(new Thyra::DefaultLinearOpSource(A.constPtr()), + true); + Teuchos::RCP > P_ptr = + Teuchos::rcp (new Thyra::DefaultPreconditioner(), true); + { + Thyra::DefaultPreconditioner + *defaultPrec = &Teuchos::dyn_cast > + (*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 + 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 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 -- 2.39.5