From d2586c0472a54c458b17a9d3f3b71b7012a8cbb4 Mon Sep 17 00:00:00 2001 From: Cristopher-Morales Date: Tue, 16 Dec 2025 15:49:08 +0100 Subject: [PATCH 1/4] initial commit --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 44 +++++++++---------------- 1 file changed, 16 insertions(+), 28 deletions(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 246da78c22b..5b5591177c2 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2799,7 +2799,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver su2double U_time_nM1[MAXNVAR], U_time_n[MAXNVAR], U_time_nP1[MAXNVAR]; su2double Volume_nM1, Volume_nP1, TimeStep; const su2double *Normal = nullptr, *GridVel_i = nullptr, *GridVel_j = nullptr; - su2double Density, Cp; + su2double Density; const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST); @@ -2807,10 +2807,9 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver const bool energy = config->GetEnergy_Equation(); const int ndim = nDim; - auto V2U = [ndim](su2double Density, su2double Cp, const su2double* V, su2double* U) { + auto V2U = [ndim](su2double Density, const su2double* V, su2double* U) { U[0] = Density; - for (int iDim = 0; iDim < ndim; iDim++) U[iDim+1] = Density*V[iDim+1]; - U[ndim+1] = Density*Cp*V[ndim+1]; + for (int iDim = 0; iDim <= ndim; iDim++) U[iDim + 1] = Density * V[iDim + 1]; }; /*--- Store the physical time step ---*/ @@ -2837,16 +2836,15 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = nodes->GetSolution_time_n(iPoint); V_time_nP1 = nodes->GetSolution(iPoint); - /*--- Access the density and Cp at this node (constant for now). ---*/ + /*--- Access the density at this node (constant for now). ---*/ Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); /*--- Compute the conservative variable vector for all time levels. ---*/ - V2U(Density, Cp, V_time_nM1, U_time_nM1); - V2U(Density, Cp, V_time_n, U_time_n); - V2U(Density, Cp, V_time_nP1, U_time_nP1); + V2U(Density, V_time_nM1, U_time_nM1); + V2U(Density, V_time_n, U_time_n); + V2U(Density, V_time_nP1, U_time_nP1); /*--- CV volume at time n+1. As we are on a static mesh, the volume of the CV will remained fixed for all time steps. ---*/ @@ -2867,13 +2865,9 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the Jacobian contribution due to the dual time source term. ---*/ if (implicit) { - su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; - - for (iDim = 0; iDim < nDim; iDim++) - Jacobian.AddVal2Diag(iPoint, iDim+1, delta); + su2double delta = (second_order ? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; - if (energy) delta *= Cp; - Jacobian.AddVal2Diag(iPoint, nDim+1, delta); + for (iDim = 0; iDim <= nDim; iDim++) Jacobian.AddVal2Diag(iPoint, iDim + 1, delta); } } END_SU2_OMP_FOR @@ -2898,8 +2892,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = nodes->GetSolution_time_n(iPoint); Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); - V2U(Density, Cp, V_time_n, U_time_n); + V2U(Density, V_time_n, U_time_n); GridVel_i = geometry->nodes->GetGridVel(iPoint); @@ -2954,8 +2947,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = nodes->GetSolution_time_n(iPoint); Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); - V2U(Density, Cp, V_time_n, U_time_n); + V2U(Density, V_time_n, U_time_n); for (iVar = 0; iVar < nVar-!energy; iVar++) LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; @@ -2981,16 +2973,15 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = nodes->GetSolution_time_n(iPoint); V_time_nP1 = nodes->GetSolution(iPoint); - /*--- Access the density and Cp at this node (constant for now). ---*/ + /*--- Access the density at this node (constant for now). ---*/ Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); /*--- Compute the conservative variable vector for all time levels. ---*/ - V2U(Density, Cp, V_time_nM1, U_time_nM1); - V2U(Density, Cp, V_time_n, U_time_n); - V2U(Density, Cp, V_time_nP1, U_time_nP1); + V2U(Density, V_time_nM1, U_time_nM1); + V2U(Density, V_time_n, U_time_n); + V2U(Density, V_time_nP1, U_time_nP1); /*--- CV volume at time n-1 and n+1. In the case of dynamically deforming grids, the volumes will change. On rigidly transforming grids, the @@ -3016,11 +3007,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver if (implicit) { su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; - for (iDim = 0; iDim < nDim; iDim++) + for (iDim = 0; iDim <= nDim; iDim++) Jacobian.AddVal2Diag(iPoint, iDim+1, delta); - - if (energy) delta *= Cp; - Jacobian.AddVal2Diag(iPoint, nDim+1, delta); } } END_SU2_OMP_FOR From 18da912a72576c7ebd4b35d5cfea45381588d768 Mon Sep 17 00:00:00 2001 From: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> Date: Tue, 16 Dec 2025 22:08:33 +0100 Subject: [PATCH 2/4] Update SU2_CFD/src/solvers/CIncEulerSolver.cpp Co-authored-by: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 5b5591177c2..f650b1cd7e3 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -3007,8 +3007,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver if (implicit) { su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; - for (iDim = 0; iDim <= nDim; iDim++) - Jacobian.AddVal2Diag(iPoint, iDim+1, delta); + for (iVar = 1; iVar < nVar; ++iVar) + Jacobian.AddVal2Diag(iPoint, iVar, delta); } } END_SU2_OMP_FOR From a322172532f08bce002ef4977b5378128d5b261a Mon Sep 17 00:00:00 2001 From: Cristopher-Morales Date: Tue, 16 Dec 2025 23:19:35 +0100 Subject: [PATCH 3/4] addressing comment --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index f650b1cd7e3..956af94ab58 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2806,10 +2806,10 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver const bool second_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); const bool energy = config->GetEnergy_Equation(); - const int ndim = nDim; - auto V2U = [ndim](su2double Density, const su2double* V, su2double* U) { + const int nvar = nVar; + auto V2U = [nvar](su2double Density, const su2double* V, su2double* U) { U[0] = Density; - for (int iDim = 0; iDim <= ndim; iDim++) U[iDim + 1] = Density * V[iDim + 1]; + for (int iVar = 1; iVar < nvar; ++iVar) U[iVar] = Density * V[iVar]; }; /*--- Store the physical time step ---*/ @@ -2867,7 +2867,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver if (implicit) { su2double delta = (second_order ? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; - for (iDim = 0; iDim <= nDim; iDim++) Jacobian.AddVal2Diag(iPoint, iDim + 1, delta); + for (iVar = 1; iVar < nVar; ++iVar) Jacobian.AddVal2Diag(iPoint, iVar, delta); } } END_SU2_OMP_FOR From f788ed7c8c6b1ca7dede97e50abf013061620e01 Mon Sep 17 00:00:00 2001 From: Cristopher-Morales Date: Wed, 17 Dec 2025 13:11:29 +0100 Subject: [PATCH 4/4] updating residuals --- TestCases/parallel_regression_AD.py | 4 ++-- TestCases/tutorials.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index b5cbd2eee1d..f5525717d2b 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -316,8 +316,8 @@ def main(): da_unsteadyCHT_cylinder.cfg_dir = "coupled_cht/disc_adj_unsteadyCHT_cylinder" da_unsteadyCHT_cylinder.cfg_file = "chtMaster.cfg" da_unsteadyCHT_cylinder.test_iter = 2 - da_unsteadyCHT_cylinder.test_vals = [-12.131633, -12.617906, -12.688933, -16.179747, -6.432277, 0.000000, 75.761000, 0.247780] - da_unsteadyCHT_cylinder.test_vals_aarch64 = [-12.131633, -12.617906, -12.688933, -16.179747, -6.432277, 0.000000, 75.761000, 0.247780] + da_unsteadyCHT_cylinder.test_vals = [-8.479629, -9.239920, -9.234868, -15.934511, -13.662012, 0.000000, 89.932000, 0.295190] + da_unsteadyCHT_cylinder.test_vals_aarch64 = [-8.479629, -9.239920, -9.234868, -15.934511, -13.662012, 0.000000, 89.932000, 0.295190] da_unsteadyCHT_cylinder.unsteady = True da_unsteadyCHT_cylinder.multizone = True test_list.append(da_unsteadyCHT_cylinder) diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index ae8df575b97..0768bc29370 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -49,7 +49,7 @@ def main(): cht_incompressible_unsteady.cfg_dir = "../Tutorials/multiphysics/unsteady_cht/" cht_incompressible_unsteady.cfg_file = "cht_2d_3cylinders.cfg" cht_incompressible_unsteady.test_iter = 2 - cht_incompressible_unsteady.test_vals = [-2.661443, -2.546289, -0.080399, -0.080399, -0.080399, -12.421980, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 238.240000] #last columns + cht_incompressible_unsteady.test_vals = [-2.661440, -2.534489, -0.080399, -0.080399, -0.080399, -12.421979, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 238.240000] #last columns cht_incompressible_unsteady.multizone = True cht_incompressible_unsteady.unsteady = True test_list.append(cht_incompressible_unsteady)