Champ_Generique_Morceau_Equation.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 <Champ_Generique_Morceau_Equation.h>
#include <Domaine_VF.h>
#include <Champ_Generique_refChamp.h>
#include <Discretisation_base.h>
#include <Operateur_base.h>
#include <Operateur.h>
#include <Synonyme_info.h>
#include <Param.h>
#include <Transport_K_Eps_base.h>
#include <Modele_turbulence_hyd_K_Eps.h>
#include <Modele_turbulence_hyd_K_Eps_Realisable.h>
#include <Modele_turbulence_hyd_RANS_K_Eps_base.h>
#include <Transport_K_Omega_base.h>
#include <Modele_turbulence_hyd_K_Omega.h>
#include <Modele_turbulence_hyd_RANS_K_Omega_base.h>
 
 
Implemente_instanciable(Champ_Generique_Morceau_Equation,"Morceau_Equation",Champ_Gen_de_Champs_Gen);
Add_synonym(Champ_Generique_Morceau_Equation,"Champ_Post_Morceau_Equation");
// XD morceau_equation champ_post_de_champs_post morceau_equation BRACE To calculate a field related to a piece of
// XD_CONT equation. For the moment, the field which can be calculated is the stability time step of an operator
// XD_CONT equation. The problem name and the unknown of the equation should be given by Source refChamp { Pb_Champ
// XD_CONT problem_name unknown_field_of_equation }
 
 
Sortie& Champ_Generique_Morceau_Equation::printOn(Sortie& s ) const
{
  return s << que_suis_je() << " " << le_nom();
}
 
//see Champ_Gen_de_Champs_Gen::readOn
Entree& Champ_Generique_Morceau_Equation::readOn(Entree& s )
{
  Champ_Gen_de_Champs_Gen::readOn(s);
  return s ;
}
 
//  type          : type of equation piece (currently "operateur")
//  numero          : number of the equation piece
//                   (operator case: 0 (diffusion) 1 (convection) 2 (gradient) 3 (divergence))
//  option          : choice of the quantity to post-process
//                     (currently "stabilite" for dt_stab, "flux_bords" for flux_bords_)
//  unite          : to impose the unit of the field
//  compo            : component number to post-process for the "flux_bords" case
//                     (if multiple components)
void Champ_Generique_Morceau_Equation::set_param(Param& param) const
{
  Champ_Gen_de_Champs_Gen::set_param(param);
  param.ajouter("type",&type_morceau_,Param::REQUIRED); // XD_ADD_P chaine
  // XD_CONT can only be operateur for equation operators.
  param.ajouter("numero",&numero_morceau_); // XD_ADD_P entier
  // XD_CONT numero will be 0 (diffusive operator) or 1 (convective operator) or 2 (gradient operator) or 3 (divergence
  // XD_CONT operator).
  param.ajouter("unite", &unite_); // XD_ADD_P chaine
  // XD_CONT will specify the field unit
  param.ajouter("option",&option_,Param::REQUIRED); // XD_ADD_P chaine(into=["stabilite","flux_bords","flux_surfacique_bords"])
  // XD_CONT option is stability for time steps or flux_bords for boundary fluxes or flux_surfacique_bords for boundary
  // XD_CONT surfacic fluxes
  param.ajouter("compo",&compo_); // XD_ADD_P entier
  // XD_CONT compo will specify the number component of the boundary flux (for boundary fluxes, in this case compo
  // XD_CONT permits to specify the number component of the boundary flux choosen).
}
 
//-Initialization of ref_eq_
//-Initialization of localisation_
//-Call to the parent class method
void Champ_Generique_Morceau_Equation::completer(const Postraitement_base& post)
{
  const Probleme_base& Pb = get_ref_pb_base();
  int numero_eq_=-1;
  if (sub_type(Champ_Generique_refChamp,get_source(0)))
    {
 
      OWN_PTR(Champ_base) espace_stockage;
      const Champ_base& mon_champ = get_source(0).get_champ(espace_stockage);
      if (sub_type(Champ_Inc_base,mon_champ))
        {
          const Champ_Inc_base& mon_champ_inc = ref_cast(Champ_Inc_base,mon_champ);
 
          {
            //We retrieve the equation even though it is not yet associated with Champ_Inc
            //We iterate over the equations of the problem and identify the one corresponding to the unknown field
 
            int nb_eq = Pb.nombre_d_equations();
            int i=0;
 
            while (i<nb_eq)
              {
                const Equation_base& eq_test = Pb.equation(i);
                if ((eq_test.inconnue().le_nom() == mon_champ_inc.le_nom()))
                  {
                    numero_eq_=i;
                    ref_eq_=Pb.equation(numero_eq_);
                    break;
                  }
 
                if (mon_champ_inc.le_nom() == "K_Eps" || mon_champ_inc.le_nom() == "K_Omega")
                  {
                    if (not sub_type(Modele_turbulence_hyd_2_eq_base, eq_test.get_modele(TURBULENCE).valeur()))
                      {
                        Cerr<<"Champ_Generique_Morceau_Equation in TriocCFD: K_Eps and K_Omega can only be used with a related RANS turbulence model"<<finl;
                        Process::exit();
                      }
                    // handles RANS fields in triocfd: search if turbulence models have an equation with this unknown name
                    const auto& le_mod_RANS = ref_cast(Modele_turbulence_hyd_2_eq_base, eq_test.get_modele(TURBULENCE).valeur());
                    const auto& rans_equation =  le_mod_RANS.get_eq_transport();
                    if ((rans_equation.inconnue().le_nom() == mon_champ_inc.le_nom()))
                      {
                        numero_eq_=i;
                        ref_eq_=rans_equation;
                        break;
                      }
                    Cerr << "Champ_Generique_Morceau_Equation in TriocCFD: K_Eps or K_Omega was probably used with the wrong RANS turbulence model." << finl;
                    Cerr << "Check your dataset consistency." << finl;
                    Process::exit();
                  }
                i++;
              }
          }
        }
    }
 
  if (numero_eq_==-1)
    {
      Cerr<<"We can apply a Champ_Generique_Morceau_Equation only to an unknown field of the problem"<<finl;
      exit();
    }
  // ref_eq_ = Pb.equation(numero_eq_); // not needed, we already did that when we had the correct equation earlier...
 
  localisation_ = morceau().get_localisation_pour_post(option_);
 
  //The parent class method is called after so that
  //the nommer_source() method has ref_eq_ initialized
  Champ_Gen_de_Champs_Gen::completer(post);
}
 
OWN_PTR(Champ_Fonc_base)& Champ_Generique_Morceau_Equation::creer_espace_stockage(const Nature_du_champ& nature,
                                                                                  const int nb_comp,
                                                                                  OWN_PTR(Champ_Fonc_base)& es_tmp) const
{
  Noms noms;
  Noms unites;
  for (int c=0; c<nb_comp; c++)
    unites.add("??");
  noms.add("bidon");
  double temps;
  temps = get_time();
  const Discretisation_base&  discr = get_discretisation();
  Motcle directive;
  if ((Motcle(option_)=="flux_bords") || (Motcle(option_)=="flux_surfacique_bords"))
    directive = "champ_face"; // To have flux_bords at faces and not at elements in VDF (ex: temperature)
  else
    directive = get_directive_pour_discr();
  const Domaine_dis_base& domaine_dis = get_ref_domaine_dis_base();
 
  discr.discretiser_champ(directive,domaine_dis,nature,noms,unites,nb_comp,temps,es_tmp);
 
  if (directive=="pression")
    {
      const  Domaine_Cl_dis_base& zcl = get_ref_zcl_dis_base();
      es_tmp->completer(zcl);
    }
 
  return es_tmp;
}
 
/*! @brief The equation piece triggers the discretisation of the storage space field and fills its value array via the calculer_pour_post(.
 *
 * ..) method
 *  Note: This procedure differs from the one applied in the other Champ_Generique_base classes, for which
 *  the filling of the espace_stockage value array is not delegated
 *
 *
 */
const Champ_base& Champ_Generique_Morceau_Equation::get_champ(OWN_PTR(Champ_base)& espace_stockage) const
{
  OWN_PTR(Champ_Fonc_base)  es_tmp;
  Nature_du_champ nature;
  int nb_comp = -1;
 
  if ((Motcle(option_)=="stabilite") || (Motcle(option_)=="flux_bords") || (Motcle(option_)=="flux_surfacique_bords"))
    {
      nature = scalaire;
      nb_comp = 1;
    }
  else
    {
      Cerr << "Error. Contact TRUST support." << finl;
      Process::exit();
    }
 
  espace_stockage = creer_espace_stockage(nature,nb_comp,es_tmp);
  morceau().calculer_pour_post(espace_stockage,option_,compo_);
  DoubleTab& es_val = espace_stockage->valeurs();
  es_val.echange_espace_virtuel();
  return espace_stockage;
}
 
const Champ_base& Champ_Generique_Morceau_Equation::get_champ_without_evaluation(OWN_PTR(Champ_base)& espace_stockage) const
{
  OWN_PTR(Champ_Fonc_base)  es_tmp;
  Nature_du_champ nature;
  int nb_comp = -1;
 
  if ((Motcle(option_)=="stabilite") || (Motcle(option_)=="flux_bords") || (Motcle(option_)=="flux_surfacique_bords"))
    {
      nature = scalaire;
      nb_comp = 1;
    }
 
  espace_stockage = creer_espace_stockage(nature,nb_comp,es_tmp);
  return espace_stockage;
}
const Noms Champ_Generique_Morceau_Equation::get_property(const Motcle& query) const
{
 
  Motcles motcles(2);
  motcles[0] = "composantes";
  motcles[1] = "unites";
 
  int rang = motcles.search(query);
  switch(rang)
    {
    case 0:
      {
        Noms compo(1);
        Nom nume((int)1);
        compo[0] = nom_post_+nume;
 
        return compo;
      }
 
    case 1:
      {
        Noms unites(1);
        if (unite_ != "??")
          {
            unites[0] = unite_;
            return unites;
          }
        if (Motcle(option_)=="stabilite")
          unites[0] = "s";
        else if (Motcle(option_).debute_par("FLUX_"))
          {
            // Very incomplete but still...:
            if (ref_eq_->inconnue().le_nom()=="vitesse")
              {
                if (numero_morceau_<2) unites[0]="N";
                else if (numero_morceau_==3) unites[0]=(dimension==2 ? "m2/s" : "m3/s");
              }
            else if (ref_eq_->inconnue().le_nom()=="temperature") unites[0]="W";
            if (Motcle(option_)==Motcle("flux_surfacique_bords"))
              unites[0]+="/m2";
          }
        return unites;
      }
 
    }
  return Champ_Gen_de_Champs_Gen::get_property(query);
}
 
Entity Champ_Generique_Morceau_Equation::get_localisation(const int index) const
{
  Entity loc;
  // For initialization
  loc = Entity::NODE;
 
  if (localisation_=="elem")
    loc = Entity::ELEMENT;
  else if (localisation_=="face")
    loc = Entity::FACE;
  else if (localisation_=="som")
    loc = Entity::NODE;
  else
    {
      Cerr << "Error of type : localisation should be specified to elem or som or face for the defined field " << nom_post_ << finl;
      throw Champ_Generique_erreur("INVALID");
    }
  return loc;
}
 
//Name the field as a source by default
//nom_eq + "_" + type_morceau + "_" + type_option
void Champ_Generique_Morceau_Equation::nommer_source()
{
  Motcle type_morceau, type_option;
  if (Motcle(type_morceau_)=="operateur")
    {
      if (numero_morceau_==0)
        type_morceau = "diffusion";
      else if (numero_morceau_==1)
        type_morceau = "convection";
      else if (numero_morceau_==2)
        type_morceau = "gradient";
      else if (numero_morceau_==3)
        type_morceau = "divergence";
      else
        {
          Cerr<<"This operateur number is not available."<<finl;
          Process::exit();
        }
    }
  else
    {
      Cerr<<"Currently the only pieces of equations considered for the postprocessing"<<finl;
      Cerr<<"are of type operator"<<finl;
      Process::exit();
    }
 
  if (Motcle(option_)=="stabilite")
    type_option = "dt";
  else if (Motcle(option_)=="flux_bords")
    {
      type_option = "flux_bords";
    }
  else if (Motcle(option_)=="flux_surfacique_bords")
    {
      type_option = "flux_surfacique_bords";
    }
  else
    {
      Cerr<<"Currently the only options considered for a piece of equation are \"stabilite\" and \"flux_bords\"."<<finl;
      Process::exit();
    }
 
 
  if (nom_post_=="??")
    {
      Nom nom_post_source;
      const Nom& nom_eq = ref_eq_->le_nom();
      nom_post_source = nom_eq;
      nom_post_source += "_";
      nom_post_source +=  type_morceau;
      nom_post_source += "_";
      nom_post_source += type_option ;
      if ((Motcle(option_)=="flux_bords") || (Motcle(option_)=="flux_surfacique_bords"))
        {
          Nom nume(compo_);
          nom_post_source += nume;
        }
      nommer(nom_post_source);
    }
 
}
 
//Returns the equation piece considered for this Champ_Generique_base
//Currently the only type of piece considered: operators
//To consider other equation pieces, type_morceau_ will need to be tested
const MorEqn& Champ_Generique_Morceau_Equation::morceau() const
{
  int nb_operateurs = ref_eq_->nombre_d_operateurs();
  if (numero_morceau_ >= nb_operateurs)
    return ref_eq_->operateur_fonctionnel(numero_morceau_-nb_operateurs).l_op_base();
 
  return ref_eq_->operateur(numero_morceau_).l_op_base();
}
 
MorEqn& Champ_Generique_Morceau_Equation::morceau()
{
  int nb_operateurs = ref_eq_->nombre_d_operateurs();
  if (numero_morceau_ >= nb_operateurs)
    return ref_eq_->operateur_fonctionnel(numero_morceau_-nb_operateurs).l_op_base();
 
  return ref_eq_->operateur(numero_morceau_).l_op_base();
}
 
const Motcle Champ_Generique_Morceau_Equation::get_directive_pour_discr() const
{
  Motcle directive;
 
  if (localisation_=="elem")
    directive = "champ_elem";
  else if (localisation_=="face")
    directive = "champ_face";
  else if (localisation_=="som")
    directive = "champ_sommets";
  else
    {
      Cerr << "Error of type : localisation should be specified to elem or som or face for the defined field " << nom_post_ << finl;
      throw Champ_Generique_erreur("INVALID");
    }
 
  return directive;
}