CSMMLab / KiT-RT / #95

Committed 28 May 2025 02:06AM UTC coverage: 46.655% (-11.1%) from 57.728%

Build # #95

Build Type

push

travis-ci

Committed by

web-flow

Commit Message

Release 2025 (#44)

* New neural models (#30)

* extend kitrt script

* added new neural models

* extend to m3 models

* added M3_2D models

* added M2 and M4 models

* configure ml optimizer for high order models

* added configuration for high order models

* add option for rotation postprcessing in neural networks. remove unneccessary python scripts

* started rotational symmetric postprocessing

* added rotational invariance postprocessing for m2 and m1 models

* fix post merge bugs

* created hohlraum mesh

* add hohlraum tes case

* add hohlraum test case

* add hohlraum cfg filees

* fixed hohlraum testcase

* add hohlraum cfg files

* changed hohlraum cfg

* changed hohlraum testcase to isotropic inflow source boundary condition

* added ghost cell bonudary for hohlraum testcase

* update readme and linesource mesh creator

* added proper scaling for linesource reference solution

* regularized newton debugging

* Data generator with reduced objective functional (#33)

* added reduced optimizer for sampling

* remove old debugging comments

* mesh acceleration (#38)

* added new ansatz (#36)

* added new ansatz

* debug new ansatz

* branch should be abandoned here

* debug new ansatz

* fix scaling error new ansatz

* fix config errors

* temporarily fixed dynamic ansatz rotation bug

* fix inheritance error for new ansatz

* mesh acceleration

* add structured hohlraum

* Mesh acc (#41)

* mesh acceleration

* added floor value for starmap moment methods

* enable accelleration

* delete minimum value starmap

* Update README.md

* Spherical harmonics nn (#40)

* added new ansatz

* debug new ansatz

* branch should be abandoned here

* debug new ansatz

* fix scaling error new ansatz

* fix config errors

* temporarily fixed dynamic ansatz rotation bug

* fix inheritance error for new ansatz

* mesh acceleration

* add structured hohlraum

* added sph... (continued)

Run Details

767 of 3019 new or added lines in 51 files covered. (25.41%)

131 existing lines in 8 files now uncovered.

4422 of 9478 relevant lines covered (46.66%)

57163.05 hits per line

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

67.89

/src/common/io.cpp

/*!
 * \file io.cpp
 * \brief Set of utility io functions for KiT-RT
 * \author  J. Kusch, S. Schotthoefer, P. Stammer,  J. Wolters, T. Xiao
 */

#include "common/io.hpp"
#include "common/config.hpp"
#include "common/mesh.hpp"
#include "common/typedef.hpp"
#include "toolboxes/errormessages.hpp"
#include "toolboxes/textprocessingtoolbox.hpp"

#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <sstream>

#ifdef IMPORT_MPI
#include <mpi.h>
#endif

#include <omp.h>
#include <sstream>
#include <vector>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkQuad.h>
#include <vtkSmartPointer.h>
#include <vtkTriangle.h>
#include <vtkUnstructuredGrid.h>
#include <vtkUnstructuredGridWriter.h>

// #include <Python.h>
// #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
// #include <numpy/arrayobject.h>

using vtkPointsSP                 = vtkSmartPointer<vtkPoints>;
using vtkUnstructuredGridSP       = vtkSmartPointer<vtkUnstructuredGrid>;
using vtkTriangleSP               = vtkSmartPointer<vtkTriangle>;
using vtkCellArraySP              = vtkSmartPointer<vtkCellArray>;
using vtkDoubleArraySP            = vtkSmartPointer<vtkDoubleArray>;
using vtkUnstructuredGridWriterSP = vtkSmartPointer<vtkUnstructuredGridWriter>;
using vtkCellDataToPointDataSP    = vtkSmartPointer<vtkCellDataToPointData>;

void ExportVTK( const std::string fileName,
                const std::vector<std::vector<std::vector<double>>>& outputFields,
                const std::vector<std::vector<std::string>>& outputFieldNames,
                const Mesh* mesh ) {
    int rank   = 0;
    int nprocs = 1;
#ifdef IMPORT_MPI
    // Initialize MPI
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
#endif
    if( rank == 0 ) {
        unsigned dim             = mesh->GetDim();
        unsigned numCells        = mesh->GetNumCells();
        unsigned numNodes        = mesh->GetNumNodes();
        auto nodes               = mesh->GetNodes();
        auto cells               = mesh->GetCells();
        unsigned numNodesPerCell = mesh->GetNumNodesPerCell();

        auto writer                 = vtkUnstructuredGridWriterSP::New();
        std::string fileNameWithExt = fileName;
        if( !TextProcessingToolbox::StringEndsWith( fileNameWithExt, ".vtk" ) ) {
            fileNameWithExt.append( ".vtk" );
        }
        writer->SetFileName( fileNameWithExt.c_str() );
        auto grid = vtkUnstructuredGridSP::New();
        auto pts  = vtkPointsSP::New();
        pts->SetDataTypeToDouble();
        pts->SetNumberOfPoints( static_cast<int>( numNodes ) );
        unsigned nodeID = 0;
        for( const auto& node : nodes ) {
            if( dim == 2 ) {
                pts->SetPoint( nodeID++, node[0], node[1], 0.0 );
            }
            else if( dim == 3 ) {
                pts->SetPoint( nodeID++, node[0], node[1], node[2] );
            }
            else {
                ErrorMessages::Error( "Unsupported dimension (d=" + std::to_string( dim ) + ")!", CURRENT_FUNCTION );
            }
        }
        vtkCellArraySP cellArray = vtkCellArraySP::New();
        if( numNodesPerCell == 3 ) {
            auto tri = vtkTriangleSP::New();
            for( unsigned i = 0; i < numCells; ++i ) {
                for( unsigned j = 0; j < numNodesPerCell; ++j ) {
                    tri->GetPointIds()->SetId( j, cells[i][j] );
                }
                cellArray->InsertNextCell( tri );
            }
        }

        if( numNodesPerCell == 4 ) {
            auto quad = vtkQuad::New();
            for( unsigned i = 0; i < numCells; ++i ) {
                for( unsigned j = 0; j < numNodesPerCell; ++j ) {
                    quad->GetPointIds()->SetId( j, cells[i][j] );
                }
                cellArray->InsertNextCell( quad );
            }
        }
        if( numNodesPerCell == 3 ) {
            grid->SetCells( VTK_TRIANGLE, cellArray );
        }
        else if( numNodesPerCell == 4 ) {
            grid->SetCells( VTK_QUAD, cellArray );
        }

        // Write the output
        auto cellData = vtkDoubleArraySP::New();

        for( unsigned idx_group = 0; idx_group < outputFields.size(); idx_group++ ) {

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

                for( unsigned idx_cell = 0; idx_cell < numCells; idx_cell++ ) {
                    cellData->InsertNextValue( outputFields[idx_group][idx_field][idx_cell] );
                }

                grid->GetCellData()->AddArray( cellData );
            }
        }

        grid->SetPoints( pts );
        grid->Squeeze();

        auto converter = vtkCellDataToPointDataSP::New();
        converter->AddInputDataObject( grid );
        converter->PassCellDataOn();
        converter->Update();

        writer->SetInputData( converter->GetOutput() );

        writer->Write();

        //  auto log = spdlog::get( "event" );
        //  log->info( "Result successfully exported to '{0}'!", fileNameWithExt );
    }
}

Mesh* LoadSU2MeshFromFile( const Config* settings ) {
    auto log = spdlog::get( "event" );
    // log->info( "| Importing mesh. This may take a while for large meshes." );
    unsigned dim;
    std::vector<Vector> nodes;
    std::vector<std::vector<unsigned>> cells;
    std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries;

    if( !std::filesystem::exists( settings->GetMeshFile() ) )
        ErrorMessages::Error( "Cannot find mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
    std::ifstream ifs( settings->GetMeshFile(), std::ios::in );
    std::string line;

    if( ifs.is_open() ) {
        while( getline( ifs, line ) ) {
            if( line.find( "NDIME", 0 ) != std::string::npos ) {
                dim = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
                // if( settings->GetDim() != dim ) {
                //     log->info( "Warning: Mesh dimension does not coinside with problem dimension! Proceed with caution!" );
                // }
                break;
            }
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        while( getline( ifs, line ) ) {
            if( line.find( "NPOIN", 0 ) != std::string::npos ) {
                unsigned numPoints = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
                nodes.resize( numPoints, Vector( dim, 0.0 ) );
                for( unsigned i = 0; i < numPoints; ++i ) {
                    getline( ifs, line );
                    std::stringstream ss;
                    ss << line;
                    unsigned id = 0;
                    Vector coords( dim, 0.0 );
                    while( !ss.eof() ) {
                        for( unsigned d = 0; d < dim; ++d ) {
                            ss >> coords[d];
                        }
                        ss >> id;
                    }

                    nodes[id] = Vector( coords );
                }
                break;
            }
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        while( getline( ifs, line ) ) {
            if( line.find( "NMARK", 0 ) != std::string::npos ) {
                unsigned numBCs = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
                boundaries.resize( numBCs );
                for( unsigned i = 0; i < numBCs; ++i ) {
                    std::string markerTag;
                    BOUNDARY_TYPE btype;
                    std::vector<unsigned> bnodes;
                    for( unsigned k = 0; k < 2; ++k ) {
                        getline( ifs, line );
                        if( line.find( "MARKER_TAG", 0 ) != std::string::npos ) {
                            markerTag    = line.substr( line.find( "=" ) + 1 );
                            auto end_pos = std::remove_if(
                                markerTag.begin(), markerTag.end(), []( char c ) { return std::isspace( static_cast<unsigned char>( c ) ); } );
                            markerTag.erase( end_pos, markerTag.end() );
                            btype = settings->GetBoundaryType( markerTag );
                            if( btype == BOUNDARY_TYPE::INVALID ) {
                                std::string errorMsg = std::string( "Invalid Boundary at marker \"" + markerTag + "\"." );
                                ErrorMessages::Error( errorMsg, CURRENT_FUNCTION );
                            }
                        }
                        else if( line.find( "MARKER_ELEMS", 0 ) != std::string::npos ) {
                            unsigned numMarkerElements = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
                            for( unsigned j = 0; j < numMarkerElements; ++j ) {
                                getline( ifs, line );
                                std::stringstream ss;
                                ss << line;
                                unsigned type = 0, id = 0;
                                while( !ss.eof() ) {
                                    ss >> type;
                                    for( unsigned d = 0; d < dim; ++d ) {
                                        ss >> id;
                                        bnodes.push_back( id );
                                    }
                                }
                            }
                        }
                        else {
                            ErrorMessages::Error( "Invalid mesh file detected! Make sure boundaries are provided.'", CURRENT_FUNCTION );
                        }
                    }
                    std::sort( bnodes.begin(), bnodes.end() );
                    auto last = std::unique( bnodes.begin(), bnodes.end() );
                    bnodes.erase( last, bnodes.end() );
                    boundaries[i] = std::make_pair( btype, bnodes );
                }
                break;
            }
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        std::vector<unsigned> numNodesPerCell;
        while( getline( ifs, line ) ) {
            if( line.find( "NELEM", 0 ) != std::string::npos ) {
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
                numNodesPerCell.resize( numCells, 0u );
                for( unsigned i = 0; i < numCells; ++i ) {
                    getline( ifs, line );
                    std::stringstream ss;
                    ss << line;
                    unsigned type = 0;
                    ss >> type;
                    switch( type ) {
                        case 5: {
                            numNodesPerCell[i] = 3;
                            break;
                        }
                        case 9: {
                            numNodesPerCell[i] = 4;
                            break;
                        }
                        default: {
                            ErrorMessages::Error( "Unsupported mesh type!'", CURRENT_FUNCTION );
                        }
                    }
                }
                break;
            }
        }
        bool mixedElementMesh = !std::equal( numNodesPerCell.begin() + 1, numNodesPerCell.end(), numNodesPerCell.begin() );
        if( mixedElementMesh ) {
            ErrorMessages::Error( "Mixed element meshes are currently not supported!'", CURRENT_FUNCTION );
        }
        ifs.clear();
        ifs.seekg( 0, std::ios::beg );
        while( getline( ifs, line ) ) {
            if( line.find( "NELEM", 0 ) != std::string::npos ) {
                unsigned numCells = static_cast<unsigned>( TextProcessingToolbox::GetTrailingNumber( line ) );
                cells.resize( numCells, std::vector<unsigned>( numNodesPerCell[0], 0u ) );
                for( unsigned i = 0; i < numCells; ++i ) {
                    getline( ifs, line );
                    std::stringstream ss;
                    ss << line;
                    unsigned type = 0, id = 0;
                    std::vector<unsigned> buffer( numNodesPerCell[0], 0u );
                    while( !ss.eof() ) {
                        ss >> type;
                        for( unsigned d = 0; d < numNodesPerCell[i]; ++d ) {
                            ss >> buffer[d];
                        }
                        ss >> id;
                    }
                    std::copy( buffer.begin(), buffer.end(), cells[id].begin() );
                }
                break;
            }
        }
    }
    else {
        ErrorMessages::Error( "Cannot open mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
    }
    ifs.close();
    // log->info( "| Mesh imported." );
    return new Mesh( settings, nodes, cells, boundaries );
}

void LoadConnectivityFromFile( const std::string inputFile,
                               std::vector<std::vector<unsigned>>& cellNeighbors,
                               std::vector<std::vector<Vector>>& cellInterfaceMidPoints,
                               std::vector<std::vector<Vector>>& cellNormals,
                               std::vector<BOUNDARY_TYPE>& cellBoundaryTypes,
                               unsigned nCells,
                               unsigned nNodesPerCell,
                               unsigned nDim ) {
    // File has nCells lines, each line is a comma separated entry containing:
    // cellNeighbors (nNodesPerCell elements),
    // cellInterfaceMidPoints (nNodesPerCell x nDim elements),
    // cellNormals (nNodesPerCell x nDim elements),
    // cellBoundaryTypes (1 element), (tranlated from  unsigned to enum BOUNDARY_TYPE)

    std::ifstream inFile( inputFile );

    if( !inFile.is_open() ) {
        ErrorMessages::Error( "Error opening connectivity file.", CURRENT_FUNCTION );
        return;
    }
    for( unsigned i = 0; i < nCells; ++i ) {
        std::string line;
        unsigned count = 1;

        std::getline( inFile, line );
        for( char ch : line ) {
            if( ch == ',' ) {
                count++;
            }
        }

        unsigned correctedNodesPerCell = nNodesPerCell;
        if( count < nNodesPerCell + nNodesPerCell * nDim * 2 + 1 ) {
            correctedNodesPerCell = nNodesPerCell - 1;
            ErrorMessages::Error( "Error opening connectivity file.", CURRENT_FUNCTION );
        }

        std::istringstream iss( line );
        // Load cellNeighbors
        cellNeighbors[i].resize( correctedNodesPerCell );

        for( unsigned j = 0; j < correctedNodesPerCell; ++j ) {
            std::getline( iss, line, ',' );
            std::istringstream converter( line );
            converter >> std::fixed >> setprecision( 15 ) >> cellNeighbors[i][j];
        }

        // Load cellInterfaceMidPoints
        cellInterfaceMidPoints[i].resize( correctedNodesPerCell );
        for( unsigned j = 0; j < correctedNodesPerCell; ++j ) {
            cellInterfaceMidPoints[i][j] = Vector( nDim, 0.0 );
            for( unsigned k = 0; k < nDim; ++k ) {
                std::getline( iss, line, ',' );
                std::istringstream converter( line );
                converter >> std::fixed >> setprecision( 15 ) >> cellInterfaceMidPoints[i][j][k];    // Replace with appropriate member of Vector
                // std::cout << std::fixed << setprecision( 15 ) << cellInterfaceMidPoints[i][j][k] << std::endl;
            }
        }
        // Load cellNormals
        cellNormals[i].resize( correctedNodesPerCell );
        for( unsigned j = 0; j < correctedNodesPerCell; ++j ) {
            cellNormals[i][j] = Vector( nDim, 0.0 );
            for( unsigned k = 0; k < nDim; ++k ) {
                std::getline( iss, line, ',' );
                std::istringstream converter( line );
                converter >> std::fixed >> setprecision( 15 ) >> cellNormals[i][j][k];    // Replace with appropriate member of Vector
            }
        }
        // Load cellBoundaryTypes
        std::getline( iss, line, ',' );
        std::istringstream converter( line );
        unsigned boundaryTypeValue;
        converter >> boundaryTypeValue;
        cellBoundaryTypes[i] = static_cast<BOUNDARY_TYPE>( boundaryTypeValue );
    }
    inFile.close();
}

void WriteConnecitivityToFile( const std::string outputFile,
                               const std::vector<std::vector<unsigned>>& cellNeighbors,
                               const std::vector<std::vector<Vector>>& cellInterfaceMidPoints,
                               const std::vector<std::vector<Vector>>& cellNormals,
                               const std::vector<BOUNDARY_TYPE>& cellBoundaryTypes,
                               unsigned nCells,
                               unsigned nDim ) {
    // File has nCells lines, each line is a comma separated entry containing:
    // cellNeighbors (nNodesPerCell elements),
    // cellInterfaceMidPoints (nNodesPerCell x nDim elements),
    // cellNormals (nNodesPerCell x nDim elements),
    // cellBoundaryTypes (1 element), (tranlated from BOUNDARY_TYPE to unsigned)

    std::ofstream outFile( outputFile );
    outFile << std::fixed << setprecision( 15 );
    // const std::size_t bufferSize = 10000;    // Adjust as needed
    // outFile.rdbuf()->pubsetbuf( 0, bufferSize );
    if( outFile.is_open() ) {
        // Write data to the file
        for( unsigned i = 0; i < nCells; ++i ) {
            // Write cellNeighbors
            for( unsigned j = 0; j < cellNeighbors[i].size(); ++j ) {
                outFile << cellNeighbors[i][j];
                if( j < cellNeighbors[i].size() - 1 ) {
                    outFile << ",";
                }
            }
            // Write cellInterfaceMidPoints
            for( unsigned j = 0; j < cellInterfaceMidPoints[i].size(); ++j ) {
                for( unsigned k = 0; k < nDim; ++k ) {
                    outFile << "," << cellInterfaceMidPoints[i][j][k];
                }
            }
            // Write cellNormals
            for( unsigned j = 0; j < cellNormals[i].size(); ++j ) {
                for( unsigned k = 0; k < nDim; ++k ) {
                    outFile << "," << cellNormals[i][j][k];
                }
            }

            // Write cellBoundaryTypes
            outFile << "," << static_cast<unsigned>( cellBoundaryTypes[i] );
            outFile << std::endl;
        }
        outFile.close();
    }
    else {
        ErrorMessages::Error( "Error opening connectivity file.", CURRENT_FUNCTION );
    }
}

void WriteRestartSolution( const std::string& baseOutputFile,
                           const std::vector<double>& solution,
                           const std::vector<double>& scalarFlux,
                           int rank,
                           int idx_iter,
                           double totalAbsorptionHohlraumCenter,
                           double totalAbsorptionHohlraumVertical,
                           double totalAbsorptionHohlraumHorizontal,
                           double totalAbsorption ) {
    // Generate filename with rank number
    std::string outputFile = baseOutputFile + "_restart_rank_" + std::to_string( rank );

    // Check if the file exists and delete it
    if( std::filesystem::exists( outputFile ) ) {
        std::filesystem::remove( outputFile );
    }

    // Open the file for binary writing
    std::ofstream outFile( outputFile, std::ios::out | std::ios::binary );
    if( !outFile ) {
        std::cerr << "Error opening restart solution file." << std::endl;
        return;
    }

    // Write the iteration number as the first item (in binary)
    outFile.write( reinterpret_cast<const char*>( &idx_iter ), sizeof( idx_iter ) );
    outFile.write( reinterpret_cast<const char*>( &totalAbsorptionHohlraumCenter ), sizeof( totalAbsorptionHohlraumCenter ) );
    outFile.write( reinterpret_cast<const char*>( &totalAbsorptionHohlraumVertical ), sizeof( totalAbsorptionHohlraumVertical ) );
    outFile.write( reinterpret_cast<const char*>( &totalAbsorptionHohlraumHorizontal ), sizeof( totalAbsorptionHohlraumHorizontal ) );
    outFile.write( reinterpret_cast<const char*>( &totalAbsorption ), sizeof( totalAbsorption ) );

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

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

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

    // Close the file
    outFile.close();
}

int LoadRestartSolution( const std::string& baseInputFile,
                         std::vector<double>& solution,
                         std::vector<double>& scalarFlux,
                         int rank,
                         unsigned long nCells,
                         double& totalAbsorptionHohlraumCenter,
                         double& totalAbsorptionHohlraumVertical,
                         double& totalAbsorptionHohlraumHorizontal,
                         double& totalAbsorption ) {
    // Generate filename with rank number
    std::string inputFile = baseInputFile + "_restart_rank_" + std::to_string( rank );

    // Open the file for binary reading
    std::ifstream inFile( inputFile, std::ios::in | std::ios::binary );
    if( !inFile ) {
        std::cerr << "Error opening restart solution file for reading." << std::endl;
        return 0;    // Signal failure
    }

    // Read the iteration number
    int idx_iter = 0;
    inFile.read( reinterpret_cast<char*>( &idx_iter ), sizeof( idx_iter ) );
    inFile.read( reinterpret_cast<char*>( &totalAbsorptionHohlraumCenter ), sizeof( totalAbsorptionHohlraumCenter ) );
    inFile.read( reinterpret_cast<char*>( &totalAbsorptionHohlraumVertical ), sizeof( totalAbsorptionHohlraumVertical ) );
    inFile.read( reinterpret_cast<char*>( &totalAbsorptionHohlraumHorizontal ), sizeof( totalAbsorptionHohlraumHorizontal ) );
    inFile.read( reinterpret_cast<char*>( &totalAbsorption ), sizeof( totalAbsorption ) );

    // Read the size of the solution and scalarFlux vectors
    size_t solution_size, scalarFlux_size;
    inFile.read( reinterpret_cast<char*>( &solution_size ), sizeof( solution_size ) );
    inFile.read( reinterpret_cast<char*>( &scalarFlux_size ), sizeof( scalarFlux_size ) );

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

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

    // Close the file
    inFile.close();

    // Verify that we have enough entries to populate scalarFlux
    if( tempData.size() < nCells ) {
        std::cerr << "Not enough data to populate scalar flux vector." << std::endl;
        return 0;    // Signal failure
    }

    // Allocate the last nCells entries to scalarFlux and the rest to solution
    solution.assign( tempData.begin(), tempData.end() - nCells );
    scalarFlux.assign( tempData.end() - nCells, tempData.end() );

    return idx_iter;    // Return the iteration index
}

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

    if( argc < 2 ) {
        std::cout << usage_help;
        exit( EXIT_FAILURE );
    }
    for( int i = 1; i < argc; i++ ) {
        std::string arg = argv[i];
        if( arg == "-h" ) {
            std::cout << usage_help;
            exit( EXIT_SUCCESS );
        }
        else {
            inputFile = std::string( argv[i] );
            std::ifstream f( inputFile );
            if( !f.is_open() ) {
                ErrorMessages::Error( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
            }
        }
    }
    return inputFile;
}

void PrintLogHeader( std::string inputFile ) {
    int nprocs = 1;
    int rank   = 0;
#ifdef IMPORT_MPI
    // Initialize MPI
    MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
    MPI_Comm_rank( MPI_COMM_WORLD, &rank );
#endif
    if( rank == 0 ) {
        auto log = spdlog::get( "event" );

        // New design
        log->info( "------------------------------------------------------------------------" );
        log->info( "|    _  _____ _____     ____ _____                                     |" );
        log->info( "|   | |/ /_ _|_   _|   |  _ \\_   _|                                    |" );
        log->info( "|   | ' / | |  | |_____| |_) || |            Version                   |" );
        log->info( "|   | . \\ | |  | |_____|  _ < | |             1.0.0                    |" );
        log->info( "|   |_|\\_\\___| |_|     |_| \\_\\|_|                                      |" );
        log->info( "|                                                                      |" );
        log->info( "------------------------------------------------------------------------" );
        log->info( "|  Copyright: MIT Licence                                              |" );
        // log->info( "|  Authors: S. Schotthöfer, J. Kusch, J. Wolters, P. Stammer ad T. Xiao |" );
        log->info( "|  Date: 21th April 2022                                               |" );
        log->info( "------------------------------------------------------------------------" );
        log->info( "|" );
        log->info( "| Git commit :\t{0}", GIT_HASH );
        log->info( "| Config file:\t{0}", inputFile );
        log->info( "| MPI Threads:\t{0}", nprocs );
        log->info( "| OMP Threads:\t{0}", omp_get_max_threads() );
        log->info( "|" );
        log->info( "-------------------------- Config File Info ----------------------------" );
        log->info( "|" );
        // print file content while omitting comments
        std::ifstream ifs( inputFile );
        if( ifs.is_open() ) {
            std::string line;
            while( !ifs.eof() ) {
                std::getline( ifs, line );
                if( line[0] != '%' ) log->info( "| {0}", line );
            }
        }
        // log->info( "------------------------------------------------------------------------" );
    }
#ifdef IMPORT_MPI
    MPI_Barrier( MPI_COMM_WORLD );
#endif
}

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

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

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

    std::string pyPath = KITRT_PYTHON_PATH;

    if( !Py_IsInitialized() ) {
        Py_InitializeEx( 0 );
        if( !Py_IsInitialized() ) {
            ErrorMessages::Error( "Python init failed!", CURRENT_FUNCTION );
        }
        PyRun_SimpleString( ( "import sys\nsys.path.append('" + pyPath + "')" ).c_str() );
    }

    PyObject *pArgs, *pReturn, *pModule, *pFunc;
    PyArrayObject* np_ret;

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

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

    pFunc = PyObject_GetAttrString( pModule, "generate" );
    if( !pFunc || !PyCallable_Check( pFunc ) ) {
        PyErr_Print();
        Py_DecRef( pModule );
        Py_DecRef( pFunc );
        ErrorMessages::Error( "'generate' is null or not callable!", CURRENT_FUNCTION );
    }

    pArgs = PyTuple_New( 2 );
    PyTuple_SetItem( pArgs, 0, reinterpret_cast<PyObject*>( image ) );
    PyTuple_SetItem( pArgs, 1, reinterpret_cast<PyObject*>( su2 ) );
    pReturn = PyObject_CallObject( pFunc, pArgs );
    np_ret  = reinterpret_cast<PyArrayObject*>( pReturn );

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

    Matrix gsImage( m, n, c_out );

    // Finalizing
    Py_DecRef( pFunc );
    Py_DecRef( pModule );
    Py_DECREF( np_ret );

    return gsImage.transpose();
}
*/

CSMMLab / KiT-RT / #95

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

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