Turbulence_paroi_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 <Turbulence_paroi_base.h>
#include <Discretisation_base.h>
#include <Schema_Temps_base.h>
#include <Probleme_base.h>
#include <Equation_base.h>
#include <EcrFicPartage.h>
#include <SFichier.h>
#include <Front_VF.h>
#include <Param.h>
#include <Noms.h>
#include <Dirichlet_paroi_defilante.h>
#include <Dirichlet_paroi_fixe.h>
 
 
Implemente_base(Turbulence_paroi_base, "Turbulence_paroi_base", Objet_U);
// XD turbulence_paroi_base objet_u turbulence_paroi_base INHERITS_BRACE Basic class for wall laws for Navier-Stokes
// XD_CONT equations.
 
// XD negligeable turbulence_paroi_base negligeable NO_BRACE Keyword to suppress the calculation of a law of the wall
// XD_CONT with a turbulence model. The wall stress is directly calculated with the derivative of the velocity, in the
// XD_CONT direction perpendicular to the wall (tau_tan /rho= nu dU/dy). NL2 Warning: This keyword is not available for
// XD_CONT k-epsilon models. In that case you must choose a wall law.
// XD negligeable_scalaire turbulence_paroi_scalaire_base negligeable_scalaire NO_BRACE Keyword to suppress the
// XD_CONT calculation of a law of the wall with a turbulence model for thermohydraulic problems. The wall stress is
// XD_CONT directly calculated with the derivative of the velocity, in the direction perpendicular to the wall.
 
Sortie& Turbulence_paroi_base::printOn(Sortie& s) const
{
  return s << que_suis_je() << " " << le_nom();
}
 
Entree& Turbulence_paroi_base::readOn(Entree& is)
{
  return is;
}
 
/*! @brief Reads the characteristics of the wall law from an input stream.
 *
 *     Format:
 *       type_de_loi_de_paroi
 *     Possible values for the wall law type:
 *       - "loi_standard_hydr"
 *       - "negligeable"
 *       - "loi_VanDriest"
 *
 * @param loi_par Wall law object to type and initialize.
 * @param mod_turb_hyd Reference to the hydraulic turbulence model.
 * @param s Input stream.
 */
void Turbulence_paroi_base::typer_lire_turbulence_paroi(OWN_PTR(Turbulence_paroi_base) &loi_par, const Modele_turbulence_hyd_base& mod_turb_hyd, Entree& s)
{
  Cerr << "Reading the wall law type " << finl;
  Motcle typ;
  s >> typ;
 
  const Equation_base& eqn = mod_turb_hyd.equation();
  if (typ == "loi_standard_hydr_scalaire" || typ == "loi_paroi_2_couches_scalaire" || typ == "negligeable_scal")
    {
      Cerr << "The data file format has changed:" << finl;
      Cerr << "Each turbulence model must have its wall law specified." << finl;
      Cerr << "For example, loi_standard_hydr is for the turbulence model of the momentum equation" << finl;
      Cerr << "and loi_standard_hydr_scalaire for the turbulence model of the energy equation." << typ << finl;
      Process::exit();
    }
  typ += "_";
 
  Nom discr = eqn.discretisation().que_suis_je();
 
  //  diffusion operators are shared by VEF and VEFP1B discretisations
  if (discr == "VEFPreP1B") discr = "VEF";
  typ += discr;
 
  Cerr << "and typing: " << typ << finl;
  loi_par.typer(typ);
  loi_par->associer_modele(mod_turb_hyd);
  loi_par->associer(eqn.domaine_dis(), eqn.domaine_Cl_dis());
}
 
void Turbulence_paroi_base::creer_champ(const Motcle& motlu)
{
  if (motlu == "u_star" && !champ_u_star_)
    {
      int nb_comp = 1;
      Noms noms(1);
      noms[0] = "u_star";
      Noms unites(1);
      unites[0] = "m2/s2";
      double temps = 0.;
      const Equation_base& equation = mon_modele_turb_hyd->equation();
      const Discretisation_base& discr = equation.probleme().discretisation();
      discr.discretiser_champ("champ_face", equation.domaine_dis(), scalaire, noms, unites, nb_comp, temps, champ_u_star_);
      champs_compris_.ajoute_champ(champ_u_star_);
    }
}
 
bool Turbulence_paroi_base::has_champ(const Motcle& nom, OBS_PTR(Champ_base)& ref_champ) const
{
  if (champ_u_star_ && (nom == champ_u_star_->le_nom()))
    {
      ref_champ = Turbulence_paroi_base::get_champ(nom);
      return true;
    }
  else
    return champs_compris_.has_champ(nom, ref_champ);
}
 
bool Turbulence_paroi_base::has_champ(const Motcle& nom) const
{
  if (champ_u_star_ && (nom == champ_u_star_->le_nom()))
    return true;
  else
    return champs_compris_.has_champ(nom);
}
 
const Champ_base& Turbulence_paroi_base::get_champ(const Motcle& nom) const
{
  if (champ_u_star_ && nom == champ_u_star_->le_nom())
    {
      // Initialise the u_star volumetric field to 0
      DoubleTab& valeurs = champ_u_star_->valeurs();
      valeurs = 0;
      const Equation_base& my_eqn = mon_modele_turb_hyd->equation();
      if (tab_u_star_.size_array() > 0)
        {
          // Loop over boundaries to retrieve u_star if tab_u_star is allocated
          int nb_front = my_eqn.domaine_dis().nb_front_Cl();
          for (int n_bord = 0; n_bord < nb_front; n_bord++)
            {
              const Cond_lim& la_cl = my_eqn.domaine_Cl_dis().les_conditions_limites(n_bord);
              const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
              int ndeb = le_bord.num_premiere_face();
              int nfin = ndeb + le_bord.nb_faces();
              for (int num_face = ndeb; num_face < nfin; num_face++)
                valeurs(num_face) = tab_u_star_(num_face);
            }
        }
      valeurs.echange_espace_virtuel();
      // Update the field time:
      champ_u_star_->mettre_a_jour(my_eqn.schema_temps().temps_courant());
      return champs_compris_.get_champ(nom);
    }
  else
    return champs_compris_.get_champ(nom);
}
 
void Turbulence_paroi_base::get_noms_champs_postraitables(Noms& nom, Option opt) const
{
  Noms noms_compris = champs_compris_.liste_noms_compris();
  noms_compris.add("u_star");
  if (opt == DESCRIPTION)
    Cerr << que_suis_je() << " : " << noms_compris << finl;
  else
    nom.add(noms_compris);
}
 
/*! @brief Opens or creates a print file for Face, uplus_, dplus_, tab_u_star, Cisaillement_paroi_.
 *
 */
void Turbulence_paroi_base::ouvrir_fichier_partage(EcrFicPartage& Ustar, const Nom& extension) const
{
 
  const Probleme_base& pb = mon_modele_turb_hyd->equation().probleme();
  const Schema_Temps_base& sch = pb.schema_temps();
 
  Nom fichier = Objet_U::nom_du_cas() + "_" + pb.le_nom() + "_" + extension + ".face";
 
  // Create the file at the first time step if no restart
  if (nb_impr_ == 0 && !pb.reprise_effectuee())
    {
      Ustar.ouvrir(fichier);
    }
  // Otherwise open in append mode
  else
    {
      Ustar.ouvrir(fichier, ios::app);
    }
 
  if (je_suis_maitre())
    {
      EcrFicPartage& fic = Ustar;
      fic << "Temps : " << sch.temps_courant();
    }
 
  nb_impr_++;
}
 
/*! @brief Opens or creates a file for printing mean values of uplus_, dplus_, tab_u_star.
 *
 */
void Turbulence_paroi_base::ouvrir_fichier_partage(EcrFicPartage& fichier, const Nom& nom_fichier, const Nom& extension) const
{
  const Probleme_base& pb = mon_modele_turb_hyd->equation().probleme();
  Nom nom_fic = nom_fichier + "." + extension;
 
  // Create the file at the first time step if no restart
  if (nb_impr0_ == 0 && !pb.reprise_effectuee())
    {
      fichier.ouvrir(nom_fic);
    }
  // Otherwise open in append mode
  else
    {
      fichier.ouvrir(nom_fic, ios::app);
    }
  nb_impr0_++;
}
 
/**
 * @brief Writes header line for u* (friction velocity) statistics file
 *
 * @param boundaries_ Flag to control boundary selection (0: all boundaries, 1: specified boundaries)
 * @param boundaries_list List of boundary names to process
 * @param nom_fichier_ Output filename
 *
 * @details Creates header line with column names for u* and d+ statistics.
 * Format: "Time Mean(u*) Mean(d+) [boundary1(u*) boundary1(d+) ...]"
 * Writes warning messages if specified boundaries are not of correct wall type
 * (Dirichlet_paroi_fixe or Dirichlet_paroi_defilante).
 */
void Turbulence_paroi_base::imprimer_premiere_ligne_ustar(int boundaries_, const LIST(Nom) &boundaries_list, const Nom& nom_fichier_) const
{
  EcrFicPartage fichier;
  ouvrir_fichier_partage(fichier, nom_fichier_, "out");
  Nom ligne, err;
  err = "";
  ligne = "# Time   \tMean(u*) \tMean(d+)";
 
  for (int n_bord = 0; n_bord < le_dom_dis_->nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
      const Nom& nom_bord = la_cl->frontiere_dis().le_nom();
      if (je_suis_maitre() && (boundaries_list.contient(nom_bord) || boundaries_list.size() == 0))
        {
          if ( sub_type(Dirichlet_paroi_fixe,la_cl.valeur()) || sub_type(Dirichlet_paroi_defilante, la_cl.valeur()))
            {
              ligne += " \t";
              ligne += nom_bord;
              ligne += "(u*)";
              ligne += " \t";
              ligne += nom_bord;
              ligne += "(d+)";
            }
          else if (boundaries_list.size() > 0)
            {
              err += "The boundary named '";
              err += nom_bord;
              err += "' is not of type Dirichlet_paroi_fixe or Dirichlet_paroi_defilante.\n";
              err += "So TRUST will not write his u_star and d_plus means.\n\n";
            }
        }
    }
  if (je_suis_maitre())
    {
      fichier << err;
      fichier << ligne;
      fichier << finl;
    }
  fichier.syncfile();
}
 
/**
 * @brief Prints mean friction velocity (u*) statistics to a file for specified boundaries
 *
 * @param os Output stream for messages
 * @param boundaries_ Flag to control boundary selection (0: all boundaries, 1: specified boundaries)
 * @param boundaries_list List of boundary names to process when boundaries_=1
 * @param nom_fichier_ Output filename
 *
 * @details Calculates and writes average u* (friction velocity) and d+ values for wall boundaries.
 * Results are written for all boundaries combined and then for each boundary separately.
 * Only processes fixed walls (Dirichlet_paroi_fixe) and moving walls (Dirichlet_paroi_defilante).
 * Output format: time mean_u* mean_d+ [boundary1_u* boundary1_d+ boundary2_u* boundary2_d+ ...]
 */
void Turbulence_paroi_base::imprimer_ustar_mean_only(Sortie& os, int boundaries_, const LIST(Nom) &boundaries_list, const Nom& nom_fichier_) const
{
  const Probleme_base& pb = mon_modele_turb_hyd->equation().probleme();
  const Schema_Temps_base& sch = pb.schema_temps();
  int ndeb, nfin, size0, num_bord;
  num_bord = 0;
 
  if (boundaries_list.size() != 0)
    {
      size0 = boundaries_list.size();
    }
  else
    {
      size0 = le_dom_dis_->nb_front_Cl();
    }
  DoubleTrav moy_bords(size0 + 1, 3);
  moy_bords = 0.;
 
  EcrFicPartage fichier;
  ouvrir_fichier_partage(fichier, nom_fichier_, "out");
 
  for (int n_bord = 0; n_bord < le_dom_dis_->nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
      if ((sub_type(Dirichlet_paroi_fixe, la_cl.valeur())) || (sub_type(Dirichlet_paroi_defilante, la_cl.valeur())))
        {
          const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
          ndeb = le_bord.num_premiere_face();
          nfin = ndeb + le_bord.nb_faces();
          if (boundaries_ == 0 || (boundaries_ == 1 && boundaries_list.contient(le_bord.le_nom())))
            {
              for (int num_face = ndeb; num_face < nfin; num_face++)
                {
                  // surface-averaged values:
                  double aire = le_dom_dis_->face_surfaces(num_face);
                  moy_bords(0, 0) += aire * tab_u_star(num_face);
                  moy_bords(0, 1) += aire;
                  moy_bords(0, 2) += aire * tab_d_plus(num_face);
                  moy_bords(num_bord + 1, 0) += aire * tab_u_star(num_face);
                  moy_bords(num_bord + 1, 1) += aire;
                  moy_bords(num_bord + 1, 2) += aire * tab_d_plus(num_face);
                }
              num_bord += 1;
            }
        }
    }
  mp_sum_for_each_item(moy_bords);
 
// print lines to the file
  if (je_suis_maitre() && moy_bords(0, 1) != 0)
    {
      fichier << sch.temps_courant() << " \t" << moy_bords(0, 0) / moy_bords(0, 1) << " \t" << moy_bords(0, 2) / moy_bords(0, 1);
    }
 
  num_bord = 0;
  for (int n_bord = 0; n_bord < le_dom_dis_->nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = le_dom_Cl_dis_->les_conditions_limites(n_bord);
      if ((sub_type(Dirichlet_paroi_fixe, la_cl.valeur())) || (sub_type(Dirichlet_paroi_defilante, la_cl.valeur())))
        {
          const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
          if (boundaries_ == 0 || (boundaries_ == 1 && boundaries_list.contient(le_bord.le_nom())))
            {
              if (je_suis_maitre())
                {
                  fichier << " \t" << moy_bords(num_bord + 1, 0) / moy_bords(num_bord + 1, 1) << " \t" << moy_bords(num_bord + 1, 2) / moy_bords(num_bord + 1, 1);
                }
              num_bord += 1;
            }
        }
    }
 
  if (je_suis_maitre())
    fichier << finl;
  fichier.syncfile();
}