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 Associates the domains to the ODE.
 *
 * @param (Domaine_dis_base& domaine) domain
 * @param (Domaine_Cl_dis_base& domaine_cl) boundary condition domain
 */
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 Returns rho with the same discretization as the velocity: one value per face in VEF.
 *
 * @return (DoubleTab&) rho discretized per 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 Computes the volumetric average of the given P1NC quantity.
 *
 * @return (DoubleTab&) rho discretized per 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);
  // The simple approach would be to compute this directly, but it is incorrect because of FAUX_EN_PERIO
  //  which counts the interlaced volumes of periodic faces twice:
  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 to average div(u) on 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();
 
  // Fill div(u) on 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 Computes the right-hand side of the continuity equation for a VEF_P1B discretization: div(U) = W = dZ/dT + U.
 *
 * grad(Z)    with Z=ln(rho)
 *
 * @return rho discretized per 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();
 
  // compute the sum of interlaced volumes around a vertex
  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__);
 
  // discretization of rho on the vertices
  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 for testing the mean of the density derivative
  // Sizing of 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);
 
  // Adding special source terms of the mass equation:
  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); // note: tab_dZ has the same size as the pressure array
    }
 
  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++)
        {
          // for now, tab_dZ has not been computed at edges
          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);
}