Consistent NEMO viscous Jacobian: distinct left block and subtractive edge assembly

.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

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);
 
 }

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

TestCases/nonequilibrium/visc_cylinder/cyl_ion_gy.cfg

@@ -57,15 +57,15 @@ TIME_DISCRE_FLOW= EULER_IMPLICIT
 
 % --------------------------- CONVERGENCE PARAMETERS --------------------------%
 %
-CONV_RESIDUAL_MINVAL= -10
+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 
+RESTART_FILENAME= restart_flow_gy_out
 TABULAR_FORMAT= TECPLOT
 CONV_FILENAME= convergence
 VOLUME_FILENAME= soln_volume

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)
 
     ##########################

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)
 
     #########################

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 <http://www.gnu.org/licenses/>.
+ */
+
+#include "catch.hpp"
+#include <algorithm>
+#include <cmath>
+#include <memory>
+#include <sstream>
+#include <vector>
+#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<CConfig> config;
+  std::unique_ptr<CGeometry> geometry;
+  CNEMONSSolver* solver{nullptr};
+
+  NEMOViscousAssemblyCase() {
+    auto origBuf = cout.rdbuf();
+    cout.rdbuf(nullptr);
+
+    stringstream ss(configOptions);
+    config = std::unique_ptr<CConfig>(new CConfig(ss, SU2_COMPONENT::SU2_CFD, false));
+
+    {
+      auto aux_geometry = std::unique_ptr<CGeometry>(new CPhysicalGeometry(config.get(), 0, 1));
+      geometry = std::unique_ptr<CGeometry>(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<CNEMONSVariable*>(solver->GetNodes());
+  REQUIRE(nodes != nullptr);
+
+  std::vector<su2double> residual_ref(nPoint * nVar, 0.0);
+  std::vector<su2double> 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);
+}