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/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..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 --------------------------%
@@ -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/parallel_regression.py b/TestCases/parallel_regression.py
index ff9b8c667061..46b0d8228763 100755
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -192,8 +192,8 @@ 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_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 = [-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)
# Viscous single wedge with Mutation++
@@ -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
@@ -220,12 +220,13 @@ 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 = [-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_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)
##########################
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index aa22e30d8ae9..7a646c7efb9f 100755
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -75,8 +75,8 @@ 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_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 = [-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)
#########################
diff --git a/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp b/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp
new file mode 100644
index 000000000000..93da083f6357
--- /dev/null
+++ b/UnitTests/SU2_CFD/nemo_viscous_assembly.cpp
@@ -0,0 +1,254 @@
+/*!
+ * \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 solver assembles edge fluxes and Jacobian blocks with the correct signs and placement",
+ "[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 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);
+
+ 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) {
+ 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. ---*/
+ 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];
+ 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)));
+ }
+ }
+ }
+
+ 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) {
+ 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]));
+ }
+ }
+ }
+
+ /*--- 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..bea3da327e01 100644
--- a/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp
+++ b/UnitTests/SU2_CFD/numerics/CNumerics_tests.cpp
@@ -26,8 +26,12 @@
*/
#include "catch.hpp"
+#include
+#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 +93,153 @@ 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());
+ }
+
+ 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);
+ /*--- 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));
+ }
+ }
+}
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'])