allow UQ + vectorization

Author: pcarruscag

SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp

@@ -36,13 +36,13 @@
  * \brief Generic factory implementation.
  */
 template<class ViscousDecorator>
-CNumericsSIMD* createNumerics(const CConfig& config, int iMesh) {
+CNumericsSIMD* createNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) {
   CNumericsSIMD* obj = nullptr;
   switch (config.GetKind_ConvNumScheme_Flow()) {
     case SPACE_UPWIND:
       switch (config.GetKind_Upwind_Flow()) {
         case ROE:
-          obj = new CRoeScheme<ViscousDecorator>(config, iMesh);
+          obj = new CRoeScheme<ViscousDecorator>(config, iMesh, turbVars);
           break;
       }
       break;
@@ -52,16 +52,16 @@ CNumericsSIMD* createNumerics(const CConfig& config, int iMesh) {
         case NO_CENTERED:
           break;
         case LAX:
-          obj = new CLaxScheme<ViscousDecorator>(config, iMesh);
+          obj = new CLaxScheme<ViscousDecorator>(config, iMesh, turbVars);
           break;
         case JST:
-          obj = new CJSTScheme<ViscousDecorator>(config, iMesh);
+          obj = new CJSTScheme<ViscousDecorator>(config, iMesh, turbVars);
           break;
         case JST_KE:
-          obj = new CJSTkeScheme<ViscousDecorator>(config, iMesh);
+          obj = new CJSTkeScheme<ViscousDecorator>(config, iMesh, turbVars);
           break;
         case JST_MAT:
-          obj = new CJSTmatScheme<ViscousDecorator>(config, iMesh);
+          obj = new CJSTmatScheme<ViscousDecorator>(config, iMesh, turbVars);
           break;
       }
       break;
@@ -74,17 +74,17 @@ CNumericsSIMD* createNumerics(const CConfig& config, int iMesh) {
  * createNumerics, which in turn instantiates the class templates of the different
  * numerical methods.
  */
-CNumericsSIMD* CNumericsSIMD::CreateNumerics(const CConfig& config, int nDim, int iMesh) {
+CNumericsSIMD* CNumericsSIMD::CreateNumerics(const CConfig& config, int nDim, int iMesh, const CVariable* turbVars) {
   if ((Double::Size < 4) && (SU2_MPI::GetRank() == MASTER_NODE)) {
     cout << "WARNING: SU2 was not compiled for an AVX-capable architecture." << endl;
   }
   CNumericsSIMD* obj = nullptr;
   if (config.GetViscous()) {
-    if (nDim == 2) obj = createNumerics<CCompressibleViscousFlux<2> >(config, iMesh);
-    if (nDim == 3) obj = createNumerics<CCompressibleViscousFlux<3> >(config, iMesh);
+    if (nDim == 2) obj = createNumerics<CCompressibleViscousFlux<2> >(config, iMesh, turbVars);
+    if (nDim == 3) obj = createNumerics<CCompressibleViscousFlux<3> >(config, iMesh, turbVars);
   } else {
-    if (nDim == 2) obj = createNumerics<CNoViscousFlux<2> >(config, iMesh);
-    if (nDim == 3) obj = createNumerics<CNoViscousFlux<3> >(config, iMesh);
+    if (nDim == 2) obj = createNumerics<CNoViscousFlux<2> >(config, iMesh, turbVars);
+    if (nDim == 3) obj = createNumerics<CNoViscousFlux<3> >(config, iMesh, turbVars);
   }
   return obj;
 }

SU2_CFD/include/numerics_simd/CNumericsSIMD.hpp

@@ -93,7 +93,8 @@ class CNumericsSIMD {
    * \param[in] config - Problem definitions.
    * \param[in] nDim - 2D or 3D.
    * \param[in] iMesh - Grid index.
+   * \param[in] turbVars - Turbulence variables.
    */
-  static CNumericsSIMD* CreateNumerics(const CConfig& config, int nDim, int iMesh);
+  static CNumericsSIMD* CreateNumerics(const CConfig& config, int nDim, int iMesh, const CVariable* turbVars = nullptr);
 
 };

SU2_CFD/include/numerics_simd/flow/diffusion/common.hpp

@@ -30,6 +30,7 @@
 #include "../../CNumericsSIMD.hpp"
 #include "../../util.hpp"
 #include "../variables.hpp"
+#include "../../../numerics/CNumerics.hpp"
 
 /*!
  * \brief Average gradients at i/j points.
@@ -65,14 +66,12 @@ FORCEINLINE void correctGradient(const PrimitiveType& V,
 }
 
 /*!
- * \brief Compute the stress tensor (using the total viscosity).
+ * \brief Compute the stress tensor.
  * \note Second viscosity term ignored.
  */
-template<size_t nVar, size_t nDim, class PrimitiveType>
-FORCEINLINE MatrixDbl<nDim> stressTensor(const PrimitiveType& V,
+template<size_t nVar, size_t nDim>
+FORCEINLINE MatrixDbl<nDim> stressTensor(Double viscosity,
                                          const MatrixDbl<nVar,nDim> grad) {
-  Double viscosity = V.laminarVisc() + V.eddyVisc();
-
   /*--- Hydrostatic term. ---*/
   Double velDiv = 0.0;
   for (size_t iDim = 0; iDim < nDim; ++iDim) {
@@ -91,6 +90,33 @@ FORCEINLINE MatrixDbl<nDim> stressTensor(const PrimitiveType& V,
   return tau;
 }
 
+/*!
+ * \brief Add perturbed stress tensor.
+ * \note Not inlined because it is not easy to vectorize properly, due to tred2 and tql2.
+ */
+template<class PrimitiveType, class MatrixType, size_t nDim, class... Ts>
+NEVERINLINE void addPerturbedRSM(const PrimitiveType& V,
+                                 const MatrixType& grad,
+                                 const Double& turb_ke,
+                                 MatrixDbl<nDim,nDim>& tau,
+                                 Ts... uq_args) {
+  /*--- Handle SIMD dimensions 1 by 1. ---*/
+  for (size_t k = 0; k < Double::Size; ++k) {
+    su2double velgrad[nDim][nDim];
+    for (size_t iVar = 0; iVar < nDim; ++iVar)
+      for (size_t iDim = 0; iDim < nDim; ++iDim)
+        velgrad[iVar][iDim] = grad(iVar+1,iDim)[k];
+
+    su2double rsm[3][3];
+    CNumerics::ComputePerturbedRSM(nDim, uq_args..., velgrad, V.density()[k],
+                                   V.eddyVisc()[k], turb_ke[k], rsm);
+
+    for (size_t iDim = 0; iDim < nDim; ++iDim)
+      for (size_t jDim = 0; jDim < nDim; ++jDim)
+        tau(iDim,jDim)[k] -= V.density()[k] * rsm[iDim][jDim];
+  }
+}
+
 /*!
  * \brief SA-QCR2000 modification of the stress tensor.
  */

SU2_CFD/include/numerics_simd/flow/diffusion/viscous_fluxes.hpp

@@ -79,19 +79,33 @@ class CCompressibleViscousFlux : public CNumericsSIMD {
   const su2double cp;
   const bool correct;
   const bool useSA_QCR;
+  const bool uq;
+  const bool uq_permute;
+  const size_t uq_eigval_comp;
+  const su2double uq_delta_b;
+  const su2double uq_urlx;
+
+  const CVariable* turbVars;
 
   /*!
    * \brief Constructor, initialize constants and booleans.
    */
   template<class... Ts>
-  CCompressibleViscousFlux(const CConfig& config, int iMesh, Ts&...) :
+  CCompressibleViscousFlux(const CConfig& config, int iMesh,
+                           const CVariable* turbVars_ = nullptr, Ts&...) :
     gamma(config.GetGamma()),
     gasConst(config.GetGas_ConstantND()),
     prandtlLam(config.GetPrandtl_Lam()),
     prandtlTurb(config.GetPrandtl_Turb()),
     cp(gamma * gasConst / (gamma - 1)),
     correct(iMesh == MESH_0),
-    useSA_QCR(config.GetQCR()) {
+    useSA_QCR(config.GetQCR()),
+    uq(config.GetUsing_UQ()),
+    uq_permute(config.GetUQ_Permute()),
+    uq_eigval_comp(config.GetEig_Val_Comp()),
+    uq_delta_b(config.GetUQ_Delta_B()),
+    uq_urlx(config.GetUQ_URLX()),
+    turbVars(turbVars_) {
   }
 
   /*!
@@ -130,12 +144,16 @@ class CCompressibleViscousFlux : public CNumericsSIMD {
     auto avgGrad = averageGradient<nPrimVarGrad,nDim>(iPoint, jPoint, gradient);
     if(correct) correctGradient(V, vector_ij, dist2_ij, avgGrad);
 
-    /// TODO: Uncertainty quantification (needs a way to access tke, maybe in ctor).
-
     /*--- Stress and heat flux tensors. ---*/
 
-    auto tau = stressTensor(avgV, avgGrad);
+    auto tau = stressTensor(avgV.laminarVisc() + (uq? Double(0.0) : avgV.eddyVisc()), avgGrad);
     if(useSA_QCR) addQCR(avgGrad, tau);
+    if(uq) {
+      Double turb_ke = 0.5*(gatherVariables(iPoint, turbVars->GetSolution()) +
+                            gatherVariables(jPoint, turbVars->GetSolution()));
+      addPerturbedRSM(avgV, avgGrad, turb_ke, tau,
+                      uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx);
+    }
 
     Double cond = cp * (avgV.laminarVisc()/prandtlLam + avgV.eddyVisc()/prandtlTurb);
     VectorDbl<nDim> heatFlux;

SU2_CFD/include/solvers/CEulerSolver.hpp

@@ -278,6 +278,13 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, COMPRESSIBLE> {
    */
   void SetCoefficient_Gradients(CConfig *config) const;
 
+  /*!
+   * \brief Instantiate a SIMD numerics object.
+   * \param[in] solvers - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void InstantiateEdgeNumerics(const CSolver* const* solvers, const CConfig* config) final;
+
 public:
   /*!
    * \brief Constructor of the class.

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -232,13 +232,18 @@ class CFVMFlowSolverBase : public CSolver {
   /*!
    * \brief Method to compute convective and viscous residual contribution using vectorized numerics.
    */
-  void EdgeFluxResidual(const CGeometry *geometry, const CConfig *config);
+  void EdgeFluxResidual(const CGeometry *geometry, const CSolver* const* solvers, const CConfig *config);
 
   /*!
    * \brief Sum the edge fluxes for each cell to populate the residual vector, only used on coarse grids.
    */
   void SumEdgeFluxes(const CGeometry* geometry);
 
+  /*!
+   * \brief Instantiate a SIMD numerics object.
+   */
+  inline virtual void InstantiateEdgeNumerics(const CSolver* const* solvers, const CConfig* config) {}
+
   /*!
    * \brief Destructor.
    */

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1287,7 +1287,12 @@ void CFVMFlowSolverBase<V, R>::BC_Custom(CGeometry* geometry, CSolver** solver_c
 }
 
 template <class V, ENUM_REGIME R>
-void CFVMFlowSolverBase<V, R>::EdgeFluxResidual(const CGeometry *geometry, const CConfig *config) {
+void CFVMFlowSolverBase<V, R>::EdgeFluxResidual(const CGeometry *geometry,
+                                                const CSolver* const* solvers,
+                                                const CConfig *config) {
+  if (!edgeNumerics) {
+    InstantiateEdgeNumerics(solvers, config);
+  }
 
   /*--- Loop over edge colors. ---*/
   for (auto color : EdgeColoring) {

SU2_CFD/include/solvers/CSolver.hpp

@@ -197,6 +197,10 @@ class CSolver {
     assert(base_nodes!=nullptr && "CSolver::base_nodes was not set properly, see brief for CSolver::SetBaseClassPointerToNodes()");
     return base_nodes;
   }
+  inline const CVariable* GetNodes() const {
+    assert(base_nodes!=nullptr && "CSolver::base_nodes was not set properly, see brief for CSolver::SetBaseClassPointerToNodes()");
+    return base_nodes;
+  }
 
   /*!
    * \brief Helper function to define the type and number of variables per point for each communication type.

SU2_CFD/include/variables/CVariable.hpp

@@ -447,7 +447,7 @@ class CVariable {
    * \brief Get the entire solution of the problem.
    * \return Reference to the solution matrix.
    */
-  inline const MatrixType& GetSolution(void) { return Solution; }
+  inline const MatrixType& GetSolution(void) const { return Solution; }
 
   /*!
    * \brief Get the solution of the problem.

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -349,26 +349,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
   /*--- Add the solver name (max 8 characters). ---*/
   SolverName = "C.FLOW";
 
-  /*--- Vectorized numerics. ---*/
-  if (config->GetUseVectorization()) {
-    const bool uncertain = config->GetUsing_UQ();
-    const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) ||
-                           (config->GetKind_FluidModel() == IDEAL_GAS);
-    const bool low_mach_corr = config->Low_Mach_Correction();
-
-    if (uncertain || !ideal_gas || low_mach_corr) {
-      SU2_MPI::Error("Some of the requested features are not yet "
-                     "supported with vectorization.", CURRENT_FUNCTION);
-    }
-
-    edgeNumerics = CNumericsSIMD::CreateNumerics(*config, nDim, iMesh);
-
-    if (!edgeNumerics) {
-      SU2_MPI::Error("The numerical scheme in use does not "
-                     "support vectorization.", CURRENT_FUNCTION);
-    }
-  }
-
   /*--- Finally, check that the static arrays will be large enough (keep this
    *    check at the bottom to make sure we consider the "final" values). ---*/
   if((nDim > MAXNDIM) || (nPrimVar > MAXNVAR) || (nSecondaryVar > MAXNVAR))
@@ -688,6 +668,31 @@ CEulerSolver::~CEulerSolver(void) {
 
 }
 
+void CEulerSolver::InstantiateEdgeNumerics(const CSolver* const* solver_container, const CConfig* config) {
+
+  SU2_OMP_MASTER {
+
+  const bool ideal_gas = (config->GetKind_FluidModel() == STANDARD_AIR) ||
+                         (config->GetKind_FluidModel() == IDEAL_GAS);
+  const bool low_mach_corr = config->Low_Mach_Correction();
+
+  if (!ideal_gas || low_mach_corr)
+    SU2_MPI::Error("Some of the requested features are not yet "
+                   "supported with vectorization.", CURRENT_FUNCTION);
+
+  if (solver_container[TURB_SOL])
+    edgeNumerics = CNumericsSIMD::CreateNumerics(*config, nDim, MGLevel, solver_container[TURB_SOL]->GetNodes());
+  else
+    edgeNumerics = CNumericsSIMD::CreateNumerics(*config, nDim, MGLevel);
+
+  if (!edgeNumerics)
+    SU2_MPI::Error("The numerical scheme in use does not "
+                   "support vectorization.", CURRENT_FUNCTION);
+
+  }
+  SU2_OMP_BARRIER
+}
+
 void CEulerSolver::InitTurboContainers(CGeometry *geometry, CConfig *config){
   unsigned short iMarker, iSpan, iVar;
   nSpanMax    = config->GetnSpanMaxAllZones();
@@ -2645,8 +2650,10 @@ void CEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container,
 void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container,
                                      CConfig *config, unsigned short iMesh, unsigned short iRKStep) {
 
-  /*--- If possible use the vectorized numerics instead. ---*/
-  if (edgeNumerics) { EdgeFluxResidual(geometry, config); return; }
+  if (config->GetUseVectorization()) {
+    EdgeFluxResidual(geometry, solver_container, config);
+    return;
+  }
 
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool jst_scheme = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0);
@@ -2732,8 +2739,10 @@ void CEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_conta
 void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_container,
                                    CNumerics **numerics_container, CConfig *config, unsigned short iMesh) {
 
-  /*--- If possible use the vectorized numerics instead. ---*/
-  if (edgeNumerics) { EdgeFluxResidual(geometry, config); return; }
+  if (config->GetUseVectorization()) {
+    EdgeFluxResidual(geometry, solver_container, config);
+    return;
+  }
 
   const auto InnerIter        = config->GetInnerIter();
   const bool implicit         = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);