Source_Con_Phase_field_Naire.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_Qdm_VDF_Phase_field.h>
#include <Navier_Stokes_phase_field.h>
#include <Source_Con_Phase_field.h>
#include <Source_Con_Phase_field_Naire.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 <Domaine_VDF.h>
#include <Milieu_base.h>
#include <SolveurSys.h>
#include <Dimension.h>
#include <EChaine.h>
#include <Param.h>
 
Implemente_instanciable(Source_Con_Phase_field_Naire,"Source_Con_Phase_field_Naire_VDF_P0_VDF",Source_Con_Phase_field);
 
// XD Source_Con_Phase_field_Naire_VDF_P0_VDF Source_Con_Phase_field Source_Con_Phase_field_Naire_VDF_P0_VDF INHERITS_BRACE Source for phase field
 
Sortie& Source_Con_Phase_field_Naire::printOn(Sortie& s) const
{
  return s << que_suis_je();
}
 
Entree& Source_Con_Phase_field_Naire::readOn(Entree& is)
{
  Cerr << "Source_Con_Phase_field_Naire::readOn" << finl;
 
  Source_Con_Phase_field::readOn(is);
  return is;
}
 
DoubleTab& Source_Con_Phase_field_Naire::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();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const int nb_comp = mil.get_fermeture().get_nb_constituants();
 
  resu = 0.;
 
  // Dans cette partie, le laplacien s'obtient comme div(grad F) ou F est un DoubleTab avec nb_comp colonnes.
  // On introduit des fonctions temporaires (temp) pour le calcul de chaque colonne avant de les regrouper dans un DoubleTab a la fin
  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)
  DoubleTab temp_resu(resu.dimension_tot(0), 1);
  temp_resu = 0;
  DoubleTab temp_F(F.dimension_tot(0), 1);
  for (int j = 0; j < nb_comp; j++)
    {
      prov_face = 0.;
      for (int i = 0; i < temp_F.dimension(0); i++)
        temp_F(i, 0) = F(i, j);
      opgrad.calculer(temp_F, prov_face);
 
      // Application solveur masse
      eq_ns.solv_masse().appliquer(prov_face);
 
      // Laplacien = Div(Grad(F))
      opdiv.calculer(prov_face, temp_resu);
      for (int i = 0; i < temp_resu.dimension(0); i++)
        resu(i, j) = temp_resu(i, 0);
    }
 
  return resu;
}
 
DoubleTab& Source_Con_Phase_field_Naire::div_kappa_grad(const DoubleTab& F, const DoubleTab& kappa_var, DoubleTab& resu) const
{
  // NOTE: Le tableau F en entrée, dont on calcule le gradient, a la forme suivante :
  // Deux colonnes, la première porte sur les éléments du volume de simulation, la deuxième porte sur le nombre d'équations de Cahn-Hilliard
 
  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 Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const int nb_comp = mil.get_fermeture().get_nb_constituants();
 
  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)
  DoubleTab temp_resu(resu.dimension_tot(0), 1);
  temp_resu = 0;
  DoubleTab temp_F(F.dimension_tot(0), 1);
  DoubleTab temp_prov_face(prov_face.dimension_tot(0), F.line_size());
  for (int j = 0; j < nb_comp; j++)
    {
      for (int i = 0; i < temp_F.dimension(0); i++)
        temp_F(i, 0) = F(i, j);
      opgrad.calculer(temp_F, prov_face);
 
      for (int i = 0; i < prov_face.dimension(0); i++)
        temp_prov_face(i, j) = prov_face(i, 0);
    }
 
  // Determine kappa_face et cface (faces interieures)
  int ndeb = domaine_VDF.premiere_face_int();
  int nbfaces = domaine_VDF.nb_faces();
  int el0, el1;
  double vol0, vol1;
  DoubleTab kappa_face(prov_face.dimension(0), nb_comp * nb_comp);
  //DoubleTab cface(prov_face.dimension(0),nb_comp);
  for (int j = 0; j < nb_comp * nb_comp; j++)
    {
      for (int fac = ndeb; fac < nbfaces; fac++)
        {
          el0 = face_voisins(fac, 0);
          el1 = face_voisins(fac, 1);
          vol0 = volumes(el0);
          vol1 = volumes(el1);
          // on calcule kappa_face comme la moyenne des kappa des voisins
          kappa_face(fac, j) = (vol0 * kappa_var(el0, j) + vol1 * kappa_var(el1, j)) / (vol0 + vol1);
        }
    }
 
  // kappa*Grad(F) - sur les faces
  DoubleTab temp_prov_face2 = temp_prov_face;
  temp_prov_face2 = 0;
  for (int j = 0; j < nb_comp; j++)
    {
      for (int i = 0; i < prov_face.dimension(0); i++)
        {
          prov_face(i, 0) = temp_prov_face(i, j);
          for (int k = 0; k < nb_comp; k++)
            temp_prov_face2(i, k) += prov_face(i, 0) * kappa_face(i, j + k * nb_comp);
        }
    }
 
  // Application solveur masse (/M)
  eq_ns.solv_masse().appliquer(temp_prov_face2);
 
  // Div(Kappa*Grad(F))
  for (int j = 0; j < nb_comp; j++)
    {
      for (int i = 0; i < prov_face.dimension(0); i++)
        prov_face(i, 0) = temp_prov_face2(i, j);
      opdiv.calculer(prov_face, temp_resu);
      for (int i = 0; i < temp_resu.dimension(0); i++)
        resu(i, j) = temp_resu(i, 0);
    }
 
  return resu;
}
 
/*! @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_Naire::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_gradC(div_alpha_rho_gradC);
   */
  DoubleTab& div_alpha_gradC = eq_c.set_div_alpha_gradC();
  // NOTE : La taille de div_alpha_gradC est la suivante :
  // Deux colonnes, une sur le nombre d'éléments et une sur le nombre d'équations de Cahn-Hilliard
  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)
    {
      Matrice_Morse matrice_diffusion_CH;
      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_);
      assembler_matrice_mobilite_explicite(matrice_diffusion_CH);
 
      // Modifie par DJ
      //---------------
      //           resolution_jfnk(c, mutilde, matrice_diffusion_CH, x1);
 
      if (getenv("TRUST_GMRES"))
        {
          // Creation  d'un solveur et typage:
          SolveurSys solveur_;
          //Nom str = "Petsc Gmres { precond diag { } seuil 1.e-6 impr }";
          Nom str = "Petsc Cholesky { impr }";
 
          EChaine chl(str);
          chl >> solveur_;
 
          //DoubleVect x,b;
          // Remplissage guess x et second membre
 
          //const int ns = 2*c.size_totale();
          int nb_elem_tot = c.size_totale();
          DoubleTab X_c(c.size_totale());
          DoubleTab X_mutilde(mutilde.size_totale());
          DoubleTab x1(nb_elem, 2);
          x1 = 0.;
          domaine_VDF.domaine().creer_tableau_elements(x1, RESIZE_OPTIONS::NOCOPY_NOINIT);
          double *x1_addr = x1.addr();
          for (int j = 0; j < nb_equation_CH_; j++)
            for (int i = 0; i < nb_elem; i++)
              {
                X_c(i + (j * nb_elem)) = c(i, j);
                X_mutilde(i + (j * nb_elem)) = mutilde(i, j);
              }
 
          for (int n_elem = 0; n_elem < c.size_totale(); n_elem++)
            {
              x1_addr[n_elem] = X_c(n_elem);
              x1_addr[n_elem + c.size_totale()] = X_mutilde(n_elem);
            }
 
          DoubleTab sm_c(nb_elem_tot);
          DoubleTab sm_mutilde(nb_elem_tot);
          DoubleTab second_membre(nb_elem_tot, 2);
          domaine_VDF.domaine().creer_tableau_elements(second_membre, RESIZE_OPTIONS::NOCOPY_NOINIT);
 
          // Assemblage du second membre
          const DoubleVect& betaMatrix = fermeture.get_betaMatrix();
          double *sm_addr = second_membre.addr();
          for (int j = 0; j < nb_equation_CH_; j++)
            for (int i = 0; i < nb_elem; i++)
              {
                sm_c(i + (j * nb_elem)) = c(i, j);
 
                sm_mutilde(i + (j * nb_elem)) = betaMatrix(j) * fermeture.call_dWdc(c(i, j));
              }
          for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
            {
              sm_addr[n_elem] = sm_c(n_elem);
              sm_addr[n_elem + nb_elem] = sm_mutilde(n_elem);
            }
 
          // Resolution:
          solveur_->reinit();
          solveur_.resoudre_systeme(matrice_diffusion_CH, second_membre, x1);
 
          // remplissage des resultats dans mutilde demi et cdemi
 
          for (int n_elem = 0; n_elem < nb_elem; n_elem++)
            {
              c_demi(n_elem) = x1_addr[n_elem];
              mutilde_demi(n_elem) = x1_addr[n_elem + nb_elem];
            }
          c_demi.echange_espace_virtuel();
          mutilde_demi.echange_espace_virtuel();
        }
      else
        resolution_jfnk(c, mutilde, matrice_diffusion_CH, c_demi, mutilde_demi);
 
      const double theta = 0.6;
      // On stocke les nouveaux c et mutilde
 
      //Cerr<<"c_demi apres resolution_jfnk"<<c_demi<<finl;
      //Cerr<<"c apres resolution_jfnk"<<c<<finl;
 
      for (int j = 0; j < nb_equation_CH_; j++)
        {
          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, j) -= (1 - theta) * c(n_elem, j);
              c_demi(n_elem, j) /= theta;
 
              // Commente par DJ
              //----------------
              //               mutilde_demi(n_elem)=x1(n_elem+nb_elem);
              //----------------
              //c=c_demi;
              //calculer_mutilde_demi(mutilde_demi, c_demi);
              //calculer_mutilde(mutilde_demi);
              //mutilde=mutilde_demi;
              mutilde_demi(n_elem, j) -= (1 - theta) * mutilde(n_elem, j);
              mutilde_demi(n_elem, j) /= theta;
            }
        }
 
      c_demi.echange_espace_virtuel();
      c.echange_espace_virtuel();
      mutilde.echange_espace_virtuel();
 
      // 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.;
      //Cerr<<"c_demi avant } fin implicitation ==1"<<c_demi<<finl;
 
    }
  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_type_kappa_auto_diffusion() == 0)
        {
          DoubleTab temp_mutilde(mutilde.dimension(0), 1);
          temp_mutilde = 0;
          DoubleTab temp_prov_elem = temp_mutilde;
 
          // laplacien(mutilde)
          laplacien(mutilde, prov_elem);
          // Cerr << "laplacien mutilde " << prov_elem<<finl;
 
          //kappa*laplacien(mutilde)
          DoubleTab temp_prov_elem2 = prov_elem;
          temp_prov_elem2 = 0;
          const int nb_comp = fermeture.get_nb_constituants();
          for (int j = 0; j < prov_elem.line_size(); j++)
            {
              for (int i = 0; i < prov_elem.dimension(0); i++)
                {
                  temp_prov_elem(i, 0) = prov_elem(i, j);
                  for (int k = 0; k < nb_comp; k++)
                    {
                      temp_prov_elem2(i, k) += temp_prov_elem(i, 0) * fermeture.get_kappaMatrix()(j + k * nb_comp);
                    }
                }
            }
          prov_elem = temp_prov_elem2;
        }
      else if (fermeture.get_type_kappa_auto_diffusion() == 1)
        {
          DoubleTab kappaMatrix_variable(prov_elem.dimension(0), c.line_size() * c.line_size()); //or (nb_elem,nb_comp*nb_comp)
          kappaMatrix_variable = fermeture.call_kappa_func_c_naire(c);
 
          //Div_kappa_grad
          div_kappa_grad(mutilde, kappaMatrix_variable, prov_elem);
 
        }
 
      // 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;
      //Cerr <<"c_demi"<<c_demi<<finl;
    }
  else
    {
      Cerr << "Type de schema errone !!" << finl;
    }
}
 
/*! @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_Naire::calculer_div_alpha_grad(const DoubleTab& c, DoubleTab& div_alpha_gradC) const
{
  // NOTE: Le tableau c donné en entrée est une concentration, ou doit être similaire à une concentration
  // Dans le cas n-aire, la concentration est un tableau à deux colonnes : une sur les éléments du domaine, et une autre sur le
  // nombre d'équations de Cahn-Hilliard.
 
  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();
 
  const int nb_comp = fermeture.get_nb_constituants();
  // Attention ici, nb_comp est le nombre de phases du système, pas la dimension du domaine ( /= nb_compo)
 
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
  prov_face = 0.;
  DoubleTab temp_prov_face2(prov_face.dimension(0), nb_comp);
  temp_prov_face2 = 0;
  // Grad(c)
  DoubleTab temp_resu(div_alpha_gradC.dimension_tot(0), 1);
  temp_resu = 0;
  DoubleTab temp_c(c.dimension_tot(0), 1);
  for (int j = 0; j < nb_comp; j++)
    {
      prov_face = 0.;
      for (int i = 0; i < temp_c.dimension(0); i++)
        {
          temp_c(i, 0) = c(i, j);
        }
      opgrad.calculer(temp_c, prov_face);
 
      // Application solveur masse
      eq_ns.solv_masse().appliquer(prov_face);
 
      // alpha*Grad(C)
      for (int i = 0; i < prov_face.dimension(0); i++)
        {
          for (int k = 0; k < nb_comp; k++)
            temp_prov_face2(i, k) += prov_face(i) * fermeture.get_alphaMatrix()(j + k * nb_comp);
          prov_face(i) = temp_prov_face2(i, j);
        }
 
      // Div(Grad(c))
      opdiv.calculer(prov_face, temp_resu);
      //opdiv.calculer(temp_prov_face2,div_alpha_gradC);
      eq_c.solv_masse().appliquer(temp_resu);
 
      for (int i = 0; i < temp_resu.dimension(0); i++)
        div_alpha_gradC(i, j) = temp_resu(i, 0);
    }
 
}
 
/*! @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_Naire::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 DoubleTab& c = eq_c.inconnue().valeurs();
//  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;
 
  // To do ?
 
}
 
/*! @brief Assemble la matrice pour le calcul du point fixe
 *
 */
void Source_Con_Phase_field_Naire::assembler_matrice_mobilite_explicite(Matrice_Morse& matrice_diffusion_CH)
{
  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;
  const double theta = 0.6;
  double kappa_naire_interpolee = 0.;
 
  if (mobilite_explicite_ == 1)
    {
      Cerr << "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_Matrix_var(nb_elem, nb_equation_CH_ * nb_equation_CH_);
      if (fermeture.get_kappa_ind() == 0)
        {
          for (int i = 0; i < nb_elem; i++)
            for (int j = 0; j < nb_equation_CH_ * nb_equation_CH_; j++)
              kappa_Matrix_var(i, j) = fermeture.get_kappaMatrix()(j);
        }
      else
        {
          kappa_Matrix_var = fermeture.call_kappa_func_c_naire(c);
        }
      //Cerr <<"kappaMatrix_Var = "<< kappa_Matrix_var<< finl;
 
      // Forme de la matrice pour un systeme -naire :
      // ici on montre uniquement pour un systeme quaternaire
      //                  ( 1 0 0 A B C )
      //                  ( 0 1 0 D E F )
      //                  ( 0 0 1 G H J )
      //                  ( K L M 1 0 0 )
      //                  ( N P Q 0 1 0 )
      //                  ( U V W 0 0 1 )
 
      // Assemblage de la moitie superieure de la matrice : 1 0 0 A B C
 
      // On dimensionne la matrice A et donc en meme temps toutes les matrices non unitaire B, C, ...(meme nombre de termes non nuls)
      // Pour cela : dimension du vecteur coeff pour chaque matrice A = (2 * dim + 1) * nb_elem - nombre_de_bords
 
      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--;
            }
        }
      // Sur la meme ligne, on a 1 matrice unitaire et nb_equation_CH* matrice non_unitaire et non nul. Pour cela,
      // On multiplie donc le nombre d'elements non nuls de la matrice A par le nb_equation_CH pour avoir le nombre d'element nnz sur la meme ligne
      // Puis on ajoute le nombre de 1 venant de la matrice unitaire sur la premiere ligne.
      // dimensionnement est le nombre d'element non nuls (nnz) dans la grosse matrice_diffusion_CH
      // Le dimensionnement total est donc le nnz de la premiere ligne multiplie par 2*nb_equation_CH
      dimensionnement *= nb_equation_CH_;
      dimensionnement += nb_elem;
      dimensionnement *= 2 * nb_equation_CH_;
 
      // Allocation des tableaux specifiques a la matrice morse
      matrice_diffusion_CH.dimensionner(nb_equation_CH_ * (2 * nb_elem), dimensionnement);
      auto& coeff = matrice_diffusion_CH.get_set_coeff();
      auto& tab2 = matrice_diffusion_CH.get_set_tab2();
      auto& tab1 = matrice_diffusion_CH.get_set_tab1();
 
      DoubleVect valeur_diag_naire(nb_equation_CH_);
      //double valeur_diag_naire;
 
      // Boucle sur le nombre de matrice dans la grosse matrice (matrice d'une matrice)
      for (int ligne = 0; ligne < nb_equation_CH_; ligne++)
        {
          // Boucle sur le nombre d'elements
          for (int elem = 0; elem < nb_elem; elem++)
            {
              valeur_diag_naire = 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) = (ligne * nb_elem) + 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++;
                }
 
              for (int ncomp = 0; ncomp < 2 * nb_compo_ - nb_faces_au_bord; ncomp++)
                {
                  min_tri = nb_elem + 1; //initialisation de min_tri pour chaque face (une sorte de reference car le numero des voisins ne doit pas depasser nb_elem+1)
                  dt = eq_c.schema_temps().pas_de_temps();    // on fait un calcul sur un demi pas de temps
                  for (int mat = 0; mat < nb_equation_CH_; mat++)
                    {
                      int j = (ligne * nb_equation_CH_) + mat;
                      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, harmonique ou arithmetique) des valeurs des mobilites
                          // des elements concernes par la face de calcul ("kappa_naire_interpolee")
 
                          //Cerr <<"f0 = "<<f0<<finl;
 
                          if (voisin != -1 && old_tri == -1)
                            {
                              dvar2 = pow((positions(elem, ori(f0)) - positions(voisin, ori(f0))), 2);
                              if (kappa_moy_ == 2)
                                kappa_naire_interpolee = pow(std::fabs(kappa_Matrix_var(elem, j) * kappa_Matrix_var(voisin, j)), 0.5);    // Moyenne geometrique
                              else if (kappa_moy_ == 1)
                                {
                                  if (kappa_Matrix_var(elem, j) == 0 || kappa_Matrix_var(voisin, j) == 0)
                                    kappa_naire_interpolee = 0;
                                  else
                                    kappa_naire_interpolee = 2. / (1. / kappa_Matrix_var(elem, j) + 1. / kappa_Matrix_var(voisin, j));  // Moyenne harmonique
                                }
                              else if (kappa_moy_ == 0)
                                kappa_naire_interpolee = (kappa_Matrix_var(elem, j) + kappa_Matrix_var(voisin, j)) / 2.;        // Moyenne arithmetique
                              valeur_diag_naire(mat) += theta * kappa_naire_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);
                            }
                        }
 
                      // 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_naire(mat);
                          tab2(compt2) = nb_elem * (mat + nb_equation_CH_) + elem; //la ou commence les matrices A, B, C, etc... donc (nb_elem*nb_equation_CH)+1 ??
                          compt++;
                          compt2++;
                        }
 
                      // On complete les tableaux avec les valeurs dans les voisins non aux bords du domaine
                      if (kappa_moy_ == 2)
                        kappa_naire_interpolee = pow(std::fabs(kappa_Matrix_var(elem, j) * kappa_Matrix_var(min_tri, j)), 0.5);    // Moyenne geometrique
                      else if (kappa_moy_ == 1)
                        {
                          if (kappa_Matrix_var(elem, j) == 0 || kappa_Matrix_var(min_tri, j) == 0)
                            kappa_naire_interpolee = 0;
                          else
                            kappa_naire_interpolee = 2. / (1. / kappa_Matrix_var(elem, j) + 1. / kappa_Matrix_var(min_tri, j));  // Moyenne harmonique
                        }
                      else if (kappa_moy_ == 0)
                        kappa_naire_interpolee = (kappa_Matrix_var(elem, j) + kappa_Matrix_var(min_tri, j)) / 2.;        // Moyenne arithmetique
 
                      coeff(compt) = -theta * kappa_naire_interpolee * (dt / dvarkeep);
                      tab2(compt2) = min_tri + nb_elem * (nb_equation_CH_ + mat);
                      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_naire(mat);
                          tab2(compt2) = elem + (nb_equation_CH_ + mat) * nb_elem;
                          compt++;
                          compt2++;
                        }
                    }
                  old_tri = min_tri;
                }
            }
 
        }
 
      // Assemblage de la moitie superieure de la matrice : B et I
      // Boucle sur le nombre de matrice dans la grosse matrice (matrice d'une matrice)
      for (int ligne = 0; ligne < nb_equation_CH_; ligne++)
        {
          // Boucle sur le nombre d'elements
          for (int elem = 0; elem < nb_elem; elem++)
            {
              valeur_diag_naire = 0.;
              old_tri = -1;
              nb_faces_au_bord = 0;
 
              // On saute de ligne en prenant la valeur de compte2 -1
              tab1(compt1) = compt2 + 1;
              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++;
                }
 
 
              for (int ncomp = 0; ncomp < 2 * nb_compo_ - nb_faces_au_bord; ncomp++)
                {
                  min_tri = nb_elem + 1; //initialisation de min_tri pour chaque face (une sorte de reference car le numero des voisins ne doit pas depasser nb_elem+1)
                  for (int mat = 0; mat < nb_equation_CH_; mat++)
                    {
                      int j = (ligne * nb_equation_CH_) + mat;
                      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, harmonique ou arithmetique) des valeurs des mobilites
                          // des elements concernes par la face de calcul ("kappa_naire_interpolee")
 
                          if (voisin != -1 && old_tri == -1)
                            {
                              dvar2 = pow((positions(elem, ori(f0)) - positions(voisin, ori(f0))), 2);
                              valeur_diag_naire(mat) += -fermeture.get_alphaMatrix()(j) / 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_naire(mat);
                          tab2(compt2) = nb_elem * (mat) + elem; //la ou commence les matrices A, B, C, etc... donc (nb_elem*nb_equation_CH)+1 ??
                          compt++;
                          compt2++;
                        }
 
                      // On complete les tableaux avec les valeurs dans les voisins non aux bords du domaine
                      coeff(compt) = fermeture.get_alphaMatrix()(j) / dvarkeep;
                      tab2(compt2) = min_tri + nb_elem * (mat);
                      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_naire(mat);
                          tab2(compt2) = elem + (mat) * nb_elem;
                          compt++;
                          compt2++;
                        }
                      // On sauve le numero d'element "minimum"
                    }
                  old_tri = min_tri;
                  //}
 
                }
              coeff(compt) = 1;
              tab2(compt2) = ((ligne + nb_equation_CH_) * nb_elem) + 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;
      matrice_diffusion_CH.clean();
    }
  else
    {
      Cerr << "STOP: argument mobilite_explicite cannot be different than 0 or 1." << finl;
      Cerr << "==> Must have an error in the code." << finl;
    }
 
}
 
/*! @brief Construire le residu du JFNK
 *
 */
void Source_Con_Phase_field_Naire::calculer_residu_jfnk(const DoubleTab& c, const Matrice_Morse& matrice_diffusion_CH, DoubleTab& v0, const DoubleTab& x1)
{
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const int nb_elem = domaine_VDF.nb_elem_tot();
  const int nb_elem_tot = c.size_totale(); // WARNING: dans le cas n-aire, size_totale comprend le nombre d'élément + le nombre de constituants
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
 
  double dt;
  const double theta = 0.6;
 
  DoubleVect second_membre(2 * nb_elem_tot);
  DoubleTab Ax1(2 * nb_elem_tot);
  DoubleTab terme_non_lin(c);
  DoubleTab x1_c(c);
 
  for (int j = 0; j < nb_equation_CH_; j++)
    for (int n_elem = 0; n_elem < nb_elem; n_elem++)
      x1_c(n_elem, j) = x1(n_elem + (j * nb_elem)); // Contient la concentration au temps intermédiaire
 
  x1_c.echange_espace_virtuel();
 
  dt = eq_c.schema_temps().pas_de_temps(); // Pas de temps
 
  // Assemblage du second membre
  terme_non_lin = fermeture.call_dWdc_naire(x1_c);
  for (int j = 0; j < nb_equation_CH_; j++)
    {
      for (int n_elem = 0; n_elem < nb_elem; n_elem++)
        {
          second_membre(n_elem + (j * nb_elem)) = c(n_elem, j);
          if (fermeture.get_type_potentiel_analytique() == 0)
            terme_non_lin(n_elem, j) *= fermeture.get_betaMatrix()(j);
          second_membre(nb_elem_tot + n_elem + (j * nb_elem)) = terme_non_lin(n_elem, j);
        }
    }
 
  if (mobilite_explicite_ == 1)
    {
      // Calcul du produit matrice / vecteur utilise
      matrice_diffusion_CH.multvect_(x1, Ax1);
      // Modifie par DJ
      //---------------
      {
        DoubleTab Ax1_c(nb_elem_tot);
        DoubleTab Ax1_mutilde(nb_elem_tot);
        for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
          {
            Ax1_c(n_elem) = Ax1(n_elem);
            Ax1_mutilde(n_elem) = Ax1(n_elem + nb_elem_tot);
          }
 
        Ax1_c.echange_espace_virtuel();
        Ax1_mutilde.echange_espace_virtuel();
 
        for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
          {
            Ax1(n_elem) = Ax1_c(n_elem);
            Ax1(n_elem + nb_elem_tot) = Ax1_mutilde(n_elem);
          }
      }
 
 
      // Calcul de v0 = [A(xn) xn - bn]
      for (int n_elem = 0; n_elem < 2 * nb_elem_tot; n_elem++)
        v0(n_elem) = (Ax1(n_elem) - second_membre(n_elem));
    }
  else if (mobilite_explicite_ == 0)
    {
 
      // Comme la matrice A n'est pas construite, on calcule directement sous forme de vecteur le produit AX
      // (X contenant la concentration et le potentiel chimique au temps intermédiaire theta)
      // Finalement le résidu évalué est H(X) = a(X) - b(X) [différence entre deux vecteurs dépendant de X]
      DoubleTab x1_mutilde(c);
      DoubleTab div_alpha_gradC(c);
      DoubleTab div_kappa_grad_mu(c);
      DoubleTab kappa_Matrix_var(nb_elem, nb_equation_CH_ * nb_equation_CH_);
      int indice_tot=0;
 
      for (int j = 0; j < nb_equation_CH_; j++)
        for (int n_elem = 0; n_elem < nb_elem; n_elem++)
          x1_mutilde(n_elem, j) = x1(nb_elem_tot + n_elem + (j * nb_elem)); // Contient le potentiel chimique au temps intermédiaire
 
      x1_mutilde.echange_espace_virtuel();
 
      if (fermeture.get_kappa_ind() == 0)
        {
          for (int i = 0; i < nb_elem; i++)
            for (int j = 0; j < nb_equation_CH_ * nb_equation_CH_; j++)
              kappa_Matrix_var(i, j) = fermeture.get_kappaMatrix()(j);
        }
      else
        {
          kappa_Matrix_var = fermeture.call_kappa_func_c_naire(x1_c);
        }
 
      calculer_div_alpha_grad(x1_c, div_alpha_gradC);
      div_kappa_grad(x1_mutilde, kappa_Matrix_var, div_kappa_grad_mu);
 
      Cerr << "La taille de mu = " << x1_mutilde.dimension(0) << ", " << x1_mutilde.dimension(1) << finl;
      Cerr << "La taille de div_alpha = " << div_alpha_gradC.dimension(0) << ", " << div_alpha_gradC.dimension(1) << finl;
      Cerr << "La taille de div_kappa = " << div_kappa_grad_mu.dimension(0) << ", " << div_kappa_grad_mu.dimension(1) << finl;
 
      for (int j = 0; j < nb_equation_CH_; j++)
        {
          for (int n_elem = 0; n_elem < nb_elem; n_elem++)
            {
              indice_tot = n_elem + (j * nb_elem);
//            Cerr << "indice n_elem = " << n_elem << finl;
//            Cerr << "indice j = " << j << finl;
//            Cerr << "indice indice_tot = " << indice_tot << finl;
              Ax1(indice_tot) = x1_c(n_elem, j) - theta * dt * div_kappa_grad_mu(n_elem, j);
              Ax1(indice_tot + nb_elem_tot) = x1_mutilde(n_elem, j) + div_alpha_gradC(n_elem, j);
            }
        }
 
      // Calcul de v0 = [a(xn) - bn]
      for (int n_elem = 0; n_elem < 2 * nb_elem_tot; n_elem++)
        v0(n_elem) = (Ax1(n_elem) - second_membre(n_elem));
    }
 
}
 
/*! @brief Construire le residu du JFNK
 *
 */
void Source_Con_Phase_field_Naire::jacobian_vect(const DoubleTab& c, const Matrice_Morse& matrice_diffusion_CH, 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, matrice_diffusion_CH, v1, x1_);
  calculer_residu_jfnk(c, matrice_diffusion_CH, v2, x1);
 
  v1 -= v2;
  v1 /= delta;
}
 
/*! @brief Algorithme résolution non linéaire par JFNK
 *
 * @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_Naire::resolution_jfnk(const DoubleTab& c, const DoubleTab& mutilde, const Matrice_Morse& matrice_diffusion_CH, DoubleTab& c_demi, DoubleTab& mutilde_demi)
 
{
  int i, j, nk, i0, im, it, ii;
  double tem = 1., res, ccos, ssin;
  const int ns = 2 * c.size_totale();
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const int nb_elem = domaine_VDF.nb_elem_tot();
 
  // Ajoute par DJ
  //--------------
  int nb_elem_tot = c.size_totale();
  DoubleTab x1(ns);
  x1 = 0.;
  for (int ncomponent = 0; ncomponent < nb_equation_CH_; ncomponent++)
    {
      for (int nelem = 0; nelem < nb_elem; nelem++)
        {
          x1(nelem + (ncomponent * nb_elem)) = c(nelem, ncomponent);
          x1(nelem + (ncomponent * nb_elem) + nb_elem_tot) = mutilde(nelem, ncomponent);
        }
    }
 
  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.;
 
  calculer_residu_jfnk(c, matrice_diffusion_CH, v0, x1);
  v0 *= -1.;
 
  res = 0.;
 
  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);
 
  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_gmres_; 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, matrice_diffusion_CH, v0, x1, v1);
          v0 = v1;
 
          // Modifie par DJ
          //---------------
          for (i = 0; i <= j; i++)
            {
              //               for (ii=0; ii<ns;ii++) h(i,j) += v0(ii) * v(ii,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);
 
              //for (ii=0; ii<ns; ii++) h(i,j) += v0(ii) * v(ii,i);
 
              h(i, j) = mp_sum(h(i, j));
              //               Debog::verifier("JFNK h(i,j) : ",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<ns;ii++) tem += v0(ii) * v0(ii) ;
 
          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);
          //           Debog::verifier("JFNK tem apres mp_sum : ",tem);
          tem = sqrt(tem);
          h(j + 1, j) = tem;
          if (tem < rtol_)
            {
              nk = j + 1;
              goto l5naire;
            }
          v0 /= tem;
        }
 
      //...Triangularisation
l5naire:
      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);
 
      calculer_residu_jfnk(c, matrice_diffusion_CH, v0, x1);
      v0 *= -1.;
 
 
      res = 0.;
 
      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);
 
      Cerr << " - At it = " << it << ", residu scalar = " << res << finl;
 
      if (res < rtol_)
        {
          // Ajoute par DJ
          //--------------
          for (int ncomponent = 0; ncomponent < nb_equation_CH_; ncomponent++)
            {
              for (int n_elem = 0; n_elem < nb_elem; n_elem++)
                {
                  c_demi(n_elem, ncomponent) = x1(n_elem + (ncomponent * nb_elem));
                  mutilde_demi(n_elem, ncomponent) = x1(n_elem + nb_elem_tot + (ncomponent * nb_elem));
                }
            }
          // L'echange espace virtuel est fait dans premier_demi_dt()
          Cerr << "Number of iterations to reach convergence : " << it + 1 << finl;
          Cerr << "" << finl;
 
          return it;
        }
      else if (it == maxit_gmres_ - 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 ncomponent = 0; ncomponent < nb_equation_CH_; ncomponent++)
            {
              for (int n_elem = 0; n_elem < nb_elem; n_elem++)
                {
                  c_demi(n_elem, ncomponent) = x1(n_elem + (ncomponent * nb_elem));
                  mutilde_demi(n_elem, ncomponent) = x1(n_elem + nb_elem_tot + (ncomponent * nb_elem));
                }
            }
 
          Cerr << "Stopped before convergence" << finl;
          Cerr << "" << finl;
        }
    }
 
  return -1;
}
 
/*! @brief Calcul de mutilde au centre des elements
 *
 * @param (mutilde) mutilde au centre des elements
 */
void Source_Con_Phase_field_Naire::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;
 
  DoubleTab potent_chimique(mutilde);
 
  potent_chimique = fermeture.call_dWdc_naire(c);
  if (fermeture.get_type_potentiel_analytique() == 0)
    {
      for (int j = 0; j < c.line_size(); j++)
        for (int i = 0; i < c.dimension(0); i++)
          potent_chimique(i, j) *= fermeture.get_betaMatrix()(j);
    }
 
  for (int j = 0; j < mutilde.line_size(); j++)
    for (int i = 0; i < mutilde.dimension(0); i++)
      {
        mutilde(i, j) += potent_chimique(i, j);
        if (mutype_ == 1) //avec Ec, n'est pas implemente pour le cas multicomposant
          Cerr << "potentiel_chimique_generalisee avec_energie_cinetique not implemented for multicomponent system" << finl;
      }
 
 
// L'espace virtuel n'est pas a jour
  assert_invalide_items_non_calcules(mutilde);
}
 
void Source_Con_Phase_field_Naire::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.;
 
  const int nb_comp = fermeture.get_nb_constituants();
 
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
  prov_face = 0.;
  DoubleTab temp_prov_face2(prov_face.dimension(0), nb_comp);
  temp_prov_face2 = 0;
  // Grad(c)
  DoubleTab temp_resu(mutilde_demi.dimension_tot(0), 1);
  temp_resu = 0;
  DoubleTab temp_c(c_demi.dimension_tot(0), 1);
  for (int j = 0; j < nb_comp; j++)
    {
      prov_face = 0.;
      for (int i = 0; i < temp_c.dimension(0); i++)
        {
          temp_c(i, 0) = c_demi(i, j);
        }
      opgrad.calculer(temp_c, prov_face);
 
      // Application solveur masse
      eq_ns.solv_masse().appliquer(prov_face);
 
      // alpha*Grad(C)
 
      for (int i = 0; i < prov_face.dimension(0); i++)
        {
          for (int k = 0; k < nb_comp; k++)
            temp_prov_face2(i, k) += prov_face(i) * fermeture.get_alphaMatrix()(j + k * nb_comp);
 
          prov_face(i) = temp_prov_face2(i, j);
        }
 
      // Div(Grad(c))
      opdiv.calculer(prov_face, temp_resu);
      //opdiv.calculer(temp_prov_face2,div_alpha_gradC);
      eq_c.solv_masse().appliquer(temp_resu);
 
      for (int i = 0; i < temp_resu.dimension(0); i++)
        {
          mutilde_demi(i, j) = temp_resu(i, 0);
        }
    }
 
  mutilde_demi *= -1.;
 
  DoubleTab potent_chimique(mutilde_demi);
 
  potent_chimique = fermeture.call_dWdc_naire(c_demi);
  if (fermeture.get_type_potentiel_analytique() == 0)
    {
      for (int j = 0; j < c_demi.line_size(); j++)
        for (int i = 0; i < c_demi.dimension(0); i++)
          potent_chimique(i, j) *= fermeture.get_betaMatrix()(j);
    }
  for (int j = 0; j < mutilde_demi.line_size(); j++)
    for (int i = 0; i < mutilde_demi.dimension(0); i++)
      {
        mutilde_demi(i, j) += potent_chimique(i, j);
        if (mutype_ == 1) //avec Ec, n'est pas implemente pour le cas multicomposant
          Cerr << "potentiel_chimique_generalisee avec_energie_cinetique not implemented for multicomponent system" << finl;
      }
 
// L'espace virtuel n'est pas a jour
  assert_invalide_items_non_calcules(mutilde_demi);
}