From 5c1f4f12ceb71a7a28576a7bc69ad385b227d8ed Mon Sep 17 00:00:00 2001 From: Jared Bellon Date: Thu, 3 Sep 2026 13:38:37 -0400 Subject: [PATCH 1/5] Return the distinct left block from the corrected NEMO viscous Jacobian and assemble it subtractively CAvgGradCorrected_NEMO::ComputeResidual returned the j-side Jacobian block in both slots, while its sibling CAvgGrad_NEMO returns the distinct blocks. CNEMONSSolver::Viscous_Residual subtracts the viscous residual at iPoint and adds it at jPoint but assembled the Jacobian with UpdateBlocks, the additive convention. With NEMO's viscous Jacobians being flux derivatives (dF/dU_i = -dF/dU_j), UpdateBlocksSub(J_i, J_j) assembles exactly dR/dU in all four blocks; the two hunks are one fix and stay in one commit. The residual is unchanged, so converged solutions do not change; the implicit operator does. Add the distinct-block unit test on a 2-D AIR-5 edge (NEMOViscousFixture) and an assembly-level test that builds a NEMO Navier-Stokes solver on a box mesh, calls Viscous_Residual and checks that the assembled residual and Jacobian blocks equal the derivative of the subtracted and added edge fluxes; both fail on develop and pass with this change. Rebaseline the regression vectors that move as a consequence: visc_cone (serial and parallel), super_cat and ion_gy (parallel). Add a regression case, ion_gy_march, that marches the ionized cylinder 100 iterations from its stored restart with the residual gate disabled: develop diverges (energy residual above 8 and a wall heat flux of order -1e10 by iteration 99) while the corrected operator holds the restart state. Give the two cylinder configurations distinct output restart names so that neither overwrites the restart the other reads. Co-Authored-By: Claude Fable 5.1 --- SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp | 2 +- SU2_CFD/src/solvers/CNEMONSSolver.cpp | 2 +- .../visc_cylinder/cyl_ion_gy.cfg | 2 +- .../visc_cylinder/cyl_ion_gy_march.cfg | 80 ++++++ TestCases/parallel_regression.py | 16 +- TestCases/serial_regression.py | 2 +- UnitTests/SU2_CFD/nemo_viscous_assembly.cpp | 246 ++++++++++++++++++ .../SU2_CFD/numerics/CNumerics_tests.cpp | 131 ++++++++++ UnitTests/meson.build | 1 + 9 files changed, 475 insertions(+), 7 deletions(-) create mode 100644 TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg create mode 100644 UnitTests/SU2_CFD/nemo_viscous_assembly.cpp diff --git a/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp b/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp index 2553fa2915d3..705c1b253adb 100644 --- a/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp +++ b/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp @@ -359,6 +359,6 @@ CNumerics::ResidualType<> CAvgGradCorrected_NEMO::ComputeResidual(const CConfig } - return ResidualType<>(Flux, Jacobian_j, Jacobian_j); + return ResidualType<>(Flux, Jacobian_i, Jacobian_j); } diff --git a/SU2_CFD/src/solvers/CNEMONSSolver.cpp b/SU2_CFD/src/solvers/CNEMONSSolver.cpp index bae5591d9460..e29934e88032 100644 --- a/SU2_CFD/src/solvers/CNEMONSSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMONSSolver.cpp @@ -238,7 +238,7 @@ void CNEMONSSolver::Viscous_Residual(CGeometry *geometry, LinSysRes.SubtractBlock(iPoint, residual); LinSysRes.AddBlock(jPoint, residual); if (implicit) { - Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); + Jacobian.UpdateBlocksSub(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j); } } } //iEdge diff --git a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg index 2f51ad7f73fc..37b682343194 100644 --- a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg +++ b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg @@ -65,7 +65,7 @@ CONV_STARTITER= 10 MESH_FILENAME= visc_cyl.su2 MESH_FORMAT= SU2 SOLUTION_FILENAME= restart_flow_gy -RESTART_FILENAME= restart_flow_gy +RESTART_FILENAME= restart_flow_gy_out TABULAR_FORMAT= TECPLOT CONV_FILENAME= convergence VOLUME_FILENAME= soln_volume diff --git a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg new file mode 100644 index 000000000000..04128bfd34eb --- /dev/null +++ b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg @@ -0,0 +1,80 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Case: ionized AIR-7 viscous cylinder, restarted from the stored solution and %% +%% marched 100 iterations with the residual gate disabled. Checks that the %% +%% implicit NEMO viscous operator holds the restart state instead of diverging. %% +%% Companion to cyl_ion_gy.cfg (same mesh and restart). %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Mach 15 viscous, ionized flow over a cylinder % +% using the Gupta-Yos transport model. % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= NEMO_NAVIER_STOKES +MATH_PROBLEM= DIRECT +RESTART_SOL= YES +AXISYMMETRIC= YES + +% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% +% +MACH_NUMBER=15 +AOA= 0.0 +SIDESLIP_ANGLE= 0.0 +FREESTREAM_PRESSURE= 21.96 +FREESTREAM_TEMPERATURE= 246.87 +FREESTREAM_TEMPERATURE_VE= 246.87 +INIT_OPTION=TD_CONDITIONS + +% ---- NONEQUILIBRIUM GAS, IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% +% +FLUID_MODEL= SU2_NONEQ +GAS_MODEL= AIR-7 +GAS_COMPOSITION= (5E-13, 0.767, 0.23, 0.001, 0.001, 0.000999999999, 5E-13) +TRANSPORT_COEFF_MODEL= GUPTA-YOS + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_ISOTHERMAL=(wall, 300.0) +MARKER_FAR= (far, exit) +MARKER_SYM= (symm) +MARKER_PLOTTING= (wall) +MARKER_MONITORING= (wall) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 3.0 +ITER= 100 +LINEAR_SOLVER= BCGSTAB +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 5 + +% -----------------------------------------------------------------------% +% +CONV_NUM_METHOD_FLOW= AUSM +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= VENKATAKRISHNAN +VENKAT_LIMITER_COEFF= 0.05 +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -30 +CONV_STARTITER= 10 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= visc_cyl.su2 +MESH_FORMAT= SU2 +SOLUTION_FILENAME= restart_flow_gy +RESTART_FILENAME= restart_flow_gy_march +TABULAR_FORMAT= TECPLOT +CONV_FILENAME= convergence +VOLUME_FILENAME= soln_volume +SURFACE_FILENAME= soln_surface +OUTPUT_WRT_FREQ= 100000 +SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY_0, RMS_DENSITY_1, RMS_DENSITY_2, RMS_DENSITY_3, RMS_DENSITY_4, RMS_DENSITY_5, RMS_DENSITY_6, RMS_ENERGY, RMS_ENERGY_VE, LIFT, DRAG, TOTAL_HEATFLUX diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ff9b8c667061..fd695f4b0df7 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -192,7 +192,7 @@ def main(): visc_cone.cfg_dir = "nonequilibrium/visc_wedge" visc_cone.cfg_file = "axi_visccone.cfg" visc_cone.test_iter = 10 - visc_cone.test_vals = [-5.215234, -5.739371, -20.559852, -20.509281, -20.408911, 1.262701, -3.205457, -0.015696, 0.093205, 32637.000000] + visc_cone.test_vals = [-5.298545, -5.823026, -20.404788, -20.318330, -20.378801, 1.067473, -3.250998, -0.015489, 0.095136, 24939.000000] visc_cone.test_vals_aarch64 = [-5.222270, -5.746525, -20.560286, -20.510152, -20.409101, 1.255758, -3.208382, -0.016014, 0.093462, 32619.000000] test_list.append(visc_cone) @@ -209,7 +209,7 @@ def main(): super_cat.cfg_dir = "nonequilibrium/visc_wedge" super_cat.cfg_file = "super_cat.cfg" super_cat.test_iter = 10 - super_cat.test_vals = [-5.232595, -5.757889, -20.641415, -20.640623, -20.541670, 1.246866, -3.205258, -0.028372, 0.250647, 32440.000000] + super_cat.test_vals = [-5.309257, -5.834048, -21.098287, -21.157699, -21.180688, 1.056908, -3.252349, -0.028039, 0.252019, 24878.000000] test_list.append(super_cat) # Viscous single wedge - partially catalytic walls @@ -225,9 +225,19 @@ def main(): ion_gy.cfg_dir = "nonequilibrium/visc_cylinder" ion_gy.cfg_file = "cyl_ion_gy.cfg" ion_gy.test_iter = 10 - ion_gy.test_vals = [-11.629873, -4.165562, -4.702662, -4.950351, -5.146155, -4.993878, -6.893332, 5.990109, 5.990004, -0.014849, 0.000000, 90090.000000] + ion_gy.test_vals = [-12.682344, -4.169796, -4.714319, -5.658742, -5.712605, -5.269125, -7.961372, 3.359972, 1.705292, -0.014849, 0.000000, 90286.000000] test_list.append(ion_gy) + # Ionized cylinder marched from its restart with the residual gate disabled: the + # corrected NEMO viscous Jacobian holds the restart state, the previous one diverged. + ion_gy_march = TestCase('ion_gy_march') + ion_gy_march.cfg_dir = "nonequilibrium/visc_cylinder" + ion_gy_march.cfg_file = "cyl_ion_gy_march.cfg" + ion_gy_march.test_iter = 99 + ion_gy_march.test_vals = [-11.662039, -4.203178, -4.868257, -5.462497, -5.232052, -4.960881, -6.951391, 4.541901, 4.552855, -0.014861, 0.000001, 90357.000000] + ion_gy_march.tol = 0.01 + test_list.append(ion_gy_march) + ########################## ### Compressible Euler ### ########################## diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index aa22e30d8ae9..0dd13519c9c3 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -75,7 +75,7 @@ def main(): visc_cone.cfg_dir = "nonequilibrium/visc_wedge" visc_cone.cfg_file = "axi_visccone.cfg" visc_cone.test_iter = 10 - visc_cone.test_vals = [-5.215230, -5.739367, -20.560781, -20.516922, -20.406516, 1.262782, -3.205476, -0.015696, 0.093206, 32641] + visc_cone.test_vals = [-5.298535, -5.823020, -20.330692, -20.368713, -20.310533, 1.067483, -3.250989, -0.015487, 0.095136, 24939] visc_cone.test_vals_aarch64 = [-5.215250, -5.739384, -20.560917, -20.517096, -20.406630, 1.262772, -3.205492, -0.015695, 0.093205, 32641.000000] test_list.append(visc_cone) diff --git a/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp b/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp new file mode 100644 index 000000000000..0fd184ff6a3a --- /dev/null +++ b/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp @@ -0,0 +1,246 @@ +/*! + * \file nemo_viscous_assembly.cpp + * \brief Unit test for the edge assembly of the NEMO viscous residual and Jacobian. + * \author J. Bellon + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include +#include +#include +#include +#include +#include "../../Common/include/geometry/CPhysicalGeometry.hpp" +#include "../../SU2_CFD/include/numerics/NEMO/NEMO_diffusion.hpp" +#include "../../SU2_CFD/include/solvers/CNEMONSSolver.hpp" +#include "../../SU2_CFD/include/variables/CNEMONSVariable.hpp" + +namespace { + +/*! + * \brief A NEMO Navier-Stokes solver on a small box mesh, built the same way + * the gradient tests build their geometry, with a non-uniform velocity field + * so that the viscous residual is not trivially zero. + */ +struct NEMOViscousAssemblyCase { + const std::string configOptions = + "SOLVER= NEMO_NAVIER_STOKES\n" + "GAS_MODEL= AIR-5\n" + "GAS_COMPOSITION= (0.77, 0.23, 0.0, 0.0, 0.0)\n" + "FLUID_MODEL= SU2_NONEQ\n" + "MESH_FORMAT= BOX\n" + "MESH_BOX_SIZE= 4,4,4\n" + "MESH_BOX_LENGTH= 1,1,1\n" + "MESH_BOX_OFFSET= 0,0,0\n" + "INIT_OPTION= TD_CONDITIONS\n" + "MACH_NUMBER= 0.3\n" + "FREESTREAM_PRESSURE= 1000.0\n" + "FREESTREAM_TEMPERATURE= 300.0\n" + "FREESTREAM_TEMPERATURE_VE= 300.0\n" + "REYNOLDS_NUMBER= 1000\n" + "KIND_TURB_MODEL= NONE\n" + "MARKER_FAR= (x_minus, x_plus, y_minus, y_plus, z_minus, z_plus)\n" + "NUM_METHOD_GRAD= GREEN_GAUSS\n" + "CONV_NUM_METHOD_FLOW= AUSM\n" + "MUSCL_FLOW= NO\n" + "TIME_DISCRE_FLOW= EULER_IMPLICIT\n"; + + std::unique_ptr config; + std::unique_ptr geometry; + CNEMONSSolver* solver{nullptr}; + + NEMOViscousAssemblyCase() { + auto origBuf = cout.rdbuf(); + cout.rdbuf(nullptr); + + stringstream ss(configOptions); + config = std::unique_ptr(new CConfig(ss, SU2_COMPONENT::SU2_CFD, false)); + + { + auto aux_geometry = std::unique_ptr(new CPhysicalGeometry(config.get(), 0, 1)); + geometry = std::unique_ptr(new CPhysicalGeometry(aux_geometry.get(), config.get())); + } + geometry->SetSendReceive(config.get()); + geometry->SetBoundaries(config.get()); + geometry->SetPoint_Connectivity(); + geometry->SetElement_Connectivity(); + geometry->SetBoundVolume(); + geometry->Check_IntElem_Orientation(config.get()); + geometry->Check_BoundElem_Orientation(config.get()); + geometry->SetEdges(); + geometry->SetVertex(config.get()); + geometry->SetControlVolume(config.get(), ALLOCATE); + geometry->SetBoundControlVolume(config.get(), ALLOCATE); + geometry->FindNormal_Neighbor(config.get()); + geometry->SetGlobal_to_Local_Point(); + geometry->PreprocessP2PComms(geometry.get(), config.get()); + + solver = new CNEMONSSolver(geometry.get(), config.get(), MESH_0); + + /*--- Scale the momentum with position so that the velocity gradients, + * and with them the viscous fluxes, are non-zero. ---*/ + const auto nDim = geometry->GetnDim(); + const auto nSpecies = config->GetnSpecies(); + auto* nodes = solver->GetNodes(); + for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) { + const auto* coord = geometry->nodes->GetCoord(iPoint); + const su2double scale = 1.0 + 0.2 * coord[0] + 0.1 * coord[1]; + for (unsigned short iDim = 0; iDim < nDim; ++iDim) { + const auto iVar = nSpecies + iDim; + nodes->SetSolution(iPoint, iVar, scale * nodes->GetSolution(iPoint, iVar)); + } + } + + cout.rdbuf(origBuf); + } + + ~NEMOViscousAssemblyCase() { delete solver; } +}; + +} // namespace + +TEST_CASE("NEMO viscous edge assembly is the derivative of the subtracted and added residual", + "[NEMO][viscous][Jacobian]") { + NEMOViscousAssemblyCase testCase; + auto* config = testCase.config.get(); + auto* geometry = testCase.geometry.get(); + auto* solver = testCase.solver; + + const auto nDim = geometry->GetnDim(); + const auto nVar = solver->GetnVar(); + const auto nPrimVar = solver->GetnPrimVar(); + const auto nPrimVarGrad = solver->GetnPrimVarGrad(); + const auto nPoint = geometry->GetnPoint(); + + /*--- Primitive variables, transport properties and gradients. ---*/ + CSolver* solver_container[MAX_SOLS] = {nullptr}; + solver_container[FLOW_SOL] = solver; + { + auto origBuf = cout.rdbuf(); + cout.rdbuf(nullptr); + solver->Preprocessing(geometry, solver_container, config, MESH_0, 0, RUNTIME_FLOW_SYS, false); + cout.rdbuf(origBuf); + } + solver->LinSysRes.SetValZero(); + solver->Jacobian.SetValZero(); + + /*--- Assemble the viscous residual and Jacobian with the solver. The + * solver's override is private, so dispatch through the public base + * interface, exactly as the integration classes do. ---*/ + CAvgGradCorrected_NEMO numerics(nDim, nVar, nPrimVar, nPrimVarGrad, config); + CNumerics* numerics_container[MAX_TERMS] = {nullptr}; + numerics_container[VISC_TERM] = &numerics; + CSolver& base = *solver; + base.Viscous_Residual(geometry, solver_container, numerics_container, config, MESH_0, 0); + + /*--- Replay every edge with an independent numerics object. The residual + * contribution F of an edge is subtracted at i and added at j, so the + * derivative of the assembled residual is -dF/dU in the i rows and +dF/dU + * in the j rows, with dF/dU_i and dF/dU_j the two blocks returned by the + * numerics. Accumulate that reference and compare it with what the solver + * assembled. ---*/ + CAvgGradCorrected_NEMO replay(nDim, nVar, nPrimVar, nPrimVarGrad, config); + auto* nodes = dynamic_cast(solver->GetNodes()); + REQUIRE(nodes != nullptr); + + std::vector residual_ref(nPoint * nVar, 0.0); + std::vector diagonal_ref(nPoint * nVar * nVar, 0.0); + su2double offdiag_error = 0.0, jacobian_scale = 0.0; + + for (auto iEdge = 0ul; iEdge < geometry->GetnEdge(); ++iEdge) { + const auto iPoint = geometry->edges->GetNode(iEdge, 0); + const auto jPoint = geometry->edges->GetNode(iEdge, 1); + + replay.SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint)); + replay.SetNormal(geometry->edges->GetNormal(iEdge)); + replay.SetConservative(nodes->GetSolution(iPoint), nodes->GetSolution(jPoint)); + replay.SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(jPoint)); + replay.SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(jPoint)); + replay.SetdPdU(nodes->GetdPdU(iPoint), nodes->GetdPdU(jPoint)); + replay.SetdTdU(nodes->GetdTdU(iPoint), nodes->GetdTdU(jPoint)); + replay.SetdTvedU(nodes->GetdTvedU(iPoint), nodes->GetdTvedU(jPoint)); + replay.SetEve(nodes->GetEve(iPoint), nodes->GetEve(jPoint)); + replay.SetCvve(nodes->GetCvve(iPoint), nodes->GetCvve(jPoint)); + replay.SetDiffusionCoeff(nodes->GetDiffusionCoeff(iPoint), nodes->GetDiffusionCoeff(jPoint)); + replay.SetLaminarViscosity(nodes->GetLaminarViscosity(iPoint), nodes->GetLaminarViscosity(jPoint)); + replay.SetEddyViscosity(nodes->GetEddyViscosity(iPoint), nodes->GetEddyViscosity(jPoint)); + replay.SetThermalConductivity(nodes->GetThermalConductivity(iPoint), nodes->GetThermalConductivity(jPoint)); + replay.SetThermalConductivity_ve(nodes->GetThermalConductivity_ve(iPoint), + nodes->GetThermalConductivity_ve(jPoint)); + + const auto edge = replay.ComputeResidual(config); + + const auto* block_ij = solver->Jacobian.GetBlock(iPoint, jPoint); + const auto* block_ji = solver->Jacobian.GetBlock(jPoint, iPoint); + REQUIRE(block_ij != nullptr); + REQUIRE(block_ji != nullptr); + + for (unsigned short iVar = 0; iVar < nVar; ++iVar) { + residual_ref[iPoint * nVar + iVar] -= edge.residual[iVar]; + residual_ref[jPoint * nVar + iVar] += edge.residual[iVar]; + + for (unsigned short jVar = 0; jVar < nVar; ++jVar) { + const su2double dFdUi = edge.jacobian_i[iVar][jVar]; + const su2double dFdUj = edge.jacobian_j[iVar][jVar]; + jacobian_scale = std::max(jacobian_scale, std::max(std::fabs(dFdUi), std::fabs(dFdUj))); + + /*--- Diagonal blocks collect every edge of a point. ---*/ + diagonal_ref[(iPoint * nVar + iVar) * nVar + jVar] -= dFdUi; + diagonal_ref[(jPoint * nVar + iVar) * nVar + jVar] += dFdUj; + + /*--- Off-diagonal blocks belong to this edge alone. ---*/ + const su2double assembled_ij = block_ij[iVar * nVar + jVar]; + const su2double assembled_ji = block_ji[iVar * nVar + jVar]; + offdiag_error = std::max(offdiag_error, std::fabs(assembled_ij - (-dFdUj))); + offdiag_error = std::max(offdiag_error, std::fabs(assembled_ji - (+dFdUi))); + } + } + } + + su2double diagonal_error = 0.0, residual_error = 0.0, residual_scale = 0.0; + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + const auto* block_ii = solver->Jacobian.GetBlock(iPoint, iPoint); + const auto* assembled_residual = solver->LinSysRes.GetBlock(iPoint); + REQUIRE(block_ii != nullptr); + for (unsigned short iVar = 0; iVar < nVar; ++iVar) { + residual_scale = std::max(residual_scale, std::fabs(residual_ref[iPoint * nVar + iVar])); + residual_error = + std::max(residual_error, std::fabs(assembled_residual[iVar] - residual_ref[iPoint * nVar + iVar])); + for (unsigned short jVar = 0; jVar < nVar; ++jVar) { + const su2double assembled = block_ii[iVar * nVar + jVar]; + diagonal_error = + std::max(diagonal_error, std::fabs(assembled - diagonal_ref[(iPoint * nVar + iVar) * nVar + jVar])); + } + } + } + + /*--- The reference must be non-trivial, otherwise the sign is untested. ---*/ + REQUIRE(jacobian_scale > 0.0); + REQUIRE(residual_scale > 0.0); + + const su2double tolerance = 1.0e-10; + CHECK(residual_error <= tolerance * residual_scale); + CHECK(offdiag_error <= tolerance * jacobian_scale); + CHECK(diagonal_error <= tolerance * jacobian_scale); +} diff --git a/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp index a4cee41e8d20..d9a32e937539 100644 --- a/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp +++ b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp @@ -26,8 +26,11 @@ */ #include "catch.hpp" +#include #include +#include #include "../../../SU2_CFD/include/numerics/CNumerics.hpp" +#include "../../../SU2_CFD/include/numerics/NEMO/NEMO_diffusion.hpp" TEST_CASE("NTS blending has a minimum of 0.05", "[Upwind/central blending]") { std::stringstream config_options; @@ -89,3 +92,131 @@ TEST_CASE("QCR2000 corrects only the turbulent stress", "[QCR]") { for (size_t jDim = 0; jDim < nDim; jDim++) REQUIRE(tau_total[iDim][jDim] == Approx(tau_lam[iDim][jDim] + tau_turb[iDim][jDim]).margin(1e-12)); } + +namespace { + +struct NEMOViscousFixture { + static constexpr unsigned short nDim = 2; + static constexpr unsigned short nSpecies = 5; + static constexpr unsigned short nVar = nSpecies + nDim + 2; + static constexpr unsigned short nPrimVar = nSpecies + nDim + 10; + static constexpr unsigned short nPrimVarGrad = nSpecies + nDim + 8; + static constexpr unsigned short T_INDEX = nSpecies; + static constexpr unsigned short TVE_INDEX = nSpecies + 1; + static constexpr unsigned short VEL_INDEX = nSpecies + 2; + static constexpr unsigned short P_INDEX = nSpecies + nDim + 2; + static constexpr unsigned short RHO_INDEX = nSpecies + nDim + 3; + static constexpr unsigned short H_INDEX = nSpecies + nDim + 4; + static constexpr unsigned short A_INDEX = nSpecies + nDim + 5; + static constexpr unsigned short RHOCVTR_INDEX = nSpecies + nDim + 6; + static constexpr unsigned short RHOCVVE_INDEX = nSpecies + nDim + 7; + + std::stringstream options; + CConfig* config = nullptr; + std::vector primitive_i = std::vector(nPrimVar, 0.0); + std::vector primitive_j = std::vector(nPrimVar, 0.0); + su2activematrix gradient_i = su2activematrix(nPrimVarGrad, nDim); + su2activematrix gradient_j = su2activematrix(nPrimVarGrad, nDim); + std::array diffusion_i{}; + std::array diffusion_j{}; + std::array eve_i{}; + std::array eve_j{}; + std::array cvve_i{}; + std::array cvve_j{}; + std::array dT_i{}; + std::array dT_j{}; + std::array dTve_i{}; + std::array dTve_j{}; + const std::array coord_i{0.0, 0.0}; + const std::array coord_j{3.0, 4.0}; + const std::array normal{0.6, 0.8}; + + NEMOViscousFixture() { + options << "SOLVER= NEMO_NAVIER_STOKES\n" + << "GAS_MODEL= AIR-5\n" + << "GAS_COMPOSITION= (0.77, 0.23, 0.0, 0.0, 0.0)\n" + << "FLUID_MODEL= SU2_NONEQ\n" + << "FROZEN_MIXTURE= YES\n" + << "TIME_DISCRE_FLOW= EULER_IMPLICIT\n" + << "CONV_NUM_METHOD_FLOW= AUSM\n"; + config = new CConfig(options, SU2_COMPONENT::SU2_CFD, false); + + primitive_i[0] = primitive_j[0] = 0.77; + primitive_i[1] = primitive_j[1] = 0.23; + primitive_i[T_INDEX] = primitive_j[T_INDEX] = 300.0; + primitive_i[TVE_INDEX] = primitive_j[TVE_INDEX] = 300.0; + primitive_i[P_INDEX] = primitive_j[P_INDEX] = 101325.0; + primitive_i[RHO_INDEX] = primitive_j[RHO_INDEX] = 1.0; + primitive_i[H_INDEX] = primitive_j[H_INDEX] = 3.0e5; + primitive_i[A_INDEX] = primitive_j[A_INDEX] = 340.0; + primitive_i[RHOCVTR_INDEX] = primitive_j[RHOCVTR_INDEX] = 700.0; + primitive_i[RHOCVVE_INDEX] = primitive_j[RHOCVVE_INDEX] = 1.0; + gradient_i = su2double(0.0); + gradient_j = su2double(0.0); + } + + ~NEMOViscousFixture() { delete config; } + + void set_common(CNumerics& numerics) { + numerics.SetCoord(coord_i.data(), coord_j.data()); + numerics.SetNormal(normal.data()); + numerics.SetPrimitive(primitive_i.data(), primitive_j.data()); + numerics.SetPrimVarGradient(CMatrixView(gradient_i), CMatrixView(gradient_j)); + numerics.SetDiffusionCoeff(diffusion_i.data(), diffusion_j.data()); + numerics.SetLaminarViscosity(2.0, 2.0); + numerics.SetEddyViscosity(0.0, 0.0); + numerics.SetThermalConductivity(0.0, 0.0); + numerics.SetThermalConductivity_ve(0.0, 0.0); + numerics.SetEve(eve_i.data(), eve_j.data()); + numerics.SetCvve(cvve_i.data(), cvve_j.data()); + numerics.SetdTdU(dT_i.data(), dT_j.data()); + numerics.SetdTvedU(dTve_i.data(), dTve_j.data()); + } + + void set_edge_linear_velocity_gradient() { + gradient_i = su2double(0.0); + gradient_j = su2double(0.0); + constexpr su2double distance_squared = 25.0; + for (unsigned short component = 0; component < nDim; ++component) { + const su2double jump = primitive_j[VEL_INDEX + component] - primitive_i[VEL_INDEX + component]; + for (unsigned short dimension = 0; dimension < nDim; ++dimension) { + const su2double edge_component = coord_j[dimension] - coord_i[dimension]; + gradient_i(VEL_INDEX + component, dimension) = jump * edge_component / distance_squared; + gradient_j(VEL_INDEX + component, dimension) = jump * edge_component / distance_squared; + } + } + } + + template + su2double directional_flux(Numerics& numerics) { + set_common(numerics); + const auto residual = numerics.ComputeResidual(config); + su2double projected = 0.0; + for (unsigned short component = 0; component < nDim; ++component) + projected += normal[component] * residual.residual[nSpecies + component]; + return projected; + } + + static su2double directional_jacobian(const su2double* const* matrix, const std::array& direction) { + su2double value = 0.0; + for (unsigned short row = 0; row < nDim; ++row) + for (unsigned short column = 0; column < nDim; ++column) + value += direction[row] * matrix[nSpecies + row][nSpecies + column] * direction[column]; + return value; + } +}; + +} // namespace + +TEST_CASE("NEMO corrected viscous residual returns distinct i and j Jacobians", "[NEMO][viscous][Jacobian]") { + NEMOViscousFixture fixture; + CAvgGradCorrected_NEMO numerics(fixture.nDim, fixture.nVar, fixture.nPrimVar, fixture.nPrimVarGrad, fixture.config); + fixture.set_common(numerics); + const auto base = numerics.ComputeResidual(fixture.config); + const su2double analytic_i = fixture.directional_jacobian(base.jacobian_i, fixture.normal); + const su2double analytic_j = fixture.directional_jacobian(base.jacobian_j, fixture.normal); + + REQUIRE(analytic_i == Approx(-8.0 / 15.0).epsilon(1.0e-12)); + REQUIRE(analytic_j == Approx(8.0 / 15.0).epsilon(1.0e-12)); + REQUIRE(analytic_i != Approx(analytic_j)); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index f8c22511b5f8..66917b9782be 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -16,6 +16,7 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', 'SU2_CFD/gradients.cpp', + 'SU2_CFD/nemo_viscous_assembly.cpp', 'SU2_CFD/windowing.cpp', 'Common/toolboxes/random_toolbox_tests.cpp', 'Common/linear_algebra/quantization_tests.cpp']) From 077f73e0654e705101a12b9817d332ddddeb6e66 Mon Sep 17 00:00:00 2001 From: Jared Bellon Date: Sat, 5 Sep 2026 16:41:46 -0400 Subject: [PATCH 2/5] Remove trailing whitespace from ion cylinder march config Keep git diff --check clean before upstream PR submission. --- TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg index 04128bfd34eb..09d9eaa50009 100644 --- a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg +++ b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg @@ -70,11 +70,11 @@ CONV_STARTITER= 10 % MESH_FILENAME= visc_cyl.su2 MESH_FORMAT= SU2 -SOLUTION_FILENAME= restart_flow_gy +SOLUTION_FILENAME= restart_flow_gy RESTART_FILENAME= restart_flow_gy_march TABULAR_FORMAT= TECPLOT CONV_FILENAME= convergence VOLUME_FILENAME= soln_volume -SURFACE_FILENAME= soln_surface +SURFACE_FILENAME= soln_surface OUTPUT_WRT_FREQ= 100000 SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY_0, RMS_DENSITY_1, RMS_DENSITY_2, RMS_DENSITY_3, RMS_DENSITY_4, RMS_DENSITY_5, RMS_DENSITY_6, RMS_ENERGY, RMS_ENERGY_VE, LIFT, DRAG, TOTAL_HEATFLUX From 8bfba2f3afd41c4ff2a7e71fb61a464e01744f64 Mon Sep 17 00:00:00 2001 From: BlueChips Date: Sat, 5 Sep 2026 20:10:58 -0400 Subject: [PATCH 3/5] Refresh native ARM64 NEMO regression vectors Update the serial and two-rank viscous-cone aarch64 expectations from GitHub Actions run 33998938091 at exact pre-update head 077f73e0654e705101a12b9817d332ddddeb6e66. Artifact digest: sha256:e7eedd1985e4e87830bdb740aa1b2a02e2f5047f0afe1a18bb88d8dce9d0b243. No x86 expectation or solver source changes. --- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index fd695f4b0df7..fcc56c2babe8 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -193,7 +193,7 @@ def main(): visc_cone.cfg_file = "axi_visccone.cfg" visc_cone.test_iter = 10 visc_cone.test_vals = [-5.298545, -5.823026, -20.404788, -20.318330, -20.378801, 1.067473, -3.250998, -0.015489, 0.095136, 24939.000000] - visc_cone.test_vals_aarch64 = [-5.222270, -5.746525, -20.560286, -20.510152, -20.409101, 1.255758, -3.208382, -0.016014, 0.093462, 32619.000000] + visc_cone.test_vals_aarch64 = [-5.298530, -5.823015, -20.404788, -20.318329, -20.378801, 1.067489, -3.250987, -0.015488, 0.095136, 24939.000000] test_list.append(visc_cone) # Viscous single wedge with Mutation++ diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 0dd13519c9c3..bc67350f59cd 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -76,7 +76,7 @@ def main(): visc_cone.cfg_file = "axi_visccone.cfg" visc_cone.test_iter = 10 visc_cone.test_vals = [-5.298535, -5.823020, -20.330692, -20.368713, -20.310533, 1.067483, -3.250989, -0.015487, 0.095136, 24939] - visc_cone.test_vals_aarch64 = [-5.215250, -5.739384, -20.560917, -20.517096, -20.406630, 1.262772, -3.205492, -0.015695, 0.093205, 32641.000000] + visc_cone.test_vals_aarch64 = [-5.298550, -5.823031, -20.330692, -20.368713, -20.310533, 1.067466, -3.251001, -0.015489, 0.095136, 24939.000000] test_list.append(visc_cone) ######################### From cd8a39d097ef15251318da727b4cf950d41cef3f Mon Sep 17 00:00:00 2001 From: BlueChips Date: Mon, 7 Sep 2026 15:41:24 -0400 Subject: [PATCH 4/5] Refine NEMO viscous Jacobian regression tests Reuse ion_gy for 100 iterations and remove the duplicate configuration. Check the momentum flux derivatives independently and clarify the assembly test coverage. Validated the focused tests and both restored-bug controls. --- .../visc_cylinder/cyl_ion_gy.cfg | 2 +- .../visc_cylinder/cyl_ion_gy_march.cfg | 80 ------------------- TestCases/parallel_regression.py | 17 +--- UnitTests/SU2_CFD/nemo_viscous_assembly.cpp | 22 +++-- .../SU2_CFD/numerics/CNumerics_tests.cpp | 51 ++++++++---- 5 files changed, 57 insertions(+), 115 deletions(-) delete mode 100644 TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg diff --git a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg index 37b682343194..f7ae13236606 100644 --- a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg +++ b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg @@ -57,7 +57,7 @@ TIME_DISCRE_FLOW= EULER_IMPLICIT % --------------------------- CONVERGENCE PARAMETERS --------------------------% % -CONV_RESIDUAL_MINVAL= -10 +CONV_RESIDUAL_MINVAL= -30 CONV_STARTITER= 10 % ------------------------- INPUT/OUTPUT INFORMATION --------------------------% diff --git a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg b/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg deleted file mode 100644 index 09d9eaa50009..000000000000 --- a/TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy_march.cfg +++ /dev/null @@ -1,80 +0,0 @@ -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%% Case: ionized AIR-7 viscous cylinder, restarted from the stored solution and %% -%% marched 100 iterations with the residual gate disabled. Checks that the %% -%% implicit NEMO viscous operator holds the restart state instead of diverging. %% -%% Companion to cyl_ion_gy.cfg (same mesh and restart). %% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% % -% SU2 configuration file % -% Case description: Mach 15 viscous, ionized flow over a cylinder % -% using the Gupta-Yos transport model. % -% % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% -% -SOLVER= NEMO_NAVIER_STOKES -MATH_PROBLEM= DIRECT -RESTART_SOL= YES -AXISYMMETRIC= YES - -% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% -% -MACH_NUMBER=15 -AOA= 0.0 -SIDESLIP_ANGLE= 0.0 -FREESTREAM_PRESSURE= 21.96 -FREESTREAM_TEMPERATURE= 246.87 -FREESTREAM_TEMPERATURE_VE= 246.87 -INIT_OPTION=TD_CONDITIONS - -% ---- NONEQUILIBRIUM GAS, IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% -% -FLUID_MODEL= SU2_NONEQ -GAS_MODEL= AIR-7 -GAS_COMPOSITION= (5E-13, 0.767, 0.23, 0.001, 0.001, 0.000999999999, 5E-13) -TRANSPORT_COEFF_MODEL= GUPTA-YOS - -% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% -% -MARKER_ISOTHERMAL=(wall, 300.0) -MARKER_FAR= (far, exit) -MARKER_SYM= (symm) -MARKER_PLOTTING= (wall) -MARKER_MONITORING= (wall) - -% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% -% -NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -CFL_NUMBER= 3.0 -ITER= 100 -LINEAR_SOLVER= BCGSTAB -LINEAR_SOLVER_ERROR= 1E-6 -LINEAR_SOLVER_ITER= 5 - -% -----------------------------------------------------------------------% -% -CONV_NUM_METHOD_FLOW= AUSM -MUSCL_FLOW= YES -SLOPE_LIMITER_FLOW= VENKATAKRISHNAN -VENKAT_LIMITER_COEFF= 0.05 -TIME_DISCRE_FLOW= EULER_IMPLICIT - -% --------------------------- CONVERGENCE PARAMETERS --------------------------% -% -CONV_RESIDUAL_MINVAL= -30 -CONV_STARTITER= 10 - -% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% -% -MESH_FILENAME= visc_cyl.su2 -MESH_FORMAT= SU2 -SOLUTION_FILENAME= restart_flow_gy -RESTART_FILENAME= restart_flow_gy_march -TABULAR_FORMAT= TECPLOT -CONV_FILENAME= convergence -VOLUME_FILENAME= soln_volume -SURFACE_FILENAME= soln_surface -OUTPUT_WRT_FREQ= 100000 -SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY_0, RMS_DENSITY_1, RMS_DENSITY_2, RMS_DENSITY_3, RMS_DENSITY_4, RMS_DENSITY_5, RMS_DENSITY_6, RMS_ENERGY, RMS_ENERGY_VE, LIFT, DRAG, TOTAL_HEATFLUX diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index fcc56c2babe8..89f2fdaa9792 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -220,24 +220,15 @@ def main(): partial_cat.test_vals = [-5.210302, -5.735065, -20.880448, -20.825971, -23.475263, 1.806201, -2.813952, -0.078400, 0.495606, 29020.000000] test_list.append(partial_cat) - # Viscous cylinder, ionization, Gupta-Yos + # Viscous cylinder, ionization, Gupta-Yos, marched 100 iterations from its restart. ion_gy = TestCase('ion_gy') ion_gy.cfg_dir = "nonequilibrium/visc_cylinder" ion_gy.cfg_file = "cyl_ion_gy.cfg" - ion_gy.test_iter = 10 - ion_gy.test_vals = [-12.682344, -4.169796, -4.714319, -5.658742, -5.712605, -5.269125, -7.961372, 3.359972, 1.705292, -0.014849, 0.000000, 90286.000000] + ion_gy.test_iter = 99 + ion_gy.test_vals = [-11.662039, -4.203178, -4.868257, -5.462497, -5.232052, -4.960881, -6.951391, 4.541901, 4.552855, -0.014861, 0.000001, 90357.000000] + ion_gy.tol = 0.01 test_list.append(ion_gy) - # Ionized cylinder marched from its restart with the residual gate disabled: the - # corrected NEMO viscous Jacobian holds the restart state, the previous one diverged. - ion_gy_march = TestCase('ion_gy_march') - ion_gy_march.cfg_dir = "nonequilibrium/visc_cylinder" - ion_gy_march.cfg_file = "cyl_ion_gy_march.cfg" - ion_gy_march.test_iter = 99 - ion_gy_march.test_vals = [-11.662039, -4.203178, -4.868257, -5.462497, -5.232052, -4.960881, -6.951391, 4.541901, 4.552855, -0.014861, 0.000001, 90357.000000] - ion_gy_march.tol = 0.01 - test_list.append(ion_gy_march) - ########################## ### Compressible Euler ### ########################## diff --git a/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp b/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp index 0fd184ff6a3a..93da083f6357 100644 --- a/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp +++ b/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp @@ -120,7 +120,7 @@ struct NEMOViscousAssemblyCase { } // namespace -TEST_CASE("NEMO viscous edge assembly is the derivative of the subtracted and added residual", +TEST_CASE("NEMO viscous solver assembles edge fluxes and Jacobian blocks with the correct signs and placement", "[NEMO][viscous][Jacobian]") { NEMOViscousAssemblyCase testCase; auto* config = testCase.config.get(); @@ -154,12 +154,11 @@ TEST_CASE("NEMO viscous edge assembly is the derivative of the subtracted and ad CSolver& base = *solver; base.Viscous_Residual(geometry, solver_container, numerics_container, config, MESH_0, 0); - /*--- Replay every edge with an independent numerics object. The residual - * contribution F of an edge is subtracted at i and added at j, so the - * derivative of the assembled residual is -dF/dU in the i rows and +dF/dU - * in the j rows, with dF/dU_i and dF/dU_j the two blocks returned by the - * numerics. Accumulate that reference and compare it with what the solver - * assembled. ---*/ + /*--- Replay every edge with a separate instance of the same numerics. + * This checks solver integration: each flux is subtracted at i and added + * at j, and the supplied i/j Jacobian blocks must have matching signs and + * positions. The reference reuses production numerics; it is not an + * independent derivative check of the full viscous Jacobian. ---*/ CAvgGradCorrected_NEMO replay(nDim, nVar, nPrimVar, nPrimVarGrad, config); auto* nodes = dynamic_cast(solver->GetNodes()); REQUIRE(nodes != nullptr); @@ -197,12 +196,15 @@ TEST_CASE("NEMO viscous edge assembly is the derivative of the subtracted and ad REQUIRE(block_ji != nullptr); for (unsigned short iVar = 0; iVar < nVar; ++iVar) { + REQUIRE(std::isfinite(SU2_TYPE::GetValue(edge.residual[iVar]))); residual_ref[iPoint * nVar + iVar] -= edge.residual[iVar]; residual_ref[jPoint * nVar + iVar] += edge.residual[iVar]; for (unsigned short jVar = 0; jVar < nVar; ++jVar) { const su2double dFdUi = edge.jacobian_i[iVar][jVar]; const su2double dFdUj = edge.jacobian_j[iVar][jVar]; + REQUIRE(std::isfinite(SU2_TYPE::GetValue(dFdUi))); + REQUIRE(std::isfinite(SU2_TYPE::GetValue(dFdUj))); jacobian_scale = std::max(jacobian_scale, std::max(std::fabs(dFdUi), std::fabs(dFdUj))); /*--- Diagonal blocks collect every edge of a point. ---*/ @@ -212,6 +214,8 @@ TEST_CASE("NEMO viscous edge assembly is the derivative of the subtracted and ad /*--- Off-diagonal blocks belong to this edge alone. ---*/ const su2double assembled_ij = block_ij[iVar * nVar + jVar]; const su2double assembled_ji = block_ji[iVar * nVar + jVar]; + REQUIRE(std::isfinite(SU2_TYPE::GetValue(assembled_ij))); + REQUIRE(std::isfinite(SU2_TYPE::GetValue(assembled_ji))); offdiag_error = std::max(offdiag_error, std::fabs(assembled_ij - (-dFdUj))); offdiag_error = std::max(offdiag_error, std::fabs(assembled_ji - (+dFdUi))); } @@ -224,11 +228,15 @@ TEST_CASE("NEMO viscous edge assembly is the derivative of the subtracted and ad const auto* assembled_residual = solver->LinSysRes.GetBlock(iPoint); REQUIRE(block_ii != nullptr); for (unsigned short iVar = 0; iVar < nVar; ++iVar) { + REQUIRE(std::isfinite(SU2_TYPE::GetValue(residual_ref[iPoint * nVar + iVar]))); + REQUIRE(std::isfinite(SU2_TYPE::GetValue(assembled_residual[iVar]))); residual_scale = std::max(residual_scale, std::fabs(residual_ref[iPoint * nVar + iVar])); residual_error = std::max(residual_error, std::fabs(assembled_residual[iVar] - residual_ref[iPoint * nVar + iVar])); for (unsigned short jVar = 0; jVar < nVar; ++jVar) { const su2double assembled = block_ii[iVar * nVar + jVar]; + REQUIRE(std::isfinite(SU2_TYPE::GetValue(assembled))); + REQUIRE(std::isfinite(SU2_TYPE::GetValue(diagonal_ref[(iPoint * nVar + iVar) * nVar + jVar]))); diagonal_error = std::max(diagonal_error, std::fabs(assembled - diagonal_ref[(iPoint * nVar + iVar) * nVar + jVar])); } diff --git a/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp index d9a32e937539..bea3da327e01 100644 --- a/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp +++ b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp @@ -27,6 +27,7 @@ #include "catch.hpp" #include +#include #include #include #include "../../../SU2_CFD/include/numerics/CNumerics.hpp" @@ -173,20 +174,6 @@ struct NEMOViscousFixture { numerics.SetdTvedU(dTve_i.data(), dTve_j.data()); } - void set_edge_linear_velocity_gradient() { - gradient_i = su2double(0.0); - gradient_j = su2double(0.0); - constexpr su2double distance_squared = 25.0; - for (unsigned short component = 0; component < nDim; ++component) { - const su2double jump = primitive_j[VEL_INDEX + component] - primitive_i[VEL_INDEX + component]; - for (unsigned short dimension = 0; dimension < nDim; ++dimension) { - const su2double edge_component = coord_j[dimension] - coord_i[dimension]; - gradient_i(VEL_INDEX + component, dimension) = jump * edge_component / distance_squared; - gradient_j(VEL_INDEX + component, dimension) = jump * edge_component / distance_squared; - } - } - } - template su2double directional_flux(Numerics& numerics) { set_common(numerics); @@ -213,10 +200,46 @@ TEST_CASE("NEMO corrected viscous residual returns distinct i and j Jacobians", CAvgGradCorrected_NEMO numerics(fixture.nDim, fixture.nVar, fixture.nPrimVar, fixture.nPrimVarGrad, fixture.config); fixture.set_common(numerics); const auto base = numerics.ComputeResidual(fixture.config); + /*--- Copy these values: later residual evaluations reuse the numerics' + * Jacobian storage. ---*/ const su2double analytic_i = fixture.directional_jacobian(base.jacobian_i, fixture.normal); const su2double analytic_j = fixture.directional_jacobian(base.jacobian_j, fixture.normal); + /*--- Independent Newtonian-stress result for this longitudinal mode: + * tau_nn = (4/3) mu du_n/dn. The unit face normal is parallel to the + * length-5 edge, with mu = 2 and rho = 1, hence the momentum derivatives + * are -(4/3) mu/(rho d) = -8/15 at i and +8/15 at j. ---*/ REQUIRE(analytic_i == Approx(-8.0 / 15.0).epsilon(1.0e-12)); REQUIRE(analytic_j == Approx(8.0 / 15.0).epsilon(1.0e-12)); REQUIRE(analytic_i != Approx(analytic_j)); + + /*--- Differentiate only the returned momentum flux, without using its + * Jacobians to construct the reference. Perturb each endpoint's momentum + * along the normal; density, thermodynamic inputs, transport and supplied + * zero gradients stay fixed. The corrected numerics reconstructs the + * edge gradient itself. This tests the controlled momentum mode, not the + * full thermochemical derivative of the approximate viscous Jacobian. ---*/ + for (unsigned short endpoint = 0; endpoint < 2; ++endpoint) { + auto& primitive = endpoint == 0 ? fixture.primitive_i : fixture.primitive_j; + const auto original = primitive; + for (const auto step : {1.0e-4, 1.0e-6}) { + INFO("endpoint = " << endpoint << ", momentum step = " << step); + for (unsigned short component = 0; component < fixture.nDim; ++component) + primitive[fixture.VEL_INDEX + component] = + original[fixture.VEL_INDEX + component] + step * fixture.normal[component] / original[fixture.RHO_INDEX]; + const su2double flux_plus = fixture.directional_flux(numerics); + + for (unsigned short component = 0; component < fixture.nDim; ++component) + primitive[fixture.VEL_INDEX + component] = + original[fixture.VEL_INDEX + component] - step * fixture.normal[component] / original[fixture.RHO_INDEX]; + const su2double flux_minus = fixture.directional_flux(numerics); + primitive = original; + + REQUIRE(std::isfinite(SU2_TYPE::GetValue(flux_plus))); + REQUIRE(std::isfinite(SU2_TYPE::GetValue(flux_minus))); + const su2double finite_difference = (flux_plus - flux_minus) / (2.0 * step); + const su2double analytic = endpoint == 0 ? analytic_i : analytic_j; + CHECK(finite_difference == Approx(analytic).epsilon(1.0e-8).margin(1.0e-10)); + } + } } From 50e3343c05105daf49ca9e69a504941489820979 Mon Sep 17 00:00:00 2001 From: BlueChips Date: Mon, 7 Sep 2026 16:29:42 -0400 Subject: [PATCH 5/5] Refresh cone references and load Mutation++ in unit CI --- .github/workflows/regression.yml | 1 + TestCases/parallel_regression.py | 2 +- TestCases/serial_regression.py | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index c4e3b652b069..76eb60ebc8f8 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -462,6 +462,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 env: OMPI_MCA_osc: pt2pt + LD_LIBRARY_PATH: /github/workspace/install/lib with: entrypoint: install/bin/${{matrix.testdriver}} - name: Post Cleanup diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 89f2fdaa9792..46b0d8228763 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -192,7 +192,7 @@ def main(): visc_cone.cfg_dir = "nonequilibrium/visc_wedge" visc_cone.cfg_file = "axi_visccone.cfg" visc_cone.test_iter = 10 - visc_cone.test_vals = [-5.298545, -5.823026, -20.404788, -20.318330, -20.378801, 1.067473, -3.250998, -0.015489, 0.095136, 24939.000000] + visc_cone.test_vals = [-5.298530, -5.823015, -20.404788, -20.318330, -20.378801, 1.067489, -3.250987, -0.015488, 0.095136, 24939.000000] visc_cone.test_vals_aarch64 = [-5.298530, -5.823015, -20.404788, -20.318329, -20.378801, 1.067489, -3.250987, -0.015488, 0.095136, 24939.000000] test_list.append(visc_cone) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index bc67350f59cd..7a646c7efb9f 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -75,7 +75,7 @@ def main(): visc_cone.cfg_dir = "nonequilibrium/visc_wedge" visc_cone.cfg_file = "axi_visccone.cfg" visc_cone.test_iter = 10 - visc_cone.test_vals = [-5.298535, -5.823020, -20.330692, -20.368713, -20.310533, 1.067483, -3.250989, -0.015487, 0.095136, 24939] + visc_cone.test_vals = [-5.298550, -5.823031, -20.330692, -20.368713, -20.310533, 1.067466, -3.251001, -0.015489, 0.095136, 24939.000000] visc_cone.test_vals_aarch64 = [-5.298550, -5.823031, -20.330692, -20.368713, -20.310533, 1.067466, -3.251001, -0.015489, 0.095136, 24939.000000] test_list.append(visc_cone)