Champ_P1iP1B_implementation.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 <Champ_P1iP1B_implementation.h>
#include <Champ_implementation_P1.h>
#include <Domaine.h>
#include <Domaine_VEF.h>
#include <Domaine_Cl_VEF.h>
#include <Debog.h>
 
#include <Matrice_Morse_Sym.h>
#include <TRUSTLists.h>
#include <SolveurSys.h>
#include <Solv_GCP.h>
#include <SSOR.h>
 
#include <Neumann_sortie_libre.h>
#include <Champ_P0_VEF.h>
#include <TRUSTTab_parts.h>
#include <Frontiere_dis_base.h>
 
#include <kokkos++.h>
#include <TRUSTArray_kokkos.tpp>
 
DoubleVect& Champ_P1iP1B_implementation::valeur_a_elem(const DoubleVect& position, DoubleVect& val, int le_poly) const
{
  int dimension=Objet_U::dimension;
  DoubleTab vals(1,1);
  IntVect les_polys(1);
  DoubleTab positions(1,dimension);
 
  les_polys(0) = le_poly;
  for (int i=0; i<dimension; i++)
    positions(0,i) = position(i);
 
  valeur_aux_elems(positions, les_polys, vals);
 
  val(0) = vals(0,0);
  return val;
}
 
double Champ_P1iP1B_implementation::valeur_a_elem_compo(const DoubleVect& position, int le_poly, int ncomp) const
{
  double val=0;
  int dimension=Objet_U::dimension;
 
  DoubleVect vals(1);
  IntVect les_polys(1);
  DoubleTab positions(1,dimension);
 
  les_polys(0) = le_poly;
  for (int i=0; i<dimension; i++)
    positions(0,i) = position(i);
 
  valeur_aux_elems_compo(positions, les_polys, vals, ncomp);
 
  val = vals(0);
 
  return val;
}
 
// Retrieves a positions array containing coordinates of points
// contained in the elements stored in the les_polys array
// Returns the val array containing the field values at those points
DoubleTab& Champ_P1iP1B_implementation::valeur_aux_elems(const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val) const
{
  // zs is properly initialized in 3D
 
  const Domaine_VEF& zvef = domaine_vef();
  const Domaine& domaine_geom = zvef.domaine();
  const DoubleTab& coord = domaine_geom.coord_sommets();
  const IntTab& sommet_poly = domaine_geom.les_elems();
  int prs=zvef.numero_premier_sommet();
  int nb_compo_=le_champ().nb_comp();
  int dimension=Objet_U::dimension;
  const int les_polys_size=les_polys.size();
 
  // TODO : FIXME
  // For FT the resize should be done in its good position and not here ...
  if (val.nb_dim() == 1 ) val.resize(les_polys_size,1);
 
  assert((val.dimension(0) == les_polys.size())||(val.dimension_tot(0) == les_polys.size()));
  assert(val.line_size() == nb_compo_);
 
  // Filter the field into the champ_filtre_ array
  champ_filtre_=filtrage(zvef,le_champ());
 
  if (nb_compo_ == 1)
    {
      // HIP Kokkos prevent from the function calls inside the Kokkos region:
      // Error: error: reference to __host__ function 'get_alphaE' in __host__ __device__ function
      bool alphaE = zvef.get_alphaE();
      bool alphaS = zvef.get_alphaS();
      CDoubleArrView champ_filtre_v = static_cast<DoubleVect&>(champ_filtre_).view_ro();
      CIntTabView sommet_poly_v = sommet_poly.view_ro();
      CDoubleTabView coord_v = coord.view_ro();
      CDoubleTabView positions_v = positions.view_ro();
      CIntArrView les_polys_v = les_polys.view_ro();
      DoubleTabView val_v = val.view_rw();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),les_polys_size, KOKKOS_LAMBDA(
                             const int rang_poly)
      {
        // Initialize
        val_v(rang_poly,0) = 0.;
 
        // Take le_poly
        int le_poly=les_polys_v(rang_poly);
        if (le_poly != -1)
          {
            // Contribution P0
            if (alphaE)
              val_v(rang_poly, 0) += champ_filtre_v(le_poly);
            // Contribution P1
            if (alphaS)
              {
                double xs = positions_v(rang_poly,0);
                double ys = positions_v(rang_poly,1);
                double zs = (dimension == 3) ? positions_v(rang_poly,2) : 0;
 
                // Compute the filtered field at point (xs,ys,zs) using P1 shape functions
                for (int i=0; i<dimension+1; i++)
                  {
                    int som = prs+sommet_poly_v(le_poly,i);
                    double champ_filtre_som = champ_filtre_v(som);
                    double li=0.;
 
                    if (dimension == 2)
                      li=coord_barycentrique_P1_triangle(sommet_poly_v,
                                                         coord_v, xs, ys,
                                                         le_poly, i);
                    else if (dimension == 3)
                      li=coord_barycentrique_P1_tetraedre(sommet_poly_v,
                                                          coord_v, xs, ys, zs,
                                                          le_poly, i);
                    val_v(rang_poly,0) += champ_filtre_som * li;
                  }
              }
          }
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
  else // nb_compo_ > 1
    {
      Cerr << "Are you really here already?" << finl;
      Process::exit();
    }
  return val;
}
 
// Same as valeur_aux_elems but for one component of a vector field
DoubleVect& Champ_P1iP1B_implementation::valeur_aux_elems_compo(const DoubleTab& positions, const IntVect& les_polys, DoubleVect& val, int ncomp) const
{
  Cerr << "Champ_P1iP1B_implementation::valeur_aux_elems_compo not implemented." << finl;
  // Factorize code with the valeur_aux_elems() method
  Process::exit();
  return val;
}
 
 
// Retrieves a domain
// Returns an array containing the field values at the domain vertices
DoubleTab& Champ_P1iP1B_implementation::valeur_aux_sommets(const Domaine& dom, DoubleTab& tab_val) const
{
  const Domaine_VEF& zvef = domaine_vef();
  int nb_compo_=le_champ().nb_comp();
  assert(nb_compo_ == tab_val.line_size());
 
  // Filter the field for post-processing (stored in the champ_filtre_ array)
  champ_filtre_=filtrage(zvef,le_champ());
 
  if (zvef.get_alphaE()) // Support P0
    {
      Champ_P0_VEF tmp;
      tmp.associer_domaine_dis_base(zvef);
      tmp.valeurs() = champ_filtre_;
      tmp.valeur_aux_sommets(zvef.domaine(), tab_val);
    }
 
  if (zvef.get_alphaS()) // Support P1
    {
      int nbs=zvef.nb_som();
      int prs=zvef.numero_premier_sommet();
      if (nb_compo_ == 1)
        {
          CDoubleArrView champ_filtre = static_cast<const DoubleVect&>(champ_filtre_).view_ro();
          DoubleTabView val = tab_val.view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nbs, KOKKOS_LAMBDA(const int num_som)
          {
            int som = prs+num_som;
            val(num_som,0) += champ_filtre(som);
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
      else // nb_compo_ > 1
        {
          Cerr << "Are you really here already?" << finl;
          Process::exit();
        }
    }
  return tab_val;
}
 
 
//
DoubleVect& Champ_P1iP1B_implementation::
valeur_aux_sommets_compo(const Domaine& dom,
                         DoubleVect& val,
                         int ncomp) const
{
  Cerr << "Champ_P1iP1B_implementation::valeur_aux_sommets_compo not implemented." << finl;
  // Factorize with the valeur_aux_sommets() method
  Process::exit();
  return val;
}
 
DoubleTab& Champ_P1iP1B_implementation::remplir_coord_noeuds(DoubleTab& coord) const
{
  const Domaine_VEF& zvef=domaine_vef();
  const Domaine& le_dom=zvef.domaine();
  const Domaine& dom=le_dom;
  const DoubleTab& coord_sommets=dom.coord_sommets();
  int dimension=Objet_U::dimension;
 
  const DoubleTab& xg = zvef.xp();
  int nbe=zvef.get_alphaE() ? zvef.nb_elem_tot() : 0;
  int nbs=zvef.get_alphaS() ? zvef.nb_som_tot() : 0;
  coord.resize(nbe+nbs,dimension);
  int k = 0;
  if (zvef.get_alphaE())
    {
      for(int i=0; i<nbe; i++)
        {
          for(int j=0; j<dimension; j++)
            coord(k,j)=xg(i,j);
          k++;
        }
    }
  if (zvef.get_alphaS())
    {
      for(int i=0; i<nbs; i++)
        {
          for(int j=0; j<dimension; j++)
            coord(k,j)=coord_sommets(i,j);
          k++;
        }
    }
  return coord;
}
 
int Champ_P1iP1B_implementation::
remplir_coord_noeuds_et_polys(DoubleTab& positions,
                              IntVect& polys) const
{
  Cerr << "Not implemented" << finl;
  assert(0);
  Process::exit();
  return 1;
}
 
void assembler(const Domaine_VEF& domaine_VEF, Matrice& matrice)
{
  matrice.typer("Matrice_Morse_Sym");
  Matrice_Morse_Sym& MatPoisson=ref_cast(Matrice_Morse_Sym, matrice.valeur());
  int nb_som_tot = domaine_VEF.domaine().nb_som_tot();
  int nb_arete_tot = domaine_VEF.domaine().nb_aretes_tot();
  int nnz=nb_som_tot+nb_arete_tot;
  const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
  const ArrOfInt& renum_arete_perio=domaine_VEF.get_renum_arete_perio();
  const Domaine& dom=domaine_VEF.domaine();
 
  IntLists voisins(nb_som_tot);
  DoubleLists coeffs(nb_som_tot);
  DoubleVect diag(nb_som_tot);
  // Loop over all non-periodic edges:
  for(int arete=0; arete<nb_arete_tot; arete++)
    {
      if (renum_arete_perio[arete]==arete)
        {
          int som1=dom.get_renum_som_perio(aretes_som(arete, 0));
          int som2=dom.get_renum_som_perio(aretes_som(arete, 1));
          if(som1>som2)
            {
              int tmp=som1;
              som1=som2;
              som2=tmp;
            }
          voisins[som1].add(som2);
          coeffs[som1].add(1);
          diag(som1)++;
          diag(som2)++;
        }
    }
  for (int i=0; i<nb_som_tot; i++)
    {
      if (diag(i)==0)
        {
          assert(i!=dom.get_renum_som_perio(i));
          diag(i)=1; // Periodic vertices
        }
    }
  MatPoisson.dimensionner(nb_som_tot, nnz) ;
  MatPoisson.remplir(voisins, coeffs, diag) ;
  MatPoisson.compacte() ;
  MatPoisson.set_est_definie(1);
  //Cerr << "Matrice P1";MatPoisson.imprimer(Cerr);
}
 
static double coeff = 4. / 5.;
static double coeff_inv = 1 / coeff;
double second_membre(const Domaine_VEF& domaine_VEF, ArrOfDouble& Pa, DoubleVect& secmem)
{
  int nb_arete_tot = domaine_VEF.domaine().nb_aretes_tot();
  const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
  const ArrOfInt& renum_arete_perio=domaine_VEF.get_renum_arete_perio();
  const Domaine& dom=domaine_VEF.domaine();
  secmem=0;
  // Loop over all non-periodic edges
  for(int arete=0; arete<nb_arete_tot; arete++)
    {
      if (renum_arete_perio[arete]==arete)
        {
          double x=Pa[arete];
          int som1=dom.get_renum_som_perio(aretes_som(arete, 0));
          int som2=dom.get_renum_som_perio(aretes_som(arete, 1));
          secmem(som1)+=x;
          secmem(som2)+=x;
        }
    }
  secmem*=coeff;
  secmem.echange_espace_virtuel();
  double norme=mp_norme_vect(secmem);
  //Cerr << "||somme Pa||= " << norme << finl;
  return norme;
}
 
void corriger(const Domaine_VEF& domaine_VEF, DoubleTab& champ_filtre_, Matrice& matrice, const int Condition_Neumann_imposee_)
{
  int nb_elem=domaine_VEF.nb_elem();
  int nb_som=domaine_VEF.nb_som();
 
  // Access arrays for Pk and Ps values
  if (!domaine_VEF.get_alphaE() || !domaine_VEF.get_alphaS())
    {
      Cerr << "Error in Champ_P1iP1B_impl.cpp: corriger(...), the field must have both vertex and element parts" << finl;
      Process::exit();
    }
  DoubleTab_parts parties_P(champ_filtre_);
  DoubleVect& Pk = parties_P[0];  // element part
  DoubleVect& Ps = parties_P[1];  // vertex part
 
  // Filtering if edge support is active
  if (domaine_VEF.get_alphaA())
    {
      if (domaine_VEF.get_renum_arete_perio().size_array()==0)
        {
          Cerr<<"We try to Champ_P1iP1B_impl::corriger but get_renum_arete_perio"<<finl;
        }
    }
  if (domaine_VEF.get_alphaA()&& domaine_VEF.get_renum_arete_perio().size_array())
    {
      champ_filtre_.ensureDataOnHost(); // We need copy H2D as the Pa part of DoubleTab_parts is computed on host
      ToDo_Kokkos("critical for P0P1Pa");
      DoubleVect& Pa = parties_P[2];  // edge part
 
      // If first pass, assemble the matrix
      if (!matrice) assembler(domaine_VEF, matrice);
 
      // Build the right-hand side
      DoubleVect secmem;
      domaine_VEF.domaine().creer_tableau_sommets(secmem);
      second_membre(domaine_VEF, Pa, secmem);
 
      // Compute the correction
      DoubleVect solution;
      solution.copy(secmem, RESIZE_OPTIONS::NOCOPY_NOINIT); // Copy of the structure
      solution = 0.;
 
      SolveurSys solveur;
      solveur.typer("Solv_GCP");
      OWN_PTR(Precond_base) p;
      p.typer("SSOR");
      ref_cast(Solv_GCP,solveur.valeur()).set_precond(p);
      solveur.nommer("Pa_filter_solver");
      solveur.resoudre_systeme(matrice, secmem, solution);
 
      // Apply periodicity to the solution:
      const Domaine& dom=domaine_VEF.domaine();
      for(int i=0; i<nb_som; i++)
        {
          int som=dom.get_renum_som_perio(i);
          if (som!=i)
            solution(i)=solution(som);
        }
 
      // Correction of Ps:
      for(int som=0; som<nb_som; som++)
        Ps(som) += solution(som);
 
      // Correction of Pa:
      int nb_arete = domaine_VEF.domaine().nb_aretes();
      const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
      const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
      for(int arete=0; arete<nb_arete; arete++)
        {
          if(!ok_arete[arete] && Pa(arete)!=0)
            {
              Cerr << "Pa(arete_superflue)!=0 in Champ_P1iP1B_implementation::corriger" << finl;
              Cerr << "Please contact TRUST support." << finl;
              Cerr << "If this is a restart from an xyz format calculation," << finl;
              Cerr << "you must overwrite the .ok_arete file of your restart calculation" << finl;
              Cerr << "with the .ok_arete file of your previous calculation so as to have the" << finl;
              Cerr << "same superfluous edges." << finl;
              Process::exit();
            }
          int som1=aretes_som(arete, 0);
          int som2=aretes_som(arete, 1);
          double sigma=solution(som1)+solution(som2);
          Pa(arete)-=sigma*coeff_inv;
        }
 
      // Correction of Pk:
      const IntTab& les_elems = domaine_VEF.domaine().les_elems();
      int nb_som_par_elem=les_elems.dimension_tot(1);
      for(int elem=0; elem<nb_elem; elem++)
        {
          double sigma=0;
          for (int som=0; som<nb_som_par_elem; som++)
            sigma+=solution(les_elems(elem,som));
          Pk(elem)-=sigma/4;
        }
    }
 
  // L2 filtering:
  double moyenne_K = mp_moyenne_vect(Pk);                // Compute the mean of the element field
  Pk -= moyenne_K;                                 // Correction of elements of champ_filtre_
  Ps += moyenne_K;                                 // Correction of vertices of champ_filtre_
 
  // Subtract the mean from Ps if there is no Neumann BC, to impose a zero mean
  // It might be better to do this in Assembleur_P_VEFPreP1B::modifier_solution
  // if possible...
  assert(Condition_Neumann_imposee_!=-1);
 
  if (!Condition_Neumann_imposee_)
    {
      // The virtual space of Ps does not need to be up to date to correctly compute the mean
      Ps -= mp_moyenne_vect(Ps); // Correction of vertices of champ_filtre_
    }
 
  champ_filtre_.echange_espace_virtuel();         // Update of virtual spaces
}
 
DoubleTab& Champ_P1iP1B_implementation::filtrage(const Domaine_VEF& zvef, const Champ_base& un_champ, bool implicitCoupling) const
{
  // Filter only if both element and vertex supports are present
  if (zvef.get_alphaE() && zvef.get_alphaS())
    {
      // No filtering if it has already been done on this field during the time step
      // and if no sub-iterations are performed during the time step
      // The condition implicitCoupling=true enables filtering and therefore updates
      // the variable "champ_filtre_" during sub-iterations.
      // Sub-iterations can be performed within a time step in the case of an implicit
      // coupling with another code, for example with a structural solver
      // in fluid-structure interaction applications
      if (un_champ.valeurs().data()==adresse_champ_filtre_ // PL: data() and not addr(). Else copyFromDevice on GPU !
          && un_champ.temps()==temps_filtrage_
          && !implicitCoupling)
        return champ_filtre_;
 
      // Copy the field to filter into the array that will hold the filtered field
      champ_filtre_=un_champ.valeurs();
      temps_filtrage_=un_champ.temps();
      adresse_champ_filtre_=un_champ.valeurs().data();
      //Cout << "Filtering the field " << un_champ.le_nom() << " at time " << un_champ.temps() << finl;
 
      // Correct the field
      corriger(zvef, champ_filtre_, matrice_filtrage_, Condition_Neumann_imposee_);
      Debog::verifier("field after filter=", champ_filtre_);
      return champ_filtre_;
    }
  else
    {
      champ_filtre_=un_champ.valeurs();
      return champ_filtre_;
    }
}
 
// Sets Condition_Neumann_imposee_ to 1 if a Neumann BC exists, 0 otherwise
// If there is no Neumann condition, a zero mean will be imposed during field filtering
// to get the same field in sequential and parallel runs
void Champ_P1iP1B_implementation::completer(const Domaine_Cl_dis_base& zcl)
{
  const Conds_lim& les_cl = zcl.les_conditions_limites();
  Condition_Neumann_imposee_ = 0;
  for(int i=0; i<les_cl.size(); i++)
    if (sub_type(Neumann_sortie_libre,les_cl[i].valeur()))
      Condition_Neumann_imposee_ = 1;
}
 
DoubleTab& Champ_P1iP1B_implementation::trace(const Frontiere_dis_base& fr, const DoubleTab& y, DoubleTab& x,int distant) const
{
  assert(distant==0);
  const Front_VF& fr_vf = ref_cast(Front_VF, fr);
  const Domaine_VEF& domaine_dis = domaine_vef();
  const DoubleTab& cg_faces = domaine_dis.xv();
  int nb_faces_fr = fr_vf.nb_faces();
  int face, elem;
 
  assert(x.dimension(0)==fr_vf.nb_faces());
 
  IntVect elem_voisins(nb_faces_fr);
  elem_voisins = -1;
  DoubleTab cg_faces_fr(nb_faces_fr,Objet_U::dimension);
  DoubleTab val_interp(nb_faces_fr, 1);
 
  // Retrieve the neighboring elements of the boundary faces
  for (int i=0; i<nb_faces_fr; i++)
    {
      face = fr_vf.num_premiere_face()+i;
      elem = domaine_dis.face_voisins(face,0);
      if (elem==-1)
        elem = domaine_dis.face_voisins(face,1);
      elem_voisins(i) = elem;
    }
 
  // Fill the array containing the coordinates
  // of the centroid of the boundary faces
  for (int num_face=0; num_face<nb_faces_fr; num_face++)
    for (int dim=0; dim<Objet_U::dimension; dim++)
      {
        face = fr_vf.num_premiere_face()+num_face;
        cg_faces_fr(num_face,dim) = cg_faces(face,dim);
      }
 
  // Interpolate field values at the centroid of boundary faces
  valeur_aux_elems(cg_faces_fr,elem_voisins,val_interp);
 
  for (int i=0; i<nb_faces_fr; i++)
    {
      x(i,0) = val_interp(i);
    }
  // Useless ? x.echange_espace_virtuel();
  return x;
}