Preconditioning for multicomponent flows

.github/workflows/regression.yml

@@ -211,7 +211,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:250717-1402
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t feature_preconditioning -c feature_restart_preconditioning -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:250717-1402
         with:

Common/include/CConfig.hpp

@@ -912,6 +912,7 @@ class CConfig {
   Initial_BCThrust,                /*!< \brief Ratio of turbulent to laminar viscosity at the actuator disk. */
   Pressure_FreeStream,             /*!< \brief Total pressure of the fluid. */
   Pressure_Thermodynamic,          /*!< \brief Thermodynamic pressure of the fluid. */
+  Standard_Ref_Temperature,        /*!< \brief Standard reference temperature for multicomponent flows. */
   Temperature_FreeStream,          /*!< \brief Total temperature of the fluid.  */
   Temperature_ve_FreeStream;       /*!< \brief Total vibrational-electronic temperature of the fluid.  */
   unsigned short wallModel_MaxIter; /*!< \brief maximum number of iterations for the Newton method for the wall model */
@@ -1951,6 +1952,12 @@ class CConfig {
    */
   su2double GetPressure_Thermodynamic(void) const { return Pressure_Thermodynamic; }
 
+  /*!
+   * \brief Get the value of the standard reference temperature for multicomponent flows.
+   * \return Standard reference temperature, Non-dimensionalized if it is needed for Non-Dimensional problems.
+   */
+  su2double GetStandard_RefTemperatureND(void) const { return Standard_Ref_Temperature / Temperature_Ref; }
+
   /*!
    * \brief Get the value of the non-dimensionalized thermodynamic pressure.
    * \return Non-dimensionalized thermodynamic pressure.

Common/src/CConfig.cpp

@@ -1229,6 +1229,9 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("GAMMA_VALUE", Gamma, 1.4);
   /*!\brief THERMODYNAMIC_PRESSURE  \n DESCRIPTION: Thermodynamics(operating) Pressure (101325 Pa), only for incompressible flows) \ingroup Config*/
   addDoubleOption("THERMODYNAMIC_PRESSURE", Pressure_Thermodynamic, 101325.0);
+  /*!\brief STANDARD_REFERENCE_TEMPERATURE  \n DESCRIPTION: Standard reference temperature (298.15K), only for
+   * multicomponent incompressible flows) \ingroup Config*/
+  addDoubleOption("STANDARD_REFERENCE_TEMPERATURE", Standard_Ref_Temperature, 298.15);
   /*!\brief CP_VALUE  \n DESCRIPTION: Specific heat at constant pressure, Cp (1004.703 J/kg*K (air), constant density incompressible fluids only) \ingroup Config*/
   addDoubleListOption("SPECIFIC_HEAT_CP", nSpecific_Heat_Cp, Specific_Heat_Cp);
   /*!\brief THERMAL_EXPANSION_COEFF  \n DESCRIPTION: Thermal expansion coefficient (0.00347 K^-1 (air), used for Boussinesq approximation for liquids/non-ideal gases) \ingroup Config*/

SU2_CFD/include/fluid/CConstantDensity.hpp

@@ -56,5 +56,15 @@ class CConstantDensity final : public CFluidModel {
        decoupled equation. Hence, we update the value.
        Note Cp = Cv, (gamma = 1).*/
     Temperature = t;
+    Enthalpy = Cp * Temperature;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using Enthalpy.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
+    Enthalpy = val_enthalpy;
+    Temperature = Enthalpy / Cp;
   }
 };

SU2_CFD/include/fluid/CFluidFlamelet.hpp

@@ -118,6 +118,13 @@ class CFluidFlamelet final : public CFluidModel {
    */
   void SetTDState_T(su2double val_temperature, const su2double* val_scalars = nullptr) override;
 
+  /*!
+   * \brief Set the thermodynamic state.
+   * \param[in] val_enthalpy - enthalpy
+   * \param[in] val_scalars - pointer to species mass fractions
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
+
   /*!
    * \brief Evaluate data-set for flamelet simulations.
    * \param[in] input_scalar - controlling variables used to interpolate manifold.

SU2_CFD/include/fluid/CFluidModel.hpp

@@ -47,6 +47,7 @@ class CLookUpTable;
 class CFluidModel {
  protected:
   su2double StaticEnergy{0.0};     /*!< \brief Internal Energy. */
+  su2double Enthalpy{0.0};         /*!<*\brief Enthalpy. */
   su2double Entropy{0.0};          /*!< \brief Entropy. */
   su2double Density{0.0};          /*!< \brief Density. */
   su2double Pressure{0.0};         /*!< \brief Pressure. */
@@ -113,6 +114,11 @@ class CFluidModel {
    */
   su2double GetStaticEnergy() const { return StaticEnergy; }
 
+  /*!
+   * \brief Get fluid enthalpy.
+   */
+  su2double GetEnthalpy() const { return Enthalpy; }
+
   /*!
    * \brief Get fluid density.
    */
@@ -186,6 +192,16 @@ class CFluidModel {
     return mass_diffusivity;
   }
 
+  /*!
+   * \brief Get heat diffusivity terms.
+   */
+  virtual void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions = nullptr) {}
+
+  /*!
+   * \brief Get gradient heat diffusivity terms.
+   */
+  virtual void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions = nullptr) {}
+
   /*!
    * \brief Get fluid pressure partial derivative.
    */
@@ -339,6 +355,13 @@ class CFluidModel {
    */
   virtual void SetTDState_T(su2double val_Temperature, const su2double* val_scalars = nullptr) {}
 
+  /*!
+   * \brief Virtual member.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   * \param[in] val_scalars - Scalar mass fractions. 
+   */
+  virtual void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) {}
+
   /*!
    * \brief Set fluid eddy viscosity provided by a turbulence model needed for computing effective thermal conductivity.
    */

SU2_CFD/include/fluid/CFluidScalar.hpp

@@ -42,6 +42,7 @@ class CFluidScalar final : public CFluidModel {
   const int n_species_mixture;            /*!< \brief Number of species in mixture. */
   su2double Gas_Constant;           /*!< \brief Specific gas constant. */
   const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */
+  const su2double Ref_Temperature;        /*!< \brief Reference temperature. */
   const su2double GasConstant_Ref;        /*!< \brief Gas constant reference needed for Nondimensional problems. */
   const su2double Prandtl_Number;         /*!< \brief Prandlt number.*/
 
@@ -91,6 +92,11 @@ class CFluidScalar final : public CFluidModel {
    */
   su2double ComputeMeanSpecificHeatCp(const su2double* val_scalars);
 
+  /*!
+   * \brief Compute Enthalpy given the temperature and scalars.
+   */
+  su2double ComputeEnthalpyFromT(const su2double val_temperature, const su2double* val_scalars);
+
   /*!
    * \brief Compute gas constant for mixture.
    */
@@ -137,9 +143,26 @@ class CFluidScalar final : public CFluidModel {
    */
   inline su2double GetMassDiffusivity(int ivar) override { return massDiffusivity[ivar]; }
 
+  /*!
+   * \brief Get enthalpy diffusivity terms.
+   */
+  void GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) override;
+
+  /*!
+   * \brief Get gradient enthalpy diffusivity terms.
+   */
+  void GetGradEnthalpyDiffusivity(su2double* grad_enthalpy_diffusions) override;
+
   /*!
    * \brief Set the Dimensionless State using Temperature.
    * \param[in] t - Temperature value at the point.
    */
   void SetTDState_T(su2double val_temperature, const su2double* val_scalars) override;
+
+  /*!
+   * \brief Virtual member.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   * \param[in] val_scalars - Scalar mass fractions. 
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
 };

SU2_CFD/include/fluid/CIncIdealGas.hpp

@@ -58,6 +58,17 @@ class CIncIdealGas final : public CFluidModel {
     /*--- The EoS only depends upon temperature. ---*/
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);
+    Enthalpy = Cp * Temperature;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using Enthalpy.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
+    Enthalpy = val_enthalpy;
+    Temperature = Enthalpy / Cp;
+    Density = Pressure / (Temperature * Gas_Constant);
   }
 
  private:

SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp

@@ -80,6 +80,25 @@ class CIncIdealGasPolynomial final : public CFluidModel {
       Cp += coeffs_[i] * t_i;
     }
     Cv = Cp / Gamma;
+    Enthalpy = Cp * Temperature;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using enthalpy.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
+    Enthalpy = val_enthalpy;
+    Temperature = Enthalpy / Cp;
+    Density = Pressure / (Temperature * Gas_Constant);
+    /* Evaluate the new Cp from the coefficients and temperature. */
+    Cp = coeffs_[0];
+    su2double t_i = 1.0;
+    for (int i = 1; i < N; ++i) {
+      t_i *= Temperature;
+      Cp += coeffs_[i] * t_i;
+    }
+    Cv = Cp / Gamma;
   }
 
  private:

SU2_CFD/include/numerics/CNumerics.hpp

@@ -129,6 +129,10 @@ class CNumerics {
   const su2double
   *ScalarVar_i,   /*!< \brief Vector of scalar variables at point i. */
   *ScalarVar_j;   /*!< \brief Vector of scalar variables at point j. */
+  su2double
+  HeatFluxDiffusion;   /*!< \brief Heat flux due to enthalpy diffusion for multicomponent. */
+  su2double
+  JacHeatFluxDiffusion;   /*!< \brief Heat flux jacobian due to enthalpy diffusion for multicomponent. */
   const su2double
   *TransVar_i,  /*!< \brief Vector of turbulent variables at point i. */
   *TransVar_j;  /*!< \brief Vector of turbulent variables at point j. */
@@ -187,6 +191,8 @@ class CNumerics {
 
   bool nemo;                      /*!< \brief Flag for NEMO problems  */
 
+  bool energy_multicomponent = false; /*!< \brief Flag for multicomponent and reacting flow  */
+
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
 
 public:
@@ -750,6 +756,20 @@ class CNumerics {
     Diffusion_Coeff_j = val_diffusioncoeff_j;
   }
 
+  /*!
+   * \brief Set the heat flux due to enthalpy diffusion
+   * \param[in] val_heatfluxdiffusion - Value of the heat flux due to enthalpy diffusion.
+   */
+  inline void SetHeatFluxDiffusion(su2double val_heatfluxdiffusion) { HeatFluxDiffusion = val_heatfluxdiffusion; }
+
+  /*!
+   * \brief Set Jacobian of the heat flux due to enthalpy diffusion
+   * \param[in] val_jacheatfluxdiffusion - Value of the heat flux jacobian due to enthalpy diffusion.
+   */
+  inline void SetJacHeatFluxDiffusion(su2double val_jacheatfluxdiffusion) {
+    JacHeatFluxDiffusion = val_jacheatfluxdiffusion;
+  }
+
   /*!
    * \brief Set the laminar viscosity.
    * \param[in] val_laminar_viscosity_i - Value of the laminar viscosity at point i.
@@ -1038,8 +1058,7 @@ class CNumerics {
    * \param[in] val_density - Value of the density.
    * \param[in] val_velocity - Pointer to the velocity.
    * \param[in] val_betainc2 - Value of the artificial compresibility factor.
-   * \param[in] val_cp - Value of the specific heat at constant pressure.
-   * \param[in] val_temperature - Value of the temperature.
+   * \param[in] val_enthalpy - Value of the enthalpy.
    * \param[in] val_dRhodT - Value of the derivative of density w.r.t. temperature.
    * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face.
    * \param[in] val_scale - Scale of the projection.
@@ -1048,8 +1067,7 @@ class CNumerics {
   void GetInviscidIncProjJac(const su2double *val_density,
                              const su2double *val_velocity,
                              const su2double *val_betainc2,
-                             const su2double *val_cp,
-                             const su2double *val_temperature,
+                             const su2double *val_enthalpy,
                              const su2double *val_dRhodT,
                              const su2double *val_normal,
                              su2double val_scale,
@@ -1060,17 +1078,15 @@ class CNumerics {
    * \param[in] val_density - Value of the density.
    * \param[in] val_velocity - Pointer to the velocity.
    * \param[in] val_betainc2 - Value of the artificial compresibility factor.
-   * \param[in] val_cp - Value of the specific heat at constant pressure.
-   * \param[in] val_temperature - Value of the temperature.
-   * \param[in] val_dRhodT - Value of the derivative of density w.r.t. temperature.
+   * \param[in] val_enthalpy - Value of the enthalpy.
+   * \param[in] val_dRhodh - Value of the derivative of density w.r.t. enthalpy.
    * \param[out] val_Precon - Pointer to the preconditioning matrix.
    */
   void GetPreconditioner(const su2double *val_density,
                          const su2double *val_velocity,
                          const su2double *val_betainc2,
-                         const su2double *val_cp,
-                         const su2double *val_temperature,
-                         const su2double *val_drhodt,
+                         const su2double *val_enthalpy,
+                         const su2double *val_drhodh,
                          su2double **val_Precon) const;
 
   /*!

SU2_CFD/include/numerics/flow/convection/centered.hpp

@@ -49,7 +49,7 @@ class CCentLaxInc_Flow final : public CNumerics {
   MeanDensity, MeanPressure,
   MeanBetaInc2, MeanEnthalpy,
   MeanCp, MeanTemperature,         /*!< \brief Mean values of primitive variables. */
-  MeandRhodT,                      /*!< \brief Derivative of density w.r.t. temperature (variable density flows). */
+  MeandRhodh,                      /*!< \brief Derivative of density w.r.t. temperature (variable density flows). */
   Param_p, Param_Kappa_0,          /*!< \brief Artificial dissipation parameters. */
   Local_Lambda_i, Local_Lambda_j,
   MeanLambda,                      /*!< \brief Local eingenvalues. */
@@ -109,7 +109,7 @@ class CCentJSTInc_Flow final : public CNumerics {
   MeanDensity, MeanPressure,
   MeanBetaInc2, MeanEnthalpy,
   MeanCp, MeanTemperature,        /*!< \brief Mean values of primitive variables. */
-  MeandRhodT,                     /*!< \brief Derivative of density w.r.t. temperature (variable density flows). */
+  MeandRhodh,                     /*!< \brief Derivative of density w.r.t. temperature (variable density flows). */
   Param_p, Param_Kappa_2,
   Param_Kappa_4,                  /*!< \brief Artificial dissipation parameters. */
   Local_Lambda_i, Local_Lambda_j,

SU2_CFD/include/numerics/flow/convection/fds.hpp

@@ -49,7 +49,7 @@ class CUpwFDSInc_Flow final : public CNumerics {
   su2double **Precon, **invPrecon_A;
   su2double Proj_ModJac_Tensor_ij, Pressure_i,
   Pressure_j, ProjVelocity,
-  MeandRhodT, dRhodT_i, dRhodT_j, /*!< \brief Derivative of density w.r.t. temperature (variable density flows). */
+  MeandRhodh, dRhodh_i, dRhodh_j, /*!< \brief Derivative of density w.r.t. enthalpy (variable density flows). */
   Temperature_i, Temperature_j,   /*!< \brief Temperature at node 0 and 1. */
   MeanDensity, MeanPressure, MeanSoundSpeed, MeanBetaInc2, MeanEnthalpy, MeanCp, MeanTemperature; /*!< \brief Mean values of primitive variables. */
   unsigned short iDim, iVar, jVar, kVar;

SU2_CFD/include/numerics/flow/flow_diffusion.hpp

@@ -290,6 +290,7 @@ class CAvgGrad_Flow final : public CAvgGrad_Base {
 class CAvgGradInc_Flow final : public CAvgGrad_Base {
 private:
   su2double Mean_Thermal_Conductivity; /*!< \brief Mean value of the effective thermal conductivity. */
+  su2double Mean_Heat_Capacity;        /*!< \brief Mean value of the heat capacity. */
   bool energy;                         /*!< \brief computation with the energy equation. */
 
   /*!

SU2_CFD/include/numerics/flow/flow_sources.hpp

@@ -126,7 +126,8 @@ class CSourceGeneralAxisymmetric_Flow final : public CSourceAxisymmetric_Flow {
 class CSourceIncAxisymmetric_Flow final : public CSourceBase_Flow {
   bool implicit, /*!< \brief Implicit calculation. */
   viscous,       /*!< \brief Viscous incompressible flows. */
-  energy;        /*!< \brief computation with the energy equation. */
+  energy,        /*!< \brief computation with the energy equation. */
+  multicomponent; /*!< \brief multicomponent incompressible flows. */
 
 public:
   /*!

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -69,7 +69,8 @@ class CFVMFlowSolverBase : public CSolver {
   su2double Mach_Inf = 0.0;          /*!< \brief Mach number at the infinity. */
   su2double Density_Inf = 0.0;       /*!< \brief Density at the infinity. */
   su2double Energy_Inf = 0.0;        /*!< \brief Energy at the infinity. */
-  su2double Temperature_Inf = 0.0;   /*!< \brief Energy at the infinity. */
+  su2double Temperature_Inf = 0.0;   /*!< \brief Temperature at the infinity. */
+  su2double Enthalpy_Inf = 0.0;      /*!< \brief Enthalpy at the infinity. */
   su2double Pressure_Inf = 0.0;      /*!< \brief Pressure at the infinity. */
   su2double* Velocity_Inf = nullptr; /*!< \brief Flow Velocity vector at the infinity. */
 

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -203,6 +203,19 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
                       CConfig *config,
                       unsigned short iMesh) final;
 
+  /*!
+   * \brief Recompute the extrapolated quantities, after MUSCL reconstruction,
+   *        in a more thermodynamically consistent way.
+   * \note This method is static to improve the chances of it being used in a
+   *       thread-safe manner.
+   * \param[in,out] fluidModel - The fluid model.
+   * \param[in] nDim - Number of physical dimensions.
+   * \param[in,out] primitive - Primitive variables.
+   * \param[in] scalar - scalar variable.
+   */
+  static void ComputeConsistentExtrapolation(CFluidModel* fluidModel, unsigned short nDim, const su2double* scalar,
+                                             su2double* primitive);
+
   /*!
    * \brief Source term integration.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CIncNSSolver.hpp

@@ -94,8 +94,19 @@ class CIncNSSolver final : public CIncEulerSolver {
    */
   void Compute_Streamwise_Periodic_Recovered_Values(CConfig *config, const CGeometry *geometry,
                                                     const unsigned short iMesh);
+  /*!
+   * \brief Compute the enthalpy due to species diffusion contribution for multicomponent and reacting flows.
+   * \param[in] iEdge - Edge for which the flux and Jacobians are to be computed.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics - Description of the numerical method.
+   * \param[in] n_species - Number of species in the mixture that a transport equation is solved.
+   * \param[in] implicit - Boolean for implicit iterations.
+   */
+  void Compute_Enthalpy_Diffusion(unsigned long iEdge, CGeometry* geometry, CSolver** solver_container,
+                                  CNumerics* numerics, const int n_species, const bool implicit);
 
-public:
+ public:
   /*!
    * \brief Constructor of the class.
    * \param[in] geometry - Geometrical definition of the problem.