Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 13 additions & 0 deletions SU2_CFD/include/solvers/CNEMOEulerSolver.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,19 @@ class CNEMOEulerSolver : public CFVMFlowSolverBase<CNEMOEulerVariable, ENUM_REGI
*/
void ComputeUnderRelaxationFactor(const CConfig *config) final;

/*!
* \brief Limit the species and total-energy updates at one point.
* \param[in] nSpecies - Number of species.
* \param[in] nVar - Number of conserved variables.
* \param[in] solution - Conserved variables at the point.
* \param[in] update - Proposed change in the conserved variables.
* \param[in] allowableRatio - Maximum allowed relative update.
* \return Under-relaxation factor for the point.
*/
static su2double ComputeUnderRelaxationFactor(unsigned short nSpecies, unsigned short nVar,
const su2double* solution, const su2double* update,
su2double allowableRatio);

/*!
* \brief Set the fluid solver nondimensionalization.
* \param[in] geometry - Geometrical definition of the problem.
Expand Down
79 changes: 37 additions & 42 deletions SU2_CFD/src/solvers/CNEMOEulerSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -936,56 +936,51 @@ void CNEMOEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**,
PrepareImplicitIteration_impl(precond, geometry, config);
}

void CNEMOEulerSolver::ComputeUnderRelaxationFactor(const CConfig *config) {
SU2_ZONE_SCOPED

/* Loop over the solution update given by relaxing the linear
system for this nonlinear iteration. */

const su2double allowableRatio = config->GetMaxUpdateFractionFlow();

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
Expand Down
28 changes: 14 additions & 14 deletions TestCases/parallel_regression.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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++
Expand All @@ -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
Expand All @@ -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)

##########################
Expand Down
8 changes: 4 additions & 4 deletions TestCases/serial_regression.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)

#########################
Expand Down
59 changes: 59 additions & 0 deletions UnitTests/SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, bootstraping the entire solver to test just this is too heavy handed, you want to have a unit test, make the logic in the loop a small re-usable function and feed it some manufactured data.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Got it, I moved the limiter logic into a small function and tested it directly with manufactured data.

Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.
*/

#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));
}
}
1 change: 1 addition & 0 deletions UnitTests/meson.build
Original file line number Diff line number Diff line change
Expand Up @@ -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'])
Expand Down
Loading