Op_Conv_EF_VEF_P1NC_Stab.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 <Op_Conv_EF_VEF_P1NC_Stab.h>
#include <Champ_P1NC.h>
#include <BilanQdmVEF.h>
#include <Sous_Domaine.h>
#include <Sous_domaine_dis_base.h>
#include <Schema_Temps_base.h>
#include <Debog.h>
#include <Porosites_champ.h>
#include <Sous_domaine_VF.h>
#include <Probleme_base.h>
#include <ArrOfBit.h>
#include <TRUSTVects.h>
#include <SFichier.h>
#include <TRUSTList.h>
#include <TRUSTTabs.h>
#include <Neumann_sortie_libre.h>
#include <Dirichlet_homogene.h>
#include <Neumann_homogene.h>
#include <Periodique.h>
#include <Symetrie.h>
#include <Echange_impose_base.h>
#include <Array_tools.h>
 
Implemente_instanciable(Op_Conv_EF_VEF_P1NC_Stab,"Op_Conv_EF_Stab_VEF_P1NC",Op_Conv_VEF_Face);
 
// XD sous_zone_valeur objet_lecture nul NO_BRACE Two words.
// XD attr sous_zone ref_sous_zone sous_zone REQ sous zone
// XD attr valeur floattant valeur REQ value
 
// XD listsous_zone_valeur listobj nul NO_BRACE sous_zone_valeur NO_COMMA List of groups of two words.
 
// XD convection_ef_stab convection_deriv ef_stab BRACE Keyword for a VEF convective scheme.
// XD attr alpha floattant alpha OPT To weight the scheme centering with the factor double (between 0 (full centered)
// XD_CONT and 1 (mix between upwind and centered), by default 1). For scalar equation, it is adviced to use alpha=1 and
// XD_CONT for the momentum equation, alpha=0.2 is adviced.
// XD attr test entier test OPT Developer option to compare old and new version of EF_stab
// XD attr tdivu rien tdivu OPT To have the convective operator calculated as div(TU)-TdivU(=UgradT).
// XD attr old rien old OPT To use old version of EF_stab scheme (default no).
// XD attr volumes_etendus rien volumes_etendus OPT Option for the scheme to use the extended volumes (default, yes).
// XD attr volumes_non_etendus rien volumes_non_etendus OPT Option for the scheme to not use the extended volumes
// XD_CONT (default, no).
// XD attr amont_sous_zone ref_sous_zone amont_sous_zone OPT Option to degenerate EF_stab scheme into Amont (upwind)
// XD_CONT scheme in the sub zone of name sz_name. The sub zone may be located arbitrarily in the domain but the more
// XD_CONT often this option will be activated in a zone where EF_stab scheme generates instabilities as for free outlet
// XD_CONT for example.
// XD attr alpha_sous_zone listsous_zone_valeur alpha_sous_zone OPT Option to change locally the alpha value on N
// XD_CONT sub-zones named sub_zone_name_I. Generally, it is used to prevent from a local divergence by increasing
// XD_CONT locally the alpha parameter.
 
Sortie& Op_Conv_EF_VEF_P1NC_Stab::printOn(Sortie& s ) const
{
  return s << que_suis_je();
}
 
Entree& Op_Conv_EF_VEF_P1NC_Stab::readOn(Entree& s )
{
  //Les mots a reconnaitre
  Motcle motlu, accouverte = "{" , accfermee = "}" ;
  Motcles les_mots(11);
  {
    les_mots[0] = "alpha";
    les_mots[1] = "test";
    les_mots[2] = "TdivU";
    les_mots[3] = "old";
    les_mots[4] = "volumes_etendus";
    les_mots[5] = "volumes_non_etendus";
    les_mots[6] = "amont_sous_domaine";
    les_mots[7] = "amont_sous_zone";
    les_mots[8] = "nouvelle_matrice_implicite";
    les_mots[9] = "alpha_sous_domaine";
    les_mots[10] = "alpha_sous_zone";
  }
 
  s >> motlu;
  if (motlu!=accouverte)
    {
      Cerr << "Error Op_Conv_EF_VEF_P1NC_Stab::readOn()" << finl;
      Cerr << "Since version 1.5.3, the syntax of the keyword EF_stab has changed." << finl;
      Cerr << "It must begin with an opening brace {" << finl;
      Cerr << "and the optional parameters are between the braces:" << finl;
      Cerr << "Convection { EF_stab } -> Convection { EF_stab { } }" << finl;
      Process::exit();
    }
  s >> motlu;
 
  while(motlu!=accfermee)
    {
      int rang=les_mots.search(motlu);
 
      switch(rang)
        {
        case 0 :
 
          s >> alpha_;
          break;
 
        case 1 :
 
          test_=1;
          break;
 
        case 2 :
 
          is_compressible_=1;
          break;
 
        case 3 :
 
          old_=1;
          break;
 
        case 4 :
 
          volumes_etendus_=1;
          break;
 
        case 5 :
 
          volumes_etendus_=0;
          break;
 
        case 6 :
        case 7 :
          sous_domaine=true;
          s >> nom_sous_domaine;
          break;
        case 8 :
          s >> new_jacobienne_;
          break;
        case 9 :
        case 10 :
          ssz_alpha=true;
          s >> nb_ssz_alpha;
          noms_ssz_alpha.dimensionner(nb_ssz_alpha);
          alpha_ssz.resize(nb_ssz_alpha);
          for (int i=0; i<nb_ssz_alpha; i++)
            {
              s>>noms_ssz_alpha[i];
              s>>alpha_ssz(i);
            }
          break;
        default :
 
          Cerr << "Error Op_Conv_EF_VEF_P1NC_Stab::readOn()" << finl;
          Cerr << "Keyword " << motlu << " not recognized." << finl;
          Cerr << "Exiting." << finl;
          Process::exit();
 
        }//end switch
 
      //Continue reading
      s >> motlu;
 
    }//end while
 
  return s ;
}
 
 
static KOKKOS_INLINE_FUNCTION double maximum(const double x,
                                             const double y)
{
  return x<y ? y : x;
}
 
static KOKKOS_INLINE_FUNCTION double maximum(const double x,
                                             const double y,
                                             const double z)
{
  return maximum(maximum(x,y),z);
}
 
static KOKKOS_INLINE_FUNCTION double minimum(const double x,
                                             const double y)
{
  return x>y ? y : x;
}
 
static KOKKOS_INLINE_FUNCTION double Dij(int elem,
                                         int face_loc_i,
                                         int face_loc_j,
                                         CDoubleTabView3 Kij)
{
  const double kij=Kij(elem,face_loc_i,face_loc_j);
  const double kji=Kij(elem,face_loc_j,face_loc_i);
  return maximum(-kij,-kji,0);
}
 
static inline double Dij(int elem,
                         int face_loc_i,
                         int face_loc_j,
                         const DoubleTab& Kij)
{
  const double kij=Kij(elem,face_loc_i,face_loc_j);
  const double kji=Kij(elem,face_loc_j,face_loc_i);
  return maximum(-kij,-kji,0);
}
 
static KOKKOS_INLINE_FUNCTION double limiteur(double r)
{
  return r<=0 ? 0 : Kokkos::fmax(Kokkos::fmin(2,r),Kokkos::fmin(1,2*r));//SuperBee
}
 
KOKKOS_INLINE_FUNCTION double formule_Id_2D(int n)
{
  return 1./3;
}
 
KOKKOS_INLINE_FUNCTION double formule_Id_3D(int n)
{
  return 0.25;
}
 
KOKKOS_INLINE_FUNCTION double formule_2D(int n)
{
  switch(n)
    {
    case 0:
      return 1./3;
    case 1:
      return 0.5;
    case 2:
      return 1.;
    default:
      Process::Kokkos_exit("Erreur Op_Conv_EF_VEF_P1NC_Stab::formule_2D() ");
    }
  return 0.;
}
 
KOKKOS_INLINE_FUNCTION double formule_3D(int n)
{
  switch(n)
    {
    case 0:
      return 0.25;
    case 1:
      return 1./3;
    case 2:
      return 0.5;
    case 3:
      return 1.;
    default:
      Process::Kokkos_exit("Erreur Op_Conv_EF_VEF_P1NC_Stab::formule_3D() ");
    }
  return 0.;
}
 
 
 
////////////////////////////////////////////////////////////////////
//
//                      Implementation of functions
//
//                   of class Op_Conv_EF_VEF_P1NC_Stab
//
////////////////////////////////////////////////////////////////////
 
void Op_Conv_EF_VEF_P1NC_Stab::reinit_conv_pour_Cl(const DoubleTab& transporte,const IntList& faces, const DoubleTabs& valeurs_faces, const DoubleTab& tab_vitesse, DoubleTab& resu) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF=la_zcl_vef.valeur();
  const DoubleTab& face_normales=domaine_VEF.face_normales();
  const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
  const int nb_comp=transporte.line_size();
  int n_bord=0, num1=0, num2=0, ind_face=0,facei=0, dim=0;
  double psc=0.;
 
  const DoubleVect& transporteV= transporte;
 
  //Loop to dimension the "faces" list
  for (n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      num1 = 0;
      num2 = le_bord.nb_faces_tot();
 
      if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
        {
          const Neumann_sortie_libre& la_sortie_libre
            = ref_cast(Neumann_sortie_libre, la_cl.valeur());
          ToDo_Kokkos("critical");
          for (ind_face=num1; ind_face<num2; ind_face++)
            {
              facei = le_bord.num_face(ind_face);
 
              psc=0.;
              for (dim=0; dim<dimension; dim++)
                psc-=tab_vitesse(facei,dim)*face_normales(facei,dim);
 
              if (psc>0)
                for (dim=0; dim<nb_comp; dim++)
                  resu(facei,dim)+=psc*(la_sortie_libre.val_ext(facei-num1,dim)-transporteV[facei*nb_comp+dim]);
 
            }
        }//end if on "Neumann_sortie_libre"
    }
}
 
void Op_Conv_EF_VEF_P1NC_Stab::calculer_coefficients_operateur_centre(DoubleTab& tab_Kij,const int nb_comp, const DoubleTab& tab_vitesse) const
{
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const int nb_faces_elem = domaine_VEF.elem_faces().line_size();
  int dim = Objet_U::dimension;
 
  assert(tab_Kij.nb_dim()==3);
  assert(tab_Kij.dimension(0)==nb_elem_tot);
  assert(tab_Kij.dimension(1)==tab_Kij.dimension(2));
  assert(tab_Kij.dimension(1)==nb_faces_elem);
 
  //
  //Compute the operator coefficients
  //
  CDoubleTabView face_normales = 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 vitesse = tab_vitesse.view_ro();
  DoubleTabView3 Kij = tab_Kij.view_rw<3>();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
                         const int elem, const int face_loci)
  {
    int face_i=elem_faces(elem,face_loci);
 
    double signei=1.0;
    if(face_voisins(face_i,0)!=elem) signei=-1.0;
 
    double psci=0;
    for(int comp=0; comp<dim; comp++)
      psci+=vitesse(face_i,comp)*face_normales(face_i,comp);
    psci*=signei;
 
    //Kij(elem,face_loci,face_loci)=0.;
    for(int face_locj=face_loci+1; face_locj<nb_faces_elem; face_locj++)
      {
        int face_j=elem_faces(elem,face_locj);
        double signej=1.0;
        if(face_voisins(face_j,0)!=elem)
          signej=-1.0;
 
        double pscj=0;
        //psci=0;
        for(int comp=0; comp<dim; comp++)
          pscj+=vitesse(face_j,comp)*face_normales(face_j,comp);
        pscj*=signej;
 
        Kokkos::atomic_add(&Kij(elem,face_loci,face_locj), -1./nb_faces_elem*pscj);
        Kokkos::atomic_add(&Kij(elem,face_loci,face_loci), +1./nb_faces_elem*pscj);
        Kokkos::atomic_add(&Kij(elem,face_locj,face_loci), -1./nb_faces_elem*psci);
        Kokkos::atomic_add(&Kij(elem,face_locj,face_locj), +1./nb_faces_elem*psci);
      }//end for on "face_locj"
  });//end for on "elem"
  end_gpu_timer(__KERNEL_NAME__);
  //
  // Correction of Kij for Dirichlet!
  //
  {
    int nb_bord=domaine_Cl_VEF.nb_cond_lim();
    for (int n_bord=0; n_bord<nb_bord; n_bord++)
      {
        const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
 
        if ( (sub_type(Dirichlet,la_cl.valeur()))
             || (sub_type(Dirichlet_homogene,la_cl.valeur()))
           )
          {
            // GF: removing the test on nb_comp because otherwise in scalar mode we always use extended volume
            if ( volumes_etendus_ )
              {
                CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
                CIntTabView elem_faces_dirichlet_v = elem_faces_dirichlet_.view_ro();
                CIntArrView elem_nb_faces_dirichlet_v = elem_nb_faces_dirichlet_.view_ro();
                CIntArrView elem_faces_frontiere_v = elem_faces_frontiere[n_bord].view_ro();
                const int elem_faces_frontiere_size = elem_faces_frontiere[n_bord].size_array();
                //
                //Modification of the matrix coefficients
                //
                Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                                     Kokkos::RangePolicy<>(0,elem_faces_frontiere_size), KOKKOS_LAMBDA(
                                       const int elem_ind)
                {
                  const int elem=elem_faces_frontiere_v(elem_ind);
                  assert(elem!=-1);
                  const int nb_faces_bord = elem_nb_faces_dirichlet_v(elem);
 
                  //
                  //Compute the weighting coefficient
                  //
 
                  const double coeff = (dim==2) ?
                                       nb_faces_bord/2. : nb_faces_bord*nb_faces_bord/6.-nb_faces_bord/3.+1./2;
 
                  //
                  //End of weighting coefficient computation
                  //
 
                  for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                    {
                      const int face_i = elem_faces(elem,face_loc_i);
 
                      /* Determine whether "face_i" is a Dirichlet face */
                      bool is_not_on_boundary = true;
                      for (int f_loc=0; f_loc<nb_faces_bord; f_loc++)
                        is_not_on_boundary&=(face_i!=elem_faces_dirichlet_v(elem,f_loc));
 
                      if (is_not_on_boundary) /* if "face_i" is not a Dirichlet face */
                        {
                          for (int face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
                            {
                              /* Modify the matrix coefficients */
                              /* according to the extended shape functions */
                              for (int f_loc=0; f_loc<nb_faces_bord; f_loc++)
                                {
                                  const int face_bord = elem_faces_dirichlet_v(elem,f_loc);
                                  const int face_loc_k = num_fac_loc(face_bord,0);
                                  assert(face_loc_k>=0);
                                  assert(face_loc_k<nb_faces_elem);
 
                                  const double kkj = Kij(elem,face_loc_k,face_loc_j);
                                  Kij(elem,face_loc_i,face_loc_j) += coeff*kkj;
                                }//end for on "f_loc"
 
                            }//end of for on "face_loc_k"
 
                        }//end of if
 
                    }//end of for on "face_loc_i"
 
 
                  //
                  //Reset to zero the coefficients associated with nodes
                  //that lie on the Dirichlet boundary
                  //
                  {
                    for (int f_loc=0; f_loc<nb_faces_bord; f_loc++)
                      {
                        const int face_bord = elem_faces_dirichlet_v(elem,f_loc);
                        const int face_loc_j = num_fac_loc(face_bord,0);
                        assert(face_loc_j>=0);
                        assert(face_loc_j<nb_faces_elem);
 
                        for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                          Kij(elem,face_loc_j,face_loc_i)=0;
                      }//end of for on "f_loc"
                  }
                  //
                  //End of reset to zero
                  //
 
                  //
                  //To recover the LED scheme
                  //
                  {
                    for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                      {
                        double sum=0.;
                        for (int face_loc_k=0; face_loc_k<nb_faces_elem; face_loc_k++)
                          {
                            sum+=Kij(elem,face_loc_i,face_loc_k);
                          }
                        //Cerr << "somme apres : " << sum << finl;
                        Kij(elem,face_loc_i,face_loc_i) -= sum;//car div(u)=0!
                      }
                  }
                  //
                  //End of LED scheme
                  //
 
                });// for face
 
                end_gpu_timer(__KERNEL_NAME__);
 
                //
                //End of matrix coefficient modification
                //
 
              }//end of if on "volumes_etendus_"
 
          }// sub_type Dirichlet
 
 
        else if (sub_type(Neumann,la_cl.valeur()) || sub_type(Neumann_homogene,la_cl.valeur()) || sub_type(Neumann_val_ext,la_cl.valeur()))
          {
            //Do nothing
          }//end of if on Neumann
 
        else if (sub_type(Symetrie,la_cl.valeur()))
          {
            //Do nothing
          }//end of if on Symetrie
 
        else if (sub_type(Periodique,la_cl.valeur()))
          {
            //Do nothing
          }//end of if on Periodique
 
        else if (sub_type(Echange_impose_base,la_cl.valeur()))
          {
            //Do nothing
          }//end of if on Echange_impose_base
 
        else
          {
            Cerr << "Error Op_Conv_EF_VEF_P1NC_Stab::calculer_coefficients_operateur_centre()" << finl;
            Cerr << "Boundary condition " << la_cl.que_suis_je() << " not implemented."   << finl;
            Cerr << "Exiting." << finl;
            Process::exit();
          }//end of else on other boundary conditions
 
      }//end of boundary conditions
  }
 
  //
  // End of Kij correction
  //
}
void Op_Conv_EF_VEF_P1NC_Stab::remplir_fluent() const
{
  // CFL computation.
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const int nb_faces = domaine_VEF.nb_faces();
  int nb_comp = Objet_U::dimension;
  CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
  CDoubleTabView tab_vitesse = vitesse_->valeurs().view_ro();
  DoubleArrView fluent = fluent_.view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       Kokkos::RangePolicy<>(0, nb_faces), KOKKOS_LAMBDA(
                         const int num_face)
  {
    double psc=0.;
    for (int i=0; i<nb_comp; i++)
      psc+=tab_vitesse(num_face,i)*face_normales(num_face,i);
    fluent(num_face)=std::fabs(psc);
  });
  end_gpu_timer(__KERNEL_NAME__);
}
 
 
DoubleTab& Op_Conv_EF_VEF_P1NC_Stab::ajouter(const DoubleTab& transporte_2,
                                             DoubleTab& resu) const
{
  DoubleTab sauv(resu);
  resu=0;
 
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
  const DoubleTab& vitesse_2=la_vitesse.valeurs();
  const DoubleVect& porosite_face = equation().milieu().porosite_face();
 
  const IntTab& elem_faces = domaine_VEF.elem_faces();
 
  const int marq=phi_u_transportant(equation());
  const int nb_faces_elem=elem_faces.dimension(1);
  assert(nb_faces_elem==(dimension+1));
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  int nb_comp=resu.line_size();
 
  DoubleTab transporte_;
  DoubleTab vitesse_face_;
  // either transporte=phi*transporte_ and vitesse=vitesse_
  // or transporte=transporte_ and vitesse=phi*vitesse_
  // depending on whether transport uses phi*u or u.
  const DoubleTab& tab_vitesse=modif_par_porosite_si_flag(vitesse_2,vitesse_face_,marq,porosite_face);
 
  // either transporte=phi*transporte_ and vitesse=vitesse_
  // or transporte=transporte_ and vitesse=phi*vitesse_
  // depending on whether transport uses phi*u or u.
  const DoubleTab& transporte=modif_par_porosite_si_flag(transporte_2,transporte_,!marq,porosite_face);
 
  DoubleTrav Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
  //Initialize the flux_bords_ array for pressure drop computation
  flux_bords_.resize(domaine_VEF.nb_faces_bord(),nb_comp);
  calculer_flux_bords(Kij,tab_vitesse,transporte);
 
  if (!old_)
    {
      calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
      if (is_compressible_) ajouter_partie_compressible(transporte,resu,tab_vitesse);
 
      ajouter_operateur_centre(Kij,transporte,resu);
 
      ajouter_diffusion(Kij,transporte,resu);
      ajouter_antidiffusion(Kij,transporte,resu);
      mettre_a_jour_pour_periodicite(resu);
 
    }
  else
    {
      assert(old_==1);
      ajouter_old(transporte,resu,tab_vitesse);
    }
 
  //To account for Neumann outflow boundary conditions
  IntList NeumannFaces;
  DoubleTabs ValeursNeumannFaces;
  reinit_conv_pour_Cl(transporte,NeumannFaces,ValeursNeumannFaces,tab_vitesse,resu);
 
  resu+=sauv;
 
  if (test_) test(transporte,resu,tab_vitesse);
  modifier_flux(*this);
  return resu;
}
 
//Porous correction: add the T*div(u) contribution
//Transported variable: T
//Transporting variable: u
//NOTE: the Kij array MUST NOT be used because by
//construction sum_{j} Kij = 0, which enforces a zero-divergence
//velocity per element — problematic in compressible flows
DoubleTab& Op_Conv_EF_VEF_P1NC_Stab::ajouter_partie_compressible(const DoubleTab& transporte,
                                                                 DoubleTab& resu, const DoubleTab& vitesse_2) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const IntTab& elem_faces=domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const DoubleTab& face_normales=domaine_VEF.face_normales();
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=elem_faces.line_size();
 
  //To account for porosity
  const int marq = phi_u_transportant(equation());
  const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
  const DoubleVect& porosite_face = equation().milieu().porosite_face();
 
  DoubleTab tab_vitesse(vitesse_->valeurs());
  DoubleTabView tab_vitesse_v = tab_vitesse.view_rw();
  CDoubleArrView porosite_face_v = porosite_face.view_ro();
 
  const int vit_size0 = tab_vitesse.dimension(0);
  const int vit_size1 = tab_vitesse.dimension(1);
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       Kokkos::MDRangePolicy<Kokkos::Rank<2>>({0, 0}, {vit_size0, vit_size1}), KOKKOS_LAMBDA(
                         const int i, const int j)
  {
    tab_vitesse_v(i,j)*=porosite_face_v(i);
  });
  end_gpu_timer(__KERNEL_NAME__);
  const int nb_comp=transporte.line_size();
  int nb_dim = dimension;
  int vol_etendus = volumes_etendus_;
 
  // read only
  CDoubleTabView face_normales_v = face_normales.view_ro();
  CIntTabView face_voisins_v = face_voisins.view_ro();
  CIntTabView elem_faces_v = elem_faces.view_ro();
  CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
  CIntArrView elem_nb_faces_dirichlet_v = elem_nb_faces_dirichlet_.view_ro();
  CDoubleArrView porosite_elem_v = porosite_elem.view_ro();
  // Write only
  DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
 
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       Kokkos::RangePolicy<>(0, nb_elem_tot), KOKKOS_LAMBDA(
                         const int elem)
  {
    //Element type: the number of Dirichlet faces
    //it contains
    int type_elem=elem_nb_faces_dirichlet_v(elem);
    double coeff;
    if (!vol_etendus)
      coeff = (nb_dim==2) ? formule_Id_2D(type_elem) : formule_Id_3D(type_elem);
    else
      coeff = (nb_dim==2) ? formule_2D(type_elem) : formule_3D(type_elem);
 
    int facei, facei_loc, dim;
 
    //Compute the divergence per element
    double div=0.;
    for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
      {
        facei=elem_faces_v(elem,facei_loc);
        double signe=(face_voisins_v(facei,0)==elem)? 1.:-1.;
 
        for (dim=0; dim<nb_dim; dim++)
          div+=signe*face_normales_v(facei,dim)*tab_vitesse_v(facei,dim);
      }
    div*=coeff;
    if (!marq) div/=porosite_elem_v(elem);
 
    //Compute the compressible part
    for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
      {
        facei=elem_faces_v(elem,facei_loc);
        for (dim=0; dim<nb_comp; dim++)
          {
            int ligne=facei*nb_comp+dim;
            double delta = div*transporteV[ligne];
            Kokkos::atomic_sub(&resuV[ligne], delta);
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  return resu;
}
 
void Op_Conv_EF_VEF_P1NC_Stab::calculer_flux_bords(const DoubleTab& Kij, const DoubleTab& tab_vitesse, const DoubleTab& transporte) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
  const int nb_comp=transporte.line_size();
 
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1 = 0;
      int num2 = le_bord.nb_faces();//only loop over real faces here
 
      if ( sub_type(Dirichlet_homogene,la_cl.valeur()) )
        {
          //Do not compute flux at Dirichlet_homogene boundaries
        }
      else if ( sub_type(Neumann,la_cl.valeur())
                || sub_type(Neumann_val_ext,la_cl.valeur())
                || sub_type(Neumann_homogene,la_cl.valeur())
                || sub_type(Symetrie,la_cl.valeur())
                || sub_type(Echange_impose_base,la_cl.valeur())
                || sub_type(Dirichlet,la_cl.valeur())
                || sub_type(Periodique,la_cl.valeur()) )
        {
          CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
          CIntArrView num_face = le_bord.num_face().view_ro();
          CDoubleTabView vitesse = tab_vitesse.view_ro();
          CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
          DoubleTabView flux_bords = flux_bords_.view_wo();
          const int nb_dim=Objet_U::dimension;
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                               Kokkos::RangePolicy<>(num1, num2), KOKKOS_LAMBDA(
                                 const int ind_face)
          {
            int facei = num_face(ind_face);
            double psc=0.;
            for (int dim=0; dim<nb_dim; dim++) //tempo fix to compile
              psc-=vitesse(facei,dim)*face_normales(facei,dim);
 
            for (int dim=0; dim<nb_comp; dim++)
              flux_bords(facei,dim)=psc*transporteV[facei*nb_comp+dim];
          });
          end_gpu_timer(__KERNEL_NAME__);
 
        }//end of if on "Neumann", "Neumann_homogene", "Symetrie", "Echange_impose_base"
      else
        {
          Cerr << "Error Op_Conv_EF_VEF_P1NC_Stab::calculer_flux_bords()" << finl;
          Cerr << "Boundary condition " << la_cl.que_suis_je() << " not implemented."   << finl;
          Cerr << "Exiting." << finl;
          Process::exit();
        }//end of else on other boundary conditions
    }
}
 
DoubleTab&
Op_Conv_EF_VEF_P1NC_Stab::ajouter_operateur_centre(const DoubleTab& tab_Kij, const DoubleTab& transporte, DoubleTab& resu) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
  const int nb_comp=transporte.line_size();
 
  CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
  CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
  CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
  DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
                         const int elem, const int facei_loc)
  {
    int facei = elem_faces(elem, facei_loc);
    for (int facej_loc = facei_loc + 1; facej_loc < nb_faces_elem; facej_loc++)
      {
        int facej = elem_faces(elem, facej_loc);
        double kij = Kij(elem, facei_loc, facej_loc);
        double kji = Kij(elem, facej_loc, facei_loc);
 
        for (int dim = 0; dim < nb_comp; dim++)
          {
            int ligne = facei * nb_comp + dim;
            int colonne = facej * nb_comp + dim;
            double delta = transporteV[colonne] - transporteV[ligne];
            Kokkos::atomic_add(&resuV[ligne], kij * delta);
            Kokkos::atomic_sub(&resuV[colonne], kji * delta);
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
  return resu;
}
 
DoubleTab&
Op_Conv_EF_VEF_P1NC_Stab::ajouter_diffusion(const DoubleTab& tab_Kij, const DoubleTab& transporte, DoubleTab& resu) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
  const int nb_comp=transporte.line_size();
 
  CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
  CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
  CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
  CDoubleArrView alpha_tab = alpha_tab_.view_ro();
  DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
                         const int elem, const int facei_loc)
  {
    int facei=elem_faces(elem,facei_loc);
    for (int facej_loc=facei_loc+1; facej_loc<nb_faces_elem; facej_loc++)
      {
        int facej=elem_faces(elem,facej_loc);
        double dij=Dij(elem,facei_loc,facej_loc,Kij);
        double coeffij=alpha_tab[facei]*dij;
        double coeffji=alpha_tab[facej]*dij;
 
        for (int dim=0; dim<nb_comp; dim++)
          {
            int ligne=facei*nb_comp+dim;
            int colonne=facej*nb_comp+dim;
            double delta=transporteV[colonne]-transporteV[ligne];
 
            //NOTE ON SIGN: here we code +div(uT)
            //NOTE: exploiting the symmetry of the operator
            Kokkos::atomic_add(&resuV[ligne], coeffij*delta);
            Kokkos::atomic_sub(&resuV[colonne], coeffji*delta);
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
  return resu;
}
 
DoubleTab&
Op_Conv_EF_VEF_P1NC_Stab::ajouter_antidiffusion(const DoubleTab& tab_Kij, const DoubleTab& transporte, DoubleTab& resu) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
  const int nb_comp=transporte.line_size();
  if (nb_comp>3) Process::exit("EF_stab is not coded for more than 3 components for array transporte.");
 
  //For the limiter
 
  CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
  CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
  CIntTabView num_fac_loc = domaine_VEF.get_num_fac_loc().view_ro();
  CDoubleArrView transporteV = static_cast<const DoubleVect&>(transporte).view_ro();
  CDoubleArrView alpha_tab = alpha_tab_.view_ro();
  CDoubleArrView beta = beta_.view_ro();
  CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
  DoubleArrView resuV = static_cast<DoubleVect&>(resu).view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
                         const int elem, const int facei_loc)
  {
    double P_plus[3],P_moins[3],Q_plus[3],Q_moins[3];
    int facei = elem_faces(elem, facei_loc);
    calculer_senseur(Kij, transporteV, nb_comp, facei, elem_faces, face_voisins, num_fac_loc, P_plus, P_moins,
                     Q_plus, Q_moins);
 
    for (int facej_loc = 0; facej_loc < nb_faces_elem; facej_loc++)
      if (facej_loc != facei_loc)
        {
          int facej = elem_faces(elem, facej_loc);
 
          double kij = Kij(elem, facei_loc, facej_loc);
          double kji = Kij(elem, facej_loc, facei_loc);
          double dij = Dij(elem, facei_loc, facej_loc, Kij);
          double lij = kij + dij;
          double lji = kji + dij;
          assert(lij >= 0);
          assert(lji >= 0);
 
          if (lij <= lji) //facei is upstream
            {
              int face_amont = facei;
              int face_aval = facej;
 
              //If lij==lji, we iterate twice in the loop
              //hence the coefficient 1/2
              double coeff = 1. * (lij < lji) + 0.5 * (lij == lji);
              assert(coeff == 1. || coeff == 0.5);
 
              // Registers for performance:
              double alpha_beta_amont = alpha_tab[face_amont] * beta[face_amont];
              double alpha_beta_aval  = alpha_tab[face_aval]  * beta[face_aval];
 
              for (int dim = 0; dim < nb_comp; dim++)
                {
                  int ligne = face_aval * nb_comp + dim;
                  int colonne = face_amont * nb_comp + dim;
 
                  double delta = transporteV[colonne] - transporteV[ligne];
 
                  //Slope limiter
                  double R;
                  if (delta >= 0.) R = (Kokkos::fabs(P_plus[dim]) < DMINFLOAT) ? 0. : Q_plus[dim] / P_plus[dim];
                  else R = (Kokkos::fabs(P_moins[dim]) < DMINFLOAT) ? 0. : Q_moins[dim] / P_moins[dim];
 
                  double limit = limiteur(R);
                  //double daij = minimum(limit * dij, lji);
                  double daij = Kokkos::fmin(limit * dij, lji);
                  assert(daij >= 0);
                  assert(daij <= lji);
 
                  double coeffij = alpha_beta_amont * daij * coeff * delta;
                  double coeffji = alpha_beta_aval  * daij * coeff * delta;
 
                  //Compute resu
                  Kokkos::atomic_add(&resuV[colonne], + coeffij);
                  Kokkos::atomic_add(&resuV[ligne],   - coeffji);
                }
            }
        }
  });
  end_gpu_timer(__KERNEL_NAME__);
  return resu;
}
 
//WARNING: assumes that the parameters P_plus, P_moins, Q_plus, Q_moins are zero on entry
KOKKOS_INLINE_FUNCTION void
Op_Conv_EF_VEF_P1NC_Stab::calculer_senseur(CDoubleTabView3 Kij, CDoubleArrView transporteV,
                                           const int nb_comp, const int face_i,
                                           CIntTabView elem_faces, CIntTabView face_voisins, CIntTabView num_fac_loc,
                                           double* P_plus, double* P_moins,
                                           double* Q_plus, double* Q_moins) const
{
  for (int i = 0; i < nb_comp; i++)
    {
      P_plus[i] = 0., P_moins[i] = 0.;
      Q_plus[i] = 0., Q_moins[i] = 0.;
    }
  const int nb_faces_elem=(int)elem_faces.extent(1);
  for (int elem_voisin=0; elem_voisin<2; elem_voisin++)
    {
      int elem = face_voisins(face_i,elem_voisin);
      if (elem!=-1)
        {
          int face_i_loc = num_fac_loc(face_i,elem_voisin);
          //Work on the faces of "elem"
          for (int face_k_loc=0; face_k_loc<nb_faces_elem; face_k_loc++)
            {
              int face_k=elem_faces(elem,face_k_loc);
              double kik=Kij(elem,face_i_loc,face_k_loc);
              //
              //Compute the intermediate variables
              //
              for (int dim=0; dim<nb_comp; dim++)
                {
                  double deltaki = transporteV[face_k*nb_comp+dim]-transporteV[face_i*nb_comp+dim];
                  //Compute P_plus and P_moins
                  //Compute Q_plus and Q_moins
                  /* Initial coding:
                     P_plus(dim)+=minimum(0.,kik)*minimum(0.,deltaki);
                     P_moins(dim)+=minimum(0.,kik)*maximum(0.,deltaki);
                     Q_plus(dim)+=maximum(0.,kik)*maximum(0.,deltaki);
                     Q_moins(dim)+=maximum(0.,kik)*minimum(0.,deltaki);
                  */
                  // Optimized coding:
                  double tmp = kik*deltaki;
                  if (kik>0)
                    {
                      if (tmp>0) Q_plus[dim]+=tmp;
                      else       Q_moins[dim]+=tmp;
                    }
                  else
                    {
                      if (tmp>0) P_plus[dim]+=tmp;
                      else       P_moins[dim]+=tmp;
                    }
                }//end of for loop over "dim"
              //
              //End of intermediate variable computation
              //
            }//end of for loop over "face_k_loc"
        }//end of if on "elem!=-1"
    }//end of for loop over "elem_voisin"
}
 
//WARNING: assumes that the parameters P_plus, P_moins, Q_plus, Q_moins are zero on entry
inline void
Op_Conv_EF_VEF_P1NC_Stab::calculer_senseur(const DoubleTab& Kij, const DoubleVect& transporteV,
                                           const int nb_comp, const int face_i,
                                           const IntTab& elem_faces, const IntTab& face_voisins, const IntTab& num_fac_loc,
                                           ArrOfDouble& P_plus, ArrOfDouble& P_moins,
                                           ArrOfDouble& Q_plus, ArrOfDouble& Q_moins) const
{
  assert(P_plus.size_array()==nb_comp);
  assert(Q_plus.size_array()==nb_comp);
  assert(P_moins.size_array()==nb_comp);
  assert(Q_moins.size_array()==nb_comp);
  const int nb_faces_elem=elem_faces.dimension(1);
  for (int elem_voisin=0; elem_voisin<2; elem_voisin++)
    {
      int elem = face_voisins(face_i,elem_voisin);
      if (elem!=-1)
        {
          int face_i_loc = num_fac_loc(face_i,elem_voisin);
          assert(face_i_loc>=0);
          assert(face_i_loc<nb_faces_elem);
          //Work on the faces of "elem"
          for (int face_k_loc=0; face_k_loc<nb_faces_elem; face_k_loc++)
            {
              int face_k=elem_faces(elem,face_k_loc);
              double kik=Kij(elem,face_i_loc,face_k_loc);
              //
              //Compute the intermediate variables
              //
              for (int dim=0; dim<nb_comp; dim++)
                {
                  double deltaki = transporteV[face_k*nb_comp+dim]-transporteV[face_i*nb_comp+dim];
                  //Compute P_plus and P_moins
                  //Compute Q_plus and Q_moins
                  /* Initial coding:
                     P_plus(dim)+=minimum(0.,kik)*minimum(0.,deltaki);
                     P_moins(dim)+=minimum(0.,kik)*maximum(0.,deltaki);
                     Q_plus(dim)+=maximum(0.,kik)*maximum(0.,deltaki);
                     Q_moins(dim)+=maximum(0.,kik)*minimum(0.,deltaki);
                  */
                  // Optimized coding:
                  double tmp = kik*deltaki;
                  if (kik>0)
                    {
                      if (tmp>0) Q_plus[dim]+=tmp;
                      else       Q_moins[dim]+=tmp;
                    }
                  else
                    {
                      if (tmp>0) P_plus[dim]+=tmp;
                      else       P_moins[dim]+=tmp;
                    }
                  assert(P_plus[dim]>=0);
                  assert(Q_plus[dim]>=0);
                  assert(P_moins[dim]<=0);
                  assert(Q_moins[dim]<=0);
                }//end of for loop over "dim"
              //
              //End of intermediate variable computation
              //
            }//end of for loop over "face_k_loc"
        }//end of if on "elem!=-1"
    }//end of for loop over "elem_voisin"
}
 
void Op_Conv_EF_VEF_P1NC_Stab::test(const DoubleTab& transporte, const DoubleTab& resu, const DoubleTab& tab_vitesse) const
{
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const int nb_faces_elem=elem_faces.dimension(1);
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
 
  DoubleTab Kij2(nb_elem_tot,nb_faces_elem, nb_faces_elem);
  Kij2=0.;
 
  DoubleTab Kij_ancien(nb_elem_tot,nb_faces_elem, nb_faces_elem);
  Kij_ancien=0.;
 
  test_difference_Kij(transporte,Kij2,Kij_ancien,tab_vitesse);
  test_difference_resu(Kij2,Kij_ancien,transporte,resu,tab_vitesse);
}
 
void Op_Conv_EF_VEF_P1NC_Stab::test_difference_Kij(const DoubleTab& transporte, DoubleTab& Kij, DoubleTab& Kij_ancien, const DoubleTab& tab_vitesse ) const
{
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const DoubleTab& face_normales=domaine_VEF.face_normales();
  const int nb_faces_elem=elem_faces.dimension(1);
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
 
  const int nb_comp = (transporte.nb_dim()!=1) ? transporte.dimension(1) : 1;
 
  calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
 
  //
  //   Compute Kij_ancien:
  //
  for(int elem0=0; elem0<nb_elem_tot; elem0++)
    {
      int face_locj=0;
      for(int face_loci=0; face_loci<nb_faces_elem; face_loci++)
        {
          int face_i0=elem_faces(elem0,face_loci);
          double signei=1.0;
          if(face_voisins(face_i0,0)!=elem0)
            signei=-1.0;
          double psci=0;
          for(int comp=0; comp<dimension; comp++)
            psci+=tab_vitesse(face_i0,comp)*face_normales(face_i0,comp);
          psci*=signei;
          //Kij_ancien(elem,face_loci,face_loci)=0.;
          for(face_locj=face_loci+1; face_locj<nb_faces_elem; face_locj++)
            {
              int face_j0=elem_faces(elem0,face_locj);
              double signej=1.0;
              if(face_voisins(face_j0,0)!=elem0)
                signej=-1.0;
 
              double pscj=0;
              //psci=0;
              for(int comp=0; comp<dimension; comp++)
                pscj+=tab_vitesse(face_j0,comp)*face_normales(face_j0,comp);
              pscj*=signej;
              Kij_ancien(elem0,face_loci,face_locj)=-1./nb_faces_elem*pscj;
              Kij_ancien(elem0,face_loci,face_loci)+=1./nb_faces_elem*pscj;
              Kij_ancien(elem0,face_locj,face_loci)=-1./nb_faces_elem*psci;
              Kij_ancien(elem0,face_locj,face_locj)+=1./nb_faces_elem*psci;
            }
        }
    }
  // Correction of Kij_ancien for Dirichlet!
  {
    int nb_bord=domaine_Cl_VEF.nb_cond_lim();
    double coeff=1./dimension;
    for (int n_bord=0; n_bord<nb_bord; n_bord++)
      {
        const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
        const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
        int nb_faces_tot=le_bord.nb_faces_tot();
 
        if ( (sub_type(Dirichlet,la_cl.valeur()))
             || (sub_type(Dirichlet_homogene,la_cl.valeur()))
           )
          {
            if ( volumes_etendus_ )
              {
                for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
                  {
                    int face = le_bord.num_face(ind_face);
                    int elem=face_voisins(face,0);
                    assert(elem!=-1);
                    int face_loc_j;
                    int face_j=-1;
                    for (face_loc_j=0; (face_loc_j<nb_faces_elem && face_j!=face); face_loc_j++)
                      {
                        face_j=elem_faces(elem,face_loc_j);
                      }
                    face_loc_j--;
                    assert(face_loc_j>=0);
                    assert(face_loc_j<nb_faces_elem);
                    assert(elem_faces(elem,face_loc_j)==face);
                    const double kjj=Kij_ancien(elem,face_loc_j,face_loc_j);
                    for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                      {
                        int face_i=elem_faces(elem,face_loc_i);
                        if(face_i!=face)
                          {
                            double& kii=Kij_ancien(elem,face_loc_i,face_loc_i);
                            const double kji=Kij_ancien(elem,face_loc_j,face_loc_i);
                            kii+=coeff*kji;
                            double& kij=Kij_ancien(elem,face_loc_i,face_loc_j);
                            kij+=coeff*kjj;
                            for (int face_loc_k=(face_loc_i+1); face_loc_k<nb_faces_elem; face_loc_k++)
                              {
                                int face_k=elem_faces(elem,face_loc_k);
                                if(face_k!=face)
                                  {
                                    double& kik=Kij_ancien(elem,face_loc_i,face_loc_k);
                                    const double kjk=Kij_ancien(elem,face_loc_j,face_loc_k);
                                    double& kki=Kij_ancien(elem,face_loc_k,face_loc_i);
                                    kik+=coeff*kjk;
                                    kki+=coeff*kji;
                                  }
                              }
                          }
                      }
                    {
                      for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                        Kij_ancien(elem,face_loc_j,face_loc_i)=0;
                    }
                    {
                      for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                        {
                          double sum=0.;
                          for (int face_loc_k=0; face_loc_k<nb_faces_elem; face_loc_k++)
                            {
                              sum+=Kij_ancien(elem,face_loc_i,face_loc_k);
                            }
                          //Cerr << "somme apres : " << sum << finl;
                          Kij_ancien(elem,face_loc_i,face_loc_i)-=sum;//car div(u)=0!
                        }
                    }
                  }// for face
 
              }//end of if on "volumes_etendus_"
 
          }// sub_type Dirichlet
      }
  }
 
  Kij-=Kij_ancien;
  const double max_kij = local_max_abs_vect(Kij);
  if (max_kij > 1.e-15)
    {
      Cerr << "Error in Kij computation: " << max_kij << finl;
      Cerr << "Exiting" << finl;
      Process::exit();
    }
 
  Kij+=Kij_ancien;
}
 
void Op_Conv_EF_VEF_P1NC_Stab::test_difference_resu(const DoubleTab& Kij, const DoubleTab& Kij_ancien,
                                                    const DoubleTab& transporte,const DoubleTab& resu, const DoubleTab& tab_vitesse) const
{
  DoubleTab resu1(resu);
  resu1=0;
 
  if (is_compressible_) ajouter_partie_compressible(transporte,resu1,tab_vitesse);
  ajouter_operateur_centre(Kij,transporte,resu1);
  ajouter_diffusion(Kij,transporte,resu1);
  ajouter_antidiffusion(Kij,transporte,resu1);
  mettre_a_jour_pour_periodicite(resu1);
 
  DoubleTab resu2(resu);
  resu2=0;
 
  //
  //   Compute resu2
  //
 
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
 
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const int nb_faces_elem=elem_faces.line_size();
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const int nb_comp=resu.line_size();
  const int nb_faces0 = transporte.dimension(0);
 
  ArrOfDouble pplusi(nb_comp);
  ArrOfDouble qplusi(nb_comp);
  ArrOfDouble pmoinsi(nb_comp);
  ArrOfDouble qmoinsi(nb_comp);
 
  for(int elem=0; elem<nb_elem_tot; elem++)
    {
      int face_locj=0;
 
      for(int face_loci=0; face_loci<nb_faces_elem; face_loci++)
        {
          int face_i0=elem_faces(elem,face_loci);
          for(int comp0=0; comp0<nb_comp; comp0++)
            resu2(face_i0,comp0)+=Kij_ancien(elem,face_loci,face_loci)*transporte(face_i0,comp0);
 
          pplusi=0.;
          qplusi=0.;
          pmoinsi=0.;
          qmoinsi=0.;
 
          int elem1=face_voisins(face_i0,0), elem2=face_voisins(face_i0,1);
          int face_loc_i=0, face_loc_j=0;
          double dT,min_dT,max_dT, K,min_K,max_K;
 
          //
          // in elem1:
          //
          while((face_loc_i<nb_faces_elem)&&(elem_faces(elem1,face_loc_i)!=face_i0))
            face_loc_i++;
          if(face_loc_i==nb_faces_elem)
            {
              //Periodique!!
              assert(elem2!=-1);
              face_loc_i=0;
              while( (face_loc_i<nb_faces_elem) &&
                     (face_voisins(elem_faces(elem1,face_loc_i),0)!=elem2) &&
                     (face_voisins(elem_faces(elem1,face_loc_i),1)!=elem2) )
                face_loc_i++;
            }
 
          assert(face_loc_i<nb_faces_elem);
          for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
            {
              int face_j=elem_faces(elem1,face_loc_j);
              if(face_j!=face_i0)
                {
                  K=Kij_ancien(elem1,face_loc_i,face_loc_j);
 
                  if(K>0.)
                    {
                      max_K=K ;
                      min_K=0.;
                    }
                  else
                    {
                      max_K=0.;
                      min_K=K ;
                    }
 
                  for(int comp0=0; comp0<nb_comp; comp0++)
                    {
                      dT =transporte(face_j,comp0);
                      dT-=transporte(face_i0,comp0);
 
                      if(dT>0.)
                        {
                          max_dT=dT ;
                          min_dT=0 ;
                        }
                      else
                        {
                          max_dT=0. ;
                          min_dT=dT ;
                        }
 
                      pplusi[comp0] +=min_K*min_dT;
                      pmoinsi[comp0]+=min_K*max_dT;
                      qplusi[comp0] +=max_K*max_dT;
                      qmoinsi[comp0]+=max_K*min_dT;
                    }
                }
            }
 
          //
          // in elem2:
          //
 
          if(elem2!=-1)
            {
              face_loc_i=0;
              while((face_loc_i<nb_faces_elem)&&(elem_faces(elem2,face_loc_i)!=face_i0))
                face_loc_i++;
              if(face_loc_i==nb_faces_elem)
                {
                  //Periodique!!
                  face_loc_i=0;
                  while( (face_loc_i<nb_faces_elem) &&
                         (face_voisins(elem_faces(elem2,face_loc_i),0)!=elem1) &&
                         (face_voisins(elem_faces(elem2,face_loc_i),1)!=elem1) )
                    face_loc_i++;
                }
              assert(face_loc_i<nb_faces_elem);
              for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
                {
                  int face_j=elem_faces(elem2,face_loc_j);
                  if(face_j!=face_i0)
                    {
                      K=Kij_ancien(elem2,face_loc_i,face_loc_j);
 
                      if(K>0.)
                        {
                          max_K=K ;
                          min_K=0.;
                        }
                      else
                        {
                          max_K=0.;
                          min_K=K ;
                        }
 
                      for(int comp0=0; comp0<nb_comp; comp0++)
                        {
                          dT =transporte(face_j,comp0);
                          dT-=transporte(face_i0,comp0);
 
                          if(dT>0.)
                            {
                              max_dT=dT ;
                              min_dT=0 ;
                            }
                          else
                            {
                              max_dT=0. ;
                              min_dT=dT ;
                            }
 
                          pplusi[comp0] +=min_K*min_dT;
                          pmoinsi[comp0]+=min_K*max_dT;
                          qplusi[comp0] +=max_K*max_dT;
                          qmoinsi[comp0]+=max_K*min_dT;
                        }
                    }
                }
            }
 
          for(face_locj=0; face_locj<nb_faces_elem; face_locj++)
            if(face_locj!=face_loci)
              {
                int face_j0=elem_faces(elem,face_locj);
                const double kij=Kij_ancien(elem,face_loci,face_locj);
                const double kji=Kij_ancien(elem,face_locj,face_loci);
                double dij=Dij(elem,face_loci,face_locj,Kij_ancien);
                double lji=kji+dij;
                double lij=kij+dij;
                assert(lij>=0);
                assert(lji>=0);
 
                for(int comp0=0; comp0<nb_comp; comp0++)
                  {
                    const double Ti=transporte(face_i0,comp0);
                    const double Tj=transporte(face_j0,comp0);
                    double deltaij=Ti-Tj;
                    double Fij=0;
                    if(lij<=lji)
                      {
                        double coef=1;
                        if (lij==lji)  coef=.5;
                        if(deltaij)
                          {
                            if(Ti >= Tj)
                              {
                                if(pplusi[comp0])
                                  {
                                    double R=qplusi[comp0]/pplusi[comp0];
                                    Fij=minimum(limiteur(R)*dij,lji);
                                  }
                              }
                            else if(pmoinsi[comp0])
                              {
                                double R=qmoinsi[comp0]/pmoinsi[comp0];
                                Fij=minimum(limiteur(R)*dij,lji);
                              }
 
                            assert(Fij*dij>=0);
                            Fij-=dij;
                            Fij*=deltaij;
                          }
                        resu2(face_i0,comp0)+=coef*(kij*Tj+Fij);
                        resu2(face_j0,comp0)+=coef*(kji*Ti-Fij);
                      }
                  }
              }
        }
    }
 
  // For periodicity
 
  int nb_bord=domaine_Cl_VEF.nb_cond_lim();
  int face;
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1 = le_bord.num_premiere_face();
      int nb_faces_b=le_bord.nb_faces();
      int num2 = num1 + nb_faces_b;
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
          int face_associee;
          IntVect fait(nb_faces_b);
          fait = 0;
 
          for (face=num1; face<num2; face++)
            {
              if (fait(face-num1) == 0)
                {
                  fait(face-num1) = 1;
                  face_associee=la_cl_perio.face_associee(face-num1);
                  fait(face_associee) = 1;
                  for (int comp=0; comp<nb_comp; comp++)
                    resu2(face_associee+num1, comp)=(resu2(face,comp)+=resu2(face_associee+num1,comp));
                }// if fait
            }// for face
        }// sub_type Perio
    }
 
  resu1-=resu2;
 
  //For Dirichlet faces, our algorithm computes no
  //value on those faces since they are overwritten
  //by the mass matrix anyway. This is not the case with the old algorithm,
  //hence the modification here to avoid seeing only the errors
  //made on Dirichlet boundaries, which have no significance
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1 = le_bord.num_premiere_face();
      int nb_faces=le_bord.nb_faces();
      int num2 = num1 + nb_faces;
 
      if (sub_type(Dirichlet,la_cl.valeur()) || sub_type(Dirichlet_homogene,la_cl.valeur()))
        {
          for (face=num1; face<num2; face++)
            for (int dim=0; dim<nb_comp; dim++)
              resu1(face,dim)=0;
        }//end of if on Dirichlet
    }//end of for on "n_bord"
 
  const double max_abs_resu1 = local_max_abs_vect(resu1);
  Journal() << "local_max_abs_vect(resu1) = " << max_abs_resu1
            << " " << equation().schema_temps().temps_courant() << finl;
 
  if (max_abs_resu1 > 1.e-15)
    {
      Cerr << "Error in resu computation: " << max_abs_resu1 << finl;
      Cerr << "Displaying boundary faces" << finl;
 
      /**************************************************/
      Cerr << "Displaying boundary faces." << finl;
 
      for (int n_bord=0; n_bord<nb_bord; n_bord++)
        {
          const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int nb_faces=le_bord.nb_faces();
          int num2 = num1 + nb_faces;
 
          if (sub_type(Periodique,la_cl.valeur()))
            {
              Cerr << "Periodic boundary: ";
              for (face=num1; face<num2; face++)
                {
                  Cerr << face << ",";
                }// for face
              Cerr << finl;
            }// sub_type Perio
 
          else if (sub_type(Dirichlet,la_cl.valeur()) || (sub_type(Dirichlet_homogene,la_cl.valeur())) )
            {
              Cerr << "Dirichlet boundary: ";
              for (face=num1; face<num2; face++)
                {
                  Cerr << face << ",";
                }// for face
              Cerr << finl;
            }
 
        }//end of for on "nbord"
      /**************************************************/
 
      /**************************************************/
      Cerr << "Display of problematic faces: " << finl;
      if (nb_comp==1)
        {
          for (int face_i=0; face_i<nb_faces0; face_i++)
            {
              if (resu1(face_i)>1.e-15)
                Cerr << face_i << "(" << face_voisins(face_i,0) << ","
                     << face_voisins(face_i,1) << ") ; ";
            }//end of for on "face_i"
        }
      else
        {
          for (int face_i=0; face_i<nb_faces0; face_i++)
            {
              Cerr << face_i << "(" << face_voisins(face_i,0) << ","
                   << face_voisins(face_i,1) << ") ";
 
              for (int dim=0; dim<nb_comp; dim++)
                {
                  Cerr  << ", resu1("
                        << face_i << "," << dim << ")= "
                        << resu1(face_i,dim);
                }
              Cerr << finl;
            }
        }
      Cerr << finl;
      /**************************************************/
 
      /**************************************************/
      resu1+=resu2;
 
      for (int n_bord=0; n_bord<nb_bord; n_bord++)
        {
          const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int nb_faces=le_bord.nb_faces();
          int num2 = num1 + nb_faces;
 
          if (sub_type(Periodique,la_cl.valeur()))
            {
              Cerr << "Display of resu values at the periodic boundary" << finl;
              for (face=num1; face<num2; face++)
                {
                  Cerr << "resu1(" << face << ") : " << resu1(face) << finl;
                  Cerr << "resu2(" << face << ") : " << resu2(face) << finl;
                }// for face
              Cerr << finl;
            }// sub_type Perio
 
        }//end of for on "nbord"
 
      /**************************************************/
 
      static int count = 0;
      count++;
      if (count==2)
        {
          Cerr << "Exiting" << finl;
          Process::exit();
        }
    }
}
 
void Op_Conv_EF_VEF_P1NC_Stab::mettre_a_jour_pour_periodicite(DoubleTab& resu) const
{
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const int nb_bord = domaine_Cl_VEF.nb_cond_lim();
  const int nb_comp = (resu.nb_dim()==1) ? 1 : resu.dimension(1);
 
  //Boundary faces
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1=0;
      int num2=le_bord.nb_faces();
 
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
          CIntArrView le_bord_num_face = static_cast<const ArrOfInt&>(le_bord.num_face()).view_ro();
          CIntArrView face_associee = static_cast<const ArrOfInt&>(la_cl_perio.face_associee()).view_ro();
          DoubleArrView resuV = static_cast<ArrOfDouble&>(resu).view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(num1, num2), KOKKOS_LAMBDA(const int ind_face)
          {
            int facei = le_bord_num_face(ind_face);
            int ind_face_associee = face_associee(ind_face);
            int faceiAss = le_bord_num_face(ind_face_associee);
 
            if (facei<faceiAss)
              for (int dim=0; dim<nb_comp; dim++)
                {
                  int ligne=facei*nb_comp+dim;
                  int ligneAss=faceiAss*nb_comp+dim;
 
                  Kokkos::atomic_add(&resuV[ligneAss],resuV[ligne]);
                  Kokkos::atomic_store(&resuV[ligne],resuV[ligneAss]);
                }
 
          });//end of for on "face_i"
          end_gpu_timer(__KERNEL_NAME__);
        }//end of if on "Periodique"
 
    }//end of for on "n_bord"
}
 
void Op_Conv_EF_VEF_P1NC_Stab::ajouter_old(const DoubleTab& transporte, DoubleTab& resu, const DoubleTab& tab_vitesse ) const
{
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  //  const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
  // const DoubleTab& vitesse=la_vitesse.valeurs();
 
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const DoubleTab& face_normales=domaine_VEF.face_normales();
  const int nb_faces_elem=elem_faces.dimension(1);
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
 
  assert(nb_faces_elem==(dimension+1));
 
  int nb_comp=resu.line_size();
 
  int face_i0, face_j0, elem0, comp0;
  DoubleTab Kij(nb_elem_tot,nb_faces_elem, nb_faces_elem);
  //
  //   Compute Kij:
  //
  for(elem0=0; elem0<nb_elem_tot; elem0++)
    {
      int face_loci=0;
      int face_locj=0;
      for(; face_loci<nb_faces_elem; face_loci++)
        {
          face_i0=elem_faces(elem0,face_loci);
          double signei=1.0;
          if(face_voisins(face_i0,0)!=elem0)
            signei=-1.0;
          double psci=0;
          for(comp0=0; comp0<dimension; comp0++)
            psci+=tab_vitesse(face_i0,comp0)*face_normales(face_i0,comp0);
          psci*=signei;
          //Kij(elem,face_loci,face_loci)=0.;
          for(face_locj=face_loci+1; face_locj<nb_faces_elem; face_locj++)
            {
              face_j0=elem_faces(elem0,face_locj);
              double signej=1.0;
              if(face_voisins(face_j0,0)!=elem0)
                signej=-1.0;
 
              double pscj=0;
              //psci=0;
              for(comp0=0; comp0<dimension; comp0++)
                pscj+=tab_vitesse(face_j0,comp0)*face_normales(face_j0,comp0);
              pscj*=signej;
              Kij(elem0,face_loci,face_locj)=-1./nb_faces_elem*pscj;
              Kij(elem0,face_loci,face_loci)+=1./nb_faces_elem*pscj;
              Kij(elem0,face_locj,face_loci)=-1./nb_faces_elem*psci;
              Kij(elem0,face_locj,face_locj)+=1./nb_faces_elem*psci;
            }
        }
    }
  //
  // Correction of Kij for Dirichlet!
  //
  {
    int nb_bord=domaine_Cl_VEF.nb_cond_lim();
    int face;
    double coeff=1./dimension;
    for (int n_bord=0; n_bord<nb_bord; n_bord++)
      {
        const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
        const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
        int nb_faces_tot=le_bord.nb_faces_tot();
        if (( (sub_type(Dirichlet,la_cl.valeur())) || (sub_type(Dirichlet_homogene,la_cl.valeur())) )
            && ( volumes_etendus_))
          {
            for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
              {
                face=le_bord.num_face(ind_face);
                int elem=face_voisins(face,0);
                assert(elem!=-1);
                int face_loc_j;
                int face_j=-1;
                for (face_loc_j=0; (face_loc_j<nb_faces_elem && face_j!=face); face_loc_j++)
                  {
                    face_j=elem_faces(elem,face_loc_j);
                  }
                face_loc_j--;
                assert(face_loc_j>=0);
                assert(face_loc_j<nb_faces_elem);
                assert(elem_faces(elem,face_loc_j)==face);
                const double kjj=Kij(elem,face_loc_j,face_loc_j);
                for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                  {
                    int face_i=elem_faces(elem,face_loc_i);
                    if(face_i!=face)
                      {
                        double& kii=Kij(elem,face_loc_i,face_loc_i);
                        const double kji=Kij(elem,face_loc_j,face_loc_i);
                        kii+=coeff*kji;
                        double& kij=Kij(elem,face_loc_i,face_loc_j);
                        kij+=coeff*kjj;
                        for (int face_loc_k=(face_loc_i+1); face_loc_k<nb_faces_elem; face_loc_k++)
                          {
                            int face_k=elem_faces(elem,face_loc_k);
                            if(face_k!=face)
                              {
                                double& kik=Kij(elem,face_loc_i,face_loc_k);
                                const double kjk=Kij(elem,face_loc_j,face_loc_k);
                                double& kki=Kij(elem,face_loc_k,face_loc_i);
                                kik+=coeff*kjk;
                                kki+=coeff*kji;
                              }
                          }
                      }
                  }
                {
                  for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                    Kij(elem,face_loc_j,face_loc_i)=0;
                }
                {
                  for (int face_loc_i=0; face_loc_i<nb_faces_elem; face_loc_i++)
                    {
                      double sum=0.;
                      for (int face_loc_k=0; face_loc_k<nb_faces_elem; face_loc_k++)
                        {
                          sum+=Kij(elem,face_loc_i,face_loc_k);
                        }
                      //Cerr << "somme apres : " << sum << finl;
                      Kij(elem,face_loc_i,face_loc_i)-=sum;//car div(u)=0!
                    }
                }
              }// for face
          }// sub_type Dirichlet
      }
  }
 
  //
  //   Compute resu
  //
 
  ArrOfDouble pplusi(nb_comp);
  ArrOfDouble qplusi(nb_comp);
  ArrOfDouble pmoinsi(nb_comp);
  ArrOfDouble qmoinsi(nb_comp);
 
  for(elem0=0; elem0<nb_elem_tot; elem0++)
    {
      int face_loci=0;
      int face_locj=0;
 
      for(; face_loci<nb_faces_elem; face_loci++)
        {
          face_i0=elem_faces(elem0,face_loci);
 
          for(comp0=0; comp0<nb_comp; comp0++)
            resu(face_i0,comp0)+=Kij(elem0,face_loci,face_loci)*transporte(face_i0,comp0);
 
          pplusi=0., qplusi=0., pmoinsi=0., qmoinsi=0.;
 
          int elem1=face_voisins(face_i0,0);
          int elem2=face_voisins(face_i0,1);
          int face_loc_i=0;
          int face_loc_j=0;
          double dT,min_dT,max_dT;
          double K,min_K,max_K;
 
          //
          // in elem1:
          //
          while((face_loc_i<nb_faces_elem)&&(elem_faces(elem1,face_loc_i)!=face_i0))
            face_loc_i++;
          if(face_loc_i==nb_faces_elem)
            {
              //Periodique!!
              assert(elem2!=-1);
              face_loc_i=0;
              while( (face_loc_i<nb_faces_elem) &&
                     (face_voisins(elem_faces(elem1,face_loc_i),0)!=elem2) &&
                     (face_voisins(elem_faces(elem1,face_loc_i),1)!=elem2) )
                face_loc_i++;
            }
          assert(face_loc_i<nb_faces_elem);
          for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
            {
              int face_j=elem_faces(elem1,face_loc_j);
              if(face_j!=face_i0)
                {
                  K=Kij(elem1,face_loc_i,face_loc_j);
 
                  if(K>0.)
                    {
                      max_K=K ;
                      min_K=0.;
                    }
                  else
                    {
                      max_K=0.;
                      min_K=K ;
                    }
 
                  for(comp0=0; comp0<nb_comp; comp0++)
                    {
                      dT =transporte(face_j,comp0);
                      dT-=transporte(face_i0,comp0);
 
                      if(dT>0.)
                        {
                          max_dT=dT ;
                          min_dT=0 ;
                        }
                      else
                        {
                          max_dT=0. ;
                          min_dT=dT ;
                        }
 
                      pplusi[comp0] +=min_K*min_dT;
                      pmoinsi[comp0]+=min_K*max_dT;
                      qplusi[comp0] +=max_K*max_dT;
                      qmoinsi[comp0]+=max_K*min_dT;
                    }
                }
            }
 
          //
          // in elem2:
          //
 
          if(elem2!=-1)
            {
              face_loc_i=0;
              while((face_loc_i<nb_faces_elem)&&(elem_faces(elem2,face_loc_i)!=face_i0))
                face_loc_i++;
              if(face_loc_i==nb_faces_elem)
                {
                  //Periodique!!
                  face_loc_i=0;
                  while( (face_loc_i<nb_faces_elem) &&
                         (face_voisins(elem_faces(elem2,face_loc_i),0)!=elem1) &&
                         (face_voisins(elem_faces(elem2,face_loc_i),1)!=elem1) )
                    face_loc_i++;
                }
              assert(face_loc_i<nb_faces_elem);
              for(face_loc_j=0; face_loc_j<nb_faces_elem; face_loc_j++)
                {
                  int face_j=elem_faces(elem2,face_loc_j);
                  if(face_j!=face_i0)
                    {
                      K=Kij(elem2,face_loc_i,face_loc_j);
 
                      if(K>0.)
                        {
                          max_K=K ;
                          min_K=0.;
                        }
                      else
                        {
                          max_K=0.;
                          min_K=K ;
                        }
 
                      for(comp0=0; comp0<nb_comp; comp0++)
                        {
                          dT =transporte(face_j,comp0);
                          dT-=transporte(face_i0,comp0);
 
                          if(dT>0.)
                            {
                              max_dT=dT ;
                              min_dT=0 ;
                            }
                          else
                            {
                              max_dT=0. ;
                              min_dT=dT ;
                            }
 
                          pplusi[comp0] +=min_K*min_dT;
                          pmoinsi[comp0]+=min_K*max_dT;
                          qplusi[comp0] +=max_K*max_dT;
                          qmoinsi[comp0]+=max_K*min_dT;
                        }
                    }
                }
            }
 
 
 
          for(face_locj=0; face_locj<nb_faces_elem; face_locj++)
            if(face_locj!=face_loci)
              {
                face_j0=elem_faces(elem0,face_locj);
                const double kij=Kij(elem0,face_loci,face_locj);
                const double kji=Kij(elem0,face_locj,face_loci);
                double dij=Dij(elem0,face_loci,face_locj,Kij);
                double lji=kji+dij;
                double lij=kij+dij;
                assert(lij>=0);
                assert(lji>=0);
 
                for(comp0=0; comp0<nb_comp; comp0++)
                  {
                    const double Ti=transporte(face_i0,comp0);
                    const double Tj=transporte(face_j0,comp0);
                    double deltaij=Ti-Tj;
                    double Fij=0;
                    if(lij<=lji)
                      {
                        double coef=1;
                        if (lij==lji)  coef=.5;
                        if(deltaij)
                          {
                            if(Ti >= Tj)
                              {
                                if(pplusi[comp0])
                                  {
                                    double R=qplusi[comp0]/pplusi[comp0];
                                    Fij=minimum(limiteur(R)*dij,lji);
                                  }
                              }
                            else if(pmoinsi[comp0])
                              {
                                double R=qmoinsi[comp0]/pmoinsi[comp0];
                                Fij=minimum(limiteur(R)*dij,lji);
                              }
                            assert(Fij*dij>=0);
                            Fij-=dij;
                            Fij*=deltaij;
                          }
                        resu(face_i0,comp0)+=coef*(kij*Tj+Fij);
                        resu(face_j0,comp0)+=coef*(kji*Ti-Fij);
                      }
                  }
              }
        }
    }
  //Apply here the treatment used for Neumann_sortie_libre boundary conditions in Op_Conv_VEF_Face
  int nb_bord=domaine_Cl_VEF.nb_cond_lim();
  int face;
  const int ncomp_ch_transporte = transporte.line_size();
 
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int nb_faces_tot=le_bord.nb_faces_tot();
 
      double psc;
      int num_face,i;
 
      if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
        {
          const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
          //const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int num2 = num1 + le_bord.nb_faces();
 
          for (num_face=num1; num_face<num2; num_face++)
            {
              psc =0;
              for (i=0; i<dimension; i++)
                psc += tab_vitesse(num_face,i)*face_normales(num_face,i);
              if (psc>0)
                {
                  for (i=0; i<ncomp_ch_transporte; i++)
                    resu(num_face,i) -= psc*transporte(num_face,i);
                }
              else
                {
                  for (i=0; i<ncomp_ch_transporte; i++)
                    resu(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
                  fluent_(num_face) -= psc;
                }
            }
        }
      // For periodicity
      else if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
          int face_associee,ind_face_associee;
          IntVect fait(nb_faces_tot);
          fait = 0;
 
          for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
            {
              face=le_bord.num_face(ind_face);
              if (fait(ind_face) == 0)
                {
                  fait(ind_face) = 1;
                  ind_face_associee=la_cl_perio.face_associee(ind_face);
                  fait(ind_face_associee) = 1;
                  face_associee=le_bord.num_face(ind_face_associee);
                  for (int comp=0; comp<nb_comp; comp++)
                    resu(face_associee, comp)=(resu(face,comp)+=resu(face_associee,comp));
                }// if fait
            }// for face
        }// sub_type Perio
    }
  /*
    ArrOfDouble bilan(nb_comp);
    BilanQdmVEF::bilan_qdm(resu, domaine_Cl_VEF, bilan);
    if(nb_comp==1)
    Cout << "Scalaire Bilan Convectif : " << bilan[0] << finl;
    else
    for (int comp=0; comp<nb_comp; comp++)
    Cout << "Vecteur Bilan Convectif " << comp << " : " << bilan[comp] << finl;
    bilan=0;
    BilanQdmVEF::bilan_energie(resu, transporte, domaine_Cl_VEF, bilan);
    if(nb_comp==1)
    Cout << "Scalaire Bilan Convectif Energie : " << bilan[0] << finl;
    else
    for (int comp=0; comp<nb_comp; comp++)
    Cout << "Vecteur Bilan Convectif Energie " << comp << " : " << bilan[comp] << finl;
    if(nb_comp==1)
    {
    Cout << "min = " << min(transporte);
    Cout << " max = " << std::max(transporte) << finl;
    }
    //  else
    //    for (int comp=0; comp<nb_comp; comp++)
    //    {
    //      Cout << "min("<<comp<<") = " << transporte.min(comp);
    //      Cout << " std::max("<<comp<<") = " << transporte.max(comp) << finl;
    //    }
    Cout << " Ratio Antidiffusion/Diffusion = " << sigma_fija/sigma_fijd << finl;
  */
}
 
//Function that initializes the attributes "elem_nb_faces_dirichlet_"
//and "elem_faces_dirichlet_"
//NOTE: "elem_nb_faces_dirichlet_" contains the number of Dirichlet faces
//for each element of the mesh
//NOTE: "elem_faces_dirichlet_" contains the global indices of Dirichlet faces
//belonging to any element of the mesh
void Op_Conv_EF_VEF_P1NC_Stab::calculer_data_pour_dirichlet()
{
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
 
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
 
  //Sizing and initialization of attributes
  //NOTE: an element cannot have more than
  //(dimension) Dirichlet faces
  elem_nb_faces_dirichlet_.resize(nb_elem_tot);
  elem_faces_dirichlet_.resize(nb_elem_tot,Objet_U::dimension);
  elem_nb_faces_dirichlet_=0;
  elem_faces_dirichlet_=-1;
  elem_faces_frontiere.dimensionner(nb_bord);
 
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int face, nb_faces_tot=le_bord.nb_faces_tot();
 
      if ( (sub_type(Dirichlet,la_cl.valeur()))
           || (sub_type(Dirichlet_homogene,la_cl.valeur()))
         )
        {
          //
          //Filling the arrays
          //
          for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
            {
              face = le_bord.num_face(ind_face);
              const int elem=face_voisins(face,0);
              assert(elem!=-1);
              elem_faces_frontiere[n_bord].append_array(elem);
 
              elem_nb_faces_dirichlet_(elem)+=1;
              assert(elem_nb_faces_dirichlet_(elem)<=Objet_U::dimension);
 
              if (elem_faces_dirichlet_(elem,0)==-1) elem_faces_dirichlet_(elem,0)=face;
              else if (elem_faces_dirichlet_(elem,1)==-1) elem_faces_dirichlet_(elem,1)=face;
              else if (Objet_U::dimension==3 && elem_faces_dirichlet_(elem,2)==-1) elem_faces_dirichlet_(elem,2)=face;
              else
                {
                  Cerr << "Error in Op_Conv_EF_VEF_P1NC_Stab::calculer_data_pour_dirichlet()" << finl;
                  Cerr << "Element number " << elem << " contains more than "
                       << Objet_U::dimension << " Dirichlet faces" << finl;
                  Cerr << "Exiting." << finl;
                  Process::exit();
                }
            }//end of for on "face"
          //
          //End of array filling
          //
 
        }//end of if on "Dirichlet"
      array_trier_retirer_doublons(elem_faces_frontiere[n_bord]);
    }//end of for on "n_bord"
}
 
void Op_Conv_EF_VEF_P1NC_Stab::completer()
{
  Op_Conv_VEF_Face::completer();
  calculer_data_pour_dirichlet();
 
  //  int nb_comp=1;
  //  if (equation().inconnue().valeurs().nb_dim()>1)
  //  nb_comp=equation().inconnue().valeurs().dimension(1);
 
  //  limiteurs_.resize(le_dom_vef->nb_faces_tot(),nb_comp);
  // limiteurs_=0.;
 
  alpha_tab_.resize_array(le_dom_vef->nb_faces_tot());
  alpha_tab_ = alpha_;
  beta_.resize_array(le_dom_vef->nb_faces_tot());
  beta_=1.;
 
  if (ssz_alpha)
    {
      for (int i=0; i<nb_ssz_alpha; i++)
        {
          OBS_PTR(Sous_domaine_VF) la_ssz;
          const Sous_Domaine& le_sous_domaine=equation().probleme().domaine().ss_domaine(noms_ssz_alpha[i]);
          const Domaine_dis_base& le_domaine_dis=le_dom_vef.valeur();
          bool trouve=false;
          for (int ssz=0; ssz<le_domaine_dis.nombre_de_sous_domaines_dis(); ssz++)
            {
              if (le_domaine_dis.sous_domaine_dis(ssz).sous_domaine().est_egal_a(le_sous_domaine))
                {
                  trouve=true;
                  la_ssz=ref_cast(Sous_domaine_VF,le_domaine_dis.sous_domaine_dis(ssz));
                }
            }
 
          if(!trouve)
            {
              Cerr << "Cannot find the discretized sub-domain associated with " << noms_ssz_alpha[i] << finl;
              Process::exit();
            }
          const Sous_domaine_VF& ssz=la_ssz.valeur();
          int nb_faces = ssz.les_faces().size();
 
          for (int face=0; face<nb_faces; face++)
            {
              int la_face=ssz.les_faces()[face];
              beta_[la_face] = 1.;
              alpha_tab_[la_face] = alpha_ssz(i);
            }
        }
    }
 
 
 
 
  if (sous_domaine)
    {
      sous_domaine=false;
      const Sous_Domaine& le_sous_domaine=equation().probleme().domaine().ss_domaine(nom_sous_domaine);
      const Domaine_dis_base& le_domaine_dis=le_dom_vef.valeur();
      for (int ssz=0; ssz<le_domaine_dis.nombre_de_sous_domaines_dis(); ssz++)
        {
          if (le_domaine_dis.sous_domaine_dis(ssz).sous_domaine().est_egal_a(le_sous_domaine))
            {
              sous_domaine=true;
              le_sous_domaine_dis=ref_cast(Sous_domaine_VF,le_domaine_dis.sous_domaine_dis(ssz));
            }
        }
 
      if(!sous_domaine)
        {
          Cerr << "Cannot find the discretized sub-domain associated with " << nom_sous_domaine << finl;
          Process::exit();
        }
 
      const Sous_domaine_VF& ssz=le_sous_domaine_dis.valeur();
      int nb_faces = ssz.les_faces().size();
 
      for (int face=0; face<nb_faces; face++)
        {
          int la_face=ssz.les_faces()[face];
          beta_[la_face] = 0.;
          alpha_tab_[la_face] = 1.;
        }
    }
 
}
 
void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution(const DoubleTab& transporte_2, Matrice_Morse& matrice) const
{
  if (new_jacobienne_==0)
    {
      Op_Conv_VEF_Face::ajouter_contribution(transporte_2, matrice) ;
      return;
    }
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
 
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
  const int nb_comp=transporte_2.line_size();
 
  DoubleTrav Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
  Kij=0.;
 
  //
  //To account for porosity
  //
  const Champ_P1NC& la_vitesse=ref_cast( Champ_P1NC, vitesse_.valeur());
  const DoubleTab& vitesse_2=la_vitesse.valeurs();
  const DoubleVect& porosite_face = equation().milieu().porosite_face();
  DoubleTrav transporte_;
  DoubleTrav vitesse_face_;
 
  // either transporte=phi*transporte_ and vitesse=vitesse_
  // or transporte=transporte_ and vitesse=phi*vitesse_
  // depending on whether transport uses phi*u or u.
  const int marq=phi_u_transportant(equation());
  const DoubleTab& transporte=modif_par_porosite_si_flag(transporte_2,transporte_,!marq,porosite_face);
  const DoubleTab& tab_vitesse=modif_par_porosite_si_flag(vitesse_2,vitesse_face_,marq,porosite_face);
 
  calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
  if (is_compressible_) ajouter_contribution_partie_compressible(transporte,tab_vitesse,matrice);
  ajouter_contribution_operateur_centre(Kij,transporte,matrice);
  ajouter_contribution_diffusion(Kij,transporte,matrice);
 
  if (test_) test_implicite();
}
 
void Op_Conv_EF_VEF_P1NC_Stab::modifier_pour_Cl (Matrice_Morse& matrice, DoubleTab& secmem) const
{
  Op_Conv_VEF_Face::modifier_pour_Cl(matrice,secmem);
}
 
void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_operateur_centre(const DoubleTab& tab_Kij, const DoubleTab& transporte, Matrice_Morse& matrice_morse) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
 
 
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=domaine_VEF.elem_faces().dimension(1);
  const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
 
  const int nb_comp=transporte.line_size();
 
  CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
  CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
  Matrice_Morse_View matrice;
  matrice.set(matrice_morse);
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
                         const int elem, const int facei_loc)
  {
    int facei = elem_faces(elem, facei_loc);
 
    for (int facej_loc = facei_loc+1; facej_loc < nb_faces_elem; facej_loc++)
      {
        int facej = elem_faces(elem, facej_loc);
 
        double kij = Kij(elem, facei_loc, facej_loc);
        double kji = Kij(elem, facej_loc, facei_loc);
 
        for (int dim = 0; dim < nb_comp; dim++)
          {
            int ligne = facei*nb_comp + dim;
            int colonne = facej*nb_comp + dim;
 
            //ATTENTION AU SIGNE : ici on code +div(uT)
            matrice.atomic_add(ligne, ligne, kij);
            matrice.atomic_add(ligne, colonne, -kij);
            matrice.atomic_add(colonne, colonne, kji);
            matrice.atomic_add(colonne, ligne, -kji);
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  //
  //For periodicity
  //
  const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1 = 0;
      int num2=le_bord.nb_faces_tot();//and not nb_faces otherwise some coefficients are missed
 
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
          int faceiAss=0,ind_faceiAss=0;
          const IntTab& face_voisins = domaine_VEF.face_voisins();
          for (int ind_face=num1; ind_face<num2; ind_face++)
            {
              ind_faceiAss=la_cl_perio.face_associee(ind_face);
 
              int facei=le_bord.num_face(ind_face);
              faceiAss=le_bord.num_face(ind_faceiAss);
 
              //To iterate over periodic faces only once
              if (facei<faceiAss)
                for (int elem_loc=0; elem_loc<2; elem_loc++)
                  {
                    int elem=face_voisins(facei,elem_loc);
                    assert(elem!=-1);
 
                    //Compute the local face index within "elem"
                    int facei_loc=num_fac_loc(facei,elem_loc);
                    int faceToComplete;
                    if (facei_loc!=-1)
                      faceToComplete=faceiAss;
                    else
                      {
                        faceToComplete=facei;
                        facei_loc=num_fac_loc(faceiAss,elem_loc);
                        assert(facei_loc!=-1);
                      }
 
                    //Compute the matrix coefficients due to "elem"
                    for (int facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
                      {
                        int facej=elem_faces(elem,facej_loc);
 
                        if (facej_loc!=facei_loc)
                          {
                            double kij=Kij(elem,facei_loc,facej_loc);
                            //double kji=Kij(elem,facej_loc,facei_loc);
 
                            for (int dim=0; dim<nb_comp; dim++)
                              {
                                int ligne=faceToComplete*nb_comp+dim;
                                int colonne=facej*nb_comp+dim;
 
                                //ATTENTION AU SIGNE : ici on code +div(uT)
                                matrice_morse(ligne,ligne)+=kij;
                                matrice_morse(ligne,colonne)-=kij;
                              }
                          }
                      }
                  }
            }
        }
    }
}
 
void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_diffusion(const DoubleTab& tab_Kij, const DoubleTab& transporte, Matrice_Morse& matrice_morse) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=domaine_VEF.elem_faces().line_size();
  const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
  const int nb_comp=transporte.line_size();
 
  CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
  CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
  CDoubleArrView alpha_tab = static_cast<const ArrOfDouble&>(alpha_tab_).view_ro();
  CDoubleTabView3 Kij = tab_Kij.view_ro<3>();
  Matrice_Morse_View matrice;
  matrice.set(matrice_morse);
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_2D({0,0}, {nb_elem_tot,nb_faces_elem}), KOKKOS_LAMBDA(
                         const int elem, const int facei_loc)
  {
    int facei = elem_faces(elem, facei_loc);
 
    for (int facej_loc = facei_loc+1; facej_loc < nb_faces_elem; facej_loc++)
      {
        int facej = elem_faces(elem, facej_loc);
 
        double dij = Dij(elem, facei_loc, facej_loc, Kij);
 
        double coeffij = alpha_tab(facei)*dij;
        double coeffji = alpha_tab(facej)*dij;
 
        for (int dim = 0; dim < nb_comp; dim++)
          {
            int ligne = facei*nb_comp + dim;
            int colonne = facej*nb_comp + dim;
 
            //NOTE ON SIGN: here we code +div(uT)
            //NOTE: exploiting the symmetry of the operator
            matrice.atomic_add(ligne, ligne, coeffij);
            matrice.atomic_add(ligne, colonne, -coeffij);
            matrice.atomic_add(colonne, colonne, coeffji);
            matrice.atomic_add(colonne, ligne, -coeffji);
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  //
  //For periodicity
  //
  const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1 = 0;
      int num2=le_bord.nb_faces_tot();//and not nb_faces() otherwise some coefficients are missed
 
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
          int faceiAss=0,ind_faceiAss=0;
 
          for (int ind_face=num1; ind_face<num2; ind_face++)
            {
              int facei=le_bord.num_face(ind_face);
              ind_faceiAss=la_cl_perio.face_associee(ind_face);
              faceiAss=le_bord.num_face(ind_faceiAss);
 
              //To iterate over periodic faces only once
              if (facei<faceiAss)
                for (int elem_loc=0; elem_loc<2; elem_loc++)
                  {
                    int elem=face_voisins(facei,elem_loc);
                    assert(elem!=-1);
 
                    //Compute the local face index within "elem"
                    int facei_loc=num_fac_loc(facei,elem_loc);
                    int faceToComplete;
                    if (facei_loc!=-1)
                      faceToComplete=faceiAss;
                    else
                      {
                        faceToComplete=facei;
                        facei_loc=num_fac_loc(faceiAss,elem_loc);
                        assert(facei_loc!=-1);
                      }
 
                    //Compute the matrix coefficients due to "elem"
                    for (int facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
                      {
                        int facej=elem_faces(elem,facej_loc);
 
                        if (facej_loc!=facei_loc)
                          {
                            double dij=Dij(elem,facei_loc,facej_loc,tab_Kij);
                            assert(dij>=0);
 
                            double coeffij=alpha_tab_[faceToComplete]*dij;
                            //double coeffji=alpha_tab_[facej]*dij;
 
                            for (int dim=0; dim<nb_comp; dim++)
                              {
                                int ligne=faceToComplete*nb_comp+dim;
                                int colonne=facej*nb_comp+dim;
 
                                //ATTENTION AU SIGNE : ici on code +div(uT)
                                matrice_morse(ligne,ligne)+=coeffij;
                                matrice_morse(ligne,colonne)-=coeffij;
                              }
                          }
                      }
                  }
            }
        }
    }
}
 
//Porous correction: add the T*div(u) contribution
//Transported variable: T
//Transporting variable: u
//NOTE: the Kij array MUST NOT be used because by
//construction sum_{j} Kij = 0, which enforces a zero-divergence
//velocity per element — problematic in compressible flows
void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_partie_compressible(const DoubleTab& transporte, const DoubleTab& vitesse_2,
                                                                        Matrice_Morse& matrice) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const IntTab& elem_faces=domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const DoubleTab& face_normales=domaine_VEF.face_normales();
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=elem_faces.line_size();
  const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
 
  //To account for porosity
  const int marq = phi_u_transportant(equation());
  const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
  const DoubleVect& porosite_face = equation().milieu().porosite_face();
 
  DoubleTrav tab_vitesse(vitesse_->valeurs());
  for (int i=0; i<tab_vitesse.dimension(0); i++)
    for (int j=0; j<tab_vitesse.line_size(); j++)
      tab_vitesse(i,j)*=porosite_face(i);
 
  const int nb_comp=transporte.line_size();
 
  double (*formule)(int);
 
  if (!volumes_etendus_)
    formule= (dimension==2) ? &formule_Id_2D : &formule_Id_3D;
  else
    formule= (dimension==2) ? &formule_2D : &formule_3D;
 
  ToDo_Kokkos("critical");
  for (int elem=0; elem<nb_elem_tot; elem++)
    {
      //Element type: the number of Dirichlet faces
      //it contains
      int type_elem=elem_nb_faces_dirichlet_(elem);
      double coeff=formule(type_elem);
 
      //Compute the divergence per element
      double div=0.;
      for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
        {
          int facei=elem_faces(elem,facei_loc);
          int signe=(face_voisins(facei,0)==elem)? 1.:-1.;
 
          for (int dim=0; dim<dimension; dim++)
            div+=signe*face_normales(facei,dim)*tab_vitesse(facei,dim);
        }
      div*=coeff;
      if (!marq) div/=porosite_elem(elem);
 
      //Compute the compressible part
      for (int facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
        {
          int facei=elem_faces(elem,facei_loc);
 
          for (int dim=0; dim<nb_comp; dim++)
            {
              int ligne=facei*nb_comp+dim;
              matrice(ligne,ligne)+=div;
            }
        }
    }
 
  //
  //For periodicity
  //
  const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1 = 0;
      int num2 = le_bord.nb_faces();//only iterate over real faces
 
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
 
          for (int ind_face=num1; ind_face<num2; ind_face++)
            {
              int facei = le_bord.num_face(ind_face);
              int ind_faceiAss = la_cl_perio.face_associee(ind_face);
              int faceiAss = le_bord.num_face(ind_faceiAss);
 
              //To iterate over periodic faces only once
              if (facei<faceiAss)
                for (int elem_loc=0; elem_loc<2; elem_loc++)
                  {
                    int elem = face_voisins(facei,elem_loc);
                    assert(elem!=-1);
 
                    //Compute the local face index within "elem"
                    int facei_loc=num_fac_loc(facei,elem_loc);
                    int faceToComplete;
                    if (facei_loc!=-1)
                      faceToComplete=faceiAss;
                    else
                      {
                        faceToComplete=facei;
                        facei_loc=num_fac_loc(faceiAss,elem_loc);
                        assert(facei_loc!=-1);
                      }
 
                    //Element type: the number of Dirichlet faces
                    //it contains
                    int type_elem=elem_nb_faces_dirichlet_(elem);
                    double coeff=formule(type_elem);
 
                    //Compute the divergence per element
                    double div=0.;
                    for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
                      {
                        facei=elem_faces(elem,facei_loc);
                        int signe=(face_voisins(facei,0)==elem)? 1.:-1.;
 
                        for (int dim=0; dim<dimension; dim++)
                          div+=signe*face_normales(facei,dim)*tab_vitesse(facei,dim);
                      }
                    div*=coeff;
                    if (!marq) div/=porosite_elem(elem);
 
                    //Compute the compressible part
                    for (int dim=0; dim<nb_comp; dim++)
                      {
                        int ligne=faceToComplete*nb_comp+dim;
                        matrice(ligne,ligne)+=div;
                      }
                  }
            }
        }
    }
}
 
void Op_Conv_EF_VEF_P1NC_Stab::ajouter_contribution_antidiffusion(const DoubleTab& Kij, const DoubleTab& transporte,
                                                                  Matrice_Morse& matrice) const
{
  const Domaine_VEF& domaine_VEF=le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const IntTab& elem_faces=domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const int nb_elem_tot=domaine_VEF.nb_elem_tot();
  const int nb_faces_elem=elem_faces.line_size();
  const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
  const int nb_comp=transporte.line_size();
 
  int elem=0, elem_loc=0, facei=0,facei_loc=0, faceiAss=0, ind_face=0,ind_faceiAss=0;
  int facej=0,facej_loc=0, ligne=0,colonne=0, dim=0, face_amont=0,face_aval=0;
  int faceToComplete=0, num1=0,num2=0, n_bord=0;
  double kij=0.,kji=0.,dij=0., lij=0.,lji=0., daij=0.;
  double delta=0., coeffij=0.,coeffji=0., coeff=0., R=0.;
 
  //For the limiter
  ArrOfDouble P_plus(nb_comp),P_moins(nb_comp);
  ArrOfDouble Q_plus(nb_comp),Q_moins(nb_comp);
  P_plus=0., P_moins=0., Q_plus=0., Q_moins=0.;
 
  const DoubleVect& transporteV = transporte;
  const ArrOfDouble& alpha_tab = alpha_tab_;
  const IntTab& num_fac_loc = domaine_VEF.get_num_fac_loc();
  ToDo_Kokkos("critical");
  for (elem=0; elem<nb_elem_tot; elem++)
    for (facei_loc=0; facei_loc<nb_faces_elem; facei_loc++)
      {
        facei=elem_faces(elem,facei_loc);
        P_plus=0., P_moins=0., Q_plus=0., Q_moins=0.;
        calculer_senseur(Kij,transporteV,nb_comp,facei,elem_faces,face_voisins,num_fac_loc,P_plus,P_moins,Q_plus,Q_moins);
        for (facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
          if (facej_loc!=facei_loc)
            {
              facej=elem_faces(elem,facej_loc);
 
              kij = Kij(elem,facei_loc,facej_loc);
              kji = Kij(elem,facej_loc,facei_loc);
              dij = Dij(elem,facei_loc,facej_loc,Kij);
              lij = kij+dij;
              lji = kji+dij;
              assert(lij>=0);
              assert(lji>=0);
 
              if (lij<=lji) //facei is upstream
                {
                  face_amont = facei;
                  face_aval = facej;
 
                  //If lij==lji, we iterate twice in the loop
                  //hence the coefficient 1/2
                  coeff = 1.*(lij<lji)+0.5*(lij==lji);
                  assert(coeff==1. || coeff==0.5);
 
                  for (dim=0; dim<nb_comp; dim++)
                    {
                      ligne=face_amont*nb_comp+dim;
                      colonne=face_aval*nb_comp+dim;
 
                      delta=transporteV[ligne]-transporteV[colonne];
 
                      //Slope limiter
                      // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/P_plus(dim);
                      // else  R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/P_moins(dim);
 
                      // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/(P_plus(dim)+DMINFLOAT);
                      // else  R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/(P_moins(dim)+DMINFLOAT);
 
                      if (delta>=0.) R=(std::fabs(P_plus[dim])<DMINFLOAT) ? 0. : Q_plus[dim]/P_plus[dim];
                      else  R=(std::fabs(P_moins[dim])<DMINFLOAT) ? 0. : Q_moins[dim]/P_moins[dim];
 
 
                      daij=minimum(limiteur(R)*dij,lji);
                      assert(daij>=0);
                      assert(daij<=lji);
                      coeffij=alpha_tab_[face_amont]*beta_[face_amont]*daij;
                      coeffji=alpha_tab_[face_aval]*beta_[face_aval]*daij;
 
                      //Compute the matrix
                      matrice(ligne,ligne)-=coeffij*coeff;
                      matrice(ligne,colonne)+=coeffij*coeff;
                      matrice(colonne,colonne)-=coeffji*coeff;
                      matrice(colonne,ligne)+=coeffji*coeff;
                    }
                }
            }
      }
 
  //
  //For periodicity
  //
  for (n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      num1 = 0;
      num2=le_bord.nb_faces();//only iterate over real faces
 
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
 
          //For the limiter
          ArrOfDouble Pj_plus(nb_comp),Pj_moins(nb_comp);
          ArrOfDouble Qj_plus(nb_comp),Qj_moins(nb_comp);
          Pj_plus=0., Pj_moins=0.;
          Qj_plus=0., Qj_moins=0.;
 
          for (ind_face=num1; ind_face<num2; ind_face++)
            {
              facei=le_bord.num_face(ind_face);
              ind_faceiAss=la_cl_perio.face_associee(ind_face);
              faceiAss=le_bord.num_face(ind_faceiAss);
 
              //To iterate over periodic faces only once
              if (facei<faceiAss)
                for (elem_loc=0; elem_loc<2; elem_loc++)
                  {
                    elem=face_voisins(facei,elem_loc);
                    assert(elem!=-1);
 
                    //Compute the local face index within "elem"
                    facei_loc=num_fac_loc(facei,elem_loc);
                    if (facei_loc!=-1)
                      faceToComplete=faceiAss;
                    else
                      {
                        faceToComplete=facei;
                        facei_loc=num_fac_loc(faceiAss,elem_loc);
                        assert(facei_loc!=-1);
                      }
 
                    //Compute the coefficient to add to the matrix
                    P_plus=0., P_moins=0.;
                    Q_plus=0., Q_moins=0.;
                    calculer_senseur(Kij,transporteV,nb_comp,faceToComplete,elem_faces,face_voisins,num_fac_loc,P_plus,P_moins,Q_plus,Q_moins);
 
                    for (facej_loc=0; facej_loc<nb_faces_elem; facej_loc++)
                      if (facej_loc!=facei_loc)
                        {
                          facej=elem_faces(elem,facej_loc);
 
                          kij = Kij(elem,facei_loc,facej_loc);
                          kji = Kij(elem,facej_loc,facei_loc);
                          dij = Dij(elem,facei_loc,facej_loc,Kij);
                          lij = kij+dij;
                          lji = kji+dij;
                          assert(lij>=0);
                          assert(lji>=0);
 
                          if (lij<=lji) //faceToComplete is upstream
                            {
                              face_amont=faceToComplete;
                              face_aval=facej;
 
                              //If lij==lji, we iterate twice in the loop
                              //hence the coefficient 1/2
                              coeff = 1.*(lij<lji)+0.5*(lij==lji);
                              assert(coeff==1. || coeff==0.5);
 
                              for (dim=0; dim<nb_comp; dim++)
                                {
                                  ligne=face_amont*nb_comp+dim;
                                  colonne=face_aval*nb_comp+dim;
                                  delta=transporteV[ligne]-transporteV[colonne];
 
                                  //Slope limiter
                                  // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/P_plus(dim);
                                  // else  R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/P_moins(dim);
 
                                  // if (delta>=0.) R=(P_plus(dim)==0.) ? 0. : Q_plus(dim)/(P_plus(dim)+DMINFLOAT);
                                  // else  R=(P_moins(dim)==0.) ? 0. : Q_moins(dim)/(P_moins(dim)+DMINFLOAT);
 
                                  if (delta>=0.) R=(std::fabs(P_plus[dim])<DMINFLOAT) ? 0. : Q_plus[dim]/P_plus[dim];
                                  else  R=(std::fabs(P_moins[dim])<DMINFLOAT) ? 0. : Q_moins[dim]/P_moins[dim];
 
                                  daij=minimum(limiteur(R)*dij,lji);
                                  assert(daij>=0);
                                  assert(daij<=lji);
                                  coeffij=alpha_tab[face_amont]*beta_[face_amont]*daij;
 
                                  //Compute the matrix
                                  matrice(ligne,ligne)-=coeffij*coeff;
                                  matrice(ligne,colonne)+=coeffij*coeff;
                                }
                            }
                          else  //faceToComplete is downstream
                            {
                              face_aval=faceToComplete;
                              face_amont=facej;
                              coeff=1.;
                              Pj_plus=0., Pj_moins=0., Qj_plus=0., Qj_moins=0.;
                              calculer_senseur(Kij,transporteV,nb_comp,facej,elem_faces,face_voisins,num_fac_loc,Pj_plus,Pj_moins,Qj_plus,Qj_moins);
 
                              for (dim=0; dim<nb_comp; dim++)
                                {
                                  ligne=face_amont*nb_comp+dim;
                                  colonne=face_aval*nb_comp+dim;
 
                                  delta=transporteV[ligne]-transporteV[colonne];
 
                                  //Slope limiter
                                  // if (delta>=0.) R=(Pj_plus(dim)==0.) ? 0. : Qj_plus(dim)/Pj_plus(dim);
                                  // else  R=(Pj_moins(dim)==0.) ? 0. : Qj_moins(dim)/Pj_moins(dim);
 
                                  // if (delta>=0.) R=(Pj_plus(dim)==0.) ? 0. : Qj_plus(dim)/(Pj_plus(dim)+DMINFLOAT);
                                  // else  R=(Pj_moins(dim)==0.) ? 0. : Qj_moins(dim)/(Pj_moins(dim)+DMINFLOAT);
 
                                  if (delta>=0.) R=(std::fabs(Pj_plus[dim])<DMINFLOAT) ? 0. : Qj_plus[dim]/Pj_plus[dim];
                                  else  R=(std::fabs(Pj_moins[dim])<DMINFLOAT) ? 0. : Qj_moins[dim]/Pj_moins[dim];
 
                                  daij=minimum(limiteur(R)*dij,lij);
                                  assert(daij>=0);
                                  assert(daij<=lij);
                                  coeffij=alpha_tab[face_aval]*beta_[face_aval]*daij;
 
                                  //Compute the matrix
                                  matrice(colonne,colonne)-=coeffij*coeff;
                                  matrice(colonne,ligne)+=coeffij*coeff;
                                }
                            }
                        }
                  }
            }
        }
    }
}
 
void Op_Conv_EF_VEF_P1NC_Stab::test_implicite() const
{
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
 
  const DoubleTab& unknown=equation().inconnue().valeurs();
  const DoubleTab& tab_vitesse=vitesse_->valeurs();
 
  DoubleTab tab_test(unknown);
  DoubleVect& test2 = tab_test;
  test2 = 0.;
 
  DoubleTab resuExp(unknown);
  DoubleVect& resu2Exp = resuExp;
  resu2Exp = 0.;
 
  DoubleTab resuImp(unknown);
  DoubleVect& resu2Imp = resuImp;
  resu2Imp = 0.;
 
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const int nb_faces_elem = domaine_VEF.domaine().nb_faces_elem();
  const int nb_faces_tot=domaine_VEF.nb_faces_tot();
  const int nb_bord=domaine_Cl_VEF.nb_cond_lim();
 
  DoubleTab Kij(nb_elem_tot,nb_faces_elem,nb_faces_elem);
  Kij=0.;
 
  const int nb_comp=unknown.line_size();
  int size=unknown.dimension(0);
  int face=0,face2=0, faceAss=0, ind_face=0,ind_faceAss=0, n_bord=0, num1=0,num2=0;
 
  SFichier testResu("test.txt");
  SFichier testMat("matrice.txt");
 
  //
  //For periodicity
  //
  IntTab faces_associees(nb_faces_tot);
  for (face=0; face<nb_faces_tot; face++)
    faces_associees(face)=face;
 
  for (n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      num1 = 0;
      num2=le_bord.nb_faces_tot();//to avoid missing any
 
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique,la_cl.valeur());
 
          for (ind_face=num1; ind_face<num2; ind_face++)
            {
              face=le_bord.num_face(ind_face);
              ind_faceAss=la_cl_perio.face_associee(ind_face);
              faceAss=le_bord.num_face(ind_faceAss);
 
              //To iterate over only half the periodic faces
              if (face<faceAss)
                {
                  faces_associees(face)=faceAss;
                  faces_associees(faceAss)=face;
                }
            }
        }
    }
  //
  //End of periodic face treatment
  //
 
  calculer_coefficients_operateur_centre(Kij,nb_comp,tab_vitesse);
 
  //
  //Build the matrix
  //
  Matrice_Morse matrice;
  dimensionner(matrice);
  if (is_compressible_)
    ajouter_contribution_partie_compressible(unknown,tab_vitesse,matrice);
  ajouter_contribution_operateur_centre(Kij,unknown,matrice);
  ajouter_contribution_diffusion(Kij,unknown,matrice);
  ajouter_contribution_antidiffusion(Kij,unknown,matrice);
  matrice.imprimer_formatte(testMat);
  //
  //End of matrix construction
  //
 
  //
  //Compute the explicit operator and comparison
  //
  for (face=0; face<size; face++)
    {
      test2[face]=1.;
      test2[faces_associees(face)]=1.;
 
      /* Compute the explicit operator */
      resuExp=0.;
      if (is_compressible_)
        ajouter_partie_compressible(tab_test,resuExp,tab_vitesse);
      ajouter_operateur_centre(Kij,tab_test,resuExp);
      ajouter_diffusion(Kij,tab_test,resuExp);
      ajouter_antidiffusion(Kij,tab_test,resuExp);
      mettre_a_jour_pour_periodicite(resuExp);
 
      /* Compute the implicit operator */
      resuImp=0.;
      matrice.ajouter_multvect_(tab_test,resuImp);
 
      /* Compute the difference */
      resuExp+=resuImp;
 
      /* Display of the difference */
      testResu<<"*************************"<<finl;
      testResu<<"Face test : "<<face<<finl;
      for (face2=0; face2<size; face2++)
        if (resu2Exp[face2]<=1.e-13)
          testResu<<face2<<" OK"<<finl;
        else
          testResu<<face2<<" residu : "<<resu2Exp[face2]<<finl;
      testResu<<"*************************"<<finl;
 
      test2[face]=0.;
      test2[faces_associees(face)]=0.;
    }
  //
  //End of explicit operator computation and comparison
  //
 
  Process::exit();
}