From 53375f99568833dd3edc582f32215ada1c7901b7 Mon Sep 17 00:00:00 2001 From: Jared Bellon Date: Thu, 3 Sep 2026 13:15:35 -0400 Subject: [PATCH 1/3] Pass the edge length, not its square, to the coarse NEMO viscous Jacobian CAvgGrad_NEMO::ComputeResidual assigned the squared norm of the edge vector to dist_ij and passed it to GetViscousProjJacs, which expects the distance. Take the square root, as CAvgGradCorrected_NEMO already does. Only the analytic viscous Jacobian changes; the residual is untouched. CDriver.cpp uses CAvgGrad_NEMO for the interior viscous term on every coarse multigrid level and for the boundary viscous term on every level including MESH_0, where it is consumed only by BC_Far_Field, BC_Supersonic_Inlet and the inherited BC_Fluid_Interface (the NEMO wall routines build their own fluxes). The defect therefore reaches every NEMO Navier-Stokes run that uses multigrid, a far-field or supersonic-inlet marker, or a fluid interface. Add a commented finite-difference unit test on a 2-D AIR-5 edge (NEMOViscousFixture) that fails on develop (8/75) and passes with this change (8/15), and rebaseline the regression vectors that move as a consequence: the machine-zero species columns of visc_cone (serial and parallel) and super_cat (parallel), and the physical residual columns rms[Rho_0..5] and rms[RhoE] of ion_gy (parallel), which shift by up to 1.05e-2 because that case has no identically-zero species. Co-Authored-By: Claude Fable 5.1 --- SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp | 2 +- TestCases/parallel_regression.py | 6 +- TestCases/serial_regression.py | 2 +- .../SU2_CFD/numerics/CNumerics_tests.cpp | 179 ++++++++++++++++++ 4 files changed, 184 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp b/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp index 2553fa2915d..4e57cdaa875 100644 --- a/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp +++ b/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp @@ -169,7 +169,7 @@ CNumerics::ResidualType<> CAvgGrad_NEMO::ComputeResidual(const CConfig *config) su2double dist_ij_2[MAXNDIM] = {0.0}; GeometryToolbox::Distance(nDim, Coord_j, Coord_i, dist_ij_2); - dist_ij = GeometryToolbox::SquaredNorm(nDim, dist_ij_2); + dist_ij = sqrt(GeometryToolbox::SquaredNorm(nDim, dist_ij_2)); for (auto iVar = 0ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < nVar; jVar++) { diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ff9b8c66706..fab53956115 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.215238, -5.739375, -20.581529, -20.580261, -20.508727, 1.262697, -3.205457, -0.015697, 0.093206, 32637.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.232595, -5.757889, -20.726945, -20.749655, -20.567151, 1.246866, -3.205258, -0.028372, 0.250647, 32440.000000] test_list.append(super_cat) # Viscous single wedge - partially catalytic walls @@ -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 = 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 = [-11.629871, -4.156741, -4.692113, -4.949447, -5.146196, -4.991664, -6.893329, 5.990109, 5.990004, -0.014849, 0.000000, 90090.000000] test_list.append(ion_gy) ########################## diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index aa22e30d8ae..c9e16d6ffdc 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.215233, -5.739371, -20.551239, -20.590906, -20.437427, 1.262779, -3.205476, -0.015697, 0.093207, 32641] 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/numerics/CNumerics_tests.cpp b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp index a4cee41e8d2..1324ed10dcf 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,179 @@ 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 { + +/*! + * \brief Fixture for the NEMO viscous numerics on a single 2-D AIR-5 edge. + * + * The edge runs from node i at (0, 0) to node j at (3, 4), so its length is 5 and its squared length + * is 25, and the unit face normal (0.6, 0.8) is the edge direction. Only the laminar viscosity (2) + * is non-zero among the transport coefficients, so the normal-projected momentum flux is a pure + * viscous stress. The analytic NEMO viscous Jacobian (CNEMONumerics::GetViscousProjJacs) builds its + * momentum block as mu / d * (theta I + n n^T / 3) with theta = |n|^2, so contracting it twice with + * the unit normal gives mu * (4 / 3) / d = 8 / 15 in magnitude on this edge: -8/15 on the i side + * and +8/15 on the j side. Passing the squared length 25 in place of the distance 5 gives 8/75 + * instead, which is what CAvgGrad_NEMO returned before the square root was added. + */ +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()); + } + + /*! + * \brief Reconstruct the gradient of a velocity field that varies linearly along the edge, so the + * nodal jump and the averaged gradient that the numerics consume describe the same profile. + */ + 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 + +/*! + * Compare the analytic jacobian_i and jacobian_j of CAvgGrad_NEMO, contracted with the unit normal, + * against a central finite difference of the normal-projected momentum flux with respect to the + * normal velocity at each node. The gradient is rebuilt after every perturbation so the flux sees + * the same linear profile that the Jacobian assumes. On the fixture edge (length 5, mu = 2) the + * exact values are -8/15 and +8/15; the un-rooted squared length gives 8/75 and fails both checks. + */ +TEST_CASE("NEMO coarse viscous Jacobian uses Euclidean edge length", "[NEMO][viscous][Jacobian]") { + NEMOViscousFixture fixture; + CAvgGrad_NEMO numerics(fixture.nDim, fixture.nVar, fixture.nPrimVar, fixture.nPrimVarGrad, fixture.config); + fixture.set_edge_linear_velocity_gradient(); + 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); + + constexpr su2double step = 1.0e-6; + for (unsigned short component = 0; component < fixture.nDim; ++component) + fixture.primitive_i[fixture.VEL_INDEX + component] += step * fixture.normal[component]; + fixture.set_edge_linear_velocity_gradient(); + const su2double plus_i = fixture.directional_flux(numerics); + for (unsigned short component = 0; component < fixture.nDim; ++component) + fixture.primitive_i[fixture.VEL_INDEX + component] -= 2.0 * step * fixture.normal[component]; + fixture.set_edge_linear_velocity_gradient(); + const su2double minus_i = fixture.directional_flux(numerics); + for (unsigned short component = 0; component < fixture.nDim; ++component) + fixture.primitive_i[fixture.VEL_INDEX + component] += step * fixture.normal[component]; + + for (unsigned short component = 0; component < fixture.nDim; ++component) + fixture.primitive_j[fixture.VEL_INDEX + component] += step * fixture.normal[component]; + fixture.set_edge_linear_velocity_gradient(); + const su2double plus_j = fixture.directional_flux(numerics); + for (unsigned short component = 0; component < fixture.nDim; ++component) + fixture.primitive_j[fixture.VEL_INDEX + component] -= 2.0 * step * fixture.normal[component]; + fixture.set_edge_linear_velocity_gradient(); + const su2double minus_j = fixture.directional_flux(numerics); + + const su2double fd_i = (plus_i - minus_i) / (2.0 * step); + const su2double fd_j = (plus_j - minus_j) / (2.0 * step); + REQUIRE(analytic_i == Approx(fd_i).epsilon(1.0e-8)); + REQUIRE(analytic_j == Approx(fd_j).epsilon(1.0e-8)); + REQUIRE(analytic_i == Approx(-8.0 / 15.0).epsilon(1.0e-12)); + REQUIRE(analytic_j == Approx(8.0 / 15.0).epsilon(1.0e-12)); +} From ab76f7f1d45b4da4922c59aaaabf628fe4f3693f Mon Sep 17 00:00:00 2001 From: Jared Bellon Date: Sat, 5 Sep 2026 18:44:27 -0400 Subject: [PATCH 2/3] Address review: use Norm and remove duplicated fixture Use GeometryToolbox::Norm directly and remove the 179-line NEMO fixture/test at maintainer request. The affected public NEMO regression cases retain path coverage. Co-Authored-By: Claude Fable 5.1 --- SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp | 2 +- .../SU2_CFD/numerics/CNumerics_tests.cpp | 179 ------------------ 2 files changed, 1 insertion(+), 180 deletions(-) diff --git a/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp b/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp index 4e57cdaa875..5a23ff0b6ce 100644 --- a/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp +++ b/SU2_CFD/src/numerics/NEMO/NEMO_diffusion.cpp @@ -169,7 +169,7 @@ CNumerics::ResidualType<> CAvgGrad_NEMO::ComputeResidual(const CConfig *config) su2double dist_ij_2[MAXNDIM] = {0.0}; GeometryToolbox::Distance(nDim, Coord_j, Coord_i, dist_ij_2); - dist_ij = sqrt(GeometryToolbox::SquaredNorm(nDim, dist_ij_2)); + dist_ij = GeometryToolbox::Norm(nDim, dist_ij_2); for (auto iVar = 0ul; iVar < nVar; iVar++) { for (auto jVar = 0ul; jVar < nVar; jVar++) { diff --git a/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp index 1324ed10dcf..a4cee41e8d2 100644 --- a/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp +++ b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp @@ -26,11 +26,8 @@ */ #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; @@ -92,179 +89,3 @@ 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 { - -/*! - * \brief Fixture for the NEMO viscous numerics on a single 2-D AIR-5 edge. - * - * The edge runs from node i at (0, 0) to node j at (3, 4), so its length is 5 and its squared length - * is 25, and the unit face normal (0.6, 0.8) is the edge direction. Only the laminar viscosity (2) - * is non-zero among the transport coefficients, so the normal-projected momentum flux is a pure - * viscous stress. The analytic NEMO viscous Jacobian (CNEMONumerics::GetViscousProjJacs) builds its - * momentum block as mu / d * (theta I + n n^T / 3) with theta = |n|^2, so contracting it twice with - * the unit normal gives mu * (4 / 3) / d = 8 / 15 in magnitude on this edge: -8/15 on the i side - * and +8/15 on the j side. Passing the squared length 25 in place of the distance 5 gives 8/75 - * instead, which is what CAvgGrad_NEMO returned before the square root was added. - */ -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()); - } - - /*! - * \brief Reconstruct the gradient of a velocity field that varies linearly along the edge, so the - * nodal jump and the averaged gradient that the numerics consume describe the same profile. - */ - 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 - -/*! - * Compare the analytic jacobian_i and jacobian_j of CAvgGrad_NEMO, contracted with the unit normal, - * against a central finite difference of the normal-projected momentum flux with respect to the - * normal velocity at each node. The gradient is rebuilt after every perturbation so the flux sees - * the same linear profile that the Jacobian assumes. On the fixture edge (length 5, mu = 2) the - * exact values are -8/15 and +8/15; the un-rooted squared length gives 8/75 and fails both checks. - */ -TEST_CASE("NEMO coarse viscous Jacobian uses Euclidean edge length", "[NEMO][viscous][Jacobian]") { - NEMOViscousFixture fixture; - CAvgGrad_NEMO numerics(fixture.nDim, fixture.nVar, fixture.nPrimVar, fixture.nPrimVarGrad, fixture.config); - fixture.set_edge_linear_velocity_gradient(); - 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); - - constexpr su2double step = 1.0e-6; - for (unsigned short component = 0; component < fixture.nDim; ++component) - fixture.primitive_i[fixture.VEL_INDEX + component] += step * fixture.normal[component]; - fixture.set_edge_linear_velocity_gradient(); - const su2double plus_i = fixture.directional_flux(numerics); - for (unsigned short component = 0; component < fixture.nDim; ++component) - fixture.primitive_i[fixture.VEL_INDEX + component] -= 2.0 * step * fixture.normal[component]; - fixture.set_edge_linear_velocity_gradient(); - const su2double minus_i = fixture.directional_flux(numerics); - for (unsigned short component = 0; component < fixture.nDim; ++component) - fixture.primitive_i[fixture.VEL_INDEX + component] += step * fixture.normal[component]; - - for (unsigned short component = 0; component < fixture.nDim; ++component) - fixture.primitive_j[fixture.VEL_INDEX + component] += step * fixture.normal[component]; - fixture.set_edge_linear_velocity_gradient(); - const su2double plus_j = fixture.directional_flux(numerics); - for (unsigned short component = 0; component < fixture.nDim; ++component) - fixture.primitive_j[fixture.VEL_INDEX + component] -= 2.0 * step * fixture.normal[component]; - fixture.set_edge_linear_velocity_gradient(); - const su2double minus_j = fixture.directional_flux(numerics); - - const su2double fd_i = (plus_i - minus_i) / (2.0 * step); - const su2double fd_j = (plus_j - minus_j) / (2.0 * step); - REQUIRE(analytic_i == Approx(fd_i).epsilon(1.0e-8)); - REQUIRE(analytic_j == Approx(fd_j).epsilon(1.0e-8)); - REQUIRE(analytic_i == Approx(-8.0 / 15.0).epsilon(1.0e-12)); - REQUIRE(analytic_j == Approx(8.0 / 15.0).epsilon(1.0e-12)); -} From b8504d8de172a20a0493391d3f3a3fd572b0cfbd Mon Sep 17 00:00:00 2001 From: BlueChips Date: Sat, 5 Sep 2026 20:09:48 -0400 Subject: [PATCH 3/3] 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 ab76f7f1d45b4da4922c59aaaabf628fe4f3693f. Artifact digest: sha256:279e65a52098a2829127441d352dd4d88d0979ae7fdb8870a9c779d076a537e6. 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 fab53956115..69a3f2fb5ae 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.215238, -5.739375, -20.581529, -20.580261, -20.508727, 1.262697, -3.205457, -0.015697, 0.093206, 32637.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.215222, -5.739359, -20.581535, -20.580239, -20.508657, 1.262709, -3.205449, -0.015696, 0.093205, 32637.000000] test_list.append(visc_cone) # Viscous single wedge with Mutation++ diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index c9e16d6ffdc..82df6232a97 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.215233, -5.739371, -20.551239, -20.590906, -20.437427, 1.262779, -3.205476, -0.015697, 0.093207, 32641] - 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.215220, -5.739356, -20.551244, -20.590883, -20.437378, 1.262790, -3.205470, -0.015696, 0.093206, 32641.000000] test_list.append(visc_cone) #########################