Schema_Implicite_Multi_TimeStep_base.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_Implicite_Multi_TimeStep_base.h>
#include <Domaine_Cl_dis_base.h>
#include <Probleme_Couple.h>
#include <communications.h>
#include <Matrice_Morse.h>
#include <Equation_base.h>
#include <Probleme_base.h>
#include <LecFicDiffuse.h>
#include <Milieu_base.h>
#include <Param.h>
 
Implemente_base(Schema_Implicite_Multi_TimeStep_base,"Schema_Implicite_Multi_TimeStep_base",Schema_Implicite_base);
 
Sortie& Schema_Implicite_Multi_TimeStep_base::printOn(Sortie& s) const
{
  return  Schema_Implicite_base::printOn(s);
}
 
Entree& Schema_Implicite_Multi_TimeStep_base::readOn(Entree& s)
{
  Schema_Implicite_base::readOn(s);
 
  if(!le_solveur)
    {
      Cerr << "A solver must be selected." << finl;
      Cerr << "Syntax : " << finl
           << "Solveur solver_name [ solver parameters ] " << finl;
      exit();
    }
  if (diffusion_implicite())
    {
      Cerr << "diffusion_implicite option cannot be used with an implicit time scheme." << finl;
      exit();
    }
 
  return s;
}
 
static void Implicite_Multi_TimeStep_calcul_fac_sec(double& residu_,double& residu_old,double& facsec_,const double facsec_max_,int& nb_ite_sans_accel_)
{
  double rap                   = 1.2;
  int    nb_ite_sans_accel_max = 20000;
 
  if (residu_old==0)
    {
      residu_old=residu_;
      nb_ite_sans_accel_=0;
    }
  else if ((residu_old>rap*residu_)||(nb_ite_sans_accel_>nb_ite_sans_accel_max))
    {
      facsec_*=sqrt(rap);
      residu_old=residu_;
      facsec_=std::min(facsec_,facsec_max_);
      nb_ite_sans_accel_=0;
    }
 
  nb_ite_sans_accel_++;
}
 
void Schema_Implicite_Multi_TimeStep_base::set_param(Param& param) const
 
{
  param.ajouter("max_iter_implicite",&nb_ite_max);
  param.ajouter("facsec_max", &facsec_max_);
  Schema_Implicite_base::set_param(param);
}
 
bool Schema_Implicite_Multi_TimeStep_base::initTimeStep(double dt)
{
  Schema_Implicite_base::initTimeStep(dt);
  residu_ = 0;
  return true;
}
////////////////////////////////
//                            //
// Schema characteristics     //
//                            //
////////////////////////////////
 
/*! @brief Returns the number of future temporal values.
 *
 * Here: n+1, so 1.
 *
 */
int Schema_Implicite_Multi_TimeStep_base::nb_valeurs_futures() const
{
  return 1 ;
}
 
/*! @brief Returns the time at the i-th future value.
 *
 * Here: t(n+1)
 *
 */
double Schema_Implicite_Multi_TimeStep_base::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_Implicite_Multi_TimeStep_base::temps_defaut() const
{
  return temps_courant()+pas_de_temps();
}
 
/////////////////////////////////////////
//                                     //
// End of schema characteristics       //
//                                     //
/////////////////////////////////////////
 
void Schema_Implicite_Multi_TimeStep_base::Initialiser_Champs(Probleme_base& pb)
{
  int i      = 0;
  int j      = 0;
 
  for(i=0; i<pb.nombre_d_equations(); i++)
    {
      if ( (pb.equation(i).que_suis_je()!="Euler") || temps_courant_==0.)
        {
          for (j=nb_valeurs_passees(); j>-1; j--)
            {
              pb.equation(i).inconnue().passe(j+1) = pb.equation(i).inconnue().passe(j);
            }
        }
    }
}
 
 
 
int Schema_Implicite_Multi_TimeStep_base::mettre_a_jour()
{
  Implicite_Multi_TimeStep_calcul_fac_sec(residu_,residu_old_,facsec_,facsec_max_,nb_ite_sans_accel_);
  return Schema_Temps_base::mettre_a_jour();
}
 
int Schema_Implicite_Multi_TimeStep_base::Iterer_Pb(Probleme_base& pb,int ite)
{
  int  i               = 0, ok = 1;
  int  ii              = 0;
  bool convergence_pb  = true;
  bool convergence_eqn = false;
  int  drap = pb.is_dilatable();
 
  for(i=0; i<pb.nombre_d_equations(); i++)
    {
      //reverse equation order if quasi-compressible medium
      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_;
 
      DoubleTab      stored_parameters;
 
      authorized_equation(eqn);
 
      // 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 equation "<<eqn.que_suis_je()<<finl;
      const DoubleTab& inut=futur;
 
      store_equation_parameters(eqn,stored_parameters);
      convergence_eqn = le_solveur->iterer_eqn(eqn,inut, present, dt_,ite, ok);
      modify_equation_parameters(eqn,stored_parameters);
 
      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 line performs:
      // Update NS (hence update unknowns)
      // Update turbulence model (k-eps)
      //   eqn.inconnue().mettre_a_jour(temps);
      //   eqn.inconnue().reculer();
      changer_temps(eqn,temps);
    }
  return (convergence_pb==true);
}
 
void Schema_Implicite_Multi_TimeStep_base::test_stationnaire(Probleme_base& pb)
{
  int i = 0;
  int j = 0;
 
  for(i=0; i<pb.nombre_d_equations(); i++)
    {
      // treatement of time convergency
      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));
      update_time_derivative(pb.equation(i),present);
      for (j=0; j<nb_valeurs_passees(); j++)
        {
          pb.equation(i).inconnue().passe(j) = pb.equation(i).inconnue().passe(j+1);
        }
    }
}
 
static void Implicite_Multi_TimeStep_recommencer_pb(int nb_valeurs_passees, Probleme_base& pb)
{
  int i = 0;
  int j = 0;
  for(i=0; i<pb.nombre_d_equations(); i++)
    {
      for (j=0; j<nb_valeurs_passees; j++)
        {
          pb.equation(i).inconnue().passe(j) = pb.equation(i).inconnue().passe(j+1);
        }
    }
}
 
 
bool Schema_Implicite_Multi_TimeStep_base::iterateTimeStep(bool& converged)
{
  Probleme_base& prob=pb_base();
 
  int convergence_pb = 0;
  Initialiser_Champs(prob);
 
  int ok=0;
  int compteur;
  while (!ok)
    {
      compteur=0;
      Cout<<" "<<finl;
      Cout<<que_suis_je()<<" : solving of problem "<<prob.que_suis_je()<<finl;
 
      while ((!convergence_pb)&&(compteur<nb_ite_max))
        {
          compteur++;
          prob.updateGivenFields();
          convergence_pb = Iterer_Pb(prob,compteur);
          Cout<<" "<<finl;;
        }
 
      if ((!convergence_pb)&&(compteur==nb_ite_max))
        {
          Cerr<<que_suis_je()<<" has not converged at t="<< temps_courant_ << " with dt =" << dt_<< finl;
 
          double rap =  sqrt(2.);
          facsec_    /= rap;
          dt_        /= rap;
          ok         =0;
 
          Cerr<<"We restart with dt_ =" << dt_ <<finl;
 
          prob.abortTimeStep();
          prob.initTimeStep(dt_);
 
          // PL fix for crash in U_in_var_impl
          const int nb_eqn=prob.nombre_d_equations();
          for(int ii=0; ii<nb_eqn; ii++)
            {
              Equation_base& eqn=prob.equation(ii);
              mettre_a_jour_equation(eqn,temps_courant_+dt_);
            }
          // Fin
 
          Implicite_Multi_TimeStep_recommencer_pb(nb_valeurs_passees(),prob);
        }
      else
        {
          Cout<<"The problem "<<prob.que_suis_je()<<" has converged after "<<compteur<<" iterations."<<finl;
          ok=1;
        }
 
    }
  // 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);
 
  converged=true;
 
  return true;
}
 
int Schema_Implicite_Multi_TimeStep_base::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 convergence_pb;
  int i;
 
  for(i=0; i<pbc.nb_problemes(); i++)
    Initialiser_Champs(ref_cast(Probleme_base,pbc.probleme(i)));
 
 
  int compteur;
 
  while (!ok)
    {
      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());
                }
 
            }
 
          for(i=0; i<pbc.nb_problemes(); i++)
            {
              pbc.probleme(i).updateGivenFields();
              convergence_pb = Iterer_Pb(ref_cast(Probleme_base,pbc.probleme(i)),compteur);
              convergence_pbc = convergence_pbc * convergence_pb ;
            }
 
        }
      if ((!convergence_pbc)&&(compteur==nb_ite_max))
        {
          Cerr<<"The process is restarted. " <<finl;
          facsec_/=2;
          dt_/=2;
          ok=0;
          for(int i2=0; i2<pbc.nb_problemes(); i2++)
            {
              ref_cast(Probleme_base,pbc.probleme(i2)).abortTimeStep();
              ref_cast(Probleme_base,pbc.probleme(i2)).initTimeStep(dt_);
              Implicite_Multi_TimeStep_recommencer_pb(nb_valeurs_passees(),ref_cast(Probleme_base,pbc.probleme(i2)));
            }
        }
      else ok=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;
}
 
void Schema_Implicite_Multi_TimeStep_base::ajouter_inertie(Matrice_Base& mat_morse,DoubleTab& secmem,const Equation_base& eqn) const
{
  if (nb_pas_dt()>nb_pas_dt_seuil())
    {
      const Equation_base&  eqn_bis   = ref_cast(Equation_base,eqn);
      int             i         = 0;
      int             offset    = 0;
      double          time_step = pas_de_temps();
 
      //Building of a time-in-ascending-order table
      DoubleTab times(nb_valeurs_temporelles_effectives());
      assert(nb_valeurs_passees()+2 == nb_valeurs_temporelles_effectives());
 
      for (i=0; i<nb_valeurs_passees(); ++i)
        {
          offset   =  nb_valeurs_passees()-i;
          times[i] =  eqn_bis.inconnue().recuperer_temps_passe(offset); //past time
        }
      times[nb_valeurs_passees()]   = eqn_bis.inconnue().recuperer_temps_futur(0); //present time
      times[nb_valeurs_passees()+1] = eqn_bis.inconnue().recuperer_temps_futur(1); //future time
 
      compute_coefficients(time_step,times);
      add_multi_timestep_data(eqn, mat_morse, secmem);
    }
  else
    {
      //Backward Euler time scheme
      Schema_Temps_base::ajouter_inertie(mat_morse,secmem,eqn);
    }
}
 
int Schema_Implicite_Multi_TimeStep_base::faire_un_pas_de_temps_eqn_base(Equation_base& eqn)
{
  //Necessary because of equations like Transport_K_Eps which must be time incremented but
  //are not equations of the solved problem and therefore are not seen by Schema_Implicite_Multi_TimeStep_base::iterateTimeStep
 
  DoubleTab& passe           = eqn.inconnue().passe();
  DoubleTab& present         = eqn.inconnue().valeurs();
  DoubleTab& futur           = eqn.inconnue().futur();
  int        compteur        = 0, ok = 1;
  int        i               = 0;
  bool       convergence_eqn = false;
 
  DoubleTab      stored_parameters;
 
  authorized_equation(eqn);
 
  //Initialize time data
  if ( (eqn.que_suis_je() != "Euler") || temps_courant_ ==0.)
    {
      for (i=nb_valeurs_passees(); i>-1; i--)
        {
          eqn.inconnue().passe(i+1) = eqn.inconnue().passe(i);
        }
    }
 
  // used for the pressure and inter-problem couplings
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
 
  Cout << finl;
  while ((!convergence_eqn)&&(compteur<nb_ite_max))
    {
      compteur++;
 
      Cout<<"==================================================================================" << finl;
      Cout<<que_suis_je()<<": Implicit iteration " << compteur << " on the "<<eqn.que_suis_je() << " equation of the problem "<< eqn.probleme().le_nom()<< " :" <<finl;
      Cout<<"==================================================================================" << finl;
 
      const DoubleTab& inut=futur;
      store_equation_parameters(eqn,stored_parameters);
      convergence_eqn=le_solveur->iterer_eqn(eqn, inut, present, dt_, compteur, ok);
      modify_equation_parameters(eqn,stored_parameters);
      futur=present;
      eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_courant()+pas_de_temps());
      present=futur;
    }
 
  //Calculation of time derivative at time n+1
  present -= passe; //"present" contains unknown at time n+1 and "passe" contains unknown at time n
  present/=dt_;
 
  update_critere_statio(present, eqn);
  update_time_derivative(eqn,present);
 
  //Update time data
  for (i=0; i<nb_valeurs_passees(); i++)
    {
      eqn.inconnue().passe(i) = eqn.inconnue().passe(i+1);
    }
 
  return 1;
}
 
int Schema_Implicite_Multi_TimeStep_base::reprendre(Entree& s)
{
  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 file opening
  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;
      while (!fichier.eof() && temps<tinit_)
        {
          facsec_lu_old = facsec_lu;
          fichier >> dt;
          fichier >> facsec_lu;
          fichier >> 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_lu;
              nb_ite_sans_accel_ = 1;
            }
          nb_ite_sans_accel_++;
        }
    }
  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;
}
 
//To get du/dt table without pressure gradient and influences of past and present times
void Schema_Implicite_Multi_TimeStep_base::modifier_second_membre_full_impl(const Equation_base& eqn, DoubleTab& secmem)
{
  Cerr<<"Error in "<<que_suis_je()<<"::modifier_second_membre_full_impl()"<<finl;
  Cerr<<que_suis_je()<<" is an implicit time scheme"<<finl;
  Cerr<<que_suis_je()<<"::modifier_second_membre_full_impl() should not be used"<<finl;
  Process::exit();
}
 
void Schema_Implicite_Multi_TimeStep_base::modifier_second_membre(const Equation_base& eqn, DoubleTab& secmem)
{
  Cerr<<"Error in "<<que_suis_je()<<"::modifier_second_membre()"<<finl;
  Cerr<<que_suis_je()<<" is an implicit time scheme"<<finl;
  Cerr<<que_suis_je()<<"::modifier_second_membre() should not be used"<<finl;
  Process::exit();
}
 
void Schema_Implicite_Multi_TimeStep_base::authorized_equation(const Equation_base& eqn)
{
  Nom  equation_name(eqn.que_suis_je());
  bool authorized =  true;
 
  authorized &= (bool) equation_name.find("FT");
  authorized &= (bool) equation_name.find("Front_Tracking");
  authorized &= (bool) equation_name.find("QC");
  authorized &= (bool) equation_name.find("Quasi_Compressible");
 
  if (!authorized)
    {
      Cerr<<"Error in "<<que_suis_je()<<"::authorized_equation()"<<finl;
      Cerr<<"Equation "<<equation_name<<" is not allowed with "<<que_suis_je()<<" time scheme"<<finl;
      Process::exit();
    }
}
 
void Schema_Implicite_Multi_TimeStep_base::authorized_equation(const Equation_base& eqn) const
{
  Nom  equation_name(eqn.que_suis_je());
  bool authorized =  true;
 
  authorized &= (bool) equation_name.find("FT");
  authorized &= (bool) equation_name.find("Front_Tracking");
  authorized &= (bool) equation_name.find("QC");
  authorized &= (bool) equation_name.find("Quasi_Compressible");
 
  if (!authorized)
    {
      Cerr<<"Error in "<<que_suis_je()<<"::authorized_equation()"<<finl;
      Cerr<<"Equation "<<equation_name<<" is not allowed with "<<que_suis_je()<<" time scheme"<<finl;
      Process::exit();
    }
}