Schema_Euler_Implicite.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 <Schema_Euler_Implicite.h>
#include <Domaine_Cl_dis_base.h>
#include <Probleme_Couple.h>
#include <communications.h>
#include <Equation_base.h>
#include <Probleme_base.h>
#include <LecFicDiffuse.h>
#include <Milieu_base.h>
#include <Param.h>
#include <string>
 
Implemente_instanciable(Schema_Euler_Implicite,"Schema_euler_implicite|Scheme_euler_implicit",Schema_Implicite_base);
 
Sortie& Schema_Euler_Implicite::printOn(Sortie& s) const
{
  return  Schema_Implicite_base::printOn(s);
}
 
Entree& Schema_Euler_Implicite::readOn(Entree& s)
{
  Schema_Implicite_base::readOn(s);
  if (facsec_max_ == DMAXFLOAT) /* facsec_max not set by the user -> ask the solver for its preferred value */
    facsec_max_ = le_solveur->get_default_facsec_max();
  dt_gf_ = le_solveur->get_default_growth_factor();
  if(!le_solveur)
    {
      Cerr << "A solver must be selected." << finl;
      Cerr << "Syntax : " << finl
           << "Solveur solver_name [ solver parameters ] " << finl;
      exit();
    }
  if (resolution_monolithique_.size() && !le_solveur->est_compatible_avec_th_mono())
    Cerr << que_suis_je() << " : deactivating thermique_monolithique for the implicit solver " << le_solveur->que_suis_je() << finl, resolution_monolithique_.clear();
  if (diffusion_implicite())
    {
      Cerr << "diffusion_implicite option cannot be used with an implicit time scheme." << finl;
      exit();
    }
  return s;
}
 
void Schema_Euler_Implicite::calcul_fac_sec(double& residu,double& residu_old,double& facsec)
{
  if (residu_old==0)
    {
      residu_old=residu;
      nb_ite_sans_accel_=0;
    }
  else if (facsec_func_)
    {
      // facsec=f(t)
      facsec_fn_.setVar(0, temps_courant());
      facsec = facsec_fn_.eval();
      facsec = std::min(facsec,facsec_max_);
    }
  else if ((residu_old>rapport_residus_*residu)||(nb_ite_sans_accel_>nb_ite_sans_accel_max_))
    {
      // Dynamic facsec=f(residual)
      facsec*=sqrt(rapport_residus_);
      residu_old=residu;
      facsec = std::min(facsec,facsec_max_);
      nb_ite_sans_accel_=0;
    }
  nb_ite_sans_accel_++;
}
 
 
int Schema_Euler_Implicite::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  if (mot == "resolution_monolithique")
    {
      Motcle m;
      is >> m;
      assert (m == "{");
      is >> m;
      while (m != "}")
        {
          if (m == "{") // block of equations solved together
            {
              std::set<std::string> bloc_domap;
              is >> m;
              while (m != "}")
                {
                  bloc_domap.insert(m.getString());
                  is >> m;
                }
              resolution_monolithique_.push_back(bloc_domap);
            }
          else
            {
              resolution_monolithique_.push_back({m.getString()});
              is >> m;
            }
        }
    }
  else if (mot == "facsec_max")  //old syntax
    is >> facsec_max_;
  else if (mot == "facsec_expert")
    lire_facsec_expert(is);
  else if (mot == "facsec_func")
    lire_facsec_func(is);
  else return Schema_Implicite_base::lire_motcle_non_standard(mot, is);
  return 1;
}
 
 
// XD facsec_expert interprete nul BRACE To parameter the safety factor for the time step during the simulation.
Entree& Schema_Euler_Implicite::lire_facsec_expert(Entree& is)
{
  Param param("facsec_expert");
  param.ajouter("facsec_ini",&facsec_);               // XD_ADD_P floattant
  // XD_CONT Initial facsec taken into account at the beginning of the simulation.
  param.ajouter("facsec_max",&facsec_max_);           // XD_ADD_P floattant
  // XD_CONT Maximum ratio allowed between time step and stability time returned by CFL condition. The initial ratio
  // XD_CONT given by facsec keyword is changed during the calculation with the implicit scheme but it couldn\'t be
  // XD_CONT higher than facsec_max value.NL2 Warning: Some implicit schemes do not permit high facsec_max, example
  // XD_CONT Schema_Adams_Moulton_order_3 needs facsec=facsec_max=1. NL2 Advice:NL2 The calculation may start with a
  // XD_CONT facsec specified by the user and increased by the algorithm up to the facsec_max limit. But the user can
  // XD_CONT also choose to specify a constant facsec (facsec_max will be set to facsec value then). Faster convergence
  // XD_CONT has been seen and depends on the kind of calculation: NL2-Hydraulic only or thermal hydraulic with forced
  // XD_CONT convection and low coupling between velocity and temperature (Boussinesq value beta low), facsec between
  // XD_CONT 20-30NL2-Thermal hydraulic with forced convection and strong coupling between velocity and temperature
  // XD_CONT (Boussinesq value beta high), facsec between 90-100 NL2-Thermohydralic with natural convection, facsec
  // XD_CONT around 300NL2 -Conduction only, facsec can be set to a very high value (1e8) as if the scheme was
  // XD_CONT unconditionally stableNL2These values can also be used as rule of thumb for initial facsec with a
  // XD_CONT facsec_max limit higher.
  param.ajouter("rapport_residus",&rapport_residus_); // XD_ADD_P floattant
  // XD_CONT Ratio between the residual at time n and the residual at time n+1 above which the facsec is increased by
  // XD_CONT multiplying by sqrt(rapport_residus) (1.2 by default).
  param.ajouter("nb_ite_sans_accel_max",&nb_ite_sans_accel_max_); // XD_ADD_P entier
  // XD_CONT Maximum number of iterations without facsec increases (20000 by default): if facsec does not increase with
  // XD_CONT the previous condition (ration between 2 consecutive residuals too high), we increase it by force after
  // XD_CONT nb_ite_sans_accel_max iterations.
 
  param.lire_avec_accolades(is);
 
  return is;
}
 
void Schema_Euler_Implicite::lire_facsec_func(Entree& is)
{
  Nom facsec_str;
 
  is >> facsec_str;
 
  facsec_fn_.setNbVar(1);
  facsec_fn_.setString(facsec_str);
  facsec_fn_.addVar("t");
  facsec_fn_.parseString();
  facsec_ = facsec_fn_.eval();
  if(facsec_str.majuscule().contient("T"))
    facsec_func_ = true;
}
 
void Schema_Euler_Implicite::set_param(Param& param) const
{
  // XD schema_euler_implicite schema_implicite_base schema_euler_implicite INHERITS_BRACE This is the Euler implicit
  // XD_CONT scheme.
  param.ajouter("max_iter_implicite",&nb_ite_max);
  param.ajouter_flag("facsec_cfl",&facsec_cfl_);    // XD_ADD_P rien
  // XD_CONT Flag to compute time step based on CFL: dt=min(facsec,facsec_max)*dt(convection)X/
  param.ajouter_non_std("facsec_max", (this)); // XD_ADD_P floattant
  // XD_CONT For old syntax, see the complete parameters of facsec for details
  param.ajouter_non_std("facsec_expert", (this)); // XD_ADD_P facsec_expert
  // XD_CONT Advanced facsec specification
  param.ajouter_non_std("facsec_func", (this)); // XD_ADD_P chaine
  // XD_CONT Advanced facsec specification as a function
  param.ajouter_non_std("resolution_monolithique", (this)); // XD_ADD_P bloc_lecture
  // XD_CONT Activate monolithic resolution for coupled problems. Solves together the equations corresponding to the
  // XD_CONT application domains in the given order. All aplication domains of the coupled equations must be given to
  // XD_CONT determine the order of resolution. If the monolithic solving is not wanted for a specific application
  // XD_CONT domain, an underscore can be added as prefix. For example, resolution_monolithique { dom1 { dom2 dom3 }
  // XD_CONT _dom4 } will solve in a single matrix the equations having dom1 as application domain, then the equations
  // XD_CONT having dom2 or dom3 as application domain in a single matrix, then the equations having dom4 as application
  // XD_CONT domain in a sequential way (not in a single matrix).
  Schema_Implicite_base::set_param(param);
}
 
 
bool Schema_Euler_Implicite::initTimeStep(double dt)
{
  Schema_Temps_base::initTimeStep(dt);
  residu_=0;
  return true;
}
////////////////////////////////
//                            //
// Schema characteristics     //
//                            //
////////////////////////////////
 
/*! @brief Returns the number of temporal values to keep.
 *
 */
int Schema_Euler_Implicite::nb_valeurs_temporelles() const
{
  return 3 ;
}
 
/*! @brief Returns the number of future temporal values.
 *
 * Here: n+1, so 1.
 *
 */
int Schema_Euler_Implicite::nb_valeurs_futures() const
{
  return 1 ;
}
 
/*! @brief Returns the time at the i-th future value.
 *
 * Here: t(n+1)
 *
 */
double Schema_Euler_Implicite::temps_futur(int i) const
{
  assert(i==1);
  return temps_courant()+pas_de_temps();
}
 
/*! @brief Returns the time that fields must return when valeurs() is called.
 *
 *     Here: t(n+1)
 *
 */
double Schema_Euler_Implicite::temps_defaut() const
{
  return temps_courant()+pas_de_temps();
}
 
/////////////////////////////////////////
//                                     //
// End of schema characteristics       //
//                                     //
/////////////////////////////////////////
 
void Schema_Euler_Implicite::Initialiser_Champs(Probleme_base& pb)
{
  int i;
  for( i=0; i<pb.nombre_d_equations(); i++)
    {
      DoubleTab& passe = pb.equation(i).inconnue().passe();
      DoubleTab& present = pb.equation(i).inconnue().valeurs();
      if ( (pb.equation(i).que_suis_je()!="Euler") || temps_courant_==0.)
        passe = present;
    }
}
 
int Schema_Euler_Implicite::mettre_a_jour()
{
  calcul_fac_sec(residu_,residu_old_,facsec_);
  return Schema_Temps_base::mettre_a_jour();
}
 
int Schema_Euler_Implicite::Iterer_Pb(Probleme_base& pb,int compteur, int& ok)
{
  int ii;
  bool convergence_pb = true;
  bool convergence_eqn = false;
  int drap = pb.is_dilatable();
  Cout << "=======================================================================================" << finl;
  Cout << "Schema_Euler_Implicite: Implicit iteration " << compteur << " on the " << pb.que_suis_je() << " problem " << pb.le_nom() << " :" << finl;
  Cout << "=======================================================================================" << finl;
  for(int i=0; i<pb.nombre_d_equations(); i++)
    {
      // reverse the order of equations for weakly compressible flow
      if ((drap)&&(i<2))
        ii = 1-i;
      else
        ii = i;
      Equation_base& eqn= pb.equation(ii);
      DoubleTab& present = eqn.inconnue().valeurs();
      DoubleTab& futur = eqn.inconnue().futur();
      double temps=temps_courant_+dt_;
 
      // imposer_cond_lim is used for the pressure and for inter-problem exchanges
      eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
      Cout<<"Solving " << eqn.que_suis_je() << " equation :" << finl;
      const DoubleTab& inut=futur;
      convergence_eqn=le_solveur->iterer_eqn(eqn, inut, present, dt_, compteur, ok);
      if (!ok) return false; // total failure -> exit without processing the following equations
      convergence_pb = convergence_pb&&convergence_eqn;
      futur=present;
      eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
      present=futur;
 
      // The following (commented-out) line performs:
      // Update NS (hence update unknowns)
      // Update turbulence model (k-eps)
      //   eqn.inconnue().mettre_a_jour(temps);
 
      //   eqn.inconnue().reculer();
      eqn.inconnue().Champ_base::changer_temps(temps);
      Cout << finl;
    }
  return (convergence_pb==true);
}
 
void Schema_Euler_Implicite::test_stationnaire(Probleme_base& pb)
{
  for(int i=0; i<pb.nombre_d_equations(); i++)
    {
      // time convergence processing
      DoubleTab& passe = pb.equation(i).inconnue().passe();
      DoubleTab& present = pb.equation(i).inconnue().valeurs();
      DoubleTab& futur = pb.equation(i).inconnue().futur();
      futur = present;
      pb.equation(i).domaine_Cl_dis().imposer_cond_lim(pb.equation(i).inconnue(),temps_courant()+pas_de_temps());
      present -= passe;
      present/=dt_;
      update_critere_statio(present, pb.equation(i));
      present = passe;
    }
}
 
void recommencer_pb(Probleme_base& pb)
{
  for(int i=0; i<pb.nombre_d_equations(); i++)
    {
      DoubleTab& passe = pb.equation(i).inconnue().passe();
      DoubleTab& present = pb.equation(i).inconnue().valeurs();
      present=passe;
    }
}
 
 
bool Schema_Euler_Implicite::iterateTimeStep(bool& converged)
{
  Probleme_base& prob=pb_base();
 
  converged = false;
  Initialiser_Champs(prob);
 
  int ok=1;
  int compteur;
  while (!converged)
    {
      compteur=0;
      Cout<<" "<<finl;
      //Cout<<"Schema_Euler_Implicite : solving of problem "<<prob.que_suis_je()<<finl;
      while (ok && !converged && compteur < nb_ite_max)
        {
          compteur++;
          prob.updateGivenFields();
          converged = Iterer_Pb(prob,compteur, ok);
          Cout<<" "<<finl;;
        }
      if (!ok || (!converged && compteur == nb_ite_max))
        {
          Cout<<"!!! Schema_Euler_Implicite has not converged at t="<< temps_courant_ << " with dt =" << dt_<< " !!!" << finl;
          converged = false;
          notify_failed_timestep(); // to propose a lower time step on the next attempt
          return false;
        }
      else
        {
          Cout<<"The "<<prob.que_suis_je()<<" problem " << prob.le_nom() << " has converged after "<<compteur<<" implicit iterations."<<finl;
        }
    }
  // modification: account for the possibility of modifying dt in Itere_Pb
  // for the compressible solver: divide the time step by 2
  // if convergence is not reached within a given number of time steps
  //  double temps = dt_ + temps_courant_;
  test_stationnaire(prob);
 
  return ok;
}
 
int Schema_Euler_Implicite::faire_un_pas_de_temps_pb_couple(Probleme_Couple& pbc, int& ok)
{
  // Modif B.M.: If the save is done between derivee_en_temps_inco and mettre_a_jour,
  //  a restarted calculation is not equivalent to an uninterrupted one
  //  (especially for front-tracking). So the save is placed at the beginning of the time step.
  //if (lsauv())
  //  for (int i=0;i<pbc.nb_problemes();i++)
  //    ref_cast(Probleme_base,pbc.probleme(i)).sauver();
 
  int convergence_pbc = 0;
  int i;
 
  for(i=0; i<pbc.nb_problemes(); i++)
    Initialiser_Champs(ref_cast(Probleme_base,pbc.probleme(i)));
 
 
  int cv = 0;
  int compteur;
 
  while (!cv)
    {
      compteur=0;
 
      while ((!convergence_pbc )&&(compteur<nb_ite_max))
        {
          compteur++;
          convergence_pbc=1;
 
 
          for(i=0; i<pbc.nb_problemes()*1; i++)
            {
              Probleme_base& pb = ref_cast(Probleme_base,pbc.probleme(i));
              const int nb_eqn=pb.nombre_d_equations();
              for(int ii=0; ii<nb_eqn; ii++)
                {
                  pb.equation(ii).domaine_Cl_dis().calculer_coeffs_echange(temps_courant_+pas_de_temps());
                }
 
            }
 
          if (resolution_monolithique_.size())
            {
              //make sure all application domains are in resolution_monolithique_
              std::set<std::string> doms_app, doms_mono;
              for(i = 0; i < pbc.nb_problemes(); i++)
                for(int j = 0; j < ref_cast(Probleme_base,pbc.probleme(i)).nombre_d_equations(); j++)
                  doms_app.insert(ref_cast(Probleme_base,pbc.probleme(i)).equation(j).domaine_application().getString());
              for (auto && s : resolution_monolithique_)
                for (auto &&d : s) doms_mono.insert((Nom(d)).getSuffix("-").getString());
 
              if (doms_mono != doms_app)
                {
                  Cerr << "Error : all the application domains should be given in the resolution_monolitique block to impose the order of resolution" << finl;
                  Cerr << "and some are missing among :";
                  for (auto &&n: doms_app) Cerr << Nom(" ") + n;
                  Cerr << finl << "(an underscore can be put at the begining of the application domains for which a standard resolution is wanted)" << finl;
                  Process::exit();
                }
 
              for (auto && s : resolution_monolithique_)
                {
                  // serach all equations of this dom app
                  LIST(OBS_PTR(Equation_base)) eqs;
                  for(i = 0; i < pbc.nb_problemes(); i++)
                    for(int j = 0; j < ref_cast(Probleme_base,pbc.probleme(i)).nombre_d_equations(); j++)
                      {
                        const Probleme_base& pb = ref_cast(Probleme_base,pbc.probleme(i));
                        const Motcle type = pb.equation(j).domaine_application();
                        if ((s.count(type.getString()) || s.count((Nom("-") + type).getString())) && !pb.equation(j).equation_non_resolue()) eqs.add(pb.equation(j));
                      }
                  if (eqs.size() == 0) continue; // equations non resolues
 
                  Cout << "RESOLUTION {";
                  for (auto &&d : s) Cout << Nom(" ") + d;
                  Cout << " }" << finl;
                  Cout << "-------------------------" << finl;
                  const bool mono = !(s.size() == 1 && Nom((*s.begin())).debute_par("-"));
                  if (mono)
                    {
                      Cout << "Resolution monolithique! the equations {";
                      for (int k = 0; k < eqs.size(); k++)
                        {
                          Cout << " " << eqs[k]->que_suis_je();
                          eqs[k]->probleme().updateGivenFields();
                        }
                      Cout << " } are solved by assembling a single matrix." << finl;
                      bool convergence_eqs = le_solveur->iterer_eqs(eqs, compteur, ok);
                      convergence_pbc = convergence_pbc && convergence_eqs;
                    }
                  else
                    {
                      for (int k = 0; k < eqs.size(); k++)
                        {
                          Equation_base& eqn = eqs[k].valeur();
                          eqn.probleme().updateGivenFields();
 
                          DoubleTab& present = eqn.inconnue().valeurs();
                          DoubleTab& futur = eqn.inconnue().futur();
                          double temps = temps_courant_ + dt_;
 
                          // imposer_cond_lim is used for the pressure and for inter-problem exchanges
                          eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
                          Cout<<"Solving " << eqn.que_suis_je() << " equation :" << finl;
                          const DoubleTab& inut=futur;
                          bool convergence_eqn=le_solveur->iterer_eqn(eqn, inut, present, dt_, compteur, ok);
                          if (!ok) break;
                          convergence_pbc = convergence_pbc && convergence_eqn;
                          futur = present;
                          eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
                          present = futur;
 
                          eqn.inconnue().Champ_base::changer_temps(temps);
                          Cout << finl;
                        }
                    }
                }
            }
          else
            {
              for(i=0; ok && i<pbc.nb_problemes(); i++)
                {
                  pbc.probleme(i).updateGivenFields();
                  int convergence_pb = Iterer_Pb(ref_cast(Probleme_base,pbc.probleme(i)),compteur, ok);
                  convergence_pbc = convergence_pbc * convergence_pb ;
                }
 
            }
        }
      if (!ok || (!convergence_pbc && compteur==nb_ite_max))
        {
          Cout << pbc.le_nom() << (ok ? " : failure" : " : non-convergence") << " at t = "<< temps_courant_ << " with dt = " << dt_<< " !!!" << finl;
          notify_failed_timestep(); // to propose a lower time step on the next attempt
          return 0;
        }
      else
        {
          Cout<<"The "<<pbc.que_suis_je()<<" problem " << pbc.le_nom() << " has converged after "<<compteur<<" implicit iterations."<<finl;
          cv=1;
        }
    }
 
  double residu_1=0;
  for(i=0; i<pbc.nb_problemes(); i++)
    {
      residu_1=residu_;
      test_stationnaire(ref_cast(Probleme_base,pbc.probleme(i)));
      residu_=std::max(residu_,residu_1);
    }
 
  return 1;
}
 
int Schema_Euler_Implicite::faire_un_pas_de_temps_eqn_base(Equation_base& eqn)
{
  DoubleTab& passe = eqn.inconnue().passe();
  DoubleTab& present = eqn.inconnue().valeurs();
  DoubleTab& futur = eqn.inconnue().futur();
  int compteur = 0, ok = 1;
  bool convergence_eqn = false;
  passe = present;
  // futur=present;
  // used for the pressure and inter-problem couplings
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
  //present=futur;
 
  compteur=0;
  Cout << finl;
  while ((!convergence_eqn)&&(compteur<nb_ite_max))
    {
      compteur++;
      Cout<<"==================================================================================" << finl;
      Cout<<"Schema_Euler_Implicite: Implicit iteration " << compteur << " on the "<<eqn.que_suis_je() << " equation of the problem "<< eqn.probleme().le_nom()<< " :" <<finl;
      Cout<<"==================================================================================" << finl;
      const DoubleTab& inut=futur;
      convergence_eqn=le_solveur->iterer_eqn(eqn, inut, present, dt_, compteur, ok);
      if (!ok) return 0; // total failure
      futur=present;
      eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
      present=futur;
    }
  present -= passe;
  present/=dt_;
 
  update_critere_statio(present, eqn);
 
 
  present = passe;
 
  return 1;
}
 
int Schema_Euler_Implicite::reprendre(Entree&)
{
  Nom nom_fichier(nom_du_cas());
  nom_fichier+=".dt_ev";
  LecFicDiffuse test;
  if (!test.ouvrir(nom_fichier))
    {
      Cerr << "*****************************************" << finl;
      Cerr << "***************** WARNING ***************" << finl;
      Cerr << "File " << nom_fichier << " does not exist." << finl;
      Cerr << "In order to restart a calculation carried out with an implicit time scheme" << finl;
      Cerr << "it is preferable to re-read the .dt_ev file to pick up some informations" << finl;
      Cerr << "and in particular the facsec of the previous calculation." << finl;
      Cerr << "TRUST will use facsec= " << facsec_ << "." << finl;
      Cerr << "Else specify the facsec wanted value in your data file." << finl;
      Cerr << "*****************************************" << finl;
      return 1;
    }
 
  // Reading facsec and residual from the .dt_ev file if it exists
  double facsec_lu=1.;
  double residu_lu=0;
  double facsec_lu_old;
 
  // Test ouverture du fichier
  if (je_suis_maitre())
    {
      EFichier fichier;
      fichier.ouvrir(nom_fichier);
      Nom chaine;
      double temps=0;
      double dt;
      std::string ligne;
      // Read the header if present
      fichier >> chaine;
      if (chaine=="#")
        {
          // Read the full line
          std::getline(fichier.get_ifstream(), ligne);
        } // Otherwise reopen
      else
        {
          fichier.ouvrir(nom_fichier);
        }
      // Search for the time step preceding tinit_
      fichier >> temps;
      double residu_old_old=0;
      while (!fichier.eof() && temps<tinit_)
        {
 
          facsec_lu_old = facsec_lu;
          fichier >> dt;
          fichier >> facsec_lu;
          fichier >> residu_lu;
          if (residu_old_==0)
            residu_old_=residu_lu;
          if (residu_old_old==0)
            residu_old_old=residu_lu;
          // Read the rest of the line
          std::getline(fichier.get_ifstream(), ligne);
          fichier >> temps;
 
          if (facsec_lu_old != facsec_lu)
            {
              residu_old_ = residu_old_old;
 
              nb_ite_sans_accel_ = 1;
            }
          residu_old_old = residu_lu;
          nb_ite_sans_accel_++;
        }
    }
// GF
  calcul_fac_sec(residu_lu,residu_old_,facsec_lu);
  envoyer_broadcast(facsec_lu, 0 /* pe source */);
  envoyer_broadcast(residu_lu, 0 /* pe source */);
  envoyer_broadcast(nb_ite_sans_accel_, 0 /* pe source */);
  envoyer_broadcast(residu_old_, 0 /* pe source */);
 
  // The read facsec is used only if it is within the bounds
  // specified in the data set.
  // Indeed, one may run an explicit calculation and then restart with an implicit scheme.
  residu_ = residu_lu;
  if (facsec_lu>=facsec_)
    {
      if (facsec_lu<=facsec_max_)
        facsec_ = facsec_lu;
      else
        facsec_ = facsec_max_;
      Cerr << "Facsec after reading in " << nom_fichier << " : " << facsec_ << finl;
      Cerr << "Residu after reading in " << nom_fichier << " : " << residu_ << finl;
    }
  else
    {
      Cerr << "The readen facsec in " << nom_fichier << " : " << facsec_lu << finl;
      Cerr << "is not used since it is lower than the facsec from the data set : " << facsec_ << finl;
    }
  return 1;
}
 
int Schema_Euler_Implicite::resolution_monolithique(const Nom& nom) const
{
  for (auto &&s : resolution_monolithique_)
    if (s.count(nom.getString())) return 1;
  return 0;
}