Implicit_steady.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 <Implicit_steady.h>
#include <Schema_Euler_Implicite_Stationnaire.h>
#include <Domaine_VF.h>
#include <Domaine_Cl_VEF.h>
#include <Navier_Stokes_std.h>
#include <EChaine.h>
#include <Debog.h>
#include <Matrice_Bloc.h>
#include <Assembleur_base.h>
#include <Schema_Temps_base.h>
#include <TRUSTTrav.h>
#include <Fluide_Quasi_Compressible.h>
#include <Dirichlet.h>
#include <Probleme_base.h>
#include <assert.h>
 
Implemente_instanciable(Implicit_steady,"Implicit_steady",Simpler);
// XD implicit_steady implicite implicit_steady BRACE this is the implicit solver using a dual time step. Remark: this
// XD_CONT solver can be used only with the Implicit_Euler_Steady_Scheme time scheme.
 
Sortie& Implicit_steady::printOn(Sortie& os ) const
{
  return Simpler::printOn(os);
}
 
Entree& Implicit_steady::readOn(Entree& is )
{
  Simpler::readOn(is);
  return is;
}
 
void test_impose_bound_cond(Equation_base& eqn,DoubleTab& current2,const char * msg,int flag)
{
  return;
  DoubleTab& present = eqn.inconnue().futur();
  DoubleTab sauv(present);
  const Schema_Temps_base& sch = eqn.probleme().schema_temps();
  eqn.domaine_Cl_dis().imposer_cond_lim(eqn.inconnue(),sch.temps_courant()+sch.pas_de_temps());
  present -= sauv;
  double ecart_max=mp_max_abs_vect(present);
  Cout<<msg <<" "<<ecart_max<<finl;
 
  if ((ecart_max>1e-10))
    abort();
  present = sauv;
}
 
 
//This function slightly modifies the "Piso :: iterate_NS with avancement_crank_=0" function in order to take into account a variable time step
//Use of a variable time step with the IMPLICITE algorithm
//replaces the global time step with a local time step
void Implicit_steady::iterer_NS(Equation_base& eqn,DoubleTab& current,DoubleTab& pression, double dt,
                                Matrice_Morse& matrice,double seuil_resol,DoubleTrav& secmem,int nb_iter,int& converge, int& ok)
{
  Schema_Temps_base& sch = eqn.probleme().schema_temps();
  Parametre_implicite& param_eqn = get_and_set_parametre_implicite(eqn);
  SolveurSys& le_solveur_ = param_eqn.solveur();
  converge = 1;
  Navier_Stokes_std& eqnNS = ref_cast(Navier_Stokes_std,eqn);
  if(!sub_type(Schema_Euler_Implicite_Stationnaire,sch))
    {
      Cerr << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << finl;
      Cerr << "Implicit_steady solveur can be used only with the Implicit_Euler_Steady_Scheme or Schema_Euler_Implicite_Stationnaire scheme!" << finl;
      Cerr << "Please, contact TRUST support." << finl;
      Process::exit();
    }
 
  DoubleTrav gradP(current);
  DoubleTrav correction_en_pression(pression);
  DoubleTrav resu(current);
 
  Operateur_Grad& gradient = eqnNS.operateur_gradient();
  Operateur_Div& divergence = eqnNS.operateur_divergence();
 
  Schema_Euler_Implicite_Stationnaire& sch_steady =
    ref_cast(Schema_Euler_Implicite_Stationnaire, sch);
 
  // --- Update dt_loc_ and decide whether to reassemble the pressure matrix ---
  // dt_loc_ is compared against the freshly computed dt_locaux_.
  // If the max relative variation exceeds seuil_variation_dt_ (default 25%),
  // dt_loc_ is updated and the pressure matrix is reassembled.
  // Otherwise dt_loc_ is kept from the previous iteration.
  // All 3 steps use the same dt_loc_, ensuring mathematical consistency:
  //   Step 1 (predictor matrix M/dt_loc + K),
  //   Step 2 (pressure matrix B*(M/dt_loc)^-1*Bt),
  //   Step 3 (velocity correction U = U* - dt_loc * grad(P')).
  DoubleVect& dt_loc = sch_steady.get_dt_loc();
  bool need_reassembly = (dt_loc.size() == 0);  // first iteration
  double max_var = 0.;
  {
    const DoubleVect dt_loc_new = sch.pas_de_temps_locaux();
    if (!need_reassembly)
      {
        for (int i = 0; i < dt_loc_new.size(); i++)
          {
            double var = std::abs(dt_loc_new[i] - dt_loc[i]) / (dt_loc[i] + 1.e-30);
            max_var = std::max(max_var, var);
          }
        max_var = Process::mp_max(max_var);
        // Reassemble only if variation > threshold AND at a reassembly interval
        need_reassembly = (max_var > seuil_variation_dt_) && (sch.nb_pas_dt() % reassembly_interval_ == 0);
      }
    if (need_reassembly)
      dt_loc = dt_loc_new;
  }
 
  if (max_var == 0. && need_reassembly)
    Cout << "Implicit_steady: first iteration  -> dt_loc_ initialised, pressure matrix reassembled" << finl;
  else
    {
      if (sch.limpr())
        Cout << "Implicit_steady: dt_loc max relative variation: " << max_var
             << (need_reassembly ? "  -> dt_loc_ updated, pressure matrix reassembled"
                 : "  -> dt_loc_ unchanged, pressure matrix reused") << finl;
    }
 
  // Step 1: Predictor — assemble K = M/dt_loc + CONV + DIFF  and  RHS = M/dt_loc*Un - Bt*Pn + Sv
  // M/dt_loc is incorporated by assembler_avec_inertie via the virtual call
  // Schema_Euler_Implicite_Stationnaire::ajouter_inertie() which uses dt_loc_.
  gradient.calculer(pression,gradP);
  resu -= gradP;
 
  eqnNS.assembler_avec_inertie(matrice,current,resu);
  le_solveur_->reinit();
  Debog::verifier("Implicit_steady::iterer_NS resu apres assembler_avec_inertie",resu);
 
  Matrice& matrice_en_pression_2 = eqnNS.matrice_pression();
  SolveurSys& solveur_pression_ = eqnNS.solveur_pression();
 
  // Solve K*U* = RHS  ->  current = U*
  Cout << "Compute U* :" << finl;
  le_solveur_.resoudre_systeme(matrice,resu,current);
  test_impose_bound_cond(eqn,current,"apres resolution ",0);
  Debog::verifier("Implicit_steady::iterer_NS current apres CL",current);
 
  current.echange_espace_virtuel();
 
  // Compute RHS of pressure equation: secmem = B*U*
  if (eqn.probleme().is_dilatable())
    {
      Cerr<<"Steady option is not compatible with the quasi/weakly compressible models !"<<finl;
      Cerr << "Please, contact TRUST support." << finl;
      Process::exit();
    }
  else
    {
      divergence.calculer(current,secmem);
    }
  secmem *= -1;
  secmem.echange_espace_virtuel();
 
  // Step 2: Pressure solve — solve (B*(M/dt_loc)^-1*Bt)*P' = B*U*
  // correction_en_pression = P'
  Cout << "Implicit_steady: solving mass equation :" << finl;
 
  // Compute M/dt_loc using dt_loc_ (same as step 1)
  DoubleVect m_dt(dt_loc);
  calcul_mat_masse_diviser_par_dt(eqnNS, m_dt, dt_loc);
 
  // Reassemble B*(M/dt_loc)^-1*Bt only if dt_loc_ was updated above
  if (need_reassembly)
    {
      eqnNS.assembleur_pression()->assembler_mat(matrice_en_pression_2, m_dt, 1, 1);
      solveur_pression_->reinit();
    }
 
  // Solve (B*(M/dt_loc)^-1*Bt)*P' = B*U*
  solveur_pression_.resoudre_systeme(matrice_en_pression_2.valeur(),
                                     secmem,correction_en_pression);
  correction_en_pression.echange_espace_virtuel();
  // Compute M^-1*Bt*P' = grad(P')
  gradient->multvect(correction_en_pression,gradP);
  eqn.solv_masse().appliquer(gradP);
  // Step 3: Velocity correction — U^{n+1} = U* - dt_loc * M^{-1} * Bt * P'
  int size=gradP.size();
 
  // Expand dt_loc (face-sized) to velocity-sized: line_size()==1 in VDF, nb_dim in VEF.
  DoubleVect dt_loc_velocity(current);
  const int nb_dim = eqnNS.inconnue().valeurs().line_size();
  assert(dt_loc.size() * nb_dim == dt_loc_velocity.size());
  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();
 
  for(int i=0; i<size; i++)
    {
      gradP.addr()[i] *=dt_loc_velocity[i];  // dt_loc * grad(P')
    }
 
  gradP.echange_espace_virtuel();
  current -= gradP;
  test_impose_bound_cond(eqn,current,"apres resolution ",0);
  // current = U^{n+1}
  current.echange_espace_virtuel();
  divergence.calculer(current,secmem);
 
  // P^{n+1} = P^n + P'
  pression += correction_en_pression;
  Debog::verifier("Implicit_steady::iterer_NS pression fin", pression);
  eqnNS.assembleur_pression()->modifier_solution(pression);
  pression.echange_espace_virtuel();
  return;
}
// Compute M/dt_loc = rho * V_entrelace * porosite / dt_loc per face.
void Implicit_steady::calcul_mat_masse_diviser_par_dt(Navier_Stokes_std& eqnNS, DoubleVect& m_dt, const DoubleVect& dt_loc)
{
  const Domaine_VF& le_dom = ref_cast(Domaine_VF, eqnNS.domaine_dis());
  DoubleVect volumes(le_dom.volumes_entrelaces());
 
  // VEF: override boundary and non-standard interior faces with CL-corrected volumes
  // (same logic as Assembleur_P_VEF::assembler). VDF has no such correction.
  if (sub_type(Domaine_Cl_VEF, eqnNS.domaine_Cl_dis()))
    {
      const DoubleVect& volumes_cl = ref_cast(Domaine_Cl_VEF, eqnNS.domaine_Cl_dis()).volumes_entrelaces_Cl();
      for (int f = 0; f < volumes_cl.size(); f++)
        if (volumes_cl(f) != 0)
          volumes(f) = volumes_cl(f);
    }
  volumes.echange_espace_virtuel();
 
  const int size = volumes.size_totale();
  const DoubleVect& porosite = eqnNS.milieu().porosite_face();
  const bool has_porosity = (porosite.size_totale() == size);
  const DoubleVect& rho = eqnNS.fluide().masse_volumique().valeurs();
  const bool rho_per_face = (rho.size_totale() == size);
  for (int face = 0; face < size; face++)
    {
      const double rho_f = rho_per_face ? rho(face) : 1.0;
      const double por_f = has_porosity ? porosite(face) : 1.0;
      m_dt(face) = volumes(face) * por_f * rho_f / dt_loc(face);
    }
  m_dt.echange_espace_virtuel();
}