CPreconditioner.hpp

Go to the documentation of this file.

/*!
 * \file CPreconditioner.hpp
 * \brief Classes related to linear preconditioner wrappers.
 *        The actual operations are currently implemented mostly by CSysMatrix.
 * \author F. Palacios, J. Hicken, T. Economon
 * \version 8.5.0 "Harrier"
 *
 * SU2 Project Website: https://su2code.github.io
 *
 * The SU2 Project is maintained by the SU2 Foundation
 * (http://su2foundation.org)
 *
 * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md)
 *
 * SU2 is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * SU2 is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
 */
 
#pragma once
 
#include <functional>
#include "../CConfig.hpp"
#include "../geometry/CGeometry.hpp"
#include "CSysVector.hpp"
#include "CSysMatrix.hpp"
 
/// \addtogroup SpLinSys
/// @{
 
/*!
 * \brief Applies a preconditioner that only has a host implementation to vectors that live
 * on the device: bring the input down, apply, put the result back.
 * \note This is what keeps LU-SGS, Linelet and PaStiX usable on the GPU path. The
 * transfers are issued by one thread with the team synchronized around them, the apply
 * itself is the normal OpenMP parallel host code.
 */
template <class ScalarType, class Apply>
inline void ApplyPreconditionerOnHost(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v, Apply&& apply) {
#ifdef SU2_ENABLE_CUDA_KERNELS
  if constexpr (su2_gpu_capable_v<ScalarType>) {
    if (VecExpr::UseDeviceExpressions()) {
      /*--- The host code must not see the device pointers of any expression it builds, so
       * the switch is flipped for the duration of the apply. It is written inside the
       * regions, by one thread, and published to the team by the trailing barrier. ---*/
      SU2_DEVICE_REGION(u.DtHTransfer(); VecExpr::SetUseDeviceExpressions(false);)
 
      apply();
 
      SU2_DEVICE_REGION(VecExpr::SetUseDeviceExpressions(true); v.HtDTransfer();)
      return;
    }
  }
#endif
  apply();
}
 
/*!
 * \brief Mirror of ApplyPreconditionerOnHost: applies a device preconditioner to host vectors.
 * \note For callers that drive the Krylov solvers themselves and so never went through
 * CSysSolve::Solve, which is what normally leaves the vectors on the device (Newton-Krylov).
 * Device expressions are on for the duration so that a nested solve also uses the device copies.
 * Only \p u is uploaded, \p v is always overwritten by the apply.
 */
template <class ScalarType, class Apply>
inline void ApplyPreconditionerOnDevice(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v, bool useCuda,
                                        Apply&& apply) {
#ifdef SU2_ENABLE_CUDA_KERNELS
  if constexpr (su2_gpu_capable_v<ScalarType>) {
    if (useCuda && !VecExpr::UseDeviceExpressions()) {
      SU2_DEVICE_REGION(u.HtDTransfer(); VecExpr::SetUseDeviceExpressions(true);)
 
      apply();
 
      SU2_DEVICE_REGION(VecExpr::SetUseDeviceExpressions(false); v.DtHTransfer();)
      return;
    }
  }
#else
  (void)useCuda;
#endif
  apply();
}
 
/*!
 * \class CPreconditioner
 * \brief Abstract base class for defining a preconditioning operation.
 * \author J. Hicken.
 *
 * See the remarks regarding the CMatrixVectorProduct class. The same
 * idea applies here to the preconditioning operation.
 */
template <class ScalarType>
class CPreconditioner {
 public:
  /*!
   * \brief Destructor of the class
   */
  virtual ~CPreconditioner() = 0;
 
  /*!
   * \brief Overload of operator (), applies the preconditioner to "u" storing the result in "v".
   */
  virtual void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const = 0;
 
  /*!
   * \brief Generic "preprocessing" hook derived classes may implement to build the preconditioner.
   */
  virtual void Build() {}
 
  /*!
   * \brief Return true to identify the identity preconditioner, may allow some solvers to be more efficient.
   */
  virtual bool IsIdentity() const { return false; }
 
  /*!
   * \brief Factory method.
   */
  static CPreconditioner* Create(ENUM_LINEAR_SOLVER_PREC kind, CSysMatrix<ScalarType>& jacobian, CGeometry* geometry,
                                 const CConfig* config);
};
template <class ScalarType>
CPreconditioner<ScalarType>::~CPreconditioner() {}
 
/*!
 * \class CIdentityPreconditioner
 * \brief No-op preconditioner used when Krylov solvers run without preconditioning.
 * \note Also serves Q_IDENTITY: Build() requests quantization of the diagonal blocks, needed by
 * the matrix-vector product shared with the Krylov solver even though this preconditioner's own
 * operation is a no-op. CSysMatrix::QuantizeDiagonalBlocks() when quantization is off.
 */
template <class ScalarType>
class CIdentityPreconditioner final : public CPreconditioner<ScalarType> {
 private:
  CSysMatrix<ScalarType>& sparse_matrix;
 
 public:
  inline explicit CIdentityPreconditioner(CSysMatrix<ScalarType>& matrix_ref) : sparse_matrix(matrix_ref) {}
 
  CIdentityPreconditioner() = delete;
 
  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override { v = u; }
 
  inline bool IsIdentity() const override { return true; }
 
  inline void Build() override { sparse_matrix.QuantizeDiagonalBlocks(); }
};
 
/*!
 * \class CJacobiPreconditioner
 * \brief Specialization of preconditioner that uses CSysMatrix class.
 */
template <class ScalarType>
class CJacobiPreconditioner final : public CPreconditioner<ScalarType> {
 private:
  CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to matrix that defines the preconditioner. */
  CGeometry* geometry;                   /*!< \brief Pointer to geometry associated with the matrix. */
  const CConfig* config;                 /*!< \brief Pointer to problem configuration. */
 
 public:
  /*!
   * \brief Constructor of the class.
   * \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
   * \param[in] geometry_ref - Geometry associated with the problem.
   * \param[in] config_ref - Config of the problem.
   */
  inline CJacobiPreconditioner(CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref, const CConfig* config_ref)
      : sparse_matrix(matrix_ref) {
    if ((geometry_ref == nullptr) || (config_ref == nullptr))
      SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
    geometry = geometry_ref;
    config = config_ref;
  }
 
  /*!
   * \note This class cannot be default constructed as that would leave us with invalid Pointers.
   */
  CJacobiPreconditioner() = delete;
 
  /*!
   * \brief operator that defines the preconditioner operation
   * \param[in] u - CSysVector that is being preconditioned
   * \param[out] v - CSysVector that is the result of the preconditioning
   */
  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
    sparse_matrix.ComputeJacobiPreconditioner(u, v, geometry, config);
  }
 
  /*!
   * \note Request the associated matrix to build the preconditioner. Also serves Q_JACOBI:
   *       BuildJacobiPreconditioner() quantizes the diagonal blocks itself when the matrix was
   *       set up for it.
   */
  inline void Build() override { sparse_matrix.BuildJacobiPreconditioner(); }
};
 
/*!
 * \class CILUPreconditioner
 * \brief Specialization of preconditioner that uses CSysMatrix class
 */
template <class ScalarType>
class CILUPreconditioner final : public CPreconditioner<ScalarType> {
 private:
  CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to matrix that defines the preconditioner. */
  CGeometry* geometry;                   /*!< \brief Pointer to geometry associated with the matrix. */
  const CConfig* config;                 /*!< \brief Pointer to problem configuration. */
 
 public:
  /*!
   * \brief Constructor of the class.
   * \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
   * \param[in] geometry_ref - Geometry associated with the problem.
   * \param[in] config_ref - Config of the problem.
   */
  inline CILUPreconditioner(CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref, const CConfig* config_ref)
      : sparse_matrix(matrix_ref) {
    if ((geometry_ref == nullptr) || (config_ref == nullptr))
      SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
    geometry = geometry_ref;
    config = config_ref;
  }
 
  /*!
   * \note This class cannot be default constructed as that would leave us with invalid Pointers.
   */
  CILUPreconditioner() = delete;
 
  /*!
   * \brief Operator that defines the preconditioner operation.
   * \param[in] u - CSysVector that is being preconditioned.
   * \param[out] v - CSysVector that is the result of the preconditioning.
   */
  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
    /*--- No host bracket, ILU has a device implementation and ComputeILUPreconditioner dispatches to it. ---*/
    sparse_matrix.ComputeILUPreconditioner(u, v, geometry, config);
  }
 
  /*!
   * \note Request the associated matrix to build the preconditioner.
   */
  inline void Build() override { sparse_matrix.BuildILUPreconditioner(); }
};
 
/*!
 * \class CLU_SGSPreconditioner
 * \brief Specialization of preconditioner that uses CSysMatrix class.
 */
template <class ScalarType>
class CLU_SGSPreconditioner final : public CPreconditioner<ScalarType> {
 private:
  CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to matrix that defines the preconditioner. */
  CGeometry* geometry;                   /*!< \brief Pointer to geometry associated with the matrix. */
  const CConfig* config;                 /*!< \brief Pointer to problem configuration. */
 
 public:
  /*!
   * \brief Constructor of the class.
   * \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
   * \param[in] geometry_ref - Geometry associated with the problem.
   * \param[in] config_ref - Config of the problem.
   */
  inline CLU_SGSPreconditioner(CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref, const CConfig* config_ref)
      : sparse_matrix(matrix_ref) {
    if ((geometry_ref == nullptr) || (config_ref == nullptr))
      SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
    geometry = geometry_ref;
    config = config_ref;
  }
 
  /*!
   * \note This class cannot be default constructed as that would leave us with invalid Pointers.
   */
  CLU_SGSPreconditioner() = delete;
 
  /*!
   * \brief operator that defines the preconditioner operation.
   * \param[in] u - CSysVector that is being preconditioned.
   * \param[out] v - CSysVector that is the result of the preconditioning.
   */
  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
    sparse_matrix.ComputeLU_SGSPreconditioner(u, v, geometry, config);
  }
 
  /*!
   * \note Also serves Q_LU_SGS: quantizes the diagonal blocks, no-op for plain LU_SGS.
   */
  inline void Build() override { sparse_matrix.BuildLU_SGSPreconditioner(); }
};
 
/*!
 * \class CLineletPreconditioner
 * \brief Specialization of preconditioner that uses CSysMatrix class.
 */
template <class ScalarType>
class CLineletPreconditioner final : public CPreconditioner<ScalarType> {
 private:
  CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to matrix that defines the preconditioner. */
  CGeometry* geometry;                   /*!< \brief Pointer to geometry associated with the matrix. */
  const CConfig* config;                 /*!< \brief Pointer to problem configuration. */
 
 public:
  /*!
   * \brief Constructor of the class.
   * \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
   * \param[in] geometry_ref - Geometry associated with the problem.
   * \param[in] config_ref - Config of the problem.
   */
  inline CLineletPreconditioner(CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref, const CConfig* config_ref)
      : sparse_matrix(matrix_ref) {
    if ((geometry_ref == nullptr) || (config_ref == nullptr))
      SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
    geometry = geometry_ref;
    config = config_ref;
  }
 
  /*!
   * \note This class cannot be default constructed as that would leave us with invalid Pointers.
   */
  CLineletPreconditioner() = delete;
 
  /*!
   * \brief Operator that defines the preconditioner operation.
   * \param[in] u - CSysVector that is being preconditioned.
   * \param[out] v - CSysVector that is the result of the preconditioning.
   */
  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputeLineletPreconditioner(u, v, geometry, config); });
  }
 
  /*!
   * \note Request the associated matrix to build the preconditioner.
   */
  inline void Build() override { sparse_matrix.BuildLineletPreconditioner(geometry, config); }
};
 
/*!
 * \class CPastixPreconditioner
 * \brief Specialization of preconditioner that uses PaStiX to factorize a CSysMatrix.
 */
template <class ScalarType>
class CPastixPreconditioner final : public CPreconditioner<ScalarType> {
 private:
  CSysMatrix<ScalarType>& sparse_matrix; /*!< \brief Pointer to the matrix. */
  CGeometry* geometry;                   /*!< \brief Geometry associated with the problem. */
  const CConfig* config;                 /*!< \brief Configuration of the problem. */
  unsigned short kind_fact;              /*!< \brief The type of factorization desired. */
 
 public:
  /*!
   * \brief Constructor of the class
   * \param[in] matrix_ref - Matrix reference that will be used to define the preconditioner.
   * \param[in] geometry_ref - Associated geometry.
   * \param[in] config_ref - Problem configuration.
   * \param[in] kind_factorization - Type of factorization required.
   */
  inline CPastixPreconditioner(CSysMatrix<ScalarType>& matrix_ref, CGeometry* geometry_ref, const CConfig* config_ref,
                               unsigned short kind_factorization)
      : sparse_matrix(matrix_ref) {
    if ((geometry_ref == nullptr) || (config_ref == nullptr))
      SU2_MPI::Error("Preconditioner needs to be built with valid references.", CURRENT_FUNCTION);
    geometry = geometry_ref;
    config = config_ref;
    kind_fact = kind_factorization;
  }
 
  /*!
   * \note This class cannot be default constructed as that would leave us with invalid Pointers.
   */
  CPastixPreconditioner() = delete;
 
  /*!
   * \brief Operator that defines the preconditioner operation.
   * \param[in] u - CSysVector that is being preconditioned.
   * \param[out] v - CSysVector that is the result of the preconditioning.
   */
  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
    ApplyPreconditionerOnHost(u, v, [&] { sparse_matrix.ComputePastixPreconditioner(u, v, geometry, config); });
  }
 
  /*!
   * \note Request the associated matrix to build the preconditioner.
   */
  inline void Build() override { sparse_matrix.BuildPastixPreconditioner(geometry, config, kind_fact); }
};
 
/*!
 * \class CAbstractPreconditioner
 * \brief Applies a std::function as the preconditioning operation.
 * \note This can be used to treat almost anything as a preconditioner.
 */
template <class ScalarType>
class CAbstractPreconditioner final : public CPreconditioner<ScalarType> {
 private:
  std::function<void(const CSysVector<ScalarType>&, CSysVector<ScalarType>&)> impl;
 
 public:
  CAbstractPreconditioner() = delete;
 
  template <class F>
  explicit CAbstractPreconditioner(const F& function) : impl(function) {}
 
  inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override { impl(u, v); }
};
 
template <class ScalarType>
CPreconditioner<ScalarType>* CPreconditioner<ScalarType>::Create(ENUM_LINEAR_SOLVER_PREC kind,
                                                                 CSysMatrix<ScalarType>& jacobian, CGeometry* geometry,
                                                                 const CConfig* config) {
  CPreconditioner<ScalarType>* prec = nullptr;
 
  switch (kind) {
    case IDENTITY:
    case Q_IDENTITY:
      prec = new CIdentityPreconditioner<ScalarType>(jacobian);
      break;
    case JACOBI:
    case Q_JACOBI:
      prec = new CJacobiPreconditioner<ScalarType>(jacobian, geometry, config);
      break;
    case LINELET:
      prec = new CLineletPreconditioner<ScalarType>(jacobian, geometry, config);
      break;
    case LU_SGS:
    case Q_LU_SGS:
      prec = new CLU_SGSPreconditioner<ScalarType>(jacobian, geometry, config);
      break;
    case ILU:
      prec = new CILUPreconditioner<ScalarType>(jacobian, geometry, config);
      break;
    case PASTIX_ILU:
    case PASTIX_LU_P:
    case PASTIX_LDLT_P:
      prec = new CPastixPreconditioner<ScalarType>(jacobian, geometry, config, kind);
      break;
  }
 
  return prec;
}
 
/// @}