From: Jan Philipp Thiele Date: Thu, 18 Jul 2024 16:42:37 +0000 (+0200) Subject: Add Tpetra wrappers for Trilinos Ifpack2 preconditioners X-Git-Tag: v9.6.0-rc1~67^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=794bf75c1fd4dea65b632de45bcf355f47e69fc5;p=dealii.git Add Tpetra wrappers for Trilinos Ifpack2 preconditioners --- diff --git a/cmake/configure/configure_20_trilinos.cmake b/cmake/configure/configure_20_trilinos.cmake index c80291be24..0c26dfa146 100644 --- a/cmake/configure/configure_20_trilinos.cmake +++ b/cmake/configure/configure_20_trilinos.cmake @@ -22,7 +22,7 @@ set(FEATURE_TRILINOS_DEPENDS MPI) # A list of optional Trilinos modules we use: # set(_deal_ii_trilinos_optional_modules - Amesos2 Belos EpetraExt Kokkos MueLu NOX ROL Sacado SEACAS Tpetra Zoltan + Amesos2 Belos EpetraExt Ifpack2 Kokkos MueLu NOX ROL Sacado SEACAS Tpetra Zoltan ) # diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 8892274e69..ebb2b05272 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -2630,3 +2630,13 @@ url = {http://sepwww.stanford.edu/oldsep/joe/Reprints/8IWSA.pdf}, } +@article{BFKY2011, + title={Multigrid smoothers for ultraparallel computing}, + author={Baker, Allison H and Falgout, Robert D and Kolev, Tzanio V and Yang, Ulrike Meier}, + journal={SIAM Journal on Scientific Computing}, + volume={33}, + number={5}, + pages={2864--2887}, + year={2011}, + publisher={SIAM} +} diff --git a/doc/news/changes/major/20240718Thiele b/doc/news/changes/major/20240718Thiele new file mode 100644 index 0000000000..670190c8ad --- /dev/null +++ b/doc/news/changes/major/20240718Thiele @@ -0,0 +1,30 @@ +New: LinearAlgebra::TpetraWrappers preconditioner interfaces +for Trilinos Ifpack2, as listed below: + +(Jan Philipp Thiele, 2024/07/18) \ No newline at end of file diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index 5546e052c8..e09209db4d 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -295,6 +295,7 @@ #cmakedefine DEAL_II_TRILINOS_WITH_AMESOS2 #cmakedefine DEAL_II_TRILINOS_WITH_BELOS #cmakedefine DEAL_II_TRILINOS_WITH_EPETRAEXT +#cmakedefine DEAL_II_TRILINOS_WITH_IFPACK2 #cmakedefine DEAL_II_TRILINOS_WITH_MUELU #cmakedefine DEAL_II_TRILINOS_WITH_NOX #cmakedefine DEAL_II_TRILINOS_WITH_ROL diff --git a/include/deal.II/lac/trilinos_tpetra_precondition.h b/include/deal.II/lac/trilinos_tpetra_precondition.h new file mode 100644 index 0000000000..9650d732bf --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_precondition.h @@ -0,0 +1,1329 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +#ifndef dealii_trilinos_tpetra_precondition_h +#define dealii_trilinos_tpetra_precondition_h + + +#include + +#include "deal.II/base/exceptions.h" +#include "deal.II/base/memory_space.h" + +#include "deal.II/lac/vector.h" + +#include + + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + +# include + +# include +# include + +# include +# include +# include +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + + /** + * The base class for all preconditioners based on Tpetra sparse matrices. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionBase : public Subscriptor + { + public: + /** + * @brief Standardized data struct to pipe additional flags to the + * preconditioner. + * + */ + struct AdditionalData + {}; + + /** + * @brief Constructor. Does not do anything. The initialize function of + * the derived classes will have to create the preconditioner from a given + * sparse matrix. + * + */ + PreconditionBase() = default; + + /** + * @brief Desctructor. + * Destroys the preconditioner, leaving an object like just after having + * called the constructor. + */ + void + clear(); + + /** + * @brief Apply the preconditioner. + * + * @param dst Input vector to apply the preconditioner to + * @param src Result vector + */ + virtual void + vmult(Vector &dst, + const Vector &src) const; + + /** + * @brief Apply the transpose preconditioner + * + * @param dst Input vector to apply the preconditioner to + * @param src Result vector + */ + virtual void + Tvmult(Vector &dst, + const Vector &src) const; + + /** + * @brief Apply the preconditioner + * + * @param dst Input vector to apply the preconditioner to + * @param src Result vector + */ + virtual void + vmult(dealii::Vector &dst, dealii::Vector &src) const; + + /** + * @brief Apply the transpose preconditioner + */ + virtual void + Tvmult(dealii::Vector &dst, dealii::Vector &src) const; + + /** + * @brief Access to underlying Trilinos data + * + * Calling this function from an uninitialized object will cause an + * exception. + */ + const TpetraTypes::LinearOperator & + trilinos_operator() const; + + /** + * @brief Access to underlying Trilinos data + * + * Calling this function from an uninitialized object will cause an + * exception. + */ + Teuchos::RCP> + trilinos_rcp() const; + + /** + * @name Partitioners + */ + /** @{ */ + + /** + * Return the partitioning of the domain space of this matrix, i.e., the + * partitioning of the vectors this matrix has to be multiplied with. + * + * @return IndexSet of the domain of this preconditioner + */ + IndexSet + locally_owned_domain_indices() const; + + /** + * Return the partitioning of the range space of this matrix, i.e., the + * partitioning of the vectors that are result from matrix-vector + * products. + * + * @return IndexSet of the range of this preconditioner + */ + IndexSet + locally_owned_range_indices() const; + + /** @} */ + + /** + * @addtogroup Exceptions + */ + /** @{ */ + /** + * The maps of the underlying matrix and one of the vectors + * do not match. + */ + 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."); + + /** + * The chosen preconditioner does not support a transposed apply. + */ + DeclExceptionMsg( + ExcTransposeNotSupported, + "The chosen preconditioner does not support transposing the matrix."); + /** @} */ + + protected: + /** + * This is a pointer to the preconditioner object that is used when + * applying the preconditioner. + */ + Teuchos::RCP> + preconditioner; + + /** + * @brief The list of preconditioner parameters. + * + * This structure is Trilinos counterpart to the AdditionalData structures + * in deal.II. Therefore any initialize will at some point pass this to + * the preconditioner. Most derived classes will handle building this list + * based on an AdditionalData object that exposes and defaults the most + * sensible parameters. But some classes will also offer full + * customization for experienced Trilinos users. + * + */ + Teuchos::ParameterList parameter_list; + }; + + + +# ifdef DEAL_II_TRILINOS_WITH_IFPACK2 + /** + * @brief Wrapper class for the IdentitySolver preconditioner of Ifpack2. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + * + * @tparam Number scalar type of the preconditioner + * @tparam MemorySpace + */ + template + class PreconditionIdentity : public PreconditionBase + { + public: + /** + * @brief Construct identity preconditioner. + * + */ + PreconditionIdentity() = default; + + /** + * @brief Initializes the preconditioner for the matrix A. + * + * Note that this only needs the matrix for information on the parallel + * distribution of the matrix and will return the original vector on + * vmult or Tvmult. + * + * @param A Matrix to base the preconditioner on. + */ + void + initialize(const SparseMatrix &A); + }; + + + + /** + * @brief The base class for all Ifpack2 preconditioners which are handled through its Factory. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionIfpackBase : public PreconditionBase + { + public: + /** + * @brief Constructor. + * + * This constructor will only save the type name given. + * The actual preconditioner has to be constructed through the + * initialize method of one of the underlying classes. + * + * @param preconditioner_type + */ + PreconditionIfpackBase(const std::string &preconditioner_type); + + /** + * Initializes the preconditioner for the matrix A based on + * the parameter_set. + * + */ + void + initialize(const SparseMatrix &A); + + /** + * The chosen preconditioner is not supported or configured with Ifpack2. + */ + DeclException1( + ExcTrilinosIpack2PreconditionerUnsupported, + std::string, + << "You tried to select the preconditioner type <" << arg1 << ">\n" + << "but this preconditioner is not supported by Trilinos/Ifpack22\n" + << "due to one of the following reasons:\n" + << "* This preconditioner does not exist\n" + << "* This preconditioner has a specialized constructor not supported by the Ifpack2 Factory.\n" + << "* This preconditioner is not (yet) supported by Trilinos/Ifpack2\n" + << "* Trilinos/Ifpack2 was not configured for its use."); + + protected: + /** + * The set preconditioner type to be handed to + * Ifpack2::Factory::create() + */ + std::string preconditioner_type; + }; + + + + /** + * @brief Wrapper to create custom Ifpack2 preconditioners. + * + * This general purpose class supports any preconditioner of Ifpack2 + * that can be constructed through its Ifpack2::Factory::create() + * function. + * + * @note This is a class for users familiar with Trilinos. + * + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionIfpack + : public PreconditionIfpackBase + { + public: + /** + * @brief Construct a new custom Ifpack2 preconditioner. + * + * The currently tested options for the type are + *
    + *
  • "RELAXATION"
    • + *
  • "CHEBYSHEV"
    • + *
  • "RILUK"
    • + *
  • "ILUT"
    • + *
  • "SCHWARZ"
    • + *
+ * + * But please refer to the User's Guide of Ifpack2 for more information. + * + * @param preconditioner_type the type based on Ifpack2 notation + */ + PreconditionIfpack(const std::string &preconditioner_type); + + /** + * @brief Set the parameter list for the preconditioner. + * + * This list will be passed to the Ifpack2 preconditioner during + * initialization. For parameter options based on the chosen + * preconditioner type see the User's Guide of Ifpack2. + * + * @param parameter_list + */ + void + set_parameter_list(Teuchos::ParameterList ¶meter_list); + }; + + + + /** + * @brief The classical Jacobi preconditioner within Ifpack2. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + * + */ + template + class PreconditionJacobi + : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * The defaults yield a textbook Jacobi preconditioner. + * + * @param omega The damping factor for the relaxation. + * @param fix_diagonal Fix small diagonal entries for the inversion? + * @param min_diagonal Threshold for fix_diagonal. + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult. + */ + AdditionalData(const double omega = 1., + const bool fix_diagonal = false, + const double min_diagonal = 0., + const int n_sweeps = 1); + + /** + * @brief Relaxation damping factor. + * + */ + double omega; + /** + * @brief Whether or not to enlarge entries below a threshold + * + * If the matrix diagonal contains entries close to zero, their + * inversion will be close to division by zero. This can be corrected by + * setting this value to true. This will result in additional + * computational work during initialization. + */ + bool fix_diagonal; + /** + * @brief Threshold below which entries will be fixed. + * + * If the threshold is zero (default) only entries which are exactly + * zero will be replaced will small nonzero values. + * + */ + double min_diagonal; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + /** + * @brief Constructor. + * + */ + PreconditionJacobi(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * @brief l1 variant of the Jacobi preconditioner. + * + * This variant adds the l1-Norm of all off-processor entries of a local row + * i to its diagonal entry d_ii to improve coupling in MPI-parallel + * applications, as described and introduced in @cite BFKY2011. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionL1Jacobi + : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * The defaults yield the l1 Jacobi preconditioner specified in the + * original publication. + * + * @param omega The damping factor for the relaxation. + * @param fix_diagonal Fix small diagonal entries for the inversion? + * @param min_diagonal Threshold for fix_diagonal. + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult. + */ + AdditionalData(const double omega = 1., + const bool fix_diagonal = false, + const double min_diagonal = 0., + const int n_sweeps = 1); + + /** + * @brief Relaxation damping factor. + * + */ + double omega; + /** + * @brief Whether or not to enlarge entries below a threshold + * + * If the matrix diagonal contains entries close to zero, their + * inversion will be close to division by zero. This can be corrected by + * setting this value to true. This will result in additional + * computational work during initialization. + */ + bool fix_diagonal; + /** + * @brief Threshold below which entries will be fixed. + * + * If the threshold is zero (default) only entries which are exactly + * zero will be replaced will small nonzero values. + * + */ + double min_diagonal; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + + /** + * @brief Constructor. + * + */ + PreconditionL1Jacobi(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * @brief l1 variant of the Gauss-Seidel preconditioner. + * + * This variant adds the l1-Norm of all off-processor entries of a single + * row i of the upper triangular part of A to its diagonal entry d_ii to + * improve coupling in MPI-parallel applications, , as described and + * introduced in @cite BFKY2011. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionL1GaussSeidel + : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * The defaults yield the l1 Gauss-Seidel preconditioner specified in + * the original publication. + * + * @param omega The damping factor for the relaxation. + * @param eta Threshold parameter for diagonal correction. + * @param fix_diagonal Fix small diagonal entries for the inversion? + * @param min_diagonal Threshold for fix_diagonal. + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult. + */ + AdditionalData(const double omega = 1, + const double eta = 1.5, + const bool fix_diagonal = false, + const double min_diagonal = 0, + const int n_sweeps = 1); + + /** + * @brief Relaxation damping factor. + * + */ + double omega; + + /** + * @brief Threshold parameter for diagonal correction. + * + * The l1 correction for a row i will only be added if it is more than + * eta times larger than the original diagonal entry. + * + */ + double eta; + /** + * @brief Whether or not to enlarge entries below a threshold + * + * If the matrix diagonal contains entries close to zero, their + * inversion will be close to division by zero. This can be corrected by + * setting this value to true. This will result in additional + * computational work during initialization. + */ + bool fix_diagonal; + /** + * @brief Threshold below which entries will be fixed. + * + * If the threshold is zero (default) only entries which are exactly + * zero will be replaced will small nonzero values. + * + */ + double min_diagonal; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + + /** + * @brief Constructor. + * + */ + PreconditionL1GaussSeidel(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * @brief The class for the SOR preconditioner within Ifpack2. + * + * If the code is executed MPI-parallel the individual processors will be + * connected with an additive Schwarz method. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionSOR : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * The defaults yield a textbook Jacobi preconditioner. + * + * @param omega The damping factor for the relaxation. + * @param overlap Overlap between processor local matrices. + * @param fix_diagonal Fix small diagonal entries for the inversion? + * @param min_diagonal Threshold for fix_diagonal. + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult. + */ + AdditionalData(const double omega = 1., + const int overlap = 0, + const bool fix_diagonal = false, + const double min_diagonal = 0., + const int n_sweeps = 1); + + /** + * @brief Relaxation damping factor. + * + */ + double omega; + /** + * @brief Overlap between processor local matrices. + * + * The amount of overlap between the processor local matrix blocks + * used in the underlying additive Schwarz method. + * + * The default will yield a block Jacobi preconditioner with each + * processor forming its own local block. + */ + int overlap; + /** + * @brief Whether or not to enlarge entries below a threshold + * + * If the matrix diagonal contains entries close to zero, their + * inversion will be close to division by zero. This can be corrected by + * setting this value to true. This will result in additional + * computational work during initialization. + */ + bool fix_diagonal; + /** + * @brief Threshold below which entries will be fixed. + * + * If the threshold is zero (default) only entries which are exactly + * zero will be replaced will small nonzero values. + * + */ + double min_diagonal; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + + /** + * @brief Constructor. + * + */ + PreconditionSOR(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * The class for the classical SSOR preconditioner within Ifpack2. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionSSOR : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * The defaults yield a textbook Jacobi preconditioner. + * + * @param omega The damping factor for the relaxation. + * @param overlap Overlap between processor local matrices. + * @param fix_diagonal Fix small diagonal entries for the inversion? + * @param min_diagonal Threshold for fix_diagonal. + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult. + */ + AdditionalData(const double omega = 1., + const int overlap = 0, + const bool fix_diagonal = false, + const double min_diagonal = 0., + const int n_sweeps = 1); + + /** + * @brief Relaxation damping factor. + * + */ + double omega; + /** + * @brief Overlap between processor local matrices. + * + * The amount of overlap between the processor local matrix blocks + * used in the underlying additive Schwarz method. + * + * The default will yield a block Jacobi preconditioner with each + * processor forming its own local block. + */ + int overlap; + /** + * @brief Whether or not to enlarge entries below a threshold + * + * If the matrix diagonal contains entries close to zero, their + * inversion will be close to division by zero. This can be corrected by + * setting this value to true. This will result in additional + * computational work during initialization. + */ + bool fix_diagonal; + /** + * @brief Threshold below which entries will be fixed. + * + * If the threshold is zero (default) only entries which are exactly + * zero will be replaced will small nonzero values. + * + */ + double min_diagonal; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + + + /** + * @brief Constructor. + * + */ + PreconditionSSOR(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * The class for the Chebyshev preconditioner within Ifpack2. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionChebyshev + : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * @param degree Degree of the Chebyshev polynomial. + * @param max_eigenvalue Upper bound for the maximum eigenvalue of the matrix. + * @param min_eigenvalue Lower bound for the minimum eigenvalue of the matrix. + * @param eigenvalue_ratio Ratio between maximum and minimum eigenvalue. + * @param min_diagonal Threshold for increasing diagonal entries. + * @param nonzero_starting Do not zero starting entries of solution. + */ + AdditionalData(const int degree = 1, + const double max_eigenvalue = 10., + const double min_eigenvalue = 1., + const double eigenvalue_ratio = 30., + const double min_diagonal = 1e-12, + const bool nonzero_starting = false); + /** + * @brief Degree of the Chebyshev polynomial. + * + * The degree directly corresponds to the number of matrix-vector + * products that have to be performed during a single application of the + * vmult() and Tvmult() operations. + * + */ + int degree; + /* + * @brief Upper bound for the maximum eigenvalue of the matrix. + * + * This needs to be set properly for the appropriate performance of this + * preconditioner. + */ + double max_eigenvalue; + /* + * @brief Lower bound for the minimum eigenvalue of the matrix. + * + */ + double min_eigenvalue; + /** + * @brief Estimated ratio between maximum and mimimum eigenvalue. + * + */ + double eigenvalue_ratio; + /** + * @brief Threshold below which entries will be fixed. + * + * If the threshold is zero (default) only entries which are exactly + * zero will be replaced will small nonzero values. + * + */ + double min_diagonal; + /** + * @brief Do not zero starting entries of solution vector. + * + * The default (false) zeroes out the entries of dst during vmult() and + * Tvmult() which is the recommended setting. + * + * However, in some situations (e.g. high-frequency error smoothing) it + * can be useful to append previous data to the Chebyshev corrections. + * The user should really know what they are doing when setting this + * flag to true. + * + */ + bool nonzero_starting; + }; + + /** + * @brief Constructor. + * + */ + PreconditionChebyshev(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * @brief The ILU/ILU(K)/RILU(K) preconditioner. + * + * The class for the modified incomplete LU factorization preconditioner + * RILUK within Ifpack2. This preconditioner works both in serial and + * parallel, depending on the matrix it is based on. + * + * In general, an incomplete factorization only has values on the sparsity + * pattern of the matrix, but one can gradually increase ilu_fill + * to eventually obtain a full LU factorization. + * + * For parallel applications this preconditioner implements an additive + * Schwarz preconditioner with RILUK as local smoothing on each processor. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionILU : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Before factorization each diagonal entry will be modified by the + * following fomula \f[a_{ii}^{new} = \alpha\sign(a_{ii})+\beta a_{ii} + * \f] with \f[\alpha\f] given by ilu_atol and \f[\beta\f] given by + * ilu_rtol. + * + * @param ilu_fill Amount of additional fill-in. + * @param ilu_atol Constant to be to each diagonal entry before factorization. + * @param ilu_rtol Factor to scale all diagonal entries by before factorization. + * @param overlapOverlap between processor local matrices. + */ + AdditionalData(const int ilu_fill = 0, + const double ilu_atol = 0., + const double ilu_rtol = 1., + const int overlap = 0); + + /** + * @brief Amount of additional fill-in. + * + * Level-of-fill to increase the sparsity pattern of the preconditioner. + * A large enough value will result in a complete LU factorization. + * + */ + int ilu_fill; + /** + * @brief Constant to be added to each diagonal entry before factorization. + * + */ + double ilu_atol; + /** + * @brief Factor to scale all diagonal entries by before factorization. + * + */ + double ilu_rtol; + /** + * @brief Overlap between processor local matrices. + * + * The amount of overlap between the processor local matrix blocks + * used in the underlying additive Schwarz method. + * + * The default will yield a block Jacobi preconditioner with each + * processor forming its own local block. + */ + int overlap; + }; + + /** + * @brief Constructor. + * + */ + PreconditionILU(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * The class for the thresholded incomplete LU (ILUT) preconditioner within + * Ifpack2. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionILUT : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Before factorization each diagonal entry will be modified by the + * following fomula \f[a_{ii}^{new} = \alpha\sign(a_{ii})+\beta a_{ii} + * \f] with \f[\alpha\f] given by ilut_atol and \f[\beta\f] given by + * ilut_rtol. + * + * @param ilut_drop Threshold for dropping entries. + * @param ilut_fill Amount of additional fill-in. + * @param ilut_atol Constant to be to each diagonal entry before factorization. + * @param ilut_rtol Factor to scale all diagonal entries by before factorization. + * @param overlapOverlap between processor local matrices. + */ + AdditionalData(const double ilut_drop = 0., + const double ilut_fill = 0., + const double ilut_atol = 0., + const double ilut_rtol = 1., + const int overlap = 0); + + /** + * @brief Threshold for dropping entries. + * + * Together with ilut_fill this controls the amount of fill-in + * and the actual values to be used or dropped. + */ + double ilut_drop; + /** + * @brief Amount of additional fill-in. + * + * Level-of-fill to increase the sparsity pattern of the preconditioner. + * A large enough value will result in a complete LU factorization. + * Note, however, that this will drastically increase the memory + * requirement, especially for 3d simulations. + * + */ + double ilut_fill; + /** + * @brief Constant to be added to each diagonal entry before factorization. + * + */ + double ilut_atol; + /** + * @brief Factor to scale all diagonal entries by before factorization. + * + */ + double ilut_rtol; + /** + * @brief Overlap between processor local matrices. + * + * The amount of overlap between the processor local matrix blocks + * used in the underlying additive Schwarz method. + * + * The default will yield a block Jacobi preconditioner with each + * processor forming its own local block. + */ + int overlap; + }; + + + + /** + * @brief Constructor. + * + */ + PreconditionILUT(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * The class for the BlockJacobi preconditioner within Ifpack2. + * + * @note This preconditioner always uses linear partitioning. + * There are other partitioners available that need additional user provided + * data, so we suggest using PreconditionIfpack with a custom + * Teuchos::ParameterList. Details on the parameters can be found in the + * User's Guide of Ifpack2. + * + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionBlockJacobi + : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * The defaults yield a textbook Jacobi preconditioner. + * + * @param n_local_parts The number of blocks per processor. + * @param omega The damping factor for the relaxation. + * @param block_overlap Amount of overlap between blocks + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult + */ + AdditionalData(const int n_local_parts = 1, + const double omega = 1., + const int block_overlap = 0, + const int n_sweeps = 1); + + /** + * @brief Number of blocks per processor. + * + * The default of 1 results in a single block per processor + * which corresponds to Preconditioner from PreconditionJacobi + */ + int n_local_parts; + /** + * @brief Relaxation damping factor. + * + */ + double omega; + /** + * @brief Amount of overlap between blocks + * + */ + int block_overlap; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + + + /** + * @brief Constructor. + * + */ + PreconditionBlockJacobi(); + + /** + * @brief Compute the preconditioner based on the given matrix and parameters. + * + * @param A The matrix to base the preconditioner on. + * @param additional_data The set of parameters to tune the preconditioner. + */ + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * @brief The class for the Block SOR preconditioner within Ifpack2. + * + * If the code is executed MPI-parallel the individual processors will be + * connected with an additive Schwarz method. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionBlockSOR + : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * + * @param n_local_parts The number of blocks per processor. + * @param omega The damping factor for the relaxation. + * @param overlap Overlap between processors + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult + */ + AdditionalData(const int n_local_parts = 1, + const double omega = 1, + const int overlap = 0, + const int n_sweeps = 1); + + /** + * @brief Number of blocks per processor. + * + * The default of 1 results in a single block per processor + * which corresponds to Preconditioner from PreconditionJacobi + */ + int n_local_parts; + /** + * @brief Relaxation damping factor. + * + */ + double omega; + /** + * @brief Overlap between processor local matrices. + * + * The amount of overlap between the processor local matrix blocks + * used in the underlying additive Schwarz method. + * + * The default will yield a block Jacobi preconditioner with each + * processor forming its own local block. + */ + int overlap; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + + PreconditionBlockSOR(); + + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; + + + + /** + * @brief The class for the Block SSOR preconditioner within Ifpack2. + * + * If the code is executed MPI-parallel the individual processors will be + * connected with an additive Schwarz method. + * + * @ingroup TpetraWrappers + * @ingroup Preconditioners + */ + template + class PreconditionBlockSSOR + : public PreconditionIfpackBase + { + public: + /** + * @brief The set of additional parameters to tune the preconditioner. + * + */ + struct AdditionalData + { + /** + * @brief Constructor. + * + * Set the parameters to be used in the preconditioner. + * + * @param n_local_parts The number of blocks per processor. + * @param omega The damping factor for the relaxation. + * @param overlap Overlap between processors + * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult + */ + AdditionalData(const int n_local_parts = 1, + const double omega = 1, + const int overlap = 1, + const int n_sweeps = 1); + /** + * @brief Number of blocks per processor. + * + * The default of 1 results in a single block per processor + * which corresponds to Preconditioner from PreconditionJacobi + */ + int n_local_parts; + /** + * @brief Relaxation damping factor. + * + */ + double omega; + /** + * @brief Overlap between processor local matrices. + * + * The amount of overlap between the processor local matrix blocks + * used in the underlying additive Schwarz method. + * + * The default will yield a block Jacobi preconditioner with each + * processor forming its own local block. + */ + int overlap; + /** + * @brief Set how often the preconditioner should be applied during vmult() or Tvmult(). + * + */ + int n_sweeps; + }; + + PreconditionBlockSSOR(); + + void + initialize(const SparseMatrix &A, + const AdditionalData &additional_data = AdditionalData()); + }; +# endif // DEAL_II_TRILINOS_WITH_IFPACK2 + + } // namespace TpetraWrappers +} // namespace LinearAlgebra + + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif diff --git a/include/deal.II/lac/trilinos_tpetra_precondition.templates.h b/include/deal.II/lac/trilinos_tpetra_precondition.templates.h new file mode 100644 index 0000000000..554f6236a3 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_precondition.templates.h @@ -0,0 +1,805 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +#ifndef dealii_trilinos_tpetra_precondition_templates_h +#define dealii_trilinos_tpetra_precondition_templates_h + + +#include + +#include "deal.II/base/exceptions.h" +#include + +#include "deal.II/lac/trilinos_tpetra_types.h" + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA +# ifdef DEAL_II_TRILINOS_WITH_IFPACK2 +# include +# include + +# include +# include +# include +# include +# include +# include +# include +# include +# include + +DEAL_II_NAMESPACE_OPEN + + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + /* ---------------------- PreconditionBase ------------------------ */ + template + void + PreconditionBase::clear() + { + preconditioner.reset(); + } + + + + template + const TpetraTypes::LinearOperator & + PreconditionBase::trilinos_operator() const + { + return *preconditioner; + } + + + + template + Teuchos::RCP> + PreconditionBase::trilinos_rcp() const + { + return preconditioner; + } + + + + template + IndexSet + PreconditionBase::locally_owned_domain_indices() const + { + return IndexSet(preconditioner->getDomainMap()); + } + + + + template + IndexSet + PreconditionBase::locally_owned_range_indices() const + { + return IndexSet(preconditioner->getRangeMap()); + } + + + + template + inline void + PreconditionBase::vmult( + Vector &dst, + const Vector &src) const + { + Assert(dst.trilinos_vector().getMap()->isSameAs( + *(preconditioner->getRangeMap())), + ExcNonMatchingMaps("dst")); + Assert(src.trilinos_rcp()->getMap()->isSameAs( + *(preconditioner->getDomainMap())), + ExcNonMatchingMaps("src")); + + preconditioner->apply(src.trilinos_vector(), dst.trilinos_vector()); + } + + + + template + inline void + PreconditionBase::Tvmult( + Vector &dst, + const Vector &src) const + { + Assert(preconditioner->hasTransposeApply(), ExcTransposeNotSupported()); + + Assert(dst.trilinos_rcp()->getMap()->isSameAs( + *(preconditioner->getDomainMap())), + ExcNonMatchingMaps("dst")); + + Assert(src.trilinos_rcp()->getMap()->isSameAs( + *(preconditioner->getRangeMap())), + ExcNonMatchingMaps("src")); + + preconditioner->apply(src.trilinos_vector(), + dst.trilinos_vector(), + Teuchos::TRANS); + } + + + + template + inline void + PreconditionBase::vmult( + dealii::Vector &dst, + dealii::Vector &src) const + { + AssertDimension(dst.size(), + preconditioner->getRangeMap()->getGlobalNumElements()); + AssertDimension(src.size(), + preconditioner->getDomainMap()->getGlobalNumElements()); + + // Construct Tpetra Vectors as views to the actual dealii::Vector data + // Since all of Tpetra is based on Kokkos::DualView this is a bit of a + // process + + // 1. Make host view with underlying dealii::Vector data. + // (This is always a n x 1 matrix since multivector expects that.) + TpetraTypes::HostViewTypeUnmanaged host_src(src.data(), + src.size(), + 1); + TpetraTypes::HostViewTypeUnmanaged host_dst(dst.data(), + dst.size(), + 1); + + // 2. Create device mirror (if Device == Host, no alloc ) + auto dev_src = Kokkos::create_mirror_view_and_copy( + typename MemorySpace::kokkos_space::execution_space(), host_src); + + auto dev_dst = Kokkos::create_mirror_view_and_copy( + typename MemorySpace::kokkos_space::execution_space(), host_dst); + + // 3. Create dual views from the previous ones + TpetraTypes::DualViewTypeUnmanaged dual_src( + dev_src, host_src); + TpetraTypes::DualViewTypeUnmanaged dual_dst( + dev_dst, host_dst); + + // 4. Create Tpetra Vector from + TpetraTypes::VectorType tpetra_dst( + preconditioner->getRangeMap(), dual_dst); + TpetraTypes::VectorType tpetra_src( + preconditioner->getDomainMap(), dual_src); + + // 5. Finally, apply our preconditioner + preconditioner->apply(tpetra_src, tpetra_dst); + } + + + + template + inline void + PreconditionBase::Tvmult( + dealii::Vector &dst, + dealii::Vector &src) const + { + AssertDimension(dst.size(), + preconditioner->getDomainMap()->getGlobalNumElements()); + AssertDimension(src.size(), + preconditioner->getRangeMap()->getGlobalNumElements()); + + Assert(preconditioner->hasTransposeApply(), ExcTransposeNotSupported()); + // Construct Tpetra Vectors as views to the actual dealii::Vector data + // Since all of Tpetra is based on Kokkos::DualView this is a bit of a + // process + + // 1. Make host view with underlying dealii::Vector data. + // (This is always a n x 1 matrix since multivector expects that.) + TpetraTypes::HostViewTypeUnmanaged host_src(src.data(), + src.size(), + 1); + TpetraTypes::HostViewTypeUnmanaged host_dst(dst.data(), + dst.size(), + 1); + + // 2. Create device mirror (if Device == Host, no alloc ) + auto dev_src = Kokkos::create_mirror_view_and_copy( + typename MemorySpace::kokkos_space::execution_space(), host_src); + + auto dev_dst = Kokkos::create_mirror_view_and_copy( + typename MemorySpace::kokkos_space::execution_space(), host_dst); + + // 3. Create dual views from the previous ones + TpetraTypes::DualViewTypeUnmanaged dual_src( + dev_src, host_src); + TpetraTypes::DualViewTypeUnmanaged dual_dst( + dev_dst, host_dst); + + // 4. Create Tpetra Vector from + TpetraTypes::VectorType tpetra_dst( + preconditioner->getRangeMap(), dual_dst); + TpetraTypes::VectorType tpetra_src( + preconditioner->getDomainMap(), dual_src); + + // 5. Finally, apply our preconditioner + preconditioner->apply(tpetra_src, tpetra_dst, Teuchos::TRANS); + } + + + + /* ---------------------- PreconditionIdentity ------------------------ */ + template + void + PreconditionIdentity::initialize( + const SparseMatrix &A) + { + Teuchos::RCP>> + ifpack2Preconditioner; + ifpack2Preconditioner = + rcp(new Ifpack2::IdentitySolver< + TpetraTypes::RowMatrixType>(A.trilinos_rcp())); + + ifpack2Preconditioner->initialize(); + ifpack2Preconditioner->compute(); + this->preconditioner = ifpack2Preconditioner; + } + + + + /* ---------------------- PreconditionIfpackBase ------------------------ */ + template + PreconditionIfpackBase::PreconditionIfpackBase( + const std::string &preconditioner_type) + : preconditioner_type(preconditioner_type) + { + Ifpack2::Factory factory; + bool supported = + factory.isSupported>( + preconditioner_type); + AssertThrow(supported, + ExcTrilinosIpack2PreconditionerUnsupported( + preconditioner_type)); + } + + + + template + void + PreconditionIfpackBase::initialize( + const SparseMatrix &A) + { + Ifpack2::Factory factory; + Teuchos::RCP> + ifpack2Preconditioner; + ifpack2Preconditioner = + factory.create(preconditioner_type, A.trilinos_rcp()); + ifpack2Preconditioner->setParameters(this->parameter_list); + ifpack2Preconditioner->initialize(); + ifpack2Preconditioner->compute(); + this->preconditioner = ifpack2Preconditioner; + } + + + + /* ---------------------- PreconditionIfpack ------------------------ */ + template + PreconditionIfpack::PreconditionIfpack( + const std::string &preconditioner_type) + : PreconditionIfpackBase(preconditioner_type) + {} + + + + template + void + PreconditionIfpack::set_parameter_list( + Teuchos::ParameterList ¶meter_list) + { + this->parameter_list.setParameters(parameter_list); + } + + + + /* ---------------------- PreconditionJacobi ------------------------ */ + template + PreconditionJacobi::AdditionalData::AdditionalData( + const double omega, + const bool fix_diagonal, + const double min_diagonal, + const int n_sweeps) + : omega(omega) + , fix_diagonal(fix_diagonal) + , min_diagonal(min_diagonal) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionJacobi::PreconditionJacobi() + : PreconditionIfpackBase("RELAXATION") + {} + + + + template + void + PreconditionJacobi::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("relaxation: type", "Jacobi"); + this->parameter_list.set("relaxation: damping factor", ad.omega); + this->parameter_list.set("relaxation: sweeps", ad.n_sweeps); + this->parameter_list.set("relaxation: fix tiny diagonal entries", + ad.fix_diagonal); + this->parameter_list.set("relaxation: min diagonal value", + ad.min_diagonal); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionL1Jacobi ------------------------ */ + template + PreconditionL1Jacobi::AdditionalData::AdditionalData( + const double omega, + const bool fix_diagonal, + const double min_diagonal, + const int n_sweeps) + : omega(omega) + , fix_diagonal(fix_diagonal) + , min_diagonal(min_diagonal) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionL1Jacobi::PreconditionL1Jacobi() + : PreconditionIfpackBase("RELAXATION") + {} + + + + template + void + PreconditionL1Jacobi::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("relaxation: type", "Jacobi"); + this->parameter_list.set("relaxation: use l1", true); + this->parameter_list.set("relaxation: damping factor", ad.omega); + this->parameter_list.set("relaxation: sweeps", ad.n_sweeps); + this->parameter_list.set("relaxation: fix tiny diagonal entries", + ad.fix_diagonal); + this->parameter_list.set("relaxation: min diagonal value", + ad.min_diagonal); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionL1GaussSeidel -------------------- */ + template + PreconditionL1GaussSeidel::AdditionalData:: + AdditionalData(const double omega, + const double eta, + const bool fix_diagonal, + const double min_diagonal, + const int n_sweeps) + : omega(omega) + , eta(eta) + , fix_diagonal(fix_diagonal) + , min_diagonal(min_diagonal) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionL1GaussSeidel::PreconditionL1GaussSeidel() + : PreconditionIfpackBase("RELAXATION") + {} + + + + template + void + PreconditionL1GaussSeidel::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("relaxation: type", "Gauss-Seidel"); + this->parameter_list.set("relaxation: use l1", true); + this->parameter_list.set("relaxation: damping factor", ad.omega); + this->parameter_list.set("relaxation: l1 eta", ad.eta); + this->parameter_list.set("relaxation: sweeps", ad.n_sweeps); + this->parameter_list.set("relaxation: fix tiny diagonal entries", + ad.fix_diagonal); + this->parameter_list.set("relaxation: min diagonal value", + ad.min_diagonal); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionSOR ------------------------ */ + template + PreconditionSOR::AdditionalData::AdditionalData( + const double omega, + const int overlap, + const bool fix_diagonal, + const double min_diagonal, + const int n_sweeps) + : omega(omega) + , overlap(overlap) + , fix_diagonal(fix_diagonal) + , min_diagonal(min_diagonal) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionSOR::PreconditionSOR() + : PreconditionIfpackBase("SCHWARZ") + {} + + + + template + void + PreconditionSOR::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("schwarz: subdomain solver name", "RELAXATION"); + { + Teuchos::ParameterList &relaxParams = + this->parameter_list.sublist("schwarz: subdomain solver parameters", + false); + relaxParams.set("relaxation: type", "Symmetric Gauss-Seidel"); + relaxParams.set("relaxation: damping factor", ad.omega); + relaxParams.set("relaxation: sweeps", ad.n_sweeps); + relaxParams.set("relaxation: fix tiny diagonal entries", + ad.fix_diagonal); + relaxParams.set("relaxation: min diagonal value", ad.min_diagonal); + } + this->parameter_list.set("schwarz: combine mode", "ADD"); + this->parameter_list.set("schwarz: overlap level", ad.overlap); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionSSOR ----------------------*/ + template + PreconditionSSOR::AdditionalData::AdditionalData( + const double omega, + const int overlap, + const bool fix_diagonal, + const double min_diagonal, + const int n_sweeps) + : omega(omega) + , overlap(overlap) + , fix_diagonal(fix_diagonal) + , min_diagonal(min_diagonal) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionSSOR::PreconditionSSOR() + : PreconditionIfpackBase("SCHWARZ") + {} + + + + template + void + PreconditionSSOR::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("schwarz: subdomain solver name", "RELAXATION"); + { + Teuchos::ParameterList &relaxParams = + this->parameter_list.sublist("schwarz: subdomain solver parameters", + false); + relaxParams.set("relaxation: type", "Symmetric Gauss-Seidel"); + relaxParams.set("relaxation: damping factor", ad.omega); + relaxParams.set("relaxation: sweeps", ad.n_sweeps); + relaxParams.set("relaxation: fix tiny diagonal entries", + ad.fix_diagonal); + relaxParams.set("relaxation: min diagonal value", ad.min_diagonal); + } + this->parameter_list.set("schwarz: combine mode", "ADD"); + this->parameter_list.set("schwarz: overlap level", ad.overlap); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionBlockJacobi ----------------------*/ + template + PreconditionBlockJacobi::AdditionalData:: + AdditionalData(const int n_local_parts, + const double omega, + const int block_overlap, + const int n_sweeps) + : n_local_parts(n_local_parts) + , omega(omega) + , block_overlap(block_overlap) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionBlockJacobi::PreconditionBlockJacobi() + : PreconditionIfpackBase("BLOCK_RELAXATION") + {} + + + + template + void + PreconditionBlockJacobi::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("relaxation: type", "Jacobi"); + this->parameter_list.set("relaxation: damping factor", ad.omega); + this->parameter_list.set("relaxation: sweeps", ad.n_sweeps); + this->parameter_list.set("partitioner: local parts", ad.n_local_parts); + this->parameter_list.set("partitioner: overlap", ad.block_overlap); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionBlockSOR ----------------------*/ + template + PreconditionBlockSOR::AdditionalData::AdditionalData( + const int n_local_parts, + const double omega, + const int overlap, + const int n_sweeps) + : n_local_parts(n_local_parts) + , omega(omega) + , overlap(overlap) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionBlockSOR::PreconditionBlockSOR() + : PreconditionIfpackBase("SCHWARZ") + {} + + + + template + void + PreconditionBlockSOR::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("schwarz: subdomain solver name", + "BLOCK_RELAXATION"); + { + Teuchos::ParameterList &relaxParams = + this->parameter_list.sublist("schwarz: subdomain solver parameters", + false); + relaxParams.set("relaxation: type", "Gauss-Seidel"); + relaxParams.set("relaxation: damping factor", ad.omega); + relaxParams.set("relaxation: sweeps", ad.n_sweeps); + relaxParams.set("partitioner: local parts", ad.n_local_parts); + } + this->parameter_list.set("schwarz: combine mode", "ADD"); + this->parameter_list.set("schwarz: overlap level", ad.overlap); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionBlockSSOR ----------------------*/ + template + PreconditionBlockSSOR::AdditionalData::AdditionalData( + const int n_local_parts, + const double omega, + const int overlap, + const int n_sweeps) + : n_local_parts(n_local_parts) + , omega(omega) + , overlap(overlap) + , n_sweeps(n_sweeps) + {} + + + + template + PreconditionBlockSSOR::PreconditionBlockSSOR() + : PreconditionIfpackBase("SCHWARZ") + {} + + + + template + void + PreconditionBlockSSOR::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("schwarz: subdomain solver name", + "BLOCK_RELAXATION"); + { + Teuchos::ParameterList &relaxParams = + this->parameter_list.sublist("schwarz: subdomain solver parameters", + false); + relaxParams.set("relaxation: type", "Symmetric Gauss-Seidel"); + relaxParams.set("relaxation: damping factor", ad.omega); + relaxParams.set("relaxation: sweeps", ad.n_sweeps); + relaxParams.set("partitioner: local parts", ad.n_local_parts); + } + this->parameter_list.set("schwarz: combine mode", "ADD"); + this->parameter_list.set("schwarz: overlap level", ad.overlap); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionChebyshev ----------------------*/ + template + PreconditionChebyshev::AdditionalData::AdditionalData( + const int degree, + const double max_eigenvalue, + const double eigenvalue_ratio, + const double min_eigenvalue, + const double min_diagonal, + const bool nonzero_starting) + : degree(degree) + , max_eigenvalue(max_eigenvalue) + , min_eigenvalue(min_eigenvalue) + , eigenvalue_ratio(eigenvalue_ratio) + , min_diagonal(min_diagonal) + , nonzero_starting(nonzero_starting) + {} + + + + template + PreconditionChebyshev::PreconditionChebyshev() + : PreconditionIfpackBase("CHEBYSHEV") + {} + + + + template + void + PreconditionChebyshev::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("chebyshev: degree", ad.degree); + this->parameter_list.set("chebyshev: max eigenvalue", ad.max_eigenvalue); + this->parameter_list.set("chebyshev: min eigenvalue", ad.min_eigenvalue); + this->parameter_list.set("chebyshev: ratio eigenvalue", + ad.eigenvalue_ratio); + this->parameter_list.set("chebyshev: min diagonal value", + ad.min_diagonal); + this->parameter_list.set("chebyshev: zero starting solution", + !ad.nonzero_starting); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionILU ----------------------*/ + template + PreconditionILU::AdditionalData::AdditionalData( + const int ilu_fill, + const double ilu_atol, + const double ilu_rtol, + const int overlap) + : ilu_fill(ilu_fill) + , ilu_atol(ilu_atol) + , ilu_rtol(ilu_rtol) + , overlap(overlap) + {} + + + + template + PreconditionILU::PreconditionILU() + : PreconditionIfpackBase("SCHWARZ") + {} + + + + template + void + PreconditionILU::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("schwarz: subdomain solver name", "RILUK"); + { + Teuchos::ParameterList &rilukParams = + this->parameter_list.sublist("schwarz: subdomain solver parameters", + false); + rilukParams.set("fact: iluk level-of-fill", ad.ilu_fill); + rilukParams.set("fact: absolute threshold", ad.ilu_atol); + rilukParams.set("fact: relative threshold", ad.ilu_rtol); + } + this->parameter_list.set("schwarz: combine mode", "ADD"); + this->parameter_list.set("schwarz: overlap level", ad.overlap); + PreconditionIfpackBase::initialize(A); + } + + + + /* ---------------------- PreconditionILUT ----------------------*/ + template + PreconditionILUT::AdditionalData::AdditionalData( + const double ilut_drop, + const double ilut_fill, + const double ilut_atol, + const double ilut_rtol, + const int overlap) + : ilut_drop(ilut_drop) + , ilut_fill(ilut_fill) + , ilut_atol(ilut_atol) + , ilut_rtol(ilut_rtol) + , overlap(overlap) + {} + + + + template + PreconditionILUT::PreconditionILUT() + : PreconditionIfpackBase("SCHWARZ") + {} + + + + template + void + PreconditionILUT::initialize( + const SparseMatrix &A, + const AdditionalData &ad) + { + this->parameter_list.set("schwarz: subdomain solver name", "ILUT"); + { + Teuchos::ParameterList &ilutParams = + this->parameter_list.sublist("schwarz: subdomain solver parameters", + false); + ilutParams.set("fact: drop tolerance", ad.ilut_drop); + ilutParams.set("fact: ilut level-of-fill", ad.ilut_fill); + ilutParams.set("fact: absolute threshold", ad.ilut_atol); + ilutParams.set("fact: relative threshold", ad.ilut_rtol); + } + this->parameter_list.set("schwarz: combine mode", "ADD"); + this->parameter_list.set("schwarz: overlap level", ad.overlap); + PreconditionIfpackBase::initialize(A); + } + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_TRILINOS_WITH_IFPACK2 +#endif // DEAL_II_TRILINOS_WITH_TPETRA + +#endif diff --git a/include/deal.II/lac/trilinos_tpetra_types.h b/include/deal.II/lac/trilinos_tpetra_types.h index 5dda027e86..bf38dc7784 100644 --- a/include/deal.II/lac/trilinos_tpetra_types.h +++ b/include/deal.II/lac/trilinos_tpetra_types.h @@ -21,10 +21,13 @@ #ifdef DEAL_II_TRILINOS_WITH_TPETRA +# include # include + // Forward declarations # ifndef DOXYGEN +# include # include # include # include @@ -32,8 +35,8 @@ # include # include # include +# include # include - # endif DEAL_II_NAMESPACE_OPEN @@ -49,6 +52,8 @@ namespace LinearAlgebra * local ordinate (processor local indices). */ using LO = int; + + /** * global ordinate (global indices). */ @@ -72,7 +77,6 @@ namespace LinearAlgebra typename MemorySpace::kokkos_space>; # endif - /** * Communication between processors. */ @@ -128,6 +132,55 @@ namespace LinearAlgebra template using MatrixType = Tpetra::CrsMatrix>; + + + /** + * Tpetra type for a parallel distributed row matrix. + */ + template + using RowMatrixType = + Tpetra::RowMatrix>; + + + /** + * @brief Typedef for the Kokkos::View type of unmanaged host memory. + * + * Unmanaged means that the view is not allocating data itself, but only + * uses an external pointer. Consequently, the memory is not freed upon + * destruction of the view. + * This is needed for shallow copies of deal.II LA structures + * to Trilinos LA structures. + * + */ + template + using HostViewTypeUnmanaged = Kokkos::View; + + /** + * @brief Typedef for the Kokkos::View type of unmanaged memory. + * + * Unmanaged means that the view is not allocating data itself, but only + * uses an external pointer. Consequently, the memory is not freed upon + * destruction of the view. + * This is needed for shallow copies of deal.II LA structures + * to Trilinos LA structures. + */ + template + using DualViewTypeUnmanaged = + Kokkos::DualView; + + /** + * Type for a Trilinos preconditioner from the Ifpack2 package. + */ + template + using Ifpack2PreconType = + Ifpack2::Preconditioner>; + } // namespace TpetraTypes } // namespace TpetraWrappers } // namespace LinearAlgebra diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index 5dec885430..f77a42be6a 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -145,6 +145,7 @@ if(DEAL_II_WITH_TRILINOS) trilinos_tpetra_block_vector.cc trilinos_tpetra_block_sparse_matrix.cc trilinos_tpetra_communication_pattern.cc + trilinos_tpetra_precondition.cc trilinos_tpetra_solver_direct.cc trilinos_tpetra_sparse_matrix.cc trilinos_tpetra_sparsity_pattern.cc diff --git a/source/lac/trilinos_tpetra_precondition.cc b/source/lac/trilinos_tpetra_precondition.cc new file mode 100644 index 0000000000..21f7ee2562 --- /dev/null +++ b/source/lac/trilinos_tpetra_precondition.cc @@ -0,0 +1,109 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +#include "deal.II/lac/trilinos_tpetra_precondition.templates.h" + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA +# ifdef DEAL_II_TRILINOS_WITH_IFPACK2 + +DEAL_II_NAMESPACE_OPEN + +# ifndef DOXYGEN +// explicit instantiations +namespace LinearAlgebra +{ + namespace TpetraWrappers + { +# ifdef HAVE_TPETRA_INST_FLOAT + template class PreconditionBase; + template class PreconditionIdentity; + template class PreconditionIfpackBase; + template class PreconditionIfpack; + template class PreconditionJacobi; + template class PreconditionL1Jacobi; + template class PreconditionL1GaussSeidel; + template class PreconditionSOR; + template class PreconditionSSOR; + template class PreconditionBlockJacobi; + template class PreconditionBlockSOR; + template class PreconditionBlockSSOR; + template class PreconditionChebyshev; + template class PreconditionILU; + template class PreconditionILUT; +# endif + +# ifdef HAVE_TPETRA_INST_DOUBLE + template class PreconditionBase; + template class PreconditionIdentity; + template class PreconditionIfpackBase; + template class PreconditionIfpack; + template class PreconditionJacobi; + template class PreconditionL1Jacobi; + template class PreconditionL1GaussSeidel; + template class PreconditionSOR; + template class PreconditionSSOR; + template class PreconditionBlockJacobi; + template class PreconditionBlockSOR; + template class PreconditionBlockSSOR; + template class PreconditionChebyshev; + template class PreconditionILU; + template class PreconditionILUT; +# endif +# ifdef DEAL_II_WITH_COMPLEX_VALUES +# ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT + template class PreconditionBase>; + template class PreconditionIdentity>; + template class PreconditionIfpackBase>; + template class PreconditionIfpack>; + template class PreconditionJacobi>; + template class PreconditionL1Jacobi>; + template class PreconditionL1GaussSeidel>; + template class PreconditionSOR>; + template class PreconditionSSOR>; + template class PreconditionBlockJacobi>; + template class PreconditionBlockSOR>; + template class PreconditionBlockSSOR>; + template class PreconditionChebyshev>; + template class PreconditionILU>; + template class PreconditionILUT>; +# endif +# ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE + template class PreconditionBase>; + template class PreconditionIdentity>; + template class PreconditionIfpackBase>; + template class PreconditionIfpack>; + template class PreconditionJacobi>; + template class PreconditionL1Jacobi>; + template class PreconditionL1GaussSeidel>; + template class PreconditionSOR>; + template class PreconditionSSOR>; + template class PreconditionBlockJacobi>; + template class PreconditionBlockSOR>; + template class PreconditionBlockSSOR>; + template class PreconditionChebyshev>; + template class PreconditionILU>; + template class PreconditionILUT>; +# endif +# endif + +# endif + + } // namespace TpetraWrappers + +} // namespace LinearAlgebra +DEAL_II_NAMESPACE_CLOSE + +# endif // DEAL_II_TRILINOS_WITH_IFPACK2 +#endif // DEAL_II_TRILINOS_WITH_TPETRA diff --git a/tests/trilinos_tpetra/add_matrices_06.cc b/tests/trilinos_tpetra/add_matrices_06.cc new file mode 100644 index 0000000000..180281b000 --- /dev/null +++ b/tests/trilinos_tpetra/add_matrices_06.cc @@ -0,0 +1,127 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2004 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// check that SparseMatrix::add(SparseMatrix) and +// SparseMatrix::copy_from(SparseMatrix) do not destroy memory when called +// from matrices with the same sparsity pattern and thus, structures depending +// on these matrices such as PreconditionJacobi do not get destroyed + +#include + +#include +#include + +#include + +#include "../tests.h" + + +void +test(LinearAlgebra::TpetraWrappers::SparseMatrix &m) +{ + LinearAlgebra::TpetraWrappers::SparseMatrix m2(m.m(), m.n(), 3U), + m3(m.m(), m.n(), 3U); + + // first set a few entries one-by-one + for (unsigned int i = 0; i < m.m(); ++i) + for (unsigned int j = 0; j < m.n(); ++j) + if ((i + 2 * j + 1) % 3 == 0 || i == j) + { + m.set(i, j, i * j * .5 + .5); + m2.set(i, j, 1.); + m3.set(i, j, -0.1); + } + + m.compress(VectorOperation::insert); + m2.compress(VectorOperation::insert); + m3.compress(VectorOperation::insert); + + deallog << "Matrix nonzeros: " << m.n_nonzero_elements() << ' ' + << m2.n_nonzero_elements() << ' ' << m3.n_nonzero_elements() + << std::endl; + + m.copy_from(m2); + m.add(0.12, m3); + + Vector src(m.m()), dst(m.m()); + src = 1.; + + LinearAlgebra::TpetraWrappers::PreconditionJacobi prec; + prec.initialize(m); + + m.vmult(dst, src); + deallog << "Vector norm: " << dst.l2_norm() << ' '; + + prec.vmult(dst, src); + deallog << dst.l2_norm() << std::endl; + + m.copy_from(m2); + m.add(0.06, m3); + + m.vmult(dst, src); + deallog << "Vector norm: " << dst.l2_norm() << ' '; + + prec.vmult(dst, src); + deallog << dst.l2_norm() << std::endl; + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + try + { + { + LinearAlgebra::TpetraWrappers::SparseMatrix m(5U, 5U, 3U); + + test(m); + } + } + catch (const std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/add_matrices_06.output b/tests/trilinos_tpetra/add_matrices_06.output new file mode 100644 index 0000000000..f1e4081029 --- /dev/null +++ b/tests/trilinos_tpetra/add_matrices_06.output @@ -0,0 +1,5 @@ + +DEAL::Matrix nonzeros: 13 13 13 +DEAL::Vector norm: 5.84509 2.26323 +DEAL::Vector norm: 5.88058 2.26323 +DEAL::OK diff --git a/tests/trilinos_tpetra/direct_solver_01.cc b/tests/trilinos_tpetra/direct_solver_01.cc index 04bf70831f..24f71ee43f 100644 --- a/tests/trilinos_tpetra/direct_solver_01.cc +++ b/tests/trilinos_tpetra/direct_solver_01.cc @@ -31,13 +31,12 @@ #include #include -#include "deal.II/lac/exceptions.h" #include #include #include #include #include -// #include +#include #include #include @@ -300,20 +299,16 @@ template void Step4::solve() { - // TODO: implement precondition SSOR Wrappers - // // Compute 'reference' solution with CG solver and SSOR preconditioner - // TrilinosWrappers::PreconditionSSOR preconditioner; - // preconditioner.initialize(system_matrix); + // Compute 'reference' solution with CG solver and SSOR preconditioner + LinearAlgebra::TpetraWrappers::PreconditionSSOR preconditioner; + preconditioner.initialize(system_matrix); LinearAlgebra::TpetraWrappers::Vector temp_solution(system_rhs); temp_solution = 0; SolverControl solver_control(1000, 1e-12); SolverCG> solver( solver_control); - solver.solve(system_matrix, - temp_solution, - system_rhs, - PreconditionIdentity()); + solver.solve(system_matrix, temp_solution, system_rhs, preconditioner); constraints.distribute(temp_solution); solution = temp_solution; @@ -323,10 +318,7 @@ Step4::solve() // do CG solve for new rhs temp_solution = 0; solution = 0; - solver.solve(system_matrix, - temp_solution, - system_rhs_two, - PreconditionIdentity()); + solver.solve(system_matrix, temp_solution, system_rhs_two, preconditioner); constraints.distribute(temp_solution); solution = temp_solution; diff --git a/tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output b/tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output index 1facc1d50a..dff97d268c 100644 --- a/tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output +++ b/tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output @@ -1,8 +1,8 @@ DEAL:cg::Starting value 21.8027 -DEAL:cg::Convergence step 109 value 0 +DEAL:cg::Convergence step 86 value 0 DEAL:cg::Starting value 0.0940512 -DEAL:cg::Convergence step 98 value 0 +DEAL:cg::Convergence step 77 value 0 DEAL:DirectKLU2::Starting value 0 DEAL:DirectKLU2::Convergence step 0 value 0 DEAL:DirectKLU2::Norm of error in direct solve 1: 0 diff --git a/tests/trilinos_tpetra/direct_solver_02.cc b/tests/trilinos_tpetra/direct_solver_02.cc index 5e01f25679..e384b2fe0b 100644 --- a/tests/trilinos_tpetra/direct_solver_02.cc +++ b/tests/trilinos_tpetra/direct_solver_02.cc @@ -35,7 +35,7 @@ #include #include #include -// #include +#include #include #include @@ -298,20 +298,16 @@ template void Step4::solve() { - // TODO: implement precondition SSOR Wrappers - // // Compute 'reference' solution with CG solver and SSOR preconditioner - // TrilinosWrappers::PreconditionSSOR preconditioner; - // preconditioner.initialize(system_matrix); + // Compute 'reference' solution with CG solver and SSOR preconditioner + LinearAlgebra::TpetraWrappers::PreconditionSSOR preconditioner; + preconditioner.initialize(system_matrix); LinearAlgebra::TpetraWrappers::Vector temp_solution(system_rhs); temp_solution = 0; SolverControl solver_control(1000, 1e-12); SolverCG> solver( solver_control); - solver.solve(system_matrix, - temp_solution, - system_rhs, - PreconditionIdentity()); + solver.solve(system_matrix, temp_solution, system_rhs, preconditioner); constraints.distribute(temp_solution); solution = temp_solution; @@ -321,10 +317,7 @@ Step4::solve() // do CG solve for new rhs temp_solution = 0; solution = 0; - solver.solve(system_matrix, - temp_solution, - system_rhs_two, - PreconditionIdentity()); + solver.solve(system_matrix, temp_solution, system_rhs_two, preconditioner); constraints.distribute(temp_solution); solution = temp_solution; diff --git a/tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output b/tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output index 2f8d30713d..12fce67a02 100644 --- a/tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output +++ b/tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output @@ -1,8 +1,8 @@ DEAL:cg::Starting value 21.8027 -DEAL:cg::Convergence step 109 value 0 +DEAL:cg::Convergence step 86 value 0 DEAL:cg::Starting value 0.0940512 -DEAL:cg::Convergence step 98 value 0 +DEAL:cg::Convergence step 77 value 0 DEAL:DirectKLU2::Starting value 0 DEAL:DirectKLU2::Convergence step 0 value 0 DEAL:DirectKLU2::Norm of error in direct solve 1: 0 diff --git a/tests/trilinos_tpetra/precondition.cc b/tests/trilinos_tpetra/precondition.cc new file mode 100644 index 0000000000..efbdcc81e3 --- /dev/null +++ b/tests/trilinos_tpetra/precondition.cc @@ -0,0 +1,584 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2013 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// solves a 2D Poisson equation for linear elements with the simple +// (Ifpack-based) Trilinos preconditioners + +#include "deal.II/base/exceptions.h" +#include "deal.II/base/logstream.h" +#include + +#include +#include + +#include +#include + +#include +#include + +#include "deal.II/lac/trilinos_tpetra_vector.h" +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + + + +template +class Step4 +{ +public: + Step4(); + void + run(); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(int cycle); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + + LinearAlgebra::TpetraWrappers::SparseMatrix system_matrix; + + LinearAlgebra::TpetraWrappers::Vector solution; + LinearAlgebra::TpetraWrappers::Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const; +}; + + + +template +double +RightHandSide::value(const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i = 0; i < dim; ++i) + return_value += 4 * std::pow(p[i], 4); + + return return_value; +} + + + +template +double +BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4() + : fe(1) + , dof_handler(triangulation) +{} + + +template +void +Step4::make_grid() +{ + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(5); +} + + + +template +void +Step4::setup_system() +{ + dof_handler.distribute_dofs(fe); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + BoundaryValues(), + constraints); + constraints.close(); + + DynamicSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, c_sparsity, constraints, false); + system_matrix.reinit(c_sparsity); + IndexSet is(dof_handler.n_dofs()); + is.add_range(0, dof_handler.n_dofs()); + solution.reinit(is); + system_rhs.reinit(is); +} + + +template +void +Step4::assemble_system() +{ + QGauss quadrature_formula(fe.degree + 1); + + const RightHandSide right_hand_side; + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += + (fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point)); + + cell_rhs(i) += + (fe_values.shape_value(i, q_point) * + right_hand_side.value(fe_values.quadrature_point(q_point)) * + fe_values.JxW(q_point)); + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + } + system_matrix.compress(VectorOperation::add); +} + + + +template +void +Step4::solve(int cycle) +{ + deallog.push(Utilities::int_to_string(dof_handler.n_dofs(), 5)); + + { + deallog.push("Identity"); + static constexpr std::array lower{{49, 100}}; + LinearAlgebra::TpetraWrappers::PreconditionIdentity preconditioner; + solution = 0; + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + Teuchos::ParameterList precon_params; + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + // Make sure the general Ifpack preconditioner works at all. + // For this we use the point relaxation (SSOR) + // with default parameters, resulting in symmetric Gauss-Seidel (SGS) + deallog.push("CustomSGS"); + static constexpr std::array lower{{49, 100}}; + LinearAlgebra::TpetraWrappers::PreconditionIfpack preconditioner( + "RELAXATION"); + solution = 0; + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + Teuchos::ParameterList precon_params; + precon_params.set("relaxation: type", "Symmetric Gauss-Seidel"); + preconditioner.set_parameter_list(precon_params); + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("Jacobi"); + static constexpr std::array lower{{49, 100}}; + LinearAlgebra::TpetraWrappers::PreconditionJacobi preconditioner; + solution = 0; + + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("l1Jacobi"); + static constexpr std::array lower{{49, 100}}; + LinearAlgebra::TpetraWrappers::PreconditionL1Jacobi preconditioner; + solution = 0; + + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("l1GaussSeidel"); + static constexpr std::array lower{{31, 62}}; + LinearAlgebra::TpetraWrappers::PreconditionL1GaussSeidel + preconditioner; + solution = 0; + + SolverControl solver_control(1000, 1e-10); + SolverBicgstab> solver( + solver_control); + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("SOR"); + static constexpr std::array lower{{31, 62}}; + LinearAlgebra::TpetraWrappers::PreconditionSOR preconditioner; + solution = 0; + + SolverControl solver_control(1000, 1e-5); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("SSOR"); + static constexpr std::array lower{{40, 77}}; + LinearAlgebra::TpetraWrappers::PreconditionSSOR preconditioner; + solution = 0; + + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("BlockJacobi"); + static constexpr std::array lower{{73, 145}}; + LinearAlgebra::TpetraWrappers::PreconditionBlockJacobi + preconditioner; + LinearAlgebra::TpetraWrappers::PreconditionBlockJacobi< + double>::AdditionalData data; + data.n_local_parts = 16; + solution = 0; + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix, data); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("BlockSOR"); + static constexpr std::array lower{{18, 37}}; + LinearAlgebra::TpetraWrappers::PreconditionBlockSOR preconditioner; + LinearAlgebra::TpetraWrappers::PreconditionBlockSOR::AdditionalData + data; + data.n_local_parts = 16; + data.omega = 0.8; + solution = 0; + SolverControl solver_control(1000, 1e-5); + SolverBicgstab> solver( + solver_control); + preconditioner.initialize(system_matrix, data); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("BlockSSOR"); + static constexpr std::array lower{{30, 59}}; + LinearAlgebra::TpetraWrappers::PreconditionBlockSSOR preconditioner; + LinearAlgebra::TpetraWrappers::PreconditionBlockSSOR::AdditionalData + data; + data.n_local_parts = 16; + data.omega = 1.2; + solution = 0; + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix, data); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + static constexpr std::array lower{{30, 56}}; + deallog.push("ILU"); + LinearAlgebra::TpetraWrappers::PreconditionILU preconditioner; + solution = 0; + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("ILUT"); + static constexpr std::array lower{{11, 19}}; + LinearAlgebra::TpetraWrappers::PreconditionILUT preconditioner; + LinearAlgebra::TpetraWrappers::PreconditionILUT::AdditionalData + data; + data.ilut_drop = 1e-6; + data.ilut_fill = 3; + solution = 0; + SolverControl solver_control(1000, 1e-5); + SolverBicgstab> solver( + solver_control); + preconditioner.initialize(system_matrix, data); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + { + deallog.push("Chebyshev"); + static constexpr std::array lower{{23, 46}}; + LinearAlgebra::TpetraWrappers::PreconditionChebyshev preconditioner; + LinearAlgebra::TpetraWrappers::PreconditionChebyshev::AdditionalData + data; + data.max_eigenvalue = 2.5; + data.degree = 3; + solution = 0; + SolverControl solver_control(1000, 1e-10); + SolverCG> solver( + solver_control); + preconditioner.initialize(system_matrix, data); + check_solver_within_range( + solver.solve(system_matrix, solution, system_rhs, preconditioner), + solver_control.last_step(), + lower[cycle], + lower[cycle] + 2); + deallog.pop(); + } + + // { + // deallog.push("Direct"); + // TrilinosWrappers::PreconditionBlockwiseDirect preconditioner; + // solution = 0; + // SolverControl solver_control(1000, 1e-10); + // SolverCG<> solver(solver_control); + // preconditioner.initialize(system_matrix); + // check_solver_within_range( + // solver.solve(system_matrix, solution, system_rhs, preconditioner), + // solver_control.last_step(), + // 1U, + // 1U); + // deallog.pop(); + // } + + { + // Make sure the general Ifpack preconditioner throws an Exception for an + // unsupported solver. + // IGLOO (incomplete global least optimal option) is of course a made up + // preconditioner name. + + deallog.push("CustomIGLOO"); + try + { + LinearAlgebra::TpetraWrappers::PreconditionIfpack + preconditioner("IGLOO"); + } + catch (dealii::ExceptionBase &exc) + { + deallog << "Error: " << std::endl; + exc.print_info(deallog.get_file_stream()); + } + deallog.pop(); + } + + deallog.pop(); +} + + + +template +void +Step4::run() +{ + for (unsigned int cycle = 0; cycle < 2; ++cycle) + { + if (cycle == 0) + make_grid(); + else + triangulation.refine_global(1); + + setup_system(); + assemble_system(); + solve(cycle); + } +} + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + try + { + Step4<2> test; + test.run(); + } + catch (const std::exception &exc) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos_tpetra/precondition.output b/tests/trilinos_tpetra/precondition.output new file mode 100644 index 0000000000..30f6004ca7 --- /dev/null +++ b/tests/trilinos_tpetra/precondition.output @@ -0,0 +1,43 @@ + +DEAL:01089:Identity::Solver stopped within 49 - 51 iterations +DEAL:01089:CustomSGS::Solver stopped after 41 iterations +DEAL:01089:Jacobi::Solver stopped within 49 - 51 iterations +DEAL:01089:l1Jacobi::Solver stopped within 49 - 51 iterations +DEAL:01089:l1GaussSeidel::Solver stopped after 48 iterations +DEAL:01089:SOR::Solver stopped after 25 iterations +DEAL:01089:SSOR::Solver stopped within 40 - 42 iterations +DEAL:01089:BlockJacobi::Solver stopped after 58 iterations +DEAL:01089:BlockSOR::Solver stopped after 16 iterations +DEAL:01089:BlockSSOR::Solver stopped after 28 iterations +DEAL:01089:ILU::Solver stopped within 30 - 32 iterations +DEAL:01089:ILUT::Solver stopped after 5 iterations +DEAL:01089:Chebyshev::Solver stopped after 38 iterations +DEAL:01089:CustomIGLOO::Error: + You tried to select the preconditioner type +but this preconditioner is not supported by Trilinos/Ifpack22 +due to one of the following reasons: +* This preconditioner does not exist +* This preconditioner has a specialized constructor not supported by the Ifpack2 Factory. +* This preconditioner is not (yet) supported by Trilinos/Ifpack2 +* Trilinos/Ifpack2 was not configured for its use. +DEAL:04225:Identity::Solver stopped within 100 - 102 iterations +DEAL:04225:CustomSGS::Solver stopped after 78 iterations +DEAL:04225:Jacobi::Solver stopped within 100 - 102 iterations +DEAL:04225:l1Jacobi::Solver stopped within 100 - 102 iterations +DEAL:04225:l1GaussSeidel::Solver stopped after 98 iterations +DEAL:04225:SOR::Solver stopped after 46 iterations +DEAL:04225:SSOR::Solver stopped within 77 - 79 iterations +DEAL:04225:BlockJacobi::Solver stopped after 82 iterations +DEAL:04225:BlockSOR::Solver stopped after 20 iterations +DEAL:04225:BlockSSOR::Solver stopped after 40 iterations +DEAL:04225:ILU::Solver stopped within 56 - 58 iterations +DEAL:04225:ILUT::Solver stopped after 10 iterations +DEAL:04225:Chebyshev::Solver stopped after 75 iterations +DEAL:04225:CustomIGLOO::Error: + You tried to select the preconditioner type +but this preconditioner is not supported by Trilinos/Ifpack22 +due to one of the following reasons: +* This preconditioner does not exist +* This preconditioner has a specialized constructor not supported by the Ifpack2 Factory. +* This preconditioner is not (yet) supported by Trilinos/Ifpack2 +* Trilinos/Ifpack2 was not configured for its use. \ No newline at end of file