Source_Transport_K_Omega_VDF_Elem.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 <Source_Transport_K_Omega_VDF_Elem.h>
#include <Modele_turbulence_hyd_K_Omega.h>
#include <Pb_Hydraulique_Turbulent.h>
#include <Discretisation_tools.h>
#include <Navier_Stokes_std.h>
#include <Milieu_base.h>
#include <TRUSTTrav.h>
#include <VDF_discretisation.h>
#include <Champ_Fonc_Face_VDF.h>
#include <K_Omega_Eps_constants.h>
#include <Fluide_base.h>
#include <Expert_mode_K_Omega.h>
#include <Dirichlet_paroi_defilante.h>
#include <Dirichlet_paroi_fixe.h>
 
Implemente_instanciable_sans_constructeur(Source_Transport_K_Omega_VDF_Elem,
                                          "Source_Transport_K_Omega_VDF_P0_VDF",
                                          Source_Transport_K_Omega_VDF_Elem_base);
 
Sortie& Source_Transport_K_Omega_VDF_Elem::printOn(Sortie& s) const { return s << que_suis_je() ; }
 
Entree& Source_Transport_K_Omega_VDF_Elem::readOn(Entree& is)
{
  Source_Transport_K_Omega_VDF_Elem_base::verifier_pb_komega(mon_equation->probleme(), que_suis_je());
  Source_Transport_K_Omega_VDF_Elem_base::readOn(is);
  discretiser();
  return is;
}
 
void Source_Transport_K_Omega_VDF_Elem::get_noms_champs_postraitables(Noms& nom,
                                                                      Option opt) const
{
  Source_Transport_K_Omega_VDF_Elem_base::get_noms_champs_postraitables(nom, opt);
  Noms noms_compris = champs_compris_.liste_noms_compris();
  noms_compris.add("grad_k_grad_omega");
 
  if (opt == DESCRIPTION)
    Cerr << que_suis_je() << " : " << noms_compris << finl;
  else
    nom.add(noms_compris);
}
void Source_Transport_K_Omega_VDF_Elem::discretiser()
{
  assert(mon_equation);
  const Discretisation_base& dis = equation().discretisation();
 
  dis.discretiser_champ("champ_elem", equation().domaine_dis(), "grad_k_grad_omega", "", 1, equation().schema_temps().temps_courant(), grad_k_omega_elem_);
  champs_compris_.ajoute_champ(grad_k_omega_elem_);
 
  dis.discretiser_champ("gradient_pression", equation().domaine_dis(), "grad_k", "", Objet_U::dimension, equation().schema_temps().temps_courant(), grad_k_face_);
  champs_compris_.ajoute_champ(grad_k_face_);
 
  dis.discretiser_champ("gradient_pression", equation().domaine_dis(), "grad_omega", "", Objet_U::dimension, equation().schema_temps().temps_courant(), grad_omega_face_);
  champs_compris_.ajoute_champ(grad_omega_face_);
 
 
  dis.discretiser_champ("champ_elem", equation().domaine_dis(), "grad_k_elem", "", Objet_U::dimension, equation().schema_temps().temps_courant(), grad_k_elem_);
  champs_compris_.ajoute_champ(grad_k_elem_);
 
  dis.discretiser_champ("champ_elem", equation().domaine_dis(), "grad_omega_elem", "", Objet_U::dimension, equation().schema_temps().temps_courant(), grad_omega_elem_);
  champs_compris_.ajoute_champ(grad_omega_elem_);
 
  dis.discretiser_champ("champ_elem", equation().domaine_dis(), "production_k", "", 1, equation().schema_temps().temps_courant(), production_k_elem_);
  champs_compris_.ajoute_champ(production_k_elem_);
}
 
void Source_Transport_K_Omega_VDF_Elem::creer_champ(const Motcle& nom)
{
  const Discretisation_base& dis = equation().discretisation();
 
  if ("production_omega" && !production_omega_elem_)
    {
      dis.discretiser_champ("champ_elem", equation().domaine_dis(), "production_omega", "", 1, equation().schema_temps().temps_courant(), production_omega_elem_);
      champs_compris_.ajoute_champ(production_omega_elem_);
    }
  else if ("dissipation_k" && !dissipation_k_elem_)
    {
      dis.discretiser_champ("champ_elem", equation().domaine_dis(), "dissipation_k", "", 1, equation().schema_temps().temps_courant(), dissipation_k_elem_);
      champs_compris_.ajoute_champ(dissipation_k_elem_);
    }
  else if ("dissipation_omega" && !dissipation_omega_elem_)
    {
      dis.discretiser_champ("champ_elem", equation().domaine_dis(), "dissipation_omega", "", 1, equation().schema_temps().temps_courant(), dissipation_omega_elem_);
      champs_compris_.ajoute_champ(dissipation_omega_elem_);
    }
  else if ("cross_diffusion_k_omega" && !cross_diffusion_k_omega_elem_)
    {
      dis.discretiser_champ("champ_elem", equation().domaine_dis(), "cross_diffusion_k_omega", "", 1, equation().schema_temps().temps_courant(), cross_diffusion_k_omega_elem_);
      champs_compris_.ajoute_champ(cross_diffusion_k_omega_elem_);
    }
  else
    Source_Transport_K_Omega_VDF_Elem_base::creer_champ(nom);
}
 
 
void Source_Transport_K_Omega_VDF_Elem::associer_pb(const Probleme_base& pb)
{
  Source_Transport_K_Omega_VDF_Elem_base::associer_pb(pb);
  eqn_K_Omega = ref_cast(Transport_K_Omega, equation());
  turbulence_model = ref_cast(Modele_turbulence_hyd_K_Omega, eqn_K_Omega->modele_turbulence());
}
 
const DoubleTab& Source_Transport_K_Omega_VDF_Elem::get_visc_turb() const
{
  return eqn_K_Omega->modele_turbulence().viscosite_turbulente().valeurs();
}
 
void Source_Transport_K_Omega_VDF_Elem::calculer_terme_production(const Champ_Face_VDF& vitesse, const DoubleTab& visco_turb, const DoubleTab& vit, DoubleVect& P, const bool& deactivate_production_limiter, const double& cst_production_limiter=0.) const
{
  const DoubleTab& K_Omega = eqn_K_Omega->inconnue().valeurs();
  if (axi) calculer_terme_production_K_Axi(le_dom_VDF.valeur(), vitesse, P, K_Omega, visco_turb);
  else calculer_terme_production_K_for_komega(le_dom_VDF.valeur(), le_dom_Cl_VDF.valeur(), P, K_Omega, vit, vitesse, visco_turb, deactivate_production_limiter, cst_production_limiter);
}
 
void Source_Transport_K_Omega_VDF_Elem::fill_resu(const DoubleVect& production_k, DoubleTab& resu) const
{
  const DoubleVect& volumes = le_dom_VDF->volumes();
  const DoubleVect& porosite_vol = le_dom_Cl_VDF->equation().milieu().porosite_elem();
  const DoubleTab& K_Omega = eqn_K_Omega->inconnue().valeurs();
  const bool is_SST_or_BSL = turbulence_model->is_SST_or_BSL();
  const bool is_SST = turbulence_model->is_SST();
  const double LeK_MIN = eqn_K_Omega->modele_turbulence().get_K_MIN();
  const DoubleTab& gradKgradOmega_elem = grad_k_omega_elem_->valeurs();
  DoubleTab production_omega_elem, dissipation_k_elem, dissipation_omega_elem, cross_diffusion_k_omega_elem;
 
  const bool has_post_production_omega = has_champ("production_omega") ;
  const bool has_post_dissipation_k = has_champ("dissipation_k") ;
  const bool has_post_dissipation_omega = has_champ("dissipation_omega") ;
  const bool has_post_cross_diffusion_k_omega = has_champ("cross_diffusion_k_omega") ;
 
  if (has_post_production_omega)
    production_omega_elem = ref_cast_non_const(DoubleTab, production_omega_elem_->valeurs());
 
  if (has_post_dissipation_k)
    dissipation_k_elem = ref_cast_non_const(DoubleTab, dissipation_k_elem_->valeurs());
 
  if (has_post_dissipation_omega)
    dissipation_omega_elem = ref_cast_non_const(DoubleTab, dissipation_omega_elem_->valeurs());
 
  if (has_post_cross_diffusion_k_omega)
    cross_diffusion_k_omega_elem = ref_cast_non_const(DoubleTab, cross_diffusion_k_omega_elem_->valeurs());
 
  const DoubleTab& enstrophy_or_strain_invariant = static_cast<const DoubleTab&>(turbulence_model->get_expert_mode().get_menter_version()==Menter_version::ORIGINAL_1994 ?
                                                                                 turbulence_model->get_enstrophy() : turbulence_model->get_strain_invariant());
  const DoubleTab& F1 = turbulence_model->get_tabF1();
  const DoubleTab& F2 = turbulence_model->get_tabF2();
  const double& gamma1 = turbulence_model->get_expert_mode().get_gamma1();
  const double& gamma2 = turbulence_model->get_expert_mode().get_gamma2();
  for (int elem = 0; elem < le_dom_VDF->nb_elem(); ++elem)
    {
      const double tke = K_Omega(elem, 0);
      const double omega = K_Omega(elem, 1);
      const double volporo = volumes(elem)*porosite_vol(elem);
      const double dissipation_k = BETA_K*tke*omega;
      resu(elem, 0) += (production_k(elem) - dissipation_k)*volporo;
 
      if (tke >= LeK_MIN)
        {
          const double cALPHA = is_SST_or_BSL
                                ? blender(gamma1, gamma2, elem)
                                : ALPHA_OMEGA;
          const double cBETA = is_SST_or_BSL
                               ? blender(BETA1, BETA2, elem)
                               : BETA_OMEGA;
 
          double cSIGMA;
          if (is_SST_or_BSL)
            cSIGMA = 2*(1 - F1(elem))*SIGMA_OMEGA2;
          else
            cSIGMA = (gradKgradOmega_elem(elem) > 0) ? 0.125 : 0.;
 
          double contrib { 0 };
          const double& CST_A1 = turbulence_model->get_CST_A1();
          const double nut = is_SST ? tke * CST_A1 /std::max(CST_A1*omega, enstrophy_or_strain_invariant[elem]*F2[elem]) : tke/omega;
          const double production_omega = cALPHA * production_k(elem) / nut; // production EB: with SST, one cannot use k/omega. See the last paragraph of the appendix of Menter, 1994.
          const double dissipation_omega = cBETA * omega * omega; // dissipation
          const double cross_diffusion_k_omega =  cSIGMA / omega * gradKgradOmega_elem(elem);
 
          contrib += production_omega; //production
          contrib -= dissipation_omega; // dissipation
          contrib += cross_diffusion_k_omega; // cross diffusion
          contrib *= volporo;
          resu(elem, 1) += contrib;
 
          if (has_post_dissipation_k)
            dissipation_k_elem(elem) = dissipation_k;
 
          if (has_post_production_omega)
            production_omega_elem(elem) = production_omega;
 
          if (has_post_dissipation_omega)
            dissipation_omega_elem(elem) = dissipation_omega;
 
          if (has_post_cross_diffusion_k_omega)
            cross_diffusion_k_omega_elem(elem) = cross_diffusion_k_omega;
        }
    }
}
 
void Source_Transport_K_Omega_VDF_Elem::ajouter_blocs(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl) const
{
  Source_Transport_K_Omega_VDF_Elem_base::ajouter_komega(secmem);
 
  const std::string& nom_inco = equation().inconnue().le_nom().getString();
  Matrice_Morse* mat = matrices.count(nom_inco) ? matrices.at(nom_inco) : nullptr;
  if(!mat) return;
 
  const DoubleTab& K_Omega = equation().inconnue().valeurs();
  const DoubleVect& porosite = le_dom_Cl_VDF->equation().milieu().porosite_elem(), &volumes = le_dom_VDF->volumes();
  const int size = K_Omega.dimension(0);
  // on implicite le -eps et le -eps^2/k
  // cAlan : impliciter omega ?
  const DoubleTab& gradKgradOmega_elem = grad_k_omega_elem_->valeurs();
  const DoubleVect& production_TKE = production_k_elem_->valeurs();
  const double& gamma1 = turbulence_model->get_expert_mode().get_gamma1();
  const double& gamma2 = turbulence_model->get_expert_mode().get_gamma2();
  for (int c = 0; c < size; ++c)
    {
      // -eps*vol  donne +vol dans la bonne case
      if (K_Omega(c, 0) > DMINFLOAT)
        {
          const double tke = K_Omega(c, 0);
          const double omega = K_Omega(c, 1);
 
          // cAlan : a adapter
          double coef_k = porosite(c)*volumes(c)*BETA_K*omega;
          (*mat)(c*2, c*2) += coef_k;
          const double cALPHA = turbulence_model->is_SST_or_BSL()
                                ? blender(gamma1, gamma2, c)
                                : ALPHA_OMEGA;
 
          const double cBETA = turbulence_model->is_SST_or_BSL()
                               ? blender(BETA1, BETA2, c)
                               : BETA_OMEGA;
 
          double cSIGMA;
          if (turbulence_model->is_SST_or_BSL())
            cSIGMA = 2*(1 - turbulence_model->get_tabF1()(c))*SIGMA_OMEGA2;
          else
            cSIGMA = (gradKgradOmega_elem(c) > 0) ? 0.125 : 0.;
 
          const double coef_omega = (-cALPHA*production_TKE(c)/tke
                                     + 2.*cBETA*omega
                                     + cSIGMA/(omega*omega)*gradKgradOmega_elem(c) ) * porosite(c)*volumes(c);
          (*mat)(c*2 + 1, c*2 + 1) += coef_omega;
        }
    }
}
 
void Source_Transport_K_Omega_VDF_Elem::compute_cross_diffusion() const
{
  const int nb_elem_tot = le_dom_VDF->nb_elem_tot();
  const int nb_elem = le_dom_VDF->nb_elem();
  DoubleTab& gradKgradOmega_elem = ref_cast_non_const(DoubleTab, grad_k_omega_elem_->valeurs());
  const DoubleTab& K_Omega = eqn_K_Omega->inconnue().valeurs();
  DoubleTab enerK; // field on elem
  DoubleTab omega; // field on elem
 
  enerK.resize(nb_elem_tot,1);
  omega.resize(nb_elem_tot,1);
  for (int elem = 0; elem < nb_elem_tot; ++elem)
    {
      enerK(elem,0) = K_Omega(elem, 0);
      omega(elem,0) = K_Omega(elem, 1);
      gradKgradOmega_elem(elem) = 0;
    }
 
  // Compute the two gradients
  const Operateur_Grad& Op_Grad_komega = eqn_K_Omega->gradient_operator_komega();
  DoubleTab& gradK_face = ref_cast_non_const(DoubleTab, grad_k_face_->valeurs());
  DoubleTab& gradOmega_face = ref_cast_non_const(DoubleTab, grad_omega_face_->valeurs());
  Op_Grad_komega.calculer(enerK, gradK_face);
  Op_Grad_komega.calculer(omega, gradOmega_face);
 
  // Interpolate from faces to elem
  const Domaine_dis_base& domaine_dis = mon_equation->inconnue().domaine_dis_base();
  DoubleTab& tab_gradK_elem = ref_cast_non_const(DoubleTab, grad_k_elem_->valeurs());
  DoubleTab& tab_gradOmega_elem = ref_cast_non_const(DoubleTab, grad_omega_elem_->valeurs());
 
  grad_k_face_->valeur_aux_centres_de_gravite(domaine_dis.domaine(), tab_gradK_elem);
  grad_omega_face_->valeur_aux_centres_de_gravite(domaine_dis.domaine(), tab_gradOmega_elem);
 
  // Dot Product gradKgradOmega
  // EB :  no need to update virtual cells
  for (int elem = 0; elem < nb_elem; ++elem)
    for (int dim=0; dim<dimension; ++dim)
      gradKgradOmega_elem(elem) += tab_gradK_elem(elem,dim) * tab_gradOmega_elem(elem,dim);
}
 
void Source_Transport_K_Omega_VDF_Elem::compute_blending_F1() const
{
  Transport_K_Omega& eq_transport = ref_cast_non_const(Transport_K_Omega, eqn_K_Omega.valeur());
  eq_transport.controler_K_Omega(); // REQUIRED for Runge_Kutta schemes
  const DoubleTab& K_Omega = eqn_K_Omega->inconnue().valeurs();
  const DoubleTab& distmin = le_dom_VDF->y_elem(); // Minimum distance to the edge
  DoubleTab& gradKgradOmega_elem = ref_cast_non_const(DoubleTab, grad_k_omega_elem_->valeurs());
  const Navier_Stokes_std& eq_ns = ref_cast(Navier_Stokes_std, eqn_K_Omega->modele_turbulence().equation());
  const double kinematic_viscosity = eq_ns.fluide().viscosite_cinematique().valeurs()(0,0);
 
  DoubleTab& tabF1 = ref_cast_non_const(DoubleTab, turbulence_model->get_tabF1());
  DoubleTab& tabF2 = ref_cast_non_const(DoubleTab, turbulence_model->get_tabF2());
  const double& cst_min_cd_komega = turbulence_model->get_expert_mode().get_cst_cd_komega();
 
  for (int elem = 0; elem < le_dom_VDF->nb_elem(); ++elem)
    {
      double const dmin = std::max(distmin(elem), 1e-12);
      double const enerK = K_Omega(elem, 0);
      double const omega = K_Omega(elem, 1);
      double const tmp1 = sqrt(enerK)/(BETA_K*omega*dmin);
      double const tmp2 = 500.0*kinematic_viscosity/(omega*dmin*dmin);
      double const maxval = std::max(2.*SIGMA_OMEGA2*gradKgradOmega_elem(elem)/omega,cst_min_cd_komega);
      double const tmp3 = 4.0*SIGMA_OMEGA2*enerK/(maxval*dmin*dmin);
 
      double const arg1 = std::min(std::max(tmp1, tmp2), tmp3); // Common name of the variable
      tabF1(elem) = std::tanh(arg1*arg1*arg1*arg1);
      double const arg2 = std::max(2.*tmp1, tmp2);
      tabF2(elem) = std::tanh(arg2*arg2);
    }
 
  // Apply the boundary conditions for F1 and F2
  for (int i = 0; i < le_dom_Cl_VDF->nb_cond_lim(); i++)
    {
      const Cond_lim_base& boundary_condition = le_dom_Cl_VDF->les_conditions_limites(i).valeur();
      const Front_VF& the_edge = ref_cast(Front_VF, boundary_condition.frontiere_dis());
 
      if (sub_type(Dirichlet_paroi_fixe, boundary_condition) ||
          sub_type(Dirichlet_paroi_defilante, boundary_condition))
        {
          int start = the_edge.num_premiere_face();
          int end = start + the_edge.nb_faces();
          for (int num_elem=start; num_elem<end; num_elem++)
            {
              tabF1(num_elem) = 1.0;
              tabF1(num_elem) = 1.0;
            }
        }
    }
  tabF1.echange_espace_virtuel();
  tabF2.echange_espace_virtuel();
}
 
double Source_Transport_K_Omega_VDF_Elem::blender(double const val1, double const val2,
                                                  int const elem) const
{
  const DoubleTab& F1 = turbulence_model->get_tabF1();
  return F1(elem)*val1 + (1 - F1(elem))*val2;
}