Heat solver using scalar framework - Part 1

Common/include/option_structure.hpp

@@ -272,22 +272,8 @@ static const MapType<std::string, MAIN_SOLVER> Solver_Map = {
   MakePair("FEM_LES", MAIN_SOLVER::FEM_LES)
   MakePair("NEMO_EULER",MAIN_SOLVER::NEMO_EULER)
   MakePair("NEMO_NAVIER_STOKES",MAIN_SOLVER::NEMO_NAVIER_STOKES)
-  MakePair("ADJ_EULER", MAIN_SOLVER::ADJ_EULER)
-  MakePair("ADJ_NAVIER_STOKES", MAIN_SOLVER::ADJ_NAVIER_STOKES)
-  MakePair("ADJ_RANS", MAIN_SOLVER::ADJ_RANS )
   MakePair("HEAT_EQUATION", MAIN_SOLVER::HEAT_EQUATION)
   MakePair("ELASTICITY", MAIN_SOLVER::FEM_ELASTICITY)
-  MakePair("DISC_ADJ_EULER", MAIN_SOLVER::DISC_ADJ_EULER)
-  MakePair("DISC_ADJ_RANS", MAIN_SOLVER::DISC_ADJ_RANS)
-  MakePair("DISC_ADJ_NAVIERSTOKES", MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES)
-  MakePair("DISC_ADJ_INC_EULER", MAIN_SOLVER::DISC_ADJ_INC_EULER)
-  MakePair("DISC_ADJ_INC_RANS", MAIN_SOLVER::DISC_ADJ_INC_RANS)
-  MakePair("DISC_ADJ_INC_NAVIERSTOKES", MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES)
-  MakePair("DISC_ADJ_HEAT_EQUATION", MAIN_SOLVER::DISC_ADJ_HEAT)
-  MakePair("DISC_ADJ_FEM_EULER", MAIN_SOLVER::DISC_ADJ_FEM_EULER)
-  MakePair("DISC_ADJ_FEM_RANS", MAIN_SOLVER::DISC_ADJ_FEM_RANS)
-  MakePair("DISC_ADJ_FEM_NS", MAIN_SOLVER::DISC_ADJ_FEM_NS)
-  MakePair("DISC_ADJ_FEM", MAIN_SOLVER::DISC_ADJ_FEM)
   MakePair("TEMPLATE_SOLVER", MAIN_SOLVER::TEMPLATE_SOLVER)
   MakePair("MULTIPHYSICS", MAIN_SOLVER::MULTIPHYSICS)
 };
@@ -689,7 +675,7 @@ enum class MIXINGVISCOSITYMODEL {
   DAVIDSON, /*!< \brief Davidson mixing viscosity model. */
 };
 static const MapType<std::string, MIXINGVISCOSITYMODEL> MixingViscosityModel_Map = {
-  MakePair("WILKE", MIXINGVISCOSITYMODEL::WILKE) 
+  MakePair("WILKE", MIXINGVISCOSITYMODEL::WILKE)
   MakePair("DAVIDSON", MIXINGVISCOSITYMODEL::DAVIDSON)
 };
 
@@ -1249,7 +1235,7 @@ inline LM_ParsedOptions ParseLMOptions(const LM_OPTIONS *LM_Options, unsigned sh
   };
 
   LMParsedOptions.LM2015 = IsPresent(LM_OPTIONS::LM2015);
-  
+
   int NFoundCorrelations = 0;
   if (IsPresent(LM_OPTIONS::MALAN)) {
     LMParsedOptions.Correlation = TURB_TRANS_CORRELATION::MALAN;

SU2_CFD/include/numerics/CNumerics.hpp

@@ -71,9 +71,6 @@ class CNumerics {
   Thermal_Conductivity_j,    /*!< \brief Thermal conductivity at point j. */
   Thermal_Conductivity_ve_i, /*!< \brief vibrational-electronic Thermal conductivity at point i. */
   Thermal_Conductivity_ve_j, /*!< \brief vibrational-electronic Thermal conductivity at point j. */
-  Thermal_Diffusivity_i,     /*!< \brief Thermal diffusivity at point i. */
-  Thermal_Diffusivity_j;     /*!< \brief Thermal diffusivity at point j. */
-  su2double
   Cp_i,               /*!< \brief Cp at point i. */
   Cp_j;               /*!< \brief Cp at point j. */
   su2double
@@ -158,7 +155,7 @@ class CNumerics {
   TurbPsi_Grad_j,  /*!< \brief Gradient of adjoint turbulent variables at point j. */
   AuxVar_Grad_i,   /*!< \brief Gradient of an auxiliary variable at point i. */
   AuxVar_Grad_j;   /*!< \brief Gradient of an auxiliary variable at point i. */
-  su2double 
+  su2double
   LocalGridLength_i; /*!< \brief Local grid length at point i. */
   const su2double *RadVar_Source;  /*!< \brief Source term from the radiative heat transfer equation. */
   const su2double
@@ -794,18 +791,6 @@ class CNumerics {
     Thermal_Conductivity_ve_j = val_thermal_conductivity_ve_j;
   }
 
-  /*!
-   * \brief Set the thermal diffusivity (translational/rotational)
-   * \param[in] val_thermal_diffusivity_i - Value of the thermal diffusivity at point i.
-   * \param[in] val_thermal_diffusivity_j - Value of the thermal diffusivity at point j.
-   * \param[in] iSpecies - Value of the species.
-   */
-  inline void SetThermalDiffusivity(su2double val_thermal_diffusivity_i,
-                                    su2double val_thermal_diffusivity_j) {
-    Thermal_Diffusivity_i = val_thermal_diffusivity_i;
-    Thermal_Diffusivity_j = val_thermal_diffusivity_j;
-  }
-
   /*!
    * \brief Set the specifc heat c_p.
    * \param[in] val_specific_heat_i - Value of the specific heat at point i.

SU2_CFD/include/numerics/heat.hpp

@@ -27,7 +27,7 @@
 
 #pragma once
 
-#include "CNumerics.hpp"
+#include "scalar/scalar_diffusion.hpp"
 
 /*!
  * \class CCentSca_Heat
@@ -100,26 +100,36 @@ class CUpwSca_Heat : public CNumerics {
  * \author O. Burghardt.
  * \version 7.5.0 "Blackbird"
  */
-class CAvgGrad_Heat : public CNumerics {
-private:
-  bool implicit, correct;
-
-public:
+class CAvgGrad_Heat final : public CAvgGrad_Scalar<CNoFlowIndices> {
+ public:
   /*!
    * \brief Constructor of the class.
    * \param[in] val_nDim - Number of dimensions of the problem.
-   * \param[in] val_nVar - Number of variables of the problem.
    * \param[in] config - Definition of the particular problem.
-   * \param[in] correct - Correct the gradient.
+   * \param[in] correct - Whether to correct the gradient.
    */
-  CAvgGrad_Heat(unsigned short val_nDim, unsigned short val_nVar, const CConfig *config, bool correct);
+  CAvgGrad_Heat(unsigned short val_nDim, const CConfig *config, bool correct)
+    : CAvgGrad_Scalar<CNoFlowIndices>(val_nDim, 1, correct, config) {}
 
+ private:
   /*!
-   * \brief Compute the viscous heat residual using an average of gradients with correction.
-   * \param[out] val_residual - Pointer to the total residual.
-   * \param[out] Jacobian_i - Jacobian of the numerical method at node i (implicit computation).
-   * \param[out] Jacobian_j - Jacobian of the numerical method at node j (implicit computation).
+   * \brief Adds extra variables to AD
+   */
+  void ExtraADPreaccIn(void) override {
+    AD::SetPreaccIn(*Diffusion_Coeff_i, *Diffusion_Coeff_j);
+  }
+
+  /*!
+   * \brief Heat-specific specific steps in the ComputeResidual method
    * \param[in] config - Definition of the particular problem.
    */
-  void ComputeResidual(su2double *val_residual, su2double **Jacobian_i, su2double **Jacobian_j, CConfig *config) override;
+  void FinishResidualCalc(const CConfig* config) override {
+    const su2double Thermal_Diffusivity_Mean = 0.5 * (*Diffusion_Coeff_i + *Diffusion_Coeff_j);
+
+    Flux[0] = Thermal_Diffusivity_Mean * Proj_Mean_GradScalarVar[0];
+
+    /*--- Use TSL for Jacobians. ---*/
+    Jacobian_i[0][0] = -Thermal_Diffusivity_Mean * proj_vector_ij;
+    Jacobian_j[0][0] = Thermal_Diffusivity_Mean * proj_vector_ij;
+  }
 };

SU2_CFD/include/numerics/scalar/scalar_diffusion.hpp

@@ -30,6 +30,18 @@
 
 #include "../CNumerics.hpp"
 
+/*!
+ * \class CNoFlowIndices
+ * \brief Dummy flow indices class to use CAvgGrad_Scalar when flow variables are not available.
+ * For example, solid heat transfer problems.
+ */
+struct CNoFlowIndices {
+  CNoFlowIndices(int, int) {}
+  inline int Density() const { return 0; }
+  inline int LaminarViscosity() const { return 0; }
+  inline int EddyViscosity() const { return 0; }
+};
+
 /*!
  * \class CAvgGrad_Scalar
  * \brief Template class for computing viscous residual of scalar values
@@ -57,7 +69,7 @@ class CAvgGrad_Scalar : public CNumerics {
   su2double* Jacobian_j[MAXNVAR];             /*!< \brief Flux Jacobian w.r.t. node j. */
   su2double JacobianBuffer[2*MAXNVAR*MAXNVAR];/*!< \brief Static storage for the two Jacobians. */
 
-  const bool correct_gradient = false, implicit = false, incompressible = false;
+  const bool correct_gradient = false, incompressible = false;
 
   /*!
    * \brief A pure virtual function; Adds any extra variables to AD
@@ -84,7 +96,6 @@ class CAvgGrad_Scalar : public CNumerics {
     : CNumerics(val_nDim, val_nVar, config),
       idx(val_nDim, config->GetnSpecies()),
       correct_gradient(correct_grad),
-      implicit(config->GetKind_TimeIntScheme_Turb() == EULER_IMPLICIT),
       incompressible(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) {
     if (nVar > MAXNVAR) {
       SU2_MPI::Error("Static arrays are too small.", CURRENT_FUNCTION);
@@ -116,15 +127,17 @@ class CAvgGrad_Scalar : public CNumerics {
       AD::SetPreaccIn(ScalarVar_i, nVar);
       AD::SetPreaccIn(ScalarVar_j, nVar);
     }
-    AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]);
-    AD::SetPreaccIn(V_j[idx.Density()], V_j[idx.LaminarViscosity()], V_j[idx.EddyViscosity()]);
-
-    Density_i = V_i[idx.Density()];
-    Density_j = V_j[idx.Density()];
-    Laminar_Viscosity_i = V_i[idx.LaminarViscosity()];
-    Laminar_Viscosity_j = V_j[idx.LaminarViscosity()];
-    Eddy_Viscosity_i = V_i[idx.EddyViscosity()];
-    Eddy_Viscosity_j = V_j[idx.EddyViscosity()];
+    if (!std::is_same<FlowIndices, CNoFlowIndices>::value) {
+      AD::SetPreaccIn(V_i[idx.Density()], V_i[idx.LaminarViscosity()], V_i[idx.EddyViscosity()]);
+      AD::SetPreaccIn(V_j[idx.Density()], V_j[idx.LaminarViscosity()], V_j[idx.EddyViscosity()]);
+
+      Density_i = V_i[idx.Density()];
+      Density_j = V_j[idx.Density()];
+      Laminar_Viscosity_i = V_i[idx.LaminarViscosity()];
+      Laminar_Viscosity_j = V_j[idx.LaminarViscosity()];
+      Eddy_Viscosity_i = V_i[idx.EddyViscosity()];
+      Eddy_Viscosity_j = V_j[idx.EddyViscosity()];
+    }
 
     ExtraADPreaccIn();
 

SU2_CFD/include/numerics/turbulent/transition/trans_diffusion.hpp

@@ -50,10 +50,9 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar<FlowIndices> {
   using Base::ScalarVar_j;
   using Base::Proj_Mean_GradScalarVar;
   using Base::proj_vector_ij;
-  using Base::implicit;
   using Base::Flux;
   using Base::Jacobian_i;
-  using Base::Jacobian_j;  
+  using Base::Jacobian_j;
 
   /*!
    * \brief Adds any extra variables to AD
@@ -65,7 +64,8 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   void FinishResidualCalc(const CConfig* config) override {
-
+    const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
+
     /*--- Compute mean effective dynamic viscosity ---*/
     const su2double diff_i_gamma = Laminar_Viscosity_i + Eddy_Viscosity_i;
     const su2double diff_j_gamma = Laminar_Viscosity_j + Eddy_Viscosity_j;
@@ -101,5 +101,5 @@ class CAvgGrad_TransLM final : public CAvgGrad_Scalar<FlowIndices> {
   CAvgGrad_TransLM(unsigned short val_nDim, unsigned short val_nVar, bool correct_grad, const CConfig* config)
     : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config){
   }
-  
+
 };

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

@@ -47,7 +47,6 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
   using Base::ScalarVar_j;
   using Base::Proj_Mean_GradScalarVar;
   using Base::proj_vector_ij;
-  using Base::implicit;
   using Base::Flux;
   using Base::Jacobian_i;
   using Base::Jacobian_j;
@@ -64,6 +63,8 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   void FinishResidualCalc(const CConfig* config) override {
+    const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
+
     /*--- Compute mean effective viscosity ---*/
 
     const su2double nu_i = Laminar_Viscosity_i/Density_i;
@@ -111,7 +112,6 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
   using Base::ScalarVar_j;
   using Base::Proj_Mean_GradScalarVar;
   using Base::proj_vector_ij;
-  using Base::implicit;
   using Base::Flux;
   using Base::Jacobian_i;
   using Base::Jacobian_j;
@@ -129,6 +129,8 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   void FinishResidualCalc(const CConfig* config) override {
+    const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
+
     /*--- Compute mean effective viscosity ---*/
 
     const su2double nu_i = Laminar_Viscosity_i/Density_i;
@@ -191,7 +193,6 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
   using Base::ScalarVar_j;
   using Base::Proj_Mean_GradScalarVar;
   using Base::proj_vector_ij;
-  using Base::implicit;
   using Base::Flux;
   using Base::Jacobian_i;
   using Base::Jacobian_j;
@@ -215,6 +216,8 @@ class CAvgGrad_TurbSST final : public CAvgGrad_Scalar<FlowIndices> {
    * \param[in] config - Definition of the particular problem.
    */
   void FinishResidualCalc(const CConfig* config) override {
+    const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT;
+
     /*--- Compute the blended constant for the viscous terms ---*/
     const su2double sigma_kine_i = F1_i*sigma_k1 + (1.0 - F1_i)*sigma_k2;
     const su2double sigma_kine_j = F1_j*sigma_k1 + (1.0 - F1_j)*sigma_k2;