Recompute density and enthalpy instead of reconstructing

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()

TestCases/cont_adj_euler/wedge/of_grad_combo.dat.ref

@@ -1,5 +1,5 @@
 VARIABLES="VARIABLE"      , "GRADIENT"      , "FINDIFF_STEP"  
-          0               , 0.00767644      , 0.0001          
-          1               , 0.00498358      , 0.0001          
-          2               , 0.00246134      , 0.0001          
-          3               , 0.000893054     , 0.0001          
+          0               , 0.00770904      , 0.0001          
+          1               , 0.00500468      , 0.0001          
+          2               , 0.00247269      , 0.0001          
+          3               , 0.000899035     , 0.0001          

TestCases/disc_adj_fsi/config.cfg

@@ -1,5 +1,5 @@
 SOLVER= MULTIPHYSICS
-MATH_PROBLEM= DISCRETE_ADJOINT
+RESTART_SOL= NO
 CONFIG_LIST=(configFlow.cfg, configFEA.cfg)
 
 MARKER_ZONE_INTERFACE = (UpperWall, UpperWallS, LowerWall, LowerWallS)

TestCases/disc_adj_fsi/configFEA.cfg

@@ -4,15 +4,11 @@
 % Author: R.Sanchez                                                            %
 % Institution: Imperial College London                                         %
 % Date: 2017.11.29                                                             %
-% File Version 8.1.0 "Harrier"                                                    %
+% File Version 8.1.0 "Harrier"                                                 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 SOLVER= ELASTICITY
 
-MATH_PROBLEM= DISCRETE_ADJOINT
-
-RESTART_SOL= NO
-
 PRESTRETCH = YES
 PRESTRETCH_FILENAME = prestretch.dat
 
@@ -27,8 +23,6 @@ REFERENCE_GEOMETRY_FORMAT = SU2
 % Consider only the surface
 REFERENCE_GEOMETRY_SURFACE = NO
 
-READ_BINARY_RESTART=NO
-
 STAT_RELAX_PARAMETER= 1.0
 BGS_RELAXATION = FIXED_PARAMETER
 PREDICTOR_ORDER=0
@@ -49,13 +43,13 @@ MARKER_CLAMPED = ( Clamped_Right, Clamped_Left )
 MARKER_FLUID_LOAD= ( LowerWallS, UpperWallS)
 
 
-LINEAR_SOLVER= CONJUGATE_GRADIENT
-LINEAR_SOLVER_PREC= JACOBI
+LINEAR_SOLVER= FGMRES
+LINEAR_SOLVER_PREC= ILU
 LINEAR_SOLVER_ERROR= 1E-9
-LINEAR_SOLVER_ITER= 50000
+LINEAR_SOLVER_ITER= 100
 
-DISCADJ_LIN_SOLVER = CONJUGATE_GRADIENT
-DISCADJ_LIN_PREC = JACOBI
+DISCADJ_LIN_SOLVER = FGMRES
+DISCADJ_LIN_PREC = ILU
 
 CONV_RESIDUAL_MINVAL= -10
 CONV_STARTITER= 10

TestCases/disc_adj_fsi/configFlow.cfg

@@ -4,20 +4,14 @@
 % Author: R.Sanchez                                                            %
 % Institution: Imperial College London                                         %
 % Date: 2017.11.29                                                             %
-% File Version 8.1.0 "Harrier"                                                     %
+% File Version 8.1.0 "Harrier"                                                 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 SOLVER= NAVIER_STOKES
 
-MATH_PROBLEM= DISCRETE_ADJOINT
-
-RESTART_SOL= NO
-
-READ_BINARY_RESTART=NO
-
 INNER_ITER= 50
 
-STAT_RELAX_PARAMETER= 1.0
+STAT_RELAX_PARAMETER= 1
 BGS_RELAXATION = FIXED_PARAMETER
 PREDICTOR_ORDER=0
 
@@ -48,7 +42,6 @@ MARKER_MONITORING= ( UpperWall, LowerWall, Wall)
 DEFORM_MESH= YES
 MARKER_DEFORM_MESH= ( UpperWall, LowerWall )
 
-
 DEFORM_STIFFNESS_TYPE = INVERSE_VOLUME
 DEFORM_POISSONS_RATIO = 1e6
 DEFORM_LINEAR_SOLVER = CONJUGATE_GRADIENT

TestCases/disc_adj_rans/naca0012/turb_NACA0012_sa.cfg

@@ -45,6 +45,7 @@ MARKER_MONITORING= ( airfoil )
 % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
 %
 NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+NUM_METHOD_GRAD_RECON= LEAST_SQUARES
 CFL_NUMBER= 10.0
 MAX_DELTA_TIME= 1E10
 CFL_ADAPT= NO
@@ -67,7 +68,6 @@ LINEAR_SOLVER_ITER= 5
 %
 CONV_NUM_METHOD_FLOW= ROE
 MUSCL_FLOW= YES
-SLOPE_LIMITER_FLOW= VENKATAKRISHNAN
 JST_SENSOR_COEFF= ( 0.5, 0.02 )
 TIME_DISCRE_FLOW= EULER_IMPLICIT
 
@@ -88,7 +88,7 @@ CONV_CAUCHY_EPS= 1E-6
 
 % ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
 %
-MESH_FILENAME= n0012_113-33.su2
+MESH_FILENAME= n0012_225-65.su2
 MESH_FORMAT= SU2
 MESH_OUT_FILENAME= mesh_out.su2
 SOLUTION_FILENAME= solution_flow_sa.dat
@@ -103,4 +103,4 @@ GRAD_OBJFUNC_FILENAME= of_grad.dat
 SURFACE_FILENAME= surface_flow
 SURFACE_ADJ_FILENAME= surface_adjoint
 OUTPUT_WRT_FREQ= 1000
-SCREEN_OUTPUT = (INNER_ITER, RMS_ADJ_DENSITY, RMS_ADJ_NU_TILDE, SENS_PRESS, SENS_AOA RMS_DENSITY RMS_NU_TILDE LIFT DRAG LINSOL_ITER LINSOL_RESIDUAL LINSOL_ITER_TURB LINSOL_RESIDUAL_TURB)
+SCREEN_OUTPUT = (INNER_ITER, RMS_ADJ_DENSITY, RMS_ADJ_NU_TILDE, SENS_PRESS, SENS_AOA RMS_DENSITY RMS_NU_TILDE LIFT DRAG LINSOL_ITER LINSOL_RESIDUAL LINSOL_ITER_TURB LINSOL_RESIDUAL_TURB)