Source_Con_Phase_field.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 <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_base(Source_Con_Phase_field,"Source_Con_Phase_field_VDF_P0_VDF",Source_Con_Phase_field_base);
 
// XD source_con_phase_field source_base source_con_phase_field BRACE Keyword to define the source term of the
// XD_CONT Cahn-Hilliard equation.
// XD attr systeme_naire systeme_naire_deriv systeme_naire OPT not_set
// XD attr temps_d_affichage entier temps_d_affichage REQ Time during the caracteristics of the problem are shown before
// XD_CONT calculation.
// XD attr moyenne_de_kappa chaine moyenne_de_kappa REQ To define how mobility kappa is calculated on faces of the mesh
// XD_CONT according to cell-centered values (chaine is arithmetique/harmonique/geometrique).
// XD attr multiplicateur_de_kappa floattant multiplicateur_de_kappa REQ To define the parameter of the mobility
// XD_CONT expression when mobility depends on C.
// XD attr couplage_NS_CH chaine couplage_NS_CH REQ Evaluating time choosen for the term source calculation into the
// XD_CONT Navier Stokes equation (chaine is mutilde(n+1/2)/mutilde(n), in order to be conservative, the first choice
// XD_CONT seems better).
// XD attr implicitation_CH chaine(into=["oui","non"]) implicitation_CH REQ To define if the Cahn-Hilliard will be
// XD_CONT solved using a implicit algorithm or not.
// XD attr seuil_residu_jfnk floattant seuil_residu_jfnk REQ Convergence threshold (an option of the Newton-Krylov
// XD_CONT method).
// XD attr dimension_espace_de_krylov entier dimension_espace_de_krylov REQ Vector numbers used in the method (an option
// XD_CONT of the Newton-Krylov method).
// XD attr nb_iterations_jfnk entier nb_iterations_jfnk REQ Maximal iteration (an option of the Newton-Krylov method).
// XD attr residu_min_jfnk floattant residu_min_jfnk REQ Minimal convergence threshold (an option of the Newton-Krylov
// XD_CONT method).
// XD attr residu_max_jfnk floattant residu_max_jfnk REQ Maximal convergence threshold (an option of the Newton-Krylov
// XD_CONT method).
 
// XD bloc_kappa_variable objet_lecture nul NO_BRACE if the parameter of the mobility, kappa, depends on C
// XD attr expr bloc_lecture expr REQ choice for kappa_variable
 
// XD bloc_potentiel_chim objet_lecture nul NO_BRACE if the chemical potential function is an univariate function
// XD attr expr bloc_lecture expr REQ choice for potentiel_chimique
 
// XD systeme_naire_deriv objet_lecture systeme_naire_deriv INHERITS_BRACE not_set
// XD systeme_naire_non systeme_naire_deriv non BRACE not_set
// XD attr alpha floattant alpha REQ Internal capillary coefficient alfa.
// XD attr beta floattant beta REQ Parameter beta of the model.
// XD attr kappa floattant kappa REQ Mobility coefficient kappa0.
// XD attr kappa_variable bloc_kappa_variable kappa_variable REQ To define a mobility which depends on concentration C.
// XD attr potentiel_chimique bloc_potentiel_chim potentiel_chimique OPT chemical potential function
 
Sortie& Source_Con_Phase_field::printOn(Sortie& s) const
{
  return s << que_suis_je();
}
 
Entree& Source_Con_Phase_field::readOn(Entree& is)
{
  Param param(que_suis_je());
  param.ajouter("Temps_d_affichage", &tpsaff_, Param::REQUIRED);
  param.ajouter("absolute_tol", &atol_, Param::REQUIRED);
  param.ajouter("relative_tol", &rtol_, Param::REQUIRED);
  param.ajouter("stagn_tol", &stol_, Param::REQUIRED);
  param.ajouter("max_it_newton", &maxit_newton_, Param::REQUIRED);
  param.ajouter("max_it_GMRES", &maxit_gmres_, Param::REQUIRED);
  param.ajouter("dimension_espace_de_krylov", &nkr_, Param::REQUIRED);
  param.ajouter("matrix_free_erel", &erel_, Param::REQUIRED);
  param.ajouter("matrix_free_umin", &umin_, Param::REQUIRED);
  param.ajouter_non_std("moyenne_de_kappa", (this), Param::REQUIRED);
  param.ajouter_non_std("couplage_NS_CH", (this), Param::REQUIRED);
  param.ajouter_non_std("implicitation_CH", (this), Param::REQUIRED);
  param.ajouter_non_std("mobilite_explicite", (this), Param::REQUIRED);
  param.ajouter_non_std("resolution_petsc", (this), Param::REQUIRED);
 
  param.lire_avec_accolades_depuis(is);
 
  if (tpsaff_ < 0. || tpsaff_ > 100.)
    {
      Cerr << "Source_Con_Phase_field::readOn: Temps_d_affichage should be in the range 0 - 100 seconds." << finl;
      Process::exit();
    }
  return is;
}
 
int Source_Con_Phase_field::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  Motcle temp;
  if (mot == "implicitation_CH")
    {
      is >> temp;
      if (temp == "oui")
        implicitation_ = 1;
      else
        implicitation_ = 0;
      return 1;
    }
  else if (mot == "mobilite_explicite")
    {
      is >> temp;
      if (temp == "oui")
        mobilite_explicite_ = 1;
      else
        mobilite_explicite_ = 0;
      return 1;
    }
  else if (mot == "resolution_petsc")
    {
      is >> temp;
      if (temp == "oui")
        with_petsc_ = 1;
      else
        with_petsc_ = 0;
      return 1;
    }
  else if (mot == "couplage_NS_CH")
    {
      is >> temp;
      if (temp == "mutilde(n)")
        couplage_ = 0;
      else
        couplage_ = 1;
      return 1;
    }
  else if (mot == "moyenne_de_kappa")
    {
      is >> temp;
      if (temp == "arithmetique")
        kappa_moy_ = 0;
      else if (temp == "harmonique")
        kappa_moy_ = 1;
      else
        kappa_moy_ = 2;
      return 1;
    }
  return -1;
}
 
void Source_Con_Phase_field::associer_pb(const Probleme_base& pb)
{
  pb_ = pb;
  Navier_Stokes_phase_field& eq_ns = ref_cast(Navier_Stokes_phase_field, pb_->equation(0));
  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());
  rho0_ = mil.rho0();
  drhodc_ = mil.drhodc();
  // Dans le modele binaire cette derivee est constante - On la calcule au debut selon que l'approx boussi soit utilisee ou non.
 
  // UPDATE Luis: C'est ici que le cas du système binaire ou du système n-aire est testé par rapport aux fermetures --------------------------
  if (mil.get_fermeture().get_type_systeme_naire() == 0)
    {
      if ((this->que_suis_je() != "Source_Con_Phase_field_Binaire_VDF_P0_VDF") && (this->que_suis_je() != "Source_Con_Phase_field_Binaire_Compact_VDF_P0_VDF"))
        {
          Cerr << "Source term is of type " << this->que_suis_je() << finl;
          Cerr << "** Error: Closure type does not fit with phase field source term... **" << finl;
          Cerr << "** Check if they are both for binary systems. **" << finl;
          Process::exit();
        }
    }
  else if ((mil.get_fermeture().get_type_systeme_naire() == 1) && (this->que_suis_je() != "Source_Con_Phase_field_Naire_VDF_P0_VDF"))
    {
      Cerr << "Source term is of type " << this->que_suis_je() << finl;
      Cerr << "** Error: Closure type does not fit with phase field source term... **" << finl;
      Cerr << "** Check if they are both for n-ary systems. **" << finl;
      Process::exit();
    }
  // -----------------------------------------------------------------------------------------------------------------------------------------
 
  boussi_ = mil.get_boussi();
  if (boussi_ != 1 && boussi_ != 0)
    {
      Cerr << "Erreur dans le choix du parametre boussi_" << finl;
      Process::exit();
    }
 
  // diff_boussi==1 si la viscosite dynamique n'est pas constante (variation de la viscosite en fonction de la concentration c)
  diff_boussi_ = mil.get_diff_boussi();
  if (diff_boussi_ != 1 && diff_boussi_ != 0)
    {
      Cerr << "Erreur dans le choix du parametre diff_boussi_" << finl;
      Process::exit();
    }
 
  g_ = mil.gravite().valeurs();
 
  // mutype==1 si c'est avec energie cinetique (dans l'expression du potentiel chimique generalise)
  mutype_ = eq_c.get_mutype_();
  if (mutype_ != 0 && mutype_ != 1)
    {
      Cerr << "Erreur dans le choix du parametre mutype_" << finl;
      Process::exit();
    }
 
  if (implicitation_ != 0 && implicitation_ != 1)
    {
      Cerr << "=================================================" << finl;
      Cerr << "Choix de l'implicitation incorrect ! (ni 0, ni 1)" << finl;
      Cerr << "=================================================" << finl;
      Process::exit();
    }
 
  //choix de la methode de calcul de la moyenne des mobilites vis a vis du gradient de mutilde dans le cas implicite
  if (implicitation_ == 1)
    {
      if (kappa_moy_ != 0 && kappa_moy_ != 1 && kappa_moy_ != 2)
        {
          Cerr << "Erreur dans le choix de la moyenne de kappa !" << finl;
          Process::exit();
        }
    }
 
  // Recapitulatif des parametres
 
  Nom choix_boussi;
  if (boussi_ == 1)
    choix_boussi = "oui";
  else
    choix_boussi = "non";
 
  Nom choix_diff_boussi;
  if (diff_boussi_ == 1)
    choix_diff_boussi = "oui";
  else
    choix_diff_boussi = "non";
  if (boussi_ == 1 && diff_boussi_ == 1)
    Cerr << "ATTENTION : les options 'approximation de Boussinesq' et" << "'approximation de Boussinesq dans le terme de diffusion' ne devrait pas etre utilisees simultanement'."
         << "(voir Source_Con_Phase_field.cpp)" << finl;
 
  const int terme_source = eq_ns.getset_terme_source();
 
  Nom choix_source;
  if (terme_source == 1)
    choix_source = "c grad(mutilde)";
  else if (terme_source == 2)
    choix_source = "c grad(laplacien(c))";
  else if (terme_source == 3)
    choix_source = "c grad(laplacien(c))-div((grad(c))^2)/2";
  else if (terme_source == 4)
    choix_source = "-laplacien(c) grad(c)";
 
  Nom choix_implicite;
  if (implicitation_ == 1)
    choix_implicite = "oui";
  else
    choix_implicite = "non";
 
  Nom type_implicite;
  type_implicite = "JFNK";
 
  Nom mutilde_d;
  if (mutype_ == 1)
    mutilde_d = "avec le terme d'Ec";
  else
    mutilde_d = "sans le terme d'Ec";
 
  Nom type_couplage;
  if (couplage_ == 1)
    type_couplage = "potentiel au temps n+1/2";
  else
    type_couplage = "potentiel au temps n";
  Nom mobilite_variable;
  Nom moyenne_kappa;
 
  if (mil.get_fermeture().get_type_kappa() == 1 || mil.get_fermeture().get_type_kappa_auto_diffusion() == 1)
    mobilite_variable = "oui";
  else
    mobilite_variable = "non";
 
  if (kappa_moy_ == 0)
    moyenne_kappa = "arithmetique (+)";
  else if (kappa_moy_ == 1)
    moyenne_kappa = "harmonique (/)";
  else if (kappa_moy_ == 2)
    moyenne_kappa = "geometrique (*)";
 
  nb_equation_CH_ = mil.get_fermeture().get_nb_constituants();
 
 
  Cerr << "" << finl;
  Cerr << "" << finl;
  Cerr << "******************************************************************************" << finl;
  Cerr << "*********************** RECAPITULATIF DU PARAMETRAGE *************************" << finl;
  Cerr << "" << finl;
  Cerr << "1) Equation de Navier-Stokes" << finl;
  Cerr << "  - approximation de Boussinesq                     : " << choix_boussi << finl;
  if (boussi_ == 0)
    {
      Cerr << "  - approximation de Boussinesq dans la diffusion   : " << choix_diff_boussi << finl;
    }
  else
    {
      Cerr << "  - masse volumique de reference                    : " << rho0_ << " kg/m3" << finl;
    }
  Cerr << "  Dans le cas ou la viscosite dynamique est variable : " << finl;
  Cerr << "  - terme source                                    : " << choix_source << finl;
  Cerr << "  - intensite du champ de gravite                   : " << norme_array(g_) << " m/s^2" << finl;
  Cerr << "" << finl;
  Cerr << "2) Equation de Cahn-Hilliard" << finl;
  Cerr << "  - discretisation implicite                        : " << choix_implicite << finl;
  if (implicitation_ == 1)
    {
      Cerr << "  - algorithme implicite                            : " << type_implicite << finl;
    }
  Cerr << "  - potentiel chimique generalise                   : " << mutilde_d << finl;
  Cerr << "  - couplage NS / CH via le potentiel chimique      : " << type_couplage << finl;
 
  Cerr << "  - mobilite variable                               : " << mobilite_variable << finl;
  if (mil.get_fermeture().get_type_kappa() == 1 || mil.get_fermeture().get_type_kappa_auto_diffusion() == 1)
    {
      Cerr << "  - type de moyenne pour la mobilite                : " << moyenne_kappa << finl;
      Cerr << "  - multiplicateur de la mobilite                   : " << mil.get_fermeture().get_mult_kappa() << finl;
    }
 
  if (implicitation_ == 1)
    {
      Cerr << "  - dans le cas de l'algorithme de JFNK ... " << finl;
      if (with_petsc_ == 1)
        {
          Cerr << "       === USING PETSc - SNES LIBRARY ===      "   << finl;
          Cerr << "        Absolute tolerance on the residual          : " << atol_ << finl;
          Cerr << "        Relative tolerance on the residual          : " << rtol_ << finl;
          Cerr << "        Stagnation tolerance on the residual        : " << rtol_ << finl;
          Cerr << "        Krylov space dimension                      : " << nkr_ << finl;
          Cerr << "        Maximum number of Newton iterations         : " << maxit_newton_ << finl;
          Cerr << "        Maximum number of GMRES iterations          : " << maxit_gmres_ << finl;
          Cerr << "        Differencing step parameter erel            : " << erel_ << finl;
          Cerr << "        Differencing step parameter umin            : " << umin_ << finl;
        }
      else
        {
          Cerr << "       === USING ADHOC JFNK SOLVER ===      "   << finl;
          Cerr << "        seuil du residu a convergence               : " << atol_ << finl;
          Cerr << "        dimension de l'espace de Krylov             : " << nkr_ << finl;
          Cerr << "        nombre d'iterations maximal de l'algorithme : " << maxit_gmres_ << finl;
          Cerr << "        residu minimum de convergence               : " << rtol_ << finl;
          Cerr << "        residu maximum de convergence               : " << stol_ << finl;
        }
 
    }
  Cerr << "" << finl;
  Cerr << "********************** Ce message persistera " << tpsaff_ << " secondes *********************" << finl;
  Cerr << "*************** Pour continuer avant les " << tpsaff_ << " secondes : Ctrl + C **************" << finl;
  Cerr << "******************************************************************************" << finl;
  Cerr << "" << finl;
  Cerr << "" << finl;
 
  Nom tps_sleep = "sleep ";
  tps_sleep += Nom(tpsaff_);
  Cerr << (int) system(tps_sleep) << finl;
}
 
void Source_Con_Phase_field::associer_domaines(const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis)
{
  le_dom_VDF_ = ref_cast(Domaine_VDF, domaine_dis);
  le_dom_Cl_VDF_ = ref_cast(Domaine_Cl_VDF, domaine_Cl_dis);
}
 
DoubleTab& Source_Con_Phase_field::ajouter(DoubleTab& resu) const
{
  resu += accr_;
  return resu;
}
 
DoubleTab& Source_Con_Phase_field::calculer(DoubleTab& resu) const
{
  resu = 0.;
  return ajouter(resu);
}
 
void Source_Con_Phase_field::mettre_a_jour(double temps)
{
  Cerr << "Temps : " << temps << " s" << finl;
  Cerr << "" << finl;
 
  Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  Champ_Fonc_base& ch = eq_c.set_mutilde_();
  ch.mettre_a_jour(temps);
}
 
void Source_Con_Phase_field::calculer_u2_elem(DoubleVect& u_carre)
{
  const Navier_Stokes_std& eq_ns = ref_cast(Navier_Stokes_std, pb_->equation(0));
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const DoubleTab& u = eq_ns.inconnue().valeurs();
  const DoubleTab& c = eq_c.inconnue().valeurs();
 
  DoubleVect u2;
 
  u_carre = c;
  u_carre = 0.;
 
  u2 = u;
  u2.carre(VECT_ALL_ITEMS);
 
  // u2 est calcule aux faces. mutildes, pression_thermo et c au centre des elements.
  // On doit donc interpoler u2 pour tout calculer au elements.
 
  // Interpolation de u2, sortie dans u_carre :
  // ------------------------------------------
 
  // Interpolation aux elements de u^2 (a partir de Discretisation_SG_VDF::calculer_tenS_capil)
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  const IntTab& face_sommets = domaine_VDF.face_sommets();
  const DoubleTab& positions = domaine_VDF.xp();
  const Domaine& domaine_geom = domaine_VDF.domaine();
  const int nb_elem = domaine_VDF.nb_elem_tot();
  int f0, f1, som0, som1;
  double psi, val0, val1;
  DoubleTab u2_elem(nb_elem, dimension);
  const int nb_compo_ = u2_elem.line_size();
 
  // Boucle sur le nombre d'elements
  for (int elem = 0; elem < nb_elem; elem++)
    {
      // Boucle sur le nombre de composantes du vecteur
      for (int ncomp = 0; ncomp < nb_compo_; ncomp++)
        {
          f0 = elem_faces(elem, ncomp);
          f1 = elem_faces(elem, dimension + ncomp);
 
          val0 = u2(f0);
          val1 = u2(f1);
 
          som0 = face_sommets(f0, 0);
          som1 = face_sommets(f1, 0);
 
          psi = (positions(elem, ncomp) - domaine_geom.coord(som0, ncomp)) / (domaine_geom.coord(som1, ncomp) - domaine_geom.coord(som0, ncomp));
 
          if (std::fabs(psi) < 1.e-12)
            u2_elem(elem, ncomp) = val0;
          else if (std::fabs(1. - psi) < 1.e-12)
            u2_elem(elem, ncomp) = val1;
          else
            u2_elem(elem, ncomp) = val0 + psi * (val1 - val0);
        }
    }
 
  // Pour chaque element, on additionne les composantes de u2_elem
  double tmp;
  for (int elem = 0; elem < nb_elem; elem++)
    {
      for (int dim = 0; dim < dimension; dim++)
        {
          tmp = u2_elem(elem, dim);
          u_carre(elem) += tmp;
        }
    }
 
  u_carre.echange_espace_virtuel();
 
}
 
/*! @brief Calcul de alpha*(Grad(C))^2 au centre des elements
 *
 * @param (DoubleTab& alpha_gradC_carre) alpha*(Grad(C))^2 au centre des elements
 */
void Source_Con_Phase_field::calculer_alpha_gradC_carre(DoubleTab& alpha_gradC_carre) const
{
  //Cette methode n'est pas utilisee dans cette version systeme multicomposant.(mr264902)
 
  const Navier_Stokes_std& eq_ns = ref_cast(Navier_Stokes_std, pb_->equation(0));
  const Convection_Diffusion_Phase_field& eq_c = ref_cast(Convection_Diffusion_Phase_field, pb_->equation(1));
  const DoubleTab& c = eq_c.inconnue().valeurs();
  const Operateur_Grad& opgrad = eq_ns.operateur_gradient();
  const Milieu_Phase_field& mil = ref_cast(Milieu_Phase_field, pb_->milieu());
 
  DoubleTab& prov_face = ref_cast_non_const(DoubleTab, prov_face_);
  if (prov_face.size() == 0)
    prov_face = eq_ns.inconnue().valeurs();
  prov_face = 0.;
 
  // On calcule Grad(c) que l'on met dans prov_face
  opgrad.calculer(c, prov_face);
  eq_ns.solv_masse().appliquer(prov_face);
 
  // On calcule (Grad(c))^2 sur les faces que l'on met dans gradc2
  DoubleVect gradc2;
  gradc2 = prov_face;
  gradc2.carre();
 
  // Interpolation aux elements de (Grad(c))^2 (a partir de Discretisation_SG_VDF::calculer_tenS_capil)
  const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  const IntTab& face_sommets = domaine_VDF.face_sommets();
  const DoubleTab& positions = domaine_VDF.xp();
  const Domaine& domaine_geom = domaine_VDF.domaine();
  const int nb_elem = domaine_VDF.nb_elem();
  int nb_compo_;
  int elem;
  int f0, f1;
  int som0, som1;
  double psi, val0, val1;
 
  DoubleTab gradc2_elem(nb_elem, dimension);  //***a dimensionner avec le nombre de composants des elements concentration
 
  // Nombre de composantes du vecteur gradc2
  if (gradc2_elem.nb_dim() == 1)
    nb_compo_ = 1;
  else
    nb_compo_ = gradc2_elem.dimension(1);
 
  // Boucle sur le nombre d'elements
  for (elem = 0; elem < nb_elem; elem++)
    {
      // Boucle sur le nombre de composantes du vecteur
      //*** composantes issues de dimension (si 2 donc x et y), faudrait-il mettre prob_dimension a la place de nb_compo et ncomp
      //*** pour eviter toute confusion avec le nombre de composants des elements c1 c2 c3 etc//***
      for (int ncomp = 0; ncomp < nb_compo_; ncomp++)
        {
          f0 = elem_faces(elem, ncomp);
          f1 = elem_faces(elem, dimension + ncomp);
 
          val0 = gradc2(f0);
          val1 = gradc2(f1);
 
          som0 = face_sommets(f0, 0);
          som1 = face_sommets(f1, 0);
 
          psi = (positions(elem, ncomp) - domaine_geom.coord(som0, ncomp)) / (domaine_geom.coord(som1, ncomp) - domaine_geom.coord(som0, ncomp));
 
          if (std::fabs(psi) < 1.e-12)
            gradc2_elem(elem, ncomp) = val0;
          else if (std::fabs(1. - psi) < 1.e-12)
            gradc2_elem(elem, ncomp) = val1;
          else
            gradc2_elem(elem, ncomp) = val0 + psi * (val1 - val0);
        }
    }
 
  // Pour chaque element, on additionne les composantes de gradc2_elem
  int dim;
  double tmp;
  for (elem = 0; elem < nb_elem; elem++)
    {
      for (dim = 0; dim < dimension; dim++)
        {
          tmp = gradc2_elem(elem, dim);
          alpha_gradC_carre(elem) += tmp;
        }
    }
 
  alpha_gradC_carre *= mil.get_fermeture().get_alpha();
  // Ajout B.M suite a plantage dans assert_espace_virtuel_vect dans op_conv...
  alpha_gradC_carre.echange_espace_virtuel();
}
 
/*! @brief Calcul de la pression thermodynamique aux elements
 *
 * @param (DoubleTab& pression_thermo) pression thermodynamique aux elements
 */
void Source_Con_Phase_field::calculer_pression_thermo(DoubleTab& pression_thermo) const
{
  //Cette methode n'est pas utilisee dans cette version systeme multicomposant.(mr264902)
 
  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 DoubleTab& c = eq_c.inconnue().valeurs();
  const DoubleTab& div_alpha_gradC = eq_c.get_div_alpha_gradC();
 
  // Recuperation de la pression en Pascal de l'etape de projection
  //---------------------------------------------------------------
  const DoubleTab& P_Pa = eq_ns.pression_pa().valeurs();
 
  const int taille = pression_thermo.size();
 
  for (int i = 0; i < taille; i++)
    pression_thermo(i) = P_Pa(i) - c(i) * div_alpha_gradC(i);
}
 
void Source_Con_Phase_field::calculer_champ_fonc_c(const double t, Champ_Don_base& champ_fonc_c, const DoubleTab& val_c) const
{
  if (sub_type(Champ_Fonc_Tabule, champ_fonc_c))
    {
      const Champ_Fonc_Tabule& ch_champ_fonc_c = ref_cast(Champ_Fonc_Tabule, champ_fonc_c);
      const Domaine_VDF& domaine_VDF = le_dom_VDF_.valeur();
      const Table& table = ch_champ_fonc_c.table();
      const int isfct = table.isfonction();
      DoubleTab& mes_valeurs = champ_fonc_c.valeurs();
      // code ci-dessous adapte de Champ_Fonc_Tabule_P0_VDF.mettre_a_jour
      if (!(val_c.nb_dim() == mes_valeurs.nb_dim()))
        {
          Cerr << "Erreur a l'initialisation d'un Champ_Fonc_Tabule" << finl;
          Cerr << "Le champ parametre et le champ a initialiser ne sont pas compatibles" << finl;
          Process::exit();
        }
      if (isfct != 2)
        {
          int nb_elems = domaine_VDF.nb_elem();
          if (val_c.line_size() == 1)
            for (int num_elem = 0; num_elem < nb_elems; num_elem++)
              mes_valeurs(num_elem, 0) = table.val(val_c(num_elem));
          else
            {
              int nb_comps = val_c.nb_dim();
              for (int num_elem = 0; num_elem < nb_elems; num_elem++)
                for (int ncomp = 0; ncomp < nb_comps; ncomp++)
                  mes_valeurs(num_elem, ncomp) = table.val(val_c(num_elem, ncomp));
            }
        }
      else
        {
          const DoubleTab& centres_de_gravites = domaine_VDF.xp();
          table.valeurs(val_c, centres_de_gravites, t, mes_valeurs);
        }
    }
}