Make the NEMO total-energy under-relaxation check reachable

SU2_CFD/src/solvers/CNEMOEulerSolver.cpp

@@ -968,12 +968,14 @@ void CNEMOEulerSolver::ComputeUnderRelaxationFactor(const CConfig *config) {
           }
         }
 
-        /*--- Energy ---*/
-        if (iVar == (nVar-2)){
-          su2double ratio = fabs(LinSysSol[index]) / (fabs(nodes->GetSolution(iPoint, iVar)) + EPS);
-          if (ratio > allowableRatio) {
-            localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
-          }
+      }
+
+      /*--- Total energy. This check must be independent of the species block:
+       * nVar-2 is nSpecies+nDim and can never satisfy iVar < nSpecies. ---*/
+      if (iVar == (nVar-2)){
+        su2double ratio = fabs(LinSysSol[index]) / (fabs(nodes->GetSolution(iPoint, iVar)) + EPS);
+        if (ratio > allowableRatio) {
+          localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
         }
       }
     }

TestCases/parallel_regression.py

@@ -129,8 +129,8 @@ def main():
     invwedge_a.cfg_dir = "nonequilibrium/invwedge"
     invwedge_a.cfg_file = "invwedge_ausm.cfg"
     invwedge_a.test_iter = 10
-    invwedge_a.test_vals = [-1.069665, -1.594428, -18.299923, -18.627315, -18.573325, 2.245732, 1.874096, 5.290295, 0.847739]
-    invwedge_a.test_vals_aarch64 = [-1.069675, -1.594438, -18.299736, -18.627126, -18.573137, 2.245721, 1.874105, 5.290285, 0.847729]
+    invwedge_a.test_vals = [-1.081764, -1.606527, -18.299923, -18.627313, -18.573325, 2.234595, 1.854912, 5.278328, 0.837083]
+    invwedge_a.test_vals_aarch64 = [-1.081764, -1.606527, -18.299923, -18.627313, -18.573325, 2.234595, 1.854912, 5.278328, 0.837083]
     test_list.append(invwedge_a)
 
     # Inviscid single wedge, ausm+-up2, implicit
@@ -147,26 +147,26 @@ def main():
     invwedge_msw.cfg_dir = "nonequilibrium/invwedge"
     invwedge_msw.cfg_file = "invwedge_msw.cfg"
     invwedge_msw.test_iter = 10
-    invwedge_msw.test_vals = [-1.212335, -1.737098, -18.301825, -18.629206, -18.575226, 2.106171, 1.651949, 5.143958, 0.704444]
-    invwedge_msw.test_vals_aarch64 = [-1.212335, -1.737098, -18.299279, -18.626656, -18.572683, 2.106171, 1.651949, 5.143958, 0.704444]
+    invwedge_msw.test_vals = [-1.206414, -1.731177, -18.301107, -18.628488, -18.574512, 2.111923, 1.660408, 5.150501, 0.710186]
+    invwedge_msw.test_vals_aarch64 = [-1.206414, -1.731177, -18.299771, -18.627181, -18.573171, 2.111923, 1.660408, 5.150501, 0.710186]
     test_list.append(invwedge_msw)
 
     # Inviscid single wedge, roe, implicit
     invwedge_roe = TestCase('invwedge_roe')
     invwedge_roe.cfg_dir = "nonequilibrium/invwedge"
     invwedge_roe.cfg_file = "invwedge_roe.cfg"
     invwedge_roe.test_iter = 10
-    invwedge_roe.test_vals = [-1.023283, -1.548046, -17.814403, -18.143369, -18.087522, 2.295025, 1.884804, 5.338440, 0.926068]
-    invwedge_roe.test_vals_aarch64 = [-1.052398, -1.577160, -17.794015, -18.122997, -18.067131, 2.266042, 1.849686, 5.304700, 0.899584]
+    invwedge_roe.test_vals = [-1.158548, -1.683311, -17.731949, -18.061692, -18.005058, 2.157602, 1.789805, 5.196306, 0.767203]
+    invwedge_roe.test_vals_aarch64 = [-1.156325, -1.681088, -17.726749, -18.056513, -17.999860, 2.159488, 1.793498, 5.198638, 0.769831]
     test_list.append(invwedge_roe)
 
     # Inviscid single wedge, lax, implicit
     invwedge_lax = TestCase('invwedge_lax')
     invwedge_lax.cfg_dir = "nonequilibrium/invwedge"
     invwedge_lax.cfg_file = "invwedge_lax.cfg"
     invwedge_lax.test_iter = 10
-    invwedge_lax.test_vals = [-0.877280, -1.402043, -32.000000, -32.000000, -24.952631, 2.451869, 1.857084, 5.486158, 1.051580]
-    invwedge_lax.test_vals_aarch64 = [-0.877280, -1.402043, -32.000000, -32.000000, -24.952631, 2.451869, 1.857084, 5.486158, 1.051580]
+    invwedge_lax.test_vals = [-0.882120, -1.406883, -32.000000, -32.000000, -24.953606, 2.447071, 1.853423, 5.480164, 1.047015]
+    invwedge_lax.test_vals_aarch64 = [-0.882120, -1.406883, -32.000000, -32.000000, -24.953606, 2.447071, 1.853423, 5.480164, 1.047015]
     test_list.append(invwedge_lax)
 
     # Inviscid single wedge, implicit, AUSM+M scheme
@@ -183,17 +183,17 @@ def main():
     invwedge_ss_inlet.cfg_dir = "nonequilibrium/invwedge"
     invwedge_ss_inlet.cfg_file = "invwedge_ss_inlet.cfg"
     invwedge_ss_inlet.test_iter = 10
-    invwedge_ss_inlet.test_vals = [-1.068634, -1.593397, -18.246265, -18.575529, -18.519338, 2.246925, 1.874200, 5.291234, 0.848731]
-    invwedge_ss_inlet.test_vals_aarch64 = [-1.068592, -1.593355, -18.250183, -18.579524, -18.523255, 2.246972, 1.874197, 5.291273, 0.848771]
+    invwedge_ss_inlet.test_vals = [-1.081061, -1.605824, -18.246258, -18.575522, -18.519330, 2.235398, 1.855022, 5.278994, 0.837755]
+    invwedge_ss_inlet.test_vals_aarch64 = [-1.081061, -1.605824, -18.246260, -18.575547, -18.519333, 2.235398, 1.855022, 5.278994, 0.837755]
     test_list.append(invwedge_ss_inlet)
 
     # Viscous single cone - axisymmetric
     visc_cone = TestCase('visc_cone')
     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.113779, -5.637767, -20.549789, -20.502278, -20.414074, 1.473388, -3.120680, -0.014621, 0.092458, 42636.000000]
+    visc_cone.test_vals_aarch64 = [-5.113765, -5.637751, -20.549769, -20.502286, -20.414073, 1.473399, -3.120671, -0.014620, 0.092456, 42636.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.139246, -5.664499, -20.627913, -20.635875, -20.535513, 1.463689, -3.125078, -0.030085, 0.252806, 42282.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.630651, -4.166100, -4.705260, -5.097030, -5.335247, -5.069078, -6.893926, 5.897432, 5.897354, -0.014849, 0.000000, 90090.000000]
     test_list.append(ion_gy)
 
     ##########################

TestCases/serial_regression.py

@@ -66,17 +66,17 @@ def main():
     invwedge.cfg_dir = "nonequilibrium/invwedge"
     invwedge.cfg_file = "invwedge_ausm.cfg"
     invwedge.test_iter = 10
-    invwedge.test_vals = [-1.073689, -1.598452, -18.299910, -18.627322, -18.573334, 2.241771, 1.868566, 5.286082, 0.843751]
-    invwedge.test_vals_aarch64 = [-1.073699, -1.598462, -18.299723, -18.627132, -18.573146, 2.241760, 1.868575, 5.286072, 0.843741]
+    invwedge.test_vals = [-1.085516, -1.610279, -18.299901, -18.627313, -18.573325, 2.230829, 1.850866, 5.274272, 0.833227]
+    invwedge.test_vals_aarch64 = [-1.085516, -1.610279, -18.299901, -18.627313, -18.573325, 2.230829, 1.850866, 5.274272, 0.833227]
     test_list.append(invwedge)
 
     # Viscous single cone - axisymmetric
     visc_cone = TestCase('visc_cone')
     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.113745, -5.637733, -20.558456, -20.511804, -20.402007, 1.473443, -3.120683, -0.014620, 0.092459, 42641]
+    visc_cone.test_vals_aarch64 = [-5.113728, -5.637715, -20.558435, -20.511813, -20.402006, 1.473454, -3.120672, -0.014619, 0.092457, 42641.000000]
     test_list.append(visc_cone)
 
     #########################

UnitTests/SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp

@@ -0,0 +1,116 @@
+/*!
+ * \file CNEMOEulerSolver_tests.cpp
+ * \brief Unit tests for the NEMO Euler solver.
+ * \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 "../../UnitQuadTestCase.hpp"
+#include "../../../SU2_CFD/include/solvers/CNEMOEulerSolver.hpp"
+
+/*!
+ * \brief Build a NEMO Euler solver on the unit box and check which conserved
+ * variables the implicit update limiter (MAX_UPDATE_FLOW) reacts to. The
+ * total-energy check used to be nested inside the species block and could
+ * never run, so the factor stayed at one for arbitrarily large energy updates.
+ */
+TEST_CASE("NEMO under-relaxation reacts to the total-energy update", "[NEMO][Solver]") {
+  UnitQuadTestCase testCase;
+  testCase.config_options =
+      "SOLVER= NEMO_EULER\n"
+      "FLUID_MODEL= SU2_NONEQ\n"
+      "GAS_MODEL= AIR-5\n"
+      "GAS_COMPOSITION= (0.77, 0.23, 0.0, 0.0, 0.0)\n"
+      "MESH_FORMAT= BOX\n"
+      "MESH_BOX_SIZE= 3,3,3\n"
+      "MESH_BOX_LENGTH= 1,1,1\n"
+      "MESH_BOX_OFFSET= 0,0,0\n"
+      "INIT_OPTION= TD_CONDITIONS\n"
+      "MACH_NUMBER= 5.0\n"
+      "FREESTREAM_PRESSURE= 101325.0\n"
+      "FREESTREAM_TEMPERATURE= 288.15\n"
+      "FREESTREAM_TEMPERATURE_VE= 288.15\n"
+      "MARKER_FAR= (x_minus, x_plus, y_minus, y_plus, z_minus, z_plus)\n"
+      "TIME_DISCRE_FLOW= EULER_IMPLICIT\n"
+      "MAX_UPDATE_FLOW= 0.2\n";
+  testCase.InitConfig();
+  testCase.InitGeometry();
+
+  CConfig* config = testCase.config.get();
+  CGeometry* geometry = testCase.geometry.get();
+
+  /*--- Construct the class under test directly instead of going through
+   * UnitQuadTestCase::InitSolver (CSolverFactory), so the test depends on
+   * CNEMOEulerSolver alone. The constructor reports to cout; silence it the
+   * way the harness helpers do. ---*/
+  std::cout.rdbuf(nullptr);
+  CNEMOEulerSolver solver(geometry, config, MESH_0);
+  std::cout.rdbuf(testCase.orig_buf);
+
+  const unsigned short nSpecies = config->GetnSpecies();
+  const unsigned short nDim = geometry->GetnDim();
+  const unsigned short nVar = solver.GetnVar();
+  REQUIRE(nVar == nSpecies + nDim + 2);
+
+  const unsigned long nPointDomain = geometry->GetnPointDomain();
+  REQUIRE(nPointDomain > 1);
+  const unsigned long iPoint = nPointDomain / 2;
+  const unsigned long jPoint = (iPoint + 1) % nPointDomain;
+
+  const su2double allowableRatio = config->GetMaxUpdateFractionFlow();
+  REQUIRE(allowableRatio == Approx(0.2));
+
+  /*--- Exceed the allowable update fraction by this factor for one variable at
+   * one point. If the limiter sees the variable, the factor drops to 1/excess. ---*/
+  const su2double excess = 4.0;
+  CVariable* nodes = solver.GetNodes();
+
+  auto factorWithUpdateOn = [&](unsigned short iVar, su2double reference) {
+    solver.LinSysSol.SetValZero();
+    solver.LinSysSol(iPoint, iVar) = excess * allowableRatio * reference;
+    solver.ComputeUnderRelaxationFactor(config);
+    return nodes->GetUnderRelaxation(iPoint);
+  };
+
+  const unsigned short rhoE = nVar - 2;
+  REQUIRE(rhoE >= nSpecies);
+
+  /*--- Control: the species check compares the summed species update with the
+   * mixture density, so the reference for the first species is the density. ---*/
+  su2double density = 0.0;
+  for (unsigned short iSpecies = 0; iSpecies < nSpecies; ++iSpecies)
+    density += fabs(nodes->GetSolution(iPoint, iSpecies));
+  CHECK(factorWithUpdateOn(0, density + EPS) == Approx(1.0 / excess).margin(1e-12));
+  CHECK(nodes->GetUnderRelaxation(jPoint) == 1.0);
+
+  /*--- The total-energy update must be limited by the same rule. ---*/
+  const su2double energy = fabs(nodes->GetSolution(iPoint, rhoE));
+  CHECK(factorWithUpdateOn(rhoE, energy + EPS) == Approx(1.0 / excess).margin(1e-12));
+  CHECK(nodes->GetUnderRelaxation(jPoint) == 1.0);
+
+  /*--- An energy update within the allowable fraction leaves the factor at one. ---*/
+  solver.LinSysSol.SetValZero();
+  solver.LinSysSol(iPoint, rhoE) = 0.5 * allowableRatio * energy;
+  solver.ComputeUnderRelaxationFactor(config);
+  CHECK(nodes->GetUnderRelaxation(iPoint) == 1.0);
+}

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/solvers/CNEMOEulerSolver_tests.cpp',
                        'SU2_CFD/windowing.cpp',
                        'Common/toolboxes/random_toolbox_tests.cpp',
                        'Common/linear_algebra/quantization_tests.cpp'])