Modele_turbulence_scal_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_scal_base.h>
#include <Modele_turbulence_hyd_base.h>
#include <EcritureLectureSpecial.h>
#include <Convection_Diffusion_std.h>
#include <Discretisation_base.h>
#include <Probleme_base.h>
#include <Domaine.h>
#include <Param.h>
 
Implemente_base(Modele_turbulence_scal_base, "Modele_turbulence_scal_base", Objet_U);
// XD turbulence_paroi_scalaire_base objet_u turbulence_paroi_scalaire_base INHERITS_BRACE Basic class for wall laws for
// XD_CONT energy equation.
 
// XD modele_turbulence_scal_base objet_u modele_turbulence_scal_base INHERITS_BRACE Basic class for turbulence model
// XD_CONT for energy equation.
 
Sortie& Modele_turbulence_scal_base::printOn(Sortie& s) const
{
  return s << que_suis_je() << " " << le_nom();
}
 
Entree& Modele_turbulence_scal_base::readOn(Entree& is)
{
  Cerr << "Reading of data for a " << que_suis_je() << " scalar turbulence model" << finl;
  Param param(que_suis_je());
  Motcle mot = "turbulence_paroi";
  set_param(param);
  param.lire_avec_accolades_depuis(is);
  const Probleme_base& le_pb = mon_equation_->probleme();
  // lp loi de paroi du modele de turbulence de l'hydraulique
  const RefObjU& modele_turbulence = le_pb.equation(0).get_modele(TURBULENCE);
  const Modele_turbulence_hyd_base& mod_turb_hydr = ref_cast(Modele_turbulence_hyd_base, modele_turbulence.valeur());
  if (!loipar_ && mod_turb_hydr.has_loi_paroi_hyd())
    {
      Cerr << finl;
      Cerr << "Mot cle " << mot << "non trouve !" << finl;
      Cerr << "On doit lire une loi de paroi dans le modele de turbulence pour un scalaire." << finl;
      Cerr << finl;
      exit();
    }
  return is;
}
 
void Modele_turbulence_scal_base::set_param(Param& param) const
{
  param.ajouter("dt_impr_nusselt", &dt_impr_nusselt_); // XD_ADD_P floattant
  // XD_CONT Keyword to print local values of Nusselt number and temperature near a wall during a turbulent calculation.
  // XD_CONT The values will be printed in the _Nusselt.face file each dt_impr_nusselt time period. The local Nusselt
  // XD_CONT expression is as follows : Nu = ((lambda+lambda_t)/lambda)*d_wall/d_eq where d_wall is the distance from
  // XD_CONT the first mesh to the wall and d_eq is given by the wall law. This option also gives the value of d_eq and
  // XD_CONT h = (lambda+lambda_t)/d_eq and the fluid temperature of the first mesh near the wall. NL2 For the Neumann
  // XD_CONT boundary conditions (flux_impose), the <<equivalent>> wall temperature given by the wall law is also
  // XD_CONT printed (Tparoi equiv.) preceded for VEF calculation by the edge temperature <<T face de bord>>.
  param.ajouter_non_std("dt_impr_nusselt_mean_only", (this)); // XD attr dt_impr_nusselt_mean_only dt_impr_nusselt_mean_only dt_impr_nusselt_mean_only OPT This keyword is used to print the mean values of Nusselt ( obtained with the wall laws) on each boundary, into a file named datafile_ProblemName_nusselt_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, then you have to specify their names.
  param.ajouter_non_std("turbulence_paroi", this); // XD_ADD_P turbulence_paroi_scalaire_base
  // XD_CONT Keyword to set the wall law.
}
int Modele_turbulence_scal_base::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  Motcle motlu;
  int retval = 1;
  if (mot == "turbulence_paroi")
    {
      Turbulence_paroi_scal_base::typer_lire_turbulence_paroi_scal(loipar_, *this, is);
      is >> loipar_.valeur();
    }
  else if (mot == "dt_impr_nusselt_mean_only")
    {
      if (loipar_->que_suis_je().contient("negligeable"))
        Process::exit("Please remove dt_impr_nusselt option if the wall law is of Negligeable type.");
      // XD dt_impr_nusselt_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() + "_nusselt_mean_only";
      Domaine& dom = equation().probleme().domaine();
      LIST(Nom) nlistbord_dom;                      //!< liste stockant tous les noms de frontiere du domaine
      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_nusselt_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_nusselt_mean_only_;
 
      is >> motlu; // boundaries ou accolade_fermee ou pasbon
      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));
                      //  verif nom bords
                      if (!nlistbord_dom.contient(boundaries_list_[i]))
                        {
                          Cerr << "Problem in the dt_impr_nusselt_mean_only instruction:" << finl;
                          Cerr << "The boundary named '" << boundaries_list_[i] << "' is not a boundary of the domain " << dom.le_nom() << "." << finl;
                          exit();
                        }
                    }
                }
              // lecture accolade fermee
              is >> motlu;
              if (motlu != accolade_fermee)
                {
                  Cerr << "Problem in the dt_impr_nusselt_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_nusselt_mean_only { dt_impr periode [boundaries nb_boundaries boundary_name1 boundary_name2 ... ] } was expected." << finl;
              exit();
            }
        }
    } // fin dt_impr_nusselt_mean_only
  else
    retval = -1;
  return retval;
}
 
/*! @brief Associe l'equation passe en parametre au modele de turbulence.
 *
 * L'equation est caste en Convection_Diffusion_std.
 *
 * @param (Equation_base& eqn) l'equation a laquelle l'objet s'associe
 */
void Modele_turbulence_scal_base::associer_eqn(const Equation_base& eqn)
{
  mon_equation_ = ref_cast(Convection_Diffusion_std, eqn);
}
 
/*! @brief NE FAIT RIEN
 *
 * @param (Domaine_dis_base&) un domaine discretise
 * @param (Domaine_Cl_dis_base&) un domaine de conditions aux limites discretisees
 */
void Modele_turbulence_scal_base::associer(const Domaine_dis_base&, const Domaine_Cl_dis_base&)
{
  ;
}
 
/*! @brief Discretise le modele de turbulence.
 *
 */
void Modele_turbulence_scal_base::discretiser()
{
  Cerr << "Turbulence scalar model discretization" << finl;
  const Schema_Temps_base& sch = mon_equation_->schema_temps();
  const Discretisation_base& dis = mon_equation_->discretisation();
  const Domaine_dis_base& z = mon_equation_->domaine_dis();
  dis.discretiser_champ("champ_elem", z, "diffusivite_turbulente", "m2/s", 1, sch.temps_courant(), diffusivite_turbulente_);
  dis.discretiser_champ("champ_elem", z, "conductivite_turbulente", "W/m/K", 1, sch.temps_courant(), conductivite_turbulente_);
  champs_compris_.ajoute_champ(diffusivite_turbulente_);
}
 
/*! @brief Complete le modele de turbulence: met a jour des references de l'objet.
 *
 */
void Modele_turbulence_scal_base::completer()
{
}
 
/*! @brief A surcharger dans les classes derivees qui possedent une ou des equations Appeler preparer_pas_de_temps sur ces equations.
 *
 */
bool Modele_turbulence_scal_base::initTimeStep(double dt)
{
  return true;
}
 
/*! @brief NE FAIT RIEN
 *
 * @return (int) renvoie toujours 1
 */
int Modele_turbulence_scal_base::preparer_calcul()
{
  if (loipar_)
    loipar_->init_lois_paroi();
  mettre_a_jour(0.);
  if (nom_fichier_!="")
    loipar_->imprimer_premiere_ligne_nusselt(boundaries_, boundaries_list_, nom_fichier_);
  return 1;
}
 
/*! @brief for PDI IO: retrieve name, type and dimensions of the fields to save/restore
 *
 */
std::vector<YAML_data> Modele_turbulence_scal_base::data_a_sauvegarder() const
{
  std::vector<YAML_data> data;
  if (loipar_)
    data = loipar_->data_a_sauvegarder();
  return data;
}
 
/*! @brief NE FAIT RIEN
 *
 * @param (Sortie&) un flot de sortie
 * @return (int) renvoie toujours 1
 */
int Modele_turbulence_scal_base::sauvegarder(Sortie& os) const
{
  if (loipar_)
    return loipar_->sauvegarder(os);
  else
    return 0;
}
 
/*! @brief Effectue l'ecriture d'une identite si cela est necessaire.
 *
 */
void Modele_turbulence_scal_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();
    }
}
 
/*! @brief NE FAIT RIEN
 *
 * @param (Entree&) un flot d'entree
 * @return (int) renvoie toujours 1
 */
int Modele_turbulence_scal_base::reprendre(Entree& is)
{
 
  if (loipar_)
    loipar_->reprendre(is);
  return 1;
}
 
void Modele_turbulence_scal_base::creer_champ(const Motcle& motlu)
{
  if (loipar_)
    loipar_->creer_champ(motlu);
}
 
bool Modele_turbulence_scal_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; /* rien trouve */
}
 
bool Modele_turbulence_scal_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; /* rien trouve */
}
 
const Champ_base& Modele_turbulence_scal_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_scal_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 Indique s'il faut imprimer ou non le Nusselt local
 *
 * @param (double temps_courant, double dt) un flot d'entree
 * @return (int) renvoie 1 si on imprime, 0 sinon
 */
int Modele_turbulence_scal_base::limpr_nusselt(double temps_courant, double temps_prec, double dt, double dt_nusselt) const
{
  const Schema_Temps_base& sch = mon_equation_->schema_temps();
  if (sch.nb_pas_dt() == 0)
    return 0;
  if (dt_nusselt <= 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())
          && dt_nusselt != DMAXFLOAT))
    return 1;
  else
    {
      // Voir Schema_Temps_base::limpr pour information sur epsilon et modf
      double i, j, epsilon = 1.e-8;
      modf(temps_courant / dt_nusselt + epsilon, &i);
      modf(temps_prec / dt_nusselt + epsilon, &j);
      return (i > j);
    }
}
 
/*! @brief Effectue l'impression si cela est necessaire.
 *
 * @return renvoie toujours 1
 */
void Modele_turbulence_scal_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 (limpr_nusselt(temps_courant, sch.temps_precedent(), dt, dt_impr_nusselt_))
    loipar_->imprimer_nusselt(os);
  if (limpr_nusselt(temps_courant, sch.temps_precedent(), dt, dt_impr_nusselt_mean_only_))
    loipar_->imprimer_nusselt_mean_only(os, boundaries_, boundaries_list_, nom_fichier_);
}