Modele_turbulence_hyd_LES_SMAGO_DYN_VDF.cpp

/****************************************************************************
* Copyright (c) 2024, 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 <Modele_turbulence_hyd_LES_SMAGO_DYN_VDF.h>
#include <Champ_Fonc_P0_VDF.h>
#include <Schema_Temps_base.h>
#include <Champ_Face_VDF.h>
#include <Domaine_Cl_VDF.h>
#include <Perf_counters.h>
#include <Equation_base.h>
#include <Domaine_VDF.h>
#include <TRUSTTrav.h>
#include <TRUSTTrav.h>
#include <SFichier.h>
#include <Domaine.h>
#include <Debog.h>
#include <Param.h>
 
Implemente_instanciable(Modele_turbulence_hyd_LES_SMAGO_DYN_VDF, "Modele_turbulence_hyd_sous_maille_SMAGO_DYN_VDF", Modele_turbulence_hyd_LES_Smago_VDF);
 
Sortie& Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::printOn(Sortie& s) const { return s << que_suis_je() << " " << le_nom(); }
 
Entree& Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::readOn(Entree& s) { return Modele_turbulence_hyd_LES_Smago_VDF::readOn(s); }
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::set_param(Param& param) const
{
  Modele_turbulence_hyd_LES_Smago_VDF::set_param(param);
  param.ajouter_non_std("stabilise", (this), Param::REQUIRED);
}
 
int Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  Motcle motlutmp;
  if (mot == "stabilise")
    {
      Motcles les_motcles(4);
      {
        les_motcles[0] = "6_points";
        les_motcles[1] = "plans_paralleles";
        les_motcles[2] = "moy_euler";
        les_motcles[3] = "moy_lagrange";
      }
 
      is >> motlutmp;
      int rang = les_motcles.search(motlutmp);
      switch(rang)
        {
        case 0:
          {
            methode_stabilise_ = "6_POINTS";
            break;
          }
        case 1:
          {
            methode_stabilise_ = "PLANS_PARALLELES";
            is >> motlutmp;
            if (motlutmp == "Nb_points")
              {
                is >> N_c_;
                calcul_tableaux_correspondance(N_c_, compt_c_, corresp_c_);
              }
            else
              {
                Cerr << "The keywords stablise plans_paralleles must be followed by the keyword Nb_points.\n" << finl;
                assert(0);
                Process::exit();
              }
            break;
          }
        case 2:
          {
            methode_stabilise_ = "MOY_EULER";
            break;
          }
        case 3:
          {
            methode_stabilise_ = "MOY_LAGRANGE";
            calc_elem_elem();
            break;
          }
        default:
          {
            Cerr << "Error while reading : " << que_suis_je() << finl;
            Cerr << motlutmp << " is not understood. Available keywords are : " << finl;
            Cerr << les_motcles;
            assert(0);
            Process::exit();
          }
        }
      return 1;
    }
  else
    return Modele_turbulence_hyd_LES_Smago_VDF::lire_motcle_non_standard(mot, is);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::associer(const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis)
{
  Modele_turbulence_hyd_LES_Smago_VDF::associer(domaine_dis, domaine_Cl_dis);
 
  Cerr << "Discretisation de coeff_field" << finl;
  Cerr << "coeff_field discretization" << finl;
  coeff_field_.typer("Champ_Fonc_P0_VDF");
  Champ_Fonc_P0_VDF& coeff = ref_cast(Champ_Fonc_P0_VDF, coeff_field_.valeur());
  coeff.associer_domaine_dis_base(le_dom_VF_.valeur());
  coeff.nommer("dynamic_coefficient");
  coeff.fixer_nb_comp(1);
  coeff.fixer_nb_valeurs_nodales(le_dom_VF_->nb_elem());
  coeff.fixer_unite("adim");
  coeff.changer_temps(0.);
 
  champs_compris_.ajoute_champ(coeff_field_);
}
 
int Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::preparer_calcul()
{
  Modele_turbulence_hyd_LES_base::preparer_calcul();
  mettre_a_jour(0.);
  return 1;
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::mettre_a_jour(double temps)
{
  statistics().begin_count(STD_COUNTERS::turbulent_viscosity, statistics().get_last_opened_counter_level()+1);
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  const Domaine_Cl_VDF& domaine_Cl_VDF = ref_cast(Domaine_Cl_VDF, le_dom_Cl_.valeur());
 
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  DoubleTab Sij_test_scale(0, dimension, dimension);
  domaine_VDF.domaine().creer_tableau_elements(Sij_test_scale);
  DoubleTrav Sij_grid_scale(Sij_test_scale);
 
  DoubleTrav S_test_scale_norme;
  domaine_VDF.domaine().creer_tableau_elements(S_test_scale_norme);
  DoubleTrav S_grid_scale_norme;
  domaine_VDF.domaine().creer_tableau_elements(S_grid_scale_norme);
 
  if (cell_cent_vel_.size()==0)
    {
      cell_cent_vel_.resize(0, dimension);
      domaine_VDF.domaine().creer_tableau_elements(cell_cent_vel_);
    }
 
 
  DoubleTab filt_vel(0, dimension);
  domaine_VDF.domaine().creer_tableau_elements(filt_vel);
  DoubleTrav Lij(Sij_test_scale);
  DoubleTrav Mij(Sij_test_scale);
  DoubleTrav l(nb_elem_tot);
 
  calculer_length_scale(l, domaine_VDF);
  calculer_cell_cent_vel(cell_cent_vel_, domaine_VDF, mon_equation_->inconnue());
  calculer_filter_field(cell_cent_vel_, filt_vel, domaine_VDF);
 
  calculer_Lij(cell_cent_vel_, filt_vel, Lij);
 
  calculer_Sij(Sij_grid_scale, domaine_VDF, domaine_Cl_VDF, mon_equation_->inconnue());
  calculer_Sij_vel_filt(filt_vel, Sij_test_scale, domaine_VDF);
 
  calculer_S_norme(Sij_grid_scale, S_grid_scale_norme, domaine_VDF.domaine().nb_elem_tot());
  calculer_S_norme(Sij_test_scale, S_test_scale_norme, domaine_VDF.domaine().nb_elem_tot());
 
  calculer_Mij(Sij_grid_scale, Sij_test_scale, l, Mij);
 
  calculer_model_coefficient(Lij, Mij);
 
  calculer_viscosite_turbulente(S_grid_scale_norme, l);
  calculer_energie_cinetique_turb(S_grid_scale_norme, l);
  loipar_->calculer_hyd(la_viscosite_turbulente_->valeurs(), energie_cinetique_turbulente().valeurs());
  limiter_viscosite_turbulente();
  controler_grandeurs_turbulentes();
  coeff_field_->changer_temps(temps);
  la_viscosite_turbulente_->valeurs().echange_espace_virtuel();
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::mettre_a_jour la_viscosite_turbulente_ after", la_viscosite_turbulente_->valeurs());
 
  statistics().end_count(STD_COUNTERS::turbulent_viscosity);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_cell_cent_vel(DoubleTab& cell_cent_vel, const Domaine_VDF& domaine_VDF, Champ_Inc_base& inco)
{
  const DoubleTab& vitesse = inco.valeurs();
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_cell_cent_vel INPUT vitesse (inco.valeurs)", vitesse);
 
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  int num0, num1, num2, num3, num4 = -1, num5 = -1;
 
  // This is to calculate the cell centered velocity values
  for (int element_number = 0; element_number < nb_elem_tot; element_number++)
    {
      num0 = elem_faces(element_number, 0);
      num1 = elem_faces(element_number, 1);
      num2 = elem_faces(element_number, 2);
      num3 = elem_faces(element_number, 3);
      if (dimension == 3)
        {
          num4 = elem_faces(element_number, 4);
          num5 = elem_faces(element_number, 5);
        }
 
      if (dimension == 2)
        {
          cell_cent_vel(element_number, 0) = 0.5 * (vitesse[num0] + vitesse[num2]);
          cell_cent_vel(element_number, 1) = 0.5 * (vitesse[num1] + vitesse[num3]);
        }
      else
        {
          cell_cent_vel(element_number, 0) = 0.5 * (vitesse[num0] + vitesse[num3]);
          cell_cent_vel(element_number, 1) = 0.5 * (vitesse[num1] + vitesse[num4]);
          cell_cent_vel(element_number, 2) = 0.5 * (vitesse[num2] + vitesse[num5]);
        }
    }
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_cell_cent_vel OUTPUT cell_cent_vel", cell_cent_vel);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_field(const DoubleTab& in_vel, DoubleTab& out_vel, const Domaine_VDF& domaine_VDF)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_field INPUT in_vel", in_vel);
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  int nb_elem = domaine_VDF.domaine().nb_elem();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  int element_number;
  int num0, num1, num2, num3, num4, num5;
  int f0, f1, f2, f3, f4, f5;
  int i;
 
  DoubleTrav temp1(out_vel);
  DoubleTrav temp2(out_vel);
 
  // This is to calculate the filtered velocity field
 
  if (dimension == 2)
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 1);
          if ((num0 == -1) || (num2 == -1))
            for (i = 0; i < dimension; i++)
              temp1(element_number, i) = in_vel(element_number, i);
          else
            {
              for (i = 0; i < dimension; i++)
                {
                  temp1(element_number, i) = 0.25 * (in_vel(num0, i) + 2.0 * in_vel(element_number, i) + in_vel(num2, i));
                }
            }
        }
      temp1.echange_espace_virtuel();
      Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_field temp1 1", temp1);
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if ((num1 == -1) || (num3 == -1))
            for (i = 0; i < dimension; i++)
              out_vel(element_number, i) = temp1(element_number, i);
          else
            {
              for (i = 0; i < dimension; i++)
                {
                  out_vel(element_number, i) = 0.25 * (temp1(num1, i) + 2.0 * temp1(element_number, i) + temp1(num3, i));
                }
            }
        }
      out_vel.echange_espace_virtuel();
      Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_field out_vel 1", out_vel);
    }
  else
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if ((num0 == -1) || (num3 == -1))
            for (i = 0; i < dimension; i++)
              temp1(element_number, i) = in_vel(element_number, i);
          else
            {
              for (i = 0; i < dimension; i++)
                {
                  temp1(element_number, i) = 0.25 * (in_vel(num0, i) + 2.0 * in_vel(element_number, i) + in_vel(num3, i));
                }
            }
        }
      temp1.echange_espace_virtuel();
      Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_field temp1 2", temp1);
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          f4 = elem_faces(element_number, 4);
          num4 = face_voisins(f4, 1);
          if ((num1 == -1) || (num4 == -1))
            for (i = 0; i < dimension; i++)
              temp2(element_number, i) = temp1(element_number, i);
          else
            {
              for (i = 0; i < dimension; i++)
                {
                  temp2(element_number, i) = 0.25 * (temp1(num1, i) + 2.0 * temp1(element_number, i) + temp1(num4, i));
                }
            }
        }
      temp2.echange_espace_virtuel();
      Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_field temp2 1", temp2);
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 0);
          f5 = elem_faces(element_number, 5);
          num5 = face_voisins(f5, 1);
          if ((num2 == -1) || (num5 == -1))
            for (i = 0; i < dimension; i++)
              out_vel(element_number, i) = temp2(element_number, i);
          else
            {
              for (i = 0; i < dimension; i++)
                {
                  out_vel(element_number, i) = 0.25 * (temp2(num2, i) + 2.0 * temp2(element_number, i) + temp2(num5, i));
                }
            }
        }
      out_vel.echange_espace_virtuel();
    }
 
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_field OUTPUT out_vel", out_vel);
 
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_filter_tensor(DoubleTab& in_vel)
{
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  int nb_elem = domaine_VDF.domaine().nb_elem();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  int element_number;
  int num0, num1, num2, num3, num4, num5;
  int f0, f1, f2, f3, f4, f5;
  int i, j;
 
  DoubleTrav temp1(in_vel);
  DoubleTrav temp2(in_vel);
 
  // This is to calculate the filtered tensor field
 
  if (dimension == 2)
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 1);
          if ((num0 == -1) || (num2 == -1))
            for (i = 0; i < dimension; i++)
              for (j = 0; j < dimension; j++)
                temp1(element_number, i, j) = in_vel(element_number, i, j);
          else
            {
              for (i = 0; i < dimension; i++)
                for (j = 0; j < dimension; j++)
                  {
                    temp1(element_number, i, j) = 0.25 * (in_vel(num0, i, j) + 2.0 * in_vel(element_number, i, j) + in_vel(num2, i, j));
                  }
            }
        }
      temp1.echange_espace_virtuel();
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if ((num1 == -1) || (num3 == -1))
            for (i = 0; i < dimension; i++)
              for (j = 0; j < dimension; j++)
                in_vel(element_number, i, j) = temp1(element_number, i, j);
          else
            {
              for (i = 0; i < dimension; i++)
                for (j = 0; j < dimension; j++)
                  {
                    in_vel(element_number, i, j) = 0.25 * (temp1(num1, i, j) + 2.0 * temp1(element_number, i, j) + temp1(num3, i, j));
                  }
            }
        }
      in_vel.echange_espace_virtuel();
    }
  else
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if ((num0 == -1) || (num3 == -1))
            for (i = 0; i < dimension; i++)
              for (j = 0; j < dimension; j++)
                temp1(element_number, i, j) = in_vel(element_number, i, j);
          else
            {
              for (i = 0; i < dimension; i++)
                for (j = 0; j < dimension; j++)
                  {
                    temp1(element_number, i, j) = 0.25 * (in_vel(num0, i, j) + 2.0 * in_vel(element_number, i, j) + in_vel(num3, i, j));
                  }
            }
        }
      temp1.echange_espace_virtuel();
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          f4 = elem_faces(element_number, 4);
          num4 = face_voisins(f4, 1);
          if ((num1 == -1) || (num4 == -1))
            for (i = 0; i < dimension; i++)
              for (j = 0; j < dimension; j++)
                temp2(element_number, i, j) = temp1(element_number, i, j);
          else
            {
              for (i = 0; i < dimension; i++)
                for (j = 0; j < dimension; j++)
                  {
                    temp2(element_number, i, j) = 0.25 * (temp1(num1, i, j) + 2.0 * temp1(element_number, i, j) + temp1(num4, i, j));
                  }
            }
        }
      temp2.echange_espace_virtuel();
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 0);
          f5 = elem_faces(element_number, 5);
          num5 = face_voisins(f5, 1);
          if ((num2 == -1) || (num5 == -1))
            for (i = 0; i < dimension; i++)
              for (j = 0; j < dimension; j++)
                in_vel(element_number, i, j) = temp2(element_number, i, j);
          else
            {
              for (i = 0; i < dimension; i++)
                for (j = 0; j < dimension; j++)
                  {
                    in_vel(element_number, i, j) = 0.25 * (temp2(num2, i, j) + 2.0 * temp2(element_number, i, j) + temp2(num5, i, j));
                  }
            }
        }
      in_vel.echange_espace_virtuel();
    }
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Lij(const DoubleTab& cell_cent_vel, const DoubleTab& filt_vel, DoubleTab& Lij)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Lij INPUT cell_cent_vel", cell_cent_vel);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Lij INPUT filt_vel", filt_vel);
 
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  int element_number;
 
  DoubleTab uij_filt(Lij);
 
  // This is to calculate the Lij term for the C coefficient
  for (element_number = 0; element_number < nb_elem_tot; element_number++)
    {
      for (int i = 0; i < dimension; i++)
        for (int j = 0; j < dimension; j++)
          uij_filt(element_number, i, j) = cell_cent_vel(element_number, i) * cell_cent_vel(element_number, j);
    }
  calculer_filter_tensor(uij_filt);
 
  for (element_number = 0; element_number < nb_elem_tot; element_number++)
    {
      for (int i = 0; i < dimension; i++)
        for (int j = 0; j < dimension; j++)
          Lij(element_number, i, j) = uij_filt(element_number, i, j) - filt_vel(element_number, i) * filt_vel(element_number, j);
    }
 
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Lij OUTPUT Lij", Lij);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Mij(const DoubleTab& Sij_grid_scale, const DoubleTab& Sij_test_scale, const DoubleVect& l, DoubleTab& Mij)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Mij INPUT Sij_grid_scale", Sij_grid_scale);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Mij INPUT Sij_test_scale", Sij_test_scale);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Mij INPUT l", l);
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
 
  DoubleTrav Sij_grid_scale_norme(nb_elem_tot);
  DoubleTrav Sij_test_scale_norme(nb_elem_tot);
  DoubleTrav S_norme_Sij_filt(Mij);
  const double alpha = sqrt(6.0);
 
  // This is to calculate the Mij term for the C coefficient
 
  calculer_S_norme(Sij_grid_scale, Sij_grid_scale_norme, domaine_VDF.domaine().nb_elem_tot());
  calculer_S_norme(Sij_test_scale, Sij_test_scale_norme, domaine_VDF.domaine().nb_elem_tot());
  for (int element_number = 0; element_number < nb_elem_tot; element_number++)
    for (int i = 0; i < dimension; i++)
      for (int j = 0; j < dimension; j++)
        S_norme_Sij_filt(element_number, i, j) = Sij_grid_scale_norme(element_number) * Sij_grid_scale(element_number, i, j);
 
  calculer_filter_tensor(S_norme_Sij_filt);
 
  for (int element_number = 0; element_number < nb_elem_tot; element_number++)
    for (int i = 0; i < dimension; i++)
      for (int j = 0; j < dimension; j++)
        Mij(element_number, i, j) = 2 * l[element_number] * l[element_number]
                                    * (S_norme_Sij_filt(element_number, i, j) - alpha * alpha * Sij_test_scale_norme(element_number) * Sij_test_scale(element_number, i, j));
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Mij OUTPUT Mij", Mij);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_model_coefficient(const DoubleTab& Lij, const DoubleTab& Mij)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_model_coefficient INPUT Lij", Lij);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_model_coefficient INPUT Mij", Mij);
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  DoubleVect& model_coeff = coeff_field_->valeurs();
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  int nb_elem = domaine_VDF.domaine().nb_elem();
  DoubleTab haut, bas;
  domaine_VDF.domaine().creer_tableau_elements(haut);
  domaine_VDF.domaine().creer_tableau_elements(bas);
  // Evaluate the dynamic model coeficient C
  for (int elem = 0; elem < nb_elem; elem++)
    {
      haut(elem) = 0.;
      bas(elem) = 0.;
      for (int i = 0; i < dimension; i++)
        for (int j = 0; j < dimension; j++)
          {
            haut(elem) += Lij(elem, i, j) * Mij(elem, i, j);
            bas(elem) += Mij(elem, i, j) * Mij(elem, i, j);
          }
    }
  haut.echange_espace_virtuel();
  bas.echange_espace_virtuel();
 
  stabilise_moyenne(haut, bas);
 
  int i = 0;
  int j = 0;
  for (int elem = 0; elem < nb_elem; elem++)
    {
      if (model_coeff(elem) < 0.0)
        {
          model_coeff(elem) = 0.0;
          i++;
        }
      else if (model_coeff(elem) > 0.5)
        {
          // model_coeff(elem) = 0.5; // if I remove this line, no more seq/par diffs
          j++;
        }
    }
  model_coeff.echange_espace_virtuel();
 
  if (i != 0)
    {
      Journal() <<"Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_model_coefficient: "
                << i<<"/"<< nb_elem_tot
                << " elements ont un coefficient inferieur a 0.0 (modele dynamique)" << finl;
    }
  if (j != 0)
    {
      Journal() <<"Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_model_coefficient: "
                << j <<"/"<< nb_elem_tot
                << " elements ont un coefficient superieur a 0.5 (modele dynamique)" << finl;
    }
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_model_coefficient OUTPUT model_coeff", model_coeff);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_length_scale(DoubleVect& l, const Domaine_VDF& domaine_VDF)
{
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  const IntVect& orientation = domaine_VDF.orientation();
  int elem;
 
  // calcul de la longueur caracteristique de l'element
 
  for (elem = 0; elem < nb_elem_tot; elem++)
    {
      double l_int;
      if (dimension == 3)
        {
          l_int = domaine_VDF.dim_elem(elem, orientation(elem_faces(elem, 0)));
          l_int *= domaine_VDF.dim_elem(elem, orientation(elem_faces(elem, 1)));
          l_int *= domaine_VDF.dim_elem(elem, orientation(elem_faces(elem, 2)));
          l(elem) = exp((log(l_int) / 3.));
        }
      else
        {
          l_int = domaine_VDF.dim_elem(elem, orientation(elem_faces(elem, 0)));
          l_int *= domaine_VDF.dim_elem(elem, orientation(elem_faces(elem, 1)));
          l[elem] = sqrt(l_int);
        }
    }
  // sometimes it is passed as ditributed, sometimes not...
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_length_scale OUTPUT l", l);
}
 
Champ_Fonc_base& Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente(const DoubleVect& SMA_barre, const DoubleVect& l)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente INPUT SMA_barre", SMA_barre);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente INPUT l", l);
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  DoubleVect& model_coeff = coeff_field_->valeurs();
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente INPUT model_coeff", model_coeff);
  double temps = mon_equation_->inconnue().temps();
  DoubleTab& visco_turb = la_viscosite_turbulente_->valeurs();
  int nb_elem = domaine_VDF.domaine().nb_elem();
 
  if (visco_turb.size() != domaine_VDF.domaine().nb_elem())
    {
      Cerr << "Size error for the array containing the values of the turbulent viscosity." << finl;
      Process::exit();
    }
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente visco_turb 0", visco_turb);
 
  for (int elem = 0; elem < nb_elem; elem++)
    visco_turb[elem] = model_coeff[elem] * l[elem] * l[elem] * sqrt(SMA_barre[elem]);
 
 
  visco_turb.echange_espace_virtuel();
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente OUTPUT visco_turb", visco_turb);
 
  la_viscosite_turbulente_->changer_temps(temps);
  return la_viscosite_turbulente_;
}
 
Champ_Fonc_base& Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente()
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_viscosite_turbulente() OUTPUT la_viscosite_turbulente_", la_viscosite_turbulente_->valeurs());
  return la_viscosite_turbulente_;
}
 
Champ_Fonc_base& Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_energie_cinetique_turb(const DoubleVect& SMA_barre, const DoubleVect& l)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_energie_cinetique_turb IN SMA_barre", SMA_barre);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_energie_cinetique_turb IN l", l);
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  DoubleVect& model_coeff = coeff_field_->valeurs();
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_energie_cinetique_turb IN model_coeff", model_coeff);
  double temps = mon_equation_->inconnue().temps();
  DoubleVect& k = energie_cinetique_turb_->valeurs();
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  constexpr double co_mu = 1. / (0.094 * 0.094);
 
  for (int elem = 0; elem < nb_elem_tot; elem++)
    {
      k[elem] = co_mu * model_coeff[elem] * model_coeff[elem] * l[elem] * l[elem] * SMA_barre[elem];
    }
 
  energie_cinetique_turb_->changer_temps(temps);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_energie_cinetique_turb OUT energie_cinetique_turb_", k);
  return energie_cinetique_turb_;
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_energie_cinetique_turb()
{
 
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::controler_grandeurs_turbulentes()
{
  static const double Cmu = CMU;
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  DoubleTab& visco_turb = la_viscosite_turbulente_->valeurs();
  DoubleVect& energie_turb = energie_cinetique_turb_->valeurs();
  double k, eps;
 
  for (int elem = 0; elem < nb_elem_tot; elem++)
    {
      k = energie_turb[elem];
      if (k <= 1.e-10)
        energie_turb[elem] = 0.;
      else
        {
          eps = Cmu * k * k / (visco_turb[elem] + DMINFLOAT);
          if (eps < 1.e-15)
            visco_turb[elem] = 0;
        }
    }
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_S_norme(const DoubleTab& Sij, DoubleVect& S_norme, int nb_elem_tot)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_S_norme INPUT Sij",Sij);
  S_norme = 0.0;
  for (int element_number = 0; element_number < nb_elem_tot; element_number++)
    for (int i = 0; i < dimension; i++)
      for (int j = 0; j < dimension; j++)
        S_norme(element_number) += Sij(element_number, i, j) * Sij(element_number, i, j);
 
  for (int element_number = 0; element_number < nb_elem_tot; element_number++)
    S_norme(element_number) = sqrt(2 * S_norme(element_number));
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_S_norme OUTPUT S_norme",S_norme);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Sij(DoubleTab& Sij, const Domaine_VDF& domaine_VDF, const Domaine_Cl_VDF& domaine_Cl_VDF, Champ_Inc_base& inco)
{
  Champ_Face_VDF& vit = ref_cast(Champ_Face_VDF, inco);
  const DoubleTab& vitesse = inco.valeurs();
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Sij INPUT vitesse",vitesse);
 
  const int nb_elem_tot = domaine_VDF.nb_elem_tot();
  const int nb_elem = domaine_VDF.nb_elem();
  assert(vitesse.line_size() == 1);
  DoubleTab gij(nb_elem_tot, dimension, dimension, vitesse.line_size());
 
  vit.calcul_duidxj(vitesse, gij, domaine_Cl_VDF);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Sij MID gij",gij);
 
  // Calcul de Sij
  for (int elem = 0; elem < nb_elem; elem++)
    for (int i = 0; i < dimension; i++)
      for (int j = 0; j < dimension; j++)
        Sij(elem, i, j) = 0.5 * (gij(elem, i, j, 0) + gij(elem, j, i, 0));
  Sij.echange_espace_virtuel();
 
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Sij OUTPUT Sij",Sij);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Sij_vel_filt(const DoubleTab& in_vel, DoubleTab& out_vel, const Domaine_VDF& domaine_VDF)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Sij_vel_filt INPUT in_vel",in_vel);
  int nb_elem = domaine_VDF.domaine().nb_elem();
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
 
  int element_number;
  int num0, num1, num2, num3, num4, num5;
  int f0, f1, f2, f3, f4, f5;
 
  if (dimension == 2)
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          if (num0 == -1)
            num0 = element_number;
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          if (num1 == -1)
            num1 = element_number;
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 1);
          if (num2 == -1)
            num2 = element_number;
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if (num3 == -1)
            num3 = element_number;
          out_vel(element_number, 0, 0) = 0.5 * ((in_vel(num2, 0) - in_vel(num0, 0)) / domaine_VDF.dim_elem(element_number, 0));
          out_vel(element_number, 0, 1) = out_vel(element_number, 1, 0) = 0.25
                                                                          * ((in_vel(num3, 0) - in_vel(num1, 0)) / domaine_VDF.dim_elem(element_number, 1) + (in_vel(num2, 1) - in_vel(num0, 1)) / domaine_VDF.dim_elem(element_number, 0));
          out_vel(element_number, 1, 1) = 0.5 * ((in_vel(num3, 1) - in_vel(num1, 1)) / domaine_VDF.dim_elem(element_number, 1));
        }
    }
  else
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          if (num0 == -1)
            num0 = element_number;
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          if (num1 == -1)
            num1 = element_number;
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 0);
          if (num2 == -1)
            num2 = element_number;
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if (num3 == -1)
            num3 = element_number;
          f4 = elem_faces(element_number, 4);
          num4 = face_voisins(f4, 1);
          if (num4 == -1)
            num4 = element_number;
          f5 = elem_faces(element_number, 5);
          num5 = face_voisins(f5, 1);
          if (num5 == -1)
            num5 = element_number;
          out_vel(element_number, 0, 0) = 0.5 * ((in_vel(num3, 0) - in_vel(num0, 0)) / domaine_VDF.dim_elem(element_number, 0));
          out_vel(element_number, 0, 1) = out_vel(element_number, 1, 0) = 0.25
                                                                          * ((in_vel(num4, 0) - in_vel(num1, 0)) / domaine_VDF.dim_elem(element_number, 1) + (in_vel(num3, 1) - in_vel(num0, 1)) / domaine_VDF.dim_elem(element_number, 0));
          out_vel(element_number, 0, 2) = out_vel(element_number, 2, 0) = 0.25
                                                                          * ((in_vel(num5, 0) - in_vel(num2, 0)) / domaine_VDF.dim_elem(element_number, 2) + (in_vel(num3, 2) - in_vel(num0, 2)) / domaine_VDF.dim_elem(element_number, 0));
          out_vel(element_number, 1, 1) = 0.5 * ((in_vel(num4, 1) - in_vel(num1, 1)) / domaine_VDF.dim_elem(element_number, 1));
          out_vel(element_number, 1, 2) = out_vel(element_number, 2, 1) = 0.25
                                                                          * ((in_vel(num5, 1) - in_vel(num2, 1)) / domaine_VDF.dim_elem(element_number, 2) + (in_vel(num4, 2) - in_vel(num1, 2)) / domaine_VDF.dim_elem(element_number, 1));
          out_vel(element_number, 2, 2) = 0.5 * ((in_vel(num5, 2) - in_vel(num2, 2)) / domaine_VDF.dim_elem(element_number, 2));
        }
    }
  out_vel.echange_espace_virtuel();
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calculer_Sij_vel_filt OUTPUT out_vel",out_vel);
}
 
//////////////////////////////////////////////////////////////////////
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne(const DoubleTab& haut, const DoubleTab& bas)
{
  if (methode_stabilise_ == "plans_paralleles")
    stabilise_moyenne_plans_paralleles(haut, bas);
  else if (methode_stabilise_ == "6_points")
    stabilise_moyenne_6_points(haut, bas);
  else if (methode_stabilise_ == "moy_euler" || methode_stabilise_ == "moy_lagrange")
    stabilise_moyenne_euler_lagrange(haut, bas);
  else
    {
      Cerr << "Problem with the definition of the stabilisation method." << finl;
      Process::exit();
    }
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_6_points(const DoubleTab& haut, const DoubleTab& bas)
{
  DoubleVect& model_coeff = coeff_field_->valeurs();
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  int nb_elem = domaine_VDF.domaine().nb_elem();
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  int element_number;
  int num0, num1, num2, num3, num4, num5;
  int f0, f1, f2, f3, f4, f5;
 
  DoubleTrav temp1(model_coeff);
  DoubleTrav temp2(model_coeff);
  // Evaluate the dynamic model coeficient C
  for (element_number = 0; element_number < nb_elem_tot; element_number++)
    {
      if (std::fabs(bas(element_number)) < 1.e-12)
        model_coeff[element_number] = 0.;
      else
        model_coeff[element_number] = haut[element_number] / bas[element_number];
    }
  // Stabilisation du modele par moyennage sur 6 noeuds
  if (dimension == 2)
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 1);
          if ((num0 == -1) || (num2 == -1))
            temp1(element_number) = 3.0 * model_coeff(element_number);
          else
            {
              temp1(element_number) = model_coeff(num0) + model_coeff(element_number) + model_coeff(num2);
            }
        }
      temp1.echange_espace_virtuel();
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if ((num1 == -1) || (num3 == -1))
            model_coeff(element_number) = temp1(element_number) / 3.0;
          else
            {
              model_coeff(element_number) = (temp1(num1) + temp1(element_number) + temp1(num3)) / 9.0;
            }
        }
    }
  else
    {
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f0 = elem_faces(element_number, 0);
          num0 = face_voisins(f0, 0);
          f3 = elem_faces(element_number, 3);
          num3 = face_voisins(f3, 1);
          if ((num0 == -1) || (num3 == -1))
            temp1(element_number) = 3.0 * model_coeff(element_number);
          else
            {
              temp1(element_number) = model_coeff(num0) + model_coeff(element_number) + model_coeff(num3);
            }
        }
      temp1.echange_espace_virtuel();
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f1 = elem_faces(element_number, 1);
          num1 = face_voisins(f1, 0);
          f4 = elem_faces(element_number, 4);
          num4 = face_voisins(f4, 1);
          if ((num1 == -1) || (num4 == -1))
            temp2(element_number) = 3.0 * temp1(element_number);
          else
            {
              temp2(element_number) = temp1(num1) + temp1(element_number) + temp1(num4);
            }
        }
      temp2.echange_espace_virtuel();
      for (element_number = 0; element_number < nb_elem; element_number++)
        {
          f2 = elem_faces(element_number, 2);
          num2 = face_voisins(f2, 0);
          f5 = elem_faces(element_number, 5);
          num5 = face_voisins(f5, 1);
          if ((num2 == -1) || (num5 == -1))
            model_coeff(element_number) = temp2(element_number) / 9.0;
          else
            {
              model_coeff(element_number) = (temp2(num2) + temp2(element_number) + temp2(num5)) / 27.0;
            }
        }
    }
  model_coeff.echange_espace_virtuel();
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_plans_paralleles(const DoubleTab& haut, const DoubleTab& bas)
{
  DoubleVect& model_coeff = coeff_field_->valeurs();
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  int nb_elem = domaine_VDF.domaine().nb_elem();
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  DoubleVect model_coeff_tmp(N_c_);
 
  // Evaluate the dynamic model coeficient C
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_plans_paralleles model_coeff avant haut/bas", model_coeff);
  for (int element_number = 0; element_number < nb_elem_tot; element_number++)
    {
      if (std::fabs(bas(element_number)) < 1.e-12)
        model_coeff[element_number] = 0.;
      else
        model_coeff[element_number] = haut[element_number] / bas[element_number];
    }
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_plans_paralleles model_coeff apres haut/bas", model_coeff);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_plans_paralleles corresp_c apres haut/bas", corresp_c_);
  for (int j = 0; j < N_c_; j++)
    model_coeff_tmp(j) = 0.0;
 
  for (int element_number = 0; element_number < nb_elem; element_number++)
    model_coeff_tmp(corresp_c_(element_number)) += model_coeff(element_number);
 
  for (int j = 0; j < N_c_; j++)
    {
      model_coeff_tmp(j) = Process::mp_sum(model_coeff_tmp(j));
      model_coeff_tmp(j) /= Process::check_int_overflow(Process::mp_sum(compt_c_(j)));
    }
  for (int element_number = 0; element_number < nb_elem; element_number++)
    model_coeff(element_number) = model_coeff_tmp(corresp_c_(element_number));
  model_coeff.echange_espace_virtuel();
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_euler_lagrange(const DoubleTab& haut, const DoubleTab& bas)
{
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_euler_lagrange INPUT haut", haut);
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_euler_lagrange INPUT bas", bas);
  static int init = 1;
  // Attention: on stocke ces deux tableaux d'un appel a l'autre:
  static DoubleTrav haut_moy;
  static DoubleTrav bas_moy;
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  DoubleVect& model_coeff = coeff_field_->valeurs();
 
 
 
  if (init == 1)
    {
      init = 0;
      haut_moy = haut;
      bas_moy = bas;
      model_coeff = 0.16 * 0.16;
    }
  else
    {
      const bool lagrange = (methode_stabilise_ == "moy_lagrange");
 
      int nb_0 = 0;
      DoubleVect l;
      domaine_VDF.domaine().creer_tableau_elements(l);
      calculer_length_scale(l, domaine_VDF);
 
      const int nmax = haut.size_totale();
      const double dt = mon_equation_->schema_temps().pas_de_temps();
 
      DoubleVect dist(8);
      IntVect num(8);
 
      for (int elem = 0; elem < nmax; elem++)
        {
          const double hmoy = haut_moy(elem);
          const double bmoy = bas_moy(elem);
          double hmoy_int, bmoy_int;
          if (lagrange)
            {
              calcul_voisins(elem, num, dist);
              interpole(num, dist, haut_moy, hmoy_int);
              interpole(num, dist, bas_moy, bmoy_int);
            }
          else
            {
              hmoy_int = hmoy;
              bmoy_int = bmoy;
            }
          const double h = haut(elem);
          const double b = bas(elem);
          const double l_elem = l(elem);
          const double a = hmoy * bmoy;
          double eps = 1.0; // default value
          if (a > 0.)
            {
              const double T = 1.5 * l_elem * pow(a, -0.125);
              const double dt_T = dt / T;
              eps = dt_T / (1. + dt_T);
            }
 
          double hmoy_new = eps * h + (1 - eps) * hmoy_int;
          if (hmoy_new < 0.)
            {
              hmoy_new = 0.0;
              nb_0++;
            }
          double bmoy_new = eps * b + (1 - eps) * bmoy_int;
          double model_c = 0.16 * 0.16; // default value
          if (bmoy_new != 0.)
            model_c = hmoy_new / bmoy_new;
 
          model_coeff(elem) = model_c;
          haut_moy(elem) = hmoy_new;
          bas_moy(elem) = bmoy_new;
 
        }
      if (nb_0 != 0)
        {
          Journal() << "Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_euler_lagrange" << finl;
          Journal() << nb_0 << " elements have a negative coefficient." << finl;
        }
    }
 
  model_coeff.echange_espace_virtuel();
  Debog::verifier("Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::stabilise_moyenne_euler_lagrange OUTPUT model_coeff", model_coeff);
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calcul_voisins(const int element_number, IntVect& num, DoubleVect& dist)
{
  double dt = mon_equation_->schema_temps().pas_de_temps();
  int i, j;
  IntTab indice(3);
 
  if (elem_elem_(element_number, 1, 1, 1) == -1)
    {
      for (j = 0; j < 8; j++)
        {
          num[j] = element_number;
          dist[j] = 1.0;
        }
    }
  else
    {
      for (i = 0; i < 3; i++)
        {
          if (cell_cent_vel_(element_number, i) > 0)
            indice(i) = 2;
          else
            indice(i) = 0;
        }
      num[0] = elem_elem_(element_number, 1, 1, 1);
      num[1] = elem_elem_(element_number, 1, 1, indice(2));
      num[2] = elem_elem_(element_number, 1, indice(1), 1);
      num[3] = elem_elem_(element_number, 1, indice(1), indice(2));
      num[4] = elem_elem_(element_number, indice(0), 1, 1);
      num[5] = elem_elem_(element_number, indice(0), 1, indice(2));
      num[6] = elem_elem_(element_number, indice(0), indice(1), 1);
      num[7] = elem_elem_(element_number, indice(0), indice(1), indice(2));
 
      const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
      const DoubleTab& xp = domaine_VDF.xp();
      DoubleTrav x(3);
      for (i = 0; i < 3; i++)
        x(i) = xp(element_number, i) - cell_cent_vel_(element_number, i) * dt;
      dist = 0.0;
      for (j = 0; j < 8; j++)
        {
          for (i = 0; i < 3; i++)
            {
              if (num[j] == -1)
                {
                  Cerr << "Problem with the algorithm Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calcul_voisins" << finl;
                  Process::exit();
                }
              dist(j) += (xp(num[j], i) - x(i)) * (xp(num[j], i) - x(i));
            }
          dist(j) = sqrt(dist(j));
        }
 
      // verification que l'on ne se trouve pas sur un element
      double d;
      d = (xp(num[0], 0) - xp(num[4], 0)) * (xp(num[0], 0) - xp(num[4], 0));
      d += (xp(num[0], 1) - xp(num[2], 1)) * (xp(num[0], 1) - xp(num[2], 1));
      d += (xp(num[0], 2) - xp(num[1], 2)) * (xp(num[0], 2) - xp(num[1], 2));
      int flag = 0;
      for (j = 0; j < 8; j++)
        if (dist(j) * dist(j) < d / 10000)
          flag = 1;
      if (flag == 1)
        for (j = 0; j < 8; j++)
          {
            num[j] = element_number;
            dist[j] = 1.0;
          }
    }
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::interpole(const IntVect& num, const DoubleVect& dist, const DoubleVect& val, double& val_int)
{
  double K = 0.0;
  val_int = 0.0;
  for (int i = 0; i < 8; i++)
    {
      K += 1 / dist(i);
      val_int += val(num(i)) / dist(i);
    }
  val_int /= K;
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calc_elem_elem()
{
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  const IntTab& face_voisins = domaine_VDF.face_voisins();
  const IntTab& elem_faces = domaine_VDF.elem_faces();
  int nb_elem_tot = domaine_VDF.domaine().nb_elem_tot();
  int element_number, f, elem;
 
  //elem_elem(element_number,i,j,k) i-x, j-y k-z et 0-avant 1-meme_position, 2-apres
 
  elem_elem_.resize(nb_elem_tot, 3, 3, 3);
 
  for (element_number = 0; element_number < nb_elem_tot; element_number++)
    {
      elem_elem_(element_number, 1, 1, 1) = element_number;  // xyz
      f = elem_faces(element_number, 0);
      elem_elem_(element_number, 0, 1, 1) = face_voisins(f, 0);  // x-
      elem = elem_elem_(element_number, 0, 1, 1);
      if (elem != -1)
        {
          f = elem_faces(elem, 1);
          elem_elem_(element_number, 0, 0, 1) = face_voisins(f, 0); // x- y-
          f = elem_faces(elem, 2);
          elem_elem_(element_number, 0, 1, 0) = face_voisins(f, 0); // x- z-
          f = elem_faces(elem, 4);
          elem_elem_(element_number, 0, 2, 1) = face_voisins(f, 1); // x- y+
          f = elem_faces(elem, 5);
          elem_elem_(element_number, 0, 1, 2) = face_voisins(f, 1); // x- z+
        }
      f = elem_faces(element_number, 3);
      elem_elem_(element_number, 2, 1, 1) = face_voisins(f, 1);  // x+
      elem = elem_elem_(element_number, 2, 1, 1);
      if (elem != -1)
        {
          f = elem_faces(elem, 1);
          elem_elem_(element_number, 2, 0, 1) = face_voisins(f, 0); // x+ y-
          f = elem_faces(elem, 2);
          elem_elem_(element_number, 2, 1, 0) = face_voisins(f, 0); // x+ z-
          f = elem_faces(elem, 4);
          elem_elem_(element_number, 2, 2, 1) = face_voisins(f, 1); // x+ y+
          f = elem_faces(elem, 5);
          elem_elem_(element_number, 2, 1, 2) = face_voisins(f, 1); // x+ z+
        }
      f = elem_faces(element_number, 1);
      elem_elem_(element_number, 1, 0, 1) = face_voisins(f, 1);  // y-
      elem = elem_elem_(element_number, 1, 0, 1);
      if (elem != -1)
        {
          f = elem_faces(elem, 2);
          elem_elem_(element_number, 1, 0, 0) = face_voisins(f, 0); // y- z-
          f = elem_faces(elem, 5);
          elem_elem_(element_number, 1, 0, 2) = face_voisins(f, 1); // y- z+
        }
      f = elem_faces(element_number, 4);
      elem_elem_(element_number, 1, 2, 1) = face_voisins(f, 1);  // y+
      elem = elem_elem_(element_number, 1, 2, 1);
      if (elem != -1)
        {
          f = elem_faces(elem, 2);
          elem_elem_(element_number, 1, 2, 0) = face_voisins(f, 0); // y+ z-
          f = elem_faces(elem, 5);
 
          elem_elem_(element_number, 1, 2, 2) = face_voisins(f, 1); // y+ z+
        }
      f = elem_faces(element_number, 2);
      elem_elem_(element_number, 1, 1, 0) = face_voisins(f, 0);  // z-
      f = elem_faces(element_number, 5);
      elem_elem_(element_number, 1, 1, 2) = face_voisins(f, 1);  // z+
      elem = elem_elem_(element_number, 1, 2, 2);
      if (elem != -1)
        {
          f = elem_faces(elem, 0);
          elem_elem_(element_number, 0, 2, 2) = face_voisins(f, 0); // x- y+ z+
          f = elem_faces(elem, 3);
          elem_elem_(element_number, 2, 2, 2) = face_voisins(f, 1); // x+ y+ z+
        }
      elem = elem_elem_(element_number, 1, 2, 0);
      if (elem != -1)
        {
          f = elem_faces(elem, 0);
          elem_elem_(element_number, 0, 2, 0) = face_voisins(f, 0); // x- y+ z-
          f = elem_faces(elem, 3);
          elem_elem_(element_number, 2, 2, 0) = face_voisins(f, 1); // x+ y+ z-
        }
      elem = elem_elem_(element_number, 1, 0, 2);
      if (elem != -1)
        {
          f = elem_faces(elem, 0);
          elem_elem_(element_number, 0, 0, 2) = face_voisins(f, 0); // x- y- z+
          f = elem_faces(elem, 3);
          elem_elem_(element_number, 2, 0, 2) = face_voisins(f, 1); // x+ y- z+
        }
      elem = elem_elem_(element_number, 1, 0, 0);
      if (elem != -1)
        {
          f = elem_faces(elem, 0);
          elem_elem_(element_number, 0, 0, 0) = face_voisins(f, 0); // x- y- z-
          f = elem_faces(elem, 3);
          elem_elem_(element_number, 2, 0, 0) = face_voisins(f, 1); // x+ y- z-
        }
 
      elem_elem_(element_number, 1, 1, 1) = 0;
 
      //Version initiale, boucles ?????
      /*
       for (int i=0; i<3; i++)
       for (int j=0; i<3; i++)
       for (int k=0; i<3; i++)
       if (elem_elem(element_number,i,j,k)==-1)
       elem_elem(element_number,1,1,1)=-1;
       */
      //A remplacer par ????
      for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
          for (int k = 0; k < 3; k++)
            if (elem_elem_(element_number, i, j, k) == -1)
              elem_elem_(element_number, 1, 1, 1) = -1;
 
    }
 
}
 
void Modele_turbulence_hyd_LES_SMAGO_DYN_VDF::calcul_tableaux_correspondance(int& N_c, IntVect& compt_c, IntVect& corresp_c)
{
  // Initialisation de : Yuv + compt_c + corresp_c
  const Domaine_VDF& domaine_VDF = ref_cast(Domaine_VDF, le_dom_VF_.valeur());
  const DoubleTab& xp = domaine_VDF.xp();
  int nb_elems = domaine_VDF.domaine().nb_elem();
  DoubleTrav Y_c;
  int num_elem, j, indic_c, trouve;
  double y;
  j = 0;
  indic_c = 0;
 
  // dimensionnement aux valeurs rentrees dans le jeu de donnees
  Y_c.resize(N_c);
  compt_c.resize(N_c);
  domaine_VDF.domaine().creer_tableau_elements(corresp_c);
 
  // initialisation
  Y_c = -100.;
  compt_c = 0;
  corresp_c = -1;
 
  // Tolerance relative pour la comparaison des coordonnees y.
  // Evite les faux plans fantomes dus aux arrondis flottants en parallele
  // sur maillages non uniformes (deadlock MPI dans stabilise_moyenne_plans_paralleles).
  const double tol_rel = 1.e-8;
 
  for (num_elem = 0; num_elem < nb_elems; num_elem++)
    {
      y = xp(num_elem, 1);
      trouve = 0;
 
      for (j = 0; j < indic_c + 1; j++)
        {
          if (j == N_c)
            {
              Cerr << "nb_points=" << N_c << " is too low  for the parallel planes number." << finl;
              Cerr << "Please increase this value." << finl;
              Process::exit();
            }
          // FIX: comparaison avec tolerance relative au lieu de l'egalite exacte
          // sur doubles, qui causait des plans fantomes en parallele sur maillages
          // non uniformes, entrainant un deadlock MPI dans mp_sum.
          const double ref = std::fabs(y) + std::fabs(Y_c[j]) + 1.e-30;
          if (std::fabs(y - Y_c[j]) < tol_rel * ref)
            {
              corresp_c[num_elem] = j;
              compt_c[j]++;
              j = indic_c + 1;
              trouve = 1;
              break;
            }
        }
      if (trouve == 0)
        {
          corresp_c[num_elem] = indic_c;
          Y_c[indic_c] = y;
          compt_c[indic_c]++;
          indic_c++;
        }
    }
  N_c = indic_c;  // nombre de y pour les elements locaux
 
  // FIX: synchronisation de N_c entre tous les processus MPI.
  // Sans cela, N_c peut differer d'un processus a l'autre (plans fantomes),
  // ce qui provoque un deadlock dans la boucle mp_sum de
  // stabilise_moyenne_plans_paralleles ou le collectif n'est pas appele
  // le meme nombre de fois par tous les processus.
  N_c = Process::mp_max(N_c);
 
  compt_c.resize(N_c);
}