Merge pull request #2356 from su2code/load_less_stuff

Recompute density and enthalpy instead of reconstructing

Author: pcarruscag

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

@@ -35,50 +35,52 @@
 /*!
  * \brief Unlimited reconstruction.
  */
-template<size_t nVar, size_t nDim, class Gradient_t>
+template<size_t nVarGrad_ = 0, size_t nVar, size_t nDim, class Gradient_t>
 FORCEINLINE void musclUnlimited(Int iPoint,
                                 const VectorDbl<nDim>& vector_ij,
                                 Double scale,
                                 const Gradient_t& gradient,
                                 VectorDbl<nVar>& vars) {
-  auto grad = gatherVariables<nVar,nDim>(iPoint, gradient);
-  for (size_t iVar = 0; iVar < nVar; ++iVar) {
+  constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : nVar;
+  auto grad = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
+  for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
     vars(iVar) += scale * dot(grad[iVar], vector_ij);
   }
 }
 
 /*!
  * \brief Limited reconstruction with point-based limiter.
  */
-template<size_t nVar, size_t nDim, class Limiter_t, class Gradient_t>
+template<size_t nVarGrad_ = 0, size_t nVar, size_t nDim, class Limiter_t, class Gradient_t>
 FORCEINLINE void musclPointLimited(Int iPoint,
                                    const VectorDbl<nDim>& vector_ij,
                                    Double scale,
                                    const Limiter_t& limiter,
                                    const Gradient_t& gradient,
                                    VectorDbl<nVar>& vars) {
-  auto lim = gatherVariables<nVar>(iPoint, limiter);
-  auto grad = gatherVariables<nVar,nDim>(iPoint, gradient);
-  for (size_t iVar = 0; iVar < nVar; ++iVar) {
+  constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : nVar;
+  auto lim = gatherVariables<nVarGrad>(iPoint, limiter);
+  auto grad = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
+  for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
     vars(iVar) += lim(iVar) * scale * dot(grad[iVar], vector_ij);
   }
 }
 
 /*!
  * \brief Limited reconstruction with edge-based limiter.
  */
-template<size_t nDim, class VarType, class Gradient_t>
+template<size_t nVarGrad_ = 0, size_t nDim, class VarType, class Gradient_t>
 FORCEINLINE void musclEdgeLimited(Int iPoint,
                                   Int jPoint,
                                   const VectorDbl<nDim>& vector_ij,
                                   const Gradient_t& gradient,
                                   CPair<VarType>& V) {
-  constexpr size_t nVar = VarType::nVar;
+  constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
 
-  auto grad_i = gatherVariables<nVar,nDim>(iPoint, gradient);
-  auto grad_j = gatherVariables<nVar,nDim>(jPoint, gradient);
+  auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
+  auto grad_j = gatherVariables<nVarGrad,nDim>(jPoint, gradient);
 
-  for (size_t iVar = 0; iVar < nVar; ++iVar) {
+  for (size_t iVar = 0; iVar < nVarGrad; ++iVar) {
     const Double proj_i = dot(grad_i[iVar], vector_ij);
     const Double proj_j = dot(grad_j[iVar], vector_ij);
     const Double delta_ij = V.j.all(iVar) - V.i.all(iVar);
@@ -93,15 +95,21 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
 
 /*!
  * \brief Retrieve primitive variables for points i/j, reconstructing them if needed.
+ * \note Density and enthalpy are recomputed from ideal gas EOS.
  * \param[in] iEdge, iPoint, jPoint - Edge and its nodes.
+ * \param[in] gamma - Heat capacity ratio.
+ * \param[in] gasConst - Specific gas constant.
  * \param[in] muscl - If true, reconstruct, else simply fetch.
+ * \param[in] limiterType - Type of flux limiter.
  * \param[in] V1st - Pair of compressible flow primitives for nodes i,j.
  * \param[in] vector_ij - Distance vector from i to j.
  * \param[in] solution - Entire solution container (a derived CVariable).
  * \return Pair of primitive variables.
  */
 template<class ReconVarType, class PrimVarType, size_t nDim, class VariableType>
 FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int jPoint,
+                                                      const su2double& gamma,
+                                                      const su2double& gasConst,
                                                       bool muscl, LIMITER limiterType,
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& vector_ij,
@@ -119,19 +127,29 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
   }
 
   if (muscl) {
+    /*--- Recompute density and enthalpy instead of reconstructing. ---*/
+    constexpr auto nVarGrad = ReconVarType::nVar - 2;
+
     switch (limiterType) {
     case LIMITER::NONE:
-      musclUnlimited(iPoint, vector_ij, 0.5, gradients, V.i.all);
-      musclUnlimited(jPoint, vector_ij,-0.5, gradients, V.j.all);
+      musclUnlimited<nVarGrad>(iPoint, vector_ij, 0.5, gradients, V.i.all);
+      musclUnlimited<nVarGrad>(jPoint, vector_ij,-0.5, gradients, V.j.all);
       break;
     case LIMITER::VAN_ALBADA_EDGE:
-      musclEdgeLimited(iPoint, jPoint, vector_ij, gradients, V);
+      musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V);
       break;
     default:
-      musclPointLimited(iPoint, vector_ij, 0.5, limiters, gradients, V.i.all);
-      musclPointLimited(jPoint, vector_ij,-0.5, limiters, gradients, V.j.all);
+      musclPointLimited<nVarGrad>(iPoint, vector_ij, 0.5, limiters, gradients, V.i.all);
+      musclPointLimited<nVarGrad>(jPoint, vector_ij,-0.5, limiters, gradients, V.j.all);
       break;
     }
+    V.i.density() = V.i.pressure() / (gasConst * V.i.temperature());
+    V.j.density() = V.j.pressure() / (gasConst * V.j.temperature());
+
+    const su2double cp = gasConst * gamma / (gamma - 1);
+    V.i.enthalpy() = cp * V.i.temperature() + 0.5 * squaredNorm<nDim>(V.i.velocity());
+    V.j.enthalpy() = cp * V.j.temperature() + 0.5 * squaredNorm<nDim>(V.j.velocity());
+
     /*--- Detect a non-physical reconstruction based on negative pressure or density. ---*/
     const Double neg_p_or_rho = fmax(fmin(V.i.pressure(), V.j.pressure()) < 0.0,
                                      fmin(V.i.density(), V.j.density()) < 0.0);

SU2_CFD/include/numerics_simd/flow/convection/roe.hpp

@@ -56,6 +56,7 @@ class CRoeBase : public Base {
 
   const su2double kappa;
   const su2double gamma;
+  const su2double gasConst;
   const su2double entropyFix;
   const bool finestGrid;
   const bool dynamicGrid;
@@ -69,6 +70,7 @@ class CRoeBase : public Base {
   CRoeBase(const CConfig& config, unsigned iMesh, Ts&... args) : Base(config, iMesh, args...),
     kappa(config.GetRoe_Kappa()),
     gamma(config.GetGamma()),
+    gasConst(config.GetGas_ConstantND()),
     entropyFix(config.GetEntropyFix_Coeff()),
     finestGrid(iMesh == MESH_0),
     dynamicGrid(config.GetDynamic_Grid()),
@@ -117,7 +119,7 @@ class CRoeBase : public Base {
     V1st.j.all = gatherVariables<nPrimVar>(jPoint, solution.GetPrimitive());
 
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-                 iEdge, iPoint, jPoint, muscl, typeLimiter, V1st, vector_ij, solution);
+        iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution);
 
     /*--- Compute conservative variables. ---*/
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1485,6 +1485,11 @@ void CFVMFlowSolverBase<V, R>::EdgeFluxResidual(const CGeometry *geometry,
                      "by the SIMD length (2, 4, or 8).", CURRENT_FUNCTION);
     }
     InstantiateEdgeNumerics(solvers, config);
+
+    /*--- The SIMD numerics do not use gradients of density and enthalpy. ---*/
+    if (!config->GetContinuous_Adjoint()) {
+      SU2_OMP_SAFE_GLOBAL_ACCESS(nPrimVarGrad = std::min<unsigned short>(nDim + 2, nPrimVarGrad);)
+    }
   }
 
   /*--- Non-physical counter. ---*/

SU2_CFD/include/variables/CEulerVariable.hpp

@@ -116,10 +116,8 @@ class CEulerVariable : public CFlowVariable {
    */
   bool SetSoundSpeed(unsigned long iPoint, su2double soundspeed2) final {
     if (soundspeed2 < 0.0) return true;
-    else {
-      Primitive(iPoint,nDim+4) = sqrt(soundspeed2);
-      return false;
-    }
+    Primitive(iPoint,nDim+4) = sqrt(soundspeed2);
+    return false;
   }
 
   /*!

SU2_CFD/src/output/CElasticityOutput.cpp

@@ -111,7 +111,6 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS
   /*--- Linear analysis: RMS of the displacements in the nDim coordinates ---*/
   /*--- Nonlinear analysis: UTOL, RTOL and DTOL (defined in the Postprocessing function) ---*/
 
-
   if (linear_analysis){
     SetHistoryOutputValue("RMS_DISP_X", log10(fea_solver->GetRes_RMS(0)));
     SetHistoryOutputValue("RMS_DISP_Y", log10(fea_solver->GetRes_RMS(1)));

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -64,6 +64,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
                          (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);
   const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING);
   const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint();
+  const bool centered = config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED;
 
   int Unst_RestartIter = 0;
   unsigned long iPoint, iMarker, counter_local = 0, counter_global = 0;
@@ -116,9 +117,12 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
 
   nDim = geometry->GetnDim();
 
-  nVar = nDim+2;
-  nPrimVar = nDim+9; nPrimVarGrad = nDim+4;
-  nSecondaryVar = nSecVar; nSecondaryVarGrad = 2;
+  nVar = nDim + 2;
+  nPrimVar = nDim + 9;
+  /*--- Centered schemes only need gradients for viscous fluxes (T and v). ---*/
+  nPrimVarGrad = nDim + (centered && !config->GetContinuous_Adjoint() ? 1 : 4);
+  nSecondaryVar = nSecVar;
+  nSecondaryVarGrad = 2;
 
   /*--- Initialize nVarGrad for deallocation ---*/
 

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -56,6 +56,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned
                     (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND));
   bool time_stepping = config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING;
   bool adjoint = (config->GetContinuous_Adjoint()) || (config->GetDiscrete_Adjoint());
+  const bool centered = config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED;
 
   /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */
   dynamic_grid = config->GetDynamic_Grid();
@@ -117,7 +118,10 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned
   nDim = geometry->GetnDim();
 
   /*--- Make sure to align the sizes with the constructor of CIncEulerVariable. ---*/
-  nVar = nDim+2; nPrimVar = nDim+9; nPrimVarGrad = nDim+4;
+  nVar = nDim + 2;
+  nPrimVar = nDim + 9;
+  /*--- Centered schemes only need gradients for viscous fluxes (T and v, but we need also to include P). ---*/
+  nPrimVarGrad = nDim + (centered ? 2 : 4);
 
   /*--- Initialize nVarGrad for deallocation ---*/
 

SU2_CFD/src/variables/CEulerVariable.cpp

@@ -30,7 +30,8 @@
 
 CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2double energy, unsigned long npoint,
                                unsigned long ndim, unsigned long nvar, const CConfig *config)
-  : CFlowVariable(npoint, ndim, nvar, ndim + 9, ndim + 4, config),
+  : CFlowVariable(npoint, ndim, nvar, ndim + 9,
+                  ndim + (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED && !config->GetContinuous_Adjoint() ? 1 : 4), config),
     indices(ndim, 0) {
 
   const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
@@ -77,7 +78,7 @@ CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2
     Grad_AuxVar.resize(nPoint, nAuxVar, nDim, 0.0);
     AuxVar.resize(nPoint, nAuxVar) = su2double(0.0);
   }
-  
+
   if (config->GetKind_FluidModel() == ENUM_FLUIDMODEL::DATADRIVEN_FLUID){
     DataDrivenFluid = true;
     DatasetExtrapolation.resize(nPoint) = 0;

SU2_CFD/src/variables/CIncEulerVariable.cpp

@@ -30,7 +30,8 @@
 
 CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *velocity, su2double temperature,
                                      unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config)
-  : CFlowVariable(npoint, ndim, nvar, ndim + 9, ndim + 4, config),
+  : CFlowVariable(npoint, ndim, nvar, ndim + 9,
+                  ndim + (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED ? 2 : 4), config),
     indices(ndim, 0) {
 
   const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||

TestCases/TestCase.py

@@ -464,6 +464,9 @@ def run_filediff(self, with_tsan=False, with_asan=False):
                 print('Ignored entries:     ' + str(ignore_counter))
                 print('Maximum difference:  ' + str(max_delta) + '%')
 
+            if not passed:
+                print(open(self.test_file).readlines())
+
         print('==================== End Test: %s ====================\n'%self.tag)
 
         sys.stdout.flush()