Champ_P1NC.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 <Dirichlet_paroi_defilante.h>
#include <Echange_externe_impose.h>
#include <Scalaire_impose_paroi.h>
#include <Dirichlet_paroi_fixe.h>
#include <Modele_turbulence_hyd_base.h>
#include <Schema_Temps_base.h>
#include <Neumann_homogene.h>
#include <Operateur_base.h>
#include <Champ_Uniforme.h>
#include <Neumann_paroi.h>
#include <distances_VEF.h>
#include <Probleme_base.h>
#include <Fluide_base.h>
#include <Periodique.h>
#include <Champ_P1NC.h>
#include <Operateur.h>
#include <TRUSTTab.h>
#include <TRUSTTrav.h>
#include <SFichier.h>
#include <Device.h>
#include <kokkos++.h>
#include <TRUSTArray_kokkos.tpp>
 
Implemente_instanciable(Champ_P1NC,"Champ_P1NC",Champ_Inc_base);
 
Sortie& Champ_P1NC::printOn(Sortie& s) const { return s << que_suis_je() << " " << le_nom(); }
 
Entree& Champ_P1NC::readOn(Entree& s)
{
  lire_donnees(s);
  return s;
}
 
// PQ : 01/10/04 : mettre_a_jour computes only once the conforming field u_bar (defined at vertices)
//                   so it can be reused in the various operators or when calling
//                   probes defined from "smooth" values (in the present case: u_bar)
void Champ_P1NC::mettre_a_jour(double un_temps)
{
  Champ_Inc_base::mettre_a_jour(un_temps);
  filtrer_L2_deja_appele_ = 0;
}
/*! @brief
 *
 */
void Champ_P1NC::abortTimeStep()
{
  ch_som_ = 0;
}
 
int Champ_P1NC::compo_normale_sortante(int num_face) const
{
  double vit_norm = 0;
  for (int ncomp = 0; ncomp < nb_comp(); ncomp++)
    vit_norm += valeurs()(num_face, ncomp) * domaine_vef().face_normales(num_face, ncomp);
  return (vit_norm > 0);
}
 
DoubleTab& Champ_P1NC::trace(const Frontiere_dis_base& fr, DoubleTab& x, double tps, int distant) const
{
  return Champ_P1NC_implementation::trace(fr, valeurs(tps), x, distant);
}
 
void Champ_P1NC::cal_rot_ordre1(DoubleTab& tab_vorticite) const
{
  const int nb_elem = domaine_vef().nb_elem(), nb_elem_tot = domaine_vef().nb_elem_tot();
  DoubleTrav tab_gradient_elem(nb_elem_tot, dimension, dimension);
  gradient(tab_gradient_elem);
  CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
  switch(dimension)
    {
    case 2:
      {
        DoubleArrView vorticite = static_cast<DoubleVect&>(tab_vorticite).view_wo();
        Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem, KOKKOS_LAMBDA(const int num_elem)
        {
          vorticite(num_elem) = gradient_elem(num_elem, 1, 0) - gradient_elem(num_elem, 0, 1);
        });
      }
      break;
    case 3:
      {
        DoubleTabView vorticite = tab_vorticite.view_wo();
        Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem, KOKKOS_LAMBDA(const int num_elem)
        {
          vorticite(num_elem, 0) = gradient_elem(num_elem, 2, 1) - gradient_elem(num_elem, 1, 2);
          vorticite(num_elem, 1) = gradient_elem(num_elem, 0, 2) - gradient_elem(num_elem, 2, 0);
          vorticite(num_elem, 2) = gradient_elem(num_elem, 1, 0) - gradient_elem(num_elem, 0, 1);
        });
      }
    }
  end_gpu_timer(__KERNEL_NAME__);
  tab_vorticite.echange_espace_virtuel();
  return;
}
 
void calculer_gradientP1NC_2D(const DoubleTab& tab_variable, const Domaine_VEF& domaine_VEF, const Domaine_Cl_VEF& domaine_Cl_VEF, DoubleTab& tab_gradient_elem)
{
  const DoubleTab& tab_face_normales = domaine_VEF.face_normales();
  const IntTab& tab_face_voisins = domaine_VEF.face_voisins();
  const DoubleVect& tab_inverse_volumes = domaine_VEF.inverse_volumes();
  const ArrOfInt& tab_est_face_bord = domaine_VEF.est_face_bord();
 
  CDoubleTabView face_normales = tab_face_normales.view_ro();
  CIntTabView face_voisins = tab_face_voisins.view_ro();
  CDoubleArrView inverse_volumes = tab_inverse_volumes.view_ro();
  CDoubleArrView variable = static_cast<const ArrOfDouble&>(tab_variable).view_ro();
  CIntArrView est_face_bord = tab_est_face_bord.view_ro();
  DoubleTabView gradient_elem = tab_gradient_elem.view_rw();
 
  int dimension = Objet_U::dimension;
  const int nb_faces_tot = domaine_VEF.nb_faces_tot();
  const int nb_elem = domaine_VEF.nb_elem_tot();
 
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       nb_faces_tot,
                       KOKKOS_LAMBDA (int fac)
  {
    int type_face = est_face_bord(fac);
 
    // Periodic boundary faces
    if (type_face == 2)
      {
        int elem1 = face_voisins(fac, 0);
        int elem2 = face_voisins(fac, 1);
 
        for (int i = 0; i < dimension; i++)
          {
            double grad = 0.5 * face_normales(fac,i) * variable(fac);
            Kokkos::atomic_add(&gradient_elem(elem1, i), +grad);
            Kokkos::atomic_add(&gradient_elem(elem2, i), -grad);
          }
      }
    // Non-periodic boundary faces
    else if (type_face == 1)
      {
        int elem1 = face_voisins(fac,0);
 
        for (int i = 0; i < dimension; i++)
          {
            double grad = face_normales(fac,i) * variable(fac);
            Kokkos::atomic_add(&gradient_elem(elem1, i), grad);
          }
      }
    // Internal faces
    else if (type_face == 0)
      {
        int elem1 = face_voisins(fac, 0);
        int elem2 = face_voisins(fac, 1);
 
        for (int i = 0; i < dimension; i++)
          {
            double grad = face_normales(fac,i) * variable(fac);
            if (elem1 >= 0) Kokkos::atomic_add(&gradient_elem(elem1, i), +grad);
            if (elem2 >= 0) Kokkos::atomic_add(&gradient_elem(elem2, i), -grad);
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  // Divide by the element volume
  auto vol_div = KOKKOS_LAMBDA (int elem, int i)
  {
    gradient_elem(elem, i) *= inverse_volumes(elem);
  };
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_2D({0,0}, {nb_elem,dimension}), vol_div);
  end_gpu_timer(__KERNEL_NAME__);
}
 
//Actually optimized kernel
template<int nbcomp_compile_time, int dims, int nbface>
void calculer_gradientP1NC_3D_template(const DoubleTab& tab_variable,
                                       const Domaine_VEF& domaine_VEF,
                                       const Domaine_Cl_VEF& domaine_Cl_VEF,
                                       DoubleTab& tab_gradient_elem)
{
  const DoubleTab& tab_face_normales = domaine_VEF.face_normales();
  const IntTab& tab_face_voisins = domaine_VEF.face_voisins();
  const DoubleVect& tab_inverse_volumes = domaine_VEF.inverse_volumes();
  const IntTab& tab_elem_faces = domaine_VEF.elem_faces();
 
  int nb_elem_tot = domaine_VEF.nb_elem_tot();
 
  CDoubleTabView face_normales = tab_face_normales.view_ro();
  CIntTabView face_voisins = tab_face_voisins.view_ro();
  CDoubleArrView inverse_volumes = tab_inverse_volumes.view_ro();
  CDoubleTabView variable = tab_variable.view_ro();
  CIntTabView elem_faces = tab_elem_faces.view_ro();
  DoubleTabView3 gradient_elem = tab_gradient_elem.view_rw<3>();
 
 
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem_tot, KOKKOS_LAMBDA (int elem)
  {
 
    double contrib_grad_loc[nbcomp_compile_time*dims]= {};
 
    for (int iface = 0; iface < nbface; iface++)
      {
        int fac = elem_faces(elem, iface);
 
        double coef = (elem == face_voisins(fac, 0) ? 1 : -1);
 
        double face_normales_loc[dims];
        double variable_loc[nbcomp_compile_time];
 
        for (int i = 0; i < dims; i++)
          {
            face_normales_loc[i] = face_normales(fac, i);
          }
 
        for (int icomp = 0; icomp < nbcomp_compile_time; icomp++)
          {
            variable_loc[icomp] = variable(fac, icomp);
          }
 
        for (int icomp = 0; icomp < nbcomp_compile_time; icomp++)
          for (int i = 0; i < dims; i++)
            {
              double grad = coef * face_normales_loc[i] * variable_loc[icomp];
              contrib_grad_loc[dims*icomp + i] += grad;
            }
      }
 
    double inverse = inverse_volumes(elem);
 
    for (int icomp = 0; icomp < nbcomp_compile_time; icomp++)
      for (int i = 0; i < dims; i++)
        {
          gradient_elem(elem, icomp, i) = contrib_grad_loc[dims*icomp + i]*inverse;
        }
 
 
  });
  end_gpu_timer(__KERNEL_NAME__);
}
 
// Helper for variadic template dispatch
template <int ...> struct TemplateIntList {};
 
// Base case: no more cases to try
template<int dims, int nbface>
void dispatch_nbcomp_impl(int, TemplateIntList<>, const DoubleTab&, const Domaine_VEF&,
                          const Domaine_Cl_VEF&, DoubleTab&)
{
  Cerr << "Error in calculer_gradientP1NC_3D: unsupported nb_comp, add more in 'dispatch_nbcomp'" << finl;
  Process::exit();
}
 
// Recursive case: try each nbcomp value
template<int dims, int nbface, int nbcomp_compile_time, int ...Rest>
void dispatch_nbcomp_impl(int nb_comp_runtime, TemplateIntList<nbcomp_compile_time, Rest...>,
                          const DoubleTab& tab_variable,
                          const Domaine_VEF& domaine_VEF,
                          const Domaine_Cl_VEF& domaine_Cl_VEF,
                          DoubleTab& tab_gradient_elem)
{
  if (nbcomp_compile_time == nb_comp_runtime)
    {
      calculer_gradientP1NC_3D_template<nbcomp_compile_time, dims, nbface>(tab_variable, domaine_VEF,
                                                                           domaine_Cl_VEF, tab_gradient_elem);
    }
  else
    {
      dispatch_nbcomp_impl<dims, nbface>(nb_comp_runtime, TemplateIntList<Rest...>(),
                                         tab_variable, domaine_VEF, domaine_Cl_VEF, tab_gradient_elem);
    }
}
 
// Wrapper for nbcomp dispatch
template<int dims, int nbface>
void dispatch_nbcomp(int nb_comp_runtime, const DoubleTab& tab_variable,
                     const Domaine_VEF& domaine_VEF,
                     const Domaine_Cl_VEF& domaine_Cl_VEF,
                     DoubleTab& tab_gradient_elem)
{
  // List all supported nbcomp values here - easy to extend!
  dispatch_nbcomp_impl<dims, nbface>(nb_comp_runtime, TemplateIntList<1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20>(),
                                     tab_variable, domaine_VEF, domaine_Cl_VEF, tab_gradient_elem);
}
 
void calculer_gradientP1NC_3D(const DoubleTab& tab_variable,
                              const Domaine_VEF& domaine_VEF,
                              const Domaine_Cl_VEF& domaine_Cl_VEF,
                              DoubleTab& tab_gradient_elem)
{
  const int nb_comp = tab_variable.line_size();
  const int dimensions = Objet_U::dimension;
  int nb_faces_elem = domaine_VEF.elem_faces().dimension(1);
 
 
  //Switch to allow compile time optimizations https://rbourgeois33.github.io./posts/post1/ namely:
  //no atomics, no redundant memory requests
  //Weird but this kernel is also called for some 2D cases..., exemple Bilans_VEF_QC_Turb_Null
  //Also weird, P1NC should always be called on tetra so we could deduce nb_faces_elem from dim (3 in 2D, 4 in 3D)
  //But in triocfd, it is called with quadrangles/hexahedron... tests Smago_hexa Smago_quadra Obstacle_Turb_quadra
 
  switch (nb_faces_elem)
    {
    case 3:
      switch (dimensions)
        {
        case 2:
          dispatch_nbcomp<2,3>(nb_comp, tab_variable,
                               domaine_VEF, domaine_Cl_VEF,
                               tab_gradient_elem);
          break;
 
        case 3:
          dispatch_nbcomp<3,3>(nb_comp, tab_variable,
                               domaine_VEF, domaine_Cl_VEF,
                               tab_gradient_elem);
          break;
 
        default:
          Cerr << "Error in calculer_gradientP1NC_3D: dimensions must be 2 or 3, got "
               << dimensions << finl;
          Process::exit();
        }
      break;
 
    case 4:
      switch (dimensions)
        {
        case 2:
          dispatch_nbcomp<2,4>(nb_comp, tab_variable,
                               domaine_VEF, domaine_Cl_VEF,
                               tab_gradient_elem);
          break;
 
        case 3:
          dispatch_nbcomp<3,4>(nb_comp, tab_variable,
                               domaine_VEF, domaine_Cl_VEF,
                               tab_gradient_elem);
          break;
 
        default:
          Cerr << "Error in calculer_gradientP1NC_3D: dimensions must be 2 or 3, got "
               << dimensions << finl;
          Process::exit();
        }
      break;
 
    case 6:
      switch (dimensions)
        {
        case 2:
          dispatch_nbcomp<2,6>(nb_comp, tab_variable,
                               domaine_VEF, domaine_Cl_VEF,
                               tab_gradient_elem);
          break;
 
        case 3:
          dispatch_nbcomp<3,6>(nb_comp, tab_variable,
                               domaine_VEF, domaine_Cl_VEF,
                               tab_gradient_elem);
          break;
 
        default:
          Cerr << "Error in calculer_gradientP1NC_3D: dimensions must be 2 or 3, got "
               << dimensions << finl;
          Process::exit();
        }
      break;
 
    default:
      Cerr << "Error in calculer_gradientP1NC_3D: nb_faces_elem must be 3, 4, or 6 got "
           << nb_faces_elem << finl;
      Process::exit();
    }
}
 
void calculer_gradientP1NC(const DoubleTab& tab_variable, const Domaine_VEF& domaine_VEF, const Domaine_Cl_VEF& domaine_Cl_VEF, DoubleTab& tab_gradient_elem)
{
  const int nb_comp = tab_variable.line_size();
  tab_gradient_elem = 0.;
 
  // Case: gradient of a vector array or scalar array tab_gradient_elem(,,)
  if (nb_comp != 1 || tab_gradient_elem.nb_dim() == 3)
    calculer_gradientP1NC_3D(tab_variable, domaine_VEF, domaine_Cl_VEF, tab_gradient_elem);
  // Case: gradient of a scalar array stored in gradient_elem(,)
  else
    calculer_gradientP1NC_2D(tab_variable, domaine_VEF, domaine_Cl_VEF, tab_gradient_elem);
}
 
void Champ_P1NC::gradient(DoubleTab& gradient_elem) const
{
  // Compute the gradient: e.g. velocity gradient for vorticity computation
  const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF, equation().domaine_Cl_dis());
  calculer_gradientP1NC(valeurs(), domaine_vef(), domaine_Cl_VEF, gradient_elem);
 
}
 
int Champ_P1NC::imprime(Sortie& os, int ncomp) const
{
  imprime_P1NC(os, ncomp);
  return 1;
}
 
Champ_base& Champ_P1NC::affecter_(const Champ_base& ch)
{
#ifdef CLEMENT
  DoubleTab& tab = valeurs();
  DoubleTab noeuds;
  remplir_coord_noeuds(noeuds);
  const Domaine_VEF& domaine_VEF = le_dom_VEF.valeur();
  const int nb_faces = domaine_VEF.nb_faces();
  const IntTab& face_sommets = domaine_VEF.face_sommets();
  const DoubleTab& sommets=domaine().coord_sommets();
  int nb_som=sommets.dimension(0);
  Cerr << "Affecter Clement : ... " << finl;
  DoubleTab valf(valeurs());
  ch.valeur_aux(noeuds, valf);
  DoubleTab vals(nb_som, nb_compo_);
  if(nb_compo_==1)
    {
      vals.resize(nb_som);
    }
  ch.valeur_aux(sommets, vals);
  for(int face=0; face<nb_faces; face++)
    {
      for(int comp=0; comp<nb_compo_; comp++)
        {
          double somme=0;
          for (int isom=0; isom<dimension; isom++)
            if(nb_compo_==1)
              somme+=vals(face_sommets(face,isom));
            else
              somme+=vals(face_sommets(face,isom), comp);
          if(nb_compo_==1)
            {
              if(dimension==2)
                tab(face)=1./6.*somme+2./3.*valf(face);
              else
                tab(face)=17./60.*somme+3./20.*valf(face);
            }
          else if(dimension==2)
            tab(face,comp)=1./6.*somme+2./3.*valf(face,comp);
          else
            tab(face,comp)=17./60.*somme+3./20.*valf(face,comp);
        }
    }
  return *this;
#else
  return Champ_Inc_base::affecter_(ch);
#endif
}
 
// Periodic BC case: ensures that values on facing periodic boundary faces are identical.
// To do so, the half-sum of the two values is used.
void Champ_P1NC::verifie_valeurs_cl()
{
  const Domaine_Cl_dis_base& zcl = equation().domaine_Cl_dis();
  int nb_cl = zcl.nb_cond_lim();
  DoubleTab& ch_tab = valeurs();
  int nb_compo = nb_comp();
  int ndeb, nfin, num_face;
 
  for (int i = 0; i < nb_cl; i++)
    {
      const Cond_lim_base& la_cl = zcl.les_conditions_limites(i).valeur();
      if (sub_type(Periodique, la_cl))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique, la_cl);
          const Front_VF& le_bord = ref_cast(Front_VF, la_cl.frontiere_dis());
          ndeb = le_bord.num_premiere_face();
          nfin = ndeb + le_bord.nb_faces();
          int voisine;
          double moy;
 
          for (num_face = ndeb; num_face < nfin; num_face++)
            {
              voisine = la_cl_perio.face_associee(num_face - ndeb) + ndeb;
              for (int comp = 0; comp < nb_compo; comp++)
                {
                  if (ch_tab(num_face, comp) != ch_tab(voisine, comp))
                    {
                      moy = 0.5 * (ch_tab(num_face, comp) + ch_tab(voisine, comp));
                      ch_tab(num_face, comp) = moy;
                      ch_tab(voisine, comp) = moy;
                    }
                }
            }
        }
    }
  ch_tab.echange_espace_virtuel();
}
 
void Champ_P1NC::calcul_critere_Q(DoubleVect& tab_Critere_Q) const
{
  const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF, equation().domaine_Cl_dis());
  const int nb_elem = domaine_vef().nb_elem(), nb_elem_tot = domaine_vef().nb_elem_tot();
  const int dim = Objet_U::dimension;
 
  DoubleTrav tab_gradient_elem(nb_elem_tot, dim, dim);
  const DoubleTab& vitesse = valeurs();
  Champ_P1NC::calcul_gradient(vitesse, tab_gradient_elem, domaine_Cl_VEF);
 
  CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
  DoubleArrView Critere_Q = tab_Critere_Q.view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_1D(0, nb_elem), KOKKOS_LAMBDA(
                         const int num_elem)
  {
    double crit = 0.;
    for (int i = 0; i < dim; i++)
      for (int j = 0; j < dim; j++)
        {
          double deriv1 = gradient_elem(num_elem, i, j);
          double deriv2 = gradient_elem(num_elem, j, i);
          crit += -0.25 * deriv1 * deriv2;
        }
    Critere_Q[num_elem] = crit;
  });
  end_gpu_timer(__KERNEL_NAME__);
}
 
void Champ_P1NC::calcul_y_plus(const Domaine_Cl_VEF& domaine_Cl_VEF, DoubleVect& y_plus) const
{
  // Initialize the y_plus field with a negative value,
  // so that when visualizing the field near the wall,
  // one only needs to remove the negative values and only
  // the values at the walls remain.
 
  int ndeb, nfin, l_unif;
  double visco0 = -1.;
  y_plus = -1.;
 
  const Equation_base& eqn_hydr = equation();
  const Fluide_base& le_fluide = ref_cast(Fluide_base, eqn_hydr.milieu());
  const Champ_Don_base& ch_visco_cin = le_fluide.viscosite_cinematique();
  const DoubleTab& tab_visco = ch_visco_cin.valeurs();
 
  if (sub_type(Champ_Uniforme, ch_visco_cin))
    {
      if (tab_visco(0, 0) > DMINFLOAT)
        visco0 = tab_visco(0, 0);
      else
        visco0 = DMINFLOAT;
      l_unif = 1;
    }
  else
    l_unif = 0;
 
  if ((!l_unif) && (tab_visco.local_min_vect() < DMINFLOAT))
    // GF we must not modify tab_visco here!
    {
      Cerr << "In Champ_P1NC::calcul_y_plus : visco = " << tab_visco.local_min_vect() << " <= 0 ? " << finl;
      exit();
    }
  // tab_visco+=DMINFLOAT;
 
  DoubleTab yplus_faces(1, 1); // will contain yplus values if available
  int yplus_already_computed = 0; // flag
 
  const RefObjU& modele_turbulence = eqn_hydr.get_modele(TURBULENCE);
  if (modele_turbulence && sub_type(Modele_turbulence_hyd_base, modele_turbulence.valeur()))
    {
      const Modele_turbulence_hyd_base& mod_turb = ref_cast(Modele_turbulence_hyd_base, modele_turbulence.valeur());
      const Turbulence_paroi_base& loipar = mod_turb.loi_paroi();
      // if( ! sub_type( Paroi_negligeable_VEF , loipar ) )
      if (loipar.use_shear())
        {
          yplus_faces.resize(domaine_vef().nb_faces_tot());
          yplus_faces.ref(loipar.tab_d_plus());
          yplus_already_computed = 1;
        }
    }
 
  for (int n_bord = 0; n_bord < domaine_vef().nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
 
      if( sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) || sub_type(Dirichlet_paroi_defilante, la_cl.valeur()) || la_cl->que_suis_je() == "Frontiere_ouverte_vitesse_imposee_ALE")
        {
          const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
          ndeb = le_bord.num_premiere_face();
          nfin = ndeb + le_bord.nb_faces();
          int dim = Objet_U::dimension;
          int nfac = domaine_vef().domaine().nb_faces_elem();
          DoubleTrav counter(y_plus.size());
          CDoubleTabView xp = domaine_vef().xp().view_ro(); // DOUBT: Why defined inside the for loop?
          CDoubleTabView xv = domaine_vef().xv().view_ro();
          CDoubleTabView face_normale = domaine_vef().face_normales().view_ro();
          CIntTabView elem_faces = domaine_vef().elem_faces().view_ro();
          CIntTabView face_voisins = domaine_vef().face_voisins().view_ro();
          CDoubleTabView vit = eqn_hydr.inconnue().valeurs().view_ro();  // DOUBT: Is this correct? Original assignment was `const Champ_P1NC& vit = *this;`
          CDoubleArrView visco = static_cast<const DoubleVect&>(tab_visco).view_ro();
          CDoubleArrView yp_faces = static_cast<const DoubleVect&>(yplus_faces).view_ro();
          DoubleArrView y_plus_view = y_plus.view_rw();
          DoubleArrView counter_view = static_cast<DoubleVect&>(counter).view_wo();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA (const int num_face)
          {
            int elem = face_voisins(num_face, 0);
            y_plus_view(elem) = 0;
          });
          end_gpu_timer(__KERNEL_NAME__);
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA (const int num_face)
          {
            int elem = face_voisins(num_face, 0);
            if (yplus_already_computed)
              {
                // y+ is only defined on faces so we take the face value to put in the element
                Kokkos::atomic_add(&y_plus_view(elem), yp_faces(num_face));
              }
            else
              {
                int num[3];
                bool ok = false;
                for (int i = 0; i < dim; i++)
                  {
                    num[i] = elem_faces(elem, i);
                    if (num[i] == num_face && !ok)
                      {
                        num[i] = elem_faces(elem, dim);
                        ok = true;
                      }
                  }
 
                double dist = distance(dim, num_face, elem, xp, xv, face_normale);
                dist *= (dim + 1.) / dim; // to convert to wall distance / midpoint of num[0]-num[1]
                double val[3];
                double norm_v = norm_vit1_lp(dim, vit, num_face, nfac, num, face_normale, val);
                double d_visco = l_unif ? visco0:visco[elem];
 
                // PQ : 01/10/03 : corrections relative to the first version
                double norm_tau = d_visco * norm_v / dist;
                double u_etoile = sqrt(norm_tau);
                Kokkos::atomic_add(&y_plus_view(elem), dist * u_etoile / d_visco);
              } // else yplus already computed
            Kokkos::atomic_add(&counter_view(elem), +1);
          }); // loop on faces
          end_gpu_timer(__KERNEL_NAME__);
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(0, y_plus.size()), KOKKOS_LAMBDA (const int elem)
          {
            if (counter_view(elem)>0) y_plus_view(elem)/=counter_view(elem);
          });
          end_gpu_timer(__KERNEL_NAME__);
 
        } // End of fixed wall
 
    } // End loop over boundary conditions
}
 
void Champ_P1NC::calcul_grad_U(const Domaine_Cl_VEF& domaine_Cl_VEF, DoubleTab& tab_grad_u) const
{
  const DoubleTab& u = valeurs();
  const int nb_elem = domaine_vef().nb_elem(), nb_elem_tot = domaine_vef().nb_elem_tot();
 
  DoubleTrav tab_gradient_elem(nb_elem_tot, dimension, dimension);
  calculer_gradientP1NC(u, domaine_vef(), domaine_Cl_VEF, tab_gradient_elem);
  int dim = dimension;
  CDoubleTabView3 gradient_elem = tab_gradient_elem.view_ro<3>();
  DoubleTabView grad_u = tab_grad_u.view_wo();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(0, nb_elem), KOKKOS_LAMBDA(const int elem)
  {
    int comp = 0;
    for (int i = 0; i < dim; i++)
      {
        for (int j = 0; j < dim; j++)
          {
            grad_u(elem, comp) = gradient_elem(elem, i, j);
            comp++;
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
}
 
void Champ_P1NC::calcul_grad_T(const Domaine_Cl_VEF& domaine_Cl_VEF, DoubleTab& grad_T) const
{
  const DoubleTab& temp = valeurs();
  calculer_gradientP1NC(temp, domaine_vef(), domaine_Cl_VEF, grad_T);
 
  // *******************************************************************
  //  PQ : 12/12/05 : to be able to visualize oscillations from
  //                    sign changes in the gradient of T
  // *******************************************************************
  /*
   const IntTab& face_voisins = domaine_VEF.face_voisins();
   const IntTab& elem_faces=domaine_VEF.elem_faces();
   int nb_faces = domaine_VEF.nb_faces();
   const int nb_elem = domaine_VEF.nb_elem_tot();
 
   int num_elem,num_face,elem0,elem1,i;
   ArrOfInt sign_opp_x(nb_faces);
   ArrOfInt sign_opp_y(nb_faces);
   sign_opp_x=0;
   sign_opp_y=0;
 
   for (num_face=0;num_face<nb_faces;num_face++)
   {
   elem0 = face_voisins(num_face,0);
   elem1 = face_voisins(num_face,1);
   if (elem1!=-1)
   {
   if  (grad_T(elem0,0)*grad_T(elem1,0)<0.) sign_opp_x(num_face)++;
   if  (grad_T(elem0,1)*grad_T(elem1,1)<0.) sign_opp_y(num_face)++;
   }
   }
 
   grad_T=0.;
 
   for (num_elem=0;num_elem<nb_elem;num_elem++)
   {
   for (i=0;i<dimension+1;i++)
   {
   grad_T(num_elem,0) +=sign_opp_x(elem_faces(num_elem,i));
   grad_T(num_elem,1) +=sign_opp_y(elem_faces(num_elem,i));
   }
   }
   */
 
}
 
// Compute the heat transfer coefficient
void Champ_P1NC::calcul_h_conv(const Domaine_Cl_VEF& domaine_Cl_VEF, DoubleTab& h_conv, int temp_ref) const
{
  const DoubleTab& temp = valeurs();
  const IntTab& face_voisins = domaine_vef().face_voisins();
  const DoubleTab& face_normale = domaine_vef().face_normales();
 
  // Retrieve the flux at boundaries
  DoubleTab Flux;
  Flux = equation().operateur(0).l_op_base().flux_bords();
 
  // Initialize to -1
  h_conv = -1.;
 
  double T_ref = double(temp_ref);
 
  for (int n_bord = 0; n_bord < domaine_vef().nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Frontiere_dis_base& la_fr = la_cl->frontiere_dis();
      const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
      int ndeb = le_bord.num_premiere_face();
      int nfin = ndeb + le_bord.nb_faces();
      double h_moy = 0.;
      // Depending on the boundary conditions
      if (sub_type(Echange_externe_impose, la_cl.valeur()))
        {
          const Champ_base& rho = equation().milieu().masse_volumique();
          const Champ_Don_base& Cp = equation().milieu().capacite_calorifique();
          int rho_uniforme = (sub_type(Champ_Uniforme,rho) ? 1 : 0);
          int cp_uniforme = (sub_type(Champ_Uniforme,Cp) ? 1 : 0);
          for (int num_face = ndeb; num_face < nfin; num_face++)
            {
              int elem = face_voisins(num_face, 0);
              double rho_cp = (rho_uniforme ? rho.valeurs()(0, 0) : (rho.nb_comp() == 1 ? rho.valeurs()(num_face) : rho.valeurs()(num_face, 0)));
              rho_cp *= (cp_uniforme ? Cp.valeurs()(0, 0) : (Cp.nb_comp() == 1 ? Cp.valeurs()(num_face) : Cp.valeurs()(num_face, 0)));
              h_conv(elem) = ref_cast(Echange_externe_impose,la_cl.valeur()).h_imp(num_face) * rho_cp;
              h_moy += h_conv(elem);
            }
        }
      else if (sub_type(Scalaire_impose_paroi,la_cl.valeur()) || sub_type(Neumann_homogene, la_cl.valeur()) || sub_type(Neumann_paroi, la_cl.valeur()))
        {
          for (int num_face = ndeb; num_face < nfin; num_face++)
            {
              int elem = face_voisins(num_face, 0);
 
              double norme_flux = 0.;
              double surf = 0.;
 
              for (int k = 0; k < Flux.dimension(1); k++)
                norme_flux += Flux(num_face, k) * Flux(num_face, k);
 
              for (int k = 0; k < Objet_U::dimension; k++)
                surf += face_normale(num_face, k) * face_normale(num_face, k);
 
              double DT = temp(num_face) - T_ref;
 
              if (DT == 0.)
                h_conv(elem) = 0.;
              else
                {
                  //This field is updated when its post-processing is requested,
                  //and consequently the fluxes have been updated and are already multiplied by rho*Cp
                  h_conv(elem) = sqrt(norme_flux / surf) / DT;
                }
              h_moy += h_conv(elem);
            }
        }
      // Optimization: combine 2 mp_sum into 1 collective call
      double nb_faces_d = static_cast<double>(le_bord.nb_faces());
      mp_sum_for_each(h_moy, nb_faces_d);
      int nb_faces = static_cast<int>(nb_faces_d);  // a single 'bord' should not have so many faces
      h_moy /= nb_faces;
      if (je_suis_maitre())
        {
          Nom fich_hconv;
          fich_hconv = "Heat_transfert_Coeff_";
          fich_hconv += la_fr.le_nom();
          fich_hconv += ".dat";
 
          SFichier fic(fich_hconv, ios::app);
          fic.setf(ios::scientific);
 
          double le_temps = equation().schema_temps().temps_courant();
          fic << le_temps << " " << h_moy << finl;
        }
    }
}
 
static double norme_L2(const DoubleTab& u, const Domaine_VEF& domaine_VEF)
{
  const int nb_faces = domaine_VEF.nb_faces();
  double norme = 0;
  int nb_compo_ = u.line_size();
 
  CDoubleArrView volumes = domaine_VEF.volumes_entrelaces().view_ro(); // DOUBT: reference?
  CDoubleTabView uview = u.view_ro(); // DOUBT: How to name this variable?
 
  Kokkos::parallel_reduce(start_gpu_timer(__KERNEL_NAME__), nb_faces,
                          KOKKOS_LAMBDA (const int i, double& s)
  {
    for (int j = 0; j < nb_compo_; j++)
      s += uview(i, j) * uview(i, j);
    s *= volumes(i);
  }, norme);
  end_gpu_timer(__KERNEL_NAME__);
 
  return sqrt(norme);
}
double Champ_P1NC::norme_L2(const Domaine& dom) const
{
  Cerr << "Champ_P1NC::norme_L2" << finl;
  {
    DoubleTab x(valeurs());
    filtrer_L2(x);
    Cerr << "Norme filtree : " << ::norme_L2(x, domaine_vef()) << finl;
  }
  const DoubleTab& u = valeurs();
  return ::norme_L2(u, domaine_vef());
}
 
double Champ_P1NC::norme_H1(const Domaine& dom) const
{
  const Domaine& mon_dom = domaine_vef().domaine();
 
  double dnorme_H1;
  // DOUBT: dimension?
  int dim = dimension;
  const int nb_elem = mon_dom.nb_elem();
  const int nb_faces_elem = mon_dom.nb_faces_elem();
  CIntTabView elem_faces = domaine_vef().elem_faces().view_ro();
  CDoubleTabView face_normales = domaine_vef().face_normales().view_ro();
  CIntTabView face_voisins = domaine_vef().face_voisins().view_ro();
  CDoubleTabView tab = valeurs().view_ro();
  CDoubleArrView volumes = domaine_vef().volumes().view_ro();
 
  // Compute the H1 norm of a P1NC unknown.
  // The algorithm accounts for the following constraints:
  // - the unknown may have several components (i.e. be a scalar or a vector)
  // - the problem dimension is arbitrary (1, 2 or 3).
  // WARNING: prisms are not supported.
  dnorme_H1 = 0.;
  Kokkos::parallel_reduce(start_gpu_timer(__KERNEL_NAME__), nb_comp(), //scalar or vector case
                          KOKKOS_LAMBDA (const int composante, double& norme_H1_comp)
  {
    for (int K = 0; K < nb_elem; K++) //loop over elements
      {
        double norme_grad_elem = 0.; //to avoid accumulation errors
        for (int i = 0; i < dim; i++) //loop over problem dimension
          {
            double int_grad_elem = 0.; //to avoid accumulation errors
            for (int face = 0; face < nb_faces_elem; face++) //loop over faces of element "K"
              {
                int face_globale = elem_faces(K, face);
 
                int_grad_elem += tab(face_globale, composante) * face_normales(face_globale, i) * oriente_normale(face_globale, K, face_voisins);
              } //end for on "face"
 
            norme_grad_elem += int_grad_elem * int_grad_elem;
          } //end for on "i"
 
        norme_H1_comp += norme_grad_elem / volumes(K);
      } //end for on "K"
 
  }, dnorme_H1); // end for on "composante"
  end_gpu_timer(__KERNEL_NAME__);
 
  return sqrt(dnorme_H1);
}
 
double Champ_P1NC::norme_L2_H1(const Domaine& dom) const
{
  double pas_de_temps = equation().schema_temps().pas_de_temps();
 
  return ((1. / pas_de_temps) * norme_L2(dom)) + norme_H1(dom);
}
 
/////////////////////////////////////////////////////////////////////////////////////////////////////
// Method that returns gij at elements from the face field.
// (gij represents the partial derivative dGi/dxj)
// From gij, one can compute Sij = 0.5(gij(i,j)+gij(j,i))
////////////////////////////////////////////////////////////////////////////////////////////////////
 
DoubleTab& Champ_P1NC::calcul_gradient(const DoubleTab& champ, DoubleTab& gij, const Domaine_Cl_VEF& domaine_Cl_VEF)
{
  const Domaine_VEF& domaine_VEF = domaine_Cl_VEF.domaine_vef();
  calculer_gradientP1NC(champ, domaine_VEF, domaine_Cl_VEF, gij);
  return gij;
 
}
 
DoubleTab& Champ_P1NC::calcul_duidxj_paroi(DoubleTab& tab_gij, const DoubleTab& tab_nu, const DoubleTab& tab_nu_turb, const DoubleTab& tab_tau_tan, const Domaine_Cl_VEF& domaine_Cl_VEF)
{
  const Domaine_VEF& domaine_VEF = domaine_Cl_VEF.domaine_vef();
 
  // Modify Sij for Dirichlet faces (Paroi_fixe) if Paroi_negligeable is not used!
  // We will have: (grad u)_f = (grad u) - (grad u . n) x n + (grad u . n)_lp x n
  //
 
  // PQ : 25/01/08 : Method cleaned up to avoid being confused by sign issues
  //                     and to know exactly what is being computed.
  //
  // PRINCIPLE:
  // (described in 2D; the extension to 3D poses no particular difficulty)
  //
  // Consider 2 frames: the global frame (x,y) and the local frame (t,n) where n is the wall normal.
  // The friction returned by the wall law corresponds to: || tau_tan || = u*^2 = nu.d(u_t)/dn
  // Let P be the change-of-basis matrix from one frame to the other:
  //
  //   P = ( tx  nx )
  //           ( ty  ny )
  //
  // The velocity gradient matrices in each frame:
  //
  //  G_(x,y) =  ( du/dx  du/dy )    and   F_(t,n) = ( du_t/dt   du_t/dn )
  //                 ( dv/dx  dv/dy )                    ( du_n/dt   du_n/dn )
  //
  //                                                                     -1
  // are related by:   G_(x,y)  =  P . F_(t,n) . P
  //
  //
  //  Thus the correction applied in F = ( du_t/dt   du_t/dn + C  )  where  C = -du_t/dn + tau_tan/nu
  //                                            ( du_n/dt   du_n/dn             )
  //
  //  maps into G as follows:
  //
  //                                              -1
  //  G*_(x,y) = G_(x,y) +  P . ( 0  C ). P   =  G_(x,y) + ( C.nx.tx  C.ny.tx )
  //                                ( 0  0 )                   ( C.nx.ty  C.ny.ty )
  //
  //
  // it is this last term that must be determined here
 
  const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
  int nb_cl = les_cl.size();
 
  // tab_tau_tan could be uninitialized at this point. In such a case,
  // we can skip the function because the value of gij won't change
  if (tab_tau_tan.dimension_tot(0) == 0) return tab_gij;
 
  for (int num_cl = 0; num_cl < nb_cl; num_cl++)
    {
      // Loop over the boundaries
      const Cond_lim& la_cl = les_cl[num_cl];
      if (sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) || sub_type(Dirichlet_paroi_defilante, la_cl.valeur()) || la_cl->que_suis_je() == "Frontiere_ouverte_vitesse_imposee_ALE")
        {
          const Front_VF& la_front_dis = ref_cast(Front_VF, la_cl->frontiere_dis());
          int ndeb = la_front_dis.num_premiere_face();
          int nfin = ndeb + la_front_dis.nb_faces();
          int dim = Objet_U::dimension;
          // Loop over the faces
          CDoubleTabView face_normale = domaine_VEF.face_normales().view_ro();
          CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
          CDoubleArrView porosite_face = domaine_Cl_VEF.equation().milieu().porosite_face().view_ro();
          CDoubleTabView tau_tan = tab_tau_tan.view_ro();
          CDoubleArrView nu = static_cast<const ArrOfDouble&>(tab_nu).view_ro();
          CDoubleArrView nu_turb = static_cast<const ArrOfDouble&>(tab_nu_turb).view_ro();
          DoubleTabView3 gij = tab_gij.view_rw<3>();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA (const int fac)
          {
            double P[3][3];
            int num1 = face_voisins(fac, 0);
            // definition of the unit vectors forming the local frame,
            // stored in the change-of-basis matrix P
            // tangential vector (aligned with the tangential velocity)
            double sum = 0.;
            for (int i = 0; i < dim; i++)
              sum += tau_tan(fac, i) * tau_tan(fac, i);
            double norme_tau_tan = sqrt(sum);
            for (int i = 0; i < dim; i++)
              P[i][0] = tau_tan(fac, i) / (norme_tau_tan + DMINFLOAT);
 
            // vector normal to the wall
            sum = 0.;
            for (int i = 0; i < dim; i++)
              sum += face_normale(fac, i) * face_normale(fac, i);
            double norme = sqrt(sum);
 
            int signe = -oriente_normale(fac, num1, face_voisins); // inward orientation
            for (int i = 0; i < dim; i++)
              P[i][1] = signe * face_normale(fac, i) / norme;
 
            // (3D) complete the basis with the second tangential vector
            if (dim == 3)
              {
                P[0][2] = P[1][0] * P[2][1] - P[2][0] * P[1][1];
                P[1][2] = P[2][0] * P[0][1] - P[0][0] * P[2][1];
                P[2][2] = P[0][0] * P[1][1] - P[1][0] * P[0][1];
              }
            //         determine the term d(u_t)/dn to be removed
            //                                                       -1
            //         term identified using the product: F =  P . G . P
            //
            double dutdn_old = 0.;
            for (int i = 0; i < dim; i++)
              for (int j = 0; j < dim; j++)
                {
                  double gij_value = Kokkos::atomic_fetch_add(&gij(num1, i, j), 0.0);
                  dutdn_old += gij_value * P[j][1] * P[i][0];
                }
 
            //         Final correction applied to matrix G
            double C = -dutdn_old + norme_tau_tan / (nu[num1] + nu_turb[num1]) * porosite_face(fac);
 
            // the division by (nu[num1]+nu_turb[num1]) is required because the diffusion operator
            // involves the product: (nu[num1]+nu_turb[num1])*g(i,j)
            for (int i = 0; i < dim; i++)
              for (int j = 0; j < dim; j++)
                Kokkos::atomic_add(&gij(num1, i, j), C * P[j][1] * P[i][0]);
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
 
  return tab_gij;
}
 
////////////////////
// Compute 2SijSij
////////////////////
DoubleVect& Champ_P1NC::calcul_S_barre(const DoubleTab& la_vitesse, DoubleVect& SMA_barre, const Domaine_Cl_VEF& domaine_Cl_VEF)
{
  const Domaine_VEF& domaine_VEF = domaine_Cl_VEF.domaine_vef();
  const int nb_elem = domaine_VEF.nb_elem();
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
 
  DoubleTrav duidxj(nb_elem_tot, dimension, dimension);
 
  Champ_P1NC::calcul_gradient(la_vitesse, duidxj, domaine_Cl_VEF);
 
  int dim = Objet_U::dimension;
  CDoubleTabView3 duidxj_v = duidxj.view_ro<3>();
  DoubleArrView SMA_barre_v = SMA_barre.view_wo();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem, KOKKOS_LAMBDA(
                         const int elem)
  {
    double temp = 0.;
    for (int i = 0; i < dim; i++)
      for (int j = 0; j < dim; j++)
        {
          double Sij = 0.5 * (duidxj_v(elem, i, j) + duidxj_v(elem, j, i));
          temp += Sij * Sij;
        }
    SMA_barre_v(elem) = 2. * temp;
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  return SMA_barre;
}
 
double Champ_P1NC::calculer_integrale_volumique(const Domaine_VEF& domaine, const DoubleVect& v, Ok_Perio ok)
{
  if (ok != FAUX_EN_PERIO)
    {
      // BM: sorry
      Cerr << "Champ_P1NC::calculer_integrale_volumique not yet correctly implemented in periodic case!" << finl;
      exit();
    }
 
  assert(v.line_size() == 1);
  const DoubleVect& volume = domaine.volumes_entrelaces();
  const double resu = mp_prodscal(v, volume);
  return resu;
}
 
int Champ_P1NC::fixer_nb_valeurs_nodales(int nb_noeuds)
{
  assert(nb_noeuds == domaine_vef().nb_faces());
  const MD_Vector& md = domaine_vef().md_vector_faces();
  creer_tableau_distribue(md);
 
  Champ_P1NC_implementation::fixer_nb_valeurs_nodales(nb_noeuds);
 
  return nb_noeuds;
}