Skip to content

Fix NEMO history MAX_/BGS_ residual groups: species-first indexing and finite maximum residuals #2882

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bellonarts wants to merge 2 commits into
base: develop
Choose a base branch
from
+265 −28
Open

Fix NEMO history MAX_/BGS_ residual groups: species-first indexing and finite maximum residuals #2882

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4751737
Fix NEMO history MAX_/BGS_ residual groups: species-first indexing an…
bellonarts Aug 22, 2026
20a6f50
Add NEMO residual output regression coverage
bellonarts Sep 5, 2026
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
70 changes: 42 additions & 28 deletions SU2_CFD/src/output/CNEMOCompOutput.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,30 +124,36 @@ void CNEMOCompOutput::SetHistoryOutputFields(CConfig *config){
/// END_GROUP

/// BEGIN_GROUP: MAX_RES, DESCRIPTION: The maximum residuals of the SOLUTION variables.
/// DESCRIPTION: Maximum residual of the density.
AddHistoryOutput("MAX_DENSITY", "max[Rho]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum square residual of the density.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the species densities.
for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
AddHistoryOutput("MAX_DENSITY_" + std::to_string(iSpecies), "max[Rho_" + std::to_string(iSpecies) + "]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the species density " + std::to_string(iSpecies) + ".", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the momentum x-component.
AddHistoryOutput("MAX_MOMENTUM-X", "max[RhoU]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum square residual of the momentum x-component.", HistoryFieldType::RESIDUAL);
AddHistoryOutput("MAX_MOMENTUM-X", "max[RhoU]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum x-component.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the momentum y-component.
AddHistoryOutput("MAX_MOMENTUM-Y", "max[RhoV]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum square residual of the momentum y-component.", HistoryFieldType::RESIDUAL);
AddHistoryOutput("MAX_MOMENTUM-Y", "max[RhoV]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the momentum y-component.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the momentum z-component.
if (nDim == 3) AddHistoryOutput("MAX_MOMENTUM-Z", "max[RhoW]", ScreenOutputFormat::FIXED,"MAX_RES", "Maximum residual of the z-component.", HistoryFieldType::RESIDUAL);
if (nDim == 3) AddHistoryOutput("MAX_MOMENTUM-Z", "max[RhoW]", ScreenOutputFormat::FIXED,"MAX_RES", "Maximum residual of the momentum z-component.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the energy.
AddHistoryOutput("MAX_ENERGY", "max[RhoE]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the energy.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the vibrational-electronic energy.
AddHistoryOutput("MAX_ENERGY_VE", "max[RhoEve]", ScreenOutputFormat::FIXED, "MAX_RES", "Maximum residual of the vibrational-electronic energy.", HistoryFieldType::RESIDUAL);
AddHistoryOutputFields_ScalarMAX_RES(config);
/// END_GROUP

/// BEGIN_GROUP: BGS_RES, DESCRIPTION: The block Gauss Seidel residuals of the SOLUTION variables.
/// DESCRIPTION: Maximum residual of the density.
AddHistoryOutput("BGS_DENSITY", "bgs[Rho]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the density.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the momentum x-component.
/// DESCRIPTION: BGS residual of the species densities.
for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
AddHistoryOutput("BGS_DENSITY_" + std::to_string(iSpecies), "bgs[Rho_" + std::to_string(iSpecies) + "]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the species density " + std::to_string(iSpecies) + ".", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: BGS residual of the momentum x-component.
AddHistoryOutput("BGS_MOMENTUM-X", "bgs[RhoU]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum x-component.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the momentum y-component.
/// DESCRIPTION: BGS residual of the momentum y-component.
AddHistoryOutput("BGS_MOMENTUM-Y", "bgs[RhoV]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum y-component.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the momentum z-component.
if (nDim == 3) AddHistoryOutput("BGS_MOMENTUM-Z", "bgs[RhoW]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the z-component.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: Maximum residual of the energy.
/// DESCRIPTION: BGS residual of the momentum z-component.
if (nDim == 3) AddHistoryOutput("BGS_MOMENTUM-Z", "bgs[RhoW]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the momentum z-component.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: BGS residual of the energy.
AddHistoryOutput("BGS_ENERGY", "bgs[RhoE]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the energy.", HistoryFieldType::RESIDUAL);
/// DESCRIPTION: BGS residual of the vibrational-electronic energy.
AddHistoryOutput("BGS_ENERGY_VE", "bgs[RhoEve]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the vibrational-electronic energy.", HistoryFieldType::RESIDUAL);
AddHistoryOutputFields_ScalarBGS_RES(config);
/// END_GROUP

Expand Down Expand Up @@ -410,24 +416,32 @@ void CNEMOCompOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSol
SetHistoryOutputValue("RMS_ENERGY", log10(NEMO_solver->GetRes_RMS(nSpecies+3)));
SetHistoryOutputValue("RMS_ENERGY_VE", log10(NEMO_solver->GetRes_RMS(nSpecies+4)));
}
SetHistoryOutputValue("MAX_DENSITY", log10(NEMO_solver->GetRes_Max(0)));
SetHistoryOutputValue("MAX_MOMENTUM-X", log10(NEMO_solver->GetRes_Max(1)));
SetHistoryOutputValue("MAX_MOMENTUM-Y", log10(NEMO_solver->GetRes_Max(2)));
if (nDim == 2)
SetHistoryOutputValue("MAX_ENERGY", log10(NEMO_solver->GetRes_Max(3)));
else {
SetHistoryOutputValue("MAX_MOMENTUM-Z", log10(NEMO_solver->GetRes_Max(3)));
SetHistoryOutputValue("MAX_ENERGY", log10(NEMO_solver->GetRes_Max(4)));
for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
SetHistoryOutputValue("MAX_DENSITY_" + std::to_string(iSpecies), log10(NEMO_solver->GetRes_Max(iSpecies)));

SetHistoryOutputValue("MAX_MOMENTUM-X", log10(NEMO_solver->GetRes_Max(nSpecies)));
SetHistoryOutputValue("MAX_MOMENTUM-Y", log10(NEMO_solver->GetRes_Max(nSpecies+1)));
if (nDim == 2){
SetHistoryOutputValue("MAX_ENERGY", log10(NEMO_solver->GetRes_Max(nSpecies+2)));
SetHistoryOutputValue("MAX_ENERGY_VE", log10(NEMO_solver->GetRes_Max(nSpecies+3)));
} else {
SetHistoryOutputValue("MAX_MOMENTUM-Z", log10(NEMO_solver->GetRes_Max(nSpecies+2)));
SetHistoryOutputValue("MAX_ENERGY", log10(NEMO_solver->GetRes_Max(nSpecies+3)));
SetHistoryOutputValue("MAX_ENERGY_VE", log10(NEMO_solver->GetRes_Max(nSpecies+4)));
}
if (multiZone){
SetHistoryOutputValue("BGS_DENSITY", log10(NEMO_solver->GetRes_BGS(0)));
SetHistoryOutputValue("BGS_MOMENTUM-X", log10(NEMO_solver->GetRes_BGS(1)));
SetHistoryOutputValue("BGS_MOMENTUM-Y", log10(NEMO_solver->GetRes_BGS(2)));
if (nDim == 2)
SetHistoryOutputValue("BGS_ENERGY", log10(NEMO_solver->GetRes_BGS(3)));
else {
SetHistoryOutputValue("BGS_MOMENTUM-Z", log10(NEMO_solver->GetRes_BGS(3)));
SetHistoryOutputValue("BGS_ENERGY", log10(NEMO_solver->GetRes_BGS(4)));
for(iSpecies = 0; iSpecies < nSpecies; iSpecies++)
SetHistoryOutputValue("BGS_DENSITY_" + std::to_string(iSpecies), log10(NEMO_solver->GetRes_BGS(iSpecies)));

SetHistoryOutputValue("BGS_MOMENTUM-X", log10(NEMO_solver->GetRes_BGS(nSpecies)));
SetHistoryOutputValue("BGS_MOMENTUM-Y", log10(NEMO_solver->GetRes_BGS(nSpecies+1)));
if (nDim == 2){
SetHistoryOutputValue("BGS_ENERGY", log10(NEMO_solver->GetRes_BGS(nSpecies+2)));
SetHistoryOutputValue("BGS_ENERGY_VE", log10(NEMO_solver->GetRes_BGS(nSpecies+3)));
} else {
SetHistoryOutputValue("BGS_MOMENTUM-Z", log10(NEMO_solver->GetRes_BGS(nSpecies+2)));
SetHistoryOutputValue("BGS_ENERGY", log10(NEMO_solver->GetRes_BGS(nSpecies+3)));
SetHistoryOutputValue("BGS_ENERGY_VE", log10(NEMO_solver->GetRes_BGS(nSpecies+4)));
}
}

Expand Down
8 changes: 8 additions & 0 deletions SU2_CFD/src/solvers/CSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2048,6 +2048,14 @@ void CSolver::SetResidual_RMS(const CGeometry *geometry, const CConfig *config,
}
}

/*--- Floor the maximum residual like the RMS residual above, so an
* exactly-zero residual (e.g. an inert species at a cold start)
* cannot become a nonfinite log10 in screen or history output. ---*/

for (unsigned short iVar = 0; iVar < nVar; iVar++) {
Residual_Max[iVar] = max(EPS*EPS, Residual_Max[iVar]);
}

}
END_SU2_OMP_SAFE_GLOBAL_ACCESS
}
Expand Down
214 changes: 214 additions & 0 deletions UnitTests/SU2_CFD/output/CNEMOCompOutput_tests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,214 @@
/*!
* \file CNEMOCompOutput_tests.cpp
* \brief Unit tests for NEMO residual history output.
* \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 <array>
#include <cmath>
#include <memory>
#include <sstream>
#include <string>
#include <vector>

#include "../../../Common/include/CConfig.hpp"
#include "../../../Common/include/geometry/CGeometry.hpp"
#include "../../../SU2_CFD/include/output/CNEMOCompOutput.hpp"
#include "../../../SU2_CFD/include/solvers/CSolver.hpp"
#include "../../../SU2_CFD/include/variables/CVariable.hpp"

namespace {

std::unique_ptr<CConfig> MakeNEMOConfig(bool multizone) {
std::stringstream options;
options << "SOLVER= NEMO_EULER\n"
<< "GAS_MODEL= AIR-5\n"
<< "GAS_COMPOSITION= (0.77, 0.23, 0.0, 0.0, 0.0)\n"
<< "FLUID_MODEL= SU2_NONEQ\n"
<< "MATH_PROBLEM= DIRECT\n"
<< "MACH_NUMBER= 5.0\n"
<< "FREESTREAM_PRESSURE= 101325.0\n"
<< "FREESTREAM_TEMPERATURE= 288.15\n"
<< "FREESTREAM_TEMPERATURE_VE= 288.15\n"
<< "COMM_LEVEL= MINIMAL\n";
auto config = std::make_unique<CConfig>(options, SU2_COMPONENT::SU2_CFD, false);
config->SetMultizone_Problem(multizone);
return config;
}

class CTestGeometry final : public CGeometry {
public:
explicit CTestGeometry(unsigned short dimension) {
nDim = dimension;
nPoint = nPointDomain = Global_nPoint = Global_nPointDomain = 1;
MGLevel = MESH_0;
}
};

class CTestSolver final : public CSolver {
private:
CVariable variables;

CVariable* GetBaseClassPointerToNodes() override { return &variables; }

public:
CTestSolver(const CConfig* config, unsigned short dimension)
: variables(0, dimension, config->GetnSpecies() + dimension + 2, config) {
nDim = dimension;
nVar = config->GetnSpecies() + dimension + 2;
nPoint = nPointDomain = 1;

Residual_RMS.resize(nVar);
Residual_Max.resize(nVar);
Residual_BGS.resize(nVar);
Residual_Max_BGS.resize(nVar);
Point_Max.resize(nVar);
Point_Max_BGS.resize(nVar);
Point_Max_Coord.resize(nVar, nDim);
Point_Max_Coord_BGS.resize(nVar, nDim);

SetBaseClassPointerToNodes();
SetCFL_Local_Stats(1.0);
SetResLinSolver(1.0);
}

void SeedResidualSums() {
for (unsigned short iVar = 0; iVar < nVar; ++iVar) {
const auto rmsExponent = static_cast<int>(iVar) + 1;
const auto maxExponent = static_cast<int>(iVar) + 22;
const auto bgsExponent = static_cast<int>(iVar) + 12;
Residual_RMS[iVar] = std::pow(10.0, -2.0 * rmsExponent);
Residual_Max[iVar] = std::pow(10.0, -maxExponent);
Residual_BGS[iVar] = std::pow(10.0, -2.0 * bgsExponent);
}

/* An exactly-zero species maximum exercises the finite output floor. */
Residual_Max[2] = 0.0;
}
};

std::vector<std::string> ExpectedResidualFields(const std::string& prefix, unsigned short dimension) {
std::vector<std::string> fields;
for (unsigned short iSpecies = 0; iSpecies < 5; ++iSpecies) {
fields.push_back(prefix + "_DENSITY_" + std::to_string(iSpecies));
}
fields.push_back(prefix + "_MOMENTUM-X");
fields.push_back(prefix + "_MOMENTUM-Y");
if (dimension == 3) fields.push_back(prefix + "_MOMENTUM-Z");
fields.push_back(prefix + "_ENERGY");
fields.push_back(prefix + "_ENERGY_VE");
return fields;
}

std::vector<std::string> ExpectedFieldNames(const std::string& prefix, unsigned short dimension) {
std::vector<std::string> names;
for (unsigned short iSpecies = 0; iSpecies < 5; ++iSpecies) {
names.push_back(prefix + "[Rho_" + std::to_string(iSpecies) + "]");
}
names.push_back(prefix + "[RhoU]");
names.push_back(prefix + "[RhoV]");
if (dimension == 3) names.push_back(prefix + "[RhoW]");
names.push_back(prefix + "[RhoE]");
names.push_back(prefix + "[RhoEve]");
return names;
}

void CheckRegisteredGroup(CNEMOCompOutput& output, const std::string& group, const std::string& prefix,
unsigned short dimension) {
const auto expectedFields = ExpectedResidualFields(group == "MAX_RES" ? "MAX" : "BGS", dimension);
const auto expectedNames = ExpectedFieldNames(prefix, dimension);
const auto registered = output.GetHistoryGroup(group);

REQUIRE(registered.size() == expectedFields.size());
const auto& fields = output.GetHistoryFields();
for (std::size_t i = 0; i < expectedFields.size(); ++i) {
INFO("history field " << expectedFields[i]);
REQUIRE(fields.count(expectedFields[i]) == 1);
CHECK(fields.at(expectedFields[i]).fieldName == expectedNames[i]);
CHECK(registered[i].fieldName == expectedNames[i]);
}
}

void CheckLoadedResiduals(unsigned short dimension) {
auto config = MakeNEMOConfig(true);
CTestGeometry geometry(dimension);
CTestSolver flow(config.get(), dimension);
flow.SeedResidualSums();

/* Exercise the actual reductions used before MAX and multizone BGS output. */
flow.SetResidual_RMS(&geometry, config.get());
flow.SetResidual_BGS(&geometry, config.get());

CNEMOCompOutput output(config.get(), dimension);
output.SetHistoryOutputFields(config.get());
std::array<CSolver*, MAX_SOLS> solvers{};
solvers[FLOW_SOL] = &flow;
solvers[MESH_SOL] = &flow;
output.LoadHistoryData(config.get(), &geometry, solvers.data());

const unsigned short nVar = config->GetnSpecies() + dimension + 2;
const auto maxFields = ExpectedResidualFields("MAX", dimension);
const auto bgsFields = ExpectedResidualFields("BGS", dimension);
REQUIRE(maxFields.size() == nVar);
REQUIRE(bgsFields.size() == nVar);

for (unsigned short iVar = 0; iVar < nVar; ++iVar) {
const auto maxValue = SU2_TYPE::GetValue(output.GetHistoryFieldValue(maxFields[iVar]));
const auto bgsValue = SU2_TYPE::GetValue(output.GetHistoryFieldValue(bgsFields[iVar]));
INFO("residual variable " << iVar);
CHECK(std::isfinite(maxValue));
CHECK(maxValue == Approx(iVar == 2 ? -32.0 : -(static_cast<double>(iVar) + 22.0)));
CHECK(bgsValue == Approx(-(static_cast<double>(iVar) + 12.0)));
}
}

} // namespace

TEST_CASE("NEMO MAX and BGS history fields use the species-first layout", "[NEMO][Output]") {
auto config = MakeNEMOConfig(false);

SECTION("2D") {
CNEMOCompOutput output(config.get(), 2);
output.SetHistoryOutputFields(config.get());
CheckRegisteredGroup(output, "MAX_RES", "max", 2);
CheckRegisteredGroup(output, "BGS_RES", "bgs", 2);
CHECK(output.GetHistoryFields().count("MAX_DENSITY") == 0);
CHECK(output.GetHistoryFields().count("BGS_DENSITY") == 0);
}

SECTION("3D") {
CNEMOCompOutput output(config.get(), 3);
output.SetHistoryOutputFields(config.get());
CheckRegisteredGroup(output, "MAX_RES", "max", 3);
CheckRegisteredGroup(output, "BGS_RES", "bgs", 3);
CHECK(output.GetHistoryFields().count("MAX_DENSITY") == 0);
CHECK(output.GetHistoryFields().count("BGS_DENSITY") == 0);
}
}

TEST_CASE("NEMO history loads finite MAX and multizone BGS residuals by species-first index", "[NEMO][Output]") {
SECTION("2D") { CheckLoadedResiduals(2); }
SECTION("3D") { CheckLoadedResiduals(3); }
}
1 change: 1 addition & 0 deletions UnitTests/meson.build
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp',
'Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp',
'SU2_CFD/numerics/CNumerics_tests.cpp',
'SU2_CFD/fluid/CFluidModel_tests.cpp',
'SU2_CFD/output/CNEMOCompOutput_tests.cpp',
'SU2_CFD/gradients.cpp',
'SU2_CFD/windowing.cpp',
'Common/toolboxes/random_toolbox_tests.cpp',
Expand Down