Axisymmetric problems : Viscous source terms and generalised convective source terms and jacobian

Common/include/option_structure.hpp

@@ -102,6 +102,7 @@ const su2double EPS = 1.0E-16;        /*!< \brief Error scale. */
 const su2double TURB_EPS = 1.0E-16;   /*!< \brief Turbulent Error scale. */
 
 const su2double ONE2 = 0.5;         /*!< \brief One divided by two. */
+const su2double ONE3 = 1.0 / 3.0;   /*!< \brief One divided by three. */
 const su2double TWO3 = 2.0 / 3.0;   /*!< \brief Two divided by three. */
 const su2double FOUR3 = 4.0 / 3.0;  /*!< \brief Four divided by three. */
 
@@ -2165,23 +2166,24 @@ enum MPI_QUANTITIES {
   MAX_EIGENVALUE       = 11,  /*!< \brief Maximum eigenvalue communication. */
   SENSOR               = 12,  /*!< \brief Dissipation sensor communication. */
   AUXVAR_GRADIENT      = 13,  /*!< \brief Auxiliary variable gradient communication. */
-  COORDINATES          = 14,  /*!< \brief Vertex coordinates communication. */
-  COORDINATES_OLD      = 15,  /*!< \brief Old vertex coordinates communication. */
-  MAX_LENGTH           = 16,  /*!< \brief Maximum length communication. */
-  GRID_VELOCITY        = 17,  /*!< \brief Grid velocity communication. */
-  CROSS_TERM           = 18,  /*!< \brief Cross term communication. */
-  CROSS_TERM_GEOMETRY  = 19,  /*!< \brief Geometric cross term communication. */
-  REF_GEOMETRY         = 20,  /*!< \brief Reference geometry communication. */
-  SOLUTION_EDDY        = 21,  /*!< \brief Turbulent solution plus eddy viscosity communication. */
-  SOLUTION_MATRIX      = 22,  /*!< \brief Matrix solution communication. */
-  SOLUTION_MATRIXTRANS = 23,  /*!< \brief Matrix transposed solution communication. */
-  NEIGHBORS            = 24,  /*!< \brief Neighbor point count communication (for JST). */
-  SOLUTION_FEA         = 25,  /*!< \brief FEA solution communication. */
-  SOLUTION_FEA_OLD     = 26,  /*!< \brief FEA solution old communication. */
-  MESH_DISPLACEMENTS   = 27,  /*!< \brief Mesh displacements at the interface. */
-  SOLUTION_TIME_N      = 28,  /*!< \brief Solution at time n. */
-  SOLUTION_TIME_N1     = 29,   /*!< \brief Solution at time n-1. */
-  PRIMITIVE            = 30   /*!< \brief Primitive solution communication. */
+  AXIAUXVAR_GRADIENT   = 14,  /*!< \brief Axisymmetric auxiliary variable gradient communication. */
+  COORDINATES          = 15,  /*!< \brief Vertex coordinates communication. */
+  COORDINATES_OLD      = 16,  /*!< \brief Old vertex coordinates communication. */
+  MAX_LENGTH           = 17,  /*!< \brief Maximum length communication. */
+  GRID_VELOCITY        = 18,  /*!< \brief Grid velocity communication. */
+  CROSS_TERM           = 19,  /*!< \brief Cross term communication. */
+  CROSS_TERM_GEOMETRY  = 20,  /*!< \brief Geometric cross term communication. */
+  REF_GEOMETRY         = 21,  /*!< \brief Reference geometry communication. */
+  SOLUTION_EDDY        = 22,  /*!< \brief Turbulent solution plus eddy viscosity communication. */
+  SOLUTION_MATRIX      = 23,  /*!< \brief Matrix solution communication. */
+  SOLUTION_MATRIXTRANS = 24,  /*!< \brief Matrix transposed solution communication. */
+  NEIGHBORS            = 25,  /*!< \brief Neighbor point count communication (for JST). */
+  SOLUTION_FEA         = 26,  /*!< \brief FEA solution communication. */
+  SOLUTION_FEA_OLD     = 27,  /*!< \brief FEA solution old communication. */
+  MESH_DISPLACEMENTS   = 28,  /*!< \brief Mesh displacements at the interface. */
+  SOLUTION_TIME_N      = 29,  /*!< \brief Solution at time n. */
+  SOLUTION_TIME_N1     = 30,  /*!< \brief Solution at time n-1. */
+  PRIMITIVE            = 31   /*!< \brief Primitive solution communication. */
 };
 
 /*!

SU2_CFD/include/numerics/CNumerics.hpp

@@ -183,6 +183,9 @@ class CNumerics {
   su2double
   *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
+  **AxiAuxVar_Grad_i,    /*!< \brief Gradient of axisymmetric auxiliary variable at point i. */
+  **AxiAuxVar_Grad_j;    /*!< \brief Gradient of axisymmetric auxiliary variable at point i. */
   const su2double *RadVar_Source;  /*!< \brief Source term from the radiative heat transfer equation. */
   su2double
   *Coord_i,      /*!< \brief Cartesians coordinates of point i. */
@@ -490,6 +493,17 @@ class CNumerics {
     AuxVar_Grad_j = val_auxvargrad_j;
   }
 
+  /*!
+   * \brief Set the gradient of the axisymmetric variables.
+   * \param[in] val_primvar_grad_i - Gradient of the primitive variable at point i.
+   * \param[in] val_primvar_grad_j - Gradient of the primitive variable at point j.
+   */
+  inline void SetAxiAuxVarGrad(su2double **val_auxvar_grad_i,
+                               su2double **val_auxvar_grad_j) {
+    AxiAuxVar_Grad_i = val_auxvar_grad_i;
+    AxiAuxVar_Grad_j = val_auxvar_grad_j;
+  }
+
   /*!
    * \brief Set the diffusion coefficient
    * \param[in] val_diffusioncoeff_i - Value of the diffusion coefficients at i.

SU2_CFD/include/numerics/flow/flow_sources.hpp

@@ -63,7 +63,16 @@ class CSourceBase_Flow : public CNumerics {
  * \ingroup SourceDiscr
  * \author F. Palacios
  */
-class CSourceAxisymmetric_Flow final : public CSourceBase_Flow {
+class CSourceAxisymmetric_Flow : public CSourceBase_Flow {
+protected:
+    bool implicit, viscous;
+    su2double yinv{0.0};
+
+  /*!
+  * \brief Diffusion residual of the axisymmetric source term.
+  */
+  void ResidualDiffusion();
+
 public:
   /*!
    * \brief Constructor of the class.
@@ -74,12 +83,31 @@ class CSourceAxisymmetric_Flow final : public CSourceBase_Flow {
   CSourceAxisymmetric_Flow(unsigned short val_nDim, unsigned short val_nVar, const CConfig* config);
 
   /*!
-   * \brief Residual of the rotational frame source term.
+   * \brief Residual of the axisymmetric source term.
    * \param[in] config - Definition of the particular problem.
    * \return Lightweight const-view of residual and Jacobian.
    */
   ResidualType<> ComputeResidual(const CConfig* config) override;
+
+};
 
+/*!
+ * \class CSourceGeneralAxisymmetric_Flow
+ * \brief Class for source term for solving axisymmetric problems for a general (non ideal) fluid.
+ * \ingroup SourceDiscr
+ * \author F. Dittmann
+ */
+class CSourceGeneralAxisymmetric_Flow final : public CSourceAxisymmetric_Flow {
+public:
+
+  using CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow;
+  /*!
+   * \brief Residual of the general axisymmetric source term.
+   * \param[in] config - Definition of the particular problem.
+   * \return Lightweight const-view of residual and Jacobian.
+   */
+  ResidualType<> ComputeResidual(const CConfig* config) override;
+
 };
 
 /*!

SU2_CFD/include/solvers/CSolver.hpp

@@ -592,6 +592,26 @@ class CSolver {
    */
   void SetAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config);
 
+  /*!
+   * \brief Compute the Green-Gauss gradient of the auxiliary variable.
+   * \param[in] geometry - Geometrical definition of the problem.
+   */
+  void SetAxiAuxVar_Gradient_GG(CGeometry *geometry, const CConfig *config);
+
+  /*!
+   * \brief Compute the Least Squares gradient of the auxiliary variable.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetAxiAuxVar_Gradient_LS(CGeometry *geometry, const CConfig *config);
+
+  /*!
+   * \brief Compute the Least Squares gradient of an auxiliar variable on the profile surface.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void SetAxiAuxVar_Surface_Gradient(CGeometry *geometry, const CConfig *config);
+
   /*!
    * \brief Add External to Solution vector.
    */

SU2_CFD/include/variables/CEulerVariable.hpp

@@ -205,6 +205,12 @@ class CEulerVariable : public CVariable {
    */
   inline su2double **GetGradient_Reconstruction(unsigned long iPoint) final { return Gradient_Reconstruction[iPoint]; }
 
+  /*!
+   * \brief Get the value of the primitive variables gradient.
+   * \return Value of the primitive variables gradient.
+   */
+  inline su2double **GetAxiAuxVarGradient(unsigned long iPoint) final { return Grad_AxiAuxVar[iPoint]; }
+
   /*!
    * \brief A virtual member.
    */

SU2_CFD/include/variables/CVariable.hpp

@@ -80,6 +80,9 @@ class CVariable {
   VectorType AuxVar;       /*!< \brief Auxiliar variable for gradient computation. */
   MatrixType Grad_AuxVar;  /*!< \brief Gradient of the auxiliar variable. */
 
+  MatrixType AxiAuxVar;             /*!< \brief Axisymmetric auxiliar variable for gradient computation. */
+  CVectorOfMatrix Grad_AxiAuxVar;   /*!< \brief Gradient of the axisymmetric variables  of the problem. */
+
   VectorType Max_Lambda_Inv;   /*!< \brief Maximun inviscid eingenvalue. */
   VectorType Max_Lambda_Visc;  /*!< \brief Maximun viscous eingenvalue. */
   VectorType Lambda;           /*!< \brief Value of the eingenvalue. */
@@ -624,6 +627,45 @@ class CVariable {
    */
   inline su2double GetAuxVarGradient(unsigned long iPoint, unsigned long iDim) const { return Grad_AuxVar(iPoint,iDim); }
 
+  /*!
+   * \brief Set axisymmetric auxiliar variables.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_auxvar - Value of the auxiliar variable.
+   */
+  inline void SetAxiAuxVar(unsigned long iPoint, const su2double *axiauxvar) {
+    for (unsigned long iVar = 0; iVar < 3; iVar++) AxiAuxVar(iPoint,iVar) = axiauxvar[iVar];
+  }
+
+  /*!
+   * \brief Get the entire axi aux vector of the problem.
+   * \return Reference to the axi aux matrix.
+   */
+  inline const MatrixType& GetAxiAuxVar(void) { return AxiAuxVar; }
+
+  /*!
+   * \brief Set value of axisymmetric auxillary gradients.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Index of the variable.
+   * \param[in] iDim - Index of the dimension.
+   * \param[in] value - Value of the gradient.
+   */
+  inline void SetAxiAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) { Grad_AxiAuxVar(iPoint,iVar,iDim) = value; }
+
+  /*!
+   * \brief Get the gradient of the axi auxilary variables.
+   * \return Reference to gradient.
+   */
+  inline CVectorOfMatrix& GetAxiAuxVarGradient(void) { return Grad_AxiAuxVar; }
+
+  /*!
+   * \brief Get the value of the axisymmetric auxilliary gradient.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Index of the variable.
+   * \param[in] iDim - Index of the dimension.
+   * \return Value of the solution gradient.
+   */
+  inline su2double GetAxiAuxVarGradient(unsigned long iPoint, unsigned long iVar, unsigned long iDim) const { return Grad_AxiAuxVar(iPoint,iVar,iDim); }
+
   /*!
    * \brief Add a value to the truncation error.
    * \param[in] iPoint - Point index.
@@ -1955,6 +1997,12 @@ class CVariable {
    */
   inline virtual CVectorOfMatrix& GetGradient_Reconstruction(void) { return Gradient; }
 
+  /*!
+   * \brief A virtual member.
+   * \return Value of the primitive variables gradient.
+   */
+  inline virtual su2double **GetAxiAuxVarGradient(unsigned long iPoint) { return nullptr; }
+
   /*!
    * \brief Set the blending function for the blending of k-w and k-eps.
    * \param[in] val_viscosity - Value of the vicosity.

SU2_CFD/src/drivers/CDriver.cpp

@@ -1865,8 +1865,10 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol
       else if (config->GetAxisymmetric() == YES) {
         if (incompressible)
           numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceIncAxisymmetric_Flow(nDim, nVar_Flow, config);
-        else
+        else if (ideal_gas)
           numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceAxisymmetric_Flow(nDim, nVar_Flow, config);
+        else
+          numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGeneralAxisymmetric_Flow(nDim, nVar_Flow, config);
       }
       else if (config->GetGravityForce() == YES) {
         numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGravity(nDim, nVar_Flow, config);