Sch_CN_EX_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 <Sch_CN_EX_iteratif.h>
#include <Probleme_base.h>
#include <Navier_Stokes_std.h>
#include <Param.h>
#include <SFichier.h>
 
Implemente_instanciable(Sch_CN_EX_iteratif,"Sch_CN_EX_iteratif",Sch_CN_iteratif);
// XD Sch_CN_EX_iteratif Sch_CN_iteratif Sch_CN_EX_iteratif INHERITS_BRACE This keyword also describes a Crank-Nicholson
// XD_CONT method of second order accuracy but here, for scalars, because of instablities encountered when dt>dt_CFL,
// XD_CONT the Crank Nicholson scheme is not applied to scalar quantities. Scalars are treated according to
// XD_CONT Euler-Explicite scheme at the end of the CN treatment for velocity flow fields (by doing p Euler explicite
// XD_CONT under-iterations at dt<=dt_CFL). Parameters are the sames (but default values may change) compare to the
// XD_CONT Sch_CN_iterative scheme plus a relaxation keyword: niter_min (2 by default), niter_max (6 by default),
// XD_CONT niter_avg (3 by default), facsec_max (20 by default), seuil (0.05 by default)
 
Sortie& Sch_CN_EX_iteratif::printOn(Sortie& s) const
{
  return  Sch_CN_iteratif::printOn(s);
}
 
 
Entree& Sch_CN_EX_iteratif::readOn(Entree& s)
{
  seuil=0.05;
  niter_min=2;
  niter_max=6;
  niter_avg=3;
  facsec_max=20;
 
  Schema_Temps_base::readOn(s) ;
 
  return s;
}
 
void Sch_CN_EX_iteratif::set_param(Param& param) const
{
  param.ajouter("omega",&omega);  // XD attr omega floattant omega OPT relaxation factor (0.1 by default)
  Sch_CN_iteratif::set_param(param);
}
 
/*! @brief Only takes into account the Navier-Stokes equation of the problem (if there is one). The other equations are not limiting: they perform as many sub-time-steps as needed.
 *
 *  Code shamelessly copied from Probleme_base::calculer_pas_de_temps()
 *
 *
 */
void Sch_CN_EX_iteratif::mettre_a_jour_dt_stab()
{
  imprimer(Cout);
  //  dt_=dt_max_; The upper bound by dt_max_ is applied in corriger_dt_calcule
  dt_stab_=DMAXFLOAT;
  const Probleme_base& prob=pb_base();
 
  for(int i=0; i<prob.nombre_d_equations(); i++)
    {
      const Equation_base& eqn=prob.equation(i);
      if (sub_type(Navier_Stokes_std,eqn))
        {
          dt_stab_=std::min(dt_stab_,eqn.calculer_pas_de_temps());
          break; // At most one Navier-Stokes equation per problem
        }
    }
}
bool Sch_CN_EX_iteratif::iterateTimeStepOnEquation(int i,bool& converged)
{
  Probleme_base& pb = pb_base();
  Equation_base& eqn = pb.equation(i);
 
  // Choose based on whether the equation is of Navier_Stokes type or not
  if (sub_type(Navier_Stokes_std,eqn))
    return iterateTimeStepOnNS(i,converged);
  else
    return iterateTimeStepOnOther(i,converged);
}
 
// Very similar to Sch_CN_iteratif::iterateTimeStepOnEquation but applies
// a relaxation factor omega between iterations.
 
bool Sch_CN_EX_iteratif::iterateTimeStepOnNS(int i,bool& converged)
{
 
  Probleme_base& pb = pb_base();
  Equation_base& eqn = pb.equation(i);
 
  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);
  DoubleTab old(intermediaire);
 
  // Impose the Dirichlet BCs at the intermediate time.
  // Indeed, the diffusion operators only use the unknown
  // and do not read the BCs.
  // This ensures in particular flux equality in a
  // coupled thermal VEF problem with Champ_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 derivative dudt for the intermediate value of the unknown.
  // Workaround: since operators take the present value by default,
  // the unknown is temporarily advanced.
 
  eqn.inconnue().avancer();
  eqn.derivee_en_temps_inco(dudt);
  eqn.inconnue().reculer();
 
  // Update the unknown values at intermediate and final times:
  // intermediaire = present + dt_intermediaire * dudt;
  // intermediaire = (1-omega)*new_intermediaire + omega*old_intermediaire
  // final =  present + dt_final * dudt;
  intermediaire = dudt;
  intermediaire*= dt_intermediaire;
  intermediaire+= present;
  intermediaire*= (1-omega);
  old*= omega;
  intermediaire+= old;
  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*=-1;
  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;
 
}
 
// Perform several explicit Euler time steps with facsec<=1.
// One step is required to land exactly on the intermediate time
// and another exactly on the final time.
 
bool Sch_CN_EX_iteratif::iterateTimeStepOnOther(int i,bool& converged)
{
  Probleme_base& pb = pb_base();
  Equation_base& eqn = pb.equation(i);
 
  double dt_eq=eqn.calculer_pas_de_temps();
 
  double temps_intermediaire=temps_futur(1);
  double temps_final=temps_futur(2);
  double dt_1=temps_intermediaire-temps_courant();
  double dt_2=temps_final-temps_intermediaire;
 
  DoubleTab& present = eqn.inconnue().valeurs();
  DoubleTab& intermediaire = eqn.inconnue().valeurs(temps_intermediaire);
  DoubleTab& final = eqn.inconnue().valeurs(temps_final);
 
  DoubleTab dIdt(present); // Pour dimensionner.
 
  // See Schema_Temps_base::limpr for information on epsilon and modf
  double nb_dt_1, nb_dt_2;
  modf(ceil(dt_1/dt_eq), &nb_dt_1); // Number of time steps in the first and second interval.
  modf(ceil(dt_2/dt_eq), &nb_dt_2);
 
  double dt_eq_1=dt_1/nb_dt_1; // Duration of time steps in the first and second interval.
  double dt_eq_2=dt_2/nb_dt_2;
 
  // Computations on the first time step
  intermediaire=present;
  for (double k=0.; k<nb_dt_1; k=k+1.)
    {
 
      // Compute the derivative dIdt for the intermediate value of the unknown.
      // Workaround: since operators take the present value by default,
      // the unknown is temporarily advanced.
      eqn.inconnue().avancer();
      eqn.derivee_en_temps_inco(dIdt);
      eqn.inconnue().reculer();
 
      // intermediaire = intermediaire + dt_eq_1 * dIdt
      dIdt*=dt_eq_1;
      intermediaire+= dIdt;
      intermediaire.echange_espace_virtuel();
 
    }
 
  // Impose the BCs at the intermediate time
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_intermediaire);
 
 
  // Computations (same) on the second time step
  final=intermediaire;
  for (double k=0.; k<nb_dt_2; k=k+1.)
    {
 
      // Compute the derivative dIdt for the final value of the unknown.
      eqn.inconnue().avancer();
      eqn.inconnue().avancer();
      eqn.derivee_en_temps_inco(dIdt);
      eqn.inconnue().reculer();
      eqn.inconnue().reculer();
 
      // intermediaire = intermediaire + dt_eq_2 * dIdt
      dIdt*=dt_eq_2;
      final+= dIdt;
      final.echange_espace_virtuel();
 
    }
  dIdt=final;
  dIdt-=present;
  dIdt/=dt_2;
  update_critere_statio(dIdt, eqn);
  // Impose the BCs at the final time
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),temps_final);
 
  converged=true;
  return true;
 
}
 
 
void Sch_CN_EX_iteratif::ajuster_facsec(type_convergence cv)
{
 
  if (facsec_!=last_facsec) // Adjust only once per initTimeStep.
    return;
 
  switch (cv)
    {
    case DIVERGENCE:
      // the divergence criterion has been reached
      facsec_=last_facsec/1.5;
      omega=0.5*(0.5+omega); // drives omega towards 0.5
      break;
    case NON_CONVERGENCE:
      // neither the divergence criterion nor the convergence criterion
      // was reached in niter_max iterations
      facsec_=last_facsec/1.5;
      omega=0.5*(0.5+omega); // (0.5+3*omega)/4; // drives omega towards 0.5, more slowly
      break;
    case CONVERGENCE_LENTE:
      // the convergence criterion was reached in more than niter_avg
      // iterations
      facsec_=last_facsec*0.9;
      if(omega>0.1)
        omega*=0.95;
      break;
    case CONVERGENCE_RAPIDE:
      // the convergence criterion was reached in fewer than niter_avg
      // iterations
      facsec_=last_facsec*1.1;
      if(omega>0.1)
        omega*=0.75;
      break;
    case CONVERGENCE_OK:
      // the convergence criterion was reached in exactly niter_avg iterations
      facsec_=last_facsec*1.01;
      break;
    default:
      break;
    }
  facsec_=std::min(facsec_,facsec_max);
}