add diffusion coefficients as volume output

Author: Bot-Enigma-0

SU2_CFD/include/numerics/CNumerics.hpp

@@ -190,6 +190,8 @@ class CNumerics {
 
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
 
+  su2double DiffCoeff_kw[2];      /*!< \brief Storage for diffusion coefficient*/
+
 public:
   /*!
    * \brief Return type used in some "ComputeResidual" overloads to give a
@@ -706,6 +708,7 @@ class CNumerics {
    * \param[in] val_CDkw_i - Value of the cross diffusion at point i.
    */
   virtual void SetCrossDiff(su2double val_CDkw_i) {/* empty */};
+  inline su2double GetDiffCoeff_kw(int index) const {return DiffCoeff_kw[index];}
 
   /*!
    * \brief Set the value of the effective intermittency for the LM model.

SU2_CFD/include/variables/CTurbVariable.hpp

@@ -43,7 +43,8 @@ class CTurbVariable : public CScalarVariable {
   static constexpr size_t MAXNVAR = 2;
   VectorType turb_index;
   VectorType intermittency;         /*!< \brief Value of the intermittency for the trans. model. */
-
+  VectorType DC_TKE;    /*!< \brief TKE Diffusion Coefficient. */ 
+  VectorType DC_Omega;  /*!< \brief Omega Diffusion Coefficient. */
   /*!
    * \brief Constructor of the class.
    * \param[in] npoint - Number of points/nodes/vertices in the domain.
@@ -100,6 +101,31 @@ class CTurbVariable : public CScalarVariable {
    */
   inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; }
 
+  /*!
+   * \brief Set the Diffusion Coefficients of TKE and omega equations.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_DC_kw - diffusion coefficient value
+   */
+  inline void SetDiffCoeff_kw(unsigned long iPoint, su2double* val_DC_kw) { 
+    DC_TKE(iPoint) = val_DC_kw[0];
+    DC_Omega(iPoint) = val_DC_kw[1];
+  }
+
+  /*!
+   * \brief Get the diffusion coefficient value of TKE.
+   * \param[in] iPoint - Point index.
+   */
+  inline su2double GetTKE_DC(unsigned long iPoint) const final { 
+    return DC_TKE(iPoint); 
+  }  
+  /*!
+   * \brief Get the diffusion coefficient value of omega.
+   * \param[in] iPoint - Point index.
+   */
+  inline su2double GetOmega_DC(unsigned long iPoint) const final { 
+    return DC_Omega(iPoint); 
+  }
+
   /*!
    * \brief Register eddy viscosity (muT) as Input or Output of an AD recording.
    * \param[in] input - Boolean whether In- or Output should be registered.

SU2_CFD/include/variables/CVariable.hpp

@@ -1664,6 +1664,15 @@ class CVariable {
    */
   inline virtual su2double GetCrossDiff(unsigned long iPoint) const { return 0.0; }
 
+  /*!
+   * \brief Get the value of the diffusion coefficient of tke and omega equations.
+   * \param[in] ipoint 
+   * \param[in] val_diff_coeff
+   */
+  inline virtual void SetDiffCoeff_kw(unsigned long iPoint, su2double* val_DiffCoeff) const { }
+  inline virtual su2double GetTKE_DC(unsigned long iPoint) const {return 0.0; }
+  inline virtual su2double GetOmega_DC(unsigned long iPoint) const {return 0.0; }
+
   /*!
    * \brief Get the value of the eddy viscosity.
    * \return the value of the eddy viscosity.

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1256,6 +1256,8 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){
     case TURB_FAMILY::KW:
       AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy");
       AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation");
+      AddVolumeOutput("DIFF_COEFF_TKE", "TKE_Diff_Coeff", "SOLUTION", "Diffusion Coefficient for TKE equation");
+      AddVolumeOutput("DIFF_COEFF_Omega", "Omega_Diff_Coeff", "SOLUTION", "Diffusion Coefficient for Omega equation");
       break;
 
     case TURB_FAMILY::NONE:
@@ -1536,6 +1538,8 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     case TURB_FAMILY::KW:
       SetVolumeOutputValue("TKE", iPoint, Node_Turb->GetSolution(iPoint, 0));
       SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1));
+      SetVolumeOutputValue("DIFF_COEFF_TKE", iPoint, Node_Turb->GetTKE_DC(iPoint));
+      SetVolumeOutputValue("DIFF_COEFF_Omega", iPoint, Node_Turb->GetOmega_DC(iPoint));
       SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0));
       SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1));
       if (limiter) {

SU2_CFD/src/solvers/CTurbSSTSolver.cpp

@@ -386,6 +386,10 @@ void CTurbSSTSolver::Viscous_Residual(const unsigned long iEdge, const CGeometry
         }
     }
   }
+  su2double DC_kw[2];
+  DC_kw[0] = numerics->GetDiffCoeff_kw(0);
+  DC_kw[1] = numerics->GetDiffCoeff_kw(1);
+  nodes->SetDiffCoeff_kw(iPoint, DC_kw);
 }
 
 void CTurbSSTSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container,

SU2_CFD/src/variables/CTurbVariable.cpp

@@ -35,6 +35,10 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned
     turb_index.resize(nPoint) = su2double(1.0);
     intermittency.resize(nPoint) = su2double(1.0);
 
+    if (TurbModelFamily(config->GetKind_Turb_Model()) == TURB_FAMILY::KW) {
+      DC_TKE.resize(nPoint) = su2double(1.0);
+      DC_Omega.resize(nPoint) = su2double(1.0);
+    }
    }
 
 void CTurbVariable::RegisterEddyViscosity(bool input) {