Piso.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 <Fluide_Dilatable_base.h>
#include <MD_Vector_composite.h>
#include <Discretisation_base.h>
#include <Navier_Stokes_IBM.h>
#include <Schema_Temps_base.h>
#include <MD_Vector_tools.h>
#include <Source_PDF_base.h>
#include <Assembleur_base.h>
#include <TRUSTTab_parts.h>
#include <Probleme_base.h>
#include <MD_Vector_std.h>
#include <Matrice_Bloc.h>
#include <TRUSTTrav.h>
#include <EChaine.h>
#include <Debog.h>
#include <Piso.h>
 
Implemente_instanciable(Piso,"Piso",Simpler);
// XD piso simpler piso INHERITS_BRACE Piso (Pressure Implicit with Split Operator) - method to solve N_S.
// XD attr seuil_convergence_implicite floattant seuil_convergence_implicite OPT Convergence criteria.
// XD attr nb_corrections_max entier nb_corrections_max OPT Maximum number of corrections performed by the PISO
// XD_CONT algorithm to achieve the projection of the velocity field. The algorithm may perform less corrections then
// XD_CONT nb_corrections_max if the accuracy of the projection is sufficient. (By default nb_corrections_max is set to
// XD_CONT 21).
 
Implemente_instanciable_sans_constructeur(Implicite,"Implicite",Piso);
// XD implicite piso implicite INHERITS_BRACE similar to PISO, but as it looks like a simplified solver, it will use
// XD_CONT fewer timesteps. But it may run faster because the pressure matrix is not re-assembled and thus provides CPU
// XD_CONT gains.
 
Implicite::Implicite()
{
  avancement_crank_ = 1;
}
 
Sortie& Piso::printOn(Sortie& os ) const { return Simpler::printOn(os); }
 
Entree& Piso::readOn(Entree& is ) { return Simpler::readOn(is); }
 
Sortie& Implicite::printOn(Sortie& os ) const { return Piso::printOn(os); }
 
Entree& Implicite::readOn(Entree& is ) { return Piso::readOn(is); }
 
 
Entree& Piso::lire(const Motcle& motlu,Entree& is)
{
  Motcles les_mots(2);
  {
    les_mots[0] = "nb_corrections_max";
    les_mots[1] = "with_sources";
  }
 
  int rang = les_mots.search(motlu);
  switch(rang)
    {
    case 0:
      {
        is >> nb_corrections_max_;
        if (nb_corrections_max_ < 2)
          {
            Cerr<<"There should be at least two corrections steps for the PISO algorithm."<<finl;
          }
        break;
      }
    case 1:
      {
        is >> with_sources_;
        break;
      }
    default :
      {
        return Simpler::lire(motlu, is);
      }
    }
  return is;
}
 
 
 
void test_imposer_cond_lim(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;
  // BM, je remplace max_abs par mp_pax_abs: du coup la methode doit etre appelee simultanement par tous les procs.
  double ecart_max=mp_max_abs_vect(present);
  Cout<<msg <<" "<<ecart_max<<finl;
 
  if ((ecart_max>1e-10))
    abort();
  present = sauv;
   */
}
 
//Entree Un ; Pn
//Sortie Un+1 = U***_k ; Pn+1 = P**_k
//n designe une etape temporelle
 
void Piso::iterer_NS(Equation_base& eqn,DoubleTab& current,DoubleTab& pression,
                     double dt,Matrice_Morse& matrice,double seuil_resol,DoubleTrav& secmem,int nb_ite,int& converge, int& ok)
{
  converge = 1;
  if (nb_ite>1) return;
  Navier_Stokes_std& eqnNS = ref_cast(Navier_Stokes_std,eqn);
  eqnNS.reassembler_pression_si_necessaire();
  const bool is_NS_IBM = sub_type(Navier_Stokes_IBM, eqnNS);
 
  eqnNS.setPressureTimeN(); //sometimes we need a second special treatement like for ALE for example
 
  if (eqnNS.discretisation().que_suis_je() == "PolyMAC_CDO")
    return iterer_NS_PolyMAC_CDO(eqnNS, current, pression, dt, matrice, ok);
 
  Parametre_implicite& param_eqn = get_and_set_parametre_implicite(eqn);
  SolveurSys& le_solveur_ = param_eqn.solveur();
 
  DoubleTrav gradP(current);
  DoubleTrav correction_en_pression(pression);
  DoubleTrav resu(current);
 
 
  int is_dilat = eqn.probleme().is_dilatable();
 
  double vitesse_norme,vitesse_norme_old ;
  double pression_norme,pression_norme_old ;
  vitesse_norme_old = 1.e10 ;
  pression_norme_old = 1.e10 ;
 
  // int deux_entrees = 0;
  //if (current.nb_dim()==2) deux_entrees = 1;
  Operateur_Grad& gradient = eqnNS.operateur_gradient();
  Operateur_Div& divergence = eqnNS.operateur_divergence();
 
  //int nb_comp = 1;
  //int nb_dim = current.nb_dim();
 
  //  if(nb_dim==2)    nb_comp = current.dimension(1);
 
  //Construction de matrice et resu
  //matrice = A[Un] = M/delta_t + CONV +DIFF
  //resu =  A[Un]Un -(A[Un]Un-Ss) + Sv -BtPn
 
  if (is_NS_IBM)
    add_penality_term(eqnNS, resu, gradP);
 
  gradient->calculer(pression,gradP);
  if (eqnNS.has_interface_blocs()) //si l'interface blocs est disponible, on l'utilise
    {
      eqnNS.assembler_blocs_avec_inertie({{ "vitesse", &matrice }}, resu);
      if (eqnNS.discretisation().is_poly_family())
        matrice.ajouter_multvect(current, resu);  //pour ne pas etre en increment
    }
  else //sinon, on passe par ajouter/contribuer
    {
      resu -= gradP;
      eqnNS.assembler_avec_inertie(matrice,current,resu);
    }
 
  le_solveur_->reinit();
 
  //Construction de matrice_en_pression_2 = BD-1Bt[Un]
  //Assemblage reeffectue seulement pour algorithme Piso (avancement_crank_==0)
  Matrice& matrice_en_pression_2 = eqnNS.matrice_pression();
  SolveurSys& solveur_pression_ = eqnNS.solveur_pression();
 
  if (avancement_crank_==0 && !with_sources_)
    {
      assembler_matrice_pression_implicite(eqnNS,matrice,matrice_en_pression_2);
      solveur_pression_->reinit();
    }
 
  //Etape predicteur
  //Resolution du systeme A[Un]U* = -BtPn + Sv + Ss
  //current = U*
  le_solveur_.resoudre_systeme(matrice,resu,current);
 
  test_imposer_cond_lim(eqn,current,"apres resolution ",0);
  current.echange_espace_virtuel();
  Debog::verifier("Piso::iterer_NS current apres solveur",current);
 
  if (with_sources_)
    {
      matrice.get_set_coeff() = 0;
      eqnNS.sources().contribuer_a_avec(current, matrice);
      eqnNS.solv_masse().ajouter_masse(dt, matrice, 1);
      assembler_matrice_pression_implicite(eqnNS,matrice,matrice_en_pression_2);
      solveur_pression_->reinit();
    }
 
  //Calcul de secmem = BU* (en incompressible) BU* -drho/dt (en quasi-compressible)
  if (is_dilat)
    {
      if (with_d_rho_dt_)
        {
          Fluide_Dilatable_base& fluide_dil = ref_cast(Fluide_Dilatable_base,eqn.milieu());
          fluide_dil.secmembre_divU_Z(secmem);
          secmem *= -1;
        }
      else   secmem = 0;
      divergence.ajouter(current,secmem);
    }
  else
    divergence.calculer(current,secmem);
  secmem *= -1;
 
  if (is_NS_IBM)
    correct_incr_pressure(eqnNS, secmem);
 
  secmem.echange_espace_virtuel();
  Debog::verifier("Piso::iterer_NS secmem",secmem);
 
  // GF il ne faut pas modifier le scd membre le terme en du/dt au bord a deja ete pris en compte dans la resolution precedente
  //  eqnNS.assembleur_pression()->modifier_secmem(secmem);
 
  //Etape de correction 1
  Cout << "Solving mass equation :" << finl;
  //Description du cas implicite
  //Resolution du systeme (BD-1Bt)P' = Bu* (D-1 = M-1 pour le cas implicite)
  //correction_en_pression = P' pour Piso et correction_en_pression = delta_t*P' pour implicite
  eqnNS.assembleur_pression()->modifier_secmem_pour_incr_p(pression, 1. / dt, secmem);
  //si la matrice varie, passage increments -> valeurs pour aider les solveurs iteratifs
  if (with_sources_)
    matrice_en_pression_2->ajouter_multvect(pression, secmem);
  solveur_pression_.resoudre_systeme(matrice_en_pression_2.valeur(),
                                     secmem,correction_en_pression);
  if (with_sources_)
    correction_en_pression -= pression;
  correction_en_pression.echange_espace_virtuel();
  Debog::verifier("Piso::iterer_NS apres correction_pression",correction_en_pression);
 
  if (avancement_crank_==1)
    {
      //calcul de Un+1
      //Calcul de Bt(delta_t*delta_P)
      gradient->multvect(correction_en_pression,gradP);
 
      if (with_sources_)
        for (int i = 0, N = gradP.nb_dim() == 2 ? gradP.dimension(1) : 1; i < gradP.dimension_tot(0); i++)
          for (int n  = 0; n < N; n++)
            gradP(i, n) = matrice.get_tab1()(N * i + n + 1) > matrice.get_tab1()(N * i + n) && matrice(N * i + n, N * i + n) ? gradP(i, n) / matrice(N * i + n, N * i + n) : 0;
      else
        eqn.solv_masse().appliquer(gradP);
 
      if (is_NS_IBM)
        correct_gradP(eqnNS, gradP);
 
      //Calcul de Un+1 = U* -delta_t*delta_P
      current -= gradP;
      eqn.solv_masse().corriger_solution(current, current); //pour PolyMAC_MPFA : sert a corriger ve
      current.echange_espace_virtuel();
      divergence.calculer(current,secmem);
 
      //Calcul de Pn+1 = Pn + (delta_t*delta_P)/delat_t
      Debog::verifier("Piso::iterer_NS correction avant dt",correction_en_pression);
      if (!with_sources_)
        correction_en_pression /= dt;
 
      if (is_NS_IBM)
        correct_pressure(eqnNS, pression, correction_en_pression);
      else
        pression += correction_en_pression;
      // </IBM>
      eqnNS.assembleur_pression()->modifier_solution(pression);
      pression.echange_espace_virtuel();
 
      Debog::verifier("Piso::iterer_NS pression finale",pression);
      Debog::verifier("Piso::iterer_NS current final",current);
      if (is_dilat)
        {
          // on redivise par rho_np_1 avant de sortir
          diviser_par_rho_np1_face(eqn,current);
        }
 
      eqnNS.updateFluidForce(current);
 
      return;
    }
 
  // calcul de la correction en vitesse premiere etape (DU' =-Bt P)
 
  //Calcul de P* = Pn + P'
  pression += correction_en_pression;
  eqnNS.assembleur_pression()->modifier_solution(pression);
  //Resolution du systeme D[Un]U' = -BtP'
  DoubleTrav correction_en_vitesse(current);
  calculer_correction_en_vitesse(correction_en_pression,gradP,correction_en_vitesse,matrice,gradient);
 
  //Calcul de U** = U* + U'
  current += correction_en_vitesse;
  test_imposer_cond_lim(eqn,current,"apres premiere correction ",0);
  Debog::verifier("Piso::iterer_NS arpes cor pression",pression);
  Debog::verifier("Piso::iterer_NS arpes cor vitesse",current);
 
  //Etape correcteur 2
  for (int compt=0; compt<nb_corrections_max_-1; compt++)
    {
      correction_en_vitesse.echange_espace_virtuel();
      //Resolution du systeme D resu = EU' + (resu2=0) pour stocker resu = D-1EU'
      DoubleTrav resu2(resu);
      int status = inverser_par_diagonale(matrice,resu2,correction_en_vitesse,resu);
 
      if (status!=0) exit();
      // calcul de P''  BD-1Bt P''= -div(D-1EU')
 
      resu.echange_espace_virtuel();
      //Calcul de B(D-1EU')
      divergence.calculer(resu,secmem);
      secmem *= -1;
      secmem.echange_espace_virtuel();
 
      //Resolution du systeme (BD-1Bt)P'' = (BD-1E)U'
      //correction_en_pression = P''
      correction_en_pression = 0;
      solveur_pression_.resoudre_systeme(matrice_en_pression_2.valeur(),
                                         secmem,correction_en_pression);
 
#ifdef DEBUG
      // Pour verifier que le calcul du gradient ne modifie pas la pression
      DoubleTrav correction_en_pression_mod(pression);
      correction_en_pression_mod = correction_en_pression;
#endif
      //Resolution du systeme D[Un]U'' = -BtP''
      //correction_en_vitesse = U''
      calculer_correction_en_vitesse(correction_en_pression,gradP,correction_en_vitesse,matrice,gradient);
 
#ifdef DEBUG
      correction_en_pression_mod -= correction_en_pression;
      assert(max_abs(correction_en_pression_mod)==0);
#endif
 
      //Calcul de U'' = U'' + D-1EU'
      correction_en_vitesse += resu;
      // ajout des increments
 
      // Optimization: combine 2 mp_norme_vect into 1 collective call
      double vitesse_carre = local_carre_norme_vect(correction_en_vitesse);
      double pression_carre = local_carre_norme_vect(correction_en_pression);
      mp_sum_for_each(vitesse_carre, pression_carre);
      vitesse_norme = sqrt(vitesse_carre);
      pression_norme = sqrt(pression_carre);
 
      if ( (vitesse_norme>vitesse_norme_old) || (pression_norme>pression_norme_old) )
        {
          Cout <<"PISO : "<< compt+1 <<" corrections to perform the projection."<< finl;
          if (is_dilat)
            {
              // on redivise par rho_np_1 avant de sortir
              diviser_par_rho_np1_face(eqn,current);
            }
          return ;
        }
 
      vitesse_norme_old = vitesse_norme;
      pression_norme_old = pression_norme;
      //Calcul de P** = P* + P''
      pression += correction_en_pression;
      eqnNS.assembleur_pression()->modifier_solution(pression);
 
      //Calcul de U*** = U** + U''
      current += correction_en_vitesse;
      test_imposer_cond_lim(eqn,current,"apres correction (int)__LINE__",0);
 
      Debog::verifier("Piso::iterer_NS apres correction pression",pression);
      Debog::verifier("Piso::iterer_NS apres correction vitesse",current);
    }
  if (is_dilat)
    {
      diviser_par_rho_np1_face(eqn,current);
      //ref_cast(Navier_Stokes_QC,eqn).impr_impl(eqnNS, Cout);
    }
  current.echange_espace_virtuel();
  // divergence.calculer(current, secmem); Cerr<<" ici DIVU  "<<mp_max_abs_vect(secmem)<<finl;;
  Cout <<"PISO : "<<nb_corrections_max_<<" corrections to perform the projection."<< finl;
}
 
//version PolyMAC_CDO de la fonction ci-dessus
void Piso::iterer_NS_PolyMAC_CDO(Navier_Stokes_std& eqn, DoubleTab& current, DoubleTab& pression, double dt, Matrice_Morse& matrice, int& ok)
{
  if (avancement_crank_ == 0)
    Process::exit("sorry, the PISO solver is not implemented with PolyMAC_CDO! Please use Implicite instead");
  Parametre_implicite& param_eqn = get_and_set_parametre_implicite(eqn);
  SolveurSys& le_solveur_ = param_eqn.solveur();
 
  Navier_Stokes_std& eqnNS = ref_cast(Navier_Stokes_std, eqn);
  Operateur_Grad& op_grad = eqnNS.operateur_gradient();
  Operateur_Div& op_div = eqnNS.operateur_divergence();
 
  DoubleTrav dv(current); //, dP(pression); //corrections en vitesse / pression
 
  /* etape de prediction : current <- v(0) ne verifiant pas div = 0 */
  DoubleTrav secmem_NS(current), v_new(current);
  op_grad.ajouter(pression, secmem_NS);
  secmem_NS *= -1;
  eqnNS.assembler_avec_inertie(matrice, current, secmem_NS);
  le_solveur_->reinit();
  le_solveur_.resoudre_systeme(matrice, secmem_NS, v_new);
  v_new.echange_espace_virtuel();
 
  Matrice& mat_press_orig = eqn.matrice_pression(), mat_press;
 
  /* etapes de correction : current <- v(i) verifiant div = 0 mais approche pour NS, pression <- p(i) correspondant a v(i) */
  for (int i = 0; i < 1; i++)
    {
      DoubleTrav sol_M(pression), secmem_M(pression);
      /* resolution en (dt * dp(i), dv(i)) */
      //second membre : divergence, NS
      DoubleTab_parts p_sec(secmem_M), p_sol(sol_M); //p_sec/sol[0] -> elements, p_sec/sol[1] -> faces
      //bloc superieur : div v(i-1)
      op_div.ajouter(v_new, p_sec[0]);
      //bloc inferieur : residu de l'etape de prediction
      p_sec[1] = 0;      //p_res[0];
 
      //matrice (sauf si avancement_crank_ == 1) : prise en compte des contributions des sources (et pas des operateurs!)
      if (avancement_crank_==0 || with_sources_)
        {
          matrice.get_set_coeff() = 0, eqn.sources().contribuer_a_avec(current, matrice);
          DoubleTrav diag(p_sec[1]);
          for (int j = 0; j < p_sec[1].dimension(0); j++)
            diag(j) = dt * matrice(j, j);
          diag.echange_espace_virtuel();
          eqn.assembleur_pression()->assembler_mat(mat_press, diag, 1, 1);
          eqn.solveur_pression()->reinit();
        }
 
      //resolution
      Cerr << "PISO : |sec dp| < " << mp_max_abs_vect(p_sec[0]) << " |sec dv| < " << mp_max_abs_vect(p_sec[1]) << finl;
      eqn.solveur_pression().resoudre_systeme(avancement_crank_ == 1 ? mat_press_orig.valeur() : mat_press.valeur(), secmem_M, sol_M);
      //mises a jour : v^(i) = v^(i-1)+dv^(i), p^(i) = p^(i-1) + dp^(i)
      eqn.assembleur_pression()->corriger_vitesses(sol_M, dv);
      Cerr << "PISO : |dp| < " << mp_max_abs_vect(p_sol[0]) / dt << " |dv| < " << mp_max_abs_vect(dv);
      v_new += dv, sol_M /= dt, pression -= sol_M;
      //
      p_sec[0] = 0, op_div.ajouter(v_new, p_sec[0]);
      Cerr << " |div v| < " << mp_max_abs_vect(p_sec[0]) << finl;
      pression.echange_espace_virtuel();
    }
  current = v_new;
}
 
void Piso::add_penality_term(Navier_Stokes_std& eqnNS, DoubleTrav& resu , DoubleTrav& gradP)
{
  // <IBM> Taking into account penality term for Immersed Boundary Method
  Navier_Stokes_IBM& eqnNS_IBM = ref_cast(Navier_Stokes_IBM, eqnNS);
  const int i_source_PDF = eqnNS_IBM.get_i_source_pdf();
  if (i_source_PDF != -1)
    {
      Source_PDF_base& src = dynamic_cast<Source_PDF_base&>((eqnNS.sources())[i_source_PDF].valeur());
      DoubleTab secmem_pdf(resu);
      src.calculer_pdf(secmem_pdf);
 
      // Terme en temps : -rho/delta_t ksi_gamma Un
      int pdf_bilan = src.get_modele().pdf_bilan();
      if (pdf_bilan == 1)
        {
          int i_traitement_special = 101;
          if (eqnNS.nombre_d_operateurs() > 1)
            {
              if (eqnNS.vitesse_pour_transport().le_nom() == "rho_u")
                i_traitement_special = 1;
            }
          DoubleTrav secmem_pdf_time(resu);
          src.calculer(secmem_pdf_time, i_traitement_special);
          secmem_pdf += secmem_pdf_time;
        }
 
      // Sauvegarde de secmem_pdf
      secmem_pdf.echange_espace_virtuel();
      src.set_sec_mem_pdf(secmem_pdf);
 
      const DoubleTab& coeff = eqnNS_IBM.get_champ_coeff_pdf_som();
      if (eqnNS_IBM.get_gradient_pression_qdm_modifie() == 1)
        {
          Cerr << "(IBM) Immersed Interface: modified pressure gradient in momentum equation." << finl;
          gradP /= coeff;
        }
      gradP.echange_espace_virtuel();
    }
}
 
void Piso::correct_gradP(Navier_Stokes_std& eqnNS, DoubleTrav& gradP)
{
  Navier_Stokes_IBM& eqnNS_IBM = ref_cast(Navier_Stokes_IBM, eqnNS);
  const int i_source_PDF = eqnNS_IBM.get_i_source_pdf();
  if ((i_source_PDF != -1) && (eqnNS_IBM.get_correction_vitesse_modifie() == 1))
    {
      Cerr << "(IBM) Immersed Interface: modified velocity correction." << finl;
      const DoubleTab& coeff = eqnNS_IBM.get_champ_coeff_pdf_som();
      gradP /= coeff;
      gradP.echange_espace_virtuel();
    }
}
 
void Piso::correct_incr_pressure(Navier_Stokes_std& eqnNS, DoubleTrav& secmem)
{
  // <IBM> Taking into account zero velocity divergence for Immersed Boundary Method
  Navier_Stokes_IBM& eqnNS_IBM = ref_cast(Navier_Stokes_IBM, eqnNS);
  const int i_source_PDF = eqnNS_IBM.get_i_source_pdf();
  if ((i_source_PDF != -1) && (eqnNS_IBM.get_correction_matrice_pression() == 1))
    {
      Cerr << "(IBM) Immersed Interface: velocity divergence is zero for crossed elements." << finl;
      const DoubleTab& coeff = eqnNS_IBM.get_champ_coeff_pdf_som();
      Source_PDF_base& src = dynamic_cast<Source_PDF_base&>((eqnNS.sources())[i_source_PDF].valeur());
      src.correct_incr_pressure(coeff, secmem);
    }
}
 
void Piso::correct_pressure(Navier_Stokes_std& eqnNS, DoubleTab& pression, DoubleTab& correction_en_pression)
{
  // <IBM> Immersed Interface: modified pressure correction.
  Navier_Stokes_IBM& eqnNS_IBM = ref_cast(Navier_Stokes_IBM, eqnNS);
  const int i_source_PDF = eqnNS_IBM.get_i_source_pdf();
  if ((i_source_PDF != -1) && (eqnNS_IBM.get_correction_pression_modifie()==1))
    {
      Cerr<<"Immersed Interface: modified pressure correction."<<finl;
      const DoubleTab& coeff = eqnNS_IBM.get_champ_coeff_pdf_som();
      const Source_PDF_base& src = dynamic_cast<Source_PDF_base&>((eqnNS.sources())[i_source_PDF].valeur());
      src.correct_pressure(coeff,pression,correction_en_pression);
    }
  else
    pression += correction_en_pression;
}