Modele_turbulence_hyd_base.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 <Modele_turbulence_hyd_base.h>
#include <EcritureLectureSpecial.h>
#include <Operateur_Diff_base.h>
#include <Operateur_Conv_base.h>
#include <Discretisation_base.h>
#include <Schema_Temps_base.h>
#include <communications.h>
#include <LecFicDiffuse.h>
#include <EcrFicPartage.h>
#include <Postraitement.h>
#include <Equation_base.h>
#include <Probleme_base.h>
#include <TRUST_2_PDI.h>
#include <Operateur.h>
#include <Domaine.h>
#include <Param.h>
#include <Debog.h>
 
Implemente_base(Modele_turbulence_hyd_base, "Modele_turbulence_hyd_base", Objet_U);
// XD modele_turbulence_hyd_deriv objet_lecture modele_turbulence_hyd_deriv INHERITS_BRACE Basic class for turbulence
// XD_CONT model for Navier-Stokes equations.
 
Sortie& Modele_turbulence_hyd_base::printOn(Sortie& s) const { return s << que_suis_je() << " " << le_nom(); }
 
/*! @brief Reads the specifications of a turbulence model from an input stream.
 *
 *     Format:
 *     Modele_turbulence type_modele
 *     {
 *       [Turbulence_paroi la_loi_de_paroi]
 *       reading block specific to the turbulence type
 *     }
 *
 * @param (Entree& is) an input stream
 * @return (Entree&) the modified input stream
 * @throws opening brace expected
 * @throws wall law incompatible with the problem
 */
Entree& Modele_turbulence_hyd_base::readOn(Entree& is)
{
  Cerr << "Reading of data for a " << que_suis_je() << " hydraulic turbulence model" << finl;
  Param param(que_suis_je());
  set_param(param);
  param.lire_avec_accolades_depuis(is);
  // Checks
  // Case where Correction_visco_turb_pour_controle_pas_de_temps_parametre is read without Correction_visco_turb_pour_controle_pas_de_temps
  if (dt_diff_sur_dt_conv_ != -1)
    calcul_borne_locale_visco_turb_ = true;
  // Case where Correction_visco_turb_pour_controle_pas_de_temps is set without Correction_visco_turb_pour_controle_pas_de_temps_parametre
  else if (calcul_borne_locale_visco_turb_)
    dt_diff_sur_dt_conv_ = 1;
  // Case where Correction_visco_turb_pour_controle_pas_de_temps_parametre is incorrectly set
  if (dt_diff_sur_dt_conv_ <= 0 && calcul_borne_locale_visco_turb_)
    {
      Cerr << "Parameter value Correction_visco_turb_pour_controle_pas_de_temps_parametre must be strictly positive." << finl;
      exit();
    }
  return is;
}
 
void Modele_turbulence_hyd_base::set_param(Param& param) const
{
  param.ajouter_non_std("turbulence_paroi", (this), Param::REQUIRED); // XD attr turbulence_paroi turbulence_paroi_base turbulence_paroi OPT Keyword to set the wall law.
  param.ajouter_non_std("dt_impr_ustar", (this));                     // XD attr dt_impr_ustar floattant dt_impr_ustar OPT This keyword is used to print the values (U +, d+, u$\star$) obtained with the wall laws into a file named datafile_ProblemName_Ustar.face and periode refers to the printing period, this value is expressed in seconds.
  param.ajouter_non_std("dt_impr_ustar_mean_only", (this));           // XD attr dt_impr_ustar_mean_only dt_impr_ustar_mean_only dt_impr_ustar_mean_only OPT This keyword is used to print the mean values of u* ( obtained with the wall laws) on each boundary, into a file named datafile_ProblemName_Ustar_mean_only.out. periode refers to the printing period, this value is expressed in seconds. If you don\'t use the optional keyword boundaries, all the boundaries will be considered. If you use it, you must specify nb_boundaries which is the number of boundaries on which you want to calculate the mean values of u*, then you have to specify their names.
  param.ajouter("nut_max", &XNUTM_);                                  // XD attr nut_max floattant nut_max OPT Upper
  // XD_CONT limitation of turbulent viscosity (default value 1.e8).
  param.ajouter_flag("Correction_visco_turb_pour_controle_pas_de_temps", &calcul_borne_locale_visco_turb_); // XD attr correction_visco_turb_pour_controle_pas_de_temps rien correction_visco_turb_pour_controle_pas_de_temps OPT Keyword to set a limitation to low time steps due to high values of turbulent viscosity. The limit for turbulent viscosity is calculated so that diffusive time-step is equal or higher than convective time-step. For a stationary flow, the correction for turbulent viscosity should apply only during the first time steps and not when permanent state is reached. To check that, we could post process the corr_visco_turb field which is the correction of turbulent viscosity: it should be 1. on the whole domain.
  param.ajouter("Correction_visco_turb_pour_controle_pas_de_temps_parametre", &dt_diff_sur_dt_conv_);       // XD attr correction_visco_turb_pour_controle_pas_de_temps_parametre floattant correction_visco_turb_pour_controle_pas_de_temps_parametre OPT Keyword to set a limitation to low time steps due to high values of turbulent viscosity. The limit for turbulent viscosity is the ratio between diffusive time-step and convective time-step is higher or equal to the given value [0-1]
  //param.ajouter_condition("not(is_read_dt_impr_ustar_mean_only_and_is_read_dt_impr_ustar)","only one of dt_impr_ustar_mean_only and dt_impr_ustar can be used");
}
 
int Modele_turbulence_hyd_base::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  Motcle motlu;
  int retval = 1;
  if (mot == "turbulence_paroi")
    {
      Turbulence_paroi_base::typer_lire_turbulence_paroi(loipar_, *this, is);
      is >> loipar_.valeur();
    }
  else if (mot == "dt_impr_ustar")
    {
      if (loipar_->que_suis_je().contient("negligeable"))
        Process::exit("Please remove dt_impr_ustar option if the wall law is of Negligeable type.");
      is >> dt_impr_ustar_;
    }
  else if (mot == "dt_impr_ustar_mean_only")
    {
      if (loipar_->que_suis_je().contient("negligeable"))
        Process::exit("Please remove dt_impr_ustar option if the wall law is of Negligeable type.");
      // XD dt_impr_ustar_mean_only objet_lecture nul BRACE not_set
      // XD attr dt_impr floattant dt_impr REQ not_set
      // XD attr boundaries listchaine boundaries OPT not_set
      Nom accolade_ouverte = "{";
      Nom accolade_fermee = "}";
      nom_fichier_ = Objet_U::nom_du_cas() + "_" + equation().probleme().le_nom() + "_ustar_mean_only";
      Domaine& dom = equation().probleme().domaine();
      LIST(Nom) nlistbord_dom;                      //!< list storing all boundary names of the domain
      int nbfr = dom.nb_front_Cl();
      for (int b = 0; b < nbfr; b++)
        {
          Frontiere& org = dom.frontiere(b);
          nlistbord_dom.add(org.le_nom());
        }
      is >> motlu;
      if (motlu != accolade_ouverte)
        {
          Cerr << motlu << " is not a keyword understood by " << que_suis_je() << " in lire_motcle_non_standard" << finl;
          Cerr << "A specification of kind : dt_impr_ustar_mean_only { dt_impr periode [boundaries nb_boundaries boundary_name1 boundary_name2 ... ] } was expected." << finl;
          exit();
        }
      is >> motlu;
      if (motlu != "dt_impr")
        {
          Cerr << "We expected dt_impr..." << finl;
          exit();
        }
      is >> dt_impr_ustar_mean_only_;
 
      is >> motlu; // boundaries, closing brace, or unknown token
      if (motlu != accolade_fermee)
        {
          if (motlu == "boundaries")
            {
              boundaries_ = 1;
              int nb_bords = 0;
              Nom nom_bord_lu;
 
              // read boundaries number
              is >> nb_bords;
              if (nb_bords != 0)
                {
                  // read boundaries
                  for (int i = 0; i < nb_bords; i++)
                    {
                      is >> nom_bord_lu;
                      boundaries_list_.add(Nom(nom_bord_lu));
                      //  check boundary name
                      if (!nlistbord_dom.contient(boundaries_list_[i]))
                        {
                          Cerr << "Problem in the dt_impr_ustar_mean_only instruction:" << finl;
                          Cerr << "The boundary named '" << boundaries_list_[i] << "' is not a boundary of the domain " << dom.le_nom() << "." << finl;
                          exit();
                        }
                    }
                }
              // read closing brace
              is >> motlu;
              if (motlu != accolade_fermee)
                {
                  Cerr << "Problem in the dt_impr_ustar_mean_only instruction:" << finl;
                  Cerr << "TRUST wants to read a '" << accolade_fermee << "' but find '" << motlu << "'!!" << finl;
                  exit();
                }
            }
          else
            {
              Cerr << motlu << " is not a keyword understood by " << que_suis_je() << " in lire_motcle_non_standard" << finl;
              Cerr << "A specification of kind : dt_impr_ustar_mean_only { dt_impr periode [boundaries nb_boundaries boundary_name1 boundary_name2 ... ] } was expected." << finl;
              exit();
            }
        }
    } // end dt_impr_ustar_mean_only
  else
    retval = -1;
 
  return retval;
}
 
/*! @brief Associates the equation passed as parameter to the turbulence model.
 *
 * @param eqn The equation to associate with this object.
 */
void Modele_turbulence_hyd_base::associer_eqn(const Equation_base& eqn)
{
  mon_equation_ = eqn;
}
 
/*! @brief Reads the domain Wall_length.xyz file to populate the wall_length field for distance_paroi post-processing.
 *
 */
void Modele_turbulence_hyd_base::lire_distance_paroi()
{
  // PQ : 25/02/04 retrieval of the wall distance from Wall_length.xyz
  DoubleTab& wall_length = wall_length_->valeurs();
  wall_length = -1.;
 
  LecFicDiffuse fic;
  Nom nom_fic = equation().probleme().domaine().le_nom() + "_Wall_length.xyz";
  fic.set_bin(1);
  if (!fic.ouvrir(nom_fic))
    {
      Cerr << " File " << nom_fic << " doesn't exist. To generate it, please, refer to html.doc (Distance_paroi) " << finl;
      exit();
    }
 
  Noms nom_paroi;
 
  fic >> nom_paroi;
 
  EcritureLectureSpecial::lecture_special(wall_length_, fic);
 
}
 
/*! @brief Discretizes the turbulence model.
 *
 */
void Modele_turbulence_hyd_base::discretiser()
{
  Cerr << "Turbulence hydraulic model discretization" << finl;
  discretiser_visc_turb(mon_equation_->schema_temps(), mon_equation_->domaine_dis(), la_viscosite_turbulente_);
  champs_compris_.ajoute_champ(la_viscosite_turbulente_);
 
  discretiser_corr_visc_turb(mon_equation_->schema_temps(), mon_equation_->domaine_dis(), corr_visco_turb_);
  champs_compris_.ajoute_champ(corr_visco_turb_);
  const Discretisation_base& dis = ref_cast(Discretisation_base, mon_equation_->discretisation());
  dis.discretiser_champ("champ_elem", mon_equation_->domaine_dis(), "distance_paroi", "m", 1, mon_equation_->schema_temps().temps_courant(), wall_length_);
  champs_compris_.ajoute_champ(wall_length_);
}
 
void Modele_turbulence_hyd_base::discretiser_visc_turb(const Schema_Temps_base& sch, Domaine_dis_base& z, OWN_PTR(Champ_Fonc_base)& ch) const
{
  int is_dilat = equation().probleme().is_dilatable();
  if (!is_dilat)
    Cerr << "Turbulent viscosity field discretization" << finl;
  else
    Cerr << "Turbulent dynamic viscosity field discretization" << finl;
  Nom nom = (is_dilat == 1) ? "viscosite_dynamique_turbulente" : "viscosite_turbulente";
  Nom unite = (is_dilat == 1) ? "kg/(m.s)" : "m2/s";
  const Discretisation_base& dis = mon_equation_->discretisation();
  dis.discretiser_champ("champ_elem", z, nom, unite, 1, sch.temps_courant(), ch);
}
 
void Modele_turbulence_hyd_base::discretiser_corr_visc_turb(const Schema_Temps_base& sch, Domaine_dis_base& z, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Turbulent viscosity correction field discretization" << finl;
  const Discretisation_base& dis = mon_equation_->discretisation();
  dis.discretiser_champ("champ_elem", z, "corr_visco_turb", "adimensionnel", 1, sch.temps_courant(), ch);
}
 
void Modele_turbulence_hyd_base::discretiser_K(const Schema_Temps_base& sch, Domaine_dis_base& z, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Kinetic turbulent energy field discretisation" << finl;
  const Discretisation_base& dis = mon_equation_->discretisation();
  dis.discretiser_champ("champ_elem", z, "K", "m2/s2", 1, sch.temps_courant(), ch);
}
 
/*! @brief Prepares the computation.
 *
 * Initializes the wall law.
 *     Fills the wall_length field for distance_paroi post-processing.
 *
 * @return Return code from Turbulence_paroi::init_lois_paroi().
 */
int Modele_turbulence_hyd_base::preparer_calcul()
{
  int res = 1;
  if (loipar_)
    res = loipar_->init_lois_paroi();
 
  bool contient_distance_paroi = false;
 
  for (auto &itr : equation().probleme().postraitements())
    if (!contient_distance_paroi)
      if (sub_type(Postraitement, itr.valeur()))
        {
          Postraitement& post = ref_cast(Postraitement, itr.valeur());
          for (int i = 0; i < post.noms_champs_a_post().size(); i++)
            {
              // The new global wall distance computation does not use the old method.
              const Nom& chmp = post.noms_champs_a_post()[i];
              if (!chmp.contient("DISTANCE_PAROI_GLOBALE") && chmp.contient("DISTANCE_PAROI"))
                {
                  lire_distance_paroi();
                  contient_distance_paroi = true;
                  break;
                }
            }
        }
  if (nom_fichier_!="")
    loipar_->imprimer_premiere_ligne_ustar(boundaries_, boundaries_list_, nom_fichier_);
  return res;
}
 
bool Modele_turbulence_hyd_base::initTimeStep(double dt)
{
  return true;
}
 
void Modele_turbulence_hyd_base::creer_champ(const Motcle& motlu)
{
  if (loipar_)
    {
      loipar_->creer_champ(motlu);
    }
}
 
bool Modele_turbulence_hyd_base::has_champ(const Motcle& nom, OBS_PTR(Champ_base)& ref_champ) const
{
  if (champs_compris_.has_champ(nom, ref_champ))
    return true;
 
  if (loipar_)
    if (loipar_->has_champ(nom, ref_champ))
      return true;
 
  return false; /* nothing found */
}
 
bool Modele_turbulence_hyd_base::has_champ(const Motcle& nom) const
{
  if (champs_compris_.has_champ(nom))
    return true;
 
  if (loipar_)
    if (loipar_->has_champ(nom))
      return true;
 
  return false; /* nothing found */
}
 
const Champ_base& Modele_turbulence_hyd_base::get_champ(const Motcle& nom) const
{
  OBS_PTR(Champ_base) ref_champ;
 
  if (champs_compris_.has_champ(nom, ref_champ))
    return ref_champ;
 
  if (loipar_)
    if (loipar_->has_champ(nom, ref_champ))
      return ref_champ;
 
  throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
}
 
void Modele_turbulence_hyd_base::get_noms_champs_postraitables(Noms& nom, Option opt) const
{
  if (opt == DESCRIPTION)
    Cerr << que_suis_je() << " : " << champs_compris_.liste_noms_compris() << finl;
  else
    nom.add(champs_compris_.liste_noms_compris());
 
  if (loipar_)
    loipar_->get_noms_champs_postraitables(nom, opt);
}
 
/*! @brief Performs printing if necessary.
 *
 * @param os Output stream.
 */
void Modele_turbulence_hyd_base::imprimer(Sortie& os) const
{
  const Schema_Temps_base& sch = mon_equation_->schema_temps();
  double temps_courant = sch.temps_courant();
  double dt = sch.pas_de_temps();
  if (loipar_ && limpr_ustar(temps_courant, sch.temps_precedent(), dt, dt_impr_ustar_))
    loipar_->imprimer_ustar(os);
  if (loipar_ && limpr_ustar(temps_courant, sch.temps_precedent(), dt, dt_impr_ustar_mean_only_))
    loipar_->imprimer_ustar_mean_only(os, boundaries_, boundaries_list_, nom_fichier_);
}
 
int Modele_turbulence_hyd_base::limpr_ustar(double temps_courant, double temps_prec, double dt, double dt_ustar) const
{
  const Schema_Temps_base& sch = mon_equation_->schema_temps();
  if (sch.nb_pas_dt() == 0)
    return 0;
  if (dt_ustar <= dt
      || ((sch.temps_cpu_max_atteint() || (!get_disable_stop() && sch.stop_lu()) || sch.temps_final_atteint() || sch.nb_pas_dt_max_atteint() || sch.nb_pas_dt() == 1 || sch.stationnaire_atteint())
          && !est_egal(dt_ustar, 1.e20)))
    return 1;
  else
    {
      // See Schema_Temps_base::limpr for information on epsilon and modf
      double i, j, epsilon = 1.e-8;
      modf(temps_courant / dt_ustar + epsilon, &i);
      modf(temps_prec / dt_ustar + epsilon, &j);
      return (i > j);
    }
}
 
 
/*! @brief for PDI IO: retrieve name, type and dimensions of the fields to save/restore
 *
 */
std::vector<YAML_data> Modele_turbulence_hyd_base::data_a_sauvegarder() const
{
  std::vector<YAML_data> data;
  if (loipar_)
    {
      data = loipar_->data_a_sauvegarder();
    }
  return data;
}
 
/*! @brief Saves the turbulence model to an output stream.
 *
 * Saves the object type.
 *
 * @param (Sortie& os) an output stream
 * @return (int) always returns 1
 */
int Modele_turbulence_hyd_base::sauvegarder(Sortie& os) const
{
  a_faire(os);
  if (loipar_)
    {
      loipar_->sauvegarder(os);
    }
  else if(!TRUST_2_PDI::is_PDI_checkpoint())  // OC: for low Re, we write negligeable — no issue doing this here
    {
      const Discretisation_base& discr = mon_equation_->discretisation();
      Nom nom_discr = discr.que_suis_je();
      Nom type = "negligeable_";
      type += nom_discr.substr_old(1, 3);
      os << type << finl;
    }
  // Verify that Modele_turbulence_hyd_base::limiter_viscosite_turbulente() was indeed called
  assert(mp_sum(borne_visco_turb_.size()) != 0);
  return 0;
}
/*! @brief Restores the wall law from a checkpoint.
 *
 * @param (Entree&) an input stream
 * @return (int) always returns 0
 */
int Modele_turbulence_hyd_base::reprendre(Entree& is)
{
  if (loipar_)
    loipar_->reprendre(is);
  return 0;
}
 
/*! @brief Writes an identity record to the output stream if required.
 *
 */
void Modele_turbulence_hyd_base::a_faire(Sortie& os) const
{
  int afaire, special;
  EcritureLectureSpecial::is_ecriture_special(special, afaire);
 
  if (afaire)
    {
      Nom mon_ident(que_suis_je());
      mon_ident += equation().probleme().domaine().le_nom();
      double temps = equation().inconnue().temps();
      mon_ident += Nom(temps, "%e");
      os << mon_ident << finl;
      os << que_suis_je() << finl;
      os.flush();
    }
}
 
void Modele_turbulence_hyd_base::limiter_viscosite_turbulente()
{
  // Initialise
  int size = viscosite_turbulente().valeurs().size();
  if (borne_visco_turb_.size() == 0)
    borne_visco_turb_ = la_viscosite_turbulente_->valeurs();
  borne_visco_turb_ = XNUTM_;
  if (calcul_borne_locale_visco_turb_ && (equation().schema_temps().nb_pas_dt() != 0 || equation().probleme().reprise_effectuee()))
    {
      // Recompute the turbulent viscosity bounds after the first time step
      const Operateur_Diff_base& op_diff = ref_cast(Operateur_Diff_base, equation().operateur(0).l_op_base());
      const Operateur_Conv_base& op_conv = ref_cast(Operateur_Conv_base, equation().operateur(1).l_op_base());
      op_diff.calculer_borne_locale(borne_visco_turb_, op_conv.dt_stab_conv(), dt_diff_sur_dt_conv_);
    }
  // Apply the turbulent viscosity bound
  int nb_elem = equation().domaine_dis().domaine().nb_elem();
  assert(nb_elem == size);
  int compt = 0;
  Debog::verifier("Modele_turbulence_hyd_base::limiter_viscosite_turbulente la_viscosite_turbulente before", la_viscosite_turbulente_->valeurs());
 
  CDoubleArrView borne_visco_turb = borne_visco_turb_.view_ro();
  DoubleArrView corr_visco_turb = static_cast<DoubleVect&>(corr_visco_turb_->valeurs()).view_wo();
  DoubleArrView visco_turb = static_cast<DoubleVect&>(la_viscosite_turbulente_->valeurs()).view_rw();
  // Fill the work arrays:
  Kokkos::parallel_reduce(start_gpu_timer(__KERNEL_NAME__),
                          Kokkos::RangePolicy<>(0, nb_elem), KOKKOS_LAMBDA(
                            const int elem, int& compt_)
  {
    if (visco_turb(elem) > borne_visco_turb(elem))
      {
        compt_++;
        corr_visco_turb(elem) = borne_visco_turb(elem) / visco_turb(elem);
        visco_turb(elem) = borne_visco_turb(elem);
      }
    else
      corr_visco_turb(elem) = 1.;
  }, Kokkos::Sum<int>(compt));
  end_gpu_timer(__KERNEL_NAME__);
  corr_visco_turb_->changer_temps(mon_equation_->inconnue().temps());
  la_viscosite_turbulente_->valeurs().echange_espace_virtuel();
  Debog::verifier("Modele_turbulence_hyd_base::limiter_viscosite_turbulente la_viscosite_turbulente after", la_viscosite_turbulente_->valeurs());
 
  // Print if required
  int imprimer_compt = 0;
  if (mon_equation_->schema_temps().temps_impr() <= mon_equation_->schema_temps().pas_de_temps())
    imprimer_compt = 1;
  else
    {
      if (mon_equation_->schema_temps().nb_pas_dt() == 0)
        imprimer_compt = 0;
      else
        {
          // See Schema_Temps_base::limpr for information on epsilon and modf
          double i, j, epsilon = 1.e-8;
          modf(mon_equation_->schema_temps().temps_courant() / mon_equation_->schema_temps().temps_impr() + epsilon, &i);
          modf((mon_equation_->schema_temps().temps_courant() - mon_equation_->schema_temps().pas_de_temps()) / mon_equation_->schema_temps().temps_impr() + epsilon, &j);
          if (i > j)
            imprimer_compt = 1;
        }
    }
  if (imprimer_compt)
    {
      // PL: optimization to avoid 2 mp_sum instead 1:
      ArrOfTID tmp(2);
      tmp[0] = compt;
      tmp[1] = size;
      mp_sum_for_each_item(tmp);
      /*      compt=mp_sum(compt);
       size=mp_sum(size); */
      double pourcent = 100.0 * (double(tmp[0]) / double(tmp[1]));
      if (je_suis_maitre() && pourcent > 0)
        Cout << "\nTurbulent viscosity has been limited on " << pourcent << " % of cells mesh." << finl;
    }
}