Solv_Petsc.h

/****************************************************************************
* Copyright (c) 2026, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#ifndef Solv_Petsc_included
#define Solv_Petsc_included
 
#include <petsc_for_kernel.h>
#include <EChaine.h>
#include <Nom.h>
#undef setbit // Sinon conflit de ArrOfBit.h avec Petsc
#include <Solv_Externe.h>
#include <ArrOfBit.h>
#ifdef PETSCKSP_H
#include <petscdm.h>
#include <TRUSTTab.h>
#endif
#include <PCShell_base.h>
#include <TRUST_Deriv.h>
 
class Matrice_Morse_Sym;
class Matrice_Morse;
 
 
enum solveur_direct_ { no, mumps, superlu_dist, petsc, umfpack, pastix, cholmod, strumpack, cli };
extern bool gmres_right_unpreconditionned;
 
/* Struct to associate Petsc preconditionner to user-provided preconditioner */
typedef struct
{
  OWN_PTR(PCShell_base) pc_shell;
} PCstruct;
 
 
class Solv_Petsc : public Solv_Externe
{
 
  Declare_instanciable_sans_constructeur_ni_destructeur(Solv_Petsc);
 
public :
 
  inline Solv_Petsc();
  inline ~Solv_Petsc() override;
  // Method which be called to create solver from chaine_lue_
  inline void create_solver()
  {
    EChaine e(chaine_lue_);
    create_solver(e);
  }
  inline int solveur_direct() const override
  {
#ifdef PETSCKSP_H
    if (SolveurPetsc_==nullptr) Process::exit("Can't call to Solv_Petsc::solveur_direct() yet.");
#endif
    return (solveur_direct_>0);
  };
  int resoudre_systeme(const Matrice_Base&, const DoubleVect&, DoubleVect& ) override;
  inline int resoudre_systeme(const Matrice_Base& M, const DoubleVect& A, DoubleVect& B, int niter_max) override
  {
    return resoudre_systeme(M,A,B);
  };
  void create_solver(Entree&);
#ifdef PETSCKSP_H
  // To switch between solver definitions in one Solv_Petsc instance:
  void reset_solver(const Nom& name)
  {
    // ToDo: regler option_prefix_ et numero_solveur dans create_solver et initialize...
    reset();
    Cout << "Setting PETSc solver: " << name << finl;
    EChaine ech(name);
    create_solver(ech);
  }
#endif
  inline void reset();
#ifdef PETSCKSP_H
  inline Solv_Petsc& operator=(const Solv_Petsc&);
  inline Solv_Petsc(const Solv_Petsc&);
 
  inline void initialize();
  inline bool gpu() const
  {
    return gpu_;
  };
  inline PCstruct& get_precond_user()
  {
    if (!pc_user_.pc_shell) create_solver();
    return pc_user_;
  }
  inline bool amgx() const
  {
    return amgx_;
  };
  inline bool amgx_initialized()
  {
    return amgx_initialized_;
  };
#if PETSC_VERSION_GE(3,24,0)
  PetscErrorCode set_convergence_test(PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void**))
#else
  PetscErrorCode set_convergence_test(PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
#endif
  {
    if (SolveurPetsc_==nullptr) create_solver();
    return KSPSetConvergenceTest(SolveurPetsc_, converge, cctx, destroy);
  }
  // Timers:
  static PetscLogStage KSPSolve_Stage_;
  static PetscLogStage Create_Stage_;
  void set_rtol(const double& rtol)
  {
    seuil_relatif_ = rtol;
    if (SolveurPetsc_!=nullptr) KSPSetTolerances(SolveurPetsc_, seuil_relatif_, seuil_, (divtol_==0 ? PETSC_DEFAULT : divtol_), nb_it_max_);
  }
#endif
 
  public_for_cuda
#ifdef PETSCKSP_H
  virtual void Update_matrix(Mat& MatricePetsc, const Matrice_Morse& mat_morse); // Fill the (previously allocated) PETSc matrix with mat_morse coefficients
#endif
 
  static int instance;               // Nombre d'instances en cours de la classe
  static int numero_solveur;         // Compte les solveurs crees et utilises pour le prefix des options
 
protected :
#ifdef PETSCKSP_H
  using ArrOfPetscInt = TRUSTArray<PetscInt, PetscInt>;
 
  bool isViennaCLVector();
  bool isKokkosVector();
  void check_aij(const Matrice_Morse&);
  void Create_DM(const DoubleVect& ); // Construit un DM (Distributed Mesh)
  virtual void Create_objects(const Matrice_Morse&, int); // Construit differents objets PETSC dont matrice
  virtual void Create_vectors(const DoubleVect&); // Construit les vecteurs Petsc x et b
  virtual void Update_vectors(const DoubleVect& secmem, DoubleVect& solution); // Remplit les vecteurs Petsc x et b
  void Create_MatricePetsc(Mat&, int, const Matrice_Morse&); // Construit et remplit une matrice Petsc depuis la matrice_morse
  virtual void Update_solution(DoubleVect& solution);
  virtual int solve(ArrOfDouble& residual); // Solve Ax=b and return residual
  virtual void finalize() {};
  virtual bool detect_new_stencil(const Matrice_Morse&);
  bool nouveau_stencil()
  {
    return nouveau_stencil_;
  }; // ToDo: Remonter dans Solveur_Sys avec nouvelle_matrice
  bool enable_ksp_view();
  bool has_option(const Nom& option, Nom& current_value);
  int add_option(const Nom& option, const double& value, int cli = 0);
  int add_option(const Nom& option, const Nom& value, int cli = 0);
  void add_amgx_option(const Nom& key, const Nom& value, const std::string& comment="");
  void add_amgx_option(const Nom& key_value);
  void SaveObjectsToFile(const DoubleVect& b, DoubleVect& x);
  void RestoreMatrixFromFile();
  PetscInt compute_nb_rows_petsc(PetscInt);
 
  // Attributes
  double seuil_;
  double seuil_relatif_;
  double divtol_;
  bool nouveau_stencil_ = true;
 
  // Objets Petsc
  Mat MatricePetsc_;
  Vec SecondMembrePetsc_;
  Vec SolutionPetsc_;       // Globale solution
  KSP SolveurPetsc_;
  PC PreconditionneurPetsc_;
  PCstruct pc_user_;   /* user-defined preconditioner context */
  DM dm_;                       //description de champs PETSC
  int preconditionnement_non_symetrique_; // Drapeau sur la symetrie de la matrice de preconditionnement
  int nb_it_max_;       // Maximal number of iterations
  int convergence_with_nb_it_max_;  // Convergence decided with nb_it_max_ specified and not by seuil threshold
  int ignore_nb_it_max_;
  int controle_residu_;         // Verification si le residu ||Ax-B||<seuil
  int block_size_;              // Block size for SBAIJ matrix
  Nom factored_matrix_;     // Deal with the A=LU factorization on disk
  int mataij_;          // Force the use of a Mataij matrix
  Nom type_pc_;         // Preconditioner type
  Nom type_ksp_;        // KSP type
  Nom option_prefix_;       // Prefix des options CLI (permet de faire plusieurs solveurs en CLI)
  Nom amgx_options_;
 
  int petsc_nb_cpus_;       // Number of CPUs used to solve the PETSc matrix
  int petsc_cpus_selection_;    // Selection of CPUs used to solve the PETSc matrix
  int different_partition_; // Different partition for the TRUST matrix and the PETSc matrix
  int petsc_decide_;        // PETSc decide of the matrix partition
  Vec LocalSolutionPetsc_;  // Local solution in case of petsc_decide_=1
  VecScatter VecScatter_;   // Scatter context needed when petsc_decide_=1 to gather values of global to local solution
#endif
 
 
  int solveur_direct_ = no;          // Pour savoir si l'on manipule un solveur direct et non iteratif
  bool gpu_ = false;                    // Utilisation des solveurs GPU de PETSc
  bool amgx_ = false;           // Utilisation des solveurs GPU de AMGX
  const Nom config();    // Nom du fichier de config eventuel
  bool amgx_initialized_ = false;   // Amgx initialise
  // Options dev:
  bool ignore_new_nonzero_ = false;
  bool rebuild_matrix_ = false;
  bool allow_realloc_ = true;
  bool mat_ignore_zero_entries_= true;
  bool reduce_ram_ = false;
  bool verbose = false; // Timing
  bool reorder_matrix_ = false;
  // Reuse preconditioner (for Inexact Newton):
  int nb_it_previous_=-1;  // Previous number of iterations
  int reuse_preconditioner_nb_it_max_=-1; // Upper limit of iterations to reuse preconditioner
#ifdef PETSCKSP_H
  IS rowperm = nullptr, colperm = nullptr;
  IS inv_rowperm = nullptr, inv_colperm = nullptr;
#endif
};
 
#define NB_IT_MAX_DEFINED 10000
 
inline Solv_Petsc::Solv_Petsc()
{
#ifdef PETSCKSP_H
  initialize();
  instance++;
  //  Journal()<<"creation solv_petsc "<<instance<<finl;
#endif
}
 
inline Solv_Petsc::~Solv_Petsc()
{
  reset();
  instance--;
  // Journal()<<" destr solv_petsc "<<instance<<finl;
}
inline void Solv_Petsc::reset()
{
#ifdef PETSCKSP_H
  if (SolveurPetsc_!=nullptr)
    {
      KSPDestroy(&SolveurPetsc_);
      finalize();
    }
  if (MatricePetsc_!=nullptr)
    {
      // Destruction des vecteurs
      VecDestroy(&SecondMembrePetsc_);
      VecDestroy(&SolutionPetsc_);
      if (LocalSolutionPetsc_!=nullptr)
        {
          VecDestroy(&LocalSolutionPetsc_);
          VecScatterDestroy(&VecScatter_);
        }
      // Destruction matrice
      MatDestroy(&MatricePetsc_);
      // Destruction DM
      if (dm_!=nullptr)
        DMDestroy(&dm_);
      // Destruction IS
      if (rowperm!=nullptr)
        ISDestroy(&rowperm), ISDestroy(&inv_rowperm);
      if (colperm!=nullptr)
        ISDestroy(&colperm), ISDestroy(&inv_colperm);
    }
  initialize();
#endif
}
#ifdef PETSCKSP_H
#define _RTOL_MIN_ 1.e-24 // Would like to set to 1.e-14 but ALE test cases in TrioCFD needs rtol very low...
inline void Solv_Petsc::initialize()
{
  matrice_symetrique_=-1;
  preconditionnement_non_symetrique_=0;
  seuil_ = 0;
  seuil_relatif_ = _RTOL_MIN_;
  divtol_ = 0;
  nouveau_stencil_ = true;
  petsc_cpus_selection_ = 0;         // By default, 0 (no selection). 1 means: first petsc_nb_cpus_ CPUs is used, 2 means: every petsc_nb_cpus_ CPUs is used
  petsc_nb_cpus_ = Process::nproc(); // By default the number of processes
  different_partition_ = 0;          // By default, same matrix partition
  petsc_decide_ = 0;                 // By default, 0. 1 is OK but does NOT improve performance (rather decrease)
  convergence_with_nb_it_max_ = 0;
  ignore_nb_it_max_ = 0;
  nb_it_max_ = NB_IT_MAX_DEFINED;
  mataij_=0;
  factored_matrix_="";
  solveur_direct_=no;
  controle_residu_=0;
  gpu_=false;
  amgx_=false;
  amgx_initialized_=false;
  block_size_=1;
  option_prefix_="??";
  dm_=nullptr;
  MatricePetsc_ = nullptr;
  SecondMembrePetsc_ = nullptr;
  SolutionPetsc_ = nullptr;
  SolveurPetsc_ = nullptr;
  LocalSolutionPetsc_ = nullptr;
  VecScatter_ = nullptr;
  // Dev:
  ignore_new_nonzero_ = false;
  rebuild_matrix_ = false;
  allow_realloc_ = true;
  mat_ignore_zero_entries_ = true;
  reorder_matrix_ = false;
  reduce_ram_ = false;
  if (instance==-1)
    {
      // First initialization:
      instance=0;
      char version[256];
      PetscGetVersion(version,256);
      Cerr << "******************************************************************************************" << finl;
      Cerr << "Commands lines possible for " << version << ":" << finl;
      Cerr << "-ksp_view : to have some informations on the solver/preconditioner used by PETSc." << finl;
      Cerr << "-info : to have among other informations on storage of matrices of PETSc." << finl;
      Cerr << "-log_summary : to have at the end of the calculation, informations about performances and memory." << finl;
      Cerr << "-log_all : trace all PETSc calls." << finl;
      Cerr << "-malloc_dump : to have at the end of the calculation the memory not deallocated by PETSc." << finl;
      Cerr << "-help : to know all the commands lines of PETSc including those related to the solver/preconditioner selected." << finl;
      Cerr << "******************************************************************************************" << finl;
      if (!disable_TU)
        Cerr << "NB: if you want to disable the wrinting of the *_petsc.TU file then specify the disable_TU flag in your datafile before reading the block of schema in time." << finl;
      else
        Cerr << "Reading of disable_TU flag => Disable the writing of the *_petsc.TU file."<< finl;
      PetscLogStageRegister("CreateStage",&Create_Stage_);
      PetscLogStageRegister("KSPSolve",&KSPSolve_Stage_);
    }
}
 
inline Solv_Petsc::Solv_Petsc(const Solv_Petsc& org):Solv_Externe::Solv_Externe()
{
  initialize();
  instance++;
  // Journal()<<"copie solv_petsc "<<instance<<finl;
  set_read_matrix(org.read_matrix());
  gpu_=org.gpu();
  amgx_=org.amgx();
  option_prefix_=org.option_prefix_;
  // on relance la lecture ....
  EChaine recup(org.get_chaine_lue());
  readOn(recup);
}
 
inline Solv_Petsc& Solv_Petsc::operator=(const Solv_Petsc&)
{
  Cerr<<"Solv_Petsc::operator=(const Solv_Petsc&) is not coded."<<finl;
  exit();
  return (*this);
}
 
#endif
 
class option
{
public:
  int defined;
  Nom name;
};
 
class option_string : public option
{
public:
  option_string(Nom n, Nom v)
  {
    defined=0;
    name=n;
    value_=v;
  }
  Nom value_;
  inline Nom& value()
  {
    return value_;
  };
};
 
class option_int : public option
{
public:
  option_int(Nom n, int v)
  {
    defined=0;
    name=n;
    value_=v;
  }
  int value_;
  inline int& value()
  {
    return value_;
  };
};
 
class option_double : public option
{
public:
  option_double(Nom n, double v)
  {
    defined=0;
    name=n;
    value_=v;
  }
  double value_;
  inline double& value()
  {
    return value_;
  };
};
 
#endif