fix old bug in SA-neg diffusion and abstract "velocity gradient"

SU2_CFD/include/numerics/turbulent/turb_diffusion.hpp

@@ -70,7 +70,6 @@ class CAvgGrad_TurbSA final : public CAvgGrad_Scalar<FlowIndices> {
     su2double nu_j = Laminar_Viscosity_j/Density_j;
     su2double nu_e = 0.5*(nu_i+nu_j+ScalarVar_i[0]+ScalarVar_j[0]);
 
-    /// TODO: I smell bug...
     Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma;
 
     /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
@@ -111,7 +110,6 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
   using Base::ScalarVar_i;
   using Base::ScalarVar_j;
   using Base::Proj_Mean_GradScalarVar;
-  using Base::Proj_Mean_GradScalarVar_Normal;
   using Base::proj_vector_ij;
   using Base::implicit;
   using Base::Flux;
@@ -150,7 +148,7 @@ class CAvgGrad_TurbSA_Neg final : public CAvgGrad_Scalar<FlowIndices> {
       nu_e = nu_ij + fn*nu_tilde_ij;
     }
 
-    Flux[0] = nu_e*Proj_Mean_GradScalarVar_Normal[0]/sigma;
+    Flux[0] = nu_e*Proj_Mean_GradScalarVar[0]/sigma;
 
     /*--- For Jacobians -> Use of TSL approx. to compute derivatives of the gradients ---*/
 

SU2_CFD/include/variables/CEulerVariable.hpp

@@ -229,6 +229,14 @@ class CEulerVariable : public CFlowVariable {
     return Primitive(iPoint,iDim+indices.Velocity());
   }
 
+  /*!
+   * \brief Get the velocity gradient.
+   * \return Value of the velocity gradient.
+   */
+  inline CMatrixView<const su2double> GetVelocityGradient(unsigned long iPoint) const final {
+    return Gradient_Primitive(iPoint, indices.Velocity());
+  }
+
   /*!
    * \brief Get the projected velocity in a unitary vector direction (compressible solver).
    * \param[in] val_vector - Direction of projection.

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -153,6 +153,14 @@ class CIncEulerVariable : public CFlowVariable {
     return Primitive(iPoint, iDim+indices.Velocity());
   }
 
+  /*!
+   * \brief Get the velocity gradient.
+   * \return Value of the velocity gradient.
+   */
+  inline CMatrixView<const su2double> GetVelocityGradient(unsigned long iPoint) const final {
+    return Gradient_Primitive(iPoint, indices.Velocity());
+  }
+
   /*!
    * \brief Get the projected velocity in a unitary vector direction (compressible solver).
    * \param[in] val_vector - Direction of projection.

SU2_CFD/include/variables/CNEMOEulerVariable.hpp

@@ -235,6 +235,14 @@ class CNEMOEulerVariable : public CFlowVariable {
    */
   inline su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const final { return Primitive(iPoint,VEL_INDEX+iDim); }
 
+  /*!
+   * \brief Get the velocity gradient.
+   * \return Value of the velocity gradient.
+   */
+  inline CMatrixView<const su2double> GetVelocityGradient(unsigned long iPoint) const final {
+    return Gradient_Primitive(iPoint, indices.Velocity());
+  }
+
   /*!
    * \brief Get the projected velocity in a unitary vector direction (compressible solver).
    * \param[in] val_vector - Direction of projection.

SU2_CFD/include/variables/CVariable.hpp

@@ -1042,6 +1042,15 @@ class CVariable {
    */
   inline virtual su2double GetVelocity(unsigned long iPoint, unsigned long iDim) const { return 0.0; }
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] iPoint - Point index.
+   * \return Value of the velocity gradient.
+   */
+  inline virtual CMatrixView<const su2double> GetVelocityGradient(unsigned long iPoint) const {
+    return CMatrixView<const su2double>();
+  }
+
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.

SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp

@@ -189,7 +189,7 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry
     su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow);
 
     su2double tau[3][3];
-    CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow,1), Viscosity);
+    CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetVelocityGradient(Point_Flow), Viscosity);
     for (auto iVar = 0u; iVar < nVar; iVar++) {
       for (auto jVar = 0u; jVar < nVar; jVar++) {
         // Viscous component in the tn vector --> Units of force (non-dimensional).

SU2_CFD/src/output/CFlowCompOutput.cpp

@@ -557,7 +557,7 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv
     } else {
       SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]);
     }
-    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1)));
+    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetVelocityGradient(iPoint)));
   }
 
   LoadCommonFVMOutputs(config, geometry, iPoint);

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -659,7 +659,7 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
     } else {
       SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]);
     }
-    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetGradient_Primitive(iPoint,1)));
+    SetVolumeOutputValue("Q_CRITERION", iPoint, GetQ_Criterion(Node_Flow->GetVelocityGradient(iPoint)));
   }
 
   // Streamwise Periodicity

SU2_CFD/src/solvers/CSolver.cpp

@@ -3984,7 +3984,7 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, const CConfig *config)
         if (viscous_flow) {
           su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint);
           su2double Tau[3][3];
-          CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetGradient_Primitive(iPoint)+1, Viscosity);
+          CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetVelocityGradient(iPoint), Viscosity);
           for (iDim = 0; iDim < nDim; iDim++) {
             auxForce[iDim] += GeometryToolbox::DotProduct(nDim, Tau[iDim], iNormal);
           }

SU2_CFD/src/solvers/CTurbSASolver.cpp

@@ -1248,7 +1248,7 @@ void CTurbSASolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contain
       visc_numerics->SetScalarVar(nodes->GetSolution(iPoint), &nu_tilde);
 
       visc_numerics->SetScalarVarGradient(nodes->GetGradient(iPoint),
-                                        nodes->GetGradient(iPoint));
+                                          nodes->GetGradient(iPoint));
 
       /*--- Compute residual, and Jacobians ---*/
 
@@ -1364,7 +1364,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
     const auto coord_i       = geometry->nodes->GetCoord(iPoint);
     const auto nNeigh        = geometry->nodes->GetnPoint(iPoint);
     const auto wallDistance  = geometry->nodes->GetWall_Distance(iPoint);
-    const auto primVarGrad   = flowNodes->GetGradient_Primitive(iPoint);
+    const auto velocityGrad  = flowNodes->GetVelocityGradient(iPoint);
     const auto vorticity     = flowNodes->GetVorticity(iPoint);
     const auto density       = flowNodes->GetDensity(iPoint);
     const auto laminarViscosity = flowNodes->GetLaminarViscosity(iPoint);
@@ -1375,7 +1375,7 @@ void CTurbSASolver::SetDES_LengthScale(CSolver **solver, CGeometry *geometry, CC
     su2double uijuij = 0.0;
     for(auto iDim = 0u; iDim < nDim; iDim++){
       for(auto jDim = 0u; jDim < nDim; jDim++){
-        uijuij += primVarGrad[1+iDim][jDim]*primVarGrad[1+iDim][jDim];
+        uijuij += pow(velocityGrad[iDim][jDim], 2);
       }
     }
     uijuij = sqrt(fabs(uijuij));