Skip to content

Bugfixes in HPC SN Solver #61

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

Merged
ScSteffen merged 12 commits into from
Feb 12, 2026
Merged

Bugfixes in HPC SN Solver #61

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
be9c48a
fix lattice perimeter definitions
ScSteffen Feb 10, 2026
a8f6dca
fix max outflow MPI
ScSteffen Feb 10, 2026
2044aa4
fix output of varAbsorption
ScSteffen Feb 10, 2026
73d98b8
fixed rms bug with mpi exec
ScSteffen Feb 10, 2026
88c19a5
fix linspace vulnerability
ScSteffen Feb 10, 2026
ecbab2a
fixed MPI flag
ScSteffen Feb 10, 2026
8044095
recomputation of scalar flux missing after heun averaging
ScSteffen Feb 10, 2026
8e6584e
more guardrails around greenprobing
ScSteffen Feb 10, 2026
189d075
made the probing blocks move with the center of the capsule
ScSteffen Feb 10, 2026
67b7f16
update restart sim time
ScSteffen Feb 10, 2026
116423b
Merge branch 'master' into bugfix_lattice_perimeter
ScSteffen Feb 11, 2026
c823526
fixed validation tests
ScSteffen Feb 11, 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
104 changes: 62 additions & 42 deletions src/solvers/snsolver_hpc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#include "quadratures/quadraturebase.hpp"
#include "solvers/snsolver_hpc.hpp"
#include "toolboxes/textprocessingtoolbox.hpp"
#include <cassert>

SNSolverHPC::SNSolverHPC( Config* settings ) {
#ifdef IMPORT_MPI
Expand Down Expand Up @@ -286,7 +287,8 @@ void SNSolverHPC::Solve() {
PrepareVolumeOutput();
DrawPreSolverOutput();
}
_curSimTime = 0.0;
// On restart, continue simulation time from the loaded iteration index.
_curSimTime = static_cast<double>( _idx_start_iter ) * _dT;
auto start = std::chrono::high_resolution_clock::now(); // Start timing

std::chrono::duration<double> duration;
Expand All @@ -311,6 +313,25 @@ void SNSolverHPC::Solve() {
_sol[Idx2D( idx_cell, idx_sys, _localNSys )] ); // Solution averaging with HEUN
}
}

// Keep scalar flux consistent with the final RK2-averaged solution used in postprocessing.
std::vector<double> temp_scalarFluxRK( _nCells, 0.0 );
#pragma omp parallel for
for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
double localScalarFlux = 0.0;
#pragma omp simd reduction( + : localScalarFlux )
for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) {
localScalarFlux += _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys];
}
temp_scalarFluxRK[idx_cell] = localScalarFlux;
}
#ifdef IMPORT_MPI
MPI_Barrier( MPI_COMM_WORLD );
MPI_Allreduce( temp_scalarFluxRK.data(), _scalarFlux.data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
MPI_Barrier( MPI_COMM_WORLD );
#else
_scalarFlux = temp_scalarFluxRK;
#endif
}
else {

Expand All @@ -334,12 +355,12 @@ void SNSolverHPC::Solve() {
PrintScreenOutput( iter );
PrintHistoryOutput( iter );
}
#ifdef BUILD_MPI
#ifdef IMPORT_MPI
MPI_Barrier( MPI_COMM_WORLD );
#endif

PrintVolumeOutput( iter );
#ifdef BUILD_MPI
#ifdef IMPORT_MPI
MPI_Barrier( MPI_COMM_WORLD );
#endif
}
Expand Down Expand Up @@ -536,8 +557,9 @@ void SNSolverHPC::FVMUpdate() {
_mass = 0.0;
_rmsFlux = 0.0;
std::vector<double> temp_scalarFlux( _nCells ); // for MPI allreduce
std::vector<double> prev_scalarFlux( _scalarFlux );

#pragma omp parallel for reduction( + : _mass, _rmsFlux )
#pragma omp parallel for reduction( + : _mass )
for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {

#pragma omp simd
Expand All @@ -556,17 +578,27 @@ void SNSolverHPC::FVMUpdate() {
localScalarFlux += _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys];
}
_mass += localScalarFlux * _areas[idx_cell];
_rmsFlux += ( localScalarFlux - _scalarFlux[idx_cell] ) * ( localScalarFlux - _scalarFlux[idx_cell] );
temp_scalarFlux[idx_cell] = localScalarFlux; // set flux
}
// MPI Allreduce: _scalarFlux
#ifdef IMPORT_MPI
MPI_Barrier( MPI_COMM_WORLD );
MPI_Allreduce( temp_scalarFlux.data(), _scalarFlux.data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
MPI_Barrier( MPI_COMM_WORLD );
if( _rank == 0 ) {
for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
double diff = _scalarFlux[idx_cell] - prev_scalarFlux[idx_cell];
_rmsFlux += diff * diff;
}
}
#endif
#ifndef IMPORT_MPI
_scalarFlux = temp_scalarFlux;
#pragma omp parallel for reduction( + : _rmsFlux )
for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
double diff = _scalarFlux[idx_cell] - prev_scalarFlux[idx_cell];
_rmsFlux += diff * diff;
}
#endif
}

Expand Down Expand Up @@ -714,6 +746,7 @@ void SNSolverHPC::IterPostprocessing() {
double tmp_curScalarOutflow = 0.0;
double tmp_curScalarOutflowPeri1 = 0.0;
double tmp_curScalarOutflowPeri2 = 0.0;
double tmp_curMaxOrdinateOutflow = 0.0;
double tmp_mass = 0.0;
double tmp_rmsFlux = 0.0;
MPI_Barrier( MPI_COMM_WORLD );
Expand All @@ -726,6 +759,9 @@ void SNSolverHPC::IterPostprocessing() {
MPI_Allreduce( &_curScalarOutflowPeri2, &tmp_curScalarOutflowPeri2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
_curScalarOutflowPeri2 = tmp_curScalarOutflowPeri2;
MPI_Barrier( MPI_COMM_WORLD );
MPI_Allreduce( &_curMaxOrdinateOutflow, &tmp_curMaxOrdinateOutflow, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
_curMaxOrdinateOutflow = tmp_curMaxOrdinateOutflow;
MPI_Barrier( MPI_COMM_WORLD );
MPI_Allreduce( &_mass, &tmp_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
_mass = tmp_mass;
MPI_Barrier( MPI_COMM_WORLD );
Expand Down Expand Up @@ -973,7 +1009,7 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) {
}
idx_field--;
break;
case VAR_ABSORPTION_GREEN: _screenOutputFields[idx_field] = _absorptionValsBlocksGreen[0]; break;
case VAR_ABSORPTION_GREEN: _screenOutputFields[idx_field] = _varAbsorptionHohlraumGreen; break;
default: ErrorMessages::Error( "Screen output group not defined!", CURRENT_FUNCTION ); break;
}
}
Expand Down Expand Up @@ -1318,7 +1354,7 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) {
_quadPts[Idx2D( idx_sys, 1, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys];
}
}
#ifdef BUILD_MPI
#ifdef IMPORT_MPI
MPI_Barrier( MPI_COMM_WORLD );
MPI_Allreduce(
_outputFields[idx_group][0].data(), _outputFields[idx_group][0].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
Expand Down Expand Up @@ -1455,28 +1491,8 @@ void SNSolverHPC::SetGhostCells() {

void SNSolverHPC::SetProbingCellsLineGreen() {

// if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) {
// double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] );
// double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] );
//
// // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen );
//
// unsigned nHorizontalProbingCells =
// (unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) );
// unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells;
//
// _probingCellsLineGreen = std::vector<unsigned>( _nProbingCellsLineGreen );
//
// // Sample points on each side of the rectangle
// std::vector<unsigned> side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells );
// std::vector<unsigned> side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, nHorizontalProbingCells );
//
// // Combine the points from each side
// _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() );
// _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() );
// }
// else
if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) {
assert( _nProbingCellsLineGreen % 2 == 0 );

std::vector<double> p1 = { _cornerUpperLeftGreen[0] + _thicknessGreen / 2.0, _cornerUpperLeftGreen[1] - _thicknessGreen / 2.0 };
std::vector<double> p2 = { _cornerUpperRightGreen[0] - _thicknessGreen / 2.0, _cornerUpperRightGreen[1] - _thicknessGreen / 2.0 };
Expand All @@ -1490,6 +1506,8 @@ void SNSolverHPC::SetProbingCellsLineGreen() {

unsigned nHorizontalProbingCells = (unsigned)std::ceil( _nProbingCellsLineGreen / 2 * pt_ratio_h );
unsigned nVerticalProbingCells = _nProbingCellsLineGreen / 2 - nHorizontalProbingCells;
assert( nHorizontalProbingCells > 1 );
assert( nVerticalProbingCells > 1 );

_probingCellsLineGreen = std::vector<unsigned>( 2 * ( nVerticalProbingCells + nHorizontalProbingCells ) );

Expand Down Expand Up @@ -1517,13 +1535,15 @@ void SNSolverHPC::SetProbingCellsLineGreen() {
std::vector<std::vector<std::vector<double>>> block_corners;

double block_size = 0.05;
const double cx = _centerGreen[0];
const double cy = _centerGreen[1];

// Loop to fill the blocks
for( int i = 0; i < 8; ++i ) { // 8 blocks in the x-direction (horizontal) (upper) (left to right)

// Top row
double x1 = -0.2 + i * block_size;
double y1 = 0.4;
double x1 = -0.2 + cx + i * block_size;
double y1 = 0.4 + cy;
double x2 = x1 + block_size;
double y2 = y1 - block_size;

Expand All @@ -1538,8 +1558,8 @@ void SNSolverHPC::SetProbingCellsLineGreen() {

for( int j = 0; j < 14; ++j ) { // 14 blocks in the y-direction (vertical)
// right column double x1 = 0.15;
double x1 = 0.15;
double y1 = 0.35 - j * block_size;
double x1 = 0.15 + cx;
double y1 = 0.35 + cy - j * block_size;
double x2 = x1 + block_size;
double y2 = y1 - block_size;

Expand All @@ -1556,8 +1576,8 @@ void SNSolverHPC::SetProbingCellsLineGreen() {

for( int i = 0; i < 8; ++i ) { // 8 blocks in the x-direction (horizontal) (lower) (right to left)
// bottom row
double x1 = 0.15 - i * block_size;
double y1 = -0.35;
double x1 = 0.15 + cx - i * block_size;
double y1 = -0.35 + cy;
double x2 = x1 + block_size;
double y2 = y1 - block_size;

Expand All @@ -1573,8 +1593,8 @@ void SNSolverHPC::SetProbingCellsLineGreen() {
for( int j = 0; j < 14; ++j ) { // 14 blocks in the y-direction (vertical) (down to up)

// left column
double x1 = -0.2;
double y1 = -0.3 + j * block_size;
double x1 = -0.2 + cx;
double y1 = -0.3 + cy + j * block_size;
double x2 = x1 + block_size;
double y2 = y1 - block_size;

Expand Down Expand Up @@ -1628,6 +1648,7 @@ std::vector<unsigned> SNSolverHPC::linspace2D( const std::vector<double>& start,
*/

std::vector<unsigned> result;
assert( num_points > 1 );
result.resize( num_points );
double stepX = ( end[0] - start[0] ) / ( num_points - 1 );
double stepY = ( end[1] - start[1] ) / ( num_points - 1 );
Expand Down Expand Up @@ -1662,23 +1683,22 @@ void SNSolverHPC::ComputeCellsPerimeterLattice() {
continue; // Skip boundary - ghost cells
}

if( abs( ( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_1 && abs( cellMids[idx_cell][0] ) < l_1 ) ||
abs( ( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_1 && abs( cellMids[idx_cell][1] ) < l_1 ) ) {
if( ( abs( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_1 && abs( cellMids[idx_cell][0] ) < l_1 ) ||
( abs( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_1 && abs( cellMids[idx_cell][1] ) < l_1 ) ) {
// neighbor is outside perimeter
_cellsLatticePerimeter1[idx_cell].push_back( idx_nbr );
_isPerimeterLatticeCell1[idx_cell] = true;
}
}
}
if( abs( cellMids[idx_cell][0] ) < l_2 && abs( cellMids[idx_cell][1] ) < l_2 && abs( cellMids[idx_cell][0] ) > l_1 &&
abs( cellMids[idx_cell][1] ) > l_1 ) {
else if( abs( cellMids[idx_cell][0] ) < l_2 && abs( cellMids[idx_cell][1] ) < l_2 ) {
// Cell is within perimeter
for( unsigned idx_nbr = 0; idx_nbr < _mesh->GetNumNodesPerCell(); ++idx_nbr ) {
if( neigbors[idx_cell][idx_nbr] == _mesh->GetNumCells() ) {
continue; // Skip boundary - ghost cells
}
if( abs( ( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_2 && abs( cellMids[idx_cell][0] ) < l_2 ) ||
abs( ( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_2 && abs( cellMids[idx_cell][1] ) < l_2 ) ) {
if( ( abs( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_2 && abs( cellMids[idx_cell][0] ) < l_2 ) ||
( abs( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_2 && abs( cellMids[idx_cell][1] ) < l_2 ) ) {
// neighbor is outside perimeter
_cellsLatticePerimeter2[idx_cell].push_back( idx_nbr );
_isPerimeterLatticeCell2[idx_cell] = true;
Expand Down
6 changes: 3 additions & 3 deletions tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200_csv_reference
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
2026-02-10 15:35:39.595600 ,Iter,Sim_time,Wall_time_[s],Mass,RMS_flux,VTK_out,CSV_out,Cur_outflow,Total_outflow,Cur_outflow_P1,Total_outflow_P1,Cur_outflow_P2,Total_outflow_P2,Max_outflow,Cur_absorption,Total_absorption,Max_absorption
2026-02-10 15:35:39.658632 ,0.000000,2.474874E-01,6.294604E-02,6.382169E+00,4.127129E+00,0.000000E+00,1.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00
2026-02-10 15:35:39.721179 ,1.000000,4.949747E-01,1.255209E-01,1.039716E+01,3.385174E+00,0.000000E+00,1.000000E+00,0.000000E+00,0.000000E+00,1.786912E-01,4.422380E-02,0.000000E+00,0.000000E+00,0.000000E+00,3.432122E-01,8.494068E-02,5.680457E-02
2026-02-10 15:35:39.781177 ,2.000000,5.000000E-01,1.855212E-01,6.939972E+00,6.744329E-02,1.000000E+00,1.000000E+00,0.000000E+00,0.000000E+00,1.850545E-01,4.515375E-02,0.000000E+00,0.000000E+00,0.000000E+00,1.885307E-01,8.588809E-02,5.680457E-02
2026-02-11 22:02:50.261192 ,0.000000,2.474874E-01,1.064435E-01,6.382169E+00,4.127129E+00,0.000000E+00,1.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00
2026-02-11 22:02:50.375590 ,1.000000,4.949747E-01,2.208639E-01,9.623368E+00,3.464680E+00,0.000000E+00,1.000000E+00,0.000000E+00,0.000000E+00,3.459946E-01,8.562928E-02,0.000000E+00,0.000000E+00,0.000000E+00,1.638836E-01,4.055913E-02,2.562927E-02
2026-02-11 22:02:50.483330 ,2.000000,5.000000E-01,3.285845E-01,6.535954E+00,6.883754E-02,1.000000E+00,1.000000E+00,0.000000E+00,0.000000E+00,3.570196E-01,8.742340E-02,0.000000E+00,0.000000E+00,0.000000E+00,1.715809E-01,4.142137E-02,2.679525E-02
Loading