fix: Fix five bugs based on reviewer feedback

Fixed three critical division by zero bugs and two moderate bugs:

1. Wall model convergence - zero-check before division
2. Wall model Gamma - protect against very small Gamma (use 1.0E-16 threshold)
3. Radiation emissivity - throw error if outside [0,1] range
4. Restart metadata - correct ITER= offset 9->5
5. Filename buffer - replace strcpy with strncpy

Changes based on reviewer feedback:
- CRadP1Solver.cpp: Error instead of silent clamp for invalid emissivity
- wall_model.cpp: Fixed Gamma logic (always negative, exp_term=0 for small |Gamma|)

Improves solver stability, correctness, and security.

Signed-off-by: shbhmexe <[email protected]>

Author: shbhmexe

Common/src/wall_model.cpp

@@ -394,20 +394,21 @@ void CWallModelLogLaw::WallShearStressAndHeatFlux(const su2double tExchange, con
     const su2double lhs = -((tExchange - TWall) * rho_wall * c_p * u_tau);
     const su2double Gamma = -(0.01 * (Pr_lam * pow(y_plus, 4.0)) / (1.0 + 5.0 * y_plus * pow(Pr_lam, 3.0)));
     const su2double rhs_1 = Pr_lam * y_plus * exp(Gamma);
-    
-    /* Protect against division by zero in exp(1/Gamma) when Gamma is very small */
-    const su2double eps_gamma = 1e-20;
+
+    /* Protect against division by zero in exp(1/Gamma) when Gamma is very small.
+     * Note: Gamma is always negative due to the minus sign in its formula.
+     * When abs(Gamma) is very small, exp(1/Gamma) approaches 0 (since 1/Gamma -> -inf).
+     */
     su2double exp_term;
-    if (abs(Gamma) > eps_gamma) {
+    if (abs(Gamma) > 1.0E-16) {
       exp_term = exp(1.0 / Gamma);
     } else {
-      /* For very small Gamma, cap the exponential to prevent overflow/division by zero */
-      exp_term = (Gamma > 0) ? exp(1.0 / eps_gamma) : exp(-1.0 / eps_gamma);
+      /* For very small |Gamma|, exp(1/Gamma) -> 0 since Gamma is negative */
+      exp_term = 0.0;
     }
-    
+
     const su2double rhs_2 =
-        (2.12 * log(1.0 + y_plus) + pow((3.85 * pow(Pr_lam, (1.0 / 3.0)) - 1.3), 2.0) + 2.12 * log(Pr_lam)) *
-        exp_term;
+        (2.12 * log(1.0 + y_plus) + pow((3.85 * pow(Pr_lam, (1.0 / 3.0)) - 1.3), 2.0) + 2.12 * log(Pr_lam)) * exp_term;
     qWall = lhs / (rhs_1 + rhs_2);
   } else {
     qWall = Wall_HeatFlux;

SU2_CFD/src/solvers/CRadP1Solver.cpp

@@ -282,18 +282,15 @@ void CRadP1Solver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_cont
   /*--- Get the specified wall emissivity from config ---*/
   Wall_Emissivity = config->GetWall_Emissivity(Marker_Tag);
 
-  /*--- Clamp emissivity to valid physical range [0,1] to prevent division by zero ---*/
-  if (Wall_Emissivity < 0.0) {
-    Wall_Emissivity = 0.0;
-  }
-  if (Wall_Emissivity > 1.0) {
-    Wall_Emissivity = 1.0;
+  /*--- Validate emissivity is in physical range [0,1] ---*/
+  if (Wall_Emissivity < 0.0 || Wall_Emissivity > 1.0) {
+    SU2_MPI::Error("Wall emissivity must be in range [0,1].", CURRENT_FUNCTION);
   }
 
   /*--- Compute the constant for the wall theta ---*/
   Theta = Wall_Emissivity / (2.0*(2.0 - Wall_Emissivity));
 
-    /*--- Retrieve the specified wall temperature ---*/
+  /*--- Retrieve the specified wall temperature ---*/
   Twall = config->GetIsothermal_Temperature(Marker_Tag)/config->GetTemperature_Ref();
 
   /*--- Loop over all of the vertices on this boundary marker ---*/
@@ -364,12 +361,9 @@ void CRadP1Solver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container,
   /*--- Get the specified wall emissivity from config ---*/
   Wall_Emissivity = config->GetWall_Emissivity(Marker_Tag);
 
-  /*--- Clamp emissivity to valid physical range [0,1] to prevent division by zero ---*/
-  if (Wall_Emissivity < 0.0) {
-    Wall_Emissivity = 0.0;
-  }
-  if (Wall_Emissivity > 1.0) {
-    Wall_Emissivity = 1.0;
+  /*--- Validate emissivity is in physical range [0,1] ---*/
+  if (Wall_Emissivity < 0.0 || Wall_Emissivity > 1.0) {
+    SU2_MPI::Error("Wall emissivity must be in range [0,1].", CURRENT_FUNCTION);
   }
 
   /*--- Compute the constant for the wall theta ---*/
@@ -446,18 +440,10 @@ void CRadP1Solver::BC_Marshak(CGeometry *geometry, CSolver **solver_container, C
   /*--- Get the specified wall emissivity from config ---*/
   Wall_Emissivity = config->GetWall_Emissivity(Marker_Tag);
 
-  /*--- Clamp emissivity to valid physical range [0,1] to prevent division by zero ---*/
-<<<<<<< HEAD
-  if (Wall_Emissivity < 0.0) {
-    Wall_Emissivity = 0.0;
-  }
-  if (Wall_Emissivity > 1.0) {
-    Wall_Emissivity = 1.0;
+  /*--- Validate emissivity is in physical range [0,1] ---*/
+  if (Wall_Emissivity < 0.0 || Wall_Emissivity > 1.0) {
+    SU2_MPI::Error("Wall emissivity must be in range [0,1].", CURRENT_FUNCTION);
   }
-=======
-  if (Wall_Emissivity < 0.0) Wall_Emissivity = 0.0;
-  if (Wall_Emissivity > 1.0) Wall_Emissivity = 1.0;
->>>>>>> 46baac1286ba393e261fef83bced67d0d4d3f269
 
   /*--- Compute the constant for the wall theta ---*/
   Theta = Wall_Emissivity / (2.0*(2.0 - Wall_Emissivity));