Source_Con_Phase_field_Binaire_Compact.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_Compact.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_Compact,"Source_Con_Phase_field_Binaire_Compact_VDF_P0_VDF",Source_Con_Phase_field);
 
 
// XD Source_Con_Phase_field_Binaire_Compact_VDF_P0_VDF Source_Con_Phase_field Source_Con_Phase_field_Binaire_Compact_VDF_P0_VDF INHERITS_BRACE Source for phase field
 
 
Sortie& Source_Con_Phase_field_Binaire_Compact::printOn(Sortie& s) const
{
  return s << que_suis_je();
}
 
Entree& Source_Con_Phase_field_Binaire_Compact::readOn(Entree& is)
{
  Cerr << "Source_Con_Phase_field_Binaire_Compact::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 FormFunction" << finl;
    }
  return is;
}
 
DoubleTab& Source_Con_Phase_field_Binaire_Compact::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_Compact::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_Compact::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_Compact::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_Compact::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_alpha_initialised_)
        {
          initialiser_matrice_alpha();
//          Cout << "Matrice avec alpha = " << finl;
//          matrice_alpha_.imprimer_formatte(Cout);
        }
 
      if (mobilite_explicite_ == 1 && !is_matrix_kappa_initialised_)
        {
          initialiser_matrice_kappa();
//          Cout << "Matrice avec kappa = " << finl;
//          matrice_kappa_.imprimer_formatte(Cout);
        }
 
      // Modifie par DJ
      //---------------
 
      // RESOLUTION JFNK
      // =========================================================
      resolution_jfnk(c, mutilde, c_demi, mutilde_demi);
      // =========================================================
 
      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;
        }
 
      // 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_Compact::initialiser_matrice_kappa()
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const DoubleTab& c = eq_c.inconnue().valeurs();
 
  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();
 
  terme_diffusif_vdf_.associer_diffusivite(ch_kappa_ou_alpha_.valeur());
  terme_diffusif_vdf_->associer_diffusivite(ch_kappa_ou_alpha_.valeur());
  terme_diffusif_vdf_->dimensionner(matrice_kappa_);
  matrice_kappa_.sort_stencil(); // XXX sort tab2 ... au cas ou ...
  update_matrice_kappa(c);
 
  is_matrix_kappa_initialised_ = true;
 
}
 
void Source_Con_Phase_field_Binaire_Compact::initialiser_matrice_alpha()
{
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
  const auto& fermeture = mil.get_fermeture();
  const DoubleTab& c = eq_c.inconnue().valeurs();
 
  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();
 
  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(matrice_alpha_);
  matrice_alpha_.sort_stencil(); // XXX sort tab2 ... au cas ou ...
  terme_diffusif_vdf_->contribuer_a_avec(c, matrice_alpha_);
 
  // Vecteur contenant la contribution des conditions limites
  vec_cl_ = c; // on initialise la structure
  vec_cl_ = 0.;
 
  terme_diffusif_vdf_->ajouter(c,vec_cl_);
  matrice_alpha_.ajouter_multvect(c, vec_cl_); // Ajout de A*Inco(n)
 
  Cout << "Vecteur de conditions limites = " << eq_c.domaine_Cl_dis();
  Cout << "Vecteur de conditions limites vec_cl = " << vec_cl_;
 
  is_matrix_alpha_initialised_ = true;
 
}
 
void Source_Con_Phase_field_Binaire_Compact::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_Compact::update_matrice_kappa(const DoubleTab& donne)
{
  const auto& fermeture = ref_cast(Milieu_Phase_field, pb_->milieu()).get_fermeture();
 
  if (is_matrix_kappa_initialised_ && fermeture.get_type_kappa() ==0)
    return;
 
  update_tab_kappa(donne);
 
  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(), matrice_kappa_);
 
}
 
/*! @brief Construire le residu du JFNK
 *
 */
 
void Source_Con_Phase_field_Binaire_Compact::calculer_residu_jfnk(const DoubleTab& c, DoubleTab& v0, const DoubleTab& c_demi)
{
  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(nb_elem);
  DoubleTab B_c_demi(nb_elem);
  DoubleTab A_B_c_demi_sec_membre(nb_elem);
 
  for (int n_elem = 0; n_elem < nb_elem; 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) = term_cin(n_elem) + fermeture.get_beta() * fermeture.call_dWdc(c_demi(n_elem));
    }
 
  // On calcule H(c_demi) = c_demi + A ( B . c_demi - beta*dW/dc(c_demi) ) - c
  // Où A = matrice_kappa et B = matrice_alpha
 
  // Calcul du premier produit matrice / vecteur B . c_demi
  matrice_alpha_.multvect_(c_demi, B_c_demi);
  // Calcul du second produit matrice / vecteur A . (B . c_demi - beta*dW/dc)
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      second_membre(n_elem) = second_membre(n_elem) - B_c_demi(n_elem);
    }
  matrice_kappa_.multvect_(second_membre, A_B_c_demi_sec_membre);
 
 
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      v0(n_elem) = c_demi(n_elem) - c(n_elem) + A_B_c_demi_sec_membre(n_elem);
    }
 
}
 
/*! @brief Construire la jacobienne multipliée par un vecteur Y
 *
 */
void Source_Con_Phase_field_Binaire_Compact::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 tabToVecBis(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 vecToTabBis(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 FormFunction de PETSc : calcule le résidu du JFNK
 *
 */
PetscErrorCode FormFunctionCompactOld(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_Compact& source = ref_cast(Source_Con_Phase_field_Binaire_Compact, eqn->sources().dernier().valeur());
  const int nb_elem = C.dimension_tot(0);
  DoubleTab c_demi(nb_elem);
 
  vecToTabBis(x,c_demi);
 
 
  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(c_demi);
 
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      second_membre(n_elem) = fermeture.get_beta() * fermeture.call_dWdc(c_demi(n_elem));
    }
 
  DoubleTrav B_c_demi(nb_elem);
  DoubleTrav A_B_c_demi_sec_membre(nb_elem);
 
  // On calcule H(c_demi) = c_demi + A ( B . c_demi - beta*dW/dc(c_demi) ) - c
  // Où A = matrice_kappa et B = matrice_alpha
 
  // Calcul du premier produit matrice / vecteur B . c_demi
  source.get_matrice_alpha().multvect_(c_demi, B_c_demi);
  // Calcul du second produit matrice / vecteur A . (B . c_demi - beta*dW/dc)
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      second_membre(n_elem) = second_membre(n_elem) - B_c_demi(n_elem);
    }
 
  source.get_matrice_kappa().multvect_(second_membre, A_B_c_demi_sec_membre);
 
  for (int n_elem = 0; n_elem < nb_elem; n_elem++)
    {
      second_membre(n_elem) = c_demi(n_elem) - C(n_elem) + A_B_c_demi_sec_membre(n_elem);
    }
 
  tabToVecBis(second_membre, f);
 
  return 0;
}
 
PetscErrorCode MySNESMonitorBis(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 LuisKSPMonitorBis(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_Compact::petsc_snes(const DoubleTab& c, const DoubleTab& mutilde, DoubleTab& c_demi, DoubleTab& mutilde_demi)
{
 
  const int nb_elem_tot = c.size_totale();
 
  Cout << "Size total = " << c.size_totale() << finl;
 
  // 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, nb_elem_tot);
  VecSetFromOptions(x);
 
  // INITIAL GUESS of X
  tabToVecBis(c, x);
 
  VecDuplicate(x, &r);
 
  // Vérification tailles
  PetscInt x_size, r_size;
  VecGetSize(x, &x_size);
  VecGetSize(r, &r_size);
  assert(x_size == nb_elem_tot);
  assert(r_size == nb_elem_tot);
 
  // ---------------------- Compute real values of tolerances ---------------------------
  // Appel explicite de ta fonction de résidu
  FormFunctionCompactOld(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 FormFunction to the SNES context ============
  SNESSetFunction(snes, r, FormFunctionCompactOld, &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, LuisKSPMonitorBis, nullptr, nullptr);
  SNESMonitorSet(snes, MySNESMonitorBis, &equation(), nullptr);
  //SNESMonitorSet(snes, MySNESMonitorBis, nullptr, nullptr);
 
  // =============================== Setting from options ===================================
  // KSPSetFromOptions(ksp);                              // Lire les options
 
  // Injecte les options comme si elles venaient de la ligne de commande
  //  PetscOptionsSetValue(nullptr, "-pc_type", "none");              // désactive le préconditionneur
  //  PetscOptionsSetValue(nullptr, "-snes_mf", nullptr);             // active matrix-free
  // PetscOptionsSetValue(nullptr, "-snes_mf_relstep", "1e-5");  // contrôle l’epsilon de l’approximation par différences finies
 
  PetscOptionsView(nullptr, PETSC_VIEWER_STDOUT_WORLD);  // Visu de la liste des options
  PetscOptionsLeft(nullptr);                             // Liste celles qui n'ont PAS été utilisées
 
  // Lecture des options ligne de commande de SNES
  SNESSetOptionsPrefix(snes, "mysolver_");
 
  SNESSetFromOptions(snes);
 
  SNESView(snes, PETSC_VIEWER_STDOUT_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(nb_elem_tot); // inconnue (vector de C et mutilde)
 
  vecToTabBis(x,solution);
 
  for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
    {
      c_demi(n_elem) = solution(n_elem);
    }
 
  DoubleTrav mu_temp(nb_elem_tot);
  DoubleTrav B_c_demi(nb_elem_tot);
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, equation().milieu());
  const auto& fermeture = mil.get_fermeture();
 
  for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
    {
      mu_temp(n_elem) = fermeture.get_beta() * fermeture.call_dWdc(c_demi(n_elem));
    }
 
  // On calcule H(c_demi) = c_demi + A ( B . c_demi - beta*dW/dc(c_demi) ) - c
  // Où A = matrice_kappa et B = matrice_alpha
 
  // Calcul du produit matrice / vecteur B . c_demi, puis soustraction du vecteur beta*dW/dc
  matrice_alpha_.multvect_(c_demi, B_c_demi);
  for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
    {
      mutilde_demi(n_elem) = mu_temp(n_elem) - B_c_demi(n_elem);
    }
 
  // 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_Compact::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 nb_elem_tot = c.size_totale();
      double temps_exec = 0.;
      int nfunc=0;
      // statistiques().set_debug_level(9);
 
      // Ajoute par DJ
      //--------------
      DoubleTrav x1(nb_elem_tot);
 
      for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
        {
          x1(n_elem) = c(n_elem);
        }
 
      DoubleTab v(nb_elem_tot, 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
 
      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 < nb_elem_tot; 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 < nb_elem_tot; 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 < nb_elem_tot; ii++)
              x1(ii) += r[i] * v(ii, i);
 
          if (mobilite_explicite_ == 0) /* implicite */
            {
              update_matrice_kappa(x1);
            }
 
          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;
 
          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_ || it == maxit_newton_ - 1)
            {
              // Ajoute par DJ
              //--------------
              for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
                {
                  c_demi(n_elem) = x1(n_elem);
                }
 
              DoubleTrav mu_temp(nb_elem_tot);
              DoubleTrav B_c_demi(nb_elem_tot);
              const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, equation().milieu());
              const auto& fermeture = mil.get_fermeture();
 
              for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
                {
                  mu_temp(n_elem) = fermeture.get_beta() * fermeture.call_dWdc(c_demi(n_elem));
                }
 
              // On calcule H(c_demi) = c_demi + A ( B . c_demi - beta*dW/dc(c_demi) ) - c
              // Où A = matrice_kappa et B = matrice_alpha
 
              // Calcul du produit matrice / vecteur B . c_demi, puis soustraction du vecteur beta*dW/dc
              matrice_alpha_.multvect_(c_demi, B_c_demi);
              for (int n_elem = 0; n_elem < nb_elem_tot; n_elem++)
                {
                  mutilde_demi(n_elem) = mu_temp(n_elem) - B_c_demi(n_elem);
                }
            }
 
          if (res < rtol_)
            {
              // 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)
            {
              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_Compact::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_Compact::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);
}