EOS_Tools_VEF.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_Fluide_Dilatable_base.h>
#include <Source_Masse_Fluide_Dilatable_base.h>
#include <Loi_Etat_Melange_GP_base.h>
#include <Fluide_Dilatable_base.h>
#include <Check_espace_virtuel.h>
#include <Schema_Temps_base.h>
#include <Porosites_champ.h>
#include <EOS_Tools_VEF.h>
#include <Domaine_VEF.h>
#include <Champ_P1NC.h>
#include <Domaine.h>
#include <Debog.h>
 
#include <kokkos++.h>
#include <TRUSTArray_kokkos.tpp>
 
Implemente_instanciable(EOS_Tools_VEF,"EOS_Tools_VEF",EOS_Tools_base);
 
Sortie& EOS_Tools_VEF::printOn(Sortie& os) const
{
  return os <<que_suis_je()<< finl;
}
 
Entree& EOS_Tools_VEF::readOn(Entree& is)
{
  return is;
}
 
/*! @brief Associe les domaines a l'EDO
 *
 * @param (Domaine_dis_base& domaine) domaine
 * @param (Domaine_Cl_dis_base& domaine_cl) domaine cl
 */
void EOS_Tools_VEF::associer_domaines(const Domaine_dis_base& dds, const Domaine_Cl_dis_base& domaine_cl)
{
  le_dom = ref_cast(Domaine_VEF,dds);
  le_dom_Cl = domaine_cl;
  if (!un_.get_md_vector())
    {
      le_dom->creer_tableau_faces(un_, RESIZE_OPTIONS::NOCOPY_NOINIT);
      mapToDevice(un_); // Copy on device this constant array
      un_ = 1.;
    }
}
 
/*! @brief Renvoie rho avec la meme discretisation que la vitesse : une valeur par face en VEF
 *
 * @return (DoubleTab&) rho discretise par face
 */
const DoubleTab& EOS_Tools_VEF::rho_discvit() const
{
  return le_fluide_->masse_volumique().valeurs();
}
 
const DoubleTab& EOS_Tools_VEF::rho_face_n() const
{
  return le_fluide_->rho_n();
}
 
const DoubleTab& EOS_Tools_VEF::rho_face_np1() const
{
  return le_fluide_->rho_np1();
}
 
/*! @brief Calcule la moyenne volumique de la grandeur P1NC donnee
 *
 * @return (DoubleTab&) rho discretise par face
 */
double EOS_Tools_VEF::moyenne_vol(const DoubleTab& tab) const
{
  assert(tab.line_size() == 1);
  double x = Champ_P1NC::calculer_integrale_volumique(le_dom.valeur(), tab, FAUX_EN_PERIO);
  // La facon simple de faire serait celle-ci, mais c'est faux a cause de FAUX_EN_PERIO
  //  qui compte deux fois les volumes entrelaces des faces periodiques:
  double y = Champ_P1NC::calculer_integrale_volumique(le_dom.valeur(), un_, FAUX_EN_PERIO);
  return x / y;
}
 
/*! @brief Div(u) is computed on faces from Div(u) on elements
 *
 */
void EOS_Tools_VEF::divu_discvit(const DoubleTab& DivVelocityElements, DoubleTab& DivVelocityFaces)
{
  ToDo_Kokkos("critical");
  assert_espace_virtuel_vect(DivVelocityElements);
  //Corrections pour moyenner div(u) sur les faces
  const DoubleVect& volumes = le_dom->volumes();
  const DoubleVect& volumes_entrelaces = le_dom->volumes_entrelaces();
  const DoubleVect& volumes_entrelaces_Cl = ref_cast(Domaine_Cl_VEF,le_dom_Cl.valeur()).volumes_entrelaces_Cl();
 
  IntTab& face_voisins = le_dom->face_voisins();
  int premiere_fac_std=le_dom->premiere_face_std();
 
  //remplissage de div(u) sur les faces
  for (int face=0; face<premiere_fac_std; face++)
    {
      int nb_comp=0;
      DivVelocityFaces(face)=0;
      for (int i=0 ; i<2 ; i++)
        {
          int elem= face_voisins(face,i);
          if (elem!=-1)
            {
              nb_comp++;
              DivVelocityFaces(face) += DivVelocityElements(elem)*(volumes_entrelaces_Cl(face)/volumes(elem));
            }
        }
      DivVelocityFaces(face) /= nb_comp;
    }
 
  int size = le_dom->nb_faces();
  for (int face=premiere_fac_std; face<size; face++)
    {
      int nb_comp=0;
      DivVelocityFaces(face)=0;
      for (int i=0 ; i<2 ; i++)
        {
          int elem= face_voisins(face,i);
          if (elem!=-1)
            {
              nb_comp++;
              DivVelocityFaces(face) += DivVelocityElements(elem)*(volumes_entrelaces(face)/volumes(elem));
            }
        }
      DivVelocityFaces(face) /= nb_comp;
    }
}
 
/*! @brief Calcule le second membre de l'equation de continuite pour une discretisation VEF_P1B: div(U) = W = dZ/dT + U.
 *
 * grad(Z)    avec Z=ln(rho)
 *
 * @return rho discretise par face
 */
void EOS_Tools_VEF::secmembre_divU_Z(DoubleTab& tab_W) const
{
  double dt = le_fluide().vitesse().equation().schema_temps().pas_de_temps();
  const IntTab& face_sommets = le_dom->face_sommets();
  int nb_elem_tot = le_dom->nb_elem_tot();
  int nb_som_tot = le_dom->domaine().nb_som_tot();
  int nb_faces_tot = le_dom->nb_faces_tot();
  const Equation_base& eq = le_fluide().vitesse().equation();
 
  const DoubleTab& tab_rhon = le_fluide().loi_etat()->rho_n();
  const DoubleTab& tab_rhonp1 = le_fluide().loi_etat()->rho_np1();
  DoubleTrav tab_rhonP1_, tab_rhonp1P1_;
  const DoubleVect& porosite_face = le_fluide().porosite_face();
  const DoubleTab& tab_rhonP1 = modif_par_porosite_si_flag(tab_rhon, tab_rhonP1_, 1, porosite_face);
  const DoubleTab& tab_rhonp1P1 = modif_par_porosite_si_flag(tab_rhonp1, tab_rhonp1P1_, 1, porosite_face);
  const DoubleVect& volumes_entrelaces = le_dom->volumes_entrelaces();
  const Domaine_VEF& zp1b = ref_cast(Domaine_VEF, le_dom.valeur());
 
  int nfe = le_dom->domaine().nb_faces_elem();
  int nsf = le_dom->nb_som_face();
  const Domaine& dom = le_dom->domaine();
 
  // calcul de la somme des volumes entrelacees autour d'un sommet
  DoubleTrav volume_int_som(nb_som_tot);
  CIntTabView face_sommets_v = face_sommets.view_ro();
  CDoubleArrView volumes_entrelaces_v = static_cast<const DoubleVect&>(volumes_entrelaces).view_ro();
  CIntArrView renum_som_perio_v = dom.get_renum_som_perio().view_ro();
  DoubleArrView volume_int_som_v = static_cast<DoubleVect&>(volume_int_som).view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_faces_tot, KOKKOS_LAMBDA(
                         const int face)
  {
    for (int som = 0; som < nsf; som++)
      {
        int som_glob = renum_som_perio_v(face_sommets_v(face, som));
        Kokkos::atomic_add(&volume_int_som_v(som_glob), volumes_entrelaces_v(face));
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  //discretisation de rho sur les sommets
  DoubleTrav tab_rhon_som(nb_som_tot);
  DoubleTrav tab_rhonp1_som(nb_som_tot);
 
  CDoubleTabView tab_rhonP1_v = tab_rhonP1.view_ro();
  CDoubleTabView tab_rhonp1P1_v = tab_rhonp1P1.view_ro();
  DoubleArrView tab_rhon_som_v = static_cast<DoubleVect&>(tab_rhon_som).view_rw();
  DoubleArrView tab_rhonp1_som_v = static_cast<DoubleVect&>(tab_rhonp1_som).view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_faces_tot, KOKKOS_LAMBDA(
                         const int face)
  {
    for (int som = 0; som < nsf; som++)
      {
        int som_glob = renum_som_perio_v(face_sommets_v(face, som));
        // double pond = 1./nsf; // version_originale
        double pond = volumes_entrelaces_v(face) / volume_int_som_v(som_glob);
        Kokkos::atomic_add(&tab_rhon_som_v(som_glob), tab_rhonP1_v(face, 0) * pond);
        Kokkos::atomic_add(&tab_rhonp1_som_v(som_glob), tab_rhonp1P1_v(face, 0) * pond);
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
//Corrections pour test de la moyenne de la derivee de la masse volumique
  // Dimensionnement de tab_dZ
  const Navier_Stokes_std& eqns = ref_cast(Navier_Stokes_std, eq);
  const DoubleVect& pression = eqns.pression().valeurs();
  Debog::verifier("EOS_Tools_VEF::secmembre_divU_Z pression=",pression);
  Debog::verifier("EOS_Tools_VEF::secmembre_divU_Z tab_rhonP1=",tab_rhonP1);
  Debog::verifier("EOS_Tools_VEF::secmembre_divU_Z tab_rhonp1P1=",tab_rhonp1P1);
  Debog::verifier("EOS_Tools_VEF::secmembre_divU_Z dt=",dt);
 
  int p_has_elem=zp1b.get_alphaE();
  int p_has_som=zp1b.get_alphaS();
  int p_has_arrete=zp1b.get_alphaA();
  int nb_ar_tot = le_dom->domaine().nb_aretes_tot();
 
  // ToDo_Kokkos, try to merge all theses kernels for efficiency:
  int decal=0;
  DoubleTrav tab_dZ = pression;
  DoubleArrView tab_dZ_v = static_cast<DoubleVect&>(tab_dZ).view_rw();
  if (p_has_elem)
    {
      CIntTabView elem_faces_v = le_dom->elem_faces().view_ro();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem_tot, KOKKOS_LAMBDA(
                             const int elem)
      {
        double rn = 0;
        double rnp1 = 0;
        for (int face = 0; face < nfe; face++)
          {
            rn += tab_rhonP1_v(elem_faces_v(elem, face), 0);
            rnp1 += tab_rhonp1P1_v(elem_faces_v(elem, face), 0);
          }
        tab_dZ_v(elem) = (rnp1 - rn) / (nfe * dt);
      });
      end_gpu_timer(__KERNEL_NAME__);
      decal += nb_elem_tot;
    }
 
  if (p_has_som)
    {
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_som_tot, KOKKOS_LAMBDA(
                             const int som)
      {
        tab_dZ_v(decal + som) = ((tab_rhonp1_som_v(som)) - (tab_rhon_som_v(som))) / dt;
      });
      end_gpu_timer(__KERNEL_NAME__);
      decal += nb_som_tot;
    }
 
  if (p_has_arrete)
    {
      for (int ar = 0; ar < nb_ar_tot; ar++)
        tab_dZ(decal + ar) = 0;
    }
  tab_dZ.echange_espace_virtuel();
  Debog::verifier("EOS_Tools_VEF::secmembre_divU_Z tab_dZ=",tab_dZ);
 
  // Ajout des termes sources speciaux de l'equation de masse:
  const bool has_mass_flux = (sub_type(Navier_Stokes_Fluide_Dilatable_base, le_fluide().vitesse().equation())) ?
                             ref_cast(Navier_Stokes_Fluide_Dilatable_base, le_fluide().vitesse().equation()).has_source_masse() : false;
 
  if (has_mass_flux)
    {
      const Source_Masse_Fluide_Dilatable_base& src_mass = ref_cast(Navier_Stokes_Fluide_Dilatable_base, le_fluide().vitesse().equation()).source_masse();
      src_mass.ajouter_projection(le_fluide(),tab_dZ); // attention : tab_dZ a taille comme pression
    }
 
  CDoubleArrView dZ = static_cast<DoubleVect&>(tab_dZ).view_ro();
  DoubleTabView W = tab_W.view_rw();
  decal=0;
  if (p_has_elem)
    {
      double coefdivelem=1;
      CDoubleArrView volumes_v = le_dom->volumes().view_ro();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem_tot, KOKKOS_LAMBDA(
                             const int elem)
      {
        W(elem, 0) = -coefdivelem * dZ(elem) * volumes_v(elem);
      });
      end_gpu_timer(__KERNEL_NAME__);
      decal+=nb_elem_tot;
    }
  if (p_has_som)
    {
      double coefdivsom=1;
      CDoubleArrView volumes_controle_v = zp1b.volume_aux_sommets().view_ro();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_som_tot, KOKKOS_LAMBDA(
                             const int som)
      {
        W(decal + som, 0) = -coefdivsom * dZ(decal + som) * volumes_controle_v(som);
      });
      end_gpu_timer(__KERNEL_NAME__);
      decal+=nb_som_tot;
    }
  if (p_has_arrete)
    {
      for (int ar = 0; ar < nb_ar_tot; ar++)
        {
          //pour l'instant on n'a pas calcule tab_dZ aux aretes
          assert(tab_dZ(decal + ar) == 0);
          tab_W(decal + ar) = 0.;
        }
    }
  tab_W.echange_espace_virtuel();
  Debog::verifier("EOS_Tools_VEF::secmembre_divU_Z tab_W=",tab_W);
}