diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp
index c1b23b4fc159..368fb4cd384b 100644
--- a/SU2_CFD/include/numerics/CNumerics.hpp
+++ b/SU2_CFD/include/numerics/CNumerics.hpp
@@ -459,6 +459,86 @@ class CNumerics {
     TurbPsi_Grad_j = val_turbpsivar_grad_j;
   }
 
+  /*!
+   * \brief Compute the mean rate of strain matrix.
+   * \details The parameter primvargrad can be e.g. PrimVar_Grad_i or Mean_GradPrimVar.
+   * \param[in] nDim - 2 or 3
+   * \param[out] rateofstrain - Rate of strain matrix
+   * \param[in] velgrad - A velocity gradient matrix.
+   * \tparam TWOINDICES_1 - any type that supports the [][] interface
+   * \tparam TWOINDICES_2 - any type that supports the [][] interface
+   */
+  template<class TWOINDICES_1, class TWOINDICES_2>
+  inline static void ComputeMeanRateOfStrainMatrix(unsigned short nDim, TWOINDICES_1& rateofstrain, const TWOINDICES_2& velgrad){
+
+    /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */
+
+    if (nDim == 3){
+      rateofstrain[0][0] = velgrad[0][0];
+      rateofstrain[1][1] = velgrad[1][1];
+      rateofstrain[2][2] = velgrad[2][2];
+      rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]);
+      rateofstrain[0][2] = 0.5 * (velgrad[0][2] + velgrad[2][0]);
+      rateofstrain[1][2] = 0.5 * (velgrad[1][2] + velgrad[2][1]);
+      rateofstrain[1][0] = rateofstrain[0][1];
+      rateofstrain[2][1] = rateofstrain[1][2];
+      rateofstrain[2][0] = rateofstrain[0][2];
+    }
+    else { // nDim==2
+      rateofstrain[0][0] = velgrad[0][0];
+      rateofstrain[1][1] = velgrad[1][1];
+      rateofstrain[2][2] = 0.0;
+      rateofstrain[0][1] = 0.5 * (velgrad[0][1] + velgrad[1][0]);
+      rateofstrain[0][2] = 0.0;
+      rateofstrain[1][2] = 0.0;
+      rateofstrain[1][0] = rateofstrain[0][1];
+      rateofstrain[2][1] = rateofstrain[1][2];
+      rateofstrain[2][0] = rateofstrain[0][2];
+    }
+  }
+
+  /*!
+   * \brief Compute the stress tensor from the velocity gradients.
+   * \details To obtain the Reynolds stress tensor +(u_i' u_j')~, divide the result
+   * of this function by (-rho). The argument density is only used if turb_ke is not 0.
+   * To select the velocity gradient components from a primitive variable gradient PrimVar_Grad_i,
+   * write PrimVar_Grad_i+1.
+   * If <code>nDim==2</code>, we use the same formula but only only access the entries [0][0]..[1][1] of
+   * stress and velgrad. If <code>reynolds3x3</code> is true, the other non-diagonal entries of stress
+   * set to zero, and <code>stress[2][2]</code> to some value.
+   * \param[in] nDim - Dimension of the flow problem, 2 or 3
+   * \param[out] stress - Stress tensor
+   * \param[in] velgrad - A velocity gradient matrix.
+   * \param[in] viscosity - Viscosity
+   * \param[in] density - Density
+   * \param[in] turb_ke - Turbulent kinetic energy, for the turbulent stress tensor
+   * \param[in] reynolds3x3 - If true, write to the third row and column of stress even if nDim==2.
+   * \tparam TWOINDICES_1 - any type that supports the [][] interface
+   * \tparam TWOINDICES_2 - any type that supports the [][] interface
+   */
+  template<class TWOINDICES_1, class TWOINDICES_2>
+  inline static void ComputeStressTensor(unsigned short nDim, TWOINDICES_1& stress, const TWOINDICES_2& velgrad,
+                                      su2double viscosity, su2double density=0.0, su2double turb_ke=0.0, bool reynolds3x3=false){
+    su2double divVel = 0;
+    for (unsigned short iDim = 0; iDim < nDim; iDim++){
+      divVel += velgrad[iDim][iDim];
+    }
+    su2double pTerm = 2./3. * (divVel * viscosity + density * turb_ke);
+
+    for (unsigned short iDim = 0; iDim < nDim; iDim++){
+      for (unsigned short jDim = 0; jDim < nDim; jDim++){
+        stress[iDim][jDim] = viscosity * (velgrad[iDim][jDim]+velgrad[jDim][iDim]);
+      }
+      stress[iDim][iDim] -= pTerm;
+    }
+
+    if(reynolds3x3 && nDim==2){ // fill the third row and column of Reynolds stress matrix
+      stress[0][2] = stress[1][2] = stress[2][0] = stress[2][1] = 0.0;
+      stress[2][2] = -pTerm;
+    }
+
+  }
+
   /*!
    * \brief Set the value of the first blending function.
    * \param[in] val_F1_i - Value of the first Menter blending function at point i.
diff --git a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp
index cf19cd0799da..c3131d48203d 100644
--- a/SU2_CFD/include/numerics/flow/flow_diffusion.hpp
+++ b/SU2_CFD/include/numerics/flow/flow_diffusion.hpp
@@ -190,12 +190,6 @@ class CAvgGrad_Base : public CNumerics {
    */
   void SetPerturbedRSM(su2double turb_ke, const CConfig* config);
 
-  /*!
-   * \brief Get the mean rate of strain matrix based on velocity gradients
-   * \param[in] S_ij
-   */
-  void GetMeanRateOfStrainMatrix(su2double **S_ij) const;
-
 public:
 
   /*!
diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
index 467b0dcb0c7a..70327066a197 100644
--- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
+++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp
@@ -339,12 +339,6 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
      */
   void SetPerturbedStrainMag(su2double turb_ke);
 
-  /*!
-   * \brief Get the mean rate of strain matrix based on velocity gradients
-   * \param[in] S_ij
-   */
-  void GetMeanRateOfStrainMatrix(su2double **S_ij);
-
 public:
   /*!
    * \brief Constructor of the class.
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index 3828b3a81237..94b8780ca034 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -2030,12 +2030,11 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
   unsigned long iVertex, iPoint, iPointNormal;
   unsigned short iMarker, iMarker_Monitoring, iDim, jDim;
   unsigned short T_INDEX = 0, TVE_INDEX = 0, VEL_INDEX = 0;
-  su2double Viscosity = 0.0, div_vel, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0,
+  su2double Viscosity = 0.0, WallDist[3] = {0.0}, Area, TauNormal, RefTemp, RefVel2 = 0.0,
             RefDensity = 0.0, GradTemperature, Density = 0.0, WallDistMod, FrictionVel, Mach2Vel, Mach_Motion,
             UnitNormal[3] = {0.0}, TauElem[3] = {0.0}, TauTangent[3] = {0.0}, Tau[3][3] = {{0.0}}, Cp,
             thermal_conductivity, thermal_conductivity_tr, thermal_conductivity_ve = 0.0,
-            MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor,
-            delta[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
+            MaxNorm = 8.0, Grad_Vel[3][3] = {{0.0}}, Grad_Temp[3] = {0.0}, AxiFactor;
   const su2double *Coord = nullptr, *Coord_Normal = nullptr, *Normal = nullptr;
   su2double **Grad_PrimVar = nullptr, dTn, dTven;
 
@@ -2190,16 +2189,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::Friction_Forces(const CGeometry* geometr
       }
 
       /*--- Evaluate Tau ---*/
-
-      div_vel = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim];
-
-      for (iDim = 0; iDim < nDim; iDim++) {
-        for (jDim = 0; jDim < nDim; jDim++) {
-          Tau[iDim][jDim] = Viscosity * (Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]) -
-                            TWO3 * Viscosity * div_vel * delta[iDim][jDim];
-        }
-      }
+      CNumerics::ComputeStressTensor(nDim, Tau, Grad_Vel, Viscosity);
 
       /*--- If necessary evaluate the QCR contribution to Tau ---*/
 
diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp
index 6ea613c5126d..8aa5415b4408 100644
--- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp
+++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp
@@ -188,22 +188,15 @@ void CFlowTractionInterface::GetDonor_Variable(CSolver *flow_solution, CGeometry
 
     su2double Viscosity = flow_nodes->GetLaminarViscosity(Point_Flow);
 
-    const su2double* const* GradVel = &flow_nodes->GetGradient_Primitive(Point_Flow)[1];
-
-    // Divergence of the velocity
-    su2double DivVel = 0.0;
-    for (auto iVar = 0u; iVar < nVar; iVar++) DivVel += GradVel[iVar][iVar];
-
+    su2double tau[3][3];
+    CNumerics::ComputeStressTensor(nVar, tau, flow_nodes->GetGradient_Primitive(Point_Flow)+1,Viscosity);
     for (auto iVar = 0u; iVar < nVar; iVar++) {
       for (auto jVar = 0u; jVar < nVar; jVar++) {
-        // Viscous stress
-        su2double delta_ij = (iVar == jVar);
-        su2double tau_ij = Viscosity*(GradVel[jVar][iVar] + GradVel[iVar][jVar] - TWO3*DivVel*delta_ij);
-
         // Viscous component in the tn vector --> Units of force (non-dimensional).
-        Donor_Variable[iVar] += tau_ij * Normal_Flow[jVar];
+        Donor_Variable[iVar] += tau[iVar][jVar] * Normal_Flow[jVar];
       }
     }
+
   }
 
   // Redimensionalize and take into account ramp transfer of the loads
diff --git a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp
index c443e35fa961..d2c61ea174be 100644
--- a/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp
+++ b/SU2_CFD/src/numerics/NEMO/CNEMONumerics.cpp
@@ -237,7 +237,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar,
   // rather than the standard V = [r1, ... , rn, T, Tve, ... ]
 
   unsigned short iSpecies, iVar, iDim, jDim;
-  su2double *Ds, *V, **GV, mu, ktr, kve, div_vel;
+  su2double *Ds, *V, **GV, mu, ktr, kve;
   su2double rho, T, Tve, RuSI, Ru;
   auto& Ms = fluidmodel->GetSpeciesMolarMass();
 
@@ -279,12 +279,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar,
   //Cpve = V[RHOCVVE_INDEX]+Ru/Mass;
   //kve += Cpve*(val_eddy_viscosity/Prandtl_Turb);
 
-  /*--- Calculate the velocity divergence ---*/
-  div_vel = 0.0;
-  for (iDim = 0 ; iDim < nDim; iDim++)
-    div_vel += GV[VEL_INDEX+iDim][iDim];
-
-
+  
   /*--- Pre-compute mixture quantities ---*/
   for (iDim = 0; iDim < nDim; iDim++) {
     Vector[iDim] = 0.0;
@@ -294,16 +289,7 @@ void CNEMONumerics::GetViscousProjFlux(su2double *val_primvar,
   }
 
   /*--- Compute the viscous stress tensor ---*/
-  for (iDim = 0; iDim < nDim; iDim++)
-    for (jDim = 0; jDim < nDim; jDim++)
-      tau[iDim][jDim] = 0.0;
-  for (iDim = 0 ; iDim < nDim; iDim++) {
-    for (jDim = 0 ; jDim < nDim; jDim++) {
-      tau[iDim][jDim] += mu * (val_gradprimvar[VEL_INDEX+jDim][iDim] +
-          val_gradprimvar[VEL_INDEX+iDim][jDim]);
-    }
-    tau[iDim][iDim] -= TWO3*mu*div_vel;
-  }
+  ComputeStressTensor(nDim,tau,val_gradprimvar+VEL_INDEX, mu);
 
   /*--- Populate entries in the viscous flux vector ---*/
   for (iDim = 0; iDim < nDim; iDim++) {
diff --git a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
index 98a50d9d9d53..5737b0f2d6ed 100644
--- a/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_diffusion.cpp
@@ -123,28 +123,22 @@ void CAvgGrad_Base::SetStressTensor(const su2double *val_primvar,
                            const su2double val_laminar_viscosity,
                            const su2double val_eddy_viscosity) {
 
-  unsigned short iDim, jDim;
   const su2double Density = val_primvar[nDim+2];
-  const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity;
-
-  su2double div_vel = 0.0;
-  for (iDim = 0 ; iDim < nDim; iDim++)
-    div_vel += val_gradprimvar[iDim+1][iDim];
 
-  /* --- If UQ methodology is used, calculate tau using the perturbed reynolds stress tensor --- */
+  /* --- If UQ methodology is used, use the perturbed Reynolds stress tensor
+   * for the turbulent part of tau. Otherwise both the laminar and turbulent
+   * parts of tau can be computed with the total viscosity. --- */
 
   if (using_uq){
-    for (iDim = 0 ; iDim < nDim; iDim++)
-      for (jDim = 0 ; jDim < nDim; jDim++)
-        tau[iDim][jDim] = val_laminar_viscosity*( val_gradprimvar[jDim+1][iDim] + val_gradprimvar[iDim+1][jDim] )
-        - TWO3*val_laminar_viscosity*div_vel*delta[iDim][jDim] - Density * MeanPerturbedRSM[iDim][jDim];
-
+    ComputeStressTensor(nDim, tau, val_gradprimvar+1, val_laminar_viscosity);  // laminar part
+    // add turbulent part which was perturbed
+    for (unsigned short iDim = 0 ; iDim < nDim; iDim++)
+      for (unsigned short jDim = 0 ; jDim < nDim; jDim++)
+        tau[iDim][jDim] += (-Density) * MeanPerturbedRSM[iDim][jDim];
   } else {
-
-    for (iDim = 0 ; iDim < nDim; iDim++)
-      for (jDim = 0 ; jDim < nDim; jDim++)
-        tau[iDim][jDim] = total_viscosity*( val_gradprimvar[jDim+1][iDim] + val_gradprimvar[iDim+1][jDim] )
-                        - TWO3*total_viscosity*div_vel*delta[iDim][jDim];
+    // compute both parts in one step
+    const su2double total_viscosity = val_laminar_viscosity + val_eddy_viscosity;
+    ComputeStressTensor(nDim, tau, val_gradprimvar+1, total_viscosity, Density, 0.0); // TODO why ignore turb_ke?
   }
 }
 
@@ -217,67 +211,14 @@ void CAvgGrad_Base::AddTauWall(const su2double *val_normal,
       tau[iDim][jDim] = tau[iDim][jDim]*(val_tau_wall/WallShearStress);
 }
 
-void CAvgGrad_Base::GetMeanRateOfStrainMatrix(su2double **S_ij) const
-{
-  /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */
-
-  if (nDim == 3){
-    S_ij[0][0] = Mean_GradPrimVar[1][0];
-    S_ij[1][1] = Mean_GradPrimVar[2][1];
-    S_ij[2][2] = Mean_GradPrimVar[3][2];
-    S_ij[0][1] = 0.5 * (Mean_GradPrimVar[1][1] + Mean_GradPrimVar[2][0]);
-    S_ij[0][2] = 0.5 * (Mean_GradPrimVar[1][2] + Mean_GradPrimVar[3][0]);
-    S_ij[1][2] = 0.5 * (Mean_GradPrimVar[2][2] + Mean_GradPrimVar[3][1]);
-    S_ij[1][0] = S_ij[0][1];
-    S_ij[2][1] = S_ij[1][2];
-    S_ij[2][0] = S_ij[0][2];
-  }
-  else {
-    S_ij[0][0] = Mean_GradPrimVar[1][0];
-    S_ij[1][1] = Mean_GradPrimVar[2][1];
-    S_ij[2][2] = 0.0;
-    S_ij[0][1] = 0.5 * (Mean_GradPrimVar[1][1] + Mean_GradPrimVar[2][0]);
-    S_ij[0][2] = 0.0;
-    S_ij[1][2] = 0.0;
-    S_ij[1][0] = S_ij[0][1];
-    S_ij[2][1] = S_ij[1][2];
-    S_ij[2][0] = S_ij[0][2];
-
-  }
-}
-
 void CAvgGrad_Base::SetReynoldsStressMatrix(su2double turb_ke){
-
-  unsigned short iDim, jDim;
-  su2double **S_ij = new su2double* [3];
-  su2double muT = Mean_Eddy_Viscosity;
-  su2double divVel = 0;
-  su2double density;
-  su2double TWO3 = 2.0/3.0;
-  density = Mean_PrimVar[nDim+2];
-
-  for (iDim = 0; iDim < 3; iDim++){
-    S_ij[iDim] = new su2double [3];
-  }
-
-  GetMeanRateOfStrainMatrix(S_ij);
-
-  /* --- Using rate of strain matrix, calculate Reynolds stress tensor --- */
-
-  for (iDim = 0; iDim < 3; iDim++){
-    divVel += S_ij[iDim][iDim];
-  }
-
-  for (iDim = 0; iDim < 3; iDim++){
-    for (jDim = 0; jDim < 3; jDim++){
-      MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim]
-      - muT / density * (2 * S_ij[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]);
+  su2double meandensity = Mean_PrimVar[nDim+2];
+  ComputeStressTensor(nDim, MeanReynoldsStress, Mean_GradPrimVar+1, Mean_Eddy_Viscosity, meandensity, turb_ke, true);
+  for(unsigned short iDim=0; iDim<3; iDim++){
+    for(unsigned short jDim=0; jDim<3; jDim++){
+      MeanReynoldsStress[iDim][jDim] /= (-meandensity);
     }
   }
-
-  for (iDim = 0; iDim < 3; iDim++)
-    delete [] S_ij[iDim];
-  delete [] S_ij;
 }
 
 void CAvgGrad_Base::SetPerturbedRSM(su2double turb_ke, const CConfig* config){
diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp
index a879554b5388..8d21f94676b0 100644
--- a/SU2_CFD/src/numerics/flow/flow_sources.cpp
+++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp
@@ -136,7 +136,7 @@ CSourceIncAxisymmetric_Flow::CSourceIncAxisymmetric_Flow(unsigned short val_nDim
 CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CConfig* config) {
 
   su2double yinv, Velocity_i[3];
-  unsigned short iDim, jDim, iVar, jVar;
+  unsigned short iDim, iVar, jVar;
 
   if (Coord_i[1] > EPS) {
 
@@ -197,21 +197,12 @@ CNumerics::ResidualType<> CSourceIncAxisymmetric_Flow::ComputeResidual(const CCo
       Eddy_Viscosity_i       = V_i[nDim+5];
       Thermal_Conductivity_i = V_i[nDim+6];
 
-      su2double total_viscosity, div_vel;
+      su2double total_viscosity;
 
       total_viscosity = (Laminar_Viscosity_i + Eddy_Viscosity_i);
 
       /*--- The full stress tensor is needed for variable density ---*/
-
-      div_vel = 0.0;
-      for (iDim = 0 ; iDim < nDim; iDim++)
-        div_vel += PrimVar_Grad_i[iDim+1][iDim];
-
-      for (iDim = 0 ; iDim < nDim; iDim++)
-        for (jDim = 0 ; jDim < nDim; jDim++)
-          tau[iDim][jDim] = (total_viscosity*(PrimVar_Grad_i[jDim+1][iDim] +
-                                              PrimVar_Grad_i[iDim+1][jDim] )
-                             -TWO3*total_viscosity*div_vel*delta[iDim][jDim]);
+      ComputeStressTensor(nDim, tau, PrimVar_Grad_i+1, total_viscosity);
 
       /*--- Viscous terms. ---*/
 
diff --git a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp
index dd10de58fb61..ea6d06a96c72 100644
--- a/SU2_CFD/src/numerics/turbulent/turb_sources.cpp
+++ b/SU2_CFD/src/numerics/turbulent/turb_sources.cpp
@@ -919,65 +919,13 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSST::ComputeResidual(const CConfi
 
 }
 
-void CSourcePieceWise_TurbSST::GetMeanRateOfStrainMatrix(su2double **S_ij)
-{
-    /* --- Calculate the rate of strain tensor, using mean velocity gradients --- */
-
-  if (nDim == 3){
-    S_ij[0][0] = PrimVar_Grad_i[1][0];
-    S_ij[1][1] = PrimVar_Grad_i[2][1];
-    S_ij[2][2] = PrimVar_Grad_i[3][2];
-    S_ij[0][1] = 0.5 * (PrimVar_Grad_i[1][1] + PrimVar_Grad_i[2][0]);
-    S_ij[0][2] = 0.5 * (PrimVar_Grad_i[1][2] + PrimVar_Grad_i[3][0]);
-    S_ij[1][2] = 0.5 * (PrimVar_Grad_i[2][2] + PrimVar_Grad_i[3][1]);
-    S_ij[1][0] = S_ij[0][1];
-    S_ij[2][1] = S_ij[1][2];
-    S_ij[2][0] = S_ij[0][2];
-  }
-  else {
-    S_ij[0][0] = PrimVar_Grad_i[1][0];
-    S_ij[1][1] = PrimVar_Grad_i[2][1];
-    S_ij[2][2] = 0.0;
-    S_ij[0][1] = 0.5 * (PrimVar_Grad_i[1][1] + PrimVar_Grad_i[2][0]);
-    S_ij[0][2] = 0.0;
-    S_ij[1][2] = 0.0;
-    S_ij[1][0] = S_ij[0][1];
-    S_ij[2][1] = S_ij[1][2];
-    S_ij[2][0] = S_ij[0][2];
-
-  }
-}
-
 void CSourcePieceWise_TurbSST::SetReynoldsStressMatrix(su2double turb_ke){
-  unsigned short iDim, jDim;
-  su2double **S_ij = new su2double* [3];
-  su2double divVel = 0;
-  su2double TWO3 = 2.0/3.0;
-
-
-
-  for (iDim = 0; iDim < 3; iDim++){
-    S_ij[iDim] = new su2double [3];
-  }
-
-  GetMeanRateOfStrainMatrix(S_ij);
-
-    /* --- Using rate of strain matrix, calculate Reynolds stress tensor --- */
-
-  for (iDim = 0; iDim < 3; iDim++){
-    divVel += S_ij[iDim][iDim];
-  }
-
-  for (iDim = 0; iDim < 3; iDim++){
-    for (jDim = 0; jDim < 3; jDim++){
-      MeanReynoldsStress[iDim][jDim] = TWO3 * turb_ke * delta3[iDim][jDim]
-      - Eddy_Viscosity_i / Density_i * (2 * S_ij[iDim][jDim] - TWO3 * divVel * delta3[iDim][jDim]);
+  ComputeStressTensor(nDim, MeanReynoldsStress, PrimVar_Grad_i+1, Eddy_Viscosity_i, Density_i, turb_ke, true);
+  for(unsigned short iDim=0; iDim<3; iDim++){
+    for(unsigned short jDim=0; jDim<3; jDim++){
+      MeanReynoldsStress[iDim][jDim] /= (-Density_i);
     }
   }
-
-  for (iDim = 0; iDim < 3; iDim++)
-    delete [] S_ij[iDim];
-  delete [] S_ij;
 }
 
 void CSourcePieceWise_TurbSST::SetPerturbedRSM(su2double turb_ke, const CConfig* config){
diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp
index 672073b26873..ec823ff9d19b 100644
--- a/SU2_CFD/src/python_wrapper_structure.cpp
+++ b/SU2_CFD/src/python_wrapper_structure.cpp
@@ -356,16 +356,9 @@ bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex)
   /*--- Parameters for the calculations ---*/
   // Pn: Pressure
   // Pinf: Pressure_infinite
-  // div_vel: Velocity divergence
-  // Dij: Dirac delta
-  su2double Pn = 0.0, div_vel = 0.0, Dij = 0.0;
+  su2double Pn = 0.0;
   su2double Viscosity = 0.0;
-  su2double Grad_Vel[3][3] = { {0.0, 0.0, 0.0} ,
-              {0.0, 0.0, 0.0} ,
-              {0.0, 0.0, 0.0} } ;
-  su2double Tau[3][3] = { {0.0, 0.0, 0.0} ,
-              {0.0, 0.0, 0.0} ,
-              {0.0, 0.0, 0.0} } ;
+  su2double Tau[3][3] = {{0.0}};
 
   su2double Pinf = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetPressure_Inf();
 
@@ -384,11 +377,6 @@ bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex)
     /*--- Get the values of pressure and viscosity ---*/
     Pn = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetPressure(iPoint);
     if (viscous_flow) {
-      for(iDim=0; iDim<nDim; iDim++) {
-        for(jDim=0; jDim<nDim; jDim++) {
-          Grad_Vel[iDim][jDim] = solver_container[ZONE_0][INST_0][FinestMesh][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint, iDim+1, jDim);
-        }
-      }
       Viscosity = solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint);
     }
 
@@ -399,23 +387,18 @@ bool CDriver::ComputeVertexForces(unsigned short iMarker, unsigned long iVertex)
 
     /*--- Calculate the viscous (shear stress) part of tn in the fluid nodes (force units) ---*/
     if ((incompressible || compressible) && viscous_flow) {
-      div_vel = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++)
-        div_vel += Grad_Vel[iDim][iDim];
-     if (incompressible) div_vel = 0.0;
-
+      CNumerics::ComputeStressTensor(nDim, Tau,
+        solver_container[ZONE_0][INST_0][FinestMesh][FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint)+1, Viscosity);
       for (iDim = 0; iDim < nDim; iDim++) {
-       for (jDim = 0 ; jDim < nDim; jDim++) {
-         Dij = 0.0; if (iDim == jDim) Dij = 1.0;
-         Tau[iDim][jDim] = Viscosity*(Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]) - TWO3*Viscosity*div_vel*Dij;
-         PyWrapNodalForce[iDim] += Tau[iDim][jDim]*Normal[jDim];
+        for (jDim = 0 ; jDim < nDim; jDim++) {
+          PyWrapNodalForce[iDim] += Tau[iDim][jDim]*Normal[jDim];
         }
       }
     }
 
     //Divide by local are in case of force density communication.
-   for(iDim = 0; iDim < nDim; iDim++) {
-     PyWrapNodalForceDensity[iDim] = PyWrapNodalForce[iDim]/Area;
+    for(iDim = 0; iDim < nDim; iDim++) {
+      PyWrapNodalForceDensity[iDim] = PyWrapNodalForce[iDim]/Area;
     }
 
     halo = false;
diff --git a/SU2_CFD/src/solvers/CAdjNSSolver.cpp b/SU2_CFD/src/solvers/CAdjNSSolver.cpp
index 7df523c9f76b..cea3942bc1fa 100644
--- a/SU2_CFD/src/solvers/CAdjNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CAdjNSSolver.cpp
@@ -604,7 +604,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
   normal_grad_psi5, normal_grad_T, sigma_partial, Laminar_Viscosity = 0.0, heat_flux_factor, temp_sens = 0.0,
   *Psi = nullptr, *U = nullptr, Enthalpy, gradPsi5_v, psi5_tau_partial, psi5_tau_grad_vel, source_v_1, Density,
   Pressure = 0.0, div_vel, val_turb_ke, vartheta, vartheta_partial, psi5_p_div_vel, Omega[3], rho_v[3] = {0.0,0.0,0.0},
-  CrossProduct[3], delta[3][3] = {{1.0, 0.0, 0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}}, r, ru, rv, rw, rE, p, T, dp_dr, dp_dru,
+  CrossProduct[3], r, ru, rv, rw, rE, p, T, dp_dr, dp_dru,
   dp_drv, dp_drw, dp_drE, dH_dr, dH_dru, dH_drv, dH_drw, dH_drE, H, D[3][3], Dd[3], Mach_Inf, eps, scale = 1.0,
   RefVel2, RefDensity, Mach2Vel, *Velocity_Inf, factor;
   const su2double* const* GridVel_Grad;
@@ -783,17 +783,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
             else
               val_turb_ke = 0.0;
 
-            div_vel = 0.0;
-            for (iDim = 0 ; iDim < nDim; iDim++) {
-              Velocity[iDim] = U[iDim+1]/Density;
-              div_vel += PrimVar_Grad[iDim+1][iDim];
-            }
-
-            for (iDim = 0 ; iDim < nDim; iDim++)
-              for (jDim = 0 ; jDim < nDim; jDim++)
-                tau[iDim][jDim] = Laminar_Viscosity*(PrimVar_Grad[jDim+1][iDim] + PrimVar_Grad[iDim+1][jDim])
-                - TWO3*Laminar_Viscosity*div_vel*delta[iDim][jDim]
-                - TWO3*Density*val_turb_ke*delta[iDim][jDim];
+            CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke);
 
             /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/
 
@@ -811,6 +801,12 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
 
             /*--- Form Sigma_Psi5v ---*/
 
+            div_vel = 0.0;
+            for (iDim = 0 ; iDim < nDim; iDim++) {
+              Velocity[iDim] = U[iDim+1]/Density;
+              div_vel += PrimVar_Grad[iDim+1][iDim];
+            }
+
             gradPsi5_v = 0.0;
             for (iDim = 0; iDim < nDim; iDim++) {
               gradPsi5_v += PsiVar_Grad[nDim+1][iDim]*Velocity[iDim];
diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp
index c8677f097906..c8add4efc55d 100644
--- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp
+++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp
@@ -1019,7 +1019,6 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry,
   su2double TauElem[3], TauTangent[3];
   su2double Tau[3][3];
   su2double TauNormal;
-  su2double div_vel=0, Delta;
 
   bool ionization = config->GetIonization();
 
@@ -1143,17 +1142,8 @@ void CNEMONSSolver::BC_Smoluchowski_Maxwell(CGeometry *geometry,
       for (iVar = 0; iVar < nVar; iVar ++)
         Res_Visc[iVar] = 0.0;
 
-      div_vel = 0.0;
-      for (iDim = 0; iDim < nDim; iDim++)
-        div_vel += Grad_PrimVar[VEL_INDEX+iDim][iDim];
-
+      CNumerics::ComputeStressTensor(nDim, Tau, Grad_PrimVar+VEL_INDEX, Viscosity);
       for (iDim = 0; iDim < nDim; iDim++) {
-        for (jDim = 0 ; jDim < nDim; jDim++) {
-          Delta = 0.0; if (iDim == jDim) Delta = 1.0;
-          Tau[iDim][jDim] = Viscosity*(Grad_PrimVar[VEL_INDEX+jDim][iDim] +
-              Grad_PrimVar[VEL_INDEX+iDim][jDim]  )
-              - TWO3*Viscosity*div_vel*Delta;
-        }
         TauElem[iDim] = 0.0;
         for (jDim = 0; jDim < nDim; jDim++)
           TauElem[iDim] += Tau[iDim][jDim]*UnitNormal[jDim];
diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp
index 418e84055b8c..f67e22a19ea6 100644
--- a/SU2_CFD/src/solvers/CNSSolver.cpp
+++ b/SU2_CFD/src/solvers/CNSSolver.cpp
@@ -471,23 +471,10 @@ void CNSSolver::AddDynamicGridResidualContribution(unsigned long iPoint, unsigne
   su2double eddy_viscosity = nodes->GetEddyViscosity(iPoint);
   su2double total_viscosity = laminar_viscosity + eddy_viscosity;
 
-  const auto Grad_Vel = &nodes->GetGradient_Primitive(iPoint)[1];
-
-  /*--- Divergence of the velocity ---*/
-
-  su2double div_vel = 0.0;
-  for (auto iDim = 0u; iDim < nDim; iDim++)
-    div_vel += Grad_Vel[iDim][iDim];
-
   /*--- Compute the viscous stress tensor ---*/
 
   su2double tau[MAXNDIM][MAXNDIM] = {{0.0}};
-  for (auto iDim = 0u; iDim < nDim; iDim++) {
-    for (auto jDim = 0u; jDim < nDim; jDim++) {
-      tau[iDim][jDim] = total_viscosity * (Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim]);
-    }
-    tau[iDim][iDim] -= TWO3*total_viscosity*div_vel;
-  }
+  CNumerics::ComputeStressTensor(nDim, tau, nodes->GetGradient_Primitive(iPoint)+1, total_viscosity);
 
   /*--- Dot product of the stress tensor with the grid velocity ---*/
 
@@ -988,21 +975,11 @@ void CNSSolver::SetTauWall_WF(CGeometry *geometry, CSolver **solver_container, C
       /*--- Compute the shear stress at the wall in the regular fashion
        by using the stress tensor on the surface ---*/
 
+      su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}, TauElem[MAXNDIM] = {0.0};
       su2double Lam_Visc_Wall = nodes->GetLaminarViscosity(iPoint);
+      CNumerics::ComputeStressTensor(nDim, tau, nodes->GetGradient_Primitive(iPoint)+1, Lam_Visc_Wall);
 
-      const auto GradVel = &nodes->GetGradient_Primitive(iPoint)[1];
-
-      su2double div_vel = 0.0;
-      for (auto iDim = 0u; iDim < nDim; iDim++)
-        div_vel += GradVel[iDim][iDim];
-
-      su2double tau[MAXNDIM][MAXNDIM] = {{0.0}}, TauElem[MAXNDIM] = {0.0};
       for (auto iDim = 0u; iDim < nDim; iDim++) {
-        for (auto jDim = 0u; jDim < nDim; jDim++) {
-          tau[iDim][jDim] = Lam_Visc_Wall * (GradVel[jDim][iDim] + GradVel[iDim][jDim]);
-        }
-        tau[iDim][iDim] -= TWO3*Lam_Visc_Wall*div_vel;
-
         TauElem[iDim] = GeometryToolbox::DotProduct(nDim, tau[iDim], UnitNormal);
       }
 
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index 8b4f08df584d..6743086dca2f 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -4142,11 +4142,7 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){
                        (config->GetKind_Solver() == DISC_ADJ_RANS));
 
   // Parameters for the calculations
-  su2double Pn = 0.0, div_vel = 0.0;
-  su2double Viscosity = 0.0;
-  su2double Tau[3][3] = {{0.0, 0.0, 0.0},{0.0, 0.0, 0.0},{0.0, 0.0, 0.0}};
-  su2double Grad_Vel[3][3] = {{0.0, 0.0, 0.0},{0.0, 0.0, 0.0},{0.0, 0.0, 0.0}};
-  su2double delta[3][3] = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{0.0, 0.0, 1.0}};
+  su2double Pn = 0.0;
   su2double auxForce[3] = {1.0, 0.0, 0.0};
 
   unsigned short iMarker;
@@ -4195,26 +4191,11 @@ void CSolver::ComputeVertexTractions(CGeometry *geometry, CConfig *config){
 
           // Calculate tn in the fluid nodes for the viscous term
           if (viscous_flow) {
-
-            Viscosity = base_nodes->GetLaminarViscosity(iPoint);
-
-            for (iDim = 0; iDim < nDim; iDim++) {
-              for (jDim = 0 ; jDim < nDim; jDim++) {
-                Grad_Vel[iDim][jDim] = base_nodes->GetGradient_Primitive(iPoint, iDim+1, jDim);
-              }
-            }
-
-            // Divergence of the velocity
-            div_vel = 0.0; for (iDim = 0; iDim < nDim; iDim++) div_vel += Grad_Vel[iDim][iDim];
-
+            su2double Viscosity = base_nodes->GetLaminarViscosity(iPoint);
+            su2double Tau[3][3];
+            CNumerics::ComputeStressTensor(nDim, Tau, base_nodes->GetGradient_Primitive(iPoint)+1, Viscosity);
             for (iDim = 0; iDim < nDim; iDim++) {
               for (jDim = 0 ; jDim < nDim; jDim++) {
-
-                // Viscous stress
-                Tau[iDim][jDim] = Viscosity*(Grad_Vel[jDim][iDim] + Grad_Vel[iDim][jDim])
-                                 - TWO3*Viscosity*div_vel*delta[iDim][jDim];
-
-                // Viscous component in the tn vector --> Units of force (non-dimensional).
                 auxForce[iDim] += Tau[iDim][jDim]*iNormal[jDim];
               }
             }
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
index 5f0ed7eed6be..054a7ba1eb48 100644
--- a/TestCases/hybrid_regression.py
+++ b/TestCases/hybrid_regression.py
@@ -232,7 +232,7 @@ def main():
     turb_naca0012_1c.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_1c.cfg_file  = "turb_NACA0012_uq_1c.cfg"
     turb_naca0012_1c.test_iter = 10
-    turb_naca0012_1c.test_vals = [-4.979389, 1.140070, 1.211965, 0.194237]
+    turb_naca0012_1c.test_vals = [-4.978913, 1.140283, 1.211887, 0.194208]
     test_list.append(turb_naca0012_1c)
 
     # NACA0012 2c
@@ -240,7 +240,7 @@ def main():
     turb_naca0012_2c.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_2c.cfg_file  = "turb_NACA0012_uq_2c.cfg"
     turb_naca0012_2c.test_iter = 10
-    turb_naca0012_2c.test_vals = [-5.484195, 0.969789, 1.310525, 0.231240]
+    turb_naca0012_2c.test_vals = [-5.484204, 0.969783, 1.311085, 0.231449]
     test_list.append(turb_naca0012_2c)
 
     # NACA0012 3c
@@ -256,7 +256,7 @@ def main():
     turb_naca0012_p1c1.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_p1c1.cfg_file  = "turb_NACA0012_uq_p1c1.cfg"
     turb_naca0012_p1c1.test_iter = 10
-    turb_naca0012_p1c1.test_vals = [-5.132081, 1.076462, 1.178093, 0.181595]
+    turb_naca0012_p1c1.test_vals = [-5.132991, 1.076082, 1.177974, 0.181556]
     test_list.append(turb_naca0012_p1c1)
 
     # NACA0012 p1c2
@@ -264,7 +264,7 @@ def main():
     turb_naca0012_p1c2.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_p1c2.cfg_file  = "turb_NACA0012_uq_p1c2.cfg"
     turb_naca0012_p1c2.test_iter = 10
-    turb_naca0012_p1c2.test_vals = [-5.556648, 0.945129, 1.240986, 0.205071]
+    turb_naca0012_p1c2.test_vals = [-5.556581, 0.945167, 1.240884, 0.205031]
     test_list.append(turb_naca0012_p1c2)
 
     ######################################
@@ -527,7 +527,7 @@ def main():
     stat_fsi.cfg_dir   = "fea_fsi/stat_fsi"
     stat_fsi.cfg_file  = "config.cfg"
     stat_fsi.test_iter = 7
-    stat_fsi.test_vals = [-3.242834, -4.866608, 0.000000, 11.000000]
+    stat_fsi.test_vals = [-3.242851, -4.866383, 0.000000, 11.000000]
     stat_fsi.multizone = True
     test_list.append(stat_fsi)
 
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index d7063c7ef4be..08bcb40c6b62 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -225,7 +225,7 @@ def main():
     poiseuille_profile.cfg_dir   = "navierstokes/poiseuille"
     poiseuille_profile.cfg_file  = "profile_poiseuille.cfg"
     poiseuille_profile.test_iter = 10
-    poiseuille_profile.test_vals = [-12.492876, -7.672445, -0.000000, 2.085796]
+    poiseuille_profile.test_vals = [-12.492859, -7.672756, -0.000000, 2.085796]
     poiseuille_profile.su2_exec  = "parallel_computation.py -f"
     poiseuille_profile.timeout   = 1600
     poiseuille_profile.tol       = 0.00001
@@ -683,7 +683,7 @@ def main():
     turb_naca0012_1c.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_1c.cfg_file  = "turb_NACA0012_uq_1c.cfg"
     turb_naca0012_1c.test_iter = 10
-    turb_naca0012_1c.test_vals = [-4.973124, 1.141759, 0.861168, 0.014208]
+    turb_naca0012_1c.test_vals = [-4.973120, 1.141760, 0.861180, 0.014211]
     turb_naca0012_1c.su2_exec  = "parallel_computation.py -f"
     turb_naca0012_1c.timeout   = 1600
     turb_naca0012_1c.tol       = 0.00001
@@ -694,7 +694,7 @@ def main():
     turb_naca0012_2c.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_2c.cfg_file  = "turb_NACA0012_uq_2c.cfg"
     turb_naca0012_2c.test_iter = 10
-    turb_naca0012_2c.test_vals = [-5.484227, 0.967174, 0.901039, 0.029201]
+    turb_naca0012_2c.test_vals = [-5.484231, 0.967172, 0.901005, 0.029186]
     turb_naca0012_2c.su2_exec  = "parallel_computation.py -f"
     turb_naca0012_2c.timeout   = 1600
     turb_naca0012_2c.tol       = 0.00001
@@ -716,7 +716,7 @@ def main():
     turb_naca0012_p1c1.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_p1c1.cfg_file  = "turb_NACA0012_uq_p1c1.cfg"
     turb_naca0012_p1c1.test_iter = 10
-    turb_naca0012_p1c1.test_vals = [-5.128931, 1.076207, 0.954587, 0.040176]
+    turb_naca0012_p1c1.test_vals = [-5.128763, 1.076245, 0.954666, 0.040210]
     turb_naca0012_p1c1.su2_exec  = "parallel_computation.py -f"
     turb_naca0012_p1c1.timeout   = 1600
     turb_naca0012_p1c1.tol       = 0.00001
@@ -727,7 +727,7 @@ def main():
     turb_naca0012_p1c2.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_p1c2.cfg_file  = "turb_NACA0012_uq_p1c2.cfg"
     turb_naca0012_p1c2.test_iter = 10
-    turb_naca0012_p1c2.test_vals = [-5.556651, 0.941872, 0.950627, 0.041946]
+    turb_naca0012_p1c2.test_vals = [-5.556632, 0.941881, 0.950638, 0.041949]
     turb_naca0012_p1c2.su2_exec  = "parallel_computation.py -f"
     turb_naca0012_p1c2.timeout   = 1600
     turb_naca0012_p1c2.tol       = 0.00001
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index f5e5d62f60ea..7a15a21d328e 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -250,7 +250,7 @@ def main():
     poiseuille_profile.cfg_dir   = "navierstokes/poiseuille"
     poiseuille_profile.cfg_file  = "profile_poiseuille.cfg"
     poiseuille_profile.test_iter = 10
-    poiseuille_profile.test_vals = [-12.494705, -7.711759, -0.000000, 2.085796] #last 4 columns
+    poiseuille_profile.test_vals = [-12.494720, -7.711373, -0.000000, 2.085796] #last 4 columns
     poiseuille_profile.su2_exec  = "SU2_CFD"
     poiseuille_profile.new_output = True
     poiseuille_profile.timeout   = 1600
@@ -806,7 +806,7 @@ def main():
     turb_naca0012_1c.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_1c.cfg_file  = "turb_NACA0012_uq_1c.cfg"
     turb_naca0012_1c.test_iter = 10
-    turb_naca0012_1c.test_vals = [-4.978069, 1.139128, 0.807062, 0.064726] #last 4 columns
+    turb_naca0012_1c.test_vals = [-4.978068, 1.139123, 0.806952, 0.064685] #last 4 columns
     turb_naca0012_1c.su2_exec  = "SU2_CFD"
     turb_naca0012_1c.new_output = True
     turb_naca0012_1c.timeout   = 1600
@@ -818,7 +818,7 @@ def main():
     turb_naca0012_2c.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_2c.cfg_file  = "turb_NACA0012_uq_2c.cfg"
     turb_naca0012_2c.test_iter = 10
-    turb_naca0012_2c.test_vals = [-5.484279, 0.967025, 0.823201, 0.070866] #last 4 columns
+    turb_naca0012_2c.test_vals = [-5.484282, 0.967023, 0.823129, 0.070840] #last 4 columns
     turb_naca0012_2c.su2_exec  = "SU2_CFD"
     turb_naca0012_2c.new_output = True
     turb_naca0012_2c.timeout   = 1600
@@ -842,7 +842,7 @@ def main():
     turb_naca0012_p1c1.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_p1c1.cfg_file  = "turb_NACA0012_uq_p1c1.cfg"
     turb_naca0012_p1c1.test_iter = 10
-    turb_naca0012_p1c1.test_vals = [-5.126796, 1.076577, 0.783116, 0.055988] #last 4 columns
+    turb_naca0012_p1c1.test_vals = [ -5.126540, 1.076620, 0.783153, 0.056001] #last 4 columns
     turb_naca0012_p1c1.su2_exec  = "SU2_CFD"
     turb_naca0012_p1c1.new_output = True
     turb_naca0012_p1c1.timeout   = 1600
@@ -854,7 +854,7 @@ def main():
     turb_naca0012_p1c2.cfg_dir   = "rans_uq/naca0012"
     turb_naca0012_p1c2.cfg_file  = "turb_NACA0012_uq_p1c2.cfg"
     turb_naca0012_p1c2.test_iter = 10
-    turb_naca0012_p1c2.test_vals = [-5.556585, 0.941677, 0.796006, 0.060817] #last 4 columns
+    turb_naca0012_p1c2.test_vals = [-5.556554, 0.941694, 0.795964, 0.060801] #last 4 columns
     turb_naca0012_p1c2.su2_exec  = "SU2_CFD"
     turb_naca0012_p1c2.new_output = True
     turb_naca0012_p1c2.timeout   = 1600