Navier_Stokes_std.h

/****************************************************************************
* 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.
*
*****************************************************************************/
 
#ifndef Navier_Stokes_std_included
#define Navier_Stokes_std_included
 
#include <Traitement_particulier_NS_base.h>
#include <Assembleur_base.h>
#include <Operateur_Grad.h>
#include <Operateur_Conv.h>
#include <Operateur_Diff.h>
#include <Operateur_Div.h>
#include <TRUST_Deriv.h>
#include <TRUST_Ref.h>
 
class Fluide_base;
 
/*! @brief Navier_Stokes_std This class carries the terms of the momentum equation
 *
 *     for a fluid without turbulence modelling.
 *     The incompressible fluid hypothesis is assumed: div U = 0.
 *     The density is considered constant (equal to rho_0) except in the
 *     gravity force term (Boussinesq hypothesis).
 *     Under these hypotheses, the following form of the Navier-Stokes equations is used:
 *        DU/dt = div(viscous term) - gradP/rho_0 + Bt(T-T0)g + other sources/rho_0
 *        div U = 0
 *     where DU/dt : material derivative of velocity
 *           rho_0 : reference density
 *           T0    : reference temperature
 *           Bt    : thermal expansion coefficient of the fluid
 *           g     : gravity vector
 *     Note: the class implementation allows neglecting
 *           certain terms of the equation (viscous term, convective term,
 *           any given source term).
 *     The unknown is the velocity field.
 *
 *     For handling special cases: addition of Traitement_particulier
 *     example: THI, channel (CA)
 *
 * @sa Equation_base Pb_Hydraulique Pb_Thermohydraulique
 */
class Navier_Stokes_std : public Equation_base
{
  Declare_instanciable_sans_constructeur(Navier_Stokes_std);
 
public :
 
  Navier_Stokes_std();
  void set_param(Param& titi) const override;
  int lire_motcle_non_standard(const Motcle&, Entree&) override;
  void associer_pb_base(const Probleme_base&) override;
  inline void associer_fluide(const Fluide_base& un_fluide ) { le_fluide = un_fluide; }
  const Milieu_base& milieu() const override;
  Milieu_base& milieu() override;
  const Fluide_base& fluide() const;
  Fluide_base& fluide();
  void associer_milieu_base(const Milieu_base& ) override;
  int nombre_d_operateurs() const override;
  int nombre_d_operateurs_tot() const override;
  const Operateur& operateur(int) const override;
  Operateur& operateur(int) override;
  const Operateur& operateur_fonctionnel(int) const override;
  Operateur& operateur_fonctionnel(int) override;
  Operateur_Div& operateur_divergence();
  const Operateur_Div& operateur_divergence() const;
  Operateur_Grad& operateur_gradient();
  const Operateur_Grad& operateur_gradient() const;
  Operateur_Diff& operateur_diff();
  const Operateur_Diff& operateur_diff() const;
  const Champ_Inc_base& inconnue() const override;
  Champ_Inc_base& inconnue() override;
  SolveurSys& solveur_pression();
  void discretiser() override;
  virtual void discretiser_vitesse();
  virtual void discretiser_grad_p();
  void completer() override;
  Entree& lire_cond_init(Entree&) override;
  bool initTimeStep(double dt) override;
  void mettre_a_jour(double temps) override;
  void abortTimeStep() override;
  void resetTime(double time) override;
  int impr(Sortie& os) const override;
 
  void dimensionner_matrice_sans_mem(Matrice_Morse& matrice) override;
 
  /*
    interface {dimensionner,assembler}_blocs
    specificites : prend en compte le gradient de pression (en dernier)
  */
  int has_interface_blocs() const override;
  void dimensionner_blocs(matrices_t matrices, const tabs_t& semi_impl = {}) const override;
  void assembler_blocs(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl = {}) const override;
 
  std::vector<YAML_data> data_a_sauvegarder() const override;
  int sauvegarder(Sortie&) const override;
  int reprendre(Entree&) override;
 
  DoubleTab& derivee_en_temps_inco(DoubleTab& ) override;
  DoubleTab& corriger_derivee_expl(DoubleTab& ) override;
  DoubleTab& corriger_derivee_impl(DoubleTab& ) override;
 
  int preparer_calcul() override;
 
  inline Matrice& matrice_pression() { return matrice_pression_; }
  inline OWN_PTR(Assembleur_base)& assembleur_pression() { return assembleur_pression_; }
 
  inline bool has_grad_P() const { return bool(gradient_P); }
  inline Champ_Inc_base& grad_P() { return gradient_P.valeur(); }
  inline const Champ_Inc_base& grad_P() const { return gradient_P.valeur(); }
  inline Champ_Inc_base& pression() { return la_pression.valeur(); }
  inline const Champ_Inc_base& pression() const { return la_pression.valeur(); }
  inline Champ_Inc_base& pression_pa() { return la_pression_en_pa.valeur(); }
  inline const Champ_Inc_base& pression_pa() const { return la_pression_en_pa.valeur(); }
  inline Champ_Inc_base& div() { return divergence_U.valeur(); }
  inline const Champ_Inc_base& div() const { return divergence_U.valeur(); }
 
  virtual const Champ_Don_base& diffusivite_pour_transport() const;
  virtual const Champ_base& diffusivite_pour_pas_de_temps() const;
  virtual const Champ_base& vitesse_pour_transport() const;
 
  // Methods of the post-processable fields interface
  /////////////////////////////////////////////////////
  void creer_champ(const Motcle& motlu) override;
  const Champ_base& get_champ(const Motcle& nom) const override;
  bool has_champ(const Motcle& nom, OBS_PTR(Champ_base) &ref_champ) const override;
  bool has_champ(const Motcle& nom) const override;
  void get_noms_champs_postraitables(Noms& nom,Option opt=NONE) const override;
  /////////////////////////////////////////////////////
 
  const Motcle& domaine_application() const override;
 
  virtual inline const Champ_Inc_base& vitesse() const { return la_vitesse.valeur(); }
  virtual inline Champ_Inc_base& vitesse() { return la_vitesse.valeur(); }
 
  virtual void projeter();
  virtual int projection_a_faire();
  virtual void sauver() const;
  virtual void calculer_la_pression_en_pa();
  virtual void calculer_pression_hydrostatique(Champ_base& pression_hydro) const;
  int verif_Cl() const override;
 
  virtual const Champ_Inc_base& rho_la_vitesse() const;
  inline Operateur_Conv& get_terme_convectif() { return terme_convectif; }
 
  virtual void updateFluidForce(DoubleTab&);
  // Retrieve true if implicit coupling with another code
  // required to perform filtering of pressure Champ_P1_isoP1Bull during sub-iterations of the implicit loop
  virtual bool getCouplingInfoForFiltering() const ;
  inline DoubleTab getPressureTimeN() {return P_n; } // Pressure at start of time step (used for FSI implicit coupling in sub-iterations)
  void setPressureTimeN();
  void update_y_plus(const DoubleTab& tab) ;
  void reassembler_pression_si_necessaire();
 
protected:
  virtual void discretiser_assembleur_pression();
  virtual void modify_initial_variable() { /* Do nothing */ }
  virtual void modify_initial_gradP( DoubleTrav& ) { /* Do nothing */ }
 
  OBS_PTR(Fluide_base) le_fluide;
 
  OWN_PTR(Champ_Inc_base) la_vitesse, la_pression, divergence_U, gradient_P, la_pression_en_pa;
  OWN_PTR(Champ_Fonc_base)  la_vorticite, grad_u, critere_Q, pression_hydrostatique_, combinaison_champ;
  OWN_PTR(Champ_Fonc_base)  distance_paroi_globale, y_plus, Reynolds_maille, Courant_maille, Taux_cisaillement;
 
  Operateur_Conv terme_convectif;
  Operateur_Diff terme_diffusif;
  Operateur_Div divergence;
  Operateur_Grad gradient;
  Matrice matrice_pression_;
  OWN_PTR(Assembleur_base) assembleur_pression_;
  SolveurSys solveur_pression_;
 
  int projection_initiale;
  double dt_projection, seuil_projection, seuil_uzawa, max_div_U, seuil_divU, raison_seuil_divU;
  mutable double cumulative_;
 
  OWN_PTR(Traitement_particulier_NS_base) le_traitement_particulier;
 
  void uzawa(const DoubleTab&, const Matrice_Base&, SolveurSys&, DoubleTab&, DoubleTab&);
  Nom chaine_champ_combi;
  int methode_calcul_pression_initiale_;
  bool postraiter_gradient_pression_sans_masse_ = false;
  // for genepi it is important to have divu=0 because of error accumulation
  // even if it is not achievable with all schemes
  int div_u_nul_et_non_dsurdt_divu_;
 
private :
  // Pressure at the beginning of the time step, useful for abortTimeStep, in particular for Piso
  // WEC: eventually, use 2 time slots for the la_pression field
  DoubleTab P_n;
  // To avoid computing an extra gradient if no post-processing of the pressure gradient
  mutable int postraitement_gradient_P_;
  double LocalFlowRateRelativeError() const;  // Estimation of a flow rate relative error
};
 
#endif /* Navier_Stokes_std_included */