Navier_Stokes_Variable_Density.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_Variable_Density.h>
#include <Fluide_Incompressible.h>
#include <Probleme_base.h>
#include <Schema_Implicite_base.h>
#include <Debog.h>
 
Implemente_instanciable( Navier_Stokes_Variable_Density, "Navier_Stokes_Variable_Density", Navier_Stokes_std ) ;
 
Sortie& Navier_Stokes_Variable_Density::printOn( Sortie& os ) const
{
  Navier_Stokes_std::printOn( os );
  return os;
}
 
Entree& Navier_Stokes_Variable_Density::readOn( Entree& is )
{
  Navier_Stokes_std::readOn( is );
  const Navier_Stokes_Variable_Density::Interface_Milieu& milieu = dynamic_cast<const Navier_Stokes_Variable_Density::Interface_Milieu&>(le_fluide.valeur());
  terme_diffusif->associer_champ_masse_volumique(milieu.masse_volumique_variable());
  return is;
}
 
void Navier_Stokes_Variable_Density::associer_sch_tps_base(const Schema_Temps_base& un_schema_en_temps)
{
  if (sub_type(Schema_Implicite_base, un_schema_en_temps))
    {
      Cerr << "implicit scheme not available (yet) with Navier_Stokes_Variable_Density variant of the Navier-Stokes equations" << finl;
      Cerr << "you should use an explicit scheme and, if you wish, activate the implicit treatment of the diffusion operator in the velocity prediction step" << finl;
      exit();
    }
  Equation_base::associer_sch_tps_base(un_schema_en_temps);
}
 
const Champ_Don_base& Navier_Stokes_Variable_Density::diffusivite_pour_transport() const
{
  return le_fluide->viscosite_dynamique();
}
 
void Navier_Stokes_Variable_Density::discretiser()
{
  Navier_Stokes_std::discretiser();
  const Domaine_dis_base& domaine_dis= le_dom_dis.valeur();
  rho_faces_.typer("Champ_Fonc_Face");
  probleme().discretisation().nommer_completer_champ_physique(domaine_dis, "masse_volumique_faces","kg/m^3", rho_faces_, probleme());
  champs_compris_.ajoute_champ(rho_faces_);
  rho_faces_->valeurs() = inconnue().valeurs();
  assembleur_pression_->set_resoudre_increment_pression(1);
  assembleur_pression_->set_resoudre_en_u(0);
}
 
int Navier_Stokes_Variable_Density::preparer_calcul()
{
  const double temps = schema_temps().temps_courant();
  sources().mettre_a_jour(temps);
  Equation_base::preparer_calcul();
 
  // begin modification of Navier_Stokes_std::preparer_calcul()
  const Navier_Stokes_Variable_Density::Interface_Milieu& milieu = dynamic_cast<const Navier_Stokes_Variable_Density::Interface_Milieu&>(le_fluide.valeur());
  const DoubleTab& tab_rho = milieu.masse_volumique_variable().valeurs();
  DoubleTab& tab_rho_faces = rho_faces_->valeurs();
  rho_aux_faces(tab_rho, tab_rho_faces);
  tab_rho_faces.echange_espace_virtuel();
  assembleur_pression_->assembler_rho_variable(matrice_pression_, rho_faces_);
  // end modification of Navier_Stokes_std::preparer_calcul()
 
  Debog::verifier("Navier_Stokes_std::preparer_calcul, la_pression av projeter", la_pression->valeurs());
  if (projection_a_faire())
    projeter();
 
  // Au cas ou une cl de pression depend de u que l'on vient de modifier
  le_dom_Cl_dis->mettre_a_jour(temps);
  Debog::verifier("Navier_Stokes_std::preparer_calcul, la_pression ap projeter", la_pression->valeurs());
 
  // Initialisation du champ de pression (resolution de Laplacien(P)=0 avec les conditions limites en pression)
  // Permet de demarrer la resolution avec une bonne approximation de la pression (important pour le Piso ou P!=0)
  if (!probleme().reprise_effectuee() && methode_calcul_pression_initiale_ != 3)
    {
      Cout << "Estimation du champ de pression au demarrage:" << finl;
      DoubleTrav secmem(la_pression->valeurs());
      DoubleTrav vpoint(gradient_P->valeurs());
      gradient.calculer(la_pression->valeurs(), gradient_P->valeurs());
      // begin modification of Navier_Stokes_std::preparer_calcul()
      gradient_P->valeurs() /= tab_rho_faces;
      // end modification of Navier_Stokes_std::preparer_calcul()
      vpoint -= gradient_P->valeurs();
 
      if (methode_calcul_pression_initiale_ >= 2)
        for (int op = 0; op < nombre_d_operateurs(); op++)
          operateur(op).ajouter(vpoint);
      if (methode_calcul_pression_initiale_ >= 1)
        {
          int mod = 0;
          if (le_schema_en_temps->pas_de_temps() == 0)
            {
              double dt = std::max(le_schema_en_temps->pas_temps_min(), calculer_pas_de_temps());
              dt = std::min(dt, le_schema_en_temps->pas_temps_max());
              le_schema_en_temps->set_dt() = (dt);
              mod = 1;
            }
          sources().ajouter(vpoint);
          if (mod)
            le_schema_en_temps->set_dt() = 0;
        }
 
      solveur_masse->appliquer(vpoint);
      vpoint.echange_espace_virtuel();
      divergence.calculer(vpoint, secmem);
      secmem *= -1;
      secmem.echange_espace_virtuel();
 
      assembleur_pression_->modifier_secmem_pour_incr_p(la_pression->valeurs(), 1, secmem);
      DoubleTrav inc_pre(la_pression->valeurs());
      solveur_pression_.resoudre_systeme(matrice_pression_.valeur(), secmem, inc_pre);
      Cerr << "Pressure increment computed successfully" << finl;
 
      // On veut que l'espace virtuel soit a jour, donc all_items
      operator_add(la_pression->valeurs(), inc_pre, VECT_ALL_ITEMS);
    }
  // Mise a jour pression
  la_pression->changer_temps(temps);
  calculer_la_pression_en_pa();
  // Calcul des forces de pression:
  gradient->calculer_flux_bords();
 
  // Calcul gradient_P (ToDo rendre coherent avec ::mettre_a_jour()):
  gradient.calculer(la_pression->valeurs(), gradient_P->valeurs());
  // begin modification of Navier_Stokes_std::preparer_calcul()
  gradient_P->valeurs() /= tab_rho_faces;
  // end modification of Navier_Stokes_std::preparer_calcul()
  gradient_P->changer_temps(temps);
 
  // Calcul divergence_U
  divergence.calculer(la_vitesse->valeurs(), divergence_U->valeurs());
  divergence_U->changer_temps(temps);
 
  if (le_traitement_particulier)
    le_traitement_particulier->preparer_calcul_particulier();
 
  Debog::verifier("Navier_Stokes_std::preparer_calcul, vitesse", inconnue());
  Debog::verifier("Navier_Stokes_std::preparer_calcul, pression", la_pression);
 
  return 1;
}
 
DoubleTab& Navier_Stokes_Variable_Density::derivee_en_temps_inco(DoubleTab& vpoint)
{
  const double dt = le_schema_en_temps->pas_de_temps();
  const Navier_Stokes_Variable_Density::Interface_Milieu& milieu = dynamic_cast<const Navier_Stokes_Variable_Density::Interface_Milieu&>(le_fluide.valeur());
  const DoubleTab& tab_rho = milieu.masse_volumique_variable().valeurs();
  DoubleTab& tab_rho_faces = rho_faces_->valeurs();
  rho_aux_faces(tab_rho, tab_rho_faces);
  tab_rho_faces.echange_espace_virtuel();
  // matrix associated with - V ∇ ⋅ ( 1 /  𝜌⁽ⁿ⁺¹⁾ ∇[] ) where V is the diagonal volume operator (solveur_masse)
  assembleur_pression_->assembler_rho_variable(matrice_pression_, rho_faces_);
  solveur_pression_->reinit();
 
  // ------------------------------
  // step 1: velocity v⋆ prediction
  // ------------------------------
  // vpoint temporary container is initialized
  vpoint = 0.;
  // vpoint -= V ( ∇ P⁽ⁿ⁾ / 𝜌⁽ⁿ⁺¹⁾ )
  gradient.calculer(la_pression->valeurs(), gradient_P->valeurs());
  gradient_P->valeurs() /= tab_rho_faces;
  vpoint -= gradient_P->valeurs();
  // vpoint += V ( F⁽ⁿ⁺¹⁾ / 𝜌⁽ⁿ⁺¹⁾ ) (les_ssources are already divided by rho)
  les_sources.ajouter(vpoint);
  // vpoint -= V ( v⁽ⁿ⁾ ⋅ ∇ v⁽ⁿ⁾ ) (terme_convectif considered on the rhs, so minus sign alread accounted for)
  DoubleTrav conv(vpoint);
  if (schema_temps().diffusion_implicite())
    {
      derivee_en_temps_conv(conv, la_vitesse->valeurs());
    }
  else
    {
      terme_convectif.calculer(la_vitesse->valeurs(), conv);
    }
  vpoint += conv;
 
  if (calculate_time_derivative())
    derivee_en_temps().valeurs() = vpoint;
  schema_temps().modifier_second_membre((*this), vpoint);
 
  if (schema_temps().diffusion_implicite())
    {
      // secmem is obtained from vpoint -> the right hand side of the equation to be solved for impliciting the diffusion term
      DoubleTrav secmem(vpoint);
      secmem = vpoint;
      solveur_masse->appliquer(secmem);
      // secmem = (- ∇ P⁽ⁿ⁾ + F⁽ⁿ⁺¹⁾ - 𝜌⁽ⁿ⁺¹⁾(v⁽ⁿ⁾ ⋅ ∇ v⁽ⁿ⁾))
      secmem *= tab_rho_faces;
      secmem.echange_espace_virtuel();
      // vpoint temporary container is set to v⁽ⁿ⁾
      vpoint = la_vitesse->valeurs();
      vpoint.echange_espace_virtuel();
      // solve (𝜌⁽ⁿ⁺¹⁾/ 𝛿t v⋆ + ∇ ⋅ 𝜏(v⋆)) = secmem
      // vpoint in input = v⁽ⁿ⁾
      // vpoint in output = (v⋆ - v⁽ⁿ⁾) / 𝛿t
      Equation_base::Gradient_conjugue_diff_impl(secmem, vpoint, tab_rho_faces.dimension(0), tab_rho_faces);
    }
  else
    {
      // vpoint += V ( ∇ ⋅ 𝜏(v⁽ⁿ⁾) ) / 𝜌⁽ⁿ⁺¹⁾
      DoubleTrav diff(vpoint);
      terme_diffusif.calculer(la_vitesse->valeurs(), diff);
      diff /= tab_rho_faces;
      vpoint += diff;
      solveur_masse->appliquer(vpoint);
      // vpoint = (v⋆ - v⁽ⁿ⁾) / 𝛿t
    }
  vpoint.echange_espace_virtuel();
 
  // -----------------------------------------------------------------------------
  // step 2: projection for calculating the pressure increment P⁽ⁿ⁺¹⁾ = P⁽ⁿ⁾ + 𝛿 P
  // -----------------------------------------------------------------------------
  // secmem temporary container = V ∇ ⋅ v⁽ⁿ⁾ / 𝛿t
  DoubleTrav secmem(la_pression->valeurs());
  divergence.calculer(la_vitesse->valeurs(), secmem);
  secmem /= dt;
  // secmem += V ∇ ⋅ (v⋆ - v⁽ⁿ⁾) / 𝛿t
  divergence.ajouter(vpoint, secmem);
  secmem *= -1;
  assembleur_pression_->modifier_secmem(secmem);
  secmem.echange_espace_virtuel();
 
  // solution of the Poisson problem - V ∇ ⋅ ( 1 /  𝜌⁽ⁿ⁺¹⁾ ∇ 𝛿 P ) = V ∇ ⋅ v⁽ⁿ⁾ / 𝛿t + V ∇ ⋅ (v⋆ - v⁽ⁿ⁾) / 𝛿t
  DoubleTrav inc_pre(la_pression->valeurs());
  solveur_pression_.resoudre_systeme(matrice_pression_.valeur(), secmem, inc_pre);
  // update pressure P⁽ⁿ⁺¹⁾ = P⁽ⁿ⁾ + 𝛿 P
  la_pression->valeurs() += inc_pre;
  assembleur_pression_->modifier_solution(la_pression->valeurs());
  la_pression->valeurs().echange_espace_virtuel();
 
  // --------------------------------------------------------------
  // step 3: velocity correction v⁽ⁿ⁺¹⁾ = v⋆ − 𝛿 t/𝜌⁽ⁿ⁺¹⁾ (∇ ⋅ 𝛿 P)
  // --------------------------------------------------------------
  // at this stage vpoint is (v⋆ - v⁽ⁿ⁾) / 𝛿t
  //
  // vpoint += ∇ ⋅ P⁽ⁿ⁾ / 𝜌⁽ⁿ⁺¹⁾
  solveur_masse->appliquer(gradient_P->valeurs());
  vpoint += gradient_P->valeurs();
  // calculate V ( ∇ ⋅ P⁽ⁿ⁺¹⁾ / 𝜌⁽ⁿ⁺¹⁾ )
  gradient.calculer(la_pression->valeurs(), gradient_P->valeurs());
  gradient_P->valeurs() /= tab_rho_faces;
  gradient_P->valeurs().echange_espace_virtuel();
  // vpoint -= ∇ ⋅ P⁽ⁿ⁺¹⁾ / 𝜌⁽ⁿ⁺¹⁾
  DoubleTrav gradP(gradient_P->valeurs());
  gradP = gradient_P->valeurs();
  solveur_masse->appliquer(gradP);
  vpoint -= gradP;
  vpoint.echange_espace_virtuel();
  // at this stage vpoint is (v⁽ⁿ⁺¹⁾ - v⁽ⁿ⁾) / 𝛿t
  return vpoint;
}
 
void Navier_Stokes_Variable_Density::calculer_la_pression_en_pa()
{
  la_pression_en_pa->valeurs() = la_pression->valeurs();
}
 
//DoubleTab& Navier_Stokes_Variable_Density::corriger_derivee_impl(DoubleTab& derivee)
//{
//  solveur_masse->set_name_of_coefficient_temporel(rho_faces_->le_nom());
//  Navier_Stokes_std::corriger_derivee_impl(derivee);
//  solveur_masse->set_name_of_coefficient_temporel("no_coeff");
//  return derivee;
//}
 
 
void Navier_Stokes_Variable_Density::modify_initial_gradP(DoubleTrav& gradP)
{
  DoubleTab& tab_rho_faces = rho_faces_->valeurs();
  operator_divide(gradP, tab_rho_faces, VECT_ALL_ITEMS);
}
 
void Navier_Stokes_Variable_Density::rho_aux_faces(const DoubleTab& tab_rho, DoubleTab& tab_rho_face)
{
  const Domaine_VF& domaineVF = ref_cast(Domaine_VF, domaine_dis());
  const IntTab& face_voisins = domaineVF.face_voisins();
  const DoubleVect& volumes = domaineVF.volumes();
  int nbfaces = domaineVF.nb_faces();
  const int nbfaces_bord = domaineVF.premiere_face_int();
  int el0, el1;
  double vol0, vol1;
  //
  for (int fac = 0; fac < nbfaces_bord; fac++)
    {
      el0 = face_voisins(fac, 0);
      if (el0 != -1)
        {
          tab_rho_face(fac) = tab_rho(el0);
        }
      else
        {
          el1 = face_voisins(fac, 1);
          tab_rho_face(fac) = tab_rho(el1);
        }
    }
 
  for (int fac = nbfaces_bord; fac < nbfaces; fac++)
    {
      el0 = face_voisins(fac, 0);
      el1 = face_voisins(fac, 1);
      vol0 = volumes(el0);
      vol1 = volumes(el1);
      tab_rho_face(fac) = (vol0 * tab_rho(el0) + vol1 * tab_rho(el1)) / (vol0 + vol1);
    }
}