diff --git a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp index e5d3c2d3c89..f1955c46a1a 100644 --- a/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CNEMOEulerSolver.hpp @@ -199,6 +199,19 @@ class CNEMOEulerSolver : public CFVMFlowSolverBaseGetMaxUpdateFractionFlow(); - - SU2_OMP_FOR_STAT(omp_chunk_size) - for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { - su2double localUnderRelaxation = 1.0; - - su2double num = 0.0; - su2double denom = 0.0; - - for (auto iVar = 0; iVar < nVar; iVar++) { - /* We impose a limit on the maximum percentage that the - density (sum of all species) and energy can change over a nonlinear iteration. */ - - const unsigned long index = iPoint * nVar + iVar; - if (iVar < config->GetnSpecies()) { - num += fabs(LinSysSol[index]); - denom += fabs(nodes->GetSolution(iPoint, iVar)); - - /*--- If final density/species, compute Under-relaxation ---*/ - if (iVar == (config ->GetnSpecies()-1)){ - su2double ratio = (num/(denom+EPS)); - if (ratio > allowableRatio) { - localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation); - } +su2double CNEMOEulerSolver::ComputeUnderRelaxationFactor(unsigned short nSpecies, unsigned short nVar, + const su2double* solution, const su2double* update, + su2double allowableRatio) { + su2double localUnderRelaxation = 1.0; + su2double num = 0.0; + su2double denom = 0.0; + + for (auto iVar = 0; iVar < nVar; iVar++) { + /*--- Limit the sum of the species updates relative to the mixture density. ---*/ + if (iVar < nSpecies) { + num += fabs(update[iVar]); + denom += fabs(solution[iVar]); + + if (iVar == nSpecies - 1) { + su2double ratio = num / (denom + EPS); + if (ratio > allowableRatio) { + localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation); } + } + } - /*--- 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 must be checked independently of the species block. ---*/ + if (iVar == nVar - 2) { + su2double ratio = fabs(update[iVar]) / (fabs(solution[iVar]) + EPS); + if (ratio > allowableRatio) { + localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation); } } + } - /* Threshold the relaxation factor in the event that there is - a very small value. This helps avoid catastrophic crashes due - to non-realizable states by canceling the update. */ + /*--- Cancel very small updates to avoid non-realizable states. ---*/ + if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0; + return localUnderRelaxation; +} - if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0; +void CNEMOEulerSolver::ComputeUnderRelaxationFactor(const CConfig *config) { + SU2_ZONE_SCOPED - /* Store the under-relaxation factor for this point. */ + const su2double allowableRatio = config->GetMaxUpdateFractionFlow(); + const unsigned short nSpecies = config->GetnSpecies(); + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + const su2double localUnderRelaxation = ComputeUnderRelaxationFactor( + nSpecies, nVar, nodes->GetSolution(iPoint), LinSysSol.GetBlock(iPoint), allowableRatio); nodes->SetUnderRelaxation(iPoint, localUnderRelaxation); } END_SU2_OMP_FOR diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index ff9b8c66706..3db05b02d75 100755 --- a/TestCases/parallel_regression.py +++ b/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,8 +147,8 @@ 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 @@ -156,8 +156,8 @@ def main(): 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 @@ -165,8 +165,8 @@ def main(): 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,8 +183,8 @@ 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 @@ -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.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) ########################## diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index aa22e30d8ae..ef98361c896 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -66,8 +66,8 @@ 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 @@ -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.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) ######################### diff --git a/UnitTests/SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp b/UnitTests/SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp new file mode 100644 index 00000000000..bbc4268f4c9 --- /dev/null +++ b/UnitTests/SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp @@ -0,0 +1,59 @@ +/*! + * \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 . + */ + +#include "catch.hpp" +#include "../../../SU2_CFD/include/solvers/CNEMOEulerSolver.hpp" + +TEST_CASE("NEMO under-relaxation limits species and total-energy updates", "[NEMO][Solver]") { + const unsigned short nSpecies = 2; + const unsigned short nVar = 6; + /*--- Two species, two momentum components, total energy, and vibrational energy. ---*/ + const su2double solution[nVar] = {0.75, 0.25, 1.0, 1.0, 10.0, 2.0}; + const su2double allowableRatio = 0.2; + + const struct { + const char* name; + su2double update[nVar]; + su2double expected; + } cases[] = { + {"zero update", {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, 1.0}, + {"small energy update", {0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, 1.0}, + {"positive energy excess", {0.0, 0.0, 0.0, 0.0, 8.0, 0.0}, 0.25}, + {"negative energy excess", {0.0, 0.0, 0.0, 0.0, -8.0, 0.0}, 0.25}, + {"opposing species updates", {0.4, -0.4, 0.0, 0.0, 0.0, 0.0}, 0.25}, + {"species sets the tighter limit", {0.8, -0.8, 0.0, 0.0, 8.0, 0.0}, 0.125}, + {"energy sets the tighter limit", {0.2, -0.2, 0.0, 0.0, 8.0, 0.0}, 0.25}, + {"momentum and vibrational energy are not limited", {0.0, 0.0, 100.0, -100.0, 0.0, 100.0}, 1.0}, + {"tiny factor cancels the update", {0.0, 0.0, 0.0, 0.0, 1e12, 0.0}, 0.0}, + }; + + for (const auto& test : cases) { + CAPTURE(test.name); + const su2double factor = + CNEMOEulerSolver::ComputeUnderRelaxationFactor(nSpecies, nVar, solution, test.update, allowableRatio); + CHECK(factor == Approx(test.expected).margin(1e-12)); + } +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index f8c22511b5f..ce909f46bd0 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/solvers/CNEMOEulerSolver_tests.cpp', 'SU2_CFD/windowing.cpp', 'Common/toolboxes/random_toolbox_tests.cpp', 'Common/linear_algebra/quantization_tests.cpp'])