use the non physical counter in SIMD numerics

RichRoos · Sep 10, 2022 · 6c53381 · 6c53381
1 parent 7e35912
commit 6c53381
Show file tree

Hide file tree

Showing 10 changed files with 107 additions and 91 deletions.
diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp
@@ -93,15 +93,15 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
 
 /*!
  * \brief Retrieve primitive variables for points i/j, reconstructing them if needed.
- * \param[in] iPoint, jPoint - Nodes of the edge.
+ * \param[in] iEdge, iPoint, jPoint - Edge and its nodes.
  * \param[in] muscl - If true, reconstruct, else simply fetch.
  * \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 iPoint, Int jPoint, bool muscl,
-                                                      LIMITER limiterType,
+FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int jPoint,
+                                                      bool muscl, LIMITER limiterType,
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& vector_ij,
                                                       const VariableType& solution) {
@@ -131,7 +131,31 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iPoint, Int jPoint, bo
       musclPointLimited(jPoint, vector_ij,-0.5, limiters, gradients, V.j.all);
       break;
     }
-    /// TODO: Extra reconstruction checks needed.
+    /*--- Detect a non-physical reconstruction based on negative pressure or density. ---*/
+    const Double neg_p_or_rho = max(min(V.i.pressure(), V.j.pressure()) < 0.0,
+                                    min(V.i.density(), V.j.density()) < 0.0);
+    /*--- Test the sign of the Roe-averaged speed of sound. ---*/
+    const Double R = sqrt(V.j.density() / V.i.density());
+    /*--- Delay dividing by R+1 until comparing enthalpy and velocity magnitude. ---*/
+    const Double enthalpy = R*V.j.enthalpy() + V.i.enthalpy();
+    Double v_squared = 0.0;
+    for (size_t iDim = 0; iDim < nDim; ++iDim) {
+      v_squared += pow(R*V.j.velocity(iDim) + V.i.velocity(iDim), 2);
+    }
+    /*--- Multiply enthalpy by R+1 since v^2 was not divided by (R+1)^2.
+     * Note: a = sqrt((gamma-1) * (H - 0.5 * v^2)) ---*/
+    const Double neg_sound_speed = enthalpy * (R+1) < 0.5 * v_squared;
+
+    /*--- Revert to first order if the state is non-physical. ---*/
+    Double bad_recon = max(neg_p_or_rho, neg_sound_speed);
+    /*--- Handle SIMD dimensions 1 by 1. ---*/
+    for (size_t k = 0; k < Double::Size; ++k) {
+      bad_recon[k] = solution.UpdateNonPhysicalEdgeCounter(iEdge[k], bad_recon[k]);
+    }
+    for (size_t iVar = 0; iVar < ReconVarType::nVar; ++iVar) {
+      V.i.all(iVar) = bad_recon * V1st.i.all(iVar) + (1-bad_recon) * V.i.all(iVar);
+      V.j.all(iVar) = bad_recon * V1st.j.all(iVar) + (1-bad_recon) * V.j.all(iVar);
+    }
   }
   return V;
 }

diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp
@@ -116,7 +116,7 @@ class CRoeBase : public Base {
     V1st.j.all = gatherVariables<nPrimVar>(jPoint, solution.GetPrimitive());
 
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-                  iPoint, jPoint, muscl, typeLimiter, V1st, vector_ij, solution);
+                 iEdge, iPoint, jPoint, muscl, typeLimiter, V1st, vector_ij, solution);
 
     /*--- Compute conservative variables. ---*/
 

diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp
@@ -43,8 +43,6 @@ class CEulerSolver : public CFVMFlowSolverBase<CEulerVariable, ENUM_REGIME::COMP
   Prandtl_Lam = 0.0,    /*!< \brief Laminar Prandtl number. */
   Prandtl_Turb = 0.0;   /*!< \brief Turbulent Prandtl number. */
 
-  su2vector<int8_t> NonPhysicalEdgeCounter;  /*!< \brief Non-physical reconstruction counter for each edge. */
-
   su2double AllBound_CEquivArea_Inv=0.0; /*!< \brief equivalent area coefficient (inviscid contribution) for all the boundaries. */
   vector<su2double> CEquivArea_Mnt;      /*!< \brief Equivalent area (inviscid contribution) for each boundary. */
   vector<su2double> CEquivArea_Inv;      /*!< \brief Equivalent area (inviscid contribution) for each boundary. */

diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
@@ -269,13 +269,50 @@ class CFVMFlowSolverBase : public CSolver {
   /*!
    * \brief Method to compute convective and viscous residual contribution using vectorized numerics.
    */
-  void EdgeFluxResidual(const CGeometry *geometry, const CSolver* const* solvers, const CConfig *config);
+  void EdgeFluxResidual(const CGeometry *geometry, const CSolver* const* solvers, 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 Sums edge fluxes (if required) and computes the global error counter.
+   * \param[in] pausePreacc - Whether preaccumulation was paused durin.
+   * \param[in] localCounter - Thread-local error counter.
+   * \param[in,out] config - Used to set the global error counter.
+   */
+  inline void FinalizeResidualComputation(const CGeometry *geometry, bool pausePreacc,
+                                          unsigned long localCounter, CConfig* config) {
+
+    /*--- Restore preaccumulation and adjoint evaluation state. ---*/
+    AD::ResumePreaccumulation(pausePreacc);
+    if (!ReducerStrategy) AD::EndNoSharedReading();
+
+    if (ReducerStrategy) {
+      SumEdgeFluxes(geometry);
+      if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) {
+        Jacobian.SetDiagonalAsColumnSum();
+      }
+    }
+
+    /*--- Warning message about non-physical reconstructions. ---*/
+    if ((MGLevel == MESH_0) && (config->GetComm_Level() == COMM_FULL)) {
+      /*--- Add counter results for all threads. ---*/
+      SU2_OMP_ATOMIC
+      ErrorCounter += localCounter;
+
+      /*--- Add counter results for all ranks. ---*/
+      BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
+      {
+        localCounter = ErrorCounter;
+        SU2_MPI::Reduce(&localCounter, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm());
+        config->SetNonphysical_Reconstr(ErrorCounter);
+      }
+      END_SU2_OMP_SAFE_GLOBAL_ACCESS
+    }
+  }
+
   /*!
    * \brief Computes and sets the required auxilliary vars (and gradients) for axisymmetric flow.
    */

diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -1575,11 +1575,17 @@ 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 CSolver* const* solvers,
-                                                const CConfig *config) {
+                                                CConfig *config) {
   if (!edgeNumerics) {
     InstantiateEdgeNumerics(solvers, config);
   }
 
+  /*--- Non-physical counter. ---*/
+  unsigned long counterLocal = 0;
+  SU2_OMP_MASTER
+  ErrorCounter = 0;
+  END_SU2_OMP_MASTER
+
   /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of
   * variables, otherwise switch to the faster adjoint evaluation mode. ---*/
   bool pausePreacc = false;
@@ -1604,20 +1610,14 @@ void CFVMFlowSolverBase<V, R>::EdgeFluxResidual(const CGeometry *geometry,
       } else {
         edgeNumerics->ComputeFlux(iEdge, *config, *geometry, *nodes, UpdateType::COLORING, mask, LinSysRes, Jacobian);
       }
+      for (auto j = 0ul; j < Double::Size; ++j) {
+        counterLocal += (nodes->NonPhysicalEdgeCounter[iEdge[j]] > 0);
+      }
     }
     END_SU2_OMP_FOR
   }
 
-  /*--- Restore preaccumulation and adjoint evaluation state. ---*/
-  AD::ResumePreaccumulation(pausePreacc);
-  if (!ReducerStrategy) AD::EndNoSharedReading();
-
-  if (ReducerStrategy) {
-    SumEdgeFluxes(geometry);
-    if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) {
-      Jacobian.SetDiagonalAsColumnSum();
-    }
-  }
+  FinalizeResidualComputation(geometry, pausePreacc, counterLocal, config);
 }
 
 template <class V, ENUM_REGIME R>

diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp
@@ -41,8 +41,6 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
   vector<CFluidModel*> FluidModel;   /*!< \brief fluid model used in the solver. */
   StreamwisePeriodicValues SPvals, SPvalsUpdated;
 
-  su2vector<int8_t> NonPhysicalEdgeCounter;  /*!< \brief Non-physical reconstruction counter for each edge. */
-
   /*!
    * \brief Preprocessing actions common to the Euler and NS solvers.
    * \param[in] geometry - Geometrical definition of the problem.

diff --git a/SU2_CFD/include/variables/CFlowVariable.hpp b/SU2_CFD/include/variables/CFlowVariable.hpp
@@ -66,6 +66,23 @@ class CFlowVariable : public CVariable {
                 unsigned long nprimvargrad, const CConfig* config);
 
  public:
+  mutable su2vector<int8_t> NonPhysicalEdgeCounter;  /*!< \brief Non-physical reconstruction counter for each edge. */
+  /*!
+   * \brief Updates the non-physical counter of an edge.
+   * \param[in] iEdge - Edge index.
+   * \param[in] isNonPhys - Should be 1 (true) if a non-physical reconstruction was occurred.
+   * \return Whether the reconstruction should be limited to first order, based on the counter.
+   */
+  template <class T>
+  inline T UpdateNonPhysicalEdgeCounter(unsigned long iEdge, const T& isNonPhys) const {
+    if (isNonPhys) {
+      /*--- Force 1st order for this edge for at least 20 iterations. ---*/
+      NonPhysicalEdgeCounter[iEdge] = 21;
+    }
+    NonPhysicalEdgeCounter[iEdge] = std::max<int8_t>(0, NonPhysicalEdgeCounter[iEdge] - 1);
+    return static_cast<T>(NonPhysicalEdgeCounter[iEdge] > 0);
+  }
+
   /*!
    * \brief Get a primitive variable.
    * \param[in] iPoint - Point index.

diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -148,7 +148,9 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
 
   Allocate(*config);
 
-  NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0;
+  if (iMesh == MESH_0) {
+    nodes->NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0;
+  }
 
   /*--- MPI + OpenMP initialization. ---*/
 
@@ -1816,13 +1818,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
 
       const bool neg_sound_speed = ((Gamma-1)*(RoeEnthalpy-0.5*sq_vel) < 0.0);
       bool bad_recon = neg_sound_speed || neg_pres_or_rho_i || neg_pres_or_rho_j;
-      if (bad_recon) {
-        /*--- Force 1st order for this edge for at least 20 iterations. ---*/
-        NonPhysicalEdgeCounter[iEdge] = 20;
-      } else if (NonPhysicalEdgeCounter[iEdge] > 0) {
-        --NonPhysicalEdgeCounter[iEdge];
-        bad_recon = true;
-      }
+      bad_recon = nodes->UpdateNonPhysicalEdgeCounter(iEdge, bad_recon);
       counter_local += bad_recon;
 
       numerics->SetPrimitive(bad_recon? V_i : Primitive_i,  bad_recon? V_j : Primitive_j);
@@ -1881,33 +1877,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
   END_SU2_OMP_FOR
   } // end color loop
 
-  /*--- Restore preaccumulation and adjoint evaluation state. ---*/
-  AD::ResumePreaccumulation(pausePreacc);
-  if (!ReducerStrategy) AD::EndNoSharedReading();
-
-  if (ReducerStrategy) {
-    SumEdgeFluxes(geometry);
-    if (implicit)
-      Jacobian.SetDiagonalAsColumnSum();
-  }
-
-  /*--- Warning message about non-physical reconstructions. ---*/
-
-  if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) {
-    /*--- Add counter results for all threads. ---*/
-    SU2_OMP_ATOMIC
-    ErrorCounter += counter_local;
-
-    /*--- Add counter results for all ranks. ---*/
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
-    {
-      counter_local = ErrorCounter;
-      SU2_MPI::Reduce(&counter_local, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm());
-      config->SetNonphysical_Reconstr(ErrorCounter);
-    }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
-  }
-
+  FinalizeResidualComputation(geometry, pausePreacc, counter_local, config);
 }
 
 void CEulerSolver::ComputeConsistentExtrapolation(CFluidModel *fluidModel, unsigned short nDim,

diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -149,7 +149,9 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned
 
   Allocate(*config);
 
-  NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0;
+  if (iMesh == MESH_0) {
+    nodes->NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0;
+  }
 
   /*--- MPI + OpenMP initialization. ---*/
 
@@ -1247,13 +1249,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
         const bool neg_density_j  = (Primitive_j[prim_idx.Density()] < 0.0);
 
         bool bad_recon = neg_temperature_i || neg_temperature_j || neg_density_i || neg_density_j;
-        if (bad_recon) {
-          /*--- Force 1st order for this edge for at least 20 iterations. ---*/
-          NonPhysicalEdgeCounter[iEdge] = 20;
-        } else if (NonPhysicalEdgeCounter[iEdge] > 0) {
-          --NonPhysicalEdgeCounter[iEdge];
-          bad_recon = true;
-        }
+        bad_recon = nodes->UpdateNonPhysicalEdgeCounter(iEdge, bad_recon);
         counter_local += bad_recon;
 
         if (bad_recon) {
@@ -1302,33 +1298,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
   END_SU2_OMP_FOR
   } // end color loop
 
-  /*--- Restore preaccumulation and adjoint evaluation state. ---*/
-  AD::ResumePreaccumulation(pausePreacc);
-  if (!ReducerStrategy) AD::EndNoSharedReading();
-
-  if (ReducerStrategy) {
-    SumEdgeFluxes(geometry);
-    if (implicit)
-      Jacobian.SetDiagonalAsColumnSum();
-  }
-
-  /*--- Warning message about non-physical reconstructions. ---*/
-
-  if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) {
-    /*--- Add counter results for all threads. ---*/
-    SU2_OMP_ATOMIC
-    ErrorCounter += counter_local;
-
-    /*--- Add counter results for all ranks. ---*/
-    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS
-    {
-      counter_local = ErrorCounter;
-      SU2_MPI::Reduce(&counter_local, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, MASTER_NODE, SU2_MPI::GetComm());
-      config->SetNonphysical_Reconstr(ErrorCounter);
-    }
-    END_SU2_OMP_SAFE_GLOBAL_ACCESS
-  }
-
+  FinalizeResidualComputation(geometry, pausePreacc, counter_local, config);
 }
 
 void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_container,

diff --git a/SU2_GEO/src/SU2_GEO.cpp b/SU2_GEO/src/SU2_GEO.cpp
@@ -677,6 +677,8 @@ int main(int argc, char *argv[]) {
 
         for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) {
 
+          Local_MoveSurface = false;
+
           switch ( config_container[ZONE_0]->GetDesign_Variable(iDV) ) {
             case FFD_CONTROL_POINT_2D : Local_MoveSurface = surface_movement->SetFFDCPChange_2D(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], FFDBox, iDV, true); break;
             case FFD_CAMBER_2D :        Local_MoveSurface = surface_movement->SetFFDCamber_2D(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], FFDBox, iDV, true); break;