Sch_CN_iteratif.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 <Navier_Stokes_std.h>
#include <Sch_CN_iteratif.h>
#include <communications.h>
#include <Probleme_base.h>
#include <Equation_base.h>
#include <SFichier.h>
#include <TRUSTTab.h>
#include <Param.h>
 
Implemente_instanciable(Sch_CN_iteratif, "Sch_CN_iteratif", Schema_Temps_base);
// XD Sch_CN_iteratif schema_temps_base Sch_CN_iteratif INHERITS_BRACE The Crank-Nicholson method of second order
// XD_CONT accuracy. A mid-point rule formulation is used (Euler-centered scheme). The basic scheme is: $$u(t+1) = u(t)
// XD_CONT + du/dt(t+1/2)*dt$$ The estimation of the time derivative du/dt at the level (t+1/2) is obtained either by
// XD_CONT iterative process. The time derivative du/dt at the level (t+1/2) is calculated iteratively with a simple
// XD_CONT under-relaxations method. Since the method is implicit, neither the cfl nor the fourier stability criteria
// XD_CONT must be respected. The time step is calculated in a way that the iterative procedure converges with the less
// XD_CONT iterations as possible. NL2 Remark : for stationary or RANS calculations, no limitation can be given for time
// XD_CONT step through high value of facsec_max parameter (for instance : facsec_max 1000). In counterpart, for LES
// XD_CONT calculations, high values of facsec_max may engender numerical instabilities.
 
 
Sortie& Sch_CN_iteratif::printOn(Sortie& s) const { return Schema_Temps_base::printOn(s); }
 
Entree& Sch_CN_iteratif::readOn(Entree& s) { return Schema_Temps_base::readOn(s); }
 
void Sch_CN_iteratif::ajuster_facsec(type_convergence cv)
{
 
  if (facsec_!=last_facsec) // Only adjust once per initTimeStep.
    return;
 
  switch (cv)
    {
    case DIVERGENCE:
      // the divergence criterion has been reached
      facsec_=last_facsec*0.8;
      break;
    case NON_CONVERGENCE:
      // neither the divergence nor the convergence criterion
      // was reached in niter_max iterations
      facsec_=last_facsec*0.9;
      break;
    case CONVERGENCE_LENTE:
      // the convergence criterion was reached in more than niter_avg
      // iterations
      facsec_=last_facsec*0.99;
      break;
    case CONVERGENCE_RAPIDE:
      // the convergence criterion was reached in fewer than niter_avg
      // iterations
      facsec_=last_facsec*1.01;
      break;
    case CONVERGENCE_OK:
      // the convergence criterion was reached in exactly niter_avg iterations
      break;
    default:
      break;
    }
  facsec_=std::min(facsec_,facsec_max);
}
void Sch_CN_iteratif::set_param(Param& param) const
{
  param.ajouter("seuil",&seuil);  // XD attr seuil floattant seuil OPT criteria for ending iterative process (Max( ||
  // XD_CONT u(p) - u(p-1)||/Max || u(p) ||) < seuil) (0.001 by default)
  param.ajouter("niter_min",&niter_min);  // XD attr niter_min entier niter_min OPT minimal number of p-iterations to
  // XD_CONT satisfy convergence criteria (2 by default)
  param.ajouter("niter_max",&niter_max);  // XD attr niter_max entier niter_max OPT number of maximum p-iterations
  // XD_CONT allowed to satisfy convergence criteria (6 by default)
  param.ajouter("niter_avg",&niter_avg);  // XD attr niter_avg entier niter_avg OPT threshold of p-iterations (3 by
  // XD_CONT default). If the number of p-iterations is greater than niter_avg, facsec is reduced, if lesser than
  // XD_CONT niter_avg, facsec is increased (but limited by the facsec_max value).
  param.ajouter("facsec_max",&facsec_max);// XD attr facsec_max floattant facsec_max OPT maximum ratio allowed between
  // XD_CONT dynamical time step returned by iterative process and stability time returned by CFL condition (2 by
  // XD_CONT default).
  Schema_Temps_base::set_param(param);
}
 
 
////////////////////////////////
//                            //
// Scheme characteristics     //
//                            //
////////////////////////////////
 
/*! @brief Returns the number of time values to keep.
 *
 * Here: n, n+1/2 and n+1, so 3.
 *
 */
int Sch_CN_iteratif::nb_valeurs_temporelles() const
{
  return 3;
}
 
/*! @brief Returns the number of future time values.
 *
 * Here: n+1/2 and n+1, so 2.
 *
 */
int Sch_CN_iteratif::nb_valeurs_futures() const
{
  return 2;
}
 
/*! @brief Returns the time at the i-th future value.
 *
 */
double Sch_CN_iteratif::temps_futur(int i) const
{
  assert(i>0 && i<=2);
  if (i==2)
    return temps_courant()+pas_de_temps();
  else
    return temps_courant()+pas_de_temps()/2;
}
 
/*! @brief Returns the time that fields should use when calling valeurs().
 *
 *     Here: t(n+1/2).
 *
 */
double Sch_CN_iteratif::temps_defaut() const
{
  return temps_futur(1);
}
 
/////////////////////////////////////////
//                                     //
// End of scheme characteristics       //
//                                     //
/////////////////////////////////////////
 
bool Sch_CN_iteratif::initTimeStep(double dt)
{
 
  if (nb_pas_dt()==0)
    mettre_a_jour_dt_stab(); // otherwise already called by validateTimeStep
  //No longer necessary since dt_ is now updated in Schema_Temps_base::computeTimeStep(bool& stop)
  //facsec_=dt/dt_;
  last_facsec=facsec_;
 
  iteration=0;
 
  return Schema_Temps_base::initTimeStep(dt);
}
 
/*! @brief
 *
 */
bool Sch_CN_iteratif::iterateTimeStep(bool& converged)
{
 
  Probleme_base& pb = pb_base();
  const int nb_eqn=pb.nombre_d_equations();
 
  SFichier fic;
  fic.ouvrir("dt_CN",ios::app);
 
  if (je_suis_maitre())
    {
      fic << "Sch_CN_iteratif : pb= " <<  pb.le_nom()
          << " present= " << temps_courant()
          << " dt= " << dt_
          << " iteration= " << iteration
          << " last_facsec= " << last_facsec;
    }
  converged=true;
  bool diverged=false;
 
  for(int i=0; i<nb_eqn && !diverged; i++)
    {
      bool cv;
      diverged = diverged || !iterateTimeStepOnEquation(i,cv);
      converged = converged && cv;
    }
 
  iteration ++;
  if (iteration<niter_min)
    converged=false; // Continue iterating anyway
 
  // Sanity check...
  assert(!(converged&&diverged));
 
  // If convergence
  if (converged)
    {
      if (iteration<niter_avg) // Convergence too fast, increase facsec
        ajuster_facsec(CONVERGENCE_RAPIDE);
      else if (iteration>niter_avg) // Convergence too slow, reduce facsec
        ajuster_facsec(CONVERGENCE_LENTE);
      else
        ajuster_facsec(CONVERGENCE_OK);
      if (je_suis_maitre())
        {
          fic << " facsec= " << facsec_
              << " result= CONVERGENCE" << finl;
        }
      fic.close();
      return true;
    }
 
  // If divergence
  else if (diverged)
    {
      ajuster_facsec(DIVERGENCE);
      if (je_suis_maitre())
        {
          fic << " facsec= " << facsec_
              << " result= DIVERGENCE" << finl;
        }
      return false;
    }
 
  // If not converged fast enough
  else if (iteration==niter_max-1)
    {
      ajuster_facsec(NON_CONVERGENCE);
      if (je_suis_maitre())
        {
          fic << " facsec= " << facsec_
              << " result= NON_CONVERGENCE" << finl;
        }
      return false;
    }
 
  // Otherwise, still converging
  if (je_suis_maitre())
    {
      fic << " facsec= " << facsec_
          << " result= CONTINUE" << finl;
    }
  return true;
}
 
 
/*! @brief Computes one iteration of the resolution on equation i.
 *
 * Computes u(n+1/2,p+1)=u(n)+f(u(n+1/2,p))*dt/2
 *  and u(n+1,p+1)=u(n)+f(u(n+1/2,p))*dt
 *  where f gives du/dt as a function of u.
 *  Returns true in converged if no change between iterations, false otherwise.
 *  Returns true if OK to continue iterating, false otherwise (diverged or too many iterations).
 *
 */
bool Sch_CN_iteratif::iterateTimeStepOnEquation(int i,bool& converged)
{
 
  Probleme_base& pb = pb_base();
  Equation_base& eqn = pb.equation(i);
  if (eqn.equation_non_resolue())
    {
      Cout<< "====================================================" << finl;
      Cout<< eqn.que_suis_je()<<" equation is not solved."<<finl;
      Cout<< "====================================================" << finl;
      // Compute the derivative once to obtain the boundary fluxes
      if (eqn.schema_temps().nb_pas_dt()==0)
        {
          DoubleTab inconnue_valeurs(eqn.inconnue().valeurs());
          eqn.derivee_en_temps_inco(inconnue_valeurs);
        }
      converged=true;
      return true;
    }
  double temps_intermediaire=temps_futur(1);
  double temps_final=temps_futur(2);
  double dt_intermediaire=temps_intermediaire-temps_courant();
  double dt_final=temps_final-temps_courant();
 
  DoubleTab& present = eqn.inconnue().valeurs();
  DoubleTab& intermediaire = eqn.inconnue().valeurs(temps_intermediaire);
  DoubleTab& final = eqn.inconnue().valeurs(temps_final);
 
  DoubleTab dudt(present);
 
  DoubleTab delta(intermediaire);
  delta*=-1;
 
  // Impose Dirichlet BCs at the intermediate time.
  // Indeed, diffusion operators only use the unknown
  // and do not read the BCs.
  // This ensures in particular equal fluxes for coupled
  // thermal VEF problems with Chap_front_contact_VEF, even before convergence.
  // WEC: /!\ the velocity at the intermediate time
  // is not necessarily divergence-free.
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_intermediaire);
 
  // Compute the time derivative dudt for the intermediate value of the unknown.
  // Trick: since operators take the present value by default,
  // we temporarily advance the unknown.
 
  eqn.inconnue().avancer();
  eqn.derivee_en_temps_inco(dudt);
  eqn.inconnue().reculer();
 
  // Update the unknown values at the intermediate and final times:
  // intermediaire = present + dt_intermediaire * dudt;
  // final =  present + dt_final * dudt;
  intermediaire = dudt;
  intermediaire*= dt_intermediaire;
  intermediaire+= present;
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_intermediaire);
  intermediaire.echange_espace_virtuel();
  final = dudt;
  final*= dt_final;
  final+= present;
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_final);
  final.echange_espace_virtuel();
 
  delta+=intermediaire; // delta = contains u(n+1/2,p+1) - u(n+1/2,p)
 
  // If the equation has diverged
  if (divergence(present,intermediaire,delta,iteration))
    {
      converged=false;
      return false;
    }
 
  // If the equation has converged
  if (convergence(present,intermediaire,delta,iteration))
    {
      converged=true;
      delta=final;
      delta-=present;
      delta/=dt_final;
      update_critere_statio(delta, eqn);
    }
  else
    {
      set_stationnaire_atteint()=0;
      converged=false;
    }
  return true;
}
 
 
// WEC: this could be coded, but it doesn't seem useful.
// It would be a loop over IterateTimeStepOnEquation
int Sch_CN_iteratif::faire_un_pas_de_temps_eqn_base(Equation_base&)
{
  Cerr << "Sch_CN_iteratif::faire_un_pas_de_temps_eqn_base non code!" << finl;
  exit();
  return 0;
}
 
 
/*! @brief Indicates whether the iterative computation has converged.
 *
 * Convergence criterion used:
 *     || u(n+1,p+1) - u(n+1,p) ||  <  seuil * || u(n+1/2,p+1) ||
 *  This is equivalent to:
 *     || u(n+1,p) - u(n+1) || < seuil * || (Id-(dt/2).(df/du))^-1 || * || u(n+1/2,p+1) ||
 *  where u(n+1) is the exact solution at n+1,
 *     df/du is the Jacobian of f(u), taken at u(n+1/2),
 *     and the norms are compatible between matrices and vectors (here the infinity norm).
 *
 * @param u0 The unknown at the beginning of the time interval u(n). Also the first estimate u(n+1/2,0).
 * @param up1 Estimate of the unknown at tn+1/2 at iteration p+1: u(n+1/2,p+1).
 * @param delta u(n+1/2,p+1) - u(n+1/2,p).
 * @return true if converged, false otherwise.
 */
bool Sch_CN_iteratif::convergence(const DoubleTab& u0, const DoubleTab& up1, const DoubleTab& delta, int p) const
{
  double a = mp_max_abs_vect(delta);
  double b = mp_max_abs_vect(up1);
  int resu = (2. * a) < (seuil * b);
  envoyer_broadcast(resu, 0); // to ensure all processes make the same decision
  return resu;
}
 
/*! @brief Indicates whether the iterative computation has diverged.
 *
 * @param u0 The unknown at the beginning of the time interval u(n). Also the first estimate u(n+1/2,0).
 * @param up1 Estimate of the unknown at tn+1/2 at iteration p+1: u(n+1/2,p+1).
 * @param delta u(n+1/2,p+1) - u(n+1/2,p).
 * @return true if diverged, false otherwise.
 */
bool Sch_CN_iteratif::divergence(const DoubleTab& u0, const DoubleTab& up1, const DoubleTab& delta, int p) const
{
  // WEC: could be improved...
  return false;
}