Axisymmetric near-axis fixes using L’Hôpital-consistent formulations

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2932,19 +2932,30 @@ su2double CFVMFlowSolverBase<V,R>::EvaluateCommonObjFunc(const CConfig& config)
 template <class V, ENUM_REGIME FlowRegime>
 void CFVMFlowSolverBase<V, FlowRegime>::ComputeAxisymmetricAuxGradients(CGeometry *geometry, const CConfig* config) {
 
+  const bool nemo = config->GetNEMOProblem();
   /*--- Loop through all points to set the auxvargrad --*/
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
     su2double yCoord          = geometry->nodes->GetCoord(iPoint, 1);
     su2double yVelocity       = nodes->GetVelocity(iPoint,1);
     su2double xVelocity       = nodes->GetVelocity(iPoint,0);
     su2double Total_Viscosity = nodes->GetLaminarViscosity(iPoint) + nodes->GetEddyViscosity(iPoint);
+    const auto VelocityGradient = nodes->GetVelocityGradient(iPoint);
 
-    if (yCoord > EPS){
+    if (yCoord > EPS) {
       su2double nu_v_on_y = Total_Viscosity*yVelocity/yCoord;
       nodes->SetAuxVar(iPoint, 0, nu_v_on_y);
       nodes->SetAuxVar(iPoint, 1, nu_v_on_y*yVelocity);
       nodes->SetAuxVar(iPoint, 2, nu_v_on_y*xVelocity);
+    } else {
+      if (!nemo) {
+        /*--- At the axis of symmetry, use L'Hôpital's rule instead of setting each to zero: lim(v/r) = dv/dr ---*/
+        su2double dv_dr = VelocityGradient(1,1);  // ∂v/∂r
+        su2double nu_dv_dr = Total_Viscosity * dv_dr;
+        nodes->SetAuxVar(iPoint, 0, nu_dv_dr);           // μ(∂v/∂r)
+        nodes->SetAuxVar(iPoint, 1, 0.0);                // μv(∂v/∂r) = 0 since v=0 at axis
+        nodes->SetAuxVar(iPoint, 2, nu_dv_dr*xVelocity); // μu(∂v/∂r)
+      }
     }
   }
   END_SU2_OMP_FOR

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -63,27 +63,64 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi
   su2double Pressure_i, Enthalpy_i, Velocity_i, sq_vel;
   unsigned short iDim, iVar, jVar;
 
-  if (Coord_i[1] > EPS) {
+  /*--- Common calculations for both branches ---*/
+  su2double rho = U_i[0];                     // density
+  su2double u = U_i[1]/U_i[0];                // u-velocity
+  su2double r = Coord_i[1];                   // radial coordinate
+  su2double dv_dr = PrimVar_Grad_i[2][1];     // ∂v/∂r (radial velocity gradient)
 
-    yinv = 1.0/Coord_i[1];
+  sq_vel = 0.0;
+  for (iDim = 0; iDim < nDim; iDim++) {
+    Velocity_i = U_i[iDim+1] / U_i[0];
+    sq_vel += Velocity_i * Velocity_i;
+  }
+  Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel);
+  Enthalpy_i = (U_i[nDim+1] + Pressure_i) / U_i[0];
+
+  /*--- Smooth blending between gradient formulation and standard formulation ---*/
+  su2double transition_width = 50.0 * EPS;  // Smooth transition over 50×EPS (much wider)
+  su2double alpha = 0.0;  // Blending factor: 0=gradient_form, 1=standard_form
+
+  if (r > transition_width) {
+    alpha = 1.0;  // Far from axis: use standard v/r formulation
+  } else if (r > EPS) {
+    // Smooth transition zone: blend between formulations (gentler slope)
+    alpha = 0.5 * (1.0 + tanh(2.0 * (r - 0.5*(EPS + transition_width))/(transition_width - EPS)));
+  } else {
+    alpha = 0.0;  // Near axis: use gradient formulation
+  }
 
-    sq_vel = 0.0;
-    for (iDim = 0; iDim < nDim; iDim++) {
-      Velocity_i = U_i[iDim+1] / U_i[0];
-      sq_vel += Velocity_i *Velocity_i;
-    }
+  /*--- Standard formulation (v/r) ---*/
+  su2double std_res[4];
+  if (r > EPS) {
+    yinv = 1.0/r;
+    std_res[0] = yinv*Volume*U_i[2];                    // ρv/r
+    std_res[1] = yinv*Volume*U_i[1]*U_i[2]/U_i[0];     // ρuv/r
+    std_res[2] = yinv*Volume*(U_i[2]*U_i[2]/U_i[0]);   // ρv²/r
+    std_res[3] = yinv*Volume*Enthalpy_i*U_i[2];        // ρHv/r
+  } else {
+    // Avoid division by zero, set to zero (will be blended out anyway)
+    std_res[0] = std_res[1] = std_res[2] = std_res[3] = 0.0;
+  }
 
-    Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel);
-    Enthalpy_i = (U_i[nDim+1] + Pressure_i) / U_i[0];
+  /*--- Gradient formulation (∂v/∂r) ---*/
+  su2double grad_res[4];
+  grad_res[0] = Volume * rho * dv_dr;              // ρ(∂v/∂r)
+  grad_res[1] = Volume * rho * u * dv_dr;          // ρu(∂v/∂r)
+  grad_res[2] = 0.0;                               // ρv(∂v/∂r) = 0 since v→0 as r→0
+  grad_res[3] = Volume * rho * Enthalpy_i * dv_dr; // ρH(∂v/∂r)
 
-    residual[0] = yinv*Volume*U_i[2];
-    residual[1] = yinv*Volume*U_i[1]*U_i[2]/U_i[0];
-    residual[2] = yinv*Volume*(U_i[2]*U_i[2]/U_i[0]);
-    residual[3] = yinv*Volume*Enthalpy_i*U_i[2];
+  /*--- Blend the two formulations ---*/
+  residual[0] = (1.0 - alpha) * grad_res[0] + alpha * std_res[0];
+  residual[1] = (1.0 - alpha) * grad_res[1] + alpha * std_res[1];
+  residual[2] = (1.0 - alpha) * grad_res[2] + alpha * std_res[2];
+  residual[3] = (1.0 - alpha) * grad_res[3] + alpha * std_res[3];
 
-    /*--- Inviscid component of the source term. ---*/
-
-    if (implicit) {
+  /*--- Jacobian calculation ---*/
+  if (implicit) {
+    if (alpha > 0.5) {
+      // Use standard Jacobian when mostly in standard formulation
+      yinv = 1.0/r;
       jacobian[0][0] = 0.0;
       jacobian[0][1] = 0.0;
       jacobian[0][2] = 1.0;
@@ -107,29 +144,18 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi
       for (iVar=0; iVar < nVar; iVar++)
         for (jVar=0; jVar < nVar; jVar++)
           jacobian[iVar][jVar] *= yinv*Volume;
-
-    }
-
-    /*--- Add the viscous terms if necessary. ---*/
-
-    if (viscous) ResidualDiffusion();
-
-  }
-
-  else {
-
-    for (iVar=0; iVar < nVar; iVar++)
-      residual[iVar] = 0.0;
-
-    if (implicit) {
+    } else {
+      // Near axis: set Jacobian to zero (gradient formulation is more complex)
       for (iVar=0; iVar < nVar; iVar++) {
         for (jVar=0; jVar < nVar; jVar++)
           jacobian[iVar][jVar] = 0.0;
       }
     }
-
   }
 
+  /*--- Add the viscous terms if necessary. ---*/
+  if (viscous) ResidualDiffusion();
+
   return ResidualType<>(residual, jacobian, nullptr);
 }
 
@@ -223,17 +249,34 @@ CNumerics::ResidualType<> CSourceGeneralAxisymmetric_Flow::ComputeResidual(const
   }
 
   else {
-
-    for (iVar=0; iVar < nVar; iVar++)
-      residual[iVar] = 0.0;
+    /*--- At the axis of symmetry, use L'Hôpital's rule: lim(v/r) = dv/dr ---*/
+    const su2double dv_dr = PrimVar_Grad_i[2][1];  // ∂v/∂r (radial velocity gradient)
+    const su2double u = U_i[1]/U_i[0];            // axial velocity u
+    const su2double rho = U_i[0];                 // density
+
+    /* Compute pressure and enthalpy consistently with the general-gas formulation. */
+    const su2double Density_i = rho;
+    const su2double Energy_i  = U_i[3]/U_i[0];
+    const su2double Pressure_i = V_j[3];
+    const su2double Enthalpy_i = Energy_i + Pressure_i/Density_i;
+
+    /*--- Apply L'Hôpital's rule to axisymmetric source terms ---*/
+    residual[0] = Volume * rho * dv_dr;               // ρ(∂v/∂r)
+    residual[1] = Volume * rho * u * dv_dr;           // ρu(∂v/∂r)
+    residual[2] = 0.0;                                // ρv(∂v/∂r) = 0 since v=0 at axis
+    residual[3] = Volume * rho * Enthalpy_i * dv_dr;  // ρH(∂v/∂r)
 
     if (implicit) {
-      for (iVar=0; iVar < nVar; iVar++) {
-        for (jVar=0; jVar < nVar; jVar++)
+      /* For now, set Jacobian to zero at axis (can be improved later to help with convergence). */
+      for (iVar = 0; iVar < nVar; iVar++) {
+        for (jVar = 0; jVar < nVar; jVar++)
           jacobian[iVar][jVar] = 0.0;
       }
     }
 
+    /*--- Add the viscous terms if necessary. ---*/
+    if (viscous) ResidualDiffusion();
+
   }
 
   return ResidualType<>(residual, jacobian, nullptr);

TestCases/hybrid_regression.py

@@ -265,7 +265,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-11.083068, -5.374686, -8.880093, -4.073548, 0.000000]
+    axi_rans_air_nozzle_restart.test_vals = [-2.663059, 2.911787, -2.522191, 1.990732, 0.000000]
     axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000]
     test_list.append(axi_rans_air_nozzle_restart)
 

TestCases/parallel_regression.py

@@ -619,7 +619,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-11.056236, -5.334116, -8.842310, -4.067917, 0.000000]
+    axi_rans_air_nozzle_restart.test_vals = [-2.662880, 2.912015, -2.712454, 1.890642, 0.000000]
     axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143310, -9.163287, -10.858232, -5.787715, 0.000000]
     axi_rans_air_nozzle_restart.tol       = 0.0001
     test_list.append(axi_rans_air_nozzle_restart)

TestCases/parallel_regression_AD.py

@@ -150,7 +150,7 @@ def main():
     discadj_axisymmetric_rans_nozzle.cfg_dir    = "axisymmetric_rans/air_nozzle"
     discadj_axisymmetric_rans_nozzle.cfg_file   = "air_nozzle_restart.cfg"
     discadj_axisymmetric_rans_nozzle.test_iter  = 10
-    discadj_axisymmetric_rans_nozzle.test_vals  = [9.909657, 5.078045, 7.129068, 2.490955]
+    discadj_axisymmetric_rans_nozzle.test_vals  = [9.911310, 5.072379, 7.118130, 2.484452]
     discadj_axisymmetric_rans_nozzle.no_restart = True
     test_list.append(discadj_axisymmetric_rans_nozzle)
 

TestCases/serial_regression.py

@@ -367,7 +367,7 @@ def main():
     axi_rans_air_nozzle_restart.cfg_dir   = "axisymmetric_rans/air_nozzle"
     axi_rans_air_nozzle_restart.cfg_file  = "air_nozzle_restart.cfg"
     axi_rans_air_nozzle_restart.test_iter = 10
-    axi_rans_air_nozzle_restart.test_vals = [-11.054281, -5.328905, -8.835591, -4.056830, 0.000000]
+    axi_rans_air_nozzle_restart.test_vals = [-2.662936, 2.911899, -2.338516, 2.161131, 0.000000]
     axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143715, -9.170705, -10.848554, -5.776746, 0.000000]
     axi_rans_air_nozzle_restart.tol       = 0.0001
     test_list.append(axi_rans_air_nozzle_restart)