Adhoc_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 <Adhoc_JFNK.h>
#include <Schema_Euler_Semi_Implicite.h>
#include <Schema_Cahn_Hilliard.h>
#include <Fermeture_Thermo_base.h>
#include <Potentiel_Chimique_base.h>
 
// XD adhoc_jfnk solveur_non_lineaire adhoc_jfnk BRACE Hand-rolled (matrix-free) Newton-Krylov GMRES non-linear solver
// XD_CONT for the Cahn-Hilliard phase-field equation, used when PETSc is not available.
Implemente_instanciable(Adhoc_JFNK,"Adhoc_JFNK",Solveur_non_lineaire);
 
Sortie& Adhoc_JFNK::printOn(Sortie& os) const
{
  os << "--------------------------------------------------" << finl;
  os << "SOLVER TYPE = " << que_suis_je() << " " << le_nom() << finl;
  os << " -- residual threshold = " << tol_ << finl ;
  os << " -- residue min = " << tol_min_ << finl;
  os << " -- residue max = " << tol_max_ << finl;
  os << " -- number of iterations = " << maxit_ << finl;
  os << " -- Krylov space dimension = " << nkr_ << finl;
  os << "--------------------------------------------------" << finl;
  return os ;
}
 
Entree& Adhoc_JFNK::readOn(Entree& is)
{
  Cerr << "Reading parameters of " << que_suis_je() << " solver ..." << finl;
  nommer("Captain Adhoc");
 
  Param param(que_suis_je());
  param.ajouter("seuil_residu", &tol_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Relative convergence threshold of the residual.
  param.ajouter("dimension_espace_de_krylov", &nkr_, Param::REQUIRED);  // XD_ADD_P entier
  // XD_CONT Dimension of the GMRES Krylov space.
  param.ajouter("nb_iterations", &maxit_, Param::REQUIRED);  // XD_ADD_P entier
  // XD_CONT Maximum number of GMRES iterations.
  param.ajouter("residu_min", &tol_min_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Minimum (absolute) convergence threshold of the residual.
  param.ajouter("residu_max", &tol_max_, Param::REQUIRED);  // XD_ADD_P floattant
  // XD_CONT Maximum (absolute) convergence threshold of the residual.
 
  param.lire_avec_accolades_depuis(is);
 
  printOn(Cout);
 
  return is;
}
 
/*! @brief Construire le jacobien
 *
 */
DoubleTab Adhoc_JFNK::jacobian_vect(const DoubleTab& v0, DoubleTab& c_theta, Cahn_Hilliard& equation)
{
  const double delta = 1.e-5;
 
  DoubleTrav c_autre_direction(c_theta);
 
  c_autre_direction = v0;
  c_autre_direction *= delta;
  c_autre_direction += c_theta; // c_autre_direction = c_theta + delta * v0
 
  DoubleTrav v1(v0);
  DoubleTab v2(v1);
  v2 = 0.;
 
  // Construction de la differentielle (dans la direction v0)
  v2 = equation.fonction_residu(c_autre_direction);
  v1 = equation.fonction_residu(c_theta);
 
  v2 -= v1;
  v2 /= delta; // v2 = [ H(c_autre_direction) - H(c_theta) ] / delta
 
  return v2;
}
 
/*! @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 Adhoc_JFNK::iterer_eqn(Equation_base& equation, const DoubleTab& inconnue, DoubleTab& result, double dt, int numero_iteration, int& ok)
{
  // Pour le moment, Adhoc_JFNK n'est utilisé que pour l'équation de Cahn_Hilliard
  if (sub_type(Cahn_Hilliard, equation))
    {
      Cout << " -------- [Adhoc_JFNK] Solving..." << finl;
    }
  else
    {
      Cerr << "[Adhoc_JFNK] Cannot solve any other equation than Cahn_Hilliard yet. Needs to be generalized." << finl;
      exit();
    }
  Cahn_Hilliard& ch = ref_cast(Cahn_Hilliard, equation);
 
  DoubleTrav c_n(inconnue);
  DoubleTrav c_theta(inconnue);
 
  c_n = inconnue;
  c_theta = result;
 
  double res;
 
  DoubleTabs v(nkr_); // Krylov vectors
  DoubleTrav h(nkr_ + 1, nkr_); // Hessenberg matrix of coefficients
  DoubleTrav r(nkr_ + 1);
  DoubleTrav v0(c_theta);
  DoubleTrav v1(c_theta);
 
  // v0 = H(c_n)
  v0 = ch.fonction_residu(c_theta);
  v0 *= -1.;
  v0.echange_espace_virtuel();
  res = 0.;
  const int nb_param = v0.line_size();
 
  // Calcul du premier résidu
  res = mp_norme_vect(v0); // Norme L2 du champ v0 = H(c_n) (on veut que cela tende vers 0)
 
  Cout << " - At it = " << 0 << ", residu scalar = " << res << finl;
 
  if (!Process::me())
    {
      SFichier fichier_residu("residu.out", ios::app);
      fichier_residu << statistics().get_time_since_last_open(STD_COUNTERS::total_execution_time)
                     << "\t" << ch.schema_temps().temps_courant() << "\t" << (double)res << finl;
    }
 
  if (res < tol_min_)
    return 0; // nothing to do
 
  tol_min_ = (tol_min_ < res * tol_) ? res * tol_ : tol_min_;
  tol_min_ = (tol_min_ < tol_max_) ? tol_min_ : tol_max_;
 
  // Ancienne façon de fixer la tolérance : voir avec Didier Jamet...
  Cout << "Stopping rule scalar : " << tol_min_ << finl;
 
  // iterations
  for (int it = 1; it < maxit_+1; it++)
    {
 
      // Dimension espace de Krylov
      int nk = nkr_;
      v0 /= res; // v0 = H(c_n) normalisé
      r = 0.;
      r[0] = res; // Stockage du résidu
      h = 0.;
 
      // Début orthogonalisation d'Arnoldi
      for (int j = 0; j < nkr_; j++)
        {
          v0.echange_espace_virtuel();
          // Remplissage des vecteurs de base de l'espace de Krylov (stockés dans v)
          v[j] = v0;
 
          // Calcul de J*v0 (approximation de la jacobienne par DF)
          v1 = jacobian_vect(v0, c_theta, ch);
          v0 = v1;
 
          // Remplissage de la matrice de Hessenberg
          for (int i = 0; i <= j; i++)
            {
              // Produit scalaire : H_ij = v_i . v_(j+1)
              DoubleTab& vvi = v[i];
              h(i,j) += mp_prodscal(v0,vvi);
 
              // v_(j+1) = v_(j+1) - H_ij v_i
              for (int ii = 0; ii < v0.dimension(0); ii++)
                for (int k = 0; k < nb_param; k++)
                  v0(ii, k) -= h(i,j) * vvi(ii, k);
 
              v0.echange_espace_virtuel();
            }
 
          double tem = mp_norme_vect(v0);
 
          // H_(j+1,j) = || v_(j+1) ||
          h(j + 1, j) = tem;
 
          // Si tem = 0, alors || v_(j+1) || = 0, donc l'espace de Krylov a cessé de grandir --> on sort
          // Happy breakdown du GMRES
          if (tem < tol_min_)
            {
              nk = j + 1;
              goto l5naire;
            }
          // Normalisation de v_(j+1)
          v0 /= tem;
        }
 
      // Résolution GMRES
l5naire:
      // Utilisation de la rotation de Givens
      // pour triangulariser la matrice de Hessenberg (se débarrasser de la diagonale inférieure)
      for (int i = 0; i < nk; i++)
        {
          int iplus1 = i + 1;
          double tem = 1. / sqrt(h(i, i) * h(i, i) + h(iplus1, i) * h(iplus1, i));
          double ccos = h(i, i) * tem;
          double ssin = -h(iplus1, i) * tem;
          for (int j = i; j < nk; j++)
            {
              tem = h(i, j);
              h(i, j) = ccos * tem - ssin * h(iplus1,j);
              h(iplus1, j) = ssin * tem + ccos * h(iplus1,j);
            }
          r[iplus1] = ssin * r[i];
          r[i] *= ccos;
        }
 
      // Solution du système linéaire
      for (int i = nk - 1; i >= 0; i--)
        {
          r[i] /= h(i, i);
          for (int i0 = i - 1; i0 >= 0; i0--)
            r[i0] -= h(i0, i) * r[i];
        }
 
      // La solution est formée d'une combinaison linéaire des vecteurs de base de l'espace de Krylov
      for(int i=0; i<nk; i++)
        {
          DoubleTrav vvi = v[i];
          vvi = v[i];
          for (int ii = 0; ii < c_theta.dimension(0); ii++)
            for (int k = 0; k < nb_param; k++)
              c_theta(ii, k) += r[i] * vvi(ii, k);
        }
      c_theta.echange_espace_virtuel();
 
      // v0 = H(c_k+1)
      v0 = ch.fonction_residu(c_theta);
      v0 *= -1.;
 
      res = 0.;
      res=mp_norme_vect(v0);
 
      Cout << " - At it = " << it << ", residu scalar = " << res << finl;
 
      double temps_exec = statistics().get_time_since_last_open(STD_COUNTERS::total_execution_time);
      if (!Process::me())
        {
          SFichier fichier_residu("residu.out", ios::app);
          fichier_residu << temps_exec << "\t" << ch.schema_temps().temps_courant() << "\t" << (double)res << finl;
        }
 
      if (res < tol_min_)
        {
          // Pour être sûr que ch.inconnue().valeurs() est mis à jour
          result = c_theta;
          return it;
        }
      else if (it == maxit_)
        {
          result = c_theta;
          // Pour être sûr que ch.inconnue().valeurs() est mis à jour
          Cerr << "Stopped before convergence" << finl;
        }
    }
 
  return true;
}