• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

CSMMLab / KiT-RT / #117

11 Feb 2026 05:12PM UTC coverage: 57.119% (+10.5%) from 46.626%
#117

push

travis-ci

web-flow
Add comprehensive unit tests for core modules (#60)

* Add 7 new test files to increase unit test coverage

Cover previously untested areas: text processing utilities, slope
limiters, entropy classes (quadratic + Maxwell-Boltzmann), upwind flux
overloads (FluxXZ, Flux1D, FluxVanLeer, vectorized), IO round-trip
(connectivity + restart solution), problem geometry (checkerboard,
lattice, line source analytical solution), expanded config parameter
parsing, and additional optimizer scenarios.

https://claude.ai/code/session_01VcLByz1wmHxA1Wev5PwVGL

* fix memory leaks in unit tests

* modify tests for memory safety

* fix memory leak

* add reference files for validation tests

* try to fix segfault

* remove compiler warnings

* attempt to fix segfaults in tests

* added cfg files for sn validation tests

* update workflow

---------

Co-authored-by: Claude <noreply@anthropic.com>

12 of 14 new or added lines in 5 files covered. (85.71%)

21 existing lines in 1 file now uncovered.

5448 of 9538 relevant lines covered (57.12%)

31303.3 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

78.85
/src/common/io.cpp
1
/*!
2
 * \file io.cpp
3
 * \brief Set of utility io functions for KiT-RT
4
 * \author  J. Kusch, S. Schotthoefer, P. Stammer,  J. Wolters, T. Xiao
5
 */
6

7
#include "common/io.hpp"
8
#include "common/config.hpp"
9
#include "common/mesh.hpp"
10
#include "common/typedef.hpp"
11
#include "toolboxes/errormessages.hpp"
12
#include "toolboxes/textprocessingtoolbox.hpp"
13

14
#include <filesystem>
15
#include <fstream>
16
#include <iomanip>
17
#include <iostream>
18
#include <sstream>
19

20
#ifdef IMPORT_MPI
21
#include <mpi.h>
22
#endif
23

24
#include <omp.h>
25
#include <sstream>
26
#include <vector>
27
#include <vtkCellArray.h>
28
#include <vtkCellData.h>
29
#include <vtkCellDataToPointData.h>
30
#include <vtkDoubleArray.h>
31
#include <vtkPointData.h>
32
#include <vtkQuad.h>
33
#include <vtkSmartPointer.h>
34
#include <vtkTriangle.h>
35
#include <vtkUnstructuredGrid.h>
36
#include <vtkUnstructuredGridWriter.h>
37

38
// #include <Python.h>
39
// #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
40
// #include <numpy/arrayobject.h>
41

42
using vtkPointsSP                 = vtkSmartPointer<vtkPoints>;
43
using vtkUnstructuredGridSP       = vtkSmartPointer<vtkUnstructuredGrid>;
44
using vtkTriangleSP               = vtkSmartPointer<vtkTriangle>;
45
using vtkCellArraySP              = vtkSmartPointer<vtkCellArray>;
46
using vtkDoubleArraySP            = vtkSmartPointer<vtkDoubleArray>;
47
using vtkUnstructuredGridWriterSP = vtkSmartPointer<vtkUnstructuredGridWriter>;
48
using vtkCellDataToPointDataSP    = vtkSmartPointer<vtkCellDataToPointData>;
49

50
void ExportVTK( const std::string fileName,
23✔
51
                const std::vector<std::vector<std::vector<double>>>& outputFields,
52
                const std::vector<std::vector<std::string>>& outputFieldNames,
53
                const Mesh* mesh ) {
54
    int rank = 0;
23✔
55
#ifdef IMPORT_MPI
56
    // Initialize MPI
57
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
58
#endif
59
    if( rank == 0 ) {
23✔
60
        unsigned dim             = mesh->GetDim();
23✔
61
        unsigned numCells        = mesh->GetNumCells();
23✔
62
        unsigned numNodes        = mesh->GetNumNodes();
23✔
63
        auto nodes               = mesh->GetNodes();
46✔
64
        auto cells               = mesh->GetCells();
46✔
65
        unsigned numNodesPerCell = mesh->GetNumNodesPerCell();
23✔
66

67
        auto writer                 = vtkUnstructuredGridWriterSP::New();
46✔
68
        std::string fileNameWithExt = fileName;
46✔
69
        if( !TextProcessingToolbox::StringEndsWith( fileNameWithExt, ".vtk" ) ) {
23✔
70
            fileNameWithExt.append( ".vtk" );
23✔
71
        }
72
        writer->SetFileName( fileNameWithExt.c_str() );
23✔
73
        auto grid = vtkUnstructuredGridSP::New();
46✔
74
        auto pts  = vtkPointsSP::New();
46✔
75
        pts->SetDataTypeToDouble();
23✔
76
        pts->SetNumberOfPoints( static_cast<int>( numNodes ) );
23✔
77
        unsigned nodeID = 0;
23✔
78
        for( const auto& node : nodes ) {
8,170✔
79
            if( dim == 2 ) {
8,147✔
80
                pts->SetPoint( nodeID++, node[0], node[1], 0.0 );
8,147✔
81
            }
82
            else if( dim == 3 ) {
×
83
                pts->SetPoint( nodeID++, node[0], node[1], node[2] );
×
84
            }
85
            else {
86
                ErrorMessages::Error( "Unsupported dimension (d=" + std::to_string( dim ) + ")!", CURRENT_FUNCTION );
×
87
            }
88
        }
89
        vtkCellArraySP cellArray = vtkCellArraySP::New();
46✔
90
        if( numNodesPerCell == 3 ) {
23✔
91
            auto tri = vtkTriangleSP::New();
38✔
92
            for( unsigned i = 0; i < numCells; ++i ) {
7,811✔
93
                for( unsigned j = 0; j < numNodesPerCell; ++j ) {
31,168✔
94
                    tri->GetPointIds()->SetId( j, cells[i][j] );
23,376✔
95
                }
96
                cellArray->InsertNextCell( tri );
7,792✔
97
            }
98
        }
99

100
        if( numNodesPerCell == 4 ) {
23✔
101
            auto quad = vtkQuad::New();
4✔
102
            for( unsigned i = 0; i < numCells; ++i ) {
3,604✔
103
                for( unsigned j = 0; j < numNodesPerCell; ++j ) {
18,000✔
104
                    quad->GetPointIds()->SetId( j, cells[i][j] );
14,400✔
105
                }
106
                cellArray->InsertNextCell( quad );
3,600✔
107
            }
108
        }
109
        if( numNodesPerCell == 3 ) {
23✔
110
            grid->SetCells( VTK_TRIANGLE, cellArray );
19✔
111
        }
112
        else if( numNodesPerCell == 4 ) {
4✔
113
            grid->SetCells( VTK_QUAD, cellArray );
4✔
114
        }
115

116
        // Write the output
117
        auto cellData = vtkDoubleArraySP::New();
46✔
118

119
        for( unsigned idx_group = 0; idx_group < outputFields.size(); idx_group++ ) {
50✔
120

121
            for( unsigned idx_field = 0; idx_field < outputFields[idx_group].size(); idx_field++ ) {    // Loop over all output fields
58✔
122
                cellData = vtkDoubleArraySP::New();
31✔
123
                cellData->SetName( outputFieldNames[idx_group][idx_field].c_str() );
31✔
124

125
                for( unsigned idx_cell = 0; idx_cell < numCells; idx_cell++ ) {
18,623✔
126
                    cellData->InsertNextValue( outputFields[idx_group][idx_field][idx_cell] );
18,592✔
127
                }
128

129
                grid->GetCellData()->AddArray( cellData );
31✔
130
            }
131
        }
132

133
        grid->SetPoints( pts );
23✔
134
        grid->Squeeze();
23✔
135

136
        auto converter = vtkCellDataToPointDataSP::New();
46✔
137
        converter->AddInputDataObject( grid );
23✔
138
        converter->PassCellDataOn();
23✔
139
        converter->Update();
23✔
140

141
        writer->SetInputData( converter->GetOutput() );
23✔
142

143
        writer->Write();
23✔
144

145
        //  auto log = spdlog::get( "event" );
146
        //  log->info( "Result successfully exported to '{0}'!", fileNameWithExt );
147
    }
148
}
23✔
149

150
Mesh* LoadSU2MeshFromFile( const Config* settings ) {
35✔
151
    auto log = spdlog::get( "event" );
105✔
152
    // log->info( "| Importing mesh. This may take a while for large meshes." );
153
    unsigned dim = settings->GetDim();
35✔
154
    bool foundDim = false;
35✔
155
    std::vector<Vector> nodes;
70✔
156
    std::vector<std::vector<unsigned>> cells;
70✔
157
    std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries;
70✔
158

159
    if( !std::filesystem::exists( settings->GetMeshFile() ) )
35✔
160
        ErrorMessages::Error( "Cannot find mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
×
161
    std::ifstream ifs( settings->GetMeshFile(), std::ios::in );
70✔
162
    std::string line;
35✔
163

164
    if( ifs.is_open() ) {
35✔
165
        while( getline( ifs, line ) ) {
35✔
166
            if( line.find( "NDIME", 0 ) != std::string::npos ) {
35✔
167
                dim = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
35✔
168
                foundDim = true;
35✔
169
                // if( settings->GetDim() != dim ) {
170
                //     log->info( "Warning: Mesh dimension does not coinside with problem dimension! Proceed with caution!" );
171
                // }
172
                break;
35✔
173
            }
174
        }
175
        if( !foundDim ) {
35✔
NEW
176
            ErrorMessages::Error( "Invalid mesh file detected! Missing NDIME entry.", CURRENT_FUNCTION );
×
177
        }
178
        ifs.clear();
35✔
179
        ifs.seekg( 0, std::ios::beg );
35✔
180
        while( getline( ifs, line ) ) {
16,309✔
181
            if( line.find( "NPOIN", 0 ) != std::string::npos ) {
16,309✔
182
                unsigned numPoints = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
35✔
183
                nodes.resize( numPoints, Vector( dim, 0.0 ) );
35✔
184
                for( unsigned i = 0; i < numPoints; ++i ) {
9,844✔
185
                    getline( ifs, line );
9,809✔
186
                    std::stringstream ss;
19,618✔
187
                    ss << line;
9,809✔
188
                    unsigned id = 0;
9,809✔
189
                    Vector coords( dim, 0.0 );
9,809✔
190
                    while( !ss.eof() ) {
19,618✔
191
                        for( unsigned d = 0; d < dim; ++d ) {
29,427✔
192
                            ss >> coords[d];
19,618✔
193
                        }
194
                        ss >> id;
9,809✔
195
                    }
196

197
                    nodes[id] = Vector( coords );
9,809✔
198
                }
199
                break;
35✔
200
            }
201
        }
202
        ifs.clear();
35✔
203
        ifs.seekg( 0, std::ios::beg );
35✔
204
        while( getline( ifs, line ) ) {
26,153✔
205
            if( line.find( "NMARK", 0 ) != std::string::npos ) {
26,153✔
206
                unsigned numBCs = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
35✔
207
                boundaries.resize( numBCs );
35✔
208
                for( unsigned i = 0; i < numBCs; ++i ) {
70✔
209
                    std::string markerTag;
70✔
210
                    BOUNDARY_TYPE btype;
211
                    std::vector<unsigned> bnodes;
35✔
212
                    for( unsigned k = 0; k < 2; ++k ) {
105✔
213
                        getline( ifs, line );
70✔
214
                        if( line.find( "MARKER_TAG", 0 ) != std::string::npos ) {
70✔
215
                            markerTag    = line.substr( line.find( "=" ) + 1 );
35✔
216
                            auto end_pos = std::remove_if(
217
                                markerTag.begin(), markerTag.end(), []( char c ) { return std::isspace( static_cast<unsigned char>( c ) ); } );
210✔
218
                            markerTag.erase( end_pos, markerTag.end() );
35✔
219
                            btype = settings->GetBoundaryType( markerTag );
35✔
220
                            if( btype == BOUNDARY_TYPE::INVALID ) {
35✔
221
                                std::string errorMsg = std::string( "Invalid Boundary at marker \"" + markerTag + "\"." );
×
222
                                ErrorMessages::Error( errorMsg, CURRENT_FUNCTION );
×
223
                            }
224
                        }
225
                        else if( line.find( "MARKER_ELEMS", 0 ) != std::string::npos ) {
35✔
226
                            unsigned numMarkerElements = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
35✔
227
                            for( unsigned j = 0; j < numMarkerElements; ++j ) {
1,579✔
228
                                getline( ifs, line );
1,544✔
229
                                std::stringstream ss;
3,088✔
230
                                ss << line;
1,544✔
231
                                unsigned type = 0, id = 0;
1,544✔
232
                                while( !ss.eof() ) {
4,552✔
233
                                    ss >> type;
3,008✔
234
                                    for( unsigned d = 0; d < dim; ++d ) {
9,024✔
235
                                        ss >> id;
6,016✔
236
                                        bnodes.push_back( id );
6,016✔
237
                                    }
238
                                }
239
                            }
240
                        }
241
                        else {
242
                            ErrorMessages::Error( "Invalid mesh file detected! Make sure boundaries are provided.'", CURRENT_FUNCTION );
×
243
                        }
244
                    }
245
                    std::sort( bnodes.begin(), bnodes.end() );
35✔
246
                    auto last = std::unique( bnodes.begin(), bnodes.end() );
35✔
247
                    bnodes.erase( last, bnodes.end() );
35✔
248
                    boundaries[i] = std::make_pair( btype, bnodes );
35✔
249
                }
250
                break;
35✔
251
            }
252
        }
253
        ifs.clear();
35✔
254
        ifs.seekg( 0, std::ios::beg );
35✔
255
        std::vector<unsigned> numNodesPerCell;
70✔
256
        while( getline( ifs, line ) ) {
70✔
257
            if( line.find( "NELEM", 0 ) != std::string::npos ) {
70✔
258
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
35✔
259
                numNodesPerCell.resize( numCells, 0u );
35✔
260
                for( unsigned i = 0; i < numCells; ++i ) {
16,239✔
261
                    getline( ifs, line );
16,204✔
262
                    std::stringstream ss;
32,408✔
263
                    ss << line;
16,204✔
264
                    unsigned type = 0;
16,204✔
265
                    ss >> type;
16,204✔
266
                    switch( type ) {
16,204✔
267
                        case 5: {
14,404✔
268
                            numNodesPerCell[i] = 3;
14,404✔
269
                            break;
14,404✔
270
                        }
271
                        case 9: {
1,800✔
272
                            numNodesPerCell[i] = 4;
1,800✔
273
                            break;
1,800✔
274
                        }
275
                        default: {
×
276
                            ErrorMessages::Error( "Unsupported mesh type!'", CURRENT_FUNCTION );
×
277
                        }
278
                    }
279
                }
280
                break;
35✔
281
            }
282
        }
283
        bool mixedElementMesh = !std::equal( numNodesPerCell.begin() + 1, numNodesPerCell.end(), numNodesPerCell.begin() );
35✔
284
        if( mixedElementMesh ) {
35✔
285
            ErrorMessages::Error( "Mixed element meshes are currently not supported!'", CURRENT_FUNCTION );
×
286
        }
287
        ifs.clear();
35✔
288
        ifs.seekg( 0, std::ios::beg );
35✔
289
        while( getline( ifs, line ) ) {
70✔
290
            if( line.find( "NELEM", 0 ) != std::string::npos ) {
70✔
291
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
35✔
292
                cells.resize( numCells, std::vector<unsigned>( numNodesPerCell[0], 0u ) );
35✔
293
                for( unsigned i = 0; i < numCells; ++i ) {
16,239✔
294
                    getline( ifs, line );
16,204✔
295
                    std::stringstream ss;
32,408✔
296
                    ss << line;
16,204✔
297
                    unsigned type = 0, id = 0;
16,204✔
298
                    std::vector<unsigned> buffer( numNodesPerCell[0], 0u );
32,408✔
299
                    while( !ss.eof() ) {
32,408✔
300
                        ss >> type;
16,204✔
301
                        for( unsigned d = 0; d < numNodesPerCell[i]; ++d ) {
66,616✔
302
                            ss >> buffer[d];
50,412✔
303
                        }
304
                        ss >> id;
16,204✔
305
                    }
306
                    std::copy( buffer.begin(), buffer.end(), cells[id].begin() );
16,204✔
307
                }
308
                break;
35✔
309
            }
310
        }
311
    }
312
    else {
313
        ErrorMessages::Error( "Cannot open mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
×
314
    }
315
    ifs.close();
35✔
316
    // log->info( "| Mesh imported." );
317
    return new Mesh( settings, nodes, cells, boundaries );
70✔
318
}
319

320
void LoadConnectivityFromFile( const std::string inputFile,
26✔
321
                               std::vector<std::vector<unsigned>>& cellNeighbors,
322
                               std::vector<std::vector<Vector>>& cellInterfaceMidPoints,
323
                               std::vector<std::vector<Vector>>& cellNormals,
324
                               std::vector<BOUNDARY_TYPE>& cellBoundaryTypes,
325
                               unsigned nCells,
326
                               unsigned nNodesPerCell,
327
                               unsigned nDim ) {
328
    // File has nCells lines, each line is a comma separated entry containing:
329
    // cellNeighbors (nNodesPerCell elements),
330
    // cellInterfaceMidPoints (nNodesPerCell x nDim elements),
331
    // cellNormals (nNodesPerCell x nDim elements),
332
    // cellBoundaryTypes (1 element), (tranlated from  unsigned to enum BOUNDARY_TYPE)
333

334
    std::ifstream inFile( inputFile );
26✔
335

336
    if( !inFile.is_open() ) {
26✔
337
        ErrorMessages::Error( "Error opening connectivity file.", CURRENT_FUNCTION );
×
338
        return;
×
339
    }
340
    for( unsigned i = 0; i < nCells; ++i ) {
12,710✔
341
        std::string line;
25,368✔
342
        unsigned count = 1;
12,684✔
343

344
        std::getline( inFile, line );
12,684✔
345
        for( char ch : line ) {
2,953,922✔
346
            if( ch == ',' ) {
2,941,238✔
347
                count++;
190,260✔
348
            }
349
        }
350

351
        unsigned correctedNodesPerCell = nNodesPerCell;
12,684✔
352
        if( count < nNodesPerCell + nNodesPerCell * nDim * 2 + 1 ) {
12,684✔
353
            correctedNodesPerCell = nNodesPerCell - 1;
×
354
            ErrorMessages::Error( "Error opening connectivity file.", CURRENT_FUNCTION );
×
355
        }
356

357
        std::istringstream iss( line );
25,368✔
358
        // Load cellNeighbors
359
        cellNeighbors[i].resize( correctedNodesPerCell );
12,684✔
360

361
        for( unsigned j = 0; j < correctedNodesPerCell; ++j ) {
50,736✔
362
            std::getline( iss, line, ',' );
38,052✔
363
            std::istringstream converter( line );
76,104✔
364
            converter >> std::fixed >> setprecision( 15 ) >> cellNeighbors[i][j];
38,052✔
365
        }
366

367
        // Load cellInterfaceMidPoints
368
        cellInterfaceMidPoints[i].resize( correctedNodesPerCell );
12,684✔
369
        for( unsigned j = 0; j < correctedNodesPerCell; ++j ) {
50,736✔
370
            cellInterfaceMidPoints[i][j] = Vector( nDim, 0.0 );
38,052✔
371
            for( unsigned k = 0; k < nDim; ++k ) {
114,156✔
372
                std::getline( iss, line, ',' );
76,104✔
373
                std::istringstream converter( line );
152,208✔
374
                converter >> std::fixed >> setprecision( 15 ) >> cellInterfaceMidPoints[i][j][k];    // Replace with appropriate member of Vector
76,104✔
375
                // std::cout << std::fixed << setprecision( 15 ) << cellInterfaceMidPoints[i][j][k] << std::endl;
376
            }
377
        }
378
        // Load cellNormals
379
        cellNormals[i].resize( correctedNodesPerCell );
12,684✔
380
        for( unsigned j = 0; j < correctedNodesPerCell; ++j ) {
50,736✔
381
            cellNormals[i][j] = Vector( nDim, 0.0 );
38,052✔
382
            for( unsigned k = 0; k < nDim; ++k ) {
114,156✔
383
                std::getline( iss, line, ',' );
76,104✔
384
                std::istringstream converter( line );
152,208✔
385
                converter >> std::fixed >> setprecision( 15 ) >> cellNormals[i][j][k];    // Replace with appropriate member of Vector
76,104✔
386
            }
387
        }
388
        // Load cellBoundaryTypes
389
        std::getline( iss, line, ',' );
12,684✔
390
        std::istringstream converter( line );
12,684✔
391
        unsigned boundaryTypeValue;
392
        converter >> boundaryTypeValue;
12,684✔
393
        cellBoundaryTypes[i] = static_cast<BOUNDARY_TYPE>( boundaryTypeValue );
12,684✔
394
    }
395
    inFile.close();
26✔
396
}
397

398
void WriteConnecitivityToFile( const std::string outputFile,
7✔
399
                               const std::vector<std::vector<unsigned>>& cellNeighbors,
400
                               const std::vector<std::vector<Vector>>& cellInterfaceMidPoints,
401
                               const std::vector<std::vector<Vector>>& cellNormals,
402
                               const std::vector<BOUNDARY_TYPE>& cellBoundaryTypes,
403
                               unsigned nCells,
404
                               unsigned nDim ) {
405
    // File has nCells lines, each line is a comma separated entry containing:
406
    // cellNeighbors (nNodesPerCell elements),
407
    // cellInterfaceMidPoints (nNodesPerCell x nDim elements),
408
    // cellNormals (nNodesPerCell x nDim elements),
409
    // cellBoundaryTypes (1 element), (tranlated from BOUNDARY_TYPE to unsigned)
410

411
    std::ofstream outFile( outputFile );
14✔
412
    outFile << std::fixed << setprecision( 15 );
7✔
413
    // const std::size_t bufferSize = 10000;    // Adjust as needed
414
    // outFile.rdbuf()->pubsetbuf( 0, bufferSize );
415
    if( outFile.is_open() ) {
7✔
416
        // Write data to the file
417
        for( unsigned i = 0; i < nCells; ++i ) {
1,185✔
418
            // Write cellNeighbors
419
            for( unsigned j = 0; j < cellNeighbors[i].size(); ++j ) {
4,712✔
420
                outFile << cellNeighbors[i][j];
3,534✔
421
                if( j < cellNeighbors[i].size() - 1 ) {
3,534✔
422
                    outFile << ",";
2,356✔
423
                }
424
            }
425
            // Write cellInterfaceMidPoints
426
            for( unsigned j = 0; j < cellInterfaceMidPoints[i].size(); ++j ) {
4,712✔
427
                for( unsigned k = 0; k < nDim; ++k ) {
10,602✔
428
                    outFile << "," << cellInterfaceMidPoints[i][j][k];
7,068✔
429
                }
430
            }
431
            // Write cellNormals
432
            for( unsigned j = 0; j < cellNormals[i].size(); ++j ) {
4,712✔
433
                for( unsigned k = 0; k < nDim; ++k ) {
10,602✔
434
                    outFile << "," << cellNormals[i][j][k];
7,068✔
435
                }
436
            }
437

438
            // Write cellBoundaryTypes
439
            outFile << "," << static_cast<unsigned>( cellBoundaryTypes[i] );
1,178✔
440
            outFile << std::endl;
1,178✔
441
        }
442
        outFile.close();
7✔
443
    }
444
    else {
445
        ErrorMessages::Error( "Error opening connectivity file.", CURRENT_FUNCTION );
×
446
    }
447
}
7✔
448

449
void WriteRestartSolution( const std::string& baseOutputFile,
4✔
450
                           const std::vector<double>& solution,
451
                           const std::vector<double>& scalarFlux,
452
                           int rank,
453
                           int idx_iter,
454
                           double totalAbsorptionHohlraumCenter,
455
                           double totalAbsorptionHohlraumVertical,
456
                           double totalAbsorptionHohlraumHorizontal,
457
                           double totalAbsorption ) {
458
    // Generate filename with rank number
459
    std::string outputFile = baseOutputFile + "_restart_rank_" + std::to_string( rank );
8✔
460

461
    // Check if the file exists and delete it
462
    if( std::filesystem::exists( outputFile ) ) {
4✔
463
        std::filesystem::remove( outputFile );
×
464
    }
465

466
    // Open the file for binary writing
467
    std::ofstream outFile( outputFile, std::ios::out | std::ios::binary );
4✔
468
    if( !outFile ) {
4✔
469
        std::cerr << "Error opening restart solution file." << std::endl;
×
470
        return;
×
471
    }
472

473
    // Write the iteration number as the first item (in binary)
474
    outFile.write( reinterpret_cast<const char*>( &idx_iter ), sizeof( idx_iter ) );
4✔
475
    outFile.write( reinterpret_cast<const char*>( &totalAbsorptionHohlraumCenter ), sizeof( totalAbsorptionHohlraumCenter ) );
4✔
476
    outFile.write( reinterpret_cast<const char*>( &totalAbsorptionHohlraumVertical ), sizeof( totalAbsorptionHohlraumVertical ) );
4✔
477
    outFile.write( reinterpret_cast<const char*>( &totalAbsorptionHohlraumHorizontal ), sizeof( totalAbsorptionHohlraumHorizontal ) );
4✔
478
    outFile.write( reinterpret_cast<const char*>( &totalAbsorption ), sizeof( totalAbsorption ) );
4✔
479

480
    // Write the size of the solution and scalarFlux vectors (optional but useful for reading)
481
    size_t solution_size   = solution.size();
4✔
482
    size_t scalarFlux_size = scalarFlux.size();
4✔
483
    outFile.write( reinterpret_cast<const char*>( &solution_size ), sizeof( solution_size ) );
4✔
484
    outFile.write( reinterpret_cast<const char*>( &scalarFlux_size ), sizeof( scalarFlux_size ) );
4✔
485

486
    // Write each element of the solution vector in binary
487
    outFile.write( reinterpret_cast<const char*>( solution.data() ), solution_size * sizeof( double ) );
4✔
488

489
    // Write each element of the scalarFlux vector in binary
490
    outFile.write( reinterpret_cast<const char*>( scalarFlux.data() ), scalarFlux_size * sizeof( double ) );
4✔
491

492
    // Close the file
493
    outFile.close();
4✔
494
}
495

496
int LoadRestartSolution( const std::string& baseInputFile,
4✔
497
                         std::vector<double>& solution,
498
                         std::vector<double>& scalarFlux,
499
                         int rank,
500
                         unsigned long nCells,
501
                         double& totalAbsorptionHohlraumCenter,
502
                         double& totalAbsorptionHohlraumVertical,
503
                         double& totalAbsorptionHohlraumHorizontal,
504
                         double& totalAbsorption ) {
505
    // Generate filename with rank number
506
    std::string inputFile = baseInputFile + "_restart_rank_" + std::to_string( rank );
12✔
507

508
    // Open the file for binary reading
509
    std::ifstream inFile( inputFile, std::ios::in | std::ios::binary );
8✔
510
    if( !inFile ) {
4✔
511
        std::cerr << "Error opening restart solution file for reading." << std::endl;
×
512
        return 0;    // Signal failure
×
513
    }
514

515
    // Read the iteration number
516
    int idx_iter = 0;
4✔
517
    inFile.read( reinterpret_cast<char*>( &idx_iter ), sizeof( idx_iter ) );
4✔
518
    inFile.read( reinterpret_cast<char*>( &totalAbsorptionHohlraumCenter ), sizeof( totalAbsorptionHohlraumCenter ) );
4✔
519
    inFile.read( reinterpret_cast<char*>( &totalAbsorptionHohlraumVertical ), sizeof( totalAbsorptionHohlraumVertical ) );
4✔
520
    inFile.read( reinterpret_cast<char*>( &totalAbsorptionHohlraumHorizontal ), sizeof( totalAbsorptionHohlraumHorizontal ) );
4✔
521
    inFile.read( reinterpret_cast<char*>( &totalAbsorption ), sizeof( totalAbsorption ) );
4✔
522

523
    // Read the size of the solution and scalarFlux vectors
524
    size_t solution_size, scalarFlux_size;
525
    inFile.read( reinterpret_cast<char*>( &solution_size ), sizeof( solution_size ) );
4✔
526
    inFile.read( reinterpret_cast<char*>( &scalarFlux_size ), sizeof( scalarFlux_size ) );
4✔
527

528
    // Temporary vector to hold the full data
529
    std::vector<double> tempData( solution_size + scalarFlux_size );
8✔
530

531
    // Read the entire data block in binary form
532
    inFile.read( reinterpret_cast<char*>( tempData.data() ), ( solution_size + scalarFlux_size ) * sizeof( double ) );
4✔
533

534
    // Close the file
535
    inFile.close();
4✔
536

537
    // Verify that we have enough entries to populate scalarFlux
538
    if( tempData.size() < nCells ) {
4✔
539
        std::cerr << "Not enough data to populate scalar flux vector." << std::endl;
×
540
        return 0;    // Signal failure
×
541
    }
542

543
    // Allocate the last nCells entries to scalarFlux and the rest to solution
544
    solution.assign( tempData.begin(), tempData.end() - nCells );
4✔
545
    scalarFlux.assign( tempData.end() - nCells, tempData.end() );
4✔
546

547
    return idx_iter;    // Return the iteration index
4✔
548
}
549

550
std::string ParseArguments( int argc, char* argv[] ) {
×
551
    std::string inputFile;
×
552
    std::string usage_help = "\nUsage: " + std::string( argv[0] ) + " inputfile\n";
×
553

554
    if( argc < 2 ) {
×
555
        std::cout << usage_help;
×
556
        exit( EXIT_FAILURE );
×
557
    }
558
    for( int i = 1; i < argc; i++ ) {
×
559
        std::string arg = argv[i];
×
560
        if( arg == "-h" ) {
×
561
            std::cout << usage_help;
×
562
            exit( EXIT_SUCCESS );
×
563
        }
564
        else {
565
            inputFile = std::string( argv[i] );
×
566
            std::ifstream f( inputFile );
×
567
            if( !f.is_open() ) {
×
568
                ErrorMessages::Error( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
×
569
            }
570
        }
571
    }
572
    return inputFile;
×
573
}
574

575
void PrintLogHeader( std::string inputFile ) {
×
576
    int nprocs = 1;
×
577
    int rank   = 0;
×
578
#ifdef IMPORT_MPI
579
    // Initialize MPI
580
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
581
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
582
#endif
583
    if( rank == 0 ) {
×
584
        auto log = spdlog::get( "event" );
×
585

586
        // New design
587
        log->info( "------------------------------------------------------------------------" );
×
588
        log->info( "|    _  _____ _____     ____ _____                                     |" );
×
589
        log->info( "|   | |/ /_ _|_   _|   |  _ \\_   _|                                    |" );
×
590
        log->info( "|   | ' / | |  | |_____| |_) || |            Version                   |" );
×
591
        log->info( "|   | . \\ | |  | |_____|  _ < | |             1.0.0                    |" );
×
592
        log->info( "|   |_|\\_\\___| |_|     |_| \\_\\|_|                                      |" );
×
593
        log->info( "|                                                                      |" );
×
594
        log->info( "------------------------------------------------------------------------" );
×
595
        log->info( "|  Copyright: MIT Licence                                              |" );
×
596
        // log->info( "|  Authors: S. Schotthöfer, J. Kusch, J. Wolters, P. Stammer ad T. Xiao |" );
597
        log->info( "|  Date: 21th April 2022                                               |" );
×
598
        log->info( "------------------------------------------------------------------------" );
×
599
        log->info( "|" );
×
600
        log->info( "| Git commit :\t{0}", GIT_HASH );
×
601
        log->info( "| Config file:\t{0}", inputFile );
×
602
        log->info( "| MPI Threads:\t{0}", nprocs );
×
603
        log->info( "| OMP Threads:\t{0}", omp_get_max_threads() );
×
604
        log->info( "|" );
×
605
        log->info( "-------------------------- Config File Info ----------------------------" );
×
606
        log->info( "|" );
×
607
        // print file content while omitting comments
608
        std::ifstream ifs( inputFile );
×
609
        if( ifs.is_open() ) {
×
610
            std::string line;
×
611
            while( !ifs.eof() ) {
×
612
                std::getline( ifs, line );
×
613
                if( line[0] != '%' ) log->info( "| {0}", line );
×
614
            }
615
        }
616
        // log->info( "------------------------------------------------------------------------" );
617
    }
618
#ifdef IMPORT_MPI
619
    MPI_Barrier( MPI_COMM_WORLD );
620
#endif
621
}
622

623
/*
624
Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ) {
625
    auto log = spdlog::get( "event" );
626

627
    // std::cout << "imageName" << imageName << "\n";
628
    // std::cout << "pyhtonPath" << KITRT_PYTHON_PATH;
629

630
    if( !std::filesystem::exists( imageName ) ) {
631
        ErrorMessages::Error( "Can not open image '" + imageName + "'!", CURRENT_FUNCTION );
632
    }
633
    std::filesystem::path outDir( std::filesystem::path( SU2Filename ).parent_path() );
634
    if( !std::filesystem::exists( outDir ) ) {
635
        ErrorMessages::Error( "Output directory '" + outDir.string() + "' does not exists!", CURRENT_FUNCTION );
636
    }
637

638
    std::string pyPath = KITRT_PYTHON_PATH;
639

640
    if( !Py_IsInitialized() ) {
641
        Py_InitializeEx( 0 );
642
        if( !Py_IsInitialized() ) {
643
            ErrorMessages::Error( "Python init failed!", CURRENT_FUNCTION );
644
        }
645
        PyRun_SimpleString( ( "import sys\nsys.path.append('" + pyPath + "')" ).c_str() );
646
    }
647

648
    PyObject *pArgs, *pReturn, *pModule, *pFunc;
649
    PyArrayObject* np_ret;
650

651
    auto image = PyUnicode_FromString( imageName.c_str() );
652
    auto su2   = PyUnicode_FromString( SU2Filename.c_str() );
653

654
    std::string moduleName = "mesh_from_image";
655
    pModule                = PyImport_ImportModule( moduleName.c_str() );
656
    if( !pModule ) {
657
        PyErr_Print();
658
        ErrorMessages::Error( "'" + moduleName + "' can not be imported!", CURRENT_FUNCTION );
659
    }
660

661
    pFunc = PyObject_GetAttrString( pModule, "generate" );
662
    if( !pFunc || !PyCallable_Check( pFunc ) ) {
663
        PyErr_Print();
664
        Py_DecRef( pModule );
665
        Py_DecRef( pFunc );
666
        ErrorMessages::Error( "'generate' is null or not callable!", CURRENT_FUNCTION );
667
    }
668

669
    pArgs = PyTuple_New( 2 );
670
    PyTuple_SetItem( pArgs, 0, reinterpret_cast<PyObject*>( image ) );
671
    PyTuple_SetItem( pArgs, 1, reinterpret_cast<PyObject*>( su2 ) );
672
    pReturn = PyObject_CallObject( pFunc, pArgs );
673
    np_ret  = reinterpret_cast<PyArrayObject*>( pReturn );
674

675
    size_t m{ static_cast<size_t>( PyArray_SHAPE( np_ret )[0] ) };
676
    size_t n{ static_cast<size_t>( PyArray_SHAPE( np_ret )[1] ) };
677
    double* c_out = reinterpret_cast<double*>( PyArray_DATA( np_ret ) );
678

679
    Matrix gsImage( m, n, c_out );
680

681
    // Finalizing
682
    Py_DecRef( pFunc );
683
    Py_DecRef( pModule );
684
    Py_DECREF( np_ret );
685

686
    return gsImage.transpose();
687
}
688
*/
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc