Skip to content

Commit 0840453

Browse files
ScSteffenScSteffen
authored
Bugfixes in HPC SN Solver (#61)
* fix lattice perimeter definitions * fix max outflow MPI * fix output of varAbsorption * fixed rms bug with mpi exec * fix linspace vulnerability * fixed MPI flag * recomputation of scalar flux missing after heun averaging * more guardrails around greenprobing * made the probing blocks move with the center of the capsule * update restart sim time * fixed validation tests
1 parent f02e8c0 commit 0840453

4 files changed

Lines changed: 75 additions & 55 deletions

File tree

‎src/solvers/snsolver_hpc.cpp‎

Lines changed: 62 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include "quadratures/quadraturebase.hpp"
1010
#include "solvers/snsolver_hpc.hpp"
1111
#include "toolboxes/textprocessingtoolbox.hpp"
12+
#include <cassert>
1213

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

292294
std::chrono::duration<double> duration;
@@ -311,6 +313,25 @@ void SNSolverHPC::Solve() {
311313
_sol[Idx2D( idx_cell, idx_sys, _localNSys )] ); // Solution averaging with HEUN
312314
}
313315
}
316+
317+
// Keep scalar flux consistent with the final RK2-averaged solution used in postprocessing.
318+
std::vector<double> temp_scalarFluxRK( _nCells, 0.0 );
319+
#pragma omp parallel for
320+
for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
321+
double localScalarFlux = 0.0;
322+
#pragma omp simd reduction( + : localScalarFlux )
323+
for( unsigned long idx_sys = 0; idx_sys < _localNSys; idx_sys++ ) {
324+
localScalarFlux += _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys];
325+
}
326+
temp_scalarFluxRK[idx_cell] = localScalarFlux;
327+
}
328+
#ifdef IMPORT_MPI
329+
MPI_Barrier( MPI_COMM_WORLD );
330+
MPI_Allreduce( temp_scalarFluxRK.data(), _scalarFlux.data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
331+
MPI_Barrier( MPI_COMM_WORLD );
332+
#else
333+
_scalarFlux = temp_scalarFluxRK;
334+
#endif
314335
}
315336
else {
316337

@@ -334,12 +355,12 @@ void SNSolverHPC::Solve() {
334355
PrintScreenOutput( iter );
335356
PrintHistoryOutput( iter );
336357
}
337-
#ifdef BUILD_MPI
358+
#ifdef IMPORT_MPI
338359
MPI_Barrier( MPI_COMM_WORLD );
339360
#endif
340361

341362
PrintVolumeOutput( iter );
342-
#ifdef BUILD_MPI
363+
#ifdef IMPORT_MPI
343364
MPI_Barrier( MPI_COMM_WORLD );
344365
#endif
345366
}
@@ -536,8 +557,9 @@ void SNSolverHPC::FVMUpdate() {
536557
_mass = 0.0;
537558
_rmsFlux = 0.0;
538559
std::vector<double> temp_scalarFlux( _nCells ); // for MPI allreduce
560+
std::vector<double> prev_scalarFlux( _scalarFlux );
539561

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

543565
#pragma omp simd
@@ -556,17 +578,27 @@ void SNSolverHPC::FVMUpdate() {
556578
localScalarFlux += _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys];
557579
}
558580
_mass += localScalarFlux * _areas[idx_cell];
559-
_rmsFlux += ( localScalarFlux - _scalarFlux[idx_cell] ) * ( localScalarFlux - _scalarFlux[idx_cell] );
560581
temp_scalarFlux[idx_cell] = localScalarFlux; // set flux
561582
}
562583
// MPI Allreduce: _scalarFlux
563584
#ifdef IMPORT_MPI
564585
MPI_Barrier( MPI_COMM_WORLD );
565586
MPI_Allreduce( temp_scalarFlux.data(), _scalarFlux.data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
566587
MPI_Barrier( MPI_COMM_WORLD );
588+
if( _rank == 0 ) {
589+
for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
590+
double diff = _scalarFlux[idx_cell] - prev_scalarFlux[idx_cell];
591+
_rmsFlux += diff * diff;
592+
}
593+
}
567594
#endif
568595
#ifndef IMPORT_MPI
569596
_scalarFlux = temp_scalarFlux;
597+
#pragma omp parallel for reduction( + : _rmsFlux )
598+
for( unsigned long idx_cell = 0; idx_cell < _nCells; ++idx_cell ) {
599+
double diff = _scalarFlux[idx_cell] - prev_scalarFlux[idx_cell];
600+
_rmsFlux += diff * diff;
601+
}
570602
#endif
571603
}
572604

@@ -714,6 +746,7 @@ void SNSolverHPC::IterPostprocessing() {
714746
double tmp_curScalarOutflow = 0.0;
715747
double tmp_curScalarOutflowPeri1 = 0.0;
716748
double tmp_curScalarOutflowPeri2 = 0.0;
749+
double tmp_curMaxOrdinateOutflow = 0.0;
717750
double tmp_mass = 0.0;
718751
double tmp_rmsFlux = 0.0;
719752
MPI_Barrier( MPI_COMM_WORLD );
@@ -726,6 +759,9 @@ void SNSolverHPC::IterPostprocessing() {
726759
MPI_Allreduce( &_curScalarOutflowPeri2, &tmp_curScalarOutflowPeri2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
727760
_curScalarOutflowPeri2 = tmp_curScalarOutflowPeri2;
728761
MPI_Barrier( MPI_COMM_WORLD );
762+
MPI_Allreduce( &_curMaxOrdinateOutflow, &tmp_curMaxOrdinateOutflow, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
763+
_curMaxOrdinateOutflow = tmp_curMaxOrdinateOutflow;
764+
MPI_Barrier( MPI_COMM_WORLD );
729765
MPI_Allreduce( &_mass, &tmp_mass, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
730766
_mass = tmp_mass;
731767
MPI_Barrier( MPI_COMM_WORLD );
@@ -973,7 +1009,7 @@ void SNSolverHPC::WriteScalarOutput( unsigned idx_iter ) {
9731009
}
9741010
idx_field--;
9751011
break;
976-
case VAR_ABSORPTION_GREEN: _screenOutputFields[idx_field] = _absorptionValsBlocksGreen[0]; break;
1012+
case VAR_ABSORPTION_GREEN: _screenOutputFields[idx_field] = _varAbsorptionHohlraumGreen; break;
9771013
default: ErrorMessages::Error( "Screen output group not defined!", CURRENT_FUNCTION ); break;
9781014
}
9791015
}
@@ -1318,7 +1354,7 @@ void SNSolverHPC::WriteVolumeOutput( unsigned idx_iter ) {
13181354
_quadPts[Idx2D( idx_sys, 1, _nDim )] * _sol[Idx2D( idx_cell, idx_sys, _localNSys )] * _quadWeights[idx_sys];
13191355
}
13201356
}
1321-
#ifdef BUILD_MPI
1357+
#ifdef IMPORT_MPI
13221358
MPI_Barrier( MPI_COMM_WORLD );
13231359
MPI_Allreduce(
13241360
_outputFields[idx_group][0].data(), _outputFields[idx_group][0].data(), _nCells, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
@@ -1455,28 +1491,8 @@ void SNSolverHPC::SetGhostCells() {
14551491

14561492
void SNSolverHPC::SetProbingCellsLineGreen() {
14571493

1458-
// if( _settings->GetProblemName() == PROBLEM_QuarterHohlraum ) {
1459-
// double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] );
1460-
// double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] );
1461-
//
1462-
// // double dx = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen );
1463-
//
1464-
// unsigned nHorizontalProbingCells =
1465-
// (unsigned)std::ceil( _nProbingCellsLineGreen * ( horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ) ) );
1466-
// unsigned nVerticalProbingCells = _nProbingCellsLineGreen - nHorizontalProbingCells;
1467-
//
1468-
// _probingCellsLineGreen = std::vector<unsigned>( _nProbingCellsLineGreen );
1469-
//
1470-
// // Sample points on each side of the rectangle
1471-
// std::vector<unsigned> side3 = linspace2D( _cornerLowerRightGreen, _cornerUpperRightGreen, nVerticalProbingCells );
1472-
// std::vector<unsigned> side4 = linspace2D( _cornerUpperRightGreen, _cornerUpperLeftGreen, nHorizontalProbingCells );
1473-
//
1474-
// // Combine the points from each side
1475-
// _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side3.begin(), side3.end() );
1476-
// _probingCellsLineGreen.insert( _probingCellsLineGreen.end(), side4.begin(), side4.end() );
1477-
// }
1478-
// else
14791494
if( _settings->GetProblemName() == PROBLEM_SymmetricHohlraum ) {
1495+
assert( _nProbingCellsLineGreen % 2 == 0 );
14801496

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

14911507
unsigned nHorizontalProbingCells = (unsigned)std::ceil( _nProbingCellsLineGreen / 2 * pt_ratio_h );
14921508
unsigned nVerticalProbingCells = _nProbingCellsLineGreen / 2 - nHorizontalProbingCells;
1509+
assert( nHorizontalProbingCells > 1 );
1510+
assert( nVerticalProbingCells > 1 );
14931511

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

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

15191537
double block_size = 0.05;
1538+
const double cx = _centerGreen[0];
1539+
const double cy = _centerGreen[1];
15201540

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

15241544
// Top row
1525-
double x1 = -0.2 + i * block_size;
1526-
double y1 = 0.4;
1545+
double x1 = -0.2 + cx + i * block_size;
1546+
double y1 = 0.4 + cy;
15271547
double x2 = x1 + block_size;
15281548
double y2 = y1 - block_size;
15291549

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

15391559
for( int j = 0; j < 14; ++j ) { // 14 blocks in the y-direction (vertical)
15401560
// right column double x1 = 0.15;
1541-
double x1 = 0.15;
1542-
double y1 = 0.35 - j * block_size;
1561+
double x1 = 0.15 + cx;
1562+
double y1 = 0.35 + cy - j * block_size;
15431563
double x2 = x1 + block_size;
15441564
double y2 = y1 - block_size;
15451565

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

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

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

15751595
// left column
1576-
double x1 = -0.2;
1577-
double y1 = -0.3 + j * block_size;
1596+
double x1 = -0.2 + cx;
1597+
double y1 = -0.3 + cy + j * block_size;
15781598
double x2 = x1 + block_size;
15791599
double y2 = y1 - block_size;
15801600

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

16301650
std::vector<unsigned> result;
1651+
assert( num_points > 1 );
16311652
result.resize( num_points );
16321653
double stepX = ( end[0] - start[0] ) / ( num_points - 1 );
16331654
double stepY = ( end[1] - start[1] ) / ( num_points - 1 );
@@ -1662,23 +1683,22 @@ void SNSolverHPC::ComputeCellsPerimeterLattice() {
16621683
continue; // Skip boundary - ghost cells
16631684
}
16641685

1665-
if( abs( ( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_1 && abs( cellMids[idx_cell][0] ) < l_1 ) ||
1666-
abs( ( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_1 && abs( cellMids[idx_cell][1] ) < l_1 ) ) {
1686+
if( ( abs( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_1 && abs( cellMids[idx_cell][0] ) < l_1 ) ||
1687+
( abs( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_1 && abs( cellMids[idx_cell][1] ) < l_1 ) ) {
16671688
// neighbor is outside perimeter
16681689
_cellsLatticePerimeter1[idx_cell].push_back( idx_nbr );
16691690
_isPerimeterLatticeCell1[idx_cell] = true;
16701691
}
16711692
}
16721693
}
1673-
if( abs( cellMids[idx_cell][0] ) < l_2 && abs( cellMids[idx_cell][1] ) < l_2 && abs( cellMids[idx_cell][0] ) > l_1 &&
1674-
abs( cellMids[idx_cell][1] ) > l_1 ) {
1694+
else if( abs( cellMids[idx_cell][0] ) < l_2 && abs( cellMids[idx_cell][1] ) < l_2 ) {
16751695
// Cell is within perimeter
16761696
for( unsigned idx_nbr = 0; idx_nbr < _mesh->GetNumNodesPerCell(); ++idx_nbr ) {
16771697
if( neigbors[idx_cell][idx_nbr] == _mesh->GetNumCells() ) {
16781698
continue; // Skip boundary - ghost cells
16791699
}
1680-
if( abs( ( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_2 && abs( cellMids[idx_cell][0] ) < l_2 ) ||
1681-
abs( ( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_2 && abs( cellMids[idx_cell][1] ) < l_2 ) ) {
1700+
if( ( abs( cellMids[neigbors[idx_cell][idx_nbr]][0] ) > l_2 && abs( cellMids[idx_cell][0] ) < l_2 ) ||
1701+
( abs( cellMids[neigbors[idx_cell][idx_nbr]][1] ) > l_2 && abs( cellMids[idx_cell][1] ) < l_2 ) ) {
16821702
// neighbor is outside perimeter
16831703
_cellsLatticePerimeter2[idx_cell].push_back( idx_nbr );
16841704
_isPerimeterLatticeCell2[idx_cell] = true;

‎tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200_csv_reference‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
11
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
2-
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
3-
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
4-
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
2+
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
3+
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
4+
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

0 commit comments

Comments
 (0)