all features supported with vectorized centered schemes

Author: pcarruscag

SU2_CFD/include/numerics_simd/CNumericsSIMD.cpp

@@ -32,43 +32,98 @@
 #include "flow/convection/centered.hpp"
 #include "flow/diffusion/viscous_fluxes.hpp"
 
+namespace {
+
 /*!
- * \brief Generic factory implementation.
+ * \brief Upwind factory implementation for ideal gas.
  */
 template<class ViscousDecorator>
+CNumericsSIMD* createUpwindIdealNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) {
+  CNumericsSIMD* obj = nullptr;
+  switch (config.GetKind_Upwind_Flow()) {
+    case ROE:
+      obj = new CRoeScheme<ViscousDecorator>(config, iMesh, turbVars);
+      break;
+  }
+  return obj;
+}
+
+/*!
+ * \brief Upwind factory implementation for real gas.
+ */
+template<class ViscousDecorator>
+CNumericsSIMD* createUpwindGeneralNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) {
+  return nullptr;
+}
+
+/*!
+ * \brief Centered factory implementation.
+ */
+template<class ViscousDecorator>
+CNumericsSIMD* createCenteredNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) {
+  CNumericsSIMD* obj = nullptr;
+  switch ((iMesh==MESH_0)? config.GetKind_Centered_Flow() : LAX) {
+    case NO_CENTERED:
+      break;
+    case LAX:
+      obj = new CLaxScheme<ViscousDecorator>(config, iMesh, turbVars);
+      break;
+    case JST:
+      obj = new CJSTScheme<ViscousDecorator>(config, iMesh, turbVars);
+      break;
+    case JST_KE:
+      obj = new CJSTkeScheme<ViscousDecorator>(config, iMesh, turbVars);
+      break;
+    case JST_MAT:
+      obj = new CJSTmatScheme<ViscousDecorator>(config, iMesh, turbVars);
+      break;
+  }
+  return obj;
+}
+
+/*!
+ * \brief Generic factory implementation.
+ */
+template<int nDim>
 CNumericsSIMD* createNumerics(const CConfig& config, int iMesh, const CVariable* turbVars) {
   CNumericsSIMD* obj = nullptr;
+  const bool ideal_gas = (config.GetKind_FluidModel() == STANDARD_AIR) ||
+                         (config.GetKind_FluidModel() == IDEAL_GAS);
+
   switch (config.GetKind_ConvNumScheme_Flow()) {
     case SPACE_UPWIND:
-      switch (config.GetKind_Upwind_Flow()) {
-        case ROE:
-          obj = new CRoeScheme<ViscousDecorator>(config, iMesh, turbVars);
-          break;
+      if (config.GetViscous()) {
+        if (ideal_gas)
+          obj = createUpwindIdealNumerics<CCompressibleViscousFlux<nDim> >(config, iMesh, turbVars);
+        else
+          obj = createUpwindGeneralNumerics<CGeneralCompressibleViscousFlux<nDim> >(config, iMesh, turbVars);
+      }
+      else {
+        if (ideal_gas)
+          obj = createUpwindIdealNumerics<CNoViscousFlux<nDim> >(config, iMesh, turbVars);
+        else
+          obj = createUpwindGeneralNumerics<CNoViscousFlux<nDim> >(config, iMesh, turbVars);
       }
       break;
 
     case SPACE_CENTERED:
-      switch ((iMesh==MESH_0)? config.GetKind_Centered_Flow() : LAX) {
-        case NO_CENTERED:
-          break;
-        case LAX:
-          obj = new CLaxScheme<ViscousDecorator>(config, iMesh, turbVars);
-          break;
-        case JST:
-          obj = new CJSTScheme<ViscousDecorator>(config, iMesh, turbVars);
-          break;
-        case JST_KE:
-          obj = new CJSTkeScheme<ViscousDecorator>(config, iMesh, turbVars);
-          break;
-        case JST_MAT:
-          obj = new CJSTmatScheme<ViscousDecorator>(config, iMesh, turbVars);
-          break;
+      if (config.GetViscous()) {
+        if (ideal_gas)
+          obj = createCenteredNumerics<CCompressibleViscousFlux<nDim> >(config, iMesh, turbVars);
+        else
+          obj = createCenteredNumerics<CGeneralCompressibleViscousFlux<nDim> >(config, iMesh, turbVars);
+      }
+      else {
+        obj = createCenteredNumerics<CNoViscousFlux<nDim> >(config, iMesh, turbVars);
       }
       break;
   }
+
   return obj;
 }
 
+} // namespace
+
 /*!
  * \brief This function instantiates both 2D and 3D versions of the implementation in
  * createNumerics, which in turn instantiates the class templates of the different
@@ -78,13 +133,8 @@ CNumericsSIMD* CNumericsSIMD::CreateNumerics(const CConfig& config, int nDim, in
   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, turbVars);
-    if (nDim == 3) obj = createNumerics<CCompressibleViscousFlux<3> >(config, iMesh, turbVars);
-  } else {
-    if (nDim == 2) obj = createNumerics<CNoViscousFlux<2> >(config, iMesh, turbVars);
-    if (nDim == 3) obj = createNumerics<CNoViscousFlux<3> >(config, iMesh, turbVars);
-  }
-  return obj;
+  if (nDim == 2) return createNumerics<2>(config, iMesh, turbVars);
+  if (nDim == 3) return createNumerics<3>(config, iMesh, turbVars);
+
+  return nullptr;
 }

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

@@ -62,14 +62,13 @@ class CNoViscousFlux : public CNumericsSIMD {
 };
 
 /*!
- * \class CCompressibleViscousFlux
- * \brief Decorator class to add viscous fluxes (compressible flow, ideal gas).
+ * \class CCompressibleViscousFluxBase
+ * \brief Decorator class to add viscous fluxes (compressible flow).
  */
-template<size_t NDIM>
-class CCompressibleViscousFlux : public CNumericsSIMD {
+template<size_t NDIM, class Derived>
+class CCompressibleViscousFluxBase : public CNumericsSIMD {
 protected:
   static constexpr size_t nDim = NDIM;
-  static constexpr size_t nPrimVar = nDim+7;
   static constexpr size_t nPrimVarGrad = nDim+1;
 
   const su2double gamma;
@@ -91,8 +90,8 @@ class CCompressibleViscousFlux : public CNumericsSIMD {
    * \brief Constructor, initialize constants and booleans.
    */
   template<class... Ts>
-  CCompressibleViscousFlux(const CConfig& config, int iMesh,
-                           const CVariable* turbVars_ = nullptr, Ts&...) :
+  CCompressibleViscousFluxBase(const CConfig& config, int iMesh,
+                               const CVariable* turbVars_ = nullptr, Ts&...) :
     gamma(config.GetGamma()),
     gasConst(config.GetGas_ConstantND()),
     prandtlLam(config.GetPrandtl_Lam()),
@@ -128,7 +127,11 @@ class CCompressibleViscousFlux : public CNumericsSIMD {
                                 MatrixDbl<nVar>& jac_i,
                                 MatrixDbl<nVar>& jac_j) const {
 
-    static_assert(PrimVarType::nVar <= nPrimVar,"");
+    static_assert(PrimVarType::nVar <= Derived::nPrimVar,"");
+
+    /*--- Pointer on which to call the "compile-time virtual" methods. ---*/
+
+    const auto derived = static_cast<const Derived*>(this);
 
     const auto& solution = static_cast<const CNSVariable&>(solution_);
     const auto& gradient = solution.GetGradient_Primitive();
@@ -155,7 +158,7 @@ class CCompressibleViscousFlux : public CNumericsSIMD {
                       uq_eigval_comp, uq_permute, uq_delta_b, uq_urlx);
     }
 
-    Double cond = cp * (avgV.laminarVisc()/prandtlLam + avgV.eddyVisc()/prandtlTurb);
+    Double cond = derived->thermalConductivity(avgV);
     VectorDbl<nDim> heatFlux;
     for (size_t iDim = 0; iDim < nDim; ++iDim) {
       heatFlux(iDim) = cond * avgGrad(0,iDim);
@@ -175,23 +178,9 @@ class CCompressibleViscousFlux : public CNumericsSIMD {
 
     Double dist_ij = sqrt(dist2_ij);
     auto dtau = stressTensorJacobian<nVar>(avgV, unitNormal, dist_ij);
-    Double contraction = 0.0;
-    for (size_t iDim = 0; iDim < nDim; ++iDim) {
-      contraction += dtau(iDim,0) * avgV.velocity(iDim);
-    }
 
     /*--- Energy flux Jacobian. ---*/
-    VectorDbl<nVar> dEdU;
-    Double vel2 = 0.5 * squaredNorm<nDim>(avgV.velocity());
-    Double phi = (gamma-1) / avgV.density();
-    Double RdTdrho = phi*vel2 - avgV.pressure() / pow(avgV.density(),2);
-    Double condOnRd = cond / (gasConst * dist_ij);
-
-    dEdU(0) = area * (condOnRd * RdTdrho - contraction);
-    for (size_t iDim = 0; iDim < nDim; ++iDim) {
-      dEdU(iDim+1) = area * (condOnRd*phi*avgV.velocity(iDim) + dtau(iDim,0));
-    }
-    dEdU(nDim+1) = area * condOnRd * phi;
+    auto dEdU = derived->energyJacobian(avgV, dtau, cond, area, dist_ij, iPoint, jPoint, solution);
 
     /*--- Update momentum and energy terms ("symmetric" part). ---*/
     for (size_t iDim = 0; iDim < nDim; ++iDim) {
@@ -252,3 +241,125 @@ class CCompressibleViscousFlux : public CNumericsSIMD {
     viscousTerms(iEdge, iPoint, jPoint, avgV, V, solution_, vector_ij, geometry, args...);
   }
 };
+
+/*!
+ * \class CCompressibleViscousFlux
+ * \brief Decorator class to add viscous fluxes (compressible flow, ideal gas).
+ */
+template<size_t NDIM>
+class CCompressibleViscousFlux : public CCompressibleViscousFluxBase<NDIM, CCompressibleViscousFlux<NDIM> > {
+public:
+  static constexpr size_t nPrimVar = NDIM+7;
+  using Base = CCompressibleViscousFluxBase<NDIM, CCompressibleViscousFlux<NDIM> >;
+  using Base::gamma;
+  using Base::gasConst;
+  using Base::prandtlLam;
+  using Base::prandtlTurb;
+  using Base::cp;
+
+  /*!
+   * \brief Constructor, initialize constants and booleans.
+   */
+  template<class... Ts>
+  CCompressibleViscousFlux(Ts&... args) : Base(args...) {}
+
+  /*!
+   * \brief Compute the thermal conductivity.
+   */
+  template<class PrimitiveType>
+  FORCEINLINE Double thermalConductivity(const PrimitiveType& V) const {
+    return cp * (V.laminarVisc()/prandtlLam + V.eddyVisc()/prandtlTurb);
+  }
+
+  /*!
+   * \brief Compute Jacobian of the energy flux, except the part due to the work of viscous forces.
+   */
+  template<size_t nVar, size_t nDim, class PrimitiveType, class... Ts>
+  FORCEINLINE VectorDbl<nVar> energyJacobian(const PrimitiveType& V,
+                                             const MatrixDbl<nDim,nVar>& dtau,
+                                             Double thermalCond,
+                                             Double area,
+                                             Double dist_ij,
+                                             Ts&... args) const {
+    Double vel2 = 0.5 * squaredNorm<nDim>(V.velocity());
+    Double phi = (gamma-1) / V.density();
+    Double RdTdrho = phi*vel2 - V.pressure() / pow(V.density(),2);
+    Double condOnRd = thermalCond / (gasConst * dist_ij);
+    Double contraction = 0.0;
+    for (size_t iDim = 0; iDim < nDim; ++iDim) {
+      contraction += dtau(iDim,0) * V.velocity(iDim);
+    }
+    VectorDbl<nVar> dEdU;
+    dEdU(0) = area * (condOnRd * RdTdrho - contraction);
+    for (size_t iDim = 0; iDim < nDim; ++iDim) {
+      dEdU(iDim+1) = area * (dtau(iDim,0) - condOnRd*phi*V.velocity(iDim));
+    }
+    dEdU(nDim+1) = area * condOnRd * phi;
+
+    return dEdU;
+  }
+};
+
+/*!
+ * \class CGeneralCompressibleViscousFlux
+ * \brief Decorator class to add viscous fluxes (compressible flow, real gas).
+ */
+template<size_t NDIM>
+class CGeneralCompressibleViscousFlux : public CCompressibleViscousFluxBase<NDIM, CGeneralCompressibleViscousFlux<NDIM> > {
+public:
+  static constexpr size_t nPrimVar = NDIM+9;
+  static constexpr size_t nSecVar = 4;
+  using Base = CCompressibleViscousFluxBase<NDIM, CGeneralCompressibleViscousFlux<NDIM> >;
+  using Base::prandtlTurb;
+
+  /*!
+   * \brief Constructor, initialize constants and booleans.
+   */
+  template<class... Ts>
+  CGeneralCompressibleViscousFlux(Ts&... args) : Base(args...) {}
+
+  /*!
+   * \brief Compute the thermal conductivity.
+   */
+  template<class PrimitiveType>
+  FORCEINLINE Double thermalConductivity(const PrimitiveType& V) const {
+    return V.thermalCond() + V.cp()*V.eddyVisc()/prandtlTurb;
+  }
+
+  /*!
+   * \brief Compute Jacobian of the energy flux, except the part due to the work of viscous forces.
+   */
+  template<size_t nVar, size_t nDim, class PrimitiveType, class VariableType>
+  FORCEINLINE VectorDbl<nVar> energyJacobian(const PrimitiveType& V,
+                                             const MatrixDbl<nDim,nVar>& dtau,
+                                             Double thermalCond,
+                                             Double area,
+                                             Double dist_ij,
+                                             Int iPoint,
+                                             Int jPoint,
+                                             const VariableType& solution) const {
+    Double vel2 = squaredNorm<nDim>(V.velocity());
+    Double contraction = 0.0;
+    for (size_t iDim = 0; iDim < nDim; ++iDim) {
+      contraction += dtau(iDim,0) * V.velocity(iDim);
+    }
+
+    auto secVar_i = gatherVariables<nSecVar>(iPoint, solution.GetSecondary());
+    auto secVar_j = gatherVariables<nSecVar>(jPoint, solution.GetSecondary());
+
+    Double dTdrho_e = 0.5 * (secVar_i(2) + secVar_j(2));
+    Double dTde_rho = 0.5 * (secVar_i(3) + secVar_j(3)) / V.density();
+
+    Double condOnDist = thermalCond / dist_ij;
+    Double dTdrho = dTdrho_e + dTde_rho*(vel2-V.enthalpy()+V.pressure()/V.density());
+
+    VectorDbl<nVar> dEdU;
+    dEdU(0) = area * (condOnDist * dTdrho - contraction);
+    for (size_t iDim = 0; iDim < nDim; ++iDim) {
+      dEdU(iDim+1) = area * (dtau(iDim,0) - condOnDist*dTde_rho*V.velocity(iDim));
+    }
+    dEdU(nDim+1) = area * condOnDist * dTde_rho;
+
+    return dEdU;
+  }
+};

SU2_CFD/include/numerics_simd/flow/variables.hpp

@@ -50,13 +50,17 @@ struct CCompressiblePrimitives {
   FORCEINLINE const Double& velocity(size_t iDim) const { return all(iDim+1); }
   FORCEINLINE const Double* velocity() const { return &velocity(0); }
 
-  /*--- Un-reconstructed variables (not allocated by default). ---*/
+  /*--- Un-reconstructed variables. ---*/
   FORCEINLINE Double& speedSound() { return all(nDim+4); }
   FORCEINLINE Double& laminarVisc() { return all(nDim+5); }
   FORCEINLINE Double& eddyVisc() { return all(nDim+6); }
+  FORCEINLINE Double& thermalCond() { return all(nDim+7); }
+  FORCEINLINE Double& cp() { return all(nDim+8); }
   FORCEINLINE const Double& speedSound() const { return all(nDim+4); }
   FORCEINLINE const Double& laminarVisc() const { return all(nDim+5); }
   FORCEINLINE const Double& eddyVisc() const { return all(nDim+6); }
+  FORCEINLINE const Double& thermalCond() const { return all(nDim+7); }
+  FORCEINLINE const Double& cp() const { return all(nDim+8); }
 };
 
 /*!