Turbulence_paroi_scal_base.cpp

/****************************************************************************
* Copyright (c) 2025, 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_scal_base.h>
#include <Modele_turbulence_hyd_base.h>
#include <Turbulence_paroi_scal_base.h>
#include <Convection_Diffusion_std.h>
#include <Discretisation_base.h>
#include <Schema_Temps_base.h>
#include <Champ_Uniforme.h>
#include <Probleme_base.h>
#include <EcrFicPartage.h>
#include <Domaine_VF.h>
#include <SFichier.h>
#include <Dirichlet_paroi_defilante.h>
#include <Dirichlet_paroi_fixe.h>
 
Implemente_base_sans_constructeur(Turbulence_paroi_scal_base, "Turbulence_paroi_scal_base", Objet_U);
 
Sortie& Turbulence_paroi_scal_base::printOn(Sortie& s) const
{
  return s << que_suis_je() << " " << le_nom();
}
 
Entree& Turbulence_paroi_scal_base::readOn(Entree& is)
{
  return is;
}
 
void Turbulence_paroi_scal_base::creer_champ(const Motcle& motlu)
{
}
 
/*! @brief Reads the characteristics of the scalar 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"
 *       - "loi_standard_hydr_scalaire"
 *
 * @param turb_par Scalar wall law object to type and initialize.
 * @param mod_turb_scal Reference to the scalar turbulence model.
 * @param s Input stream.
 */
void Turbulence_paroi_scal_base::typer_lire_turbulence_paroi_scal(OWN_PTR(Turbulence_paroi_scal_base)& turb_par, const Modele_turbulence_scal_base& mod_turb_scal, Entree& s)
{
  Cerr << "Reading the wall law type " << finl;
  Motcle typ;
  s >> typ;
 
  const Equation_base& eqn = mod_turb_scal.equation().probleme().equation(0); // hydraulic equation
  const RefObjU& modele_turbulence = eqn.get_modele(TURBULENCE);
  const Modele_turbulence_hyd_base& mod_turb_hydr = ref_cast(Modele_turbulence_hyd_base, modele_turbulence.valeur());
  const Turbulence_paroi_base& loi = mod_turb_hydr.loi_paroi();
 
  if (typ != "negligeable_scalaire")
    if ((loi.que_suis_je() == "negligeable_VDF") || (loi.que_suis_je() == "negligeable_VEF"))
      {
        Cerr << "The wall law of type " << typ << " chosen for the scalar is not compatible with" << finl;
        Cerr << "the wall law of type " << loi.que_suis_je() << " chosen for the hydraulics" << finl;
        Cerr << "Use the type 'negligeable_scalaire' for the scalar or use a non-negligible wall law" << finl;
        Cerr << "for the hydraulics" << finl;
        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;
 
  if (typ == "loi_analytique_scalaire_VDF")
    {
      Cerr << "The scalar wall law of type loi_analytique_scalaire" << finl;
      Cerr << "can only be used with a VEF discretization" << finl;
      exit();
    }
 
  Cerr << "and typing: " << typ << finl;
  turb_par.typer(typ);
  turb_par->associer_modele(mod_turb_scal);
  turb_par->associer(eqn.domaine_dis(), eqn.domaine_Cl_dis());
}
 
const Champ_base& Turbulence_paroi_scal_base::get_champ(const Motcle& nom) const
{
  return champs_compris_.get_champ(nom);
}
 
bool Turbulence_paroi_scal_base::has_champ(const Motcle& nom, OBS_PTR(Champ_base)& ref_champ) const
{
  return champs_compris_.has_champ(nom, ref_champ);
}
bool Turbulence_paroi_scal_base::has_champ(const Motcle& nom) const
{
  return champs_compris_.has_champ(nom);
}
 
void Turbulence_paroi_scal_base::get_noms_champs_postraitables(Noms& nom, Option opt) const
{
  if (opt == DESCRIPTION)
    Cerr << "Turbulence_paroi_scal_base : " << champs_compris_.liste_noms_compris() << finl;
  else
    nom.add(champs_compris_.liste_noms_compris());
 
}
/*! @brief Opens or creates a print file for Face, d_eq, local Nu, h.
 *
 * @param Nusselt Shared output file object.
 * @param extension File extension suffix.
 */
void Turbulence_paroi_scal_base::ouvrir_fichier_partage(EcrFicPartage& Nusselt, const Nom& extension) const
{
  const Probleme_base& pb = mon_modele_turb_scal->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 there is no restart
  if (nb_impr_ == 0 && !pb.reprise_effectuee())
    Nusselt.ouvrir(fichier);
  else // Otherwise open in append mode
    Nusselt.ouvrir(fichier, ios::app);
 
  if (je_suis_maitre())
    {
      EcrFicPartage& fic = Nusselt;
      fic << "Temps : " << sch.temps_courant();
    }
  nb_impr_++;
}
 
/*! @brief Opens or creates a file for printing mean values.
 *
 * @param fichier Shared output file object.
 * @param nom_fichier Base filename.
 * @param extension File extension suffix.
 */
void Turbulence_paroi_scal_base::ouvrir_fichier_partage(EcrFicPartage& fichier, const Nom& nom_fichier, const Nom& extension) const
{
  const Probleme_base& pb = mon_modele_turb_scal->equation().probleme();
  Nom nom_fic = nom_fichier + "." + extension;
 
  // Create the file at the first time step if there is 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 Give a boolean indicating if we need to use equivant distance by default we consider that we use the equivalent distance
 *
 * @return (boolean)
 */
bool Turbulence_paroi_scal_base::use_equivalent_distance() const
{
  return true;
}
 
/**
 * @brief Writes header line for Nusselt number and heat transfer 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 thermal statistics including:
 * - Characteristic distance (m)
 * - Local Nusselt number
 * - Convective heat transfer coefficient (W/m²/K)
 * - Near-wall fluid temperature (K)
 * - Boundary face temperature (K)
 * - Equivalent wall temperature (K)
 *
 * Format: "Time Mean(field1) Mean(field2) ... [boundary1(field1) boundary1(field2) ...]"
 * Writes warning messages if specified boundaries are not valid wall types.
 */
void Turbulence_paroi_scal_base::imprimer_premiere_ligne_nusselt(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 = "";
  Noms fields(nb_fields_);
  fields[0] = "dist. carac. (m)";
  fields[1] = "Nusselt (local)";
  fields[2] = "h (Conv. W/m2/K)";
  fields[3] = "Tf cote paroi (K)";
  fields[4] = "T face de bord (K)";
  fields[5] = "Tparoi equiv.(K)";
  for (int i=0; i<nb_fields_; i++)
    {
      ligne+="# ";
      ligne+="field";
      ligne += (Nom)i;
      ligne += " : ";
      ligne += fields[i];
      ligne += "\n";
    }
  ligne += "# Time  ";
  for (int i=0; i<nb_fields_; i++)
    {
      ligne += " \tMean(";
      //ligne += fields[i];
      ligne += "field";
      ligne += (Nom)i;
      ligne += ")";
    }
  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()))
            {
              for (int i=0; i<nb_fields_; i++)
                {
                  ligne += " \t";
                  ligne += nom_bord;
                  ligne += "(";
                  //ligne += fields[i];
                  ligne += "field";
                  ligne += (Nom)i;
                  ligne += ")";
                }
            }
          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 Nusselt mean values.\n\n";
            }
        }
    }
  if (je_suis_maitre())
    {
      fichier << err;
      fichier << ligne;
      fichier << finl;
    }
  fichier.syncfile();
}
 
/**
 * @brief Prints mean Nusselt number and heat transfer 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 values for multiple thermal parameters:
 * - Characteristic distance
 * - Nusselt number
 * - Heat transfer coefficient
 * - Near-wall temperatures
 * - Wall temperatures
 *
 * 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_values] [boundary1_values] [boundary2_values] ...
 */
void Turbulence_paroi_scal_base::imprimer_nusselt_mean_only(Sortie& os, int boundaries_, const LIST(Nom) &boundaries_list, const Nom& nom_fichier_) const
{
  compute_nusselt();
  const Probleme_base& pb = mon_modele_turb_scal->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, nb_fields_+1);
  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);
                  for (int i=0; i<nb_fields_; i++)
                    {
                      moy_bords(0, i) += aire * tab_(num_face, i);
                      moy_bords(num_bord + 1, i) += aire * tab_(num_face, i);
                    }
                  moy_bords(0, nb_fields_) += aire;
                  moy_bords(num_bord + 1, nb_fields_) += aire;
                }
              num_bord += 1;
            }
        }
    }
  mp_sum_for_each_item(moy_bords);
 
// print lines to the file
  if (je_suis_maitre())
    {
      fichier << sch.temps_courant();
      if (moy_bords(0, nb_fields_) != 0)
        for (int i=0; i<nb_fields_; i++)
          fichier << " \t" << moy_bords(0, i) / moy_bords(0, nb_fields_);
    }
 
  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())
                {
                  for (int i=0; i<nb_fields_; i++)
                    fichier << " \t" << moy_bords(num_bord + 1, i) / moy_bords(num_bord + 1, nb_fields_);
                }
              num_bord += 1;
            }
        }
    }
 
  if (je_suis_maitre())
    fichier << finl;
  fichier.syncfile();
}