edge-based strategy for reconstructions

pcarruscag · pcarruscag · commit c523946c6120 · 2022-07-02T15:58:45.000-07:00
diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp
@@ -43,6 +43,8 @@ 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/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp
@@ -41,6 +41,8 @@ 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/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp
@@ -148,6 +148,8 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config,
 
   Allocate(*config);
 
+  NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0;
+
   /*--- MPI + OpenMP initialization. ---*/
 
   HybridParallelInitialization(*config, *geometry);
@@ -1829,32 +1831,19 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
       }
       su2double RoeEnthalpy = (R * Primitive_j[prim_idx.Enthalpy()] + Primitive_i[prim_idx.Enthalpy()]) / (R+1);
 
-      bool neg_sound_speed = ((Gamma-1)*(RoeEnthalpy-0.5*sq_vel) < 0.0);
-
-      bool bad_i = neg_sound_speed || neg_pres_or_rho_i;
-      bool bad_j = neg_sound_speed || neg_pres_or_rho_j;
-
-      auto update_nonphysical = [&](){
-        nodes->SetNon_Physical(iPoint, bad_i);
-        nodes->SetNon_Physical(jPoint, bad_j);
-
-        /*--- Get updated state, in case the point recovered after the set. ---*/
-        bad_i = nodes->GetNon_Physical(iPoint);
-        bad_j = nodes->GetNon_Physical(jPoint);
-      };
-
-      if(ReducerStrategy){
-        SU2_OMP_CRITICAL
-        update_nonphysical();
-        END_SU2_OMP_CRITICAL
-      } else {
-        update_nonphysical();
+      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;
       }
+      counter_local += bad_recon;
 
-      counter_local += bad_i+bad_j;
-
-      numerics->SetPrimitive(bad_i? V_i : Primitive_i,  bad_j? V_j : Primitive_j);
-      numerics->SetSecondary(bad_i? S_i : Secondary_i,  bad_j? S_j : Secondary_j);
+      numerics->SetPrimitive(bad_recon? V_i : Primitive_i,  bad_recon? V_j : Primitive_j);
+      numerics->SetSecondary(bad_recon? S_i : Secondary_i,  bad_recon? S_j : Secondary_j);
 
     }
 
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -149,6 +149,8 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned
 
   Allocate(*config);
 
+  NonPhysicalEdgeCounter.resize(geometry->GetnEdge()) = 0;
+
   /*--- MPI + OpenMP initialization. ---*/
 
   HybridParallelInitialization(*config, *geometry);
@@ -471,7 +473,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
 
       case INC_IDEAL_GAS:
         fluidModel = new CIncIdealGas(Specific_Heat_CpND, Gas_ConstantND, Pressure_ThermodynamicND);
-        fluidModel->SetTDState_T(Temperature_FreeStreamND);        
+        fluidModel->SetTDState_T(Temperature_FreeStreamND);
         break;
 
       case FLUID_MIXTURE:
@@ -490,7 +492,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
         }
         fluidModel->SetTDState_T(Temperature_FreeStreamND);
         break;
- 
+
     }
 
     if (viscous) {
@@ -1254,36 +1256,27 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
        checked. Pressure is the dynamic pressure (can be negative). ---*/
 
       if (config->GetEnergy_Equation()) {
-        bool neg_temperature_i = (Primitive_i[prim_idx.Temperature()] < 0.0);
-        bool neg_temperature_j = (Primitive_j[prim_idx.Temperature()] < 0.0);
-
-        bool neg_density_i  = (Primitive_i[prim_idx.Density()] < 0.0);
-        bool neg_density_j  = (Primitive_j[prim_idx.Density()] < 0.0);
-
-        auto update_nonphysical = [&](){
-          nodes->SetNon_Physical(iPoint, neg_density_i || neg_temperature_i);
-          nodes->SetNon_Physical(jPoint, neg_density_j || neg_temperature_j);
-
-          /* Lastly, check for existing first-order points still active from previous iterations. */
+        const bool neg_temperature_i = (Primitive_i[prim_idx.Temperature()] < 0.0);
+        const bool neg_temperature_j = (Primitive_j[prim_idx.Temperature()] < 0.0);
+
+        const bool neg_density_i  = (Primitive_i[prim_idx.Density()] < 0.0);
+        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;
+        }
+        counter_local += bad_recon;
 
-          if (nodes->GetNon_Physical(iPoint)) {
-            counter_local++;
-            for (iVar = 0; iVar < nPrimVar; iVar++)
-              Primitive_i[iVar] = V_i[iVar];
-          }
-          if (nodes->GetNon_Physical(jPoint)) {
-            counter_local++;
-            for (iVar = 0; iVar < nPrimVar; iVar++)
-              Primitive_j[iVar] = V_j[iVar];
+        if (bad_recon) {
+          for (iVar = 0; iVar < nPrimVar; iVar++) {
+            Primitive_i[iVar] = V_i[iVar];
+            Primitive_j[iVar] = V_j[iVar];
           }
-        };
-
-        if(ReducerStrategy){
-          SU2_OMP_CRITICAL
-          update_nonphysical();
-          END_SU2_OMP_CRITICAL
-        } else {
-          update_nonphysical();
         }
       }
 
@@ -2235,7 +2228,7 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
     Flow_Dir = Inlet_FlowDir[val_marker][iVertex];
     Flow_Dir_Mag = GeometryToolbox::Norm(nDim, Flow_Dir);
 
-    /*--- Store the unit flow direction vector. 
+    /*--- Store the unit flow direction vector.
      If requested, use the local boundary normal (negative),
      instead of the prescribed flow direction in the config. ---*/
 
@@ -2352,7 +2345,7 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
         V_inlet[iDim+prim_idx.Velocity()] = nodes->GetVelocity(iPoint,iDim);
       /* pressure obtained from interior */
       V_inlet[prim_idx.Pressure()] = nodes->GetPressure(iPoint);
-    } 
+    }
 
     /*--- Access density at the node. This is either constant by
       construction, or will be set fixed implicitly by the temperature
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
@@ -1055,7 +1055,7 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry,
 #ifdef HAVE_MPI
       /*--- Once we have recv'd a message, get the source rank. ---*/
       int ind;
-      SU2_OMP_SAFE_GLOBAL_ACCESS(SU2_MPI::Waitany(geometry->nPeriodicRecv, eometry->req_PeriodicRecv, &ind, &status);)
+      SU2_OMP_SAFE_GLOBAL_ACCESS(SU2_MPI::Waitany(geometry->nPeriodicRecv, geometry->req_PeriodicRecv, &ind, &status);)
       source = status.MPI_SOURCE;
 #else
       /*--- For serial calculations, we know the rank. ---*/
