Schema_Euler_Implicite_Stationnaire.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 <Equation_base.h>
#include <Probleme_base.h>
#include <Probleme_Couple.h>
#include <Milieu_base.h>
#include <Param.h>
#include <Debog.h>
#include <LecFicDiffuse.h>
#include <communications.h>
#include <string>
#include <Schema_Euler_Implicite_Stationnaire.h>
 
Implemente_instanciable(Schema_Euler_Implicite_Stationnaire,"Schema_Euler_Implicite_Stationnaire|Implicit_Euler_Steady_Scheme",Schema_Euler_Implicite);
// XD Implicit_Euler_Steady_Scheme schema_implicite_base schema_euler_implicite_stationnaire BRACE This is the Implicit
// XD_CONT Euler scheme using a dual time step procedure (using local and global dt) for steady problems. The only
// XD_CONT possible solver choice for this scheme is the implicit_steady solver. The stationarity criterion checks all
// XD_CONT equations (velocity, turbulence, temperature, etc.) using the global seuil_statio threshold. For RANS
// XD_CONT turbulence models (k-epsilon, k-omega), turbulence variables may exhibit persistent low-amplitude
// XD_CONT oscillations that prevent the stationarity criterion from being reached. In such cases, it is recommended to
// XD_CONT set nb_pas_dt_max as the primary stopping criterion and monitor residuals manually. The local time step used
// XD_CONT as pseudo-inertia is dt_loc = steady_security_factor * facsec * dt_CFL. When facsec_max is set to a finite
// XD_CONT value, facsec grows dynamically (based on residual reduction) up to facsec_max, accelerating convergence
// XD_CONT similarly to the standard implicit scheme. Without facsec_max, facsec stays fixed at its initial value
// XD_CONT (default 1). The pseudo-time displayed in outputs (TIME in .dt_ev) equals the iteration count (dt=1 per
// XD_CONT step), making it comparable across meshes and convection schemes.
// XD attr facsec_max floattant facsec_max OPT Maximum ratio allowed between local pseudo-time step and stability time
// XD_CONT returned by the CFL condition; default is 1 (facsec stays fixed). Setting a value greater than 1 enables
// XD_CONT dynamic facsec growth.
 
 
Sortie& Schema_Euler_Implicite_Stationnaire::printOn(Sortie& s) const
{
  return  Schema_Implicite_base::printOn(s);
}
 
 
Entree& Schema_Euler_Implicite_Stationnaire::readOn(Entree& s)
{
  // Set default facsec_max to 1 for the steady scheme (parent default is DMAXFLOAT).
  // With facsec_max = 1, facsec stays fixed at 1 and calcul_fac_sec is never called.
  // The user can override this by explicitly setting facsec_max > 1 in the data file.
  facsec_max_ = 1.0;
  Schema_Euler_Implicite::readOn(s);
  if(!le_solveur)
    {
      Cerr << "A solver must be selected." << finl;
      Cerr << "Syntax : " << finl
           << "Solveur solver_name [ solver parameters ] " << finl;
      Process::exit();
    }
  if (diffusion_implicite())
    {
      Cerr << "diffusion_implicite option cannot be used with an implicit time scheme." << finl;
      Process::exit();
    }
 
  if (le_solveur->le_nom()!="Implicit_steady")
    {
      Cerr << "You can only select the Implicit_steady solver" << finl;
      Cerr << "Syntax : " << finl
           << "Solveur Implicit_steady [ solver parameters ] " << finl;
      Process::exit();
    }
 
  return s;
}
 
void Schema_Euler_Implicite_Stationnaire::set_param(Param& param) const
{
  Schema_Euler_Implicite::set_param(param);
 
 
}
 
bool Schema_Euler_Implicite_Stationnaire::isStationary() const
{
  return Schema_Temps_base::isStationary();
}
 
 
int Schema_Euler_Implicite_Stationnaire::mettre_a_jour()
{
  // Skip Schema_Euler_Implicite::mettre_a_jour() intentionally to avoid
  // calling calcul_fac_sec() unconditionally.
  // facsec_ acts as a dynamic CFL multiplier for dt_loc_ (see calculer_pas_de_temps_local_pb).
  // Default facsec_max = 1: facsec_ stays fixed at 1 (no growth).
  // Growth is enabled only when the user sets facsec_max > 1.
  if (facsec_max_ > 1.0)
    calcul_fac_sec(residu_, residu_old_, facsec_);
  return Schema_Temps_base::mettre_a_jour();
}
 
int Schema_Euler_Implicite_Stationnaire::reprendre(Entree& is)
{
  // Restore parent data
  int ok = Schema_Euler_Implicite::reprendre(is);
  // Force pressure matrix reassembly after restart —
  // dt_loc_ from previous run is no longer valid
  dt_loc_.resize(0);
  return ok;
}
 
void Schema_Euler_Implicite_Stationnaire::ajouter_inertie(Matrice_Base& mat_morse,DoubleTab& secmem,const Equation_base& eqn) const
{
  // Implicit Euler inertia: adds M/dt_loc to the matrix and M/dt_loc*U^n to the RHS
  DoubleVect dt_loc = dt_loc_;
 
  // No penalty on Dirichlet or symmetry boundary faces
  int pen=0;
  if (dt_loc.size()!= secmem.size()) // VEF velocity: size = nb_dim*nb_faces, while dt_loc has size nb_faces
    {
      int nb_dim=eqn.inconnue().dimension;
      if(dt_loc.size()* nb_dim == secmem.size())
        {
          DoubleVect  dt_loc_velocity(secmem);
          int j = 0;
          for(int i=0; i<dt_loc_velocity.size(); i++)
            {
              if(i != 0 && i%nb_dim == 0) j++;
              dt_loc_velocity[i]= dt_loc[j];
            }
          dt_loc_velocity.echange_espace_virtuel();
          eqn.solv_masse().ajouter_masse_dt_local(dt_loc_velocity,secmem,eqn.inconnue().passe(),pen);
          eqn.solv_masse().ajouter_masse_dt_local(dt_loc_velocity,mat_morse,pen);
        }
      else if (dt_loc.size()*2==secmem.size())
        {
          // k-eps: 2 DOFs per face
          DoubleVect  dt_loc_velocity(secmem);
          int j = 0;
          for(int i=0; i<dt_loc_velocity.size(); i++)
            {
              if(i != 0 && i%2 == 0) j++;
              dt_loc_velocity[i]= dt_loc[j];
            }
          dt_loc_velocity.echange_espace_virtuel();
          eqn.solv_masse().ajouter_masse_dt_local(dt_loc_velocity,secmem,eqn.inconnue().passe(),pen);
          eqn.solv_masse().ajouter_masse_dt_local(dt_loc_velocity,mat_morse,pen);
        }
      else
        {
          // dt_loc is face-sized (nb_faces), but secmem has a different size.
          // Known compatible cases:
          //   - VEF velocity:  secmem.size() == dt_loc.size() * nb_dim  (nb_dim components per face)
          //   - VEF k-eps/k-omega: secmem.size() == dt_loc.size() * 2   (k and eps/omega face-centered)
          // A mismatch that falls here almost always means the scheme is used with VDF
          // and a scalar equation (RANS k-eps, k-omega, temperature ...) whose unknowns
          // are element-centered (nb_elem DOFs) instead of face-centered (nb_faces DOFs).
          // Implicit_Euler_Steady_Scheme / Implicit_steady is only compatible with VEF
          // discretization for multi-equation problems (turbulence, thermohydraulics).
          Cerr << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << finl;
          Cerr << "Implicit_Euler_Steady_Scheme: incompatible equation detected." << finl;
          Cerr << "  Equation      : " << eqn.que_suis_je() << finl;
          Cerr << "  dt_loc.size() = " << dt_loc.size()   << "  (face-sized, from NS convection operator)" << finl;
          Cerr << "  secmem.size() = " << secmem.size()   << "  (unknown DOF count for this equation)" << finl;
          Cerr << finl;
          Cerr << "This scheme is ONLY compatible with VEF discretization when used" << finl;
          Cerr << "with RANS turbulence (k-eps, k-omega) or temperature equations." << finl;
          Cerr << "In VDF, turbulence and scalar unknowns are element-centered (nb_elem)" << finl;
          Cerr << "while dt_loc is face-centered (nb_faces), which are incompatible." << finl;
          Cerr << "Please switch to a VEF discretization or use a different time scheme." << finl;
          Cerr << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << finl;
          Process::exit();
        }
    }
  else
    {
      eqn.solv_masse().ajouter_masse_dt_local(dt_loc,secmem,eqn.inconnue().passe(),pen);
      eqn.solv_masse().ajouter_masse_dt_local(dt_loc,mat_morse,pen);
    }
}
 
void Schema_Euler_Implicite_Stationnaire::calculer_pas_de_temps_local_pb()
{
 
  pb_base().equation(0).calculer_pas_de_temps_locaux(dt_locaux_);
  for(int i=1; i< pb_base().nombre_d_equations(); i++)
    {
      DoubleTab dt_temp (dt_locaux_);
      pb_base().equation(i).calculer_pas_de_temps_locaux(dt_temp);
      for (int count=0; count<dt_locaux_.size(); count++)
        {
          if (dt_locaux_[count] > dt_temp[count])
            dt_locaux_[count] =  dt_temp[count];
        }
    }
 
  dt_locaux_ *= steady_security_factor_ * facsec_;
}
void Schema_Euler_Implicite_Stationnaire::initialize()
{
  // VDF is not supported for multi-equation or turbulent problems:
  //  - Thermohydraulics: nombre_d_equations() > 1 (NS + temperature)
  //  - RANS turbulence:  nombre_d_equations() == 1 BUT the NS equation is of type
  //    Navier_Stokes_Turbulent, which internally drives a transport equation (k-eps,
  //    k-omega ...) that Schema_Euler_Implicite solves via faire_un_pas_de_temps_eqn_base.
  //    The local dt_loc coupling is not validated for these cases in VDF and diverges.
  if (pb_base().discretisation().is_vdf())
    {
      bool unsupported = (pb_base().nombre_d_equations() > 1);
      if (!unsupported)
        unsupported = pb_base().equation(0).que_suis_je().contient("Turbulent");
      if (unsupported)
        {
          Cerr << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << finl;
          Cerr << "Implicit_Euler_Steady_Scheme is NOT compatible with VDF" << finl;
          Cerr << "for turbulent or multi-equation problems." << finl;
          Cerr << "  Discretization : VDF" << finl;
          Cerr << "  Equation(0)    : " << pb_base().equation(0).que_suis_je() << finl;
          Cerr << "  Nb equations   : " << pb_base().nombre_d_equations() << finl;
          Cerr << "RANS turbulence (k-eps, k-omega) and temperature equations are" << finl;
          Cerr << "not supported with this scheme in VDF: the local time-step" << finl;
          Cerr << "coupling diverges due to inconsistent DOF layout." << finl;
          Cerr << "Please switch to VEF discretization or use a standard time scheme." << finl;
          Cerr << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << finl;
          Process::exit();
        }
    }
  Schema_Euler_Implicite::initialize();
  calculer_pas_de_temps_local_pb();
  // Force pressure matrix reassembly at first iteration
  dt_loc_.resize(0);
}
 
void Schema_Euler_Implicite_Stationnaire::mettre_a_jour_dt_stab()
{
  imprimer(Cout);
  if(nb_pas_dt_ <1)
    dt_stab_=pb_base().calculer_pas_de_temps();
  calculer_pas_de_temps_local_pb();
  max_dt_loc_= dt_locaux_.mp_max_abs_vect();
 
  // Print dt_loc statistics for the first 5 steps to help the user
  // assess mesh heterogeneity and check local time step stability.
  // Skipped afterwards to avoid costly MPI reductions at every step.
  if (nb_pas_dt_ < 5)
    {
      double dt_loc_min = dt_locaux_.mp_min_vect();
      double dt_loc_max = max_dt_loc_;
      double dt_loc_sum = mp_somme_vect(dt_locaux_);
      int    dt_loc_n   = (int)mp_sum(dt_locaux_.size_reelle());
      Cout << "dt_loc stats: min=" << dt_loc_min
           << "  max=" << dt_loc_max
           << "  mean=" << dt_loc_sum / dt_loc_n
           << "  ratio_max/min=" << dt_loc_max / (dt_loc_min + 1.e-300) << finl;
    }
}
bool Schema_Euler_Implicite_Stationnaire::iterateTimeStep(bool& converged)
{
  bool ok = Schema_Euler_Implicite::iterateTimeStep(converged);
 
  // After the parent's iterateTimeStep(), test_stationnaire() has run with
  // divisor dt_ (the CFL-based pseudo-timestep), giving a wrong absolute criterion.
  // We recompute the stationarity check at every step with the relative criterion:
  //   ||U^{n+1} - U^n|| / ||U^n||  < seuil_statio
  //
  // Post-condition left by test_stationnaire():
  //   valeurs() = U^n   (restored by "present = passe")
  //   futur()   = U^{n+1}
  {
    set_stationnaire_atteint() = 1;
    residu_ = 0.;
    for (int i = 0; i < pb_base().nombre_d_equations(); i++)
      {
        Equation_base& eqn = pb_base().equation(i);
        const DoubleTab& un   = eqn.inconnue().valeurs();  // U^n
        const DoubleTab& unp1 = eqn.inconnue().futur();    // U^{n+1}
        const double norm_un = mp_norme_vect(un);
        DoubleTab tab_critere(unp1);
        tab_critere -= un;
        tab_critere /= std::max(norm_un, 1.e-10);
        update_critere_statio(tab_critere, eqn);
      }
  }
 
  return ok;
}
 
bool Schema_Euler_Implicite_Stationnaire::corriger_dt_calcule(double& dt_global) const
{
  // Replicate the RAM print from Schema_Temps_base::corriger_dt_calcule,
  if (limpr())
    {
      Cout << finl;
      imprimer_ram_totale();
    }
 
  // Set the global pseudo-timestep to 1 so that temps_courant_ = nb_pas_dt.
  // This makes the TIME axis in outputs and .dt_ev equal to the iteration count,
  // which is scheme- and mesh-independent and allows meaningful cross-case comparison.
  dt_global = 1.0;
  return 1;
}