cleanup code

Author: pcarruscag

Common/include/linear_algebra/CSysMatrix.hpp

@@ -553,6 +553,22 @@ class CSysMatrix {
     SetBlock<MatrixType, false>(block_i, block_j, val_block, alpha);
   }
 
+  /*!
+   * \brief Returns the 4 blocks ii, ij, ji, jj used by "UpdateBlocks".
+   * \note This method assumes an FVM-type sparse pattern.
+   * \param[in] edge - Index of edge that connects iPoint and jPoint.
+   * \param[in] iPoint - Row to which we add the blocks.
+   * \param[in] jPoint - Row from which we subtract the blocks.
+   * \param[out] bii, bij, bji, bjj - Blocks of the matrix.
+   */
+  inline void GetBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, ScalarType*& bii,
+                        ScalarType*& bij, ScalarType*& bji, ScalarType*& bjj) {
+    bii = &matrix[dia_ptr[iPoint] * nVar * nEqn];
+    bjj = &matrix[dia_ptr[jPoint] * nVar * nEqn];
+    bij = &matrix[edge_ptr(iEdge, 0) * nVar * nEqn];
+    bji = &matrix[edge_ptr(iEdge, 1) * nVar * nEqn];
+  }
+
   /*!
    * \brief Update 4 blocks ii, ij, ji, jj (add to i* sub from j*).
    * \note This method assumes an FVM-type sparse pattern.
@@ -566,10 +582,8 @@ class CSysMatrix {
   template <class MatrixType, class OtherType = ScalarType>
   inline void UpdateBlocks(unsigned long iEdge, unsigned long iPoint, unsigned long jPoint, const MatrixType& block_i,
                            const MatrixType& block_j, OtherType scale = 1) {
-    ScalarType* bii = &matrix[dia_ptr[iPoint] * nVar * nEqn];
-    ScalarType* bjj = &matrix[dia_ptr[jPoint] * nVar * nEqn];
-    ScalarType* bij = &matrix[edge_ptr(iEdge, 0) * nVar * nEqn];
-    ScalarType* bji = &matrix[edge_ptr(iEdge, 1) * nVar * nEqn];
+    ScalarType *bii, *bij, *bji, *bjj;
+    GetBlocks(iEdge, iPoint, jPoint, bii, bij, bji, bjj);
 
     unsigned long iVar, jVar, offset = 0;
 

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

@@ -54,6 +54,8 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
   const su2double sigma = 2.0/3.0;
   const su2double cb2 = 0.622;
 
+  const bool use_accurate_jacobians;
+
   /*!
    * \brief Adds any extra variables to AD
    */
@@ -68,41 +70,39 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
 
     /*--- Compute mean effective viscosity ---*/
 
-    /* First Term */
-    su2double nu_i = Laminar_Viscosity_i/Density_i;
-    su2double nu_j = Laminar_Viscosity_j/Density_j;
-    su2double nu_e = 0.5*(nu_i + nu_j + (1+cb2)*(ScalarVar_i[0] + ScalarVar_j[0])); 
-    su2double term_1 = nu_e; 
+    /*--- First Term. Normal diffusion, and conservative part of the quadratic diffusion.
+     * ||grad nu_t||^2 = div(nu_t grad nu_t) - nu_t div grad nu_t ---*/
+    const su2double nu_i = Laminar_Viscosity_i/Density_i;
+    const su2double nu_j = Laminar_Viscosity_j/Density_j;
+    const su2double nu_e = 0.5 * (nu_i + nu_j + (1 + cb2) * (ScalarVar_i[0] + ScalarVar_j[0]));
+    const su2double term_1 = nu_e;
 
-    /* Second Term */
-    su2double nu_tilde_i = ScalarVar_i[0];
-    su2double term_2 = cb2 * nu_tilde_i;
+    /* Second Term (quadratic diffusion, non conservative). */
+    const su2double nu_tilde_i = ScalarVar_i[0];
+    const su2double term_2 = cb2 * nu_tilde_i;
 
-    Flux[0] = (term_1 - term_2)*Proj_Mean_GradScalarVar[0]/sigma;
-    su2double dTerm1_dnut_i = (1+cb2)*0.5;
-    su2double dTerm1_dnut_j = (1+cb2)*0.5;
+    const su2double diffusion_coefficient = term_1 - term_2;
+    Flux[0] = diffusion_coefficient * Proj_Mean_GradScalarVar[0] / sigma;
 
-    su2double dTerm2_dnut_i = cb2;
-    su2double dTerm2_dnut_j = 0.0;
+    if (implicit) {
+      /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
+      Jacobian_i[0][0] = -diffusion_coefficient * proj_vector_ij / sigma;
+      Jacobian_j[0][0] = diffusion_coefficient * proj_vector_ij / sigma;
 
-    su2double dDC_dnut_i = dTerm1_dnut_i - dTerm2_dnut_i;
-    su2double dDC_dnut_j = dTerm1_dnut_j - dTerm2_dnut_j;
+      if (use_accurate_jacobians) {
+        /*--- The diffusion coefficient is also a function of nu_t. ---*/
+        const su2double dTerm1_dnut_i = (1 + cb2) * 0.5;
+        const su2double dTerm1_dnut_j = (1 + cb2) * 0.5;
 
-    su2double dGrad_dnut_i = -proj_vector_ij;
-    su2double dGrad_dnut_j =  proj_vector_ij;
+        const su2double dTerm2_dnut_i = cb2;
+        const su2double dTerm2_dnut_j = 0.0;
 
-    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
-    su2double diffusion_coefficient = term_1 - term_2;
-    bool UseExactJacobians = config->GetUse_Accurate_Turb_Jacobians();
+        const su2double dDC_dnut_i = dTerm1_dnut_i - dTerm2_dnut_i;
+        const su2double dDC_dnut_j = dTerm1_dnut_j - dTerm2_dnut_j;
 
-    if (implicit) {
-      if (UseExactJacobians) {
-        Jacobian_i[0][0] = (dDC_dnut_i*Proj_Mean_GradScalarVar[0] + dGrad_dnut_i*diffusion_coefficient)/sigma; //exact
-        Jacobian_j[0][0] = (dDC_dnut_j*Proj_Mean_GradScalarVar[0] + dGrad_dnut_j*diffusion_coefficient)/sigma; //exact
-      } else {
-        Jacobian_i[0][0] = diffusion_coefficient*-proj_vector_ij/sigma; //frozen diffusion coefficients
-        Jacobian_j[0][0] = diffusion_coefficient* proj_vector_ij/sigma; //frozen diffusion coefficients
-      }    
+        Jacobian_i[0][0] += dDC_dnut_i * Proj_Mean_GradScalarVar[0] / sigma;
+        Jacobian_j[0][0] += dDC_dnut_j * Proj_Mean_GradScalarVar[0] / sigma;
+      }
     }
   }
 
@@ -116,7 +116,8 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
    */
   CAvgGrad_TurbSA(unsigned short val_nDim, unsigned short val_nVar,
                   bool correct_grad, const CConfig* config)
-    : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config) {}
+    : CAvgGrad_Scalar<FlowIndices>(val_nDim, val_nVar, correct_grad, config),
+      use_accurate_jacobians(config->GetUse_Accurate_Turb_Jacobians()) {}
 };
 
 /*!
@@ -162,42 +163,40 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
     const su2double nu_i = Laminar_Viscosity_i/Density_i;
     const su2double nu_j = Laminar_Viscosity_j/Density_j;
 
-    const su2double nu_ij = 0.5*(nu_i+nu_j);
+    const su2double nu_ij = 0.5 * (nu_i + nu_j);
     const su2double nu_tilde_i = ScalarVar_i[0];
     const su2double nu_tilde_j = ScalarVar_j[0];
-    const su2double nu_tilde_ij = 0.5*(nu_tilde_i + nu_tilde_j);
+    const su2double nu_tilde_ij = 0.5 * (nu_tilde_i + nu_tilde_j);
 
-    /* Following Diskin's implementation from 10.2514/1.J064629, they propose a new fn function
-       to be evaluated at the cell to maintain positivity in the diffusion coefficient, which is 
-       used in both terms. The new fn term averaged across the face reverts to the original fn
-       function. */ 
+    /*--- Following Diskin's implementation from 10.2514/1.J064629, they propose a new fn function
+     * to be evaluated at the cell to maintain positivity in the diffusion coefficient, which is
+     * used in both terms. The new fn term averaged across the face reverts to the original fn
+     * function. ---*/
 
     /*--- Second Term (LHS) ---*/
-    su2double zeta_i = ((1 + cb2)*nu_tilde_ij - cb2*nu_tilde_i)/nu_ij;
-    su2double fn_i;
-    if (zeta_i >= 0.0) {
-      fn_i = 1.0;  
-    } else {
-      fn_i = (cn1 + pow(zeta_i,3.0))/(cn1 - pow(zeta_i,3.0));
+    const su2double zeta_i = ((1 + cb2) * nu_tilde_ij - cb2 * nu_tilde_i) / nu_ij;
+    su2double fn_i = 1.0;
+    if (zeta_i < 0.0) {
+      fn_i = (cn1 + pow(zeta_i,3)) / (cn1 - pow(zeta_i,3));
     }
 
-    su2double term_1 = (nu_ij + (1 + cb2)*nu_tilde_ij*fn_i);
-    su2double term_2 = cb2*nu_tilde_i*fn_i;
-    Flux[0] = (term_1 - term_2)*Proj_Mean_GradScalarVar[0]/sigma;
+    const su2double term_1 = (nu_ij + (1 + cb2) * nu_tilde_ij * fn_i);
+    const su2double term_2 = cb2 * nu_tilde_i * fn_i;
+    Flux[0] = (term_1 - term_2) * Proj_Mean_GradScalarVar[0] / sigma;
 
-    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients 
+    /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients
     * Exact Jacobians were tested on multiple cases but resulted in divergence of all
     * simulations, hence only frozen diffusion coefficient (approximate) Jacobians are used. ---*/
 
-    su2double diffusion_coefficient = (term_1 - term_2);
+    if (implicit) {
+      const su2double diffusion_coefficient = (term_1 - term_2);
 
-    su2double dGrad_dnut_i = -proj_vector_ij;
-    su2double dGrad_dnut_j =  proj_vector_ij;
+      const su2double dGrad_dnut_i = -proj_vector_ij;
+      const su2double dGrad_dnut_j = proj_vector_ij;
 
-    if (implicit) { 
-        Jacobian_i[0][0] = diffusion_coefficient * dGrad_dnut_i/sigma;  //frozen diffusion
-        Jacobian_j[0][0] = diffusion_coefficient * dGrad_dnut_j/sigma;  //frozen diffusion
-      }
+      Jacobian_i[0][0] = diffusion_coefficient * dGrad_dnut_i / sigma;
+      Jacobian_j[0][0] = diffusion_coefficient * dGrad_dnut_j / sigma;
+    }
   }
 
 public:

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -86,7 +86,7 @@ class CSourceBase_TurbSA : public CNumerics {
    */
   inline void ResidualAxisymmetricDiffusion(su2double sigma) {
     if (Coord_i[1] < EPS) return;
-    
+
     const su2double yinv = 1.0 / Coord_i[1];
     const su2double& nue = ScalarVar_i[0];
 
@@ -250,7 +250,7 @@ class CSourceBase_TurbSA : public CNumerics {
       Residual = (Production - Destruction) * Volume;
 
       if (axisymmetric) ResidualAxisymmetricDiffusion(var.sigma);
-      
+
       Jacobian_i[0] *= Volume;
     }
 
@@ -424,6 +424,7 @@ struct Edw {
 
 /*!
  * \brief SA source terms classes: production and destruction term and their derivative.
+ * \note Quadratic diffusion is included in the viscous fluxes.
  * \ingroup SourceDiscr
  * \param[in] nue: SA variable.
  * \param[in] var: Common SA variables struct.
@@ -550,7 +551,7 @@ class CCompressibilityCorrection final : public ParentClass {
       const su2double v = V_i[idx.Velocity() + 1];
 
       const su2double d_axiCorrection = 2.0 * c5 * nue * pow(v * yinv / sound_speed, 2) * Volume;
-      const su2double axiCorrection = 0.5 * nue * d_axiCorrection; 
+      const su2double axiCorrection = 0.5 * nue * d_axiCorrection;
 
       this->Residual -= axiCorrection;
       this->Jacobian_i[0] -= d_axiCorrection;

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -77,7 +77,7 @@ class CScalarSolver : public CSolver {
 
   /*--- Edge fluxes for reducer strategy (see the notes in CEulerSolver.hpp). ---*/
   CSysVector<su2double> EdgeFluxes; /*!< \brief Flux across each edge. */
-  CSysVector<su2double> EdgeFluxes_Diff; /*!< \brief Flux difference across ij and ji for non-conservative discretisation. */
+  CSysVector<su2double> EdgeFluxesDiff; /*!< \brief Flux difference between ij and ji for non-conservative discretisation. */
 
   /*!
    * \brief The highest level in the variable hierarchy this solver can safely use.

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -341,6 +341,8 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
 
 template <class VariableType>
 void CScalarSolver<VariableType>::SumEdgeFluxes(CGeometry* geometry) {
+  const bool nonConservative = EdgeFluxesDiff.GetLocSize() > 0;
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
     LinSysRes.SetBlock_Zero(iPoint);
@@ -350,8 +352,8 @@ void CScalarSolver<VariableType>::SumEdgeFluxes(CGeometry* geometry) {
         LinSysRes.AddBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
       else {
         LinSysRes.SubtractBlock(iPoint, EdgeFluxes.GetBlock(iEdge));
-        if (EdgeFluxes_Diff.GetLocSize() > 0) {
-          LinSysRes.SubtractBlock(iPoint, EdgeFluxes_Diff.GetBlock(iEdge));
+        if (nonConservative) {
+          LinSysRes.SubtractBlock(iPoint, EdgeFluxesDiff.GetBlock(iEdge));
         }
       }
     }

SU2_CFD/include/solvers/CSolver.hpp

@@ -191,12 +191,12 @@ class CSolver {
   CSysVector<su2double> LinSysSol;    /*!< \brief vector to store iterative solution of implicit linear system. */
   CSysVector<su2double> LinSysRes;    /*!< \brief vector to store iterative residual of implicit linear system. */
 #ifndef CODI_FORWARD_TYPE
-  CSysMatrix<su2mixedfloat> Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */
-  CSysSolve<su2mixedfloat>  System;   /*!< \brief Linear solver/smoother. */
+  using JacobianScalarType = su2mixedfloat;
 #else
-  CSysMatrix<su2double> Jacobian;
-  CSysSolve<su2double>  System;
+  using JacobianScalarType = su2double;
 #endif
+  CSysMatrix<JacobianScalarType> Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */
+  CSysSolve<JacobianScalarType> System;    /*!< \brief Linear solver/smoother. */
 
   CSysVector<su2double> OutputVariables;    /*!< \brief vector to store the extra variables to be written. */
   string* OutputHeadingNames;               /*!< \brief vector of strings to store the headings for the exra variables */