Modele_turbulence_hyd_K_Omega.cpp

/****************************************************************************
* Copyright (c) 2023, 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 <Modifier_pour_fluide_dilatable.h>
#include <Modele_turbulence_hyd_K_Omega.h>
#include <Modele_turbulence_scal_base.h>
#include <Navier_Stokes_Turbulent.h>
#include <Schema_Temps_base.h>
#include <Champ_Inc_P0_base.h>
#include <communications.h>
#include <Champ_Uniforme.h>
#include <TRUSTTab_parts.h>
#include <Probleme_base.h>
#include <Perf_counters.h>
#include <Fluide_base.h>
#include <TRUSTTrav.h>
#include <Param.h>
#include <Debog.h>
#include <Domaine_VDF.h>
#include <Domaine_VEF.h>
#include <VDF_discretisation.h>
#include <VEF_discretisation.h>
#include <Expert_mode_K_Omega.h>
 
Implemente_instanciable(Modele_turbulence_hyd_K_Omega, "Modele_turbulence_hyd_K_Omega", Modele_turbulence_hyd_RANS_K_Omega_base);
// XD k_omega mod_turb_hyd_rans k_omega INHERITS_BRACE Turbulence model (k-omega).
 
Sortie& Modele_turbulence_hyd_K_Omega::printOn(Sortie& s) const
{
  return s << que_suis_je() << " " << le_nom();
}
 
Entree& Modele_turbulence_hyd_K_Omega::readOn(Entree& s)
{
  Modele_turbulence_hyd_RANS_K_Omega_base::readOn(s);
  is_SST_ = false;
  is_BSL_ = false;
  Navier_Stokes_Turbulent& ns_equation = ref_cast(Navier_Stokes_Turbulent, equation());
  if (model_variant_ == "SST")
    {
      Cout << "SST model: initialize wall distances." << "\n";
      ns_equation.creer_champ("distance_paroi_globale");
      is_SST_ = true;
      is_BSL_ = true;
    }
  else if (model_variant_ == "BSL")
    {
      Cout << "BSL model: initialize wall distances." << "\n";
      ns_equation.creer_champ("distance_paroi_globale");
      is_BSL_ = true;
    }
  else if (model_variant_ != "STD" && model_variant_ != "")
    {
      Cerr << "Error in Modele_turbulence_hyd_K_Omega::readOn()" << finl;
      Cerr << "Unknown model_variant_ : " << model_variant_ << finl;
      Cerr << "Valid options are: STD (default), SST (Shear Stress Tensor)  or BSL (Baseline)." << finl;
      Process::exit();
    }
  return s;
}
 
void Modele_turbulence_hyd_K_Omega::set_param(Param& param) const
{
  Modele_turbulence_hyd_RANS_K_Omega_base::set_param(param);
  param.ajouter("PRANDTL_K|Sigma_K", &Sigma_K_); // XD_ADD_P double
  // XD_CONT Prandtl_K|Sigma_K (default value 1./2.). See https://turbmodels.larc.nasa.gov/wilcox.html
  param.ajouter("PRANDTL_Omega|Sigma_Omega", &Sigma_Omega_); // XD_ADD_P double
  // XD_CONT Prandtl_Omega|Sigma_Omega (default value 1./2.).
  param.ajouter("model_variant", &model_variant_, Param::OPTIONAL); // XD_ADD_P chaine
  // XD_CONT Model variant for k-omega (default value STD)
  param.ajouter("Sigma_K1", &Sigma_K1_); // XD_ADD_P double
  // XD_CONT Sigma_K1 for SST model (default value 0.85). See https://turbmodels.larc.nasa.gov/sst.html
  param.ajouter("Sigma_K2", &Sigma_K2_); // XD_ADD_P double
  // XD_CONT Sigma_K2 for SST model (default value 1.).
  param.ajouter("Sigma_Omega1", &Sigma_OMEGA1_); // XD_ADD_P double
  // XD_CONT Sigma_Omega1 for SST model (default value 1./2.).
  param.ajouter("Sigma_Omega2", &Sigma_OMEGA2_); // XD_ADD_P double
  // XD_CONT Prandt_Omega2 for SST model (default value 0.856).
  param.ajouter("expert_mode", &expert_mode_); // XD_ADD_P bloc_lecture
  // XD_CONT not_set
}
 
 
void Modele_turbulence_hyd_K_Omega::completer()
{
  Modele_turbulence_hyd_RANS_K_Omega_base::completer();
}
 
/*! @brief Calcule la viscosite turbulente au temps demande.
 *
 * @param (double temps) le temps auquel il faut calculer la viscosite
 * @return (Champ_Fonc_base&) la viscosite turbulente au temps demande
 * @throws erreur de taille de visco_turb_K_omega
 */
Champ_Fonc_base& Modele_turbulence_hyd_K_Omega::calculer_viscosite_turbulente(double temps)
{
  const Champ_base& chK_Omega = get_eq_transport().inconnue();
  const Nom type = chK_Omega.que_suis_je();
  const DoubleTab& tab_K_Omega = chK_Omega.valeurs();
  DoubleTab& visco_turb = la_viscosite_turbulente_->valeurs();
 
  // K_Omega(i, 0) = K on node i
  // K_Omega(i, 1) = Omega on node i
 
  int komega_size_real = tab_K_Omega.dimension(0);
  if (komega_size_real < 0)
    {
      // Check field K_Omega type
      if (!sub_type(Champ_Inc_P0_base, chK_Omega))
        {
          Cerr << "Unsupported K_Omega field in "
               "Modele_turbulence_hyd_K_Omega::calculer_viscosite_turbulente"
               << finl;
          Process::exit();
        }
      komega_size_real = get_eq_transport().domaine_dis().domaine().nb_elem();
    }
 
  // cAlan, le 20/01/2023 : sortir cette partie et en faire une fonction à
  // part ? dans le cas d'une domaine nulle on doit effectuer le
  // dimensionnement
  Debog::verifier("Modele_turbulence_hyd_K_Omega::calculer_viscosite_turbulente la_viscosite_turbulente before",
                  la_viscosite_turbulente_->valeurs());
  Debog::verifier("Modele_turbulence_hyd_K_Omega::calculer_viscosite_turbulente tab_K_Omega",
                  tab_K_Omega);
  double non_prepare = (visco_turb.size() == komega_size_real) ? 0 : 1;
  non_prepare = mp_max(non_prepare);
 
  if (non_prepare == 1)
    {
      if (!visco_turb_au_format_K_Omega_)
        {
          // Create field visco_turb_au_format_K_Omega_
          visco_turb_au_format_K_Omega_.typer(type);
          complete_viscosity_field(komega_size_real, get_eq_transport().domaine_dis(),visco_turb_au_format_K_Omega_);
        }
      DoubleTab& visco_turb_K_Omega = visco_turb_au_format_K_Omega_->valeurs();
      if (visco_turb_K_Omega.size() != komega_size_real)
        {
          Cerr << "visco_turb_K_Omega size is " << visco_turb_K_Omega.size()
               << " instead of " << komega_size_real << finl;
          Process::exit();
        }
 
      // A la fin de cette boucle, le tableau visco_turb_K_Omega contient les valeurs de la viscosite turbulente
      // au centre des faces du maillage.
      fill_turbulent_viscosity_tab(tab_K_Omega, visco_turb_K_Omega);
 
      // On connait donc la viscosite turbulente au centre des faces de chaque element
      // On cherche maintenant a interpoler cette viscosite turbulente au centre des elements.
      la_viscosite_turbulente_->affecter(visco_turb_au_format_K_Omega_.valeur());
      Debog::verifier("Modele_turbulence_hyd_K_Omega::calculer_viscosite_turbulente visco_turb_au_format_K_Omega", visco_turb_au_format_K_Omega_.valeur());
    }
  else
    fill_turbulent_viscosity_tab(tab_K_Omega, visco_turb);
 
  la_viscosite_turbulente_->changer_temps(temps);
  Debog::verifier("Modele_turbulence_hyd_K_Omega::calculer_viscosite_turbulente la_viscosite_turbulente after",
                  visco_turb);
  return la_viscosite_turbulente_;
}
 
/*! @brief Fill the turbulent viscosity table depending on the model.
 *
 * If omega is lower that the threshold OMEGA_MIN_, turbulent viscosity is null.
 * If the model variant is SST, a special inverse time scale is computed in place of omega.
 *
 * @param[in] (tab_K_Omega_dim0) the dimension(0) of the K_Omega table, i.e. the number of face or elements.
 * @param[in] (tab_K_Omega&) the table containing the turbulent kinetic and the inverse time scale omega
 * @param[in] (turbulent_viscosity&) the turbulent viscosity table
 */
void Modele_turbulence_hyd_K_Omega::fill_turbulent_viscosity_tab(const DoubleTab& tab_K_Omega,
                                                                 DoubleTab& tab_turbulent_viscosity)
{
  const double OMEGA_MIN = OMEGA_MIN_;
  bool is_sst = is_SST();
 
  CDoubleTabView K_Omega = tab_K_Omega.view_ro();
  DoubleArrView turbulent_viscosity = static_cast<ArrOfDouble&>(tab_turbulent_viscosity).view_wo();
 
  if (is_sst)
    {
 
      CDoubleArrView val = static_cast<const ArrOfDouble&>(get_expert_mode().get_menter_version()==Menter_version::ORIGINAL_1994 ?
                                                           enstrophy_->valeurs() : strain_invariant_ ->valeurs()).view_ro();
      CDoubleArrView F2 = static_cast<const ArrOfDouble&>(tabF2_->valeurs()).view_ro();
 
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), tab_K_Omega.dimension(0), KOKKOS_LAMBDA(const int i)
      {
        const double k = K_Omega(i, 0);
        const double omega = K_Omega(i, 1);
        turbulent_viscosity[i] = omega <= OMEGA_MIN ? 0 : k * CST_A1 / Kokkos::max(CST_A1*omega, val[i]*F2[i]);
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
  else
    {
 
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), tab_K_Omega.dimension(0), KOKKOS_LAMBDA(const int i)
      {
        const double k = K_Omega(i, 0);
        const double omega = K_Omega(i, 1);
        turbulent_viscosity[i] = omega <= OMEGA_MIN ? 0 :  k/omega;
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
 
}
 
 
void Modele_turbulence_hyd_K_Omega::get_noms_champs_postraitables(Noms& nom,
                                                                  Option opt) const
{
  Modele_turbulence_hyd_RANS_K_Omega_base::get_noms_champs_postraitables(nom, opt);
  Noms noms_compris = champs_compris_.liste_noms_compris();
  noms_compris.add("tab_F1");
  noms_compris.add("tab_F2");
  noms_compris.add("enstrophy");
  noms_compris.add("strain_invariant");
 
  if (opt == DESCRIPTION)
    Cerr << " Modele_turbulence_hyd_K_Omega : " << noms_compris << finl;
  else
    nom.add(noms_compris);
}
 
void Modele_turbulence_hyd_K_Omega::discretiser()
{
  Modele_turbulence_hyd_RANS_K_Omega_base::discretiser();
  const Discretisation_base& dis = ref_cast(Discretisation_base, mon_equation_->discretisation());
  Motcle loc="";
  if (sub_type(Domaine_VDF,get_eq_transport().domaine_dis()))
    {
      //const VDF_discretisation& disc = ref_cast(VDF_discretisation, equation().discretisation());
      loc="champ_elem";
    }
  else if (sub_type(Domaine_VEF,get_eq_transport().domaine_dis()))
    {
      loc="champ_face";
    }
  else
    {
      Cerr << "Modele_turbulence_hyd_K_Omega::creer_champ: "
           << "The domain must be VDF or VEF." << finl;
      Process::exit();
    }
 
  Noms noms(1), unites(1);
  if (!tabF1_)
    {
      Nom name="tab_F1";
      dis.discretiser_champ(loc, equation().domaine_dis(), name, "", 1, mon_equation_->schema_temps().temps_courant(), tabF1_);
      champs_compris_.ajoute_champ(tabF1_);
    }
  if (!tabF2_)
    {
      Nom name="tab_F2";
      dis.discretiser_champ(loc, equation().domaine_dis(), name, "", 1, mon_equation_->schema_temps().temps_courant(), tabF2_);
      champs_compris_.ajoute_champ(tabF2_);
    }
  if (!enstrophy_)
    {
      Nom name="enstrophy";
      dis.discretiser_champ(loc, equation().domaine_dis(), name, "", 1, mon_equation_->schema_temps().temps_courant(), enstrophy_);
      champs_compris_.ajoute_champ(enstrophy_);
    }
  if (!strain_invariant_)
    {
      Nom name="strain_invariant";
      dis.discretiser_champ(loc, equation().domaine_dis(), name, "", 1, mon_equation_->schema_temps().temps_courant(), strain_invariant_);
      champs_compris_.ajoute_champ(strain_invariant_);
    }
}
 
 
/*! @brief Prepare the computation of the k-omega model.
 *
 * Initialise tables if model variant is SST and compute the viscosity limit.
 *
 */
int Modele_turbulence_hyd_K_Omega::preparer_calcul()
{
  get_set_eq_transport().preparer_calcul();
  Modele_turbulence_hyd_RANS_K_Omega_base::preparer_calcul();
  calculate_limit_viscosity<MODELE_TYPE::K_OMEGA>(get_set_eq_K_Omega(), -123. /* unused */ );
  Debog::verifier("Modele_turbulence_hyd_K_Omega::preparer_calcul la_viscosite_turbulente", la_viscosite_turbulente_->valeurs());
  return 1;
}
 
bool Modele_turbulence_hyd_K_Omega::initTimeStep(double dt)
{
  return get_set_eq_transport().initTimeStep(dt);
}
 
/*! @brief Performs a time update of the turbulence model.
 *
 * Update the transport equation of k and omega, computes the wall function and the turbulence
 * viscosity.
 *
 * @param[in] (double temps) new time
 */
void Modele_turbulence_hyd_K_Omega::mettre_a_jour(double temps)
{
  Schema_Temps_base& sch = get_set_eq_transport().schema_temps();
  get_set_eq_transport().domaine_Cl_dis().mettre_a_jour(temps);
  if (!get_set_eq_transport().equation_non_resolue())
    {
      statistics().end_count(STD_COUNTERS::update_variables);
      get_set_eq_transport().controler_K_Omega(); // for Euler_Explicit scheme
      sch.faire_un_pas_de_temps_eqn_base(get_set_eq_transport());
      get_set_eq_transport().controler_K_Omega(); // for RUNGE_Kutta schemes
      statistics().begin_count(STD_COUNTERS::update_variables, statistics().get_last_opened_counter_level()+1);
    }
  get_set_eq_transport().mettre_a_jour(temps);
 
  statistics().begin_count(STD_COUNTERS::turbulent_viscosity, statistics().get_last_opened_counter_level()+1);
  Debog::verifier("Modele_turbulence_hyd_K_Omega::mettre_a_jour la_viscosite_turbulente before", la_viscosite_turbulente_->valeurs());
  calculate_limit_viscosity<MODELE_TYPE::K_OMEGA>(get_set_eq_K_Omega(), -123. /* unused */ );
  Debog::verifier("Modele_turbulence_hyd_K_Omega::mettre_a_jour la_viscosite_turbulente after", la_viscosite_turbulente_->valeurs());
  statistics().end_count(STD_COUNTERS::turbulent_viscosity);
}