Source_Con_Phase_field_Binaire.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 <Convection_Diffusion_Phase_field.h>
#include <Source_Con_Phase_field_Binaire.h>
#include <Source_Qdm_VDF_Phase_field.h>
#include <Navier_Stokes_phase_field.h>
#include <Source_Con_Phase_field.h>
#include <Check_espace_virtuel.h>
#include <MD_Vector_composite.h>
#include <Milieu_Phase_field.h>
#include <Champ_Fonc_Tabule.h>
#include <Domaine_Cl_VDF.h>
#include <TRUSTTab_parts.h>
#include <Champ_Uniforme.h>
#include <Equation_base.h>
#include <Probleme_base.h>
#include <Matrix_tools.h>
#include <Domaine_VDF.h>
#include <Milieu_base.h>
#include <SolveurSys.h>
#include <Dimension.h>
#ifdef PETSCKSP_H
#include <petscsnes.h>
#endif
#include <EChaine.h>
#include <Param.h>
#include <SFichier.h>
#include <Perf_counters.h>
 
 
Implemente_instanciable(Source_Con_Phase_field_Binaire,"Source_Con_Phase_field_Binaire_VDF_P0_VDF",Source_Con_Phase_field);
 
// XD Source_Con_Phase_field_Binaire_VDF_P0_VDF Source_Con_Phase_field Source_Con_Phase_field_Binaire_VDF_P0_VDF INHERITS_BRACE Source for phase field
 
 
Sortie& Source_Con_Phase_field_Binaire::printOn(Sortie& s) const
{
  return s << que_suis_je();
}
 
Entree& Source_Con_Phase_field_Binaire::readOn(Entree& is)
{
  Cerr << "Source_Con_Phase_field_Binaire::readOn" << finl;
  Source_Con_Phase_field::readOn(is);
  if (!Process::me())
    {
      SFichier fichier_residu("residu.txt");
      fichier_residu << "# SAUVEGARDE DU RESIDU #" << finl;
      fichier_residu << "# Time \t Time simu \t Residu" << finl;
 
      SFichier fichier_eval("eval.txt");
      fichier_eval << "# Time simu \t Number eval FormFunctionOld" << finl;
    }
  return is;
}
 
DoubleTab& Source_Con_Phase_field_Binaire::laplacien(const DoubleTab& F, DoubleTab& resu) const
{
  const Navier_Stokes_std& eq_ns = ref_cast(Navier_Stokes_std, pb_->equation(0));
  const Operateur_Grad& opgrad = eq_ns.operateur_gradient();
  const Operateur_Div& opdiv = eq_ns.operateur_divergence();
 
  resu = 0.;
 
  // Dans cette partie, le laplacien est obtenu comme div.(grad F) avec
  // une application du solveur masse entre l'etape de gradient et de divergence
 
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
 
  prov_face = 0.;
  // Grad(F)
  opgrad.calculer(F, prov_face);
  // Application solveur masse
  eq_ns.solv_masse().appliquer(prov_face);
 
  // Div(M*Grad(F))
  opdiv.calculer(prov_face, resu);
 
  return resu;
}
 
DoubleTab& Source_Con_Phase_field_Binaire::div_kappa_grad(const DoubleTab& F, const DoubleTab& kappa_var, DoubleTab& resu) const
{
  const Navier_Stokes_std& eq_ns = ref_cast(Navier_Stokes_std, pb_->equation(0));
  const Operateur_Grad& opgrad = eq_ns.operateur_gradient();
  const Operateur_Div& opdiv = eq_ns.operateur_divergence();
 
  // const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  // const DoubleTab& c = eq_c.inconnue().valeurs();
 
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  const DoubleVect& volumes = domaine_VDF.volumes();
 
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
  prov_face = 0.;
  resu = 0.;
 
  // Grad(F)
  opgrad.calculer(F, prov_face);
 
  // M*Kappa*Grad(F)
  int ndeb = domaine_VDF.premiere_face_int();
  int nbfaces = domaine_VDF.nb_faces();
  int el0, el1;
  double kappa_face, vol0, vol1; // cface
 
  for (int fac = ndeb; fac < nbfaces; fac++)
    {
      el0 = face_voisins(fac, 0);
      el1 = face_voisins(fac, 1);
      vol0 = volumes(el0);
      vol1 = volumes(el1);
      kappa_face = (vol0 * kappa_var(el0) + vol1 * kappa_var(el1)) / (vol0 + vol1);
      // cface = (vol0 * c(el0) + vol1 * c(el1)) / (vol0 + vol1);
      // prov_face(fac) = kappa_face * mobilite(cface) * prov_face(fac);
      prov_face(fac) = kappa_face * prov_face(fac); // Changement de Luis: mobilite = 1 de toute façon dans le code initial
    }
 
  // Application solveur masse (/M)
  eq_ns.solv_masse().appliquer(prov_face);
 
  // Div(Kappa*Grad(F))
  opdiv.calculer(prov_face, resu);
  //Cerr <<"div_kappa_grad"<<resu<<finl;
 
  return resu;
}
 
/*! @brief Calcul de Div(alpha*rho*Grad((C)) au centre des elements
 *
 * @param (DoubleTab& div_alpha_gradC) Div(alpha*Grad((C)) au centre des elements
 */
void Source_Con_Phase_field_Binaire::calculer_div_alpha_grad(const DoubleTab& c, DoubleTab& div_alpha_gradC) const
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const Navier_Stokes_std& eq_ns = ref_cast(Navier_Stokes_std, pb_->equation(0));
  const Operateur_Grad& opgrad = eq_ns.operateur_gradient();
  const Operateur_Div& opdiv = eq_ns.operateur_divergence();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
 
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
  prov_face = 0.;
 
  // Grad(C)
  opgrad.calculer(c, prov_face);
  eq_ns.solv_masse().appliquer(prov_face);
  prov_face *= fermeture.get_alpha();
 
  // Div(alpha*Grad(c))
  opdiv.calculer(prov_face, div_alpha_gradC);
  eq_c.solv_masse().appliquer(div_alpha_gradC);
  //Cerr <<"div_alpha_gradC"<<div_alpha_gradC<<finl;
 
}
 
/*! @brief Calcul de Div(alpha*rho*Grad((C)) au centre des elements
 *
 * @param (DoubleTab& div_alpha_rho_gradC) Div(alpha*rho*Grad((C)) au centre des elements
 */
void Source_Con_Phase_field_Binaire::calculer_div_alpha_rho_grad(const DoubleTab& c, DoubleTab& div_alpha_rho_gradC) const
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const Navier_Stokes_phase_field& eq_ns = ref_cast(Navier_Stokes_phase_field, pb_->equation(0));
  const Operateur_Grad& opgrad = eq_ns.operateur_gradient();
  const Operateur_Div& opdiv = eq_ns.operateur_divergence();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
 
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  const DoubleVect& volumes = domaine_VDF.volumes();
  const int ndeb = domaine_VDF.premiere_face_int();
  const int nbfaces = domaine_VDF.nb_faces();
 
  int el0, el1;
  double vol0, vol1;
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
  prov_face = 0.;
 
  // Grad(C)
  opgrad.calculer(c, prov_face);
  eq_ns.solv_masse().appliquer(prov_face);
 
  const DoubleTab& rhoPF = mil.rho().valeurs();
  double rho_face;
 
  if (boussi_ == 1)
    prov_face *= fermeture.get_alpha() * rho0_; // Cas approximation de Boussinesq
  else if (boussi_ == 0)
    {
      for (int fac = ndeb; fac < nbfaces; fac++)
        {
          el0 = face_voisins(fac, 0);
          el1 = face_voisins(fac, 1);
          vol0 = volumes(el0);
          vol1 = volumes(el1);
 
          rho_face = (vol0 * rhoPF(el0) + vol1 * rhoPF(el1)) / (vol0 + vol1);
          prov_face(fac) *= fermeture.get_alpha() * rho_face;
        }
    }
 
  // Div(alpha*rho*Grad(c))
  opdiv.calculer(prov_face, div_alpha_rho_gradC);
  eq_c.solv_masse().appliquer(div_alpha_rho_gradC);
 
}
 
/*! @brief Calcule le premier demi pas de temps dans le cas implicite Calcule le pas de temps dans le cas explicite
 *
 */
void Source_Con_Phase_field_Binaire::premier_demi_dt()
{
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const int nb_elem = domaine_VDF.nb_elem_tot();
 
  Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  DoubleTab& c = eq_c.inconnue().valeurs();
 
  DoubleTab& mutilde = eq_c.set_mutilde_().valeurs();
  DoubleTab& mutilde_demi = eq_c.set_mutilde_demi();
  DoubleTab& c_demi = eq_c.set_c_demi();
 
  Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  Fermeture_Phase_field_base& fermeture = mil.get_fermeture();
 
  // Utilise dans l'ancien modele de Didier Jamet
 
  /* commente par mr264902 car pas utilise dans la version presente
   //DoubleTab& alpha_gradC_carre=eq_c.set_alpha_gradC_carre();
   //calculer_alpha_gradC_carre(alpha_gradC_carre);
   //DoubleTab& div_alpha_rho_gradC=eq_c.set_div_alpha_rho_gradC();
   //
   //calculer_div_alpha_rho_grad(c, div_alpha_rho_gradC);
   */
  DoubleTab& div_alpha_gradC = eq_c.set_div_alpha_gradC();
  calculer_div_alpha_grad(c, div_alpha_gradC);
  // commente par mr264902 car pas utilise dans la version presente
  //DoubleTab& pression_thermo=eq_c.set_pression_thermo();
  //calculer_pression_thermo(pression_thermo);
 
  accr_ = c; // on copie la structure // !!
  c_demi = c; // on copie la structure // !!
 
  accr_ = 0.; // initialise a 0
  c_demi = 0.; // initialise a 0
 
  if (implicitation_ == 1)
    {
      mutilde_demi = c; // on copie la structure // !!
      mutilde_demi = 0.; // initialise a 0
 
      // Pour le JFNK ----------------------------
 
      // Commente par DJ
      //----------------
      //       DoubleTab x1(2*nb_elem);
      //       x1=0.;
      //       if(gmres_==1)
      //         {
 
      //           // Ecriture de x1(0)
      //           for(int n_elem=0;n_elem<nb_elem;n_elem++)
      //             {
      //               x1(n_elem)=c(n_elem);
      //               x1(n_elem+nb_elem)=mutilde(n_elem);
      //             }
      //         }
      //----------------
      // ---------------------------------------------
 
      calculer_u2_elem(u_carre_);
 
      // only once !
      if (!is_matrix_initialised_)
        initialiser_matrice_bloc();
 
      assembler_ancienne_matrice();
 
      // Modifie par DJ
      //---------------
 
      // RESOLUTION JFNK
      // =========================================================
      resolution_jfnk(c, mutilde, c_demi, mutilde_demi);
      // =========================================================
      Cout << "C à la sortie de JFNK = " << c_demi << finl;
      Cout << "mutilde à la sortie de JFNK = " << mutilde_demi << finl;
 
      Cout << "C^n = " << c << finl;
      Cout << "mutilde^n = " << mutilde << finl;
 
      const double theta = 0.6;
      // On stocke les nouveaux c et mutilde
 
      for (int n_elem = 0; n_elem < nb_elem; n_elem++)
        {
          // Commente par DJ
          //----------------
          //               c_demi(n_elem)=x1(n_elem);
          //----------------
          c_demi(n_elem) -= (1 - theta) * c(n_elem);
          c_demi(n_elem) /= theta;
 
          // Commente par DJ
          //----------------
          //               mutilde_demi(n_elem)=x1(n_elem+nb_elem);
          //----------------
          //c=c_demi;
 
          //---------- Calcul de mutilde_demi = dw/dc (c_demi)-div_alpha_gradC(c_demi) ------------
          //calculer_mutilde_demi(mutilde_demi, c_demi);
          //calculer_mutilde(mutilde_demi);
          //mutilde=mutilde_demi;
          mutilde_demi(n_elem) -= (1 - theta) * mutilde(n_elem);
          mutilde_demi(n_elem) /= theta;
        }
 
      Cout << " -------------------------------------- " << finl;
      Cout << "C^n+1 = " << c_demi << finl;
      Cout << "mutilde^n+1 = " << mutilde_demi << finl;
 
      // ATTENTION : A VERIFIER
      //=======================
      //       c_demi.echange_espace_virtuel();
      //       mutilde_demi.echange_espace_virtuel();
      //       c.echange_espace_virtuel();
      //       mutilde.echange_espace_virtuel();
 
      // Mise a jour
 
      if (couplage_ == 0)
        {
          // Si traitement explicite de mutilde dans le cas implicite
          calculer_mutilde(mutilde);
        }
      else if (couplage_ == 1)
        {
          // Si traitement implicite de mutilde dans le cas implicite
          mutilde = mutilde_demi;
        }
      mutilde.echange_espace_virtuel();
 
      // Calcul de l'accroissement entre n et n+1/2
      // Utile pour pouvoir utiliser n'importe quel dt
      accr_ = 0.;
 
    }
  else if (implicitation_ == 0)
    {
      // Cerr<<"Schema explicite"<<finl;
      calculer_u2_elem(u_carre_);
      calculer_mutilde(mutilde);
      mutilde.echange_espace_virtuel();
 
      DoubleTab& prov_elem = prov_elem_;
      if (prov_elem.size() == 0)
        prov_elem = mutilde;
 
      if (fermeture.get_kappa_ind() == 1)
        {
          DoubleTab kappa_var(prov_elem);
          kappa_var = 0.;
          for (int ikappa = 0; ikappa < nb_elem; ikappa++)
            kappa_var(ikappa) = fermeture.call_kappa_func_c(c(ikappa));
          // Div(Kappa*Grad(mutilde))
          div_kappa_grad(mutilde, kappa_var, prov_elem);
        }
      else if (fermeture.get_kappa_ind() == 0)
        {
          // Kappa*Laplacien(mutilde)
          laplacien(mutilde, prov_elem);
          prov_elem *= fermeture.get_kappa();
        }
 
      // Pour equation Allen-Cahn (kappa constant)
      //------------------------------------------
      //   prov_elem=F;
      //   prov_elem*=-kappa;
 
      accr_ += prov_elem;
      //eq_c.solv_masse().appliquer(prov_elem);
      //Cerr<<"masse de la force sur c "<<prov_elem.max_abs()<<finl;
 
      // Utile pour pouvoir utiliser n'importe quel dt - voir Schema_Phase_Field
      c_demi = c;
    }
  else
    {
      Cerr << "Type de schema errone !!" << finl;
    }
}
 
void Source_Con_Phase_field_Binaire::initialiser_matrice_bloc()
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const DoubleTab& c = eq_c.inconnue().valeurs();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
  const int nb_elem = le_dom_VDF_->nb_elem_tot();
 
  kappa_tab_ = c; // on initialise la structure
 
  /*
   * Typage d'une operateur utile
   */
 
  ch_kappa_ou_alpha_ = eq_c.get_mutilde_(); // on copie la structure
  ch_kappa_ou_alpha_->nommer("ch_kappa_ou_alpha");
 
  terme_diffusif_vdf_.associer_eqn(eq_c);
  terme_diffusif_vdf_.associer_diffusivite(ch_kappa_ou_alpha_.valeur());
  EChaine diff( "{ }");
  diff >> terme_diffusif_vdf_;
  terme_diffusif_vdf_.associer_diffusivite_pour_pas_de_temps(ch_kappa_ou_alpha_.valeur());
  terme_diffusif_vdf_->associer(le_dom_VDF_.valeur(), eq_c.domaine_Cl_dis(), eq_c.inconnue());
  terme_diffusif_vdf_->completer();
 
  Cerr << "------------ INITIALISER MATRICE ---------------" << finl;
  Cerr << "le_dom_dCl_dis = " << eq_c.domaine_Cl_dis() << finl;
 
 
  matrice_diffusion_CH_.dimensionner(2, 2);
 
  for (int i = 0; i < 2; i++)
    for (int j = 0; j < 2; j++)
      {
        matrice_diffusion_CH_.get_bloc(i, j).typer("Matrice_Morse");
        Matrice_Morse& mat = ref_cast(Matrice_Morse, matrice_diffusion_CH_.get_bloc(i, j).valeur());
 
        if (i == j) /* Si diagonale => identite */
          {
            Stencil stencil(0, 2);
 
            for (int e = 0; e < nb_elem; e++)
              stencil.append_line(e, e);
 
            Matrix_tools::allocate_morse_matrix(nb_elem, nb_elem, stencil, mat);
            mat.sort_stencil(); // XXX sort tab2 ... au cas ou ...
 
            for (int e = 0; e < nb_elem; e++)
              mat(e, e) = 1.0;
          }
        else
          {
            if (i == 1) /* Matrice de alpha */
              {
                ch_kappa_ou_alpha_->valeurs() = -fermeture.get_alpha();
 
                /*
                * XXX Elie Saikali : attention la y a pas le solveur de masse qui va diviser par V apres ....
                */
                const DoubleVect& vol = le_dom_VDF_->volumes();
                tab_divide_any_shape(ch_kappa_ou_alpha_->valeurs(), vol);
                terme_diffusif_vdf_.associer_diffusivite(ch_kappa_ou_alpha_.valeur());
                terme_diffusif_vdf_->associer_diffusivite(ch_kappa_ou_alpha_.valeur());
                terme_diffusif_vdf_->dimensionner(mat);
                mat.sort_stencil(); // XXX sort tab2 ... au cas ou ...
                terme_diffusif_vdf_->contribuer_a_avec(eq_c.inconnue().valeurs(), mat);
              }
            else /* j == 1 ==> Matrice de kappa */
              {
                terme_diffusif_vdf_.associer_diffusivite(ch_kappa_ou_alpha_.valeur());
                terme_diffusif_vdf_->associer_diffusivite(ch_kappa_ou_alpha_.valeur());
                terme_diffusif_vdf_->dimensionner(mat);
                mat.sort_stencil(); // XXX sort tab2 ... au cas ou ...
                update_bloc_kappa(c);
              }
          }
      }
//
  Cout << "MATRICE BLOCS =========== " << finl;
  matrice_diffusion_CH_.imprimer_formatte(Cout);
 
  is_matrix_initialised_ = true;
 
}
 
void Source_Con_Phase_field_Binaire::update_tab_kappa(const DoubleTab& donne)
{
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
  const int nb_elem = le_dom_VDF_->nb_elem_tot();
 
  if (fermeture.get_kappa_ind() == 0)
    kappa_tab_ = fermeture.get_kappa();
  else
    {
      for (int elem = 0; elem < nb_elem; elem++)
        kappa_tab_(elem) = fermeture.call_kappa_func_c(donne(elem));
    }
}
 
void Source_Con_Phase_field_Binaire::update_bloc_kappa(const DoubleTab& donne)
{
  const auto& fermeture = ref_cast(Milieu_Phase_field, pb_->milieu()).get_fermeture();
 
  if (is_matrix_initialised_ && fermeture.get_type_kappa() ==0)
    return;
 
  update_tab_kappa(donne);
  Matrice_Morse& mat_kappa = ref_cast(Matrice_Morse, matrice_diffusion_CH_.get_bloc(0, 1).valeur());
 
 
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  double dt = eq_c.schema_temps().pas_de_temps();    // on fait un calcul sur un demi pas de temps
  const double theta = 0.6;
 
  ch_kappa_ou_alpha_->valeurs() = kappa_tab_;
  ch_kappa_ou_alpha_->valeurs() *= (theta * dt);
  const DoubleVect& vol = le_dom_VDF_->volumes();
  tab_divide_any_shape(ch_kappa_ou_alpha_->valeurs(), vol);
 
  terme_diffusif_vdf_.associer_diffusivite(ch_kappa_ou_alpha_.valeur());
  terme_diffusif_vdf_->associer_diffusivite(ch_kappa_ou_alpha_.valeur());
  terme_diffusif_vdf_->contribuer_a_avec(eq_c.inconnue().valeurs(), mat_kappa);
}
 
/*! @brief Ancien assemblage de matrice "à la main" : ne prend pas en compte les différents
 * types de schémas numériques (que VDF, pas de VEF) et les différentes CL.
 *
 */
 
void Source_Con_Phase_field_Binaire::assembler_ancienne_matrice()
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const DoubleTab& c = eq_c.inconnue().valeurs();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
 
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  const DoubleTab& positions = domaine_VDF.xp();
  const IntVect& ori = domaine_VDF.orientation();
 
  int compt = 0;
  int compt1 = 0;
  int compt2 = 0;
 
  const int nb_elem = domaine_VDF.nb_elem_tot();
  int nb_compo_;
  int f0;
  int min_tri;
  int old_tri;
  int voisin;
  int nb_faces_au_bord;
  int dimensionnement;
  double dvar2 = 0.;
  double dvarkeep = 0.;
  double dt;
  double valeur_diag = 0.;
  const double theta = 0.6;
  double kappa_interpolee = 0.;
 
  if (mobilite_explicite_ == 1)
    {
      Cout << "Explicit mobility: Building matrix A used in JFNK..." << finl;
      //   Cerr<<"Assemblage de la matrice du point fixe"<<finl;
 
      // Dimension du pb
      nb_compo_ = dimension;
 
      // Forme de la matrice : ( I A )
      //                       ( B I )
 
      // Assemblage de la moitie superieure de la matrice : I et A
 
      // On dimensionne la matrice A et donc en meme temps B (meme nombre de termes non nuls)
      // Pour cela : dimension du vecteur coeff = (2 * dim + 1) * nb_elem - nombre_de_bords
 
      DoubleTab kappa_var(c);
 
      if (fermeture.get_kappa_ind() == 0)
        kappa_var = fermeture.get_kappa();
      else
        {
          for (int elem = 0; elem < nb_elem; elem++)
            kappa_var(elem) = fermeture.call_kappa_func_c(c(elem));
        }
 
      dimensionnement = (2 * nb_compo_ + 1) * nb_elem;
 
      for (int elem = 0; elem < nb_elem; elem++)
        {
          for (int ncomp = 0; ncomp < 2 * nb_compo_; ncomp++)
            {
              f0 = elem_faces(elem, ncomp);
              voisin = face_voisins(f0, 0);
              if (face_voisins(f0, 0) == elem)
                voisin = face_voisins(f0, 1);
 
              if (voisin == -1)
                dimensionnement--;
            }
        }
      // On ajoute le nombre de 1 venant de I, ceci valant pour I et A. Le dimensionnement total est donc le double
      dimensionnement += nb_elem;
      dimensionnement *= 2;
 
      // Allocation des tableaux specifiques a la matrice morse
      ancienne_matrice_.dimensionner(2 * nb_elem, dimensionnement);
      auto& coeff = ancienne_matrice_.get_set_coeff();
      auto& tab2 = ancienne_matrice_.get_set_tab2();
      auto& tab1 = ancienne_matrice_.get_set_tab1();
      //Cerr <<"ancienne_matrice_ = "<<ancienne_matrice_<<finl;
 
      // Boucle sur le nombre d'elements
      for (int elem = 0; elem < nb_elem; elem++)
        {
          valeur_diag = 0.;
          old_tri = -1;
          nb_faces_au_bord = 0;
 
          // On ajoute le 1 de la diagonale de la sous-matrice I du haut a gauche
          coeff(compt) = 1;
          tab2(compt2) = elem;
          compt++;
          compt2++;
          tab1(compt1) = compt2;
          compt1++;
 
          // Calcul du nombre de bords
          for (int ncomp = 0; ncomp < 2 * nb_compo_; ncomp++)
            {
              f0 = elem_faces(elem, ncomp);
              voisin = face_voisins(f0, 0);
              if (face_voisins(f0, 0) == elem)
                voisin = face_voisins(f0, 1);
 
              // On connait le voisin associe a la face en cours. On regarde s'il est au bord
              if (voisin == -1)
                nb_faces_au_bord++;
            }
 
          // Boucle sur les faces
          for (int ncomp = 0; ncomp < 2 * nb_compo_ - nb_faces_au_bord; ncomp++)
            {
 
              min_tri = nb_elem + 1;
              dt = eq_c.schema_temps().pas_de_temps();    // on fait un calcul sur un demi pas de temps
 
              for (int ncomp_tri = 0; ncomp_tri < 2 * nb_compo_; ncomp_tri++)
                {
                  f0 = elem_faces(elem, ncomp_tri);
 
                  // Voisin associe a la face - On s'assure de traiter un voisin et pas l'element resident
                  voisin = face_voisins(f0, 0);
                  if (face_voisins(f0, 0) == elem)
                    voisin = face_voisins(f0, 1);
 
                  // On va compter le nombre de voisins existants (faces non aux bords de l'element elem)
                  // Ceci sert a calculer la contribution au terme relatif a l'element elem
                  // Cas kappa variable : on utilise la moyenne geometrique des valeurs des mobilites
                  // des elements concernes par la face de calcul ("kappa_interpolee")
                  if (voisin != -1 && old_tri == -1)
                    {
                      dvar2 = pow((positions(elem, ori(f0)) - positions(voisin, ori(f0))), 2);
                      if (kappa_moy_ == 2)
                        kappa_interpolee = pow(std::fabs(kappa_var(elem) * kappa_var(voisin)), 0.5);    // Moyenne geometrique
                      else if (kappa_moy_ == 1)
                        {
                          if (kappa_var(elem) == 0 || kappa_var(voisin) == 0)
                            kappa_interpolee = 0;
                          else
                            kappa_interpolee = 2. / (1. / kappa_var(elem) + 1. / kappa_var(voisin));  // Moyenne harmonique
                        }
                      else if (kappa_moy_ == 0)
                        kappa_interpolee = (kappa_var(elem) + kappa_var(voisin)) / 2.;        // Moyenne arithmetique
 
                      valeur_diag += theta * kappa_interpolee * (dt / dvar2);
                    }
 
                  // Si le voisin a un numero plus petit que l'ancien minimum, on ne le prend pas en compte
                  // En effet, le nouveau minimum est un numero d'element qui lui sera superieur
                  // Rq : voisin>old_tri => voisin != -1 donc pour une face au bord le if est false
                  if (voisin > old_tri)
                    {
                      min_tri = std::min(min_tri, voisin);
                      if (min_tri == voisin)
                        dvarkeep = pow((positions(elem, ori(f0)) - positions(voisin, ori(f0))), 2);
                    }
                }
 
              //Cerr <<" min_tri = std::min(min_tri,voisin) if voisin>old_tri ====== "<<min_tri<<finl;
              // Si l'ancien minimum est inferieur a elem et le nouveau superieur, il faut traiter elem juste avant le nouveau voisin
              if (old_tri < elem && min_tri > elem)
                {
                  coeff(compt) = valeur_diag;
                  tab2(compt2) = elem + nb_elem;
                  compt++;
                  compt2++;
                }
 
              // On complete les tableaux avec les valeurs dans les voisins non aux bords du domaine
              if (kappa_moy_ == 2)
                kappa_interpolee = pow(std::fabs(kappa_var(elem) * kappa_var(min_tri)), 0.5);    // Moyenne geometrique
              else if (kappa_moy_ == 1)
                {
                  if (kappa_var(elem) == 0 || kappa_var(min_tri) == 0)
                    kappa_interpolee = 0;
                  else
                    kappa_interpolee = 2. / (1. / kappa_var(elem) + 1. / kappa_var(min_tri));  // Moyenne harmonique
                }
              else if (kappa_moy_ == 0)
                kappa_interpolee = (kappa_var(elem) + kappa_var(min_tri)) / 2.;        // Moyenne arithmetique
 
              coeff(compt) = -theta * kappa_interpolee * (dt / dvarkeep);
              tab2(compt2) = min_tri + nb_elem;
              compt++;
              compt2++;
 
              // Si on est a la fin et que elem n'a pas encore ete mis dans coeff, il faut le faire avant la fin de la boucle sur les faces de l'element
              // Ceci arrive si elem > tous les numeros d'elements de ses voisins
              if (ncomp == 2 * nb_compo_ - nb_faces_au_bord - 1 && min_tri < elem)
                {
                  coeff(compt) = valeur_diag;
                  tab2(compt2) = elem + nb_elem;
                  compt++;
                  compt2++;
                }
 
              // On sauve le numero d'element "minimum"
              old_tri = min_tri;
            }
 
        }
 
      // Assemblage de la moitie superieure de la matrice : B et I
 
      // Boucle sur le nombre d'elements
      for (int elem = 0; elem < nb_elem; elem++)
        {
          valeur_diag = 0;
          old_tri = -1;
          nb_faces_au_bord = 0;
 
          // Calcul du nombre de bords
          for (int ncomp = 0; ncomp < 2 * nb_compo_; ncomp++)
            {
              f0 = elem_faces(elem, ncomp);
              voisin = face_voisins(f0, 0);
              if (face_voisins(f0, 0) == elem)
                voisin = face_voisins(f0, 1);
 
              // On connait le voisin associe a la face en cours. On regarde s'il est au bord
              if (voisin == -1)
                nb_faces_au_bord++;
            }
 
          // Boucle sur les faces
          for (int ncomp = 0; ncomp < 2 * nb_compo_ - nb_faces_au_bord; ncomp++)
            {
 
              min_tri = nb_elem + 1;
 
              for (int ncomp_tri = 0; ncomp_tri < 2 * nb_compo_; ncomp_tri++)
                {
                  f0 = elem_faces(elem, ncomp_tri);
 
                  // Voisin associe a la face - On s'assure de traiter un voisin et pas l'element resident
                  voisin = face_voisins(f0, 0);
                  if (face_voisins(f0, 0) == elem)
                    voisin = face_voisins(f0, 1);
 
                  // On va compter le nombre de voisins existants (faces non aux bords de l'element elem)
                  // Ceci sert a calculer la contribution au terme relatif a l'element elem
                  if (voisin != -1 && old_tri == -1)
                    {
                      dvar2 = pow((positions(elem, ori(f0)) - positions(voisin, ori(f0))), 2);
                      valeur_diag += -fermeture.get_alpha() / dvar2;
                    }
 
                  // Si le voisin a un numero plus petit que l'ancien minimum, on ne le prend pas en compte
                  // En effet, le nouveau minimum est un numero d'element qui lui sera superieur
                  // Rq : voisin>old_tri => voisin != -1 donc pour une face au bord le if est false
                  if (voisin > old_tri)
                    {
                      min_tri = std::min(min_tri, voisin);
                      if (min_tri == voisin)
                        dvarkeep = pow((positions(elem, ori(f0)) - positions(voisin, ori(f0))), 2);
                    }
 
                }
 
              // Si l'ancien minimum est inferieur a elem et le nouveau superieur, il faut traiter elem juste avant le nouveau voisin
              if (old_tri < elem && min_tri > elem)
                {
                  coeff(compt) = valeur_diag;
                  tab2(compt2) = elem;
                  compt++;
                  compt2++;
                  if (old_tri == -1)
                    {
                      // On ajoute un terme a tab1 dans le cas ou l'element diagonal est le premier terme
                      tab1(compt1) = compt2;
                      compt1++;
                      old_tri = elem;
                    }
                }
 
              // On complete les tableaux avec les valeurs dans les voisins non aux bords du domaine
              coeff(compt) = fermeture.get_alpha() / dvarkeep;
              tab2(compt2) = min_tri;
              compt++;
              compt2++;
              if (old_tri == -1)
                {
                  // On ajoute un terme a tab1 dans le cas ou le premier terme est extradiagonal
                  tab1(compt1) = compt2;
                  compt1++;
                }
 
              // Si on est a la fin et que elem n'a pas encore ete mis dans coeff, il faut le faire avant la fin de la boucle sur les faces de l'element
              // Ceci arrive si elem > tous les numeros d'elements de ses voisins
              // Bien sur dans ce cas la, l'element diagonal est forcement le dernier mais pas le premier - il n'existe pas de maillage avec des mailles seules sans voisins
              if (ncomp == 2 * nb_compo_ - nb_faces_au_bord - 1 && min_tri < elem)
                {
                  coeff(compt) = valeur_diag;
                  tab2(compt2) = elem;
                  compt++;
                  compt2++;
                }
 
              // On sauve le numero d'element "minimum"
              old_tri = min_tri;
 
            }
 
          // On ajoute le 1 de la diagonale de la sous-matrice I du bas a droite
          coeff(compt) = 1;
          tab2(compt2) = elem + nb_elem;
          compt++;
          compt2++;
 
        }
 
      // cf Mat_Morse.h/.cpp, tab1() et tab2() sont a utiliser au sens Fortran
      // EXPLICATION :
      // coeff stocke les valeurs des coefficients
      // tab2 stocke les numeros de colonnes dans la matrice de ces coefficients, au sens FORTRAN (i-e 1 <= tab2[i] <= n)
      // tab1 stocke le rang de l'element de tab2() pour lequel on change de ligne. Ainsi pour tab1[j], on est a la i-eme ligne avec :
      //      ( tab1[i] <= j < tab1[i+1] ) qui implique qu'on considere toujours la i-eme ligne,
      //      et des que cela n'est plus respecte, le passage a la ligne suivante est effectue
 
      // De par l'algorithme, tab1() est deja au sens FORTRAN. Reste a le faire pour tab2()
      // De plus on doit mettre le dernier terme de tab1() a dimensionnement+1
      tab2 += 1;
      tab1(compt1) = dimensionnement + 1;
      compt1 += 1;
 
      if (compt != dimensionnement)
        {
          Cerr << "Erreur lors du calcul de la matrice du point fixe : nombre d'elements non nuls calcules different du nombre d'elements non nuls prevus" << finl;
          Process::exit();
        }
    }
  else if (mobilite_explicite_ == 0)
    {
      Cerr << "Implicit mobility: the matrix A used in JFNK is not implemented." << finl;
      ancienne_matrice_.clean();
    }
  else
    {
      Cerr << "STOP: argument mobilite_explicite cannot be different than 0 or 1." << finl;
      Cerr << "==> Must have an error in the code." << finl;
    }
 
//  Cerr << "---------------------- Ancienne matrice = " << finl;
//  ancienne_matrice_.imprimer_formatte(Cerr);
 
}
 
/*! @brief Construire le residu du JFNK
 *
 */
 
void Source_Con_Phase_field_Binaire::calculer_residu_jfnk(const DoubleTab& c, DoubleTab& v0, const DoubleTab& x1)
{
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
  const int nb_elem = le_dom_VDF_->nb_elem_tot();
 
  DoubleVect term_cin(nb_elem);
  // Ceci correspond a u2/2*drhodc=mutilde_dyn-mutilde
 
  DoubleVect second_membre(2 * nb_elem);
  DoubleTab Ax1(2 * nb_elem);
 
  // Assemblage du second membre
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      second_membre(n_elem) = c(n_elem);
 
      term_cin(n_elem) = 0.;
      if (mutype_ == 1)
        term_cin(n_elem) += (0.5 * u_carre_(n_elem)) * drhodc(n_elem);
 
      second_membre(n_elem + nb_elem) = term_cin(n_elem) + fermeture.get_beta() * fermeture.call_dWdc(x1(n_elem));
    }
 
  // Calcul du produit matrice / vecteur utilise
  matrice_diffusion_CH_.multvect_(x1, Ax1);
  // ancienne_matrice_.multvect_(x1, Ax1);
 
  // Modifie par DJ
  //---------------
  DoubleTab Ax1_c(c);
  DoubleTab Ax1_mutilde(c);
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      Ax1_c(n_elem) = Ax1(n_elem);
      Ax1_mutilde(n_elem) = Ax1(n_elem + nb_elem);
    }
 
  Ax1_c.echange_espace_virtuel();
  Ax1_mutilde.echange_espace_virtuel();
 
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      Ax1(n_elem) = Ax1_c(n_elem);
      Ax1(n_elem + nb_elem) = Ax1_mutilde(n_elem);
    }
 
  // Calcul de v0 = [A(xn) xn - bn]
  for (int n_elem = 0; n_elem < 2 * nb_elem; n_elem++)
    v0(n_elem) = (Ax1(n_elem) - second_membre(n_elem));
 
}
 
/*! @brief Construire la jacobienne multipliée par un vecteur Y
 *
 */
void Source_Con_Phase_field_Binaire::jacobian_vect(const DoubleTab& c, const DoubleTab& v0, const DoubleTab& x1, DoubleTab& v1)
{
  const double delta = 1.e-5;
 
  DoubleTab v2(v1);
  DoubleTab x1_(v0);
 
  v1 = 0.;
 
  x1_ = v0;
  x1_ *= delta;
  x1_ += x1;
 
  // Construction de la differentielle (dans la direction v0)
  calculer_residu_jfnk(c, v1, x1_);
  calculer_residu_jfnk(c, v2, x1);
 
  v1 -= v2;
  v1 /= delta;
}
 
#ifdef PETSCKSP_H
// Convertir DoubleTab -> PETSc Vec
void tabToVecOld(const DoubleTab& tab, Vec& f)
{
  PetscScalar *f_array;
  PetscInt size;
  VecGetArray(f, &f_array);
  VecGetSize(f, &size);
  assert(size == tab.dimension_tot(0));
  for (int i = 0; i < tab.dimension_tot(0); i++)
    f_array[i] = tab(i);
  VecRestoreArray(f, &f_array);
}
 
// Convertir PETSc Vec -> DoubleTab
void vecToTabOld(const Vec& v, DoubleTab& result)
{
  PetscInt size;
  const PetscScalar *array;
  VecGetSize(v, &size);
  VecGetArrayRead(v, &array);
  assert((int)size == result.dimension_tot(0));
  for (PetscInt i = 0; i < size; ++i)
    result((int)i) = array[i];
  VecRestoreArrayRead(v, &array);
}
 
/*! @brief ForFormFunctionOld de PETSc : calcule le résidu du JFNK
 *
 */
//PetscErrorCode FormFunctionOld(SNES snes, Vec x, Vec f, void* ctx)
//{
//  Equation_base* eqn = static_cast<Equation_base*>(ctx);
//  Source_Con_Phase_field_Binaire& source = ref_cast(Source_Con_Phase_field_Binaire, eqn->sources().dernier().valeur());
//
//  const int nb_elem = eqn->inconnue().valeurs().dimension_tot(0);
//  DoubleTab x_tab(2 * nb_elem);
//  vecToTabOld(x, x_tab);
//  DoubleTab c(nb_elem), mutilde(nb_elem);
//  DoubleTab res(2 * nb_elem);
//
//  for (int i = 0; i < nb_elem; ++i)
//    {
//      c(i) = x_tab(i);
//      mutilde(i) = x_tab(i + nb_elem);
//    }
//
//  source.calculer_residu_jfnk(c, res, x_tab);
//
//  tabToVecOld(res, f);
//
//  // Cerr << "F(X) = " << res << finl;
//
//  return 0;
//}
 
PetscErrorCode FormFunctionOld(SNES snes, Vec x, Vec f, void* ctx)
{
  Equation_base* eqn = static_cast<Equation_base*>(ctx);
  Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, *eqn);
  const DoubleTab& C = eq_c.inconnue().valeurs();
 
  Source_Con_Phase_field_Binaire& source = ref_cast(Source_Con_Phase_field_Binaire, eqn->sources().dernier().valeur());
  const int nb_elem = C.dimension_tot(0);
  DoubleTab x1(2*nb_elem);
  vecToTabOld(x,x1);
 
//  Cerr << "&&&&& D'abord Vec X = " << finl;
//  VecView(x, PETSC_VIEWER_STDOUT_WORLD);
 
//  Cerr << "Et X1 = " << x1 << finl;
 
  DoubleTrav Ax1(x1);
 
//  Cerr << "Ax1 avant calcul = " << Ax1 << finl;
 
  source.get_matrice_diffusion_CH().multvect_(x1, Ax1);
 
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, eq_c.milieu());
  const auto& fermeture = mil.get_fermeture();
  // Assemblage du second membre
  DoubleTrav second_membre(x1);
 
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      second_membre(n_elem) = C(n_elem);
      second_membre(n_elem + nb_elem) = fermeture.get_beta() * fermeture.call_dWdc(x1(n_elem));
    }
//  Cerr << "second_membre = " << second_membre << finl;
 
  second_membre -= Ax1;
 
  tabToVecOld(second_membre, f);
 
  //  Cerr << "F =" << finl;
//  VecView(f, PETSC_VIEWER_STDOUT_WORLD);
 
//  Cerr << "F(X) = " << second_membre << finl;
 
  return 0;
}
 
PetscErrorCode MySNESMonitor(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);
 
  Cout << "[SNES] Iteration " << its << " : ||F(x)|| = " << norm << finl;
 
  if (!Process::me())
    {
      SFichier fichier_residu("residu.txt", ios::app);
      fichier_residu << temps_exec << "\t" << eqn->schema_temps().temps_courant() << "\t" << (double)norm << finl;
    }
 
 
//  Vec x;
//  SNESGetSolution(snes, &x);
//
//  Vec r;
//  SNESGetFunction(snes, &r, nullptr, nullptr);  // récupère F(x) courant
// Facultatif : afficher x
//  Cerr << "Solution Iteration " << its << " : x = " << finl;
//  VecView(x, PETSC_VIEWER_STDOUT_WORLD);
//  Cerr << "Résidu Iteration " << its << " : r = " << finl;
//  VecView(r, PETSC_VIEWER_STDOUT_WORLD);
 
  return 0;
}
 
PetscErrorCode LuisKSPMonitor(KSP ksp, PetscInt its, PetscReal rnorm, void* ctx)
{
  Cout << "  [KSP] Iteration " << its << " : residu = " << rnorm << finl;
  //  Vec x;
  //  KSPBuildSolution(ksp, nullptr, &x);  // récupère x courant (solution approx.)
 
//  VecView(x, PETSC_VIEWER_STDOUT_WORLD);  // affiche le contenu de x
 
  return 0;
}
 
int Source_Con_Phase_field_Binaire::petsc_snes(const DoubleTab& c, const DoubleTab& mutilde, DoubleTab& c_demi, DoubleTab& mutilde_demi)
{
 
  const int ns = 2 * c.size_totale();
  const int nb_elem_tot = c.size_totale();
  DoubleTrav x1(ns); // inconnue (vector de C et mutilde)
  DoubleTrav v0(ns); // tableau TRUST du résidu
 
  for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
    {
      x1(n_elem) = c(n_elem);
      x1(n_elem + nb_elem_tot) = mutilde(n_elem);
    }
 
  // Variables SNES ########################################
  SNES            snes;     // Contexte SNES
  Vec             x, r;     // Vecteurs : solution, résidu
  KSP ksp;
  Mat J;
  PC pc;
  // #######################################################
 
  // INITIALISATION DU VECTEUR D'INCONNUES X
  VecCreate(PETSC_COMM_WORLD, &x);
  VecSetSizes(x, PETSC_DECIDE, ns);
  VecSetFromOptions(x);
 
  // INITIAL GUESS of X
  tabToVecOld(x1, x);
 
  VecDuplicate(x, &r);
 
  // Vérification tailles
  PetscInt x_size, r_size;
  VecGetSize(x, &x_size);
  VecGetSize(r, &r_size);
  assert(x_size == ns);
  assert(r_size == ns);
 
  // ---------------------- Compute real values of tolerances ---------------------------
  // Appel explicite de ta fonction de résidu
  FormFunctionOld(nullptr, x, r, &equation());
 
  // 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 FormFunctionOld to the SNES context ============
  SNESSetFunction(snes, r, FormFunctionOld, &equation());
  // ======================================================================
 
  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);
 
  // 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, LuisKSPMonitor, nullptr, nullptr);
  SNESMonitorSet(snes, MySNESMonitor, &equation(), 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_WORLD);
 
  // ############################### SOLVING ################################################
  SNESSolve(snes, nullptr, x);                                                          //  #
  // ########################################################################################
 
  PetscInt nfunc;
  SNESGetNumberFunctionEvals(snes, &nfunc);
 
  if (!Process::me())
    {
      SFichier fichier_eval("eval.txt", ios::app);
      fichier_eval << equation().schema_temps().temps_courant() << "\t" << nfunc << finl;
    }
 
 
  const char *reasonstr;
  SNESGetConvergedReasonString(snes, &reasonstr);
  Cout << "SNES converged for this reason: " << reasonstr << finl;
 
  DoubleTrav solution(ns); // inconnue (vector de C et mutilde)
 
  vecToTabOld(x,solution);
 
  for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
    {
      c_demi(n_elem) = solution(n_elem);
      mutilde_demi(n_elem) = solution(n_elem + nb_elem_tot);
    }
 
  // Nettoyage
  VecDestroy(&x);
  VecDestroy(&r);
  MatDestroy(&J);
  SNESDestroy(&snes);
 
  return 1;
}
#endif
 
/*! @brief Algorithme JFNK (ancienne version de D. Jamet et nouvelle version de PETSc)
 *
 * @param (mutilde) inverse A x = b
 * @return (int)
 *
 * Solves Ax = b by GMRES with nit iterations (reinitializations) and nkr Krilov vectors.
 * Copied from Sch_Crank_Nicholson.cpp (scalar)//
 *
 */
int Source_Con_Phase_field_Binaire::resolution_jfnk(const DoubleTab& c, const DoubleTab& mutilde, DoubleTab& c_demi, DoubleTab& mutilde_demi)
{
#ifdef PETSCKSP_H
  if (with_petsc_ == 1)
    return petsc_snes(c, mutilde, c_demi, mutilde_demi);
  else if (with_petsc_ == 0)
    {
      int i, j, nk, i0, im, it, ii;
      double tem = 1., res, ccos, ssin;
      const int ns = 2 * c.size_totale();
      const int nb_elem_tot = c.size_totale();
      double temps_exec = 0.;
      int nfunc=0;
//      statistiques().set_debug_level(9);
 
//      Cout << "MATRICE BLOCS M =========== " << finl;
//      matrice_diffusion_CH_.get_bloc(0,1)->imprimer_formatte(Cout);
//      Cout << "MATRICE BLOCS K =========== " << finl;
//      matrice_diffusion_CH_.get_bloc(1,0)->imprimer_formatte(Cout);
 
      // Ajoute par DJ
      //--------------
      DoubleTrav x1(ns);
      DoubleTrav c_update(nb_elem_tot);
 
      for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
        {
          x1(n_elem) = c(n_elem);
          x1(n_elem + nb_elem_tot) = mutilde(n_elem);
        }
 
      Cout << "Vecteur x1 = " << x1 << finl;
 
      DoubleTab v(ns, nkr_);                         // Krilov vectors
      DoubleTab h(nkr_ + 1, nkr_);                // Heisenberg maatrix of coefficients
      DoubleVect r(nkr_ + 1);
      DoubleTab v0(x1);
      DoubleTab v1(x1);
 
      // Initialisation
      v = 0.;
      v0 = 0.;
      v1 = 0.;
 
      // v0 = -1.*calculer_residu_jfnk(eqn, x1); // DJ
 
      calculer_residu_jfnk(c, v0, x1);
      nfunc++;
      v0 *= -1.; // b - AX
 
      Cout << "Residu v0 = " << v0 << finl;
 
      res = 0.;
      // Modifie par DJ
      //---------------
      for (ii = 0; ii < c.size(); ii++) // Partie c
        res += v0(ii) * v0(ii);
      for (ii = c.size_totale(); ii < c.size_totale() + c.size(); ii++) // Partie mutilde
        res += v0(ii) * v0(ii);
 
      res = mp_sum(res);
      res = sqrt(res);
 
      // Cout << "Premier vecteur v0 = " << v0 << finl;
      Cout << "Norme du premier résidu = " << res << finl;
 
      if (res < rtol_)
        return 0; // nothing to do
 
      rtol_ = (rtol_ < res * atol_) ? res * atol_ : rtol_;
      rtol_ = (rtol_ < stol_) ? rtol_ : stol_;
 
      Cout << "Source Concentration Phase Field - JFNK" << finl;
      Cout << "Stopping rule scalar : " << rtol_ << finl;
 
      // iterations
      for (it = 0; it < maxit_newton_; it++)
        {
          nk = nkr_;
 
          //...Orthogonalisation of Arnoldi
          v0 /= res;
          r = 0.;
          r[0] = res;
          h = 0.;
 
          for (j = 0; j < nkr_; j++)
            {
              for (ii = 0; ii < ns; ii++)
                v(ii, j) = v0(ii);
 
              jacobian_vect(c, v0, x1, v1);
              nfunc++;
              nfunc++;
              v0 = v1;
 
              // Modifie par DJ
              //---------------
              for (i = 0; i <= j; i++)
                {
                  for (ii = 0; ii < c.size(); ii++)
                    h(i, j) += v0(ii) * v(ii, i);
                  for (ii = c.size_totale(); ii < c.size_totale() + c.size(); ii++)
                    h(i, j) += v0(ii) * v(ii, i);
                  h(i, j) = mp_sum(h(i, j));
 
                  for (ii = 0; ii < ns; ii++)
                    v0(ii) -= h(i, j) * v(ii, i);
                }
              // tem = sqrt(v0 * v0);
              tem = 0.;
              // Modifie par DJ
              //---------------
              for (ii = 0; ii < c.size(); ii++)
                tem += v0(ii) * v0(ii);
              for (ii = c.size_totale(); ii < c.size_totale() + c.size(); ii++)
                tem += v0(ii) * v0(ii);
 
              tem = mp_sum(tem);
              tem = sqrt(tem);
              h(j + 1, j) = tem;
              if (tem < rtol_)
                {
                  nk = j + 1;
                  goto l5;
                }
              v0 /= tem;
            }
 
          //...Triangularisation QR de la matrice de Hessenberg
l5:
          for (i = 0; i < nk; i++)
            {
              im = i + 1;
              tem = 1. / sqrt(h(i, i) * h(i, i) + h(im, i) * h(im, i));
              ccos = h(i, i) * tem;
              ssin = -h(im, i) * tem;
              for (j = i; j < nk; j++)
                {
                  tem = h(i, j);
                  h(i, j) = ccos * tem - ssin * h(im, j);
                  h(im, j) = ssin * tem + ccos * h(im, j);
                }
              r[im] = ssin * r[i];
              r[i] *= ccos;
            }
 
          //...Solution of linear system
          for (i = nk - 1; i >= 0; i--)
            {
              r[i] /= h(i, i);
              for (i0 = i - 1; i0 >= 0; i0--)
                r[i0] -= h(i0, i) * r[i];
            }
          for (i = 0; i < nk; i++)
            for (ii = 0; ii < ns; ii++)
              x1(ii) += r[i] * v(ii, i);
 
          if (mobilite_explicite_ == 0) /* implicite */
            {
              for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
                c_update(n_elem) = x1(n_elem);
 
              update_bloc_kappa(c_update);
            }
 
          calculer_residu_jfnk(c, v0, x1);
          nfunc++;
          v0 *= -1.;
 
          res = 0.;
          // Modifie par DJ
          //---------------
          for (ii = 0; ii < c.size(); ii++)
            res += v0(ii) * v0(ii);
          for (ii = c.size_totale(); ii < c.size_totale() + c.size(); ii++)
            res += v0(ii) * v0(ii);
 
          res = mp_sum(res);
          res = sqrt(res);
 
          Cout << " - At it = " << it << ", residu scalar = " << res; // << finl;
          // Cout << " and vector R = " << v0 << finl;
 
          temps_exec = statistics().get_time_since_last_open(STD_COUNTERS::total_execution_time);
 
          if (!Process::me())
            {
              SFichier fichier_residu("residu.txt", ios::app);
              fichier_residu << temps_exec << "\t" << equation().schema_temps().temps_courant() << "\t" << res << finl;
            }
 
          if (res < rtol_)
            {
              // Ajoute par DJ
              //--------------
              for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
                {
                  c_demi(n_elem) = x1(n_elem);
                  mutilde_demi(n_elem) = x1(n_elem + nb_elem_tot);
                }
              // L'echange espace virtuel est fait dans premier_demi_dt()
              Cout << "Number of iterations to reach convergence : " << it + 1 << finl;
              Cout << "" << finl;
 
              return it;
            }
          else if (it == maxit_newton_ - 1)
            {
              // Ajoute par DJ
              //--------------
              // On fait le choix de mettre a jour c_demi et mutilde_demi meme s'il n'y a pas convergence...
              for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
                {
                  c_demi(n_elem) = x1(n_elem);
                  mutilde_demi(n_elem) = x1(n_elem + nb_elem_tot);
                }
 
              Cout << "Stopped before convergence" << finl;
              Cout << "" << finl;
 
            }
        }
      if (!Process::me())
        {
          SFichier fichier_eval("eval.txt", ios::app);
          fichier_eval << equation().schema_temps().temps_courant() << "\t" << nfunc << finl;
        }
    }
  else
    {
      Cerr << "Error: Phase_field - Bad input for choosing between resolution with PETSc or not." << finl;
      Process::exit();
    }
#endif
  return -1;
 
}
 
/*! @brief Calcul de mutilde au centre des elements
 *
 * @param (mutilde) mutilde au centre des elements
 */
void Source_Con_Phase_field_Binaire::calculer_mutilde(DoubleTab& mutilde) const
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const DoubleTab& c = eq_c.inconnue().valeurs();
  const DoubleTab& div_alpha_gradC = eq_c.get_div_alpha_gradC();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
 
  // Calcul de mutilde
 
  mutilde = div_alpha_gradC;
  mutilde *= -1.;
 
  //Cerr << "mutilde  = -div_alpha_gradC " << mutilde << finl;
 
  const int taille = mutilde.size();
  for (int i = 0; i < taille; i++)
    {
      mutilde(i) += fermeture.get_beta() * fermeture.call_dWdc(c(i));
      if (mutype_ == 1) //avec Ec
        mutilde(i) += (0.5 * u_carre_(i)) * drhodc(i);
    }
 
// L'espace virtuel n'est pas a jour
  assert_invalide_items_non_calcules(mutilde);
 
}
 
void Source_Con_Phase_field_Binaire::calculer_mutilde_demi(DoubleTab& mutilde_demi, DoubleTab& c_demi) const
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const Navier_Stokes_std& eq_ns = ref_cast(Navier_Stokes_std, pb_->equation(0));
  const Operateur_Grad& opgrad = eq_ns.operateur_gradient();
  const Operateur_Div& opdiv = eq_ns.operateur_divergence();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
 
  mutilde_demi = 0.;
 
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
  prov_face = 0.;
 
  // Grad(C)
  opgrad.calculer(c_demi, prov_face);
  eq_ns.solv_masse().appliquer(prov_face);
  prov_face *= fermeture.get_alpha();
 
  // Div(alpha*Grad(c))
  opdiv.calculer(prov_face, mutilde_demi);
  eq_c.solv_masse().appliquer(mutilde_demi);
  //Cerr <<"div_alpha_gradC"<<div_alpha_gradC<<finl;
 
  mutilde_demi *= -1.;
 
  const int taille = mutilde_demi.size();
  for (int i = 0; i < taille; i++)
    {
      mutilde_demi(i) += fermeture.get_beta() * fermeture.call_dWdc(c_demi(i));
      if (mutype_ == 1) //avec Ec
        {
          mutilde_demi(i) += (0.5 * u_carre_(i)) * drhodc(i);
        }
    }
 
// L'espace virtuel n'est pas a jour
  assert_invalide_items_non_calcules(mutilde_demi);
}