Schema_Temps_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_Temps_base.h>
#include <EcrFicCollecte.h>
#include <communications.h>
#include <Probleme_base.h>
#include <Matrice_Morse.h> // necessary for visual
#include <LecFicDiffuse.h>
#include <Equation_base.h>
#include <TRUST_2_PDI.h>
#include <sys/stat.h>
#include <SFichier.h>
#include <Param.h>
#include <Debog.h>
#include <cfloat>
#include <DeviceMemory.h>
#include <Perf_counters.h>
 
// XD dt_start class_generic dt_start NO_BRACE not_set
// XD dt_calc_dt_calc dt_start dt_calc NO_BRACE The time step at first iteration is calculated in agreement with CFL
// XD_CONT condition.
// XD dt_calc_dt_min dt_start dt_min NO_BRACE The first iteration is based on dt_min.
 
// XD dt_calc_dt_fixe dt_start dt_fixe NO_BRACE The first time step is fixed by the user (recommended when resuming
// XD_CONT calculation with Crank Nicholson temporal scheme to ensure continuity).
// XD attr value floattant value REQ first time step.
 
 
Implemente_base_sans_constructeur(Schema_Temps_base,"Schema_Temps_base",Objet_U);
// XD schema_temps_base objet_u schema_temps_base INHERITS_BRACE Basic class for time schemes. This scheme will be
// XD_CONT associated with a problem and the equations of this problem.
/* Attributes further down in the cpp: */
 
/*! @brief Default constructor of a time scheme.
 *
 * Initializes different class members.
 *
 */
Schema_Temps_base::Schema_Temps_base() :
  norm_residu_("max")
{
}
 
void Schema_Temps_base::initialize()
{
  nb_impr_= 0;
  nb_pas_dt_ = 0;
 
  // GF re-calling calculer_pas_de_temps because some implicit time schemes use it for dimensioning (in particular ovap)
  pb_base().calculer_pas_de_temps();
  // calling it a second time for alternating schemes....
  pb_base().calculer_pas_de_temps();
 
  // Outputs:
  bool init = true;
  write_dt_ev(init);
  write_progress(init);
 
  if ( nb_pas_dt_ == 0 && ( ( mode_dt_start_ == 0. && est_egal(tinit_,0.)) || mode_dt_start_ == -1.) )
    {
      //We divide by facsec_ because multiplication by facsec_ is done in corriger_dt_calcule()
      //regardless of the dt_ initialization mode
      dt_=dt_min_/facsec_;
    }
  else if (mode_dt_start_ > 0.)
    {
      dt_ = mode_dt_start_ / facsec_;
    }
  else
    {
      mettre_a_jour_dt_stab();
      dt_=dt_stab_;
    }
  dt_stab_=dt_;
  dt_failed_ = DBL_MAX;
}
 
double Schema_Temps_base::computeTimeStep(bool& is_stop) const
{
  //re-evaluation of dt_max_ if it is a function of time
  if (dt_max_str_ != Nom())
    {
      dt_max_fn_.setVar(0, temps_courant());
      dt_max_ = dt_max_fn_.eval();
    }
  is_stop=false;
  // First correct the time step
  double dt = dt_stab_;
  if (!corriger_dt_calcule(dt)) // Change le contenu de dt
    is_stop=true;
  dt = std::min(dt, dt_failed_ / sqrt(2)); //to force a dt decrease in case of failure at the previous resolution
  if (temps_courant_ > temps_precedent_)
    dt = std::min(dt, (temps_courant_ - temps_precedent_) * dt_gf_); //to avoid increasing dt too fast (like facsec)
 
  if (limpr() || (nb_pas_dt_ == 0))
    Cout << "Time step finally used to solve the next time step (taking into account facsec) : " << dt << " s." << finl;
 
  // If dt has been reduced after a failed step, enforce dt_min check again.
  if ((dt - dt_min_) / (dt + DMINFLOAT) < -1.e-6 && !adapt_dt_tmax_)
    {
      Cerr << "---------------------------------------------------------" << finl;
      Cerr << "Problem with the time step " << dt << " which is less than dt_min " << dt_min_ << finl;
      Cerr << "Lower dt_min value or check why the time step decreases..." << finl;
      Cerr << "Results are saved to provide help." << finl;
      Cerr << "---------------------------------------------------------" << finl;
      Probleme_base& pb = ref_cast_non_const(Probleme_base, mon_probleme.valeur());
      pb.postraiter(1);
      pb.sauver();
      Process::exit();
    }
 
  // Immediate update of dt_ attribute so that pas_de_temps()
  // is always up to date (especially at post-processing time)
  // This was not the case for versions <= 1.6.3 and fields depending
  // on dt were not correct if facsec<>1
  //  dt_ = dt;
 
  if (this->stop())
    {
      is_stop=true;
      //return 0; Why ?
      return dt;
    }
  else
    return dt;
}
 
bool Schema_Temps_base::initTimeStep(double dt)
{
 
  dt_=dt;
  stationnaire_atteint_ = -1;
  residu_=0;
 
  return true;
}
 
/*! @brief Calculate the U(n+1) unknown for each equation (if solved) of the problem with the selected time scheme
 *
 */
bool Schema_Temps_base::iterateTimeStep(bool& converged)
{
  // Loop on the equations of the problem:
  for(int i=0; i<pb_base().nombre_d_equations(); i++)
    {
      Equation_base& equation=pb_base().equation(i);
      // Calculate (if equation is solved) the unknown equation U(n+1)
      // according to the selected time scheme:
      if (!equation.equation_non_resolue())
        faire_un_pas_de_temps_eqn_base(equation);
      else
        {
          Cout<< "====================================================" << finl;
          Cout<< equation.que_suis_je()<<" equation is not solved."<<finl;
          Cout<< "====================================================" << finl;
          // Compute the derivative once to obtain the boundary fluxes
          if (equation.schema_temps().nb_pas_dt()==0)
            {
              DoubleTab inconnue_valeurs(equation.inconnue().valeurs());
              equation.derivee_en_temps_inco(inconnue_valeurs);
            }
        }
    }
  converged=true;
  return true;
}
 
/*! @brief Returns 1 if there is a need to perform a print (cf dt_impr) Returns 0 otherwise
 *
 * @return (int) 1 if there is a need to perform a print 0 otherwise
 */
int Schema_Temps_base::limpr() const
{
  if (dt_impr_<=0)
    return 0; // Negative value for dt_impr : we never print
  if (nb_pas_dt_==0)
    return 0; // No solve cause 0 time step : we print nothing
  if (dt_impr_<=dt_)
    return 1; // We print every time step if dt_impr is lower than dt
  if (nb_pas_dt_<=1)
    return 1; // We print the first time step
  if (tmax_<=temps_courant_ || nb_pas_dt_max_<=nb_pas_dt_ || stationnaires_atteints_)
    return 1; // We print the last time step
 
  // 26/01/2010: We now use the function "modf(operation, &integer_part)" from math.h
  // This function decomposes the result of "operation" into an integer part and a decimal part.
  // This avoids converting a double to int with all the associated risks.
  // ex: modf(9/7,&i) gives 1.000000 for the integer part and 0.285714 for the decimal part
  //
  // epsilon ensures that the result of the operation is independent of machine precision:
  // ex: 0.99999999/1.00000000 can give 0.99999999 or 1.00000000 depending on the machine
  // whereas
  // ex: 0.99999999/1.00000000 + 1.e-8 always gives 1.00000000.
  double i, j, epsilon = 1.e-8;
  modf(temps_courant_/dt_impr_ + epsilon, &i);
  modf(temps_precedent_/dt_impr_ + epsilon, &j);
  return ( i>j );
}
 
void Schema_Temps_base::validateTimeStep()
{
  statistics().begin_count(STD_COUNTERS::update_variables,statistics().get_last_opened_counter_level()+1);
  // Update the problem:
  Probleme_base& problem=pb_base();
  problem.mettre_a_jour(temps_courant_+dt_);
 
  // Outputs
  bool init=false;
  write_dt_ev(init);
  write_progress(init);
 
  // Update time scheme:
  mettre_a_jour();
  statistics().end_count(STD_COUNTERS::update_variables);
  dt_failed_ = DBL_MAX;
}
void Schema_Temps_base::terminate()
{
  if (je_suis_maitre() && !disable_progress())
    progress_<< (int)100<< finl;
}
 
/*! @brief Returns 1 if during the last time step, the problem has not evolved.
 *
 * @param (temps) the time to reach
 * @return (1 if the problem has not evolved, 0 otherwise.)
 */
bool Schema_Temps_base::isStationary() const
{
  if (stationnaire_atteint_ && (nb_pas_dt_ >= 2))
    {
      Cerr << "---------------------------------------------------------"
           << finl
           << "The problem " << pb_base().le_nom()
           << " has reached the steady state"
           << finl << finl;
      return true;
    }
  return false;
}
 
void Schema_Temps_base::abortTimeStep()
{
  Probleme_base& prob=pb_base();
  for(int i=0; i<prob.nombre_d_equations(); i++)
    prob.equation(i).abortTimeStep();
 
  dt_ = prob.calculer_pas_de_temps();
}
 
void Schema_Temps_base::resetTime(double time)
{
  //
  // It happens here: -)
  //
  temps_courant_ = time;
}
 
void Schema_Temps_base::associer_pb(const Probleme_base& un_probleme)
{
  mon_probleme=un_probleme;
}
 
 
 
void Schema_Temps_base::set_param(Param& param) const
{
  param.ajouter("tinit",&tinit_); // XD_ADD_P double
  // XD_CONT Value of initial calculation time (0 by default).
  param.ajouter( "tmax",&tmax_); // XD_ADD_P double
  // XD_CONT Time during which the calculation will be stopped (1e30s by default).
  param.ajouter_non_std( "tcpumax",(this)); // XD_ADD_P double
  // XD_CONT CPU time limit (must be specified in hours) for which the calculation is stopped (1e30s by default).
  param.ajouter( "dt_min",&dt_min_); // XD_ADD_P double
  // XD_CONT Minimum calculation time step (1e-16s by default).
  param.ajouter( "dt_max",&dt_max_str_); // XD_ADD_P chaine
  // XD_CONT Maximum calculation time step as function of time (1e30s by default).
  param.ajouter( "dt_sauv",&dt_sauv_); // XD_ADD_P double
  // XD_CONT Save time step value (1e30s by default). Every dt_sauv, fields are saved in the .sauv file. The file
  // XD_CONT contains all the information saved over time. If this instruction is not entered, results are saved only
  // XD_CONT upon calculation completion. To disable the writing of the .sauv files, you must specify 0. Note that
  // XD_CONT dt_sauv is in terms of physical time (not cpu time).
  param.ajouter( "nb_sauv_max",&nb_sauv_max_); // XD_ADD_P entier
  // XD_CONT Maximum number of timesteps that will be stored in backup file (10 by default). This value is only useful
  // XD_CONT when doing a complete backup of the calculation with parallel PDI (as it needs to allocate the proper
  // XD_CONT amount of dataspace in advance). If this number is reached (ie we already stored the data of nb_sauv_max
  // XD_CONT timesteps in the file), the next checkpoints will overwrite the first ones
  param.ajouter( "dt_impr",&dt_impr_); // XD_ADD_P double
  // XD_CONT Scheme parameter printing time step in time (1e30s by default). The time steps and the flux balances are
  // XD_CONT printed (incorporated onto every side of processed domains) into the .out file.
  param.ajouter_non_std("facsec",(this)); // XD_ADD_P chaine
  // XD_CONT Value assigned to the safety factor for the time step (1. by default). It can also be a function of time.
  // XD_CONT The time step calculated is multiplied by the safety factor. The first thing to try when a calculation does
  // XD_CONT not converge with an explicit time scheme is to reduce the facsec to 0.5. NL2 Warning: Some schemes needs a
  // XD_CONT facsec lower than 1 (0.5 is a good start), for example Schema_Adams_Bashforth_order_3.
  param.ajouter( "seuil_statio",&seuil_statio_); // XD_ADD_P double
  // XD_CONT Value of the convergence threshold (1e-12 by default). Problems using this type of time scheme converge
  // XD_CONT when the derivatives dGi/dt NL1 of all the unknown transported values Gi have a combined absolute value
  // XD_CONT less than this value. This is the keyword used to set the permanent rating threshold.
  param.ajouter_non_std("residuals", (this));    // XD_ADD_P residuals
  // XD_CONT To specify how the residuals will be computed (default max norm, possible to choose L2-norm instead).
  param.ajouter( "diffusion_implicite",&ind_diff_impl_); // XD_ADD_P entier
  // XD_CONT Keyword to make the diffusive term in the Navier-Stokes equations implicit (in this case, it should be set
  // XD_CONT to 1). The stability time step is then only based on the convection time step (dt=facsec*dt_convection).
  // XD_CONT Thus, in some circumstances, an important gain is achieved with respect to the time step (large diffusion
  // XD_CONT with respect to convection on tightened meshes). Caution: It is however recommended that the user avoids
  // XD_CONT exceeding the convection time step by selecting a too large facsec value. Start with a facsec value of 1
  // XD_CONT and then increase it gradually if you wish to accelerate calculation. In addition, for a natural convection
  // XD_CONT calculation with a zero initial velocity, in the first time step, the convection time is infinite and
  // XD_CONT therefore dt=facsec*dt_max.
  param.ajouter( "seuil_diffusion_implicite",&seuil_diff_impl_); // XD_ADD_P double
  // XD_CONT This keyword changes the default value (1e-6) of convergency criteria for the resolution by conjugate
  // XD_CONT gradient used for implicit diffusion.
  param.ajouter( "impr_diffusion_implicite",&impr_diff_impl_); // XD_ADD_P entier
  // XD_CONT Unactivate (default) or not the printing of the convergence during the resolution of the conjugate
  // XD_CONT gradient.
  param.ajouter( "impr_extremums",&impr_extremums_); // XD_ADD_P entier
  // XD_CONT Print unknowns extremas
  param.ajouter( "no_error_if_not_converged_diffusion_implicite",&no_error_if_not_converged_diff_impl_); // XD_ADD_P entier
  // XD_CONT not_set
  param.ajouter( "no_conv_subiteration_diffusion_implicite",&no_conv_subiteration_diff_impl_); // XD_ADD_P entier
  // XD_CONT not_set
  param.ajouter_non_std( "dt_start",(this)); // XD_ADD_P dt_start
  // XD_CONT dt_start dt_min : the first iteration is based on dt_min. NL2 dt_start dt_calc : the time step at first
  // XD_CONT iteration is calculated in agreement with CFL condition. NL2 dt_start dt_fixe value : the first time step
  // XD_CONT is fixed by the user (recommended when resuming calculation with Crank Nicholson temporal scheme to ensure
  // XD_CONT continuity). NL2 By default, the first iteration is based on dt_calc.
  param.ajouter_non_std( "nb_pas_dt_max",(this)); // XD_ADD_P entier
  // XD_CONT Maximum number of calculation time steps (1e9 by default).
  param.ajouter( "niter_max_diffusion_implicite",&niter_max_diff_impl_); // XD_ADD_P entier
  // XD_CONT This keyword changes the default value (number of unknowns) of the maximal iterations number in the
  // XD_CONT conjugate gradient method used for implicit diffusion.
  param.ajouter( "precision_impr",&precision_impr_); // XD_ADD_P entier
  // XD_CONT Optional keyword to define the digit number for flux values printed into .out files (by default 3).
  param.ajouter_non_std( "periode_sauvegarde_securite_en_heures",(this)); // XD_ADD_P double
  // XD_CONT To change the default period (23 hours) between the save of the fields in .sauv file.
  param.ajouter_non_std( "no_check_disk_space",(this)); // XD_ADD_P flag
  // XD_CONT To disable the check of the available amount of disk space during the calculation.
  param.ajouter_flag( "disable_progress",&disable_progress_); // XD_ADD_P flag
  // XD_CONT To disable the writing of the .progress file.
  param.ajouter_flag( "disable_dt_ev",&disable_dt_ev_); // XD_ADD_P flag
  // XD_CONT To disable the writing of the .dt_ev file.
  param.ajouter_flag("adapt_dt_tmax", &adapt_dt_tmax_); // XD_ADD_P flag
  // XD_CONT Use to adapt final dt when approaching tmax.
  param.ajouter( "gnuplot_header",&gnuplot_header_); // XD_ADD_P entier
  // XD_CONT Optional keyword to modify the header of the .out files. Allows to use the column title instead of columns
  // XD_CONT number.
 
  // XD residuals interprete nul BRACE To specify how the residuals will be computed.
  // XD attr norm chaine(into=["L2","max"]) norm OPT allows to choose the norm we want to use (max norm by default).
  // XD_CONT Possible to specify L2-norm.
// XD attr relative chaine(into=["0","1","2"]) relative OPT This is the old keyword seuil_statio_relatif_deconseille. If
// XD_CONT it is set to 1, it will normalize the residuals with the residuals of the first 5 timesteps (default is 0).
// XD_CONT if set to 2, residual will be computed as R/(max-min).
}
 
/*! @brief Overrides Objet_U::printOn(Sortie&): prints the time scheme to an output stream.
 *
 *     !! Note: this is not symmetric with reading !!
 *     Writes the various parameters of the time scheme.
 *
 * @param (Sortie& os) the output stream
 * @return (Sortie&) the modified output stream
 */
Sortie& Schema_Temps_base::printOn(Sortie& os) const
{
  os << "dt " << dt_ << finl;
  os << "temps_courant " << temps_courant_ << finl ;
  os << "tinit " << tinit_ << finl;
  os << "tmax " << tmax_ << finl ;
  os << "tcpumax " << tcpumax_ << finl ;
  os << "nb_pas_dt " << nb_pas_dt_ << finl;
  os << "nb_pas_dt_max " <<  nb_pas_dt_max_ << finl;
  os << "dt_min " << dt_min_ << finl;
  os << "dt_max " << dt_max_ << finl;
  os << "facsec " << facsec_ << finl;
  os << "seuil_statio" << seuil_statio_ << finl;
  os << "seuil_statio_relatif_deconseille" << seuil_statio_relatif_deconseille_ << finl;
  os << "norm_residu" << norm_residu_ << finl;
  os << "dt_sauv " << dt_sauv_ << finl;
  os << "nb_sauv_max " << nb_sauv_max_ << finl;
  os << "limite_cpu_sans_sauvegarde " << limite_cpu_sans_sauvegarde_ << finl;
  os << "dt_impr " << dt_impr_ << finl;
  os << "precision_impr " << precision_impr_ << finl;
  os << "stationnaire_atteint " << stationnaire_atteint_ << finl;
  os << "diffusion_implicite " << ind_diff_impl_ << finl ;
  os << "seuil_diffusion_implicite " << seuil_diff_impl_ << finl ;
  os << "impr_diffusion_implicite " << impr_diff_impl_ << finl ;
  os << "impr_extremums " << impr_extremums_ << finl ;
  os << "niter_max_diffusion_implicite " << niter_max_diff_impl_ << finl ;
  os << "no_file_allocation " << file_allocation_ << finl ;
  os << "disable_progress " << int(disable_progress_) << finl ; // TODO allow outputs of bool in Sortie
  os << "disable_dt_ev " << int(disable_dt_ev_) << finl ;
  os << "adapt_dt_tmax " << int(adapt_dt_tmax_) << finl ;
  os << "fin " << finl;
  return os ;
}
 
 
/*! @brief Reads a time scheme from an input stream.
 *
 * The expected read format is:
 *      {
 *       [Motcle valeur_reelle]
 *      }
 *      The keywords can be:
 *      tinit, tmax, nb_pas_dt_max, dt_min, dt_max,
 *      dt_sauv, dt_impr, facsec, seuil_statio,
 *
 * @param (Entree& is) the input stream
 * @return (Entree&) the modified input stream
 * @throws opening brace expected
 * @throws unknown keyword at this location
 */
Entree& Schema_Temps_base::readOn(Entree& is)
{
  Cerr<<"Reading of data for a "<<que_suis_je()<<" time scheme"<<finl;
  Param param(que_suis_je());
  set_param(param);
  param.lire_avec_accolades_depuis(is);
  if (dt_min_==0)
    {
      Cerr<<"dt_min has not been read or is set to 0."<<finl;
      Cerr<<"Assign a value strictly positive to dt_min."<<finl;
      exit();
    }
  temps_courant_=tinit_;
  temps_precedent_=tinit_;
  lu_=true;
  if (dt_sauv_ <= 0.)
    Cerr << "NO next backup, by security, because dt_sauv = " << dt_sauv_ << finl;
  else
    Cerr << "The next backup, by security, will take place after " << limite_cpu_sans_sauvegarde_/3600 << " hours of calculation." << finl;
 
  if (dt_max_str_ != Nom())
    {
      dt_max_fn_.setNbVar(1);
      dt_max_fn_.setString(dt_max_str_);
      dt_max_fn_.addVar("t");
      dt_max_fn_.parseString();
      dt_max_ = dt_max_fn_.eval();
    }
  return is ;
}
 
int Schema_Temps_base::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  Motcle motlu;
  int retval = 1;
  if (mot=="dt_start")
    {
      is>>motlu;
      Motcles les_mots2(4);
      les_mots2[0]="dt_std";
      les_mots2[1]="dt_min";
      les_mots2[2]="dt_calc";
      les_mots2[3]="dt_fixe";
      int rang2=les_mots2.search(motlu);
      switch(rang2)
        {
        case 0:
        case 1:
        case 2:
          {
            mode_dt_start_ = (double)(-rang2);
            break;
          }
        case 3:
          {
            is >> mode_dt_start_;
            break;
          }
        default:
          {
            Cerr<<" We do not understand "<<motlu <<"in Schema_Temps_base::lire_motcle_non_standard"<<finl;
            Cerr<<"keywords understood "<<les_mots2<<finl;
            exit();
          }
        }
      return 1;
    }
  else if (mot=="nb_pas_dt_max")
    lire_nb_pas_dt_max(is);
  else if (mot=="periode_sauvegarde_securite_en_heures")
    lire_periode_sauvegarde_securite_en_heures(is);
  else if (mot=="tcpumax")
    lire_temps_cpu_max(is);
  else if (mot=="no_check_disk_space")
    file_allocation_=0;
  else if (mot == "residuals")
    lire_residuals(is);
  else if(mot == "facsec")
    is >> facsec_;
  else
    retval = -1;
 
  return retval;
}
 
/*! @brief Reads the maximum number of time steps.
 *
 */
Entree& Schema_Temps_base::lire_nb_pas_dt_max(Entree& is)
{
  Cerr << "Reading of the maximum number of time steps"  << finl;
  is >> nb_pas_dt_max_;
  return is;
}
 
Entree& Schema_Temps_base::lire_periode_sauvegarde_securite_en_heures(Entree& is)
{
  Cerr << "Reading of the safety backup period in hours"  << finl;
  is >> periode_cpu_sans_sauvegarde_;
  periode_cpu_sans_sauvegarde_*=3600; // Conversion to seconds
  limite_cpu_sans_sauvegarde_ = periode_cpu_sans_sauvegarde_;
  return is;
}
 
Entree& Schema_Temps_base::lire_temps_cpu_max(Entree& is)
{
  Cerr << "Reading the max cpu time allowed"  << finl;
  is >> tcpumax_;
  tcpumax_*=3600; // Conversion to seconds
  return is;
}
 
 
Entree& Schema_Temps_base::lire_residuals(Entree& is)
{
  Motcle m;
  is >> m;
  assert (m == "{");
  is >> m;
  while (m != "}")
    {
      Motcles residuals_mots(2);
      residuals_mots[0]="relative";
      residuals_mots[1]="norm";
      int res_rang=residuals_mots.search(m);
      switch(res_rang)
        {
        case 0:
          is >> seuil_statio_relatif_deconseille_;
          break;
        case 1:
          is >> norm_residu_;
          break;
        default :
          {
            Cerr<<" We do not understand "<<m <<"in Schema_Temps_base::lire_motcle_non_standard"<<finl;
            Cerr<<" keywords understood "<<residuals_mots<<finl;
            exit();
          }
        }
      is >> m;
    }
  return is;
}
 
 
/*! @brief Prints the time step number, the time step value, and the current time.
 *
 * @param (Sortie& os) the output stream
 * @return (int) always returns 0
 */
int Schema_Temps_base::impr(Sortie& os) const
{
  os << finl;
  os << "-------------------------------------------------------------------" << finl;
  os << "We finished treating the time step number "<< nb_pas_dt_ << " , for the time scheme ..." << finl
     << "   stable dt used = " << dt_ << finl /* dt_stab_ perhaps? it is the same I think */
     << "   time achieved (in seconds) = " << temps_courant_ << finl;
  return 0;
}
 
/*! @brief Prints the time step number, the time step value, and the current time.
 *
 * @param (Sortie& os) the output stream
 * @return (int) always returns 0
 */
int Schema_Temps_base::impr(Sortie& os,Probleme_base& pb) const
{
  os << finl;
  os << " Problem : " << pb.le_nom() << finl;
  os << "We treat the time step number "<< nb_pas_dt_ << finl;
  os << "   dt = " << dt_ << finl;
  os << "   time = " << temps_courant_ << finl;
  os << "   time scheme : " << pb.schema_temps() << finl;
  return 0;
}
int Schema_Temps_base::impr(Sortie& os,const Probleme_base& pb) const
{
  os << finl;
  os << " Problem : " << pb.le_nom() << finl;
  os << "We treat the time step number "<< nb_pas_dt_ << finl;
  os << "   dt = " << dt_ << finl;
  os << "   time = " << temps_courant_ << finl;
  os << "   time scheme : " << pb.schema_temps() << finl;
  return 0;
}
 
extern "C" {
  int ccc_tremain(double*);
}
/*! @brief Updates the current time (t+=dt) and the number of time steps performed (nb_pas_dt_++).
 *
 * @return (int) always returns 1
 */
int Schema_Temps_base::mettre_a_jour()
{
  temps_precedent_ = temps_courant_;
  temps_courant_ += dt_;
  nb_pas_dt_++;
  // Compute next time step stability:
  statistics().end_count(STD_COUNTERS::update_variables,0,0);
  mettre_a_jour_dt_stab();
  statistics().begin_count(STD_COUNTERS::update_variables,statistics().get_last_opened_counter_level()+1);
  assert_parallel(dt_stab_);
  assert_parallel(temps_courant_);
  if (!ind_tps_final_atteint)
    ind_tps_final_atteint = (temps_courant_ >= tmax_);
 
  if (!ind_nb_pas_dt_max_atteint)
    ind_nb_pas_dt_max_atteint = (nb_pas_dt_ >= nb_pas_dt_max_);
 
  if (!ind_temps_cpu_max_atteint)
    ind_temps_cpu_max_atteint = (temps_cpu_ecoule_ >= tcpumax_);
 
  // Increment limite_cpu_sans_sauvegarde_
  if ( temps_cpu_ecoule_ > limite_cpu_sans_sauvegarde_ )
    {
      limite_cpu_sans_sauvegarde_ += periode_cpu_sans_sauvegarde_;
      //Finally, we double the limit, so by default security backup at 10h, 20h, 40h, 80h, 160h, ...
      //limite_cpu_sans_sauvegarde_ = 2 * limite_cpu_sans_sauvegarde_;
      if (dt_sauv_ <= 0.)
        Cerr << "NO next backup, by security, because dt_sauv = " << dt_sauv_ << finl;
      else
        Cerr << "The next backup, by security, will take place after " << limite_cpu_sans_sauvegarde_/3600 << " hours of calculation." << finl;
    }
// GF to ensure that all processes have the same elapsed time
  if (je_suis_maitre())
    {
      temps_cpu_ecoule_ = statistics().get_time_since_last_open(STD_COUNTERS::total_execution_time);
    }
 
  envoyer_broadcast(temps_cpu_ecoule_,0);
 
 
#ifdef LIBCCC_USER
  if (je_suis_maitre() && limpr())
    {
      Cout << "[CCRT] ";
      double second_remain;
      int error = ccc_tremain(&second_remain);
      if(!error)
        {
          int hour_remain = (int)(second_remain/3600);
          second_remain-=hour_remain*3600;
          int minute_remain = (int)(second_remain/60);
          second_remain-=minute_remain*60;
          Cout << hour_remain <<"h"<<minute_remain<<"mn"<<second_remain<<"s before job is killed on CCRT." << finl;
        }
      else
        Cout << "Error." << finl;
    }
#endif
  return 1;
}
// Makes a protection backup on 24h runs on CCRT machines
int Schema_Temps_base::lsauv() const
{
  if (dt_sauv_ <= 0.)
    return 0;
  else
    {
      if ( temps_cpu_ecoule_ > limite_cpu_sans_sauvegarde_ )
        {
          Cerr << "After these " << limite_cpu_sans_sauvegarde_/3600 << " hours of calculation, a backup will be made by security." << finl;
          return 1;
        }
      else
        {
          if (dt_sauv_<=dt_)
            return 1;
          else if (tmax_<=temps_courant_ || nb_pas_dt_max_<=nb_pas_dt_ || stationnaires_atteints_ || ind_temps_cpu_max_atteint || stop_lu())
            return 1;
          else
            {
              // See Schema_Temps_base::limpr for information on epsilon and modf
              double i, j, epsilon = 1.e-8;
              modf(temps_courant_/dt_sauv_ + epsilon, &i);
              modf(temps_precedent_/dt_sauv_ + epsilon, &j);
              return ( i>j ) ;
            }
        }
    }
}
 
void Schema_Temps_base::mettre_a_jour_dt_stab()
{
  DeviceMemory::nb_pas_dt_ = nb_pas_dt_;
  imprimer(Cout);
  dt_stab_=pb_base().calculer_pas_de_temps();
}
/*! @brief Prints the time step to an output stream if appropriate.
 *
 * @param (Sortie& os) the output stream
 */
void Schema_Temps_base::imprimer(Sortie& os) const
{
  if (limpr())
    {
      nb_impr_++;
      if (je_suis_maitre() && schema_impr())
        {
          impr(os);
        }
    }
}
 
/*! @brief Prints the time step to an output stream if appropriate.
 *
 * @param (Sortie& os) the output stream
 */
void Schema_Temps_base::imprimer(Sortie& os, Probleme_base& pb) const
{
  if (limpr() && je_suis_maitre())
    impr(os,pb);
}
void Schema_Temps_base::imprimer(Sortie& os, const Probleme_base& pb) const
{
  if (limpr() && je_suis_maitre())
    impr(os,pb);
}
 
/*! @brief Saves the current time and the number of time steps to an output stream.
 *
 * @param (Sortie& os) the output stream for the backup
 */
int Schema_Temps_base::sauvegarder(Sortie& os) const
{
  int bytes = 0;
  if(TRUST_2_PDI::is_PDI_checkpoint())
    {
      int simple_checkpoint = ref_cast_non_const(Probleme_base, mon_probleme.valeur()).is_sauvegarde_simple();
      int i = 0;
      int i_max = 1;
      if(!simple_checkpoint)
        {
          if(nb_sauv_ > nb_sauv_max_)
            Cerr << "WARNING ! Overwriting the first " << nb_sauv_max_ << " backups..." << finl;
          i =  nb_sauv_ % nb_sauv_max_;
          i_max = nb_sauv_max_;
        }
      TRUST_2_PDI pdi_interface;
      pdi_interface.TRUST_start_sharing("iter", &i);
      pdi_interface.TRUST_start_sharing("nb_iter_max", &i_max);
      pdi_interface.TRUST_start_sharing("temps", &temps_courant_);
      if(Process::node_master())
        pdi_interface.trigger("time_scheme");
      pdi_interface.stop_sharing_last_variable(); //temps
      pdi_interface.stop_sharing_last_variable(); //nb_iter_max
      pdi_interface.stop_sharing_last_variable(); //iter
 
      bytes = 8+4; // one double and one int
    }
  nb_sauv_++;
  return bytes;
}
 
/*! @brief Restarts (reads) the current time and the number of time steps performed from an input stream.
 *
 * @param (Entree& is) the input stream
 * @return (int) always returns 1
 */
int Schema_Temps_base::reprendre(Entree& is)
{
  return 1;
}
 
/*! @brief Returns 1 if the .stop file contains a 1, returns 0 otherwise.
 *
 * @return (int) 1 if the .stop file contains 1, 0 otherwise
 */
int Schema_Temps_base::stop_lu() const
{
  int stop_lu_l = 0;
  if (!get_disable_stop())
    {
      Nom nomfic(nom_du_cas());
      nomfic += ".stop";
      if (nb_pas_dt_ < 1)
        {
          if (je_suis_maitre())
            {
              SFichier ficstop(nomfic);
              ficstop << stop_lu_l;
            }
        }
      else
        {
          LecFicDiffuse ficstop(nomfic);
          ficstop >> stop_lu_l;
        }
    }
  return stop_lu_l;
}
 
/*! @brief Corrects the computed time step passed as parameter and verifies it is not "too" small (< dt_min_).
 *
 *     The correction is:
 *        delta_t = min((security factor * dt_calc), dt_max)
 *     And verifies that delta_t is "sufficiently" larger than dt_min_.
 *
 * @param (double& dt_calc) the computed time step to verify
 * @throws if the computed time step is less than dt_min
 */
bool Schema_Temps_base::corriger_dt_calcule(double& dt_calc) const
{
  // Print the RAM
  if (limpr())
    {
      Cout << finl;
      imprimer_ram_totale();
    }
 
  // proposed time step = stability time step (dt_stab) * security factor (facsec)
  const double dt_propose = dt_stab_ * facsec_;
 
  // Compute the time step dt as the minimal value between dt_max_ and proposed time step
  double dt = std::min(dt_max_, dt_propose);
 
  // if adapt_dt_tmax option is active ... cathare-style
  bool adapt_dt_tmax = false;
  if (adapt_dt_tmax_)
    {
      if (temps_courant_ + dt > tmax_)
        {
          dt = tmax_ - temps_courant_;
          adapt_dt_tmax = true;
        }
      else if (temps_courant_ + 2. * dt > tmax_)
        {
          dt = 0.5 * (tmax_ - temps_courant_);
          adapt_dt_tmax = true;
        }
    }
 
  dt_calc = dt;
 
  if ((dt - dt_min_) / (dt + DMINFLOAT) < -1.e-6 && !adapt_dt_tmax)
    {
      // Calculation stops if time step dt is less than dt_min
      Cerr << "---------------------------------------------------------" << finl;
      Cerr << "Problem with the time step " << dt << " which is less than dt_min " << dt_min_ << finl;
      Cerr << "Lower dt_min value or check why the time step decreases..." << finl;
      Cerr << "Results are saved to provide help." << finl;
      Cerr << "---------------------------------------------------------" << finl;
      Probleme_base& pb = ref_cast_non_const(Probleme_base, mon_probleme.valeur());
      pb.postraiter(1);
      pb.sauver();
      Process::exit();
    }
 
  return true;
}
 
 
/*! @brief Returns 1 if it is necessary to stop the computation for various reasons:
 *
 *         - the final time has been reached
 *         - the maximum number of time steps has been exceeded
 *         - the steady state has been reached
 *         - file-based stop indicator (see int Schema_Temps_base::stop_lu())
 *     Returns 0 otherwise.
 *
 * @return (int) 1 if the computation should stop, 0 otherwise.
 */
int Schema_Temps_base::stop() const
{
  if (temps_final_atteint())
    {
      Cerr << "---------------------------------------------------------"
           << finl
           << "The problem " << pb_base().le_nom()
           << " wants to stop : final time reached"
           << finl << finl;
      return 1;
    }
 
  if ((nb_pas_dt_max_atteint()) || nb_pas_dt_max_==0)
    {
      Cerr << "---------------------------------------------------------"
           << finl
           << "The problem " << pb_base().le_nom()
           << " wants to stop : the maximum number of time steps reached"
           << finl << finl;
      return 1;
    }
 
  if (temps_cpu_max_atteint())
    {
      Cerr << "---------------------------------------------------------"
           << finl
           << "The problem " << pb_base().le_nom()
           << " wants to stop : max cpu time reached"
           << finl << finl;
      return 1;
    }
 
  if (!get_disable_stop())
    {
      if (stop_lu())
        {
          Cerr << "---------------------------------------------------------"
               << finl
               << "The problem " << pb_base().le_nom()
               << " wants to stop : stop file detected"
               << finl << finl;
          return 1;
        }
    }
  return 0;
}
 
 
Probleme_base& Schema_Temps_base::pb_base()
{
  return mon_probleme.valeur();
}
 
const Probleme_base& Schema_Temps_base::pb_base() const
{
  return mon_probleme.valeur();
}
 
void Schema_Temps_base::ajouter_inertie(Matrice_Base& mat_morse,DoubleTab& secmem,const Equation_base& eqn) const
{
  // implementation for euler_implicite
  double dt=pas_de_temps();
  // We do not penalize the matrices and the right-hand sides even in the
  // Dirichlet, symmetry cases
  int pen=0;
  eqn.solv_masse().ajouter_masse(dt,mat_morse,pen); //important order for PolyMAC_MPFA
  const bool use_old_volumes = eqn.domaine_dis().domaine().deformable();
  eqn.solv_masse().ajouter_masse(dt, secmem, eqn.inconnue().passe(), pen, use_old_volumes);
}
 
void Schema_Temps_base::ajouter_blocs(matrices_t matrices, DoubleTab& secmem, const Equation_base& eqn, const tabs_t& semi_impl) const
{
  eqn.solv_masse().ajouter_blocs(matrices, secmem, pas_de_temps(), semi_impl, 1);
}
 
/*! @brief Updates stationnaire_atteint_ and residu_ (criterion: residu_ < seuil_statio_).
 *
 * @param (double& t) the new current time
 */
void Schema_Temps_base::update_critere_statio(const DoubleTab& tab_critere, Equation_base& equation)
{
  DoubleVect& residu_equation = equation.get_residu();
  // Depending on the size of the residu_equation array
  // we take the max of tab_critere over the whole array or by column
  int size = residu_equation.size_array();
  if (size==0)
    {
      Cerr << "Error in Schema_Temps_base::update_critere_statio" << finl;
      Cerr << "Array residu_equation has a null size" << finl;
      for (int i=0; i<pb_base().nombre_d_equations(); i++)
        if (pb_base().equation(i).get_residu().size_array()==0)
          {
            Cerr << "The equation of kind " << pb_base().equation(i).que_suis_je() << " didn't size the array." << finl;
            Cerr << "Please, contact TRUST support." << finl;
            exit();
          }
    }
  else if (size==1)
    {
      if(norm_residu_ == "max")
        residu_equation(0) = mp_max_abs_vect(tab_critere);
      else if (norm_residu_ == "L2" || norm_residu_ == "l2")
        residu_equation(0) = mp_norme_vect(tab_critere);
      else
        {
          Cerr << "Schema_Temps_base::update_critere_statio : only norm max and norm L2 are allowed to compute residuals ("
               << norm_residu_ << " not understood)" << finl;
          Process::exit();
        }
    }
  else
    {
      if(norm_residu_ == "max")
        mp_max_abs_tab(tab_critere, residu_equation);
      else if (norm_residu_ == "L2" || norm_residu_ == "l2")
        mp_norme_tab(tab_critere, residu_equation);
      else
        {
          Cerr << "Schema_Temps_base::update_critere_statio : only norm max and norm L2 are allowed to compute residuals ("
               << norm_residu_ << " not understood)" << finl;
          Process::exit();
        }
    }
  equation.set_residuals(tab_critere);
  // Compute residu_initial_equation over the first 5 time steps
  if (seuil_statio_relatif_deconseille_ == 1)
    {
      DoubleVect& residu_initial_equation = equation.residu_initial();
      if (nb_pas_dt()<6)
        {
          for (int i=0; i<size; i++)
            residu_initial_equation(i) = residu_equation(i);
          // We no longer take the max since it can be large at the first time step:
          // residu_initial_equation(i) = std::max(residu_initial_equation(i), residu_equation(i));
        }
      // Normalize residu_equation by residu_initial_equation
      for (int i=0; i<size; i++)
        if (residu_initial_equation(i)>0)
          residu_equation(i) /= residu_initial_equation(i);
    }
  else if (seuil_statio_relatif_deconseille_ == 2)
    {
      const double max_var = mp_max_abs_vect(equation.inconnue().futur());
      const double min_var = mp_min_abs_vect(equation.inconnue().futur());
      residu_equation /= max_var - min_var + 1e-2;
    }
  if (!equation.disable_equation_residual())
    {
      //double equation_residual = mp_max_abs_vect(residu_equation);
      double equation_residual = local_max_abs_vect(residu_equation);
      residu_ = std::max(residu_, equation_residual);
      set_stationnaire_atteint() *= (equation_residual < seuil_statio_);
    }
}
 
// Print the current time taking into account dt
// for the precision of the output (used in operators)
void Schema_Temps_base::imprimer_temps_courant(SFichier& os) const
{
  int precision_actuelle=os.get_precision();
  int precision_temps;
  if (!est_egal(temps_courant(),0.))
    precision_temps=std::max( precision_actuelle, (int)(2+log10(1/std::fabs(pas_de_temps()))+(int)(log10(std::fabs(temps_courant())))) );
  else
    precision_temps=std::max( precision_actuelle, (int)(2+log10(1/std::fabs(pas_de_temps()))) );
  os.precision(precision_temps);
  os << temps_courant();
  os.precision(precision_actuelle);
}
 
/*! @brief Writes the .progress file (estimated remaining CPU time).
 *
 */
void Schema_Temps_base::write_progress(bool init)
{
  if (je_suis_maitre() && !disable_progress())
    {
      if (init)
        {
          if (!progress_.is_open())
            progress_.ouvrir(nom_du_cas() + ".progress");
        }
      else
        {
          if (schema_impr())
            {
              if ((residu_ > 0) && (residu_old_slope_ > 0))
                cumul_slope_ += (log(residu_) - log(residu_old_slope_)) / dt_;
              residu_old_slope_ = residu_;
            }
          if (schema_impr() && (nb_pas_dt() > 0 && pas_de_temps() > 0))
            {
              // Compute the average CPU time per time step; drawback: it can vary strongly over time in case of divergence or acceleration
              // But Statistiques does not allow retrieving the CPU time of the last time step (last_time called here would return the CPU time since the start of the time step)
              double nb_pas_selon_tmax = (temps_max() - temps_courant()) / pas_de_temps();
              double nb_pas_selon_nb_pas_dt_max = nb_pas_dt_max() - nb_pas_dt();
              double nb_pas_avant_fin = std::min(nb_pas_selon_tmax, nb_pas_selon_nb_pas_dt_max);
              //double seconds_to_finish  = nb_pas_avant_fin * cpu_per_timestep;
              double dpercent = (1. - nb_pas_avant_fin /
                                 (nb_pas_avant_fin +
                                  nb_pas_dt()));    // works even if it is tmax that limits
              // if the slope is >0 we are diverging ....
              if ((seuil_statio_ > 0) && (cumul_slope_ < -1e-20) &&
                  seuil_statio_ <
                  residu_) // in a pb_couple one of the problems may have (seuil_statio_ > residu)
                {
                  double distance = (-log(residu_ + 1e-20) + log(seuil_statio_)) / (cumul_slope_) * nb_pas_dt();
 
                  //Cerr<<distance<<" DDDDDD"<<finl;
                  double dpercent2 = temps_courant() / (temps_courant() + distance);
                  dpercent = std::max(dpercent, dpercent2);
                }
              int percent = int(dpercent * 100);
 
              if (limpr())
                {
                  double seconds_to_finish = statistics().get_time_since_last_open(STD_COUNTERS::total_execution_time) * (1. - dpercent) / dpercent;
                  int integer_limit = (int) (pow(2.0, (double) ((sizeof(int) * 8) - 1)) - 1);
                  if (seconds_to_finish < integer_limit)
                    {
                      int h = int(seconds_to_finish / 3600);
                      int mn = int((seconds_to_finish - 3600 * h) / 60);
                      int s = int(seconds_to_finish - 3600 * h - 60 * mn);
 
                      Cout << finl << "Estimated CPU time to finish the run (according to "
                           << (nb_pas_selon_tmax < nb_pas_selon_nb_pas_dt_max ? "tmax" : "nb_pas_dt_max")
                           << " value) : ";
                      if (seconds_to_finish < 1)
                        Cout << seconds_to_finish << " s";
                      else
                        Cout << h << "h" << mn << "mn" << s << "s";
 
                      Cout << ". Progress: " << (percent) << finl;
                    }
                }
              if (!disable_progress() && Process::je_suis_maitre())
                progress_ << (percent) << finl;
            }
        }
    }
}
 
/*! @brief Writes the .dt_ev file.
 *
 */
static SFichier dt_ev_;
static bool header_complete = false;
void open_dt_ev(IOS_OPEN_MODE mode)
{
  if (!dt_ev_.is_open())
    {
      Nom fichier(Objet_U::nom_du_cas() + ".dt_ev");
      dt_ev_.ouvrir(fichier, mode);
      dt_ev_.setf(ios::scientific);
    }
}
void Schema_Temps_base::write_dt_ev(bool init)
{
  if (je_suis_maitre() && !disable_dt_ev())
    {
      if (init)
        {
          Nom fichier(nom_du_cas() + ".dt_ev");
          struct stat f;
          // Initialize the .dt_ev file if it does not exist or if this is a fresh computation without restart
          if ((nb_pas_dt_ == 0) && ((stat(fichier, &f)) || !(pb_base().reprise_effectuee() == 1)))
            {
              if (schema_impr())
                {
                  open_dt_ev(ios::out);
                  dt_ev_ << "# temps\t\t dt\t\t facsec\t\t residu=max|Ri|\t dt_stab\t ";
                }
              open_dt_ev(ios::app);
              for (int i = 0; i < pb_base().nombre_d_equations(); i++)
                dt_ev_ << pb_base().equation(i).expression_residu();
            }
          else
            {
              if (schema_impr())
                open_dt_ev(ios::app);
            }
        }
      else
        {
          open_dt_ev(ios::app);
          if (schema_impr())
            {
              dt_ev_ << finl << temps_courant_ << "\t " << dt_ << "\t " << facsec_ << "\t " << residu_ << "\t "
                     << dt_stab_ << "\t ";
              header_complete = true;
            }
          for (int i = 0; i < pb_base().nombre_d_equations(); i++)
            pb_base().equation(i).imprime_residu(dt_ev_);
        }
    }
}
 
/*! @brief Closes the .dt_ev file.
 *
 */
void Schema_Temps_base::finir() const
{
  if (je_suis_maitre() && dt_ev_.is_open())
    {
      if (header_complete) dt_ev_ << finl;
      dt_ev_.close();
      header_complete = false;
    }
}