SNES_JFNK.cpp

/****************************************************************************
* Copyright (c) 2021, 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.
*
*****************************************************************************/
 
#include <SNES_JFNK.h>
#include <Cahn_Hilliard.h>
#include <Schema_Cahn_Hilliard.h>
#include <Debog.h>
#include <EcrFicPartageBin.h>
#include <petsc_for_kernel.h>
#ifdef PETSCKSP_H
// Include PETSC library
#include <petscsnes.h>
#endif
 
// XD snes_jfnk solveur_non_lineaire snes_jfnk BRACE Matrix-free Newton-Krylov (JFNK) non-linear solver built on the
// XD_CONT PETSc SNES framework, used by the Cahn-Hilliard phase-field time schemes.
Implemente_instanciable(SNES_JFNK,"SNES_JFNK",Solveur_non_lineaire);
 
Sortie& SNES_JFNK::printOn(Sortie& os) const
{
  os << "--------------------------------------------------" << finl;
  os << "SOLVER TYPE = " << que_suis_je() << " " << le_nom() << finl;
  os << " -- absolute tolerance = " << atol_ << finl ;
  os << " -- relative tolerance = " << rtol_ << finl;
  os << " -- stagnation tolerance = " << stol_ << finl;
  os << " -- max number of iterations of Newton = " << maxit_newton_ << finl;
  os << " -- max number of iterations of GMRES = " << maxit_gmres_ << finl;
  os << " -- dimension of Krylov space = " << nkr_ << finl;
  os << " -- Matrix Free: <<erel>> parameter to compute finite difference step = " << erel_ << finl;
  os << " -- Matrix Free: <<umin>> parameter to compute finite difference step = " << umin_ << finl;
  os << "--------------------------------------------------" << finl;
  return os ;
}
 
Entree& SNES_JFNK::readOn(Entree& is)
{
  Cerr << "Reading parameters of " << que_suis_je() << " solver ..." << finl;
  nommer("SNES");
 
  Param param(que_suis_je());
  param.ajouter("absolute_tol", &atol_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Absolute residual tolerance of the Newton iterations.
  param.ajouter("relative_tol", &rtol_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Relative residual tolerance of the Newton iterations.
  param.ajouter("stagn_tol", &stol_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Stagnation (step-size) tolerance of the Newton iterations.
  param.ajouter("max_it_newton", &maxit_newton_, Param::REQUIRED);  // XD_ADD_P entier
  // XD_CONT Maximum number of Newton iterations.
  param.ajouter("max_it_GMRES", &maxit_gmres_, Param::REQUIRED);  // XD_ADD_P entier
  // XD_CONT Maximum number of GMRES iterations of the inner linear solve.
  param.ajouter("dimension_espace_de_krylov", &nkr_, Param::REQUIRED);  // XD_ADD_P entier
  // XD_CONT Dimension (restart) of the GMRES Krylov space.
  param.ajouter("matrix_free_erel", &erel_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Relative step used to build the matrix-free finite-difference Jacobian.
  param.ajouter("matrix_free_umin", &umin_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Minimum step used to build the matrix-free finite-difference Jacobian.
 
  param.lire_avec_accolades_depuis(is);
 
  printOn(Cout);
 
  return is;
}
 
/*! @brief FormFunction de PETSc : calcule le résidu du JFNK
 *
 *  Le vecteur v en entrée est uniquement la concentration
 */
#ifdef PETSCKSP_H
PetscErrorCode FormFunction(SNES snes, Vec v, Vec f, void* ctx)
{
  std::pair<Equation_base&, ArrOfTID&>* pair = static_cast< std::pair<Equation_base&, ArrOfTID&> * >(ctx);
 
  Cahn_Hilliard& ch = ref_cast(Cahn_Hilliard, pair->first);
  auto& ix = pair->second;
 
  DoubleTrav c_theta(ch.inconnue().valeurs());
  DoubleTrav residu(c_theta);
 
  VecGetValues(v, ix.size_array(), (PetscInt*)ix.addr(), c_theta.addr());
 
  c_theta.echange_espace_virtuel(); // Virt values are zeros ...
 
  residu = ch.fonction_residu(c_theta);
 
  VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
  VecSetValues(f, ix.size_array(), (PetscInt*)ix.addr(), residu.addr(), INSERT_VALUES);
 
  VecAssemblyBegin(f);
  VecAssemblyEnd(f);
 
  return 0;
}
 
PetscErrorCode SNESMonitor(SNES snes, PetscInt its, PetscReal norm, void* ctx)
{
  Equation_base* eqn = static_cast<Equation_base*>(ctx);
 
  // statistiques().set_debug_level(9);
  double temps_exec = statistics().get_time_since_last_open(STD_COUNTERS::total_execution_time);
 
  Cerr << "[SNES] Iteration " << its << " : ||F(x)|| = " << norm << finl;
 
  if (!Process::me())
    {
      SFichier fichier_residu("residu.out", ios::app);
      fichier_residu << temps_exec << "\t" << eqn->schema_temps().temps_courant() << "\t" << (double)norm << finl;
    }
 
  return 0;
}
 
PetscErrorCode KSPMonitor(KSP ksp, PetscInt its, PetscReal rnorm, void* ctx)
{
  Cerr << "  [KSP] Iteration " << its << " : residu = " << rnorm << finl;
 
  return 0;
}
#endif
/*! @brief Permet de résoudre l'équation non linéaire H(c) = c_theta - c_n - theta * dt * M^-1 * D * mu_theta = 0
 *
 * inconnue = passe
 * result = present
 */
bool SNES_JFNK::iterer_eqn(Equation_base& equation, const DoubleTab& inconnue, DoubleTab& result, double dt, int numero_iteration, int& ok)
{
#ifdef PETSCKSP_H
  Debog::verifier("SNES_JFNK::iterer_eqn avant result ", result);
  // Pour le moment, Adhoc_JFNK n'est utilisé que pour l'équation de Cahn_Hilliard
  if (sub_type(Cahn_Hilliard, equation))
    {
      Cout << " -------- [SNES_JFNK] Solving..." << finl;
    }
  else
    {
      Cerr << "[SNES_JFNK] Cannot solve any other equation than Cahn_Hilliard yet. Needs to be generalized." << finl;
      exit();
    }
  Cahn_Hilliard& ch = ref_cast(Cahn_Hilliard, equation);
 
  if (ix.size_array() == 0)
    construit_renum(inconnue);
 
  DoubleTrav residu(result);
  residu = inconnue;
 
  // Variables SNES ########################################
  SNES            snes;     // Contexte SNES
  Vec             v, r;     // Vecteurs : solution, résidu
  KSP ksp;
  Mat J;
  PC pc;
  // #######################################################
 
  // INITIALISATION DU VECTEUR D'INCONNUES v
  VecCreate(PETSC_COMM_WORLD, &v);
  VecSetSizes(v, nb_rows_, nb_rows_tot_);
  VecSetFromOptions(v);
 
  // INITIAL GUESS of v
  VecSetOption(v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
  VecSetValues(v, ix.size_array(), (PetscInt*)ix.addr(), residu.addr(), INSERT_VALUES);
  VecAssemblyBegin(v);
  VecAssemblyEnd(v);
 
  VecDuplicate(v, &r);
 
  // Vérification tailles
  PetscInt v_size, r_size;
  VecGetSize(v, &v_size);
  VecGetSize(r, &r_size);
//  assert(v_size == taille_vec && r_size == taille_vec);
 
  // ---------------------- Compute real values of tolerances ---------------------------
  // Appel explicite de ta fonction de résidu
  std::pair<Equation_base&, ArrOfTID&> pair = { ch, ix };
  FormFunction(nullptr, v, r, &pair);
 
  // Affiche le résidu (par exemple sa norme)
  PetscReal norm;
  VecNorm(r, NORM_2, &norm);
  PetscPrintf(PETSC_COMM_WORLD, "Résidu F(x) = %g\n", (double)norm);
 
  // ------------------------------------ SNES ------------------------------------------
  SNESCreate(PETSC_COMM_WORLD, &snes);
 
  SNESSetType(snes, SNESNEWTONLS);
 
  // ========== Associate the FormFunction to the SNES context ============
  SNESSetFunction(snes, r, FormFunction, &pair);
  // ======================================================================
 
  SNESLineSearch linesearch;
  SNESGetLineSearch(snes, &linesearch);                      // Newton Line Search ne marche pas dans notre cas
  // C'est pour ça qu'on le désactive ici
  SNESLineSearchSetType(linesearch, SNESLINESEARCHNONE);     // Désactive la recherche de pas
 
  // Checking if linesearch is desactivated
  const char* ls_type;
  SNESLineSearchGetType(linesearch, &ls_type);
  PetscPrintf(PETSC_COMM_WORLD, "Type of Line Search used: %s\n", ls_type); // Printing linesearch type
 
  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  // Crée la matrice matrix-free à partir du SNES
  MatCreateSNESMF(snes, &J);
 
  // Fixe epsilon sur la vraie matrice MF
  // MatMFFDSetBase(J, delta);
  MatMFFDSetFunctionError(J, erel_);
  MatMFFDSetType(J, "ds");  // default: ds (difference scaled), alternative: wp (Walker–Pernice)
  MatMFFDDSSetUmin(J, umin_);
 
  SNESSetJacobian(snes, J, J, MatMFFDComputeJacobian, nullptr);
  // SNESSetJacobian(snes, J, nullptr, nullptr, nullptr);
 
  MatView(J, PETSC_VIEWER_STDOUT_WORLD);
  SNESView(snes, PETSC_VIEWER_STDOUT_WORLD);
  PetscInt local_size, global_size;
  VecGetLocalSize(v, &local_size);
  VecGetSize(v, &global_size);
 
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
                          "Proc %d : local=%d, global=%d\n",
                          (int)PetscGlobalRank, (int)local_size, (int)global_size);
  PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
                          "Proc %d : nb_rows=%d, nb_rows_tot=%d\n",
                          Process::me(), nb_rows_, (int)nb_rows_tot_);
  PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
  // Pas de Jacobien explicite => JFNK par défaut : SNES construit la matrice par DF entre colonnes
  // SNESSetJacobian(snes, nullptr, nullptr, SNESComputeJacobianDefault, nullptr);
  // SNESSetUseMatrixFree(snes, PETSC_TRUE, PETSC_TRUE); // Impose l'option Matrix-Free
 
  // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
  // Tolerance for Newton
  SNESSetTolerances(snes, atol_, rtol_, stol_, maxit_newton_, 100000);
 
  // TRES IMPORTANT : Force à continuer l'algorithme de Newton même si la méthode de Krylov n'a pas convergé
  SNESSetForceIteration(snes, PETSC_TRUE);
  SNESSetMaxLinearSolveFailures(snes, 200);
 
  // Vérification
  PetscBool force_it;
  SNESGetForceIteration(snes, &force_it);
  PetscPrintf(PETSC_COMM_WORLD, "Force iteration: %s\n", force_it ? "ON" : "OFF");
 
  // --------------------------------------- KSP ---------------------------------------------
  SNESGetKSP(snes, &ksp);
  KSPSetType(ksp, KSPGMRES);  // GMRES
  KSPGMRESSetRestart(ksp, nkr_);   // GMRES(m) avec m la dimension de l'espace de Krylov définie par l'utilisateur
  // KSPSetType(ksp, KSPCG); //  CG
  // KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); // XXX Luis : PAS BON DU TOUT
  // PetscOptionsSetValue(nullptr, "-ksp_gmres_orthogonalization", "classical"); // Facultatif : dans notre config, on ne peut pas changer le type d'orthogonalisation
  KSPSetTolerances(ksp, rtol_, atol_, 5000., maxit_gmres_);
 
  // -------------------------------- Préconditionneur --------------------------------------
  KSPGetPC(ksp, &pc);
  PCSetType(pc, PCNONE); // No preconditionning
  // PCSetFromOptions(pc);
 
  // Monitoring
  // KSPMonitorSet(ksp, KSPMonitor, nullptr, nullptr);
  SNESMonitorSet(snes, SNESMonitor, &ch, nullptr);
  //SNESMonitorSet(snes, MySNESMonitorBis, nullptr, nullptr);
 
  // =============================== Setting from options ===================================
  // KSPSetFromOptions(ksp);                              // Lire les options
 
  // Injecte les options comme si elles venaient de la ligne de commande
  //  PetscOptionsSetValue(nullptr, "-pc_type", "none");              // désactive le préconditionneur
  //  PetscOptionsSetValue(nullptr, "-snes_mf", nullptr);             // active matrix-free
  // PetscOptionsSetValue(nullptr, "-snes_mf_relstep", "1e-5");  // contrôle l’epsilon de l’approximation par différences finies
 
  PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_WORLD);  // Visu de la liste des options
  PetscOptionsLeft(nullptr);                             // Liste celles qui n'ont PAS été utilisées
 
  // Lecture des options ligne de commande de SNES
  SNESSetOptionsPrefix(snes, "mysolver_");
 
  SNESSetFromOptions(snes);
 
  SNESView(snes, PETSC_VIEWER_STDOUT_SELF);
 
  // ############################### SOLVING ################################################
  SNESSolve(snes, nullptr, v);                                                          //  #
  // ########################################################################################
 
 
  const char *reasonstr;
  SNESGetConvergedReasonString(snes, &reasonstr);
  Cout << "SNES converged for this reason: " << reasonstr << finl;
 
  residu = 0;
  VecGetValues(v, ix.size_array(), (PetscInt*)ix.addr(), residu.addr());
 
  result = residu;
  result.echange_espace_virtuel();
  Debog::verifier("SNES_JFNK::iterer_eqn apres result ", result);
 
  // Nettoyage
  VecDestroy(&v);
  VecDestroy(&r);
  MatDestroy(&J);
  SNESDestroy(&snes);
 
#endif
  return true;
}