Echange_contact_Correlation_VDF.cpp

/****************************************************************************
* Copyright (c) 2026, 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 <Echange_contact_Correlation_VDF.h>
#include <Domaine_Cl_dis_base.h>
#include <Champ_front_calc.h>
#include <communications.h>
#include <Champ_Uniforme.h>
#include <Probleme_base.h>
#include <Milieu_base.h>
#include <Domaine_VDF.h>
#include <Conduction.h>
#include <Solv_TDMA.h>
#include <SFichier.h>
#include <Param.h>
 
Implemente_instanciable(Echange_contact_Correlation_VDF,"Paroi_Echange_contact_Correlation_VDF",Echange_global_impose);
// XD paroi_echange_contact_correlation_vdf condlim_base paroi_echange_contact_correlation_vdf BRACE Class to define a
// XD_CONT thermohydraulic 1D model which will apply to a boundary of 2D or 3D domain. NL2 Warning : For parallel
// XD_CONT calculation, the only possible partition will be according the axis of the model with the keyword Tranche.
 
Sortie& Echange_contact_Correlation_VDF::printOn(Sortie& s ) const { return s << que_suis_je() << finl; }
 
Entree& Echange_contact_Correlation_VDF::readOn(Entree& is )
{
  if (app_domains.size() == 0) app_domains = { Motcle("Thermique") };
 
  Param param(que_suis_je());
  Reprise_temperature=false;
  dt_impr = 1e10;
  set_param(param);
  param.lire_avec_accolades_depuis(is);
 
  le_champ_front.typer("Champ_front_fonc");
  T_ext().fixer_nb_comp(1);
  h_imp_.typer("Champ_front_fonc");
  h_imp_->fixer_nb_comp(1);
  return is;
}
 
void Echange_contact_Correlation_VDF::set_param(Param& param) const
{
  param.ajouter("dir",&dir); // XD_ADD_P entier
  // XD_CONT Direction (0 : axis X, 1 : axis Y, 2 : axis Z) of the 1D model.
  param.ajouter_condition("(value_of_dir_ge_0)_AND_(value_of_dir_le_2)", "The direction must be 0, 1 or 2 in Echange_contact_Correlation_VDF");
  param.ajouter("Tinf",&Tinf); // XD_ADD_P floattant
  // XD_CONT Inlet fluid temperature of the 1D model (oC or K).
  param.ajouter("Tsup",&Tsup); // XD_ADD_P floattant
  // XD_CONT Outlet fluid temperature of the 1D model (oC or K).
  param.ajouter_non_std("lambda",(this)); // XD_ADD_P chaine
  // XD_CONT Thermal conductivity of the fluid (W.m-1.K-1).
  param.ajouter_non_std("rho",(this)); // XD_ADD_P chaine
  // XD_CONT Mass density of the fluid (kg.m-3) which may be a function of the temperature T.
  param.ajouter("dt_impr",&dt_impr); // XD_ADD_P floattant
  // XD_CONT Printing period in name_of_data_file_time.dat files of the 1D model results.
  param.ajouter("Cp",&Cp); // XD_ADD_P floattant
  // XD_CONT Calorific capacity value at a constant pressure of the fluid (J.kg-1.K-1).
  param.ajouter_non_std("mu",(this)); // XD_ADD_P chaine
  // XD_CONT Dynamic viscosity of the fluid (kg.m-1.s-1) which may be a function of thetemperature T.
  param.ajouter("debit",&debit); // XD_ADD_P floattant
  // XD_CONT Surface flow rate (kg.s-1.m-2) of the fluid into the channel.
  param.ajouter("Dh",&diam); // XD_ADD_P floattant
  // XD_CONT Hydraulic diameter may be a function f(x) with x position along the 1D axis (xinf <= x <= xsup)
  param.ajouter_non_std("volume",(this)); // XD_ADD_P chaine
  // XD_CONT Exact volume of the 1D domain (m3) which may be a function of the hydraulic diameter (Dh) and the lateral
  // XD_CONT surface (S) of the meshed boundary.
  param.ajouter_non_std("Nu",(this)); // XD_ADD_P chaine
  // XD_CONT Nusselt number which may be a function of the Reynolds number (Re) and the Prandtl number (Pr).
  // Added by Cyril MALOD (14/09/2006)
  param.ajouter_flag("Reprise_correlation",&Reprise_temperature); // XD_ADD_P rien
  // XD_CONT Keyword in the case of a resuming calculation with this correlation.
}
 
int Echange_contact_Correlation_VDF::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  int retval = 1;
 
  if (mot=="lambda")
    {
      Nom tmp;
      is >> tmp;
      lambda_T.setNbVar(1);
      lambda_T.setString(tmp);
      lambda_T.addVar("T");
      lambda_T.parseString();
    }
  else if (mot=="rho")
    {
      Nom tmp;
      is >> tmp;
      rho_T.setNbVar(1);
      rho_T.setString(tmp);
      rho_T.addVar("T");
      rho_T.parseString();
    }
  else if (mot=="mu")
    {
      Nom tmp;
      is >> tmp;
      mu_T.setNbVar(1);
      mu_T.setString(tmp);
      mu_T.addVar("T");
      mu_T.parseString();
    }
  else if (mot=="volume")
    {
      Nom tmp;
      is >> tmp;
      fct_vol.setNbVar(2);
      fct_vol.setString(tmp);
      fct_vol.addVar("Dh");
      fct_vol.addVar("S");
      fct_vol.parseString();
    }
  else if (mot=="Nu")
    {
      Nom tmp;
      is >> tmp;
      fct_Nu.setNbVar(2);
      fct_Nu.setString(tmp);
      fct_Nu.addVar("Re");
      fct_Nu.addVar("Pr");
      fct_Nu.parseString();
    }
  else retval = -1;
 
  return retval;
}
 
 
/**
 * @brief Computes the local heat exchange coefficient in the solid cell.
 */
void Echange_contact_Correlation_VDF::calculer_h_mon_pb(DoubleTab& tab)
{
  // necessarily local
  const Equation_base& mon_eqn = domaine_Cl_dis().equation();
  const Milieu_base& mon_milieu = mon_eqn.milieu();
  const Domaine_VDF& ma_zvdf = ref_cast(Domaine_VDF,domaine_Cl_dis().domaine_dis());
  const Front_VF& ma_front_vf = ref_cast(Front_VF,frontiere_dis());
  calculer_h_solide(tab,mon_eqn,ma_zvdf,ma_front_vf,mon_milieu);
}
 
 
 
/**
 * @brief Completes and initializes the class attributes.
 */
void Echange_contact_Correlation_VDF::completer()
{
  Echange_impose_base::completer();
  trier_coord();
  T=std::min(Tinf,Tsup);
  calculer_prop_physique();
  if (T(0)==Tinf)
    U=debit/rho(0);
  else U = debit/rho(N-1);
 
  for (int i=0; i<N; i++)
    h_correlation(i) = calculer_coefficient_echange(i);
 
  DoubleTab& mon_h= h_imp_->valeurs();
  // calculer_h_mon_pb(mon_h);
 
 
  const Front_VF& ma_front_vf = ref_cast(Front_VF,frontiere_dis());
  const int nb_faces_bord = ma_front_vf.nb_faces();
  DoubleTab& Text_valeurs = T_ext().valeurs();
  Text_valeurs.resize(nb_faces_bord,1);
 
  mon_h.resize(nb_faces_bord,1);
 
  for (int ii=0; ii<nb_faces_bord; ii++)
    {
      mon_h(ii,0) = h_correlation(correspondance_solide_fluide(ii));
      Text_valeurs(ii,0) = T(correspondance_solide_fluide(ii));
    }
 
  init_tab_echange();
 
  //Update because there is no dependency on data from another problem
  //in the case of coupled problems.
  //The completer() method may need to be replaced by the initialiser() method.
  Equation_base& mon_eqn = domaine_Cl_dis().equation();
  const double temps = mon_eqn.schema_temps().temps_courant();
  mettre_a_jour(temps);
}
 
/**
 * @brief Initializes the tab_ech array for parallel execution.
 */
void Echange_contact_Correlation_VDF::init_tab_echange()
{
  DoubleVect y_envoye(2);
  DoubleVect y_recu(2);
  y_envoye(0) = coord(1); // y_min
  y_envoye(1) = coord(N-2); // y_max
 
 
  const Joints& joints = domaine_Cl_dis().domaine().faces_joint();
  const int nb_voisins = joints.size();
 
  const int ME=Process::me();
  tab_ech.resize(Process::nproc(),Process::nproc());
  tab_ech=-1;
  for (int i=0; i<nb_voisins; i++)
    {
      int mon_voisin = joints[i].PEvoisin();
      Cerr << " Echange_contact_Correlation_VDF::init_tab_echange() To check for blocking send" << finl;
      exit();
      envoyer(y_envoye,ME,mon_voisin,ME);
 
      recevoir(y_recu,mon_voisin,ME,mon_voisin);
      if (y_recu(1)<coord(1))
        {
          tab_ech(mon_voisin,ME)=1; // mon_voisin m'enverra le max
        }
      else if (y_recu(0)>coord(N-2))
        {
          tab_ech(mon_voisin,ME)=0; // mon_voisin m'enverra le min
        }
    }
}
 
 
/**
 * @brief Computes rho, mu and lambda of the fluid at the current temperature.
 */
void Echange_contact_Correlation_VDF::calculer_prop_physique()
{
  for (int i=0; i<N; i++)
    {
      rho_T.setVar("T",T(i));
      mu_T.setVar("T",T(i));
      lambda_T.setVar("T",T(i));
      rho(i) = rho_T.eval();
      mu(i) = mu_T.eval();
      lambda(i) = lambda_T.eval();
    }
}
 
/**
 * @brief Computes the heat exchange coefficient using the correlation entered in the data file.
 */
double Echange_contact_Correlation_VDF::calculer_coefficient_echange(int i)
{
  double Re,Pr;
  Re = std::fabs(getQh()*getDh()/getMu(i));
  Pr = getMu(i)*getCp()/getLambda(i);
  fct_Nu.setVar("Re",Re);
  fct_Nu.setVar("Pr",Pr);
  return fct_Nu.eval()*getLambda(i)/getDh();
}
 
/**
 * @brief Computes the volumetric power source term in the 1D fluid energy equation.
 */
void Echange_contact_Correlation_VDF::calculer_Q()
{
  const Domaine_VDF& ma_zvdf = ref_cast(Domaine_VDF,domaine_Cl_dis().domaine_dis());
  const Front_VF& ma_front_vf = ref_cast(Front_VF,frontiere_dis());
  const int ndeb = ma_front_vf.num_premiere_face();
  const int nb_faces_bord = ma_front_vf.nb_faces();
  const DoubleVect& surfaces = ma_zvdf.face_surfaces();
  const IntTab& face_voisins = ma_zvdf.face_voisins();
  Champ_Inc_base& inco = domaine_Cl_dis().equation().inconnue();
  DoubleTab& Ts = inco.valeurs();
  DoubleTab& mon_h= h_imp_->valeurs();
 
 
  Qvol=0.;
  for(int iface=0; iface<nb_faces_bord; iface++)
    {
      int corresp = correspondance_solide_fluide(iface);
      int elem =  face_voisins(ndeb+iface,0);
      if(elem == -1)
        elem =  face_voisins(ndeb+iface,1);
      Qvol(corresp)+=mon_h(iface,0)*(Ts(elem)-T(corresp))*surfaces(ndeb+iface);
    }
  for (int i=0; i<N; i++)
    Qvol(i)/=vol(i);
}
 
 
 
void Echange_contact_Correlation_VDF::trier_coord()
{
  const Domaine_VDF& ma_zvdf = ref_cast(Domaine_VDF,domaine_Cl_dis().domaine_dis());
  const DoubleVect& surfaces = ma_zvdf.face_surfaces();
  const Front_VF& ma_front_vf = ref_cast(Front_VF,frontiere_dis());
 
  const int ndeb = ma_front_vf.num_premiere_face();
  const int nb_faces_bord = ma_front_vf.nb_faces();
 
  IntVect face_triee(nb_faces_bord);
  const DoubleTab& xv = ma_zvdf.xv();
  int i, j ,tmp;
 
  // Sort faces by increasing coordinate (along the coordinate corresponding to the tube direction)
  if (nb_faces_bord>1)
    {
      for(i=0; i<nb_faces_bord; i++)
        face_triee(i) = ndeb+i;
 
      for(i=1; i<nb_faces_bord; i++)
        {
          j=i-1;
          tmp = face_triee(i);
          while ( j>=0 && xv(face_triee(j),dir)>xv(tmp,dir))
            {
              face_triee(j+1) = face_triee(j);
              j--;
            }
          face_triee(j+1) = tmp;
        }
    }
  // Determine the number of 1D cells in the fictitious fluid
  // and fill the solid-to-fluid element correspondence array
  correspondance_solide_fluide.resize(nb_faces_bord);
  N=3; // add the boundary nodes
  correspondance_solide_fluide(face_triee(0)-ndeb) = 1;
  for(i=1; i<nb_faces_bord; i++)
    {
      if (xv(face_triee(i),dir) != xv(face_triee(i-1),dir))
        N++;
      correspondance_solide_fluide(face_triee(i)-ndeb) = N-2;
    }
 
  U.resize(N);
  T.resize(N);
  Qvol.resize(N);
  rho.resize(N);
  mu.resize(N);
  lambda.resize(N);
  vol.resize(N);
  coord.resize(N);
  h_correlation.resize(N);
 
 
 
  // Compute slice volumes
  DoubleVect surf(N);
  surf=0.;
  vol=0.;
  surf(0) = surf(N-1) = 0.; // not used
  for(int iface=0; iface<nb_faces_bord; iface++)
    {
      int corresp = correspondance_solide_fluide(iface);
      surf(corresp)+=surfaces(ndeb+iface);
      coord(corresp) = xv(ndeb+iface,dir);
    }
  coord(0) = coord(1)-0.5*(coord(2)-coord(1)); // the first mesh point
  coord(N-1) = coord(N-2)+0.5*(coord(N-2)-coord(N-3)); // the last mesh point
  for (i=1; i<N-1; i++)
    {
      vol(i)=volume(surf(i),diam);
    }
  vol(0) = vol(N-1) = 1.; // not used, but avoids a division by zero in the Qvol computation.
}
 
 
 
/**
 * @brief Computes the volume of a slice with lateral surface area s and hydraulic diameter d.
 * @param s Lateral surface area of the slice.
 * @param d Hydraulic diameter.
 * @return Volume of the slice.
 */
double Echange_contact_Correlation_VDF::volume(double s, double d)
{
  fct_vol.setVar("Dh",d);
  fct_vol.setVar("S",s);
 
  return fct_vol.eval();
}
 
 
/**
 * @brief Computes the boundary conditions to apply to the energy equation.
 */
void Echange_contact_Correlation_VDF::calculer_CL()
{
  const int ME=Process::me();
 
  const Joints& joints = domaine_Cl_dis().domaine().faces_joint();
  const int nb_voisins = joints.size();
 
  DoubleVect les_cl_envoyees(2);
  DoubleVect les_cl_recues(2);
  les_cl_envoyees(0) = T(1);
  les_cl_envoyees(1) = T(N-2);
 
  T_CL0=Tinf;
  T_CL1=Tsup;
  for (int i=0; i<nb_voisins; i++)
    {
      int mon_voisin = joints[i].PEvoisin();
 
      Cerr << " Echange_contact_Correlation_VDF::calculer_CL() To check for blocking send" << finl;
      exit();
      envoyer(les_cl_envoyees,ME,mon_voisin,ME);
 
      recevoir(les_cl_recues,mon_voisin,ME,mon_voisin);
      if (tab_ech(mon_voisin,ME) == 1)
        {
          T_CL0 = 0.5*(les_cl_recues(1)+T(1));
          T(0) = T_CL0;
        }
      else if (tab_ech(mon_voisin,ME) == 0)
        {
          T_CL1 = 0.5*(les_cl_recues(0)+T(N-2));
          T(N-1) = T_CL1;
        }
      else
        {
          Cerr << "Communication error between processors " <<ME<< " and " <<mon_voisin << finl;
          Cerr << "in Echange_contact_Correlation_VDF::calculer_CL" << finl;
          if (Process::is_parallel()) Cerr << "Verify that your partition is in slices along the axis of the 1D model." << finl;
          exit();
 
        }
    }
}
 
 
/**
 * @brief Computes the velocity from mass conservation.
 */
void Echange_contact_Correlation_VDF::calculer_Vitesse()
{
  for (int i=0; i<N; i++)
    U(i) = debit/rho(i);
}
 
 
/**
 * @brief Computes the 1D temperature in the fluid by solving energy conservation.
 */
void Echange_contact_Correlation_VDF::calculer_Tfluide()
{
  const Equation_base& mon_eqn = domaine_Cl_dis().equation();
  const double dt = mon_eqn.schema_temps().pas_de_temps();
  DoubleVect ma(N); // diagonal
  DoubleVect mb(N-1); // sub-diagonal
  DoubleVect mc(N-1); // super-diagonal
  DoubleVect sm(N);
  const int sgn = (debit>0) ? 1 : -1; // upwind scheme for transport
 
 
  ma(0) = 1.;
  ma(N-1) = 1.;
  sm(0) = T_CL0;
  sm(N-1) = T_CL1;
  mc(0) = 0.;
  mb(N-2) = 0.;
  int i;
  for (i=1; i<N-1; i++)
    {
      const double dtrhoCp = dt/rho(i)/Cp;
      const double dz1 = coord(i+1)-coord(i);
      const double dz2 = coord(i)-coord(i-1);
      const double l1 = 0.5*(lambda(i+1)+lambda(i));
      const double l2 = 0.5*(lambda(i)+lambda(i-1));
 
      ma(i)=1+dtrhoCp*(sgn*Cp*debit/dz1+2*(l1/dz1+l2/dz2)/(dz1+dz2));
 
      sm(i) = T(i)+dtrhoCp*Qvol(i);
 
      mc(i) = dtrhoCp*(0.5*(1-sgn)*Cp*debit/dz1-2*l1/dz1/(dz1+dz2));
 
      mb(i-1) = dtrhoCp*(-2*l2/dz2/(dz1+dz2));
      mb(i-1) += dtrhoCp*(-0.5*(1+sgn)*Cp*debit/dz1);
    }
 
 
 
 
  DoubleVect xx(N);
  xx=0.;
  Solv_TDMA::resoudre(ma,mb,mc,sm,xx,N);
  for (i=0; i<N; i++)
    T(i) = xx(i);
}
 
 
 
 
/**
 * @brief Updates the velocity, fluid temperature and heat exchange coefficient.
 */
void Echange_contact_Correlation_VDF::mettre_a_jour(double temps)
{
 
  // Name of the backup files
  Nom Fichier_sauv_nom=domaine_Cl_dis().equation().probleme().le_nom();
  Fichier_sauv_nom+="_";
  Fichier_sauv_nom+=frontiere_dis().frontiere().le_nom();
  Fichier_sauv_nom+=".sauv";
 
  // Resume operation for the fluid temperature field
  const int ME = Process::me();
  const int nbproc = Process::nproc();
  FILE *Fichier_sauv;
  if (Reprise_temperature)
    {
      if (nbproc>1)
        {
          Cerr << "To check for blocking send" << finl;
          exit();
          envoyer(T,ME,0,ME);                                                                // send T to know its array dimension
          if (je_suis_maitre())
            {
              double temps_lu=0., Temperature_reprise=0.;
              char tmp_lu;
              //int code;
              Fichier_sauv=fopen(Fichier_sauv_nom,"r");                                // open the backup file in read-only mode
              if (Fichier_sauv!=nullptr)
                {
                  if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                                // read the first word "Temps" (discarded)
                  if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                                // read the second word "=" (discarded)
                  if (!fscanf(Fichier_sauv,"%lf",&temps_lu)) exit();                                // read the time and store it in temps_lu
                  if (temps_lu==temps)
                    {
                      for (int p=0; p<nbproc; p++)                                      // loop over the number of processors
                        {
                          Cerr << "Resuming temperature field (" << Fichier_sauv_nom << ") for the correlation on processor " << p << "...";
                          DoubleVect T_tmp;
                          recevoir(T_tmp,p,0,p);                                        // retrieve T to know its array dimension at the master processor level
                          for (int i=0; i<T_tmp.size(); i++)                      // loop over the number of temperatures to define
                            {
                              if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                // read the word "T(i)" (discarded)
                              if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                // read the word "=" (discarded)
                              if (!fscanf(Fichier_sauv,"%lf",&Temperature_reprise)) exit();// read the temperature value and store it in Temperature_reprise
                              T_tmp(i) = Temperature_reprise;                        // assign the computation variable T(i) to Temperature_reprise
                            }
                          envoyer(T_tmp,0,p,0);                                        // send T_tmp to all processors
                          Cerr << "   OK!" << finl;
                        }
                      fclose(Fichier_sauv);
                    }
                  else
                    {
                      Cerr << "\nThe time indicated in the file \""<< Fichier_sauv_nom << "\" "<< temps_lu << finl;
                      Cerr << "differs from the restart time of the computation "<< temps << " !!!" << finl << finl;
                      Cerr << "You cannot resume with this correlation. Remove" << finl;
                      Cerr << "the keyword \"Reprise\" from the correlation. The transient phase" << finl;
                      Cerr << "of the computation will be wrong, but the steady state will be correct." << finl << finl;
                      fclose(Fichier_sauv);
                      exit();
                    }
                }
              else
                {
                  Cerr << "\nThe resume cannot be performed: the file \"" << Fichier_sauv_nom << "\" is missing !!!" << finl << finl;
                  exit();
                }
            }
          recevoir(T,0,ME,0);                                                                // retrieve T=T_tmp at each processor level
        }
      else
        {
          Fichier_sauv=fopen(Fichier_sauv_nom,"r");
          double temps_lu=0., Temperature_reprise=0.;
          char tmp_lu;
 
          if (Fichier_sauv!=nullptr)
            {
              if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                        // read the first word "Temps" (discarded)
              if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                        // read the second word "=" (discarded)
              if (!fscanf(Fichier_sauv,"%lf",&temps_lu)) exit();                        // read the time and store it in temps_lu
              if (temps_lu==temps)
                {
                  Cerr << "Resuming temperature field (" << Fichier_sauv_nom << ") for the correlation...";
                  for (int i=0; i<N; i++)
                    {
                      if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                // read the word "T(i)" (discarded)
                      if (!fscanf(Fichier_sauv,"%s",&tmp_lu)) exit();                // read the word "=" (discarded)
                      if (!fscanf(Fichier_sauv,"%lf",&Temperature_reprise)) exit();// read the temperature value and store it in Temperature_reprise
                      T(i) = Temperature_reprise;                        // assign the computation variable T(i) to Temperature_reprise
                    }
                  fclose(Fichier_sauv);
                  Cerr << "   OK!" << finl;
                }
              else
                {
                  Cerr << "\nThe time indicated in the file \""<< Fichier_sauv_nom << "\" "<< temps_lu << finl;
                  Cerr << "differs from the restart time of the computation "<< temps << " !!!" << finl << finl;
                  Cerr << "You cannot resume with this correlation. Remove" << finl;
                  Cerr << "the keyword \"Reprise\" from the correlation. The transient phase" << finl;
                  Cerr << "of the computation will be wrong, but the steady state will be correct." << finl << finl;
                  fclose(Fichier_sauv);
                  exit();
                }
            }
          else
            {
              Cerr << "\nThe resume cannot be performed: the file \"" << Fichier_sauv_nom << "\" is missing !!!" << finl << finl;
              exit();
            }
        }
      Reprise_temperature=false;
    }
 
  calculer_CL();
  calculer_prop_physique();
  calculer_Q();
  calculer_Vitesse();
  calculer_Tfluide();
 
 
 
  for (int i=0; i<N; i++)
    h_correlation(i) = calculer_coefficient_echange(i);
 
 
 
 
  DoubleTab& mon_h= h_imp_->valeurs();
  calculer_h_mon_pb(mon_h);
  const int taille=mon_h.dimension(0);
  DoubleTab& Text_valeurs = T_ext().valeurs();
 
  for (int ii=0; ii<taille; ii++)
    {
      mon_h(ii,0) = 1./(1./mon_h(ii,0)+1./h_correlation(correspondance_solide_fluide(ii)));
      Text_valeurs(ii,0) = T(correspondance_solide_fluide(ii));
    }
 
  const Equation_base& mon_eqn = domaine_Cl_dis().equation();
 
  if (limpr(temps,mon_eqn.schema_temps().pas_de_temps())) imprimer(temps);
 
  Echange_global_impose::mettre_a_jour(temps);
 
  // Save operation for the fluid temperature field for a future restart
  if (nbproc>1)
    {
      envoyer(T,ME,0,ME);                                                                // send T to the master processor
      if (je_suis_maitre())
        {
          Fichier_sauv=fopen(Fichier_sauv_nom,"w");                // open the backup file in write mode
          fprintf(Fichier_sauv,"Temps\t=\t%f\n",temps);                                // print the time
          int j=0;
          for (int p=0; p<nbproc; p++)                                              // loop over the number of processors
            {
              DoubleVect T_tmp;
              recevoir(T_tmp,p,0,p);                                                // retrieve temperatures from each processor
              for (int i=0; i<T_tmp.size(); i++)
                {
                  fprintf(Fichier_sauv,"T(%i)\t=\t%f\n",(int)j,T_tmp(i));        // print the temperatures in the backup file
                  j=j+1;
                }
            }
          fclose(Fichier_sauv);
        }
    }
  else
    {
      Fichier_sauv=fopen(Fichier_sauv_nom,"w");
      fprintf(Fichier_sauv,"Temps\t=\t%f\n",temps);
      for (int i=0; i<N; i++)
        {
          fprintf(Fichier_sauv,"T(%i)\t=\t%f\n",(int)i,T(i));
        }
      fclose(Fichier_sauv);
    }
}
 
void Echange_contact_Correlation_VDF::calculer_h_solide(DoubleTab& tab,const Equation_base& une_eqn,const Domaine_VDF& zvdf_2,const Front_VF& front_vf,const    Milieu_base& le_milieu)
{
  DoubleVect e;
  const IntTab& face_voisins = zvdf_2.face_voisins();
  int i;
  int nb_comp = le_milieu.conductivite().nb_comp();
  int ndeb = front_vf.num_premiere_face();
  int nfin = ndeb + front_vf.nb_faces();
 
  e.resize(front_vf.nb_faces());
 
  for (int face=ndeb; face<nfin; face++)
    e(face-ndeb) = zvdf_2.dist_norm_bord(face);
 
  // Compute tab = 1/(e/lambda + 1/h_paroi) =1/(e/lambda+invhparoi)
  if(!sub_type(Champ_Uniforme,le_milieu.conductivite()))
    {
      //Cerr << "raccord local homogene et conductivite non uniforme" << finl;
      const DoubleTab& tab_lambda = le_milieu.conductivite().valeurs();
      for (int face=ndeb; face<nfin; face++)
        {
          int elem = face_voisins(face,0);
          if (elem == -1)
            elem = face_voisins(face,1);
          for(i=0; i<nb_comp; i++)
            {
              assert(le_milieu.conductivite().valeurs()(elem,i)!=0.);
              tab(face-ndeb,i) = tab_lambda(elem,i)/e(face-ndeb);
            }
        }
    }
  else  // the conductivity is a uniform OWN_PTR(Champ_base)
    {
      for (int face=ndeb; face<nfin; face++)
        {
          for(i=0; i<nb_comp; i++)
            {
              assert(le_milieu.conductivite().valeurs()(0,i)!=0.);
              tab(face-ndeb,i) = le_milieu.conductivite().valeurs()(0,i)/e(face-ndeb);
            }
        }
    }
}
 
 
/**
 * @brief Tests whether printing is requested.
 */
int Echange_contact_Correlation_VDF::limpr(double temps_courant,double dt) const
{
  const Schema_Temps_base& sch = domaine_Cl_dis().equation().schema_temps();
  // Slightly convoluted test, but it works !!!! (CM 05/07/2007)
  if (dt_impr<=dt || ((sch.temps_max()<=temps_courant || sch.nb_pas_dt_max()<=(sch.nb_pas_dt()+1) || (temps_courant!=sch.temps_courant() && sch.nb_pas_dt()==0)) && dt_impr!=1e10))
    return 1;
  else
    {
      // See Schema_Temps_base::limpr for information on epsilon and modf
      static const double epsilon = 1.e-9;
      double i, j;
      modf(temps_courant/dt_impr + epsilon, &i);
      modf((temps_courant-dt)/dt_impr + epsilon, &j);
      return ( i>j );
    }
}
 
 
/**
 * @brief Prints the results.
 */
void Echange_contact_Correlation_VDF::imprimer(double temps) const
{
  const int ME = Process::me();
  const int nbproc = Process::nproc();
 
  if (nbproc>1)
    {
      envoyer(coord,ME,0,ME);
      envoyer(T,ME,0,ME);
      envoyer(U,ME,0,ME);
      envoyer(h_correlation,ME,0,ME);
      envoyer(rho,ME,0,ME);
      envoyer(mu,ME,0,ME);
      envoyer(lambda,ME,0,ME);
      envoyer(Qvol,ME,0,ME);
      envoyer(vol,ME,0,ME);
 
      if (je_suis_maitre())
        {
          Nom nom_bord=frontiere_dis().frontiere().le_nom();
          nom_bord+="_";
          nom_bord+=Nom(temps);
          nom_bord+=".dat";
          SFichier fic(nom_bord);
          fic.precision(domaine_Cl_dis().equation().schema_temps().precision_impr());
          fic.setf(ios::scientific);
          double Qt=0.;
          fic << "# X                T                U                h                rho                mu                lambda                        Q[W]" << finl;
 
          for (int p=0; p<nbproc; p++)
            {
              DoubleVect coord_tmp;
              DoubleVect T_tmp;
              DoubleVect U_tmp;
              DoubleVect h_tmp;
              DoubleVect rho_tmp;
              DoubleVect mu_tmp;
              DoubleVect lambda_tmp;
              DoubleVect Qvol_tmp;
              DoubleVect vol_tmp;
 
              recevoir(coord_tmp,p,0,p);
              recevoir(T_tmp,p,0,p);
              recevoir(U_tmp,p,0,p);
              recevoir(h_tmp,p,0,p);
              recevoir(rho_tmp,p,0,p);
              recevoir(mu_tmp,p,0,p);
              recevoir(lambda_tmp,p,0,p);
              recevoir(Qvol_tmp,p,0,p);
              recevoir(vol_tmp,p,0,p);
 
              for (int i =0; i<coord_tmp.size(); i++)
                {
                  fic << coord_tmp(i) << " \t" << T_tmp(i) << " \t" << U_tmp(i) << " \t" << h_tmp(i) << " \t" << rho_tmp(i);
                  fic << " \t" << mu_tmp(i) << " \t" << lambda_tmp(i) << " \t" << Qvol_tmp(i)*vol_tmp(i) << finl;
                  Qt+=Qvol_tmp(i)*vol_tmp(i);
                }
 
            }
 
 
          fic << "# Q total[W] = " << Qt << finl;
          fic << finl;
        }
    }
  else
    {
      Nom nom_bord=frontiere_dis().frontiere().le_nom();
      nom_bord+="_";
      nom_bord+=Nom(temps);
      nom_bord+=".dat";
      SFichier fic(nom_bord);
      fic.precision(domaine_Cl_dis().equation().schema_temps().precision_impr());
      fic.setf(ios::scientific);
      double Qt=0.;
      fic << "# X                 T                 U                 h                 rho                 mu                 lambda                  Q[W]" << finl;
      for (int i =0; i<N; i++)
        {
          fic << coord(i) << " \t" << T(i) << " \t" << U(i) << " \t" << h_correlation(i) << " \t" << rho(i);
          fic << " \t" << mu(i) << " \t" << lambda(i) << " \t" << Qvol(i)*vol(i) << finl;
          Qt+=Qvol(i)*vol(i);
        }
      fic << "# Q total[W] = " << Qt << finl;
      fic << finl;
 
    }
}