From f4467cf14e92d5bc6049e86c409c6eb81adf933c Mon Sep 17 00:00:00 2001 From: Kyle Hanquist Date: Fri, 31 Oct 2025 15:33:53 -0400 Subject: [PATCH 01/11] Trying to fix axisymmetry issue seen in heat flux --- .../include/solvers/CFVMFlowSolverBase.inl | 8 ++++++ SU2_CFD/src/numerics/flow/flow_sources.cpp | 27 +++++++++++++++---- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index aea0dc1d6a26..1a132f5008e1 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2939,12 +2939,20 @@ void CFVMFlowSolverBase::ComputeAxisymmetricAuxGradients(CGeometr 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){ 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 { + /*--- 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 diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 8cf9642b4251..db5028873e8c 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -223,17 +223,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). */ + 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); From 1b53ab2f06445c2472580d0ebc777f17c5dcc9b6 Mon Sep 17 00:00:00 2001 From: Kyle Hanquist Date: Fri, 31 Oct 2025 15:45:10 -0400 Subject: [PATCH 02/11] Clarifying the change --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index db5028873e8c..af61d6eb84a7 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -241,7 +241,7 @@ CNumerics::ResidualType<> CSourceGeneralAxisymmetric_Flow::ComputeResidual(const residual[3] = Volume * rho * Enthalpy_i * dv_dr; // ρH(∂v/∂r) if (implicit) { - /* For now, set Jacobian to zero at axis (can be improved later). */ + /* 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; From a26ee2a574dd889319a29d2b53a36a7f95401083 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Wed, 21 Jan 2026 19:15:49 -0500 Subject: [PATCH 03/11] Axisymmetric source term: blend dv/dr with v/r near axis --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 93 ++++++++++++++-------- 1 file changed, 60 insertions(+), 33 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index af61d6eb84a7..c4d178929390 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -63,27 +63,65 @@ 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) { - - 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; - } + /*--- Common calculations for both branches ---*/ + su2double rho = U_i[0]; // density + su2double u = U_i[1]/U_i[0]; // u-velocity + su2double v = U_i[2]/U_i[0]; // v-velocity + su2double r = Coord_i[1]; // radial coordinate + su2double dv_dr = PrimVar_Grad_i[2][1]; // ∂v/∂r (radial velocity gradient) + + 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 + } - 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]; + /*--- 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; + } - 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]; + /*--- 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) - /*--- Inviscid component of the source term. ---*/ + /*--- 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]; - 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 +145,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); } From 80027c21991a2904bf6c65fc7b5ebb59f2bf53c9 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Wed, 8 Apr 2026 20:22:34 -0400 Subject: [PATCH 04/11] Remove redundant velocity definition --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index c4d178929390..f07475628b4e 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -66,7 +66,6 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi /*--- Common calculations for both branches ---*/ su2double rho = U_i[0]; // density su2double u = U_i[1]/U_i[0]; // u-velocity - su2double v = U_i[2]/U_i[0]; // v-velocity su2double r = Coord_i[1]; // radial coordinate su2double dv_dr = PrimVar_Grad_i[2][1]; // ∂v/∂r (radial velocity gradient) From 060b82321fff15463ac29c5c3d5ae5542df6553a Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Tue, 22 Sep 2026 14:57:31 -0400 Subject: [PATCH 05/11] Update regression values for axisymmetric source term --- TestCases/hybrid_regression.py | 2 +- TestCases/parallel_regression.py | 4 ++-- TestCases/parallel_regression_AD.py | 2 +- TestCases/serial_regression.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 149f4a921607..f63748ba9283 100644 --- a/TestCases/hybrid_regression.py +++ b/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) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 36766da22528..85e6d3f5d64d 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -225,7 +225,7 @@ def main(): ion_gy.cfg_dir = "nonequilibrium/visc_cylinder" ion_gy.cfg_file = "cyl_ion_gy.cfg" ion_gy.test_iter = 99 - ion_gy.test_vals = [-12.191649, -4.245292, -4.904190, -5.585803, -5.472455, -5.057619, -7.442352, 3.429183, 3.433581, -0.014861, 0.000001, 90357.000000] + ion_gy.test_vals = [-11.673961, -4.180906, -4.838570, -5.457224, -5.150665, -4.930094, -6.962117, 4.562850, 4.581001, -0.014893, -0.000032, 90471.000000] ion_gy.tol = 0.01 test_list.append(ion_gy) @@ -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) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index 0b5d8724f62b..26a32c3e2762 100644 --- a/TestCases/parallel_regression_AD.py +++ b/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) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index ab1de7fe1431..435e32bf6744 100755 --- a/TestCases/serial_regression.py +++ b/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) From cd725f1834bf1a8f405844f4f52bc13eec1ebfb3 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Wed, 23 Sep 2026 12:59:27 -0400 Subject: [PATCH 06/11] Restore ion_gy regression values --- TestCases/parallel_regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 85e6d3f5d64d..e85b5b7be660 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -225,7 +225,7 @@ def main(): ion_gy.cfg_dir = "nonequilibrium/visc_cylinder" ion_gy.cfg_file = "cyl_ion_gy.cfg" ion_gy.test_iter = 99 - ion_gy.test_vals = [-11.673961, -4.180906, -4.838570, -5.457224, -5.150665, -4.930094, -6.962117, 4.562850, 4.581001, -0.014893, -0.000032, 90471.000000] + ion_gy.test_vals = [-12.191649, -4.245292, -4.904190, -5.585803, -5.472455, -5.057619, -7.442352, 3.429183, 3.433581, -0.014861, 0.000001, 90357.000000] ion_gy.tol = 0.01 test_list.append(ion_gy) From b765070ffcca4261aecceda124323f7eeceb363e Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Thu, 24 Sep 2026 19:03:13 -0400 Subject: [PATCH 07/11] Exclude NEMO from axisymmetric auxiliary correction --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 1a132f5008e1..cb44757c6280 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -2932,6 +2932,7 @@ su2double CFVMFlowSolverBase::EvaluateCommonObjFunc(const CConfig& config) template void CFVMFlowSolverBase::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++) { @@ -2941,18 +2942,20 @@ void CFVMFlowSolverBase::ComputeAxisymmetricAuxGradients(CGeometr 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 { - /*--- 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) + 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 From 503e6fdfc635872389aa735114b418a97efce3fd Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Fri, 25 Sep 2026 15:08:48 -0400 Subject: [PATCH 08/11] Remove trailing whitespace --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index f07475628b4e..c9921be6fe14 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -68,7 +68,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi 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) - + sq_vel = 0.0; for (iDim = 0; iDim < nDim; iDim++) { Velocity_i = U_i[iDim+1] / U_i[0]; @@ -80,7 +80,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi /*--- 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) { @@ -257,7 +257,7 @@ CNumerics::ResidualType<> CSourceGeneralAxisymmetric_Flow::ComputeResidual(const /* 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 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 ---*/ From acbc66af38f33055f4bfc4ab1ae0406b2c8c56cb Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Mon, 28 Sep 2026 15:37:21 -0400 Subject: [PATCH 09/11] Simplify inverse radius handling --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index c9921be6fe14..c1158cc91f1f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -92,8 +92,8 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi /*--- Standard formulation (v/r) ---*/ su2double std_res[4]; + yinv = (r > EPS) ? 1.0/r : 0.0; 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 @@ -120,7 +120,6 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi 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; From d3260087b023893caf0672bf78debe328f9e6a31 Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Tue, 29 Sep 2026 14:24:25 -0400 Subject: [PATCH 10/11] Update axisymmetric species regression values --- TestCases/serial_regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 435e32bf6744..c34fe6f47225 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -377,7 +377,7 @@ def main(): axi_rans_air_nozzle_species.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_species.cfg_file = "air_nozzle_species.cfg" axi_rans_air_nozzle_species.test_iter = 10 - axi_rans_air_nozzle_species.test_vals = [-1.690665, 3.882506, -2.928702, 5.760933, -3.144560, 0.000000] + axi_rans_air_nozzle_species.test_vals = [-1.676190, 3.896581, -2.912396, 5.762790, -3.565909, 0.000000] axi_rans_air_nozzle_species.tol = 0.0001 test_list.append(axi_rans_air_nozzle_species) From 4c1f9b2b2531baa5beef4372242b8edddac9c38a Mon Sep 17 00:00:00 2001 From: Raghava Davuluri Date: Tue, 29 Sep 2026 15:21:15 -0400 Subject: [PATCH 11/11] Fix inverse radius initialization for AD builds --- SU2_CFD/src/numerics/flow/flow_sources.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index c1158cc91f1f..e3797c7c5a5f 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -92,8 +92,9 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi /*--- Standard formulation (v/r) ---*/ su2double std_res[4]; - yinv = (r > EPS) ? 1.0/r : 0.0; + yinv = 0.0; 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