diff --git a/.github/workflows/c-cpp.yml b/.github/workflows/c-cpp.yml index b300b81a..be28f9aa 100644 --- a/.github/workflows/c-cpp.yml +++ b/.github/workflows/c-cpp.yml @@ -14,6 +14,10 @@ jobs: env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }} + OMP_NUM_THREADS: 1 + OMP_DYNAMIC: FALSE + OPENBLAS_NUM_THREADS: 1 + MKL_NUM_THREADS: 1 steps: - uses: actions/checkout@v2 @@ -30,8 +34,6 @@ jobs: cd build ./unit_tests - OMP_NUM_THREADS=1 ./unit_tests -r junit > unit_tests_report.xml - - name: Code coverage run: | - cpp-coveralls -r . -b "build/" -i "src/" -i "include/" --exclude "ext/" --gcov-options "\-lp" --verbose \ No newline at end of file + cpp-coveralls -r . -b "build/" -i "src/" -i "include/" --exclude "ext/" --gcov-options "\-lp" --verbose diff --git a/CMakeLists.txt b/CMakeLists.txt index 58470146..bf1fee28 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,7 +45,7 @@ include_directories( ${LAPACK_INCLUDE_DIR} ) find_package( BLAS ) -find_package( VTK COMPONENTS vtkIOGeometry vtkFiltersCore REQUIRED ) +find_package( VTK COMPONENTS IOGeometry FiltersCore REQUIRED ) add_compile_definitions( KITRT_PYTHON_PATH="${CMAKE_SOURCE_DIR}/python" ) add_compile_definitions( BLAZE_USE_SHARED_MEMORY_PARALLELIZATION=0 ) diff --git a/include/common/mesh.hpp b/include/common/mesh.hpp index 92f42f41..6ae11686 100644 --- a/include/common/mesh.hpp +++ b/include/common/mesh.hpp @@ -133,7 +133,7 @@ class Mesh * @return cell_idxs: std::vector */ std::vector GetCellsofBall( const double x, const double y, const double r ) const; - /*! @brief Returns index of cells contained in the rectangle with specified corner coordinates/* + /*! @brief Returns index of cells contained in the rectangle with specified corner coordinates * @return cell_idxs: std::vector */ std::vector GetCellsofRectangle( const std::vector>& cornercoordinates ) const; diff --git a/include/optimizers/neuralnetworkoptimizer.hpp b/include/optimizers/neuralnetworkoptimizer.hpp index cad5ca9c..e44805b3 100644 --- a/include/optimizers/neuralnetworkoptimizer.hpp +++ b/include/optimizers/neuralnetworkoptimizer.hpp @@ -71,15 +71,15 @@ class NeuralNetworkOptimizer : public OptimizerBase inline ~NeuralNetworkOptimizer(); - inline void Solve( Vector& alpha, const Vector& u, const VectorVector& moments, unsigned idx_cell = 0 ) override{}; + inline void Solve( Vector&, const Vector&, const VectorVector&, unsigned = 0 ) override{}; - inline void SolveMultiCell( VectorVector& alpha, const VectorVector& u, const VectorVector& moments, Vector& alpha_norms ) override{}; + inline void SolveMultiCell( VectorVector&, const VectorVector&, const VectorVector&, Vector& ) override{}; /*! @brief Reconstruct the moment sol from the Lagrange multiplier alpha * @param sol moment vector * @param alpha Lagrange multipliers * @param moments Moment basis */ - inline void ReconstructMoments( Vector& sol, const Vector& alpha, const VectorVector& moments ) override{}; + inline void ReconstructMoments( Vector&, const Vector&, const VectorVector& ) override{}; }; #endif #endif // MLOPTIMIZER_H diff --git a/include/solvers/snsolver_hpc.hpp b/include/solvers/snsolver_hpc.hpp index 410c5086..da61e072 100644 --- a/include/solvers/snsolver_hpc.hpp +++ b/include/solvers/snsolver_hpc.hpp @@ -158,7 +158,8 @@ class SNSolverHPC void IterPostprocessing(); - void SetGhostCells(); /*!< @brief Sets vector of ghost cells for + /*! @brief Sets vector of ghost cell values according to boundary conditions */ + void SetGhostCells(); // Helper /*! @brief ComputeTimeStep calculates the maximal stable time step using the diff --git a/src/common/config.cpp b/src/common/config.cpp index 01f87ea6..246efdb3 100644 --- a/src/common/config.cpp +++ b/src/common/config.cpp @@ -530,7 +530,7 @@ void Config::SetPostprocessing() { _dataDir = std::filesystem::path( _inputDir ).append( _dataDir ).lexically_normal(); // create directories if they dont exist - if( !std::filesystem::exists( _outputDir ) ) std::filesystem::create_directory( _outputDir ); + if( !std::filesystem::exists( _outputDir ) ) std::filesystem::create_directories( _outputDir ); // init logger InitLogger(); @@ -1249,7 +1249,7 @@ void Config::InitLogger() { // create log dir if not existent if( !std::filesystem::exists( _logDir ) ) { - std::filesystem::create_directory( _logDir ); + std::filesystem::create_directories( _logDir ); } if( !spdlog::get( "event" ) ) { @@ -1398,4 +1398,4 @@ template K Config::findKey( const std::map& myMap } // If the value is not found, you can return a default value or throw an exception throw std::out_of_range( "Value not found in the map" ); -} \ No newline at end of file +} diff --git a/src/common/io.cpp b/src/common/io.cpp index aa10c89b..a8042c51 100644 --- a/src/common/io.cpp +++ b/src/common/io.cpp @@ -51,11 +51,9 @@ void ExportVTK( const std::string fileName, const std::vector>>& outputFields, const std::vector>& outputFieldNames, const Mesh* mesh ) { - int rank = 0; - 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 ) { @@ -152,7 +150,8 @@ void ExportVTK( const std::string fileName, Mesh* LoadSU2MeshFromFile( const Config* settings ) { auto log = spdlog::get( "event" ); // log->info( "| Importing mesh. This may take a while for large meshes." ); - unsigned dim; + unsigned dim = settings->GetDim(); + bool foundDim = false; std::vector nodes; std::vector> cells; std::vector>> boundaries; @@ -166,12 +165,16 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) { while( getline( ifs, line ) ) { if( line.find( "NDIME", 0 ) != std::string::npos ) { dim = static_cast( TextProcessingToolbox::GetTrailingNumber( line ) ); + foundDim = true; // if( settings->GetDim() != dim ) { // log->info( "Warning: Mesh dimension does not coinside with problem dimension! Proceed with caution!" ); // } break; } } + if( !foundDim ) { + ErrorMessages::Error( "Invalid mesh file detected! Missing NDIME entry.", CURRENT_FUNCTION ); + } ifs.clear(); ifs.seekg( 0, std::ios::beg ); while( getline( ifs, line ) ) { @@ -682,4 +685,4 @@ Matrix createSU2MeshFromImage( std::string imageName, std::string SU2Filename ) return gsImage.transpose(); } -*/ \ No newline at end of file +*/ diff --git a/src/common/mesh.cpp b/src/common/mesh.cpp index ed7ef3c3..594f1ae3 100644 --- a/src/common/mesh.cpp +++ b/src/common/mesh.cpp @@ -26,10 +26,8 @@ Mesh::Mesh( const Config* settings, else { ErrorMessages::Error( "Unsupported mesh dimension!", CURRENT_FUNCTION ); } - int nprocs = 1; int rank = 0; #ifdef IMPORT_MPI - MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); #endif if( rank == 0 ) { @@ -93,18 +91,13 @@ Mesh::Mesh( const Config* settings, Mesh::~Mesh() {} void Mesh::ComputeConnectivity() { - int nprocs = 1; int rank = 0; #ifdef IMPORT_MPI - MPI_Comm_size( MPI_COMM_WORLD, &nprocs ); MPI_Comm_rank( MPI_COMM_WORLD, &rank ); #endif - unsigned comm_size = 1; // No MPI implementation right now // determine number/chunk size and indices of cells treated by each mpi thread (deactivated for now) unsigned chunkSize = _numCells; // std::ceil( static_cast( _numCells ) / static_cast( comm_size ) ); - unsigned mpiCellStart = 0; // comm_rank * chunkSize; - unsigned mpiCellEnd = _numCells; // std::min( ( comm_rank + 1 ) * chunkSize, _numCells ); // 'flat' vectors are a flattened representation of the neighbors numCells nested vectors; easier for MPI // 'part' vectors store information for each single MPI thread diff --git a/src/datagenerator/datageneratorbase.cpp b/src/datagenerator/datageneratorbase.cpp index d8b40dfc..5197e5be 100644 --- a/src/datagenerator/datageneratorbase.cpp +++ b/src/datagenerator/datageneratorbase.cpp @@ -85,8 +85,10 @@ DataGeneratorBase::DataGeneratorBase( Config* settings ) { } DataGeneratorBase::~DataGeneratorBase() { - delete _quadrature; delete _entropy; + delete _optimizer; + delete _basisGenerator; + delete _quadrature; } DataGeneratorBase* DataGeneratorBase::Create( Config* settings ) { diff --git a/src/solvers/mnsolver_normalized.cpp b/src/solvers/mnsolver_normalized.cpp index fd1ef480..9a6e4986 100644 --- a/src/solvers/mnsolver_normalized.cpp +++ b/src/solvers/mnsolver_normalized.cpp @@ -29,7 +29,7 @@ MNSolverNormalized::MNSolverNormalized( Config* settings ) : MNSolver( settings _sol2 = VectorVector( _nCells, Vector( _nSystem, 0.0 ) ); } -MNSolverNormalized::~MNSolverNormalized() {} +MNSolverNormalized::~MNSolverNormalized() { delete _optimizer2; } void MNSolverNormalized::IterPreprocessing( unsigned /*idx_pseudotime*/ ) { Vector alpha_norm_per_cell( _nCells, 0 ); // ONLY FOR DEBUGGING! THIS SLOWS DOWN THE CODE diff --git a/src/solvers/snsolver_hpc.cpp b/src/solvers/snsolver_hpc.cpp index 73b20619..a2756c6e 100644 --- a/src/solvers/snsolver_hpc.cpp +++ b/src/solvers/snsolver_hpc.cpp @@ -38,7 +38,7 @@ SNSolverHPC::SNSolverHPC( Config* settings ) { _nq = static_cast( quad->GetNq() ); _nSys = _nq; - if( _numProcs > _nq ) { + if( static_cast( _numProcs ) > _nq ) { ErrorMessages::Error( "The number of processors must be less than or equal to the number of quadrature points.", CURRENT_FUNCTION ); } @@ -287,7 +287,7 @@ void SNSolverHPC::Solve() { DrawPreSolverOutput(); } _curSimTime = 0.0; - auto start = std::chrono::high_resolution_clock::now(); // Start timing + auto start = std::chrono::high_resolution_clock::now(); // Start timing std::chrono::duration duration; // Loop over energies (pseudo-time of continuous slowing down approach) @@ -1487,7 +1487,6 @@ void SNSolverHPC::SetProbingCellsLineGreen() { double horizontalLineWidth = std::abs( p1[0] - p2[0] ); double pt_ratio_h = horizontalLineWidth / ( horizontalLineWidth + verticalLineWidth ); - double pt_ratio_v = verticalLineWidth / ( horizontalLineWidth + verticalLineWidth ); unsigned nHorizontalProbingCells = (unsigned)std::ceil( _nProbingCellsLineGreen / 2 * pt_ratio_h ); unsigned nVerticalProbingCells = _nProbingCellsLineGreen / 2 - nHorizontalProbingCells; @@ -1591,20 +1590,13 @@ void SNSolverHPC::SetProbingCellsLineGreen() { } // Compute the probing cells for each block - for( int i = 0; i < _nProbingCellsBlocksGreen; i++ ) { + for( unsigned i = 0; i < _nProbingCellsBlocksGreen; i++ ) { _probingCellsBlocksGreen.push_back( _mesh->GetCellsofRectangle( block_corners[i] ) ); } } } void SNSolverHPC::ComputeQOIsGreenProbingLine() { - double verticalLineWidth = std::abs( _cornerUpperLeftGreen[1] - _cornerLowerLeftGreen[1] - _thicknessGreen ); - double horizontalLineWidth = std::abs( _cornerUpperLeftGreen[0] - _cornerUpperRightGreen[0] - _thicknessGreen ); - - double dl = 2 * ( horizontalLineWidth + verticalLineWidth ) / ( (double)_nProbingCellsLineGreen ); - double area = dl * _thicknessGreen; - double l_max = _nProbingCellsLineGreen * dl; - #pragma omp parallel for for( unsigned i = 0; i < _nProbingCellsLineGreen; i++ ) { // Loop over probing cells _absorptionValsLineSegment[i] = @@ -1694,4 +1686,4 @@ void SNSolverHPC::ComputeCellsPerimeterLattice() { } } } -} \ No newline at end of file +} diff --git a/tests/input/unit_tests/result/logs/2026-02-10_15:15:38 b/tests/input/unit_tests/result/logs/2026-02-10_15:15:38 new file mode 100644 index 00000000..e69de29b diff --git a/tests/input/unit_tests/result/logs/2026-02-10_15:15:38.csv b/tests/input/unit_tests/result/logs/2026-02-10_15:15:38.csv new file mode 100644 index 00000000..e69de29b diff --git a/tests/input/unit_tests/result/logs/2026-02-10_15:15:381.csv b/tests/input/unit_tests/result/logs/2026-02-10_15:15:381.csv new file mode 100644 index 00000000..e69de29b diff --git a/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200.cfg b/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200.cfg new file mode 100644 index 00000000..a2808cf7 --- /dev/null +++ b/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200.cfg @@ -0,0 +1,45 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Lattice Validation SN-HPC % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ---- File specifications ---- +% +OUTPUT_DIR = ../../../result +OUTPUT_FILE = lattice_hpc_200 +LOG_DIR = ../../../result/logs +LOG_FILE = lattice_hpc_200_output +MESH_FILE = lattice_hpc_200.mesh +FORCE_CONNECTIVITY_RECOMPUTE = YES +% +% ---- Problem specifications ---- +% +PROBLEM = LATTICE +SPATIAL_DIM = 2 +SOURCE_MAGNITUDE = 1.0 +% +% ---- Solver specifications ---- +% +HPC_SOLVER = YES +SOLVER = SN_SOLVER +CFL_NUMBER = 0.5 +TIME_FINAL = 0.5 +RECONS_ORDER = 2 +TIME_INTEGRATION_ORDER = 2 +% +% ---- Boundary Conditions ---- +% +BC_NEUMANN = ( void ) +% +% ---- Quadrature ---- +% +QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED_2D +QUAD_ORDER = 6 +% +% ----- Output ---- +% +VOLUME_OUTPUT = (MINIMAL) +VOLUME_OUTPUT_FREQUENCY = 0 +SCREEN_OUTPUT = (ITER, WALL_TIME, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT) +SCREEN_OUTPUT_FREQUENCY = 1 +HISTORY_OUTPUT = (ITER, SIM_TIME, WALL_TIME, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT, CUR_OUTFLOW, TOTAL_OUTFLOW, CUR_OUTFLOW_P1, TOTAL_OUTFLOW_P1, CUR_OUTFLOW_P2, TOTAL_OUTFLOW_P2, MAX_OUTFLOW, CUR_PARTICLE_ABSORPTION, TOTAL_PARTICLE_ABSORPTION, MAX_PARTICLE_ABSORPTION) +HISTORY_OUTPUT_FREQUENCY = 1 diff --git a/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200.mesh b/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200.mesh new file mode 100644 index 00000000..1bc14db2 --- /dev/null +++ b/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200.mesh @@ -0,0 +1,367 @@ +NDIME= 2 +NELEM= 200 +5 0 1 12 0 +5 0 12 11 1 +5 1 2 13 2 +5 1 13 12 3 +5 2 3 14 4 +5 2 14 13 5 +5 3 4 15 6 +5 3 15 14 7 +5 4 5 16 8 +5 4 16 15 9 +5 5 6 17 10 +5 5 17 16 11 +5 6 7 18 12 +5 6 18 17 13 +5 7 8 19 14 +5 7 19 18 15 +5 8 9 20 16 +5 8 20 19 17 +5 9 10 21 18 +5 9 21 20 19 +5 11 12 23 20 +5 11 23 22 21 +5 12 13 24 22 +5 12 24 23 23 +5 13 14 25 24 +5 13 25 24 25 +5 14 15 26 26 +5 14 26 25 27 +5 15 16 27 28 +5 15 27 26 29 +5 16 17 28 30 +5 16 28 27 31 +5 17 18 29 32 +5 17 29 28 33 +5 18 19 30 34 +5 18 30 29 35 +5 19 20 31 36 +5 19 31 30 37 +5 20 21 32 38 +5 20 32 31 39 +5 22 23 34 40 +5 22 34 33 41 +5 23 24 35 42 +5 23 35 34 43 +5 24 25 36 44 +5 24 36 35 45 +5 25 26 37 46 +5 25 37 36 47 +5 26 27 38 48 +5 26 38 37 49 +5 27 28 39 50 +5 27 39 38 51 +5 28 29 40 52 +5 28 40 39 53 +5 29 30 41 54 +5 29 41 40 55 +5 30 31 42 56 +5 30 42 41 57 +5 31 32 43 58 +5 31 43 42 59 +5 33 34 45 60 +5 33 45 44 61 +5 34 35 46 62 +5 34 46 45 63 +5 35 36 47 64 +5 35 47 46 65 +5 36 37 48 66 +5 36 48 47 67 +5 37 38 49 68 +5 37 49 48 69 +5 38 39 50 70 +5 38 50 49 71 +5 39 40 51 72 +5 39 51 50 73 +5 40 41 52 74 +5 40 52 51 75 +5 41 42 53 76 +5 41 53 52 77 +5 42 43 54 78 +5 42 54 53 79 +5 44 45 56 80 +5 44 56 55 81 +5 45 46 57 82 +5 45 57 56 83 +5 46 47 58 84 +5 46 58 57 85 +5 47 48 59 86 +5 47 59 58 87 +5 48 49 60 88 +5 48 60 59 89 +5 49 50 61 90 +5 49 61 60 91 +5 50 51 62 92 +5 50 62 61 93 +5 51 52 63 94 +5 51 63 62 95 +5 52 53 64 96 +5 52 64 63 97 +5 53 54 65 98 +5 53 65 64 99 +5 55 56 67 100 +5 55 67 66 101 +5 56 57 68 102 +5 56 68 67 103 +5 57 58 69 104 +5 57 69 68 105 +5 58 59 70 106 +5 58 70 69 107 +5 59 60 71 108 +5 59 71 70 109 +5 60 61 72 110 +5 60 72 71 111 +5 61 62 73 112 +5 61 73 72 113 +5 62 63 74 114 +5 62 74 73 115 +5 63 64 75 116 +5 63 75 74 117 +5 64 65 76 118 +5 64 76 75 119 +5 66 67 78 120 +5 66 78 77 121 +5 67 68 79 122 +5 67 79 78 123 +5 68 69 80 124 +5 68 80 79 125 +5 69 70 81 126 +5 69 81 80 127 +5 70 71 82 128 +5 70 82 81 129 +5 71 72 83 130 +5 71 83 82 131 +5 72 73 84 132 +5 72 84 83 133 +5 73 74 85 134 +5 73 85 84 135 +5 74 75 86 136 +5 74 86 85 137 +5 75 76 87 138 +5 75 87 86 139 +5 77 78 89 140 +5 77 89 88 141 +5 78 79 90 142 +5 78 90 89 143 +5 79 80 91 144 +5 79 91 90 145 +5 80 81 92 146 +5 80 92 91 147 +5 81 82 93 148 +5 81 93 92 149 +5 82 83 94 150 +5 82 94 93 151 +5 83 84 95 152 +5 83 95 94 153 +5 84 85 96 154 +5 84 96 95 155 +5 85 86 97 156 +5 85 97 96 157 +5 86 87 98 158 +5 86 98 97 159 +5 88 89 100 160 +5 88 100 99 161 +5 89 90 101 162 +5 89 101 100 163 +5 90 91 102 164 +5 90 102 101 165 +5 91 92 103 166 +5 91 103 102 167 +5 92 93 104 168 +5 92 104 103 169 +5 93 94 105 170 +5 93 105 104 171 +5 94 95 106 172 +5 94 106 105 173 +5 95 96 107 174 +5 95 107 106 175 +5 96 97 108 176 +5 96 108 107 177 +5 97 98 109 178 +5 97 109 108 179 +5 99 100 111 180 +5 99 111 110 181 +5 100 101 112 182 +5 100 112 111 183 +5 101 102 113 184 +5 101 113 112 185 +5 102 103 114 186 +5 102 114 113 187 +5 103 104 115 188 +5 103 115 114 189 +5 104 105 116 190 +5 104 116 115 191 +5 105 106 117 192 +5 105 117 116 193 +5 106 107 118 194 +5 106 118 117 195 +5 107 108 119 196 +5 107 119 118 197 +5 108 109 120 198 +5 108 120 119 199 +NPOIN= 121 +-3.5 -3.5 0 +-2.8 -3.5 1 +-2.1 -3.5 2 +-1.4 -3.5 3 +-0.7 -3.5 4 +0 -3.5 5 +0.7 -3.5 6 +1.4 -3.5 7 +2.1 -3.5 8 +2.8 -3.5 9 +3.5 -3.5 10 +-3.5 -2.8 11 +-2.8 -2.8 12 +-2.1 -2.8 13 +-1.4 -2.8 14 +-0.7 -2.8 15 +0 -2.8 16 +0.7 -2.8 17 +1.4 -2.8 18 +2.1 -2.8 19 +2.8 -2.8 20 +3.5 -2.8 21 +-3.5 -2.1 22 +-2.8 -2.1 23 +-2.1 -2.1 24 +-1.4 -2.1 25 +-0.7 -2.1 26 +0 -2.1 27 +0.7 -2.1 28 +1.4 -2.1 29 +2.1 -2.1 30 +2.8 -2.1 31 +3.5 -2.1 32 +-3.5 -1.4 33 +-2.8 -1.4 34 +-2.1 -1.4 35 +-1.4 -1.4 36 +-0.7 -1.4 37 +0 -1.4 38 +0.7 -1.4 39 +1.4 -1.4 40 +2.1 -1.4 41 +2.8 -1.4 42 +3.5 -1.4 43 +-3.5 -0.7 44 +-2.8 -0.7 45 +-2.1 -0.7 46 +-1.4 -0.7 47 +-0.7 -0.7 48 +0 -0.7 49 +0.7 -0.7 50 +1.4 -0.7 51 +2.1 -0.7 52 +2.8 -0.7 53 +3.5 -0.7 54 +-3.5 0 55 +-2.8 0 56 +-2.1 0 57 +-1.4 0 58 +-0.7 0 59 +0 0 60 +0.7 0 61 +1.4 0 62 +2.1 0 63 +2.8 0 64 +3.5 0 65 +-3.5 0.7 66 +-2.8 0.7 67 +-2.1 0.7 68 +-1.4 0.7 69 +-0.7 0.7 70 +0 0.7 71 +0.7 0.7 72 +1.4 0.7 73 +2.1 0.7 74 +2.8 0.7 75 +3.5 0.7 76 +-3.5 1.4 77 +-2.8 1.4 78 +-2.1 1.4 79 +-1.4 1.4 80 +-0.7 1.4 81 +0 1.4 82 +0.7 1.4 83 +1.4 1.4 84 +2.1 1.4 85 +2.8 1.4 86 +3.5 1.4 87 +-3.5 2.1 88 +-2.8 2.1 89 +-2.1 2.1 90 +-1.4 2.1 91 +-0.7 2.1 92 +0 2.1 93 +0.7 2.1 94 +1.4 2.1 95 +2.1 2.1 96 +2.8 2.1 97 +3.5 2.1 98 +-3.5 2.8 99 +-2.8 2.8 100 +-2.1 2.8 101 +-1.4 2.8 102 +-0.7 2.8 103 +0 2.8 104 +0.7 2.8 105 +1.4 2.8 106 +2.1 2.8 107 +2.8 2.8 108 +3.5 2.8 109 +-3.5 3.5 110 +-2.8 3.5 111 +-2.1 3.5 112 +-1.4 3.5 113 +-0.7 3.5 114 +0 3.5 115 +0.7 3.5 116 +1.4 3.5 117 +2.1 3.5 118 +2.8 3.5 119 +3.5 3.5 120 +NMARK= 1 +MARKER_TAG= void +MARKER_ELEMS= 40 +3 0 1 +3 1 2 +3 2 3 +3 3 4 +3 4 5 +3 5 6 +3 6 7 +3 7 8 +3 8 9 +3 9 10 +3 10 21 +3 21 32 +3 32 43 +3 43 54 +3 54 65 +3 65 76 +3 76 87 +3 87 98 +3 98 109 +3 109 120 +3 120 119 +3 119 118 +3 118 117 +3 117 116 +3 116 115 +3 115 114 +3 114 113 +3 113 112 +3 112 111 +3 111 110 +3 110 99 +3 99 88 +3 88 77 +3 77 66 +3 66 55 +3 55 44 +3 44 33 +3 33 22 +3 22 11 +3 11 0 diff --git a/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200_csv_reference b/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200_csv_reference new file mode 100644 index 00000000..ef356ba0 --- /dev/null +++ b/tests/input/validation_tests/SN_solver_hpc/lattice_hpc_200_csv_reference @@ -0,0 +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 diff --git a/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200.cfg b/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200.cfg new file mode 100644 index 00000000..49d83aeb --- /dev/null +++ b/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200.cfg @@ -0,0 +1,45 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Symmetric Hohlraum Validation SN-HPC% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ---- File specifications ---- +% +OUTPUT_DIR = ../../../result +OUTPUT_FILE = symmetric_hohlraum_hpc_200 +LOG_DIR = ../../../result/logs +LOG_FILE = symmetric_hohlraum_hpc_200_output +MESH_FILE = symmetric_hohlraum_hpc_200.mesh +FORCE_CONNECTIVITY_RECOMPUTE = YES +% +% ---- Problem specifications ---- +% +PROBLEM = SYMMETRIC_HOHLRAUM +SPATIAL_DIM = 2 +N_SAMPLING_PTS_LINE_GREEN = 20 +% +% ---- Solver specifications ---- +% +HPC_SOLVER = YES +SOLVER = SN_SOLVER +CFL_NUMBER = 0.5 +TIME_FINAL = 0.2 +RECONS_ORDER = 2 +TIME_INTEGRATION_ORDER = 2 +% +% ---- Boundary Conditions ---- +% +BC_NEUMANN = ( void ) +% +% ---- Quadrature ---- +% +QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED_2D +QUAD_ORDER = 6 +% +% ----- Output ---- +% +VOLUME_OUTPUT = (MINIMAL) +VOLUME_OUTPUT_FREQUENCY = 0 +SCREEN_OUTPUT = (ITER, WALL_TIME, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT) +SCREEN_OUTPUT_FREQUENCY = 1 +HISTORY_OUTPUT = (ITER, SIM_TIME, WALL_TIME, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT, CUR_OUTFLOW, TOTAL_OUTFLOW, MAX_OUTFLOW, TOTAL_PARTICLE_ABSORPTION, TOTAL_PARTICLE_ABSORPTION_CENTER, TOTAL_PARTICLE_ABSORPTION_VERTICAL, TOTAL_PARTICLE_ABSORPTION_HORIZONTAL, PROBE_MOMENT_TIME_TRACE, VAR_ABSORPTION_GREEN, ABSORPTION_GREEN_BLOCK, ABSORPTION_GREEN_LINE) +HISTORY_OUTPUT_FREQUENCY = 1 diff --git a/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200.mesh b/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200.mesh new file mode 100644 index 00000000..1595178b --- /dev/null +++ b/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200.mesh @@ -0,0 +1,367 @@ +NDIME= 2 +NELEM= 200 +5 0 1 12 0 +5 0 12 11 1 +5 1 2 13 2 +5 1 13 12 3 +5 2 3 14 4 +5 2 14 13 5 +5 3 4 15 6 +5 3 15 14 7 +5 4 5 16 8 +5 4 16 15 9 +5 5 6 17 10 +5 5 17 16 11 +5 6 7 18 12 +5 6 18 17 13 +5 7 8 19 14 +5 7 19 18 15 +5 8 9 20 16 +5 8 20 19 17 +5 9 10 21 18 +5 9 21 20 19 +5 11 12 23 20 +5 11 23 22 21 +5 12 13 24 22 +5 12 24 23 23 +5 13 14 25 24 +5 13 25 24 25 +5 14 15 26 26 +5 14 26 25 27 +5 15 16 27 28 +5 15 27 26 29 +5 16 17 28 30 +5 16 28 27 31 +5 17 18 29 32 +5 17 29 28 33 +5 18 19 30 34 +5 18 30 29 35 +5 19 20 31 36 +5 19 31 30 37 +5 20 21 32 38 +5 20 32 31 39 +5 22 23 34 40 +5 22 34 33 41 +5 23 24 35 42 +5 23 35 34 43 +5 24 25 36 44 +5 24 36 35 45 +5 25 26 37 46 +5 25 37 36 47 +5 26 27 38 48 +5 26 38 37 49 +5 27 28 39 50 +5 27 39 38 51 +5 28 29 40 52 +5 28 40 39 53 +5 29 30 41 54 +5 29 41 40 55 +5 30 31 42 56 +5 30 42 41 57 +5 31 32 43 58 +5 31 43 42 59 +5 33 34 45 60 +5 33 45 44 61 +5 34 35 46 62 +5 34 46 45 63 +5 35 36 47 64 +5 35 47 46 65 +5 36 37 48 66 +5 36 48 47 67 +5 37 38 49 68 +5 37 49 48 69 +5 38 39 50 70 +5 38 50 49 71 +5 39 40 51 72 +5 39 51 50 73 +5 40 41 52 74 +5 40 52 51 75 +5 41 42 53 76 +5 41 53 52 77 +5 42 43 54 78 +5 42 54 53 79 +5 44 45 56 80 +5 44 56 55 81 +5 45 46 57 82 +5 45 57 56 83 +5 46 47 58 84 +5 46 58 57 85 +5 47 48 59 86 +5 47 59 58 87 +5 48 49 60 88 +5 48 60 59 89 +5 49 50 61 90 +5 49 61 60 91 +5 50 51 62 92 +5 50 62 61 93 +5 51 52 63 94 +5 51 63 62 95 +5 52 53 64 96 +5 52 64 63 97 +5 53 54 65 98 +5 53 65 64 99 +5 55 56 67 100 +5 55 67 66 101 +5 56 57 68 102 +5 56 68 67 103 +5 57 58 69 104 +5 57 69 68 105 +5 58 59 70 106 +5 58 70 69 107 +5 59 60 71 108 +5 59 71 70 109 +5 60 61 72 110 +5 60 72 71 111 +5 61 62 73 112 +5 61 73 72 113 +5 62 63 74 114 +5 62 74 73 115 +5 63 64 75 116 +5 63 75 74 117 +5 64 65 76 118 +5 64 76 75 119 +5 66 67 78 120 +5 66 78 77 121 +5 67 68 79 122 +5 67 79 78 123 +5 68 69 80 124 +5 68 80 79 125 +5 69 70 81 126 +5 69 81 80 127 +5 70 71 82 128 +5 70 82 81 129 +5 71 72 83 130 +5 71 83 82 131 +5 72 73 84 132 +5 72 84 83 133 +5 73 74 85 134 +5 73 85 84 135 +5 74 75 86 136 +5 74 86 85 137 +5 75 76 87 138 +5 75 87 86 139 +5 77 78 89 140 +5 77 89 88 141 +5 78 79 90 142 +5 78 90 89 143 +5 79 80 91 144 +5 79 91 90 145 +5 80 81 92 146 +5 80 92 91 147 +5 81 82 93 148 +5 81 93 92 149 +5 82 83 94 150 +5 82 94 93 151 +5 83 84 95 152 +5 83 95 94 153 +5 84 85 96 154 +5 84 96 95 155 +5 85 86 97 156 +5 85 97 96 157 +5 86 87 98 158 +5 86 98 97 159 +5 88 89 100 160 +5 88 100 99 161 +5 89 90 101 162 +5 89 101 100 163 +5 90 91 102 164 +5 90 102 101 165 +5 91 92 103 166 +5 91 103 102 167 +5 92 93 104 168 +5 92 104 103 169 +5 93 94 105 170 +5 93 105 104 171 +5 94 95 106 172 +5 94 106 105 173 +5 95 96 107 174 +5 95 107 106 175 +5 96 97 108 176 +5 96 108 107 177 +5 97 98 109 178 +5 97 109 108 179 +5 99 100 111 180 +5 99 111 110 181 +5 100 101 112 182 +5 100 112 111 183 +5 101 102 113 184 +5 101 113 112 185 +5 102 103 114 186 +5 102 114 113 187 +5 103 104 115 188 +5 103 115 114 189 +5 104 105 116 190 +5 104 116 115 191 +5 105 106 117 192 +5 105 117 116 193 +5 106 107 118 194 +5 106 118 117 195 +5 107 108 119 196 +5 107 119 118 197 +5 108 109 120 198 +5 108 120 119 199 +NPOIN= 121 +-0.65 -0.65 0 +-0.52 -0.65 1 +-0.39 -0.65 2 +-0.26 -0.65 3 +-0.13 -0.65 4 +0 -0.65 5 +0.13 -0.65 6 +0.26 -0.65 7 +0.39 -0.65 8 +0.52 -0.65 9 +0.65 -0.65 10 +-0.65 -0.52 11 +-0.52 -0.52 12 +-0.39 -0.52 13 +-0.26 -0.52 14 +-0.13 -0.52 15 +0 -0.52 16 +0.13 -0.52 17 +0.26 -0.52 18 +0.39 -0.52 19 +0.52 -0.52 20 +0.65 -0.52 21 +-0.65 -0.39 22 +-0.52 -0.39 23 +-0.39 -0.39 24 +-0.26 -0.39 25 +-0.13 -0.39 26 +0 -0.39 27 +0.13 -0.39 28 +0.26 -0.39 29 +0.39 -0.39 30 +0.52 -0.39 31 +0.65 -0.39 32 +-0.65 -0.26 33 +-0.52 -0.26 34 +-0.39 -0.26 35 +-0.26 -0.26 36 +-0.13 -0.26 37 +0 -0.26 38 +0.13 -0.26 39 +0.26 -0.26 40 +0.39 -0.26 41 +0.52 -0.26 42 +0.65 -0.26 43 +-0.65 -0.13 44 +-0.52 -0.13 45 +-0.39 -0.13 46 +-0.26 -0.13 47 +-0.13 -0.13 48 +0 -0.13 49 +0.13 -0.13 50 +0.26 -0.13 51 +0.39 -0.13 52 +0.52 -0.13 53 +0.65 -0.13 54 +-0.65 0 55 +-0.52 0 56 +-0.39 0 57 +-0.26 0 58 +-0.13 0 59 +0 0 60 +0.13 0 61 +0.26 0 62 +0.39 0 63 +0.52 0 64 +0.65 0 65 +-0.65 0.13 66 +-0.52 0.13 67 +-0.39 0.13 68 +-0.26 0.13 69 +-0.13 0.13 70 +0 0.13 71 +0.13 0.13 72 +0.26 0.13 73 +0.39 0.13 74 +0.52 0.13 75 +0.65 0.13 76 +-0.65 0.26 77 +-0.52 0.26 78 +-0.39 0.26 79 +-0.26 0.26 80 +-0.13 0.26 81 +0 0.26 82 +0.13 0.26 83 +0.26 0.26 84 +0.39 0.26 85 +0.52 0.26 86 +0.65 0.26 87 +-0.65 0.39 88 +-0.52 0.39 89 +-0.39 0.39 90 +-0.26 0.39 91 +-0.13 0.39 92 +0 0.39 93 +0.13 0.39 94 +0.26 0.39 95 +0.39 0.39 96 +0.52 0.39 97 +0.65 0.39 98 +-0.65 0.52 99 +-0.52 0.52 100 +-0.39 0.52 101 +-0.26 0.52 102 +-0.13 0.52 103 +0 0.52 104 +0.13 0.52 105 +0.26 0.52 106 +0.39 0.52 107 +0.52 0.52 108 +0.65 0.52 109 +-0.65 0.65 110 +-0.52 0.65 111 +-0.39 0.65 112 +-0.26 0.65 113 +-0.13 0.65 114 +0 0.65 115 +0.13 0.65 116 +0.26 0.65 117 +0.39 0.65 118 +0.52 0.65 119 +0.65 0.65 120 +NMARK= 1 +MARKER_TAG= void +MARKER_ELEMS= 40 +3 0 1 +3 1 2 +3 2 3 +3 3 4 +3 4 5 +3 5 6 +3 6 7 +3 7 8 +3 8 9 +3 9 10 +3 10 21 +3 21 32 +3 32 43 +3 43 54 +3 54 65 +3 65 76 +3 76 87 +3 87 98 +3 98 109 +3 109 120 +3 120 119 +3 119 118 +3 118 117 +3 117 116 +3 116 115 +3 115 114 +3 114 113 +3 113 112 +3 112 111 +3 111 110 +3 110 99 +3 99 88 +3 88 77 +3 77 66 +3 66 55 +3 55 44 +3 44 33 +3 33 22 +3 22 11 +3 11 0 diff --git a/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200_csv_reference b/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200_csv_reference new file mode 100644 index 00000000..6846f20b --- /dev/null +++ b/tests/input/validation_tests/SN_solver_hpc/symmetric_hohlraum_hpc_200_csv_reference @@ -0,0 +1,6 @@ +2026-02-10 15:35:51.136388 ,Iter,Sim_time,Wall_time_[s],Mass,RMS_flux,VTK_out,CSV_out,Cur_outflow,Total_outflow,Max_outflow,Total_absorption,Cumulated_absorption_center,Cumulated_absorption_vertical_wall,Cumulated_absorption_horizontal_wall,Var. absorption green,Probe 0 u_0,Probe 0 u_1,Probe 0 u_2,Probe 1 u_0,Probe 1 u_1,Probe 1 u_2,Probe 2 u_0,Probe 2 u_1,Probe 2 u_2,Probe 3 u_0,Probe 3 u_1,Probe 3 u_2,Probe Green Line 0,Probe Green Line 1,Probe Green Line 2,Probe Green Line 3,Probe Green Line 4,Probe Green Line 5,Probe Green Line 6,Probe Green Line 7,Probe Green Line 8,Probe Green Line 9,Probe Green Line 10,Probe Green Line 11,Probe Green Line 12,Probe Green Line 13,Probe Green Line 14,Probe Green Line 15,Probe Green Line 16,Probe Green Line 17,Probe Green Line 18,Probe Green Line 19,Probe Green Block 0,Probe Green Block 1,Probe Green Block 2,Probe Green Block 3,Probe Green Block 4,Probe Green Block 5,Probe Green Block 6,Probe Green Block 7,Probe Green Block 8,Probe Green Block 9,Probe Green Block 10,Probe Green Block 11,Probe Green Block 12,Probe Green Block 13,Probe Green Block 14,Probe Green Block 15,Probe Green Block 16,Probe Green Block 17,Probe Green Block 18,Probe Green Block 19,Probe Green Block 20,Probe Green Block 21,Probe Green Block 22,Probe Green Block 23,Probe Green Block 24,Probe Green Block 25,Probe Green Block 26,Probe Green Block 27,Probe Green Block 28,Probe Green Block 29,Probe Green Block 30,Probe Green Block 31,Probe Green Block 32,Probe Green Block 33,Probe Green Block 34,Probe Green Block 35,Probe Green Block 36,Probe Green Block 37,Probe Green Block 38,Probe Green Block 39,Probe Green Block 40,Probe Green Block 41,Probe Green Block 42,Probe Green Block 43 +2026-02-10 15:35:51.140259 ,0.000000,4.596194E-02,3.629437E-03,3.912731E-01,3.889031E+00,0.000000E+00,1.000000E+00,1.699372E-01,7.810644E-03,2.896596E-01,9.441258E-02,0.000000E+00,0.000000E+00,9.441258E-02,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,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,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,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,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,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,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,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:51.141164 ,1.000000,9.192388E-02,4.564381E-03,6.199409E-01,3.864600E+00,0.000000E+00,1.000000E+00,4.601696E-01,2.896093E-02,8.046266E-01,2.878594E-01,0.000000E+00,0.000000E+00,2.878594E-01,0.000000E+00,3.055669E+00,2.371090E+00,-1.322582E-01,1.120635E+00,-9.104930E-01,2.929389E-01,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,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,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,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,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,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,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:51.142161 ,2.000000,1.378858E-01,5.536966E-03,8.721190E-01,4.519524E+00,0.000000E+00,1.000000E+00,8.639430E-01,6.866942E-02,1.565070E+00,6.283514E-01,0.000000E+00,0.000000E+00,6.283514E-01,0.000000E+00,9.563148E+00,7.246645E+00,-5.433469E-02,4.386915E+00,-3.513465E+00,5.851197E-01,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,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,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,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,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,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,1.612150E-02,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,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00 +2026-02-10 15:35:51.143099 ,3.000000,1.838478E-01,6.507687E-03,1.182900E+00,6.059086E+00,0.000000E+00,1.000000E+00,1.475481E+00,1.364854E-01,2.641718E+00,1.204017E+00,8.438287E-05,0.000000E+00,1.203932E+00,4.895838E-07,1.719370E+01,1.271281E+01,1.233221E-01,9.616643E+00,-7.516527E+00,7.126574E-01,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,1.454970E-02,1.454970E-02,2.052099E-02,1.723574E-02,1.723630E-02,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,1.454970E-02,1.454970E-02,2.052099E-02,1.723574E-02,1.723630E-02,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,1.229449E-04,0.000000E+00,1.734024E-04,0.000000E+00,0.000000E+00,1.456420E-04,0.000000E+00,1.456468E-04,0.000000E+00,0.000000E+00,1.456467E-04,0.000000E+00,0.000000E+00,1.454678E-04,3.315284E-06,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,6.069284E-02,1.229449E-04,0.000000E+00,1.734024E-04,0.000000E+00,0.000000E+00,1.456420E-04,0.000000E+00,1.456468E-04,0.000000E+00,0.000000E+00,1.456467E-04,0.000000E+00,0.000000E+00,1.456341E-04 +2026-02-10 15:35:51.146269 ,4.000000,2.000000E-01,9.664172E-03,9.150041E-01,3.761988E+00,1.000000E+00,1.000000E+00,9.972836E-01,1.525938E-01,2.641718E+00,1.247551E+00,1.397201E-04,0.000000E+00,1.247411E+00,4.707055E-07,2.004454E+01,1.470103E+01,9.776727E-02,1.146703E+01,-8.877042E+00,6.309063E-01,6.206449E-04,5.566883E-04,-1.933037E-04,1.092471E-03,8.122029E-04,-5.300880E-04,0.000000E+00,1.207213E-02,5.702114E-03,1.204194E-02,1.204194E-02,1.827069E-02,1.504029E-02,1.506527E-02,0.000000E+00,0.000000E+00,0.000000E+00,1.207369E-02,5.702114E-03,1.204194E-02,1.204194E-02,1.827069E-02,1.504029E-02,1.506528E-02,0.000000E+00,0.000000E+00,0.000000E+00,0.000000E+00,2.300753E-04,0.000000E+00,1.746656E-04,0.000000E+00,0.000000E+00,0.000000E+00,1.017544E-04,0.000000E+00,1.543873E-04,0.000000E+00,0.000000E+00,1.270905E-04,0.000000E+00,1.273016E-04,0.000000E+00,0.000000E+00,1.273329E-04,0.000000E+00,0.000000E+00,1.268685E-04,9.557303E-05,0.000000E+00,1.770442E-06,0.000000E+00,1.753167E-04,0.000000E+00,0.000000E+00,2.014592E-02,1.017544E-04,0.000000E+00,1.543873E-04,0.000000E+00,0.000000E+00,1.270905E-04,0.000000E+00,1.273016E-04,0.000000E+00,0.000000E+00,1.273303E-04,0.000000E+00,0.000000E+00,1.304508E-04 diff --git a/tests/test_cases.cpp b/tests/test_cases.cpp index 37329c85..dfd82ddb 100644 --- a/tests/test_cases.cpp +++ b/tests/test_cases.cpp @@ -1,3 +1,8 @@ +#include +#include +#include +#include +#include #include #include #include @@ -7,6 +12,7 @@ #include "common/config.hpp" #include "datagenerator/datageneratorbase.hpp" +#include "solvers/snsolver_hpc.hpp" #include "solvers/solverbase.hpp" using vtkUnstructuredGridReaderSP = vtkSmartPointer; @@ -31,6 +37,20 @@ std::vector readVTKFile( std::string filename ) { return data; } +void requireAndFlushTestLoggers() { + auto log = spdlog::get( "event" ); + auto logCSV = spdlog::get( "tabular" ); + REQUIRE( log ); + REQUIRE( logCSV ); + log->flush(); + logCSV->flush(); +} + +void resetTestLoggers() { + // Config::InitLogger uses spdlog::flush_every; shutdown avoids flusher-thread races between test cases. + spdlog::shutdown(); +} + TEST_CASE( "SN_SOLVER", "[validation_tests]" ) { std::string sn_fileDir = "input/validation_tests/SN_solver/"; SECTION( "checkerboard" ) { @@ -294,6 +314,8 @@ TEST_CASE( "CSD_PN_SOLVER", "[validation_tests]" ) { if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false; } REQUIRE( errorWithinBounds ); + delete solver; + delete config; } SECTION( "starmap validation 2nd order" ) { @@ -315,6 +337,8 @@ TEST_CASE( "CSD_PN_SOLVER", "[validation_tests]" ) { if( std::fabs( test[i] - reference[i] ) > eps ) errorWithinBounds = false; } REQUIRE( errorWithinBounds ); + delete solver; + delete config; } } @@ -371,9 +395,189 @@ void tokenize( std::string const& str, const char delim, std::vector( token[parsedChars] ) ) ) { + parsedChars++; + } + return parsedChars == token.size(); + } + catch( ... ) { + return false; + } +} + +std::string findNewestCSVFile( const std::string& directory, const std::string& prefix ) { + namespace fs = std::filesystem; + + if( !fs::exists( directory ) ) { + return ""; + } + + bool found = false; + fs::file_time_type newestWrite = fs::file_time_type::min(); + fs::path newestPath = ""; + const bool enforcePrefix = !prefix.empty(); + auto selectNewestMatchingFileWithPredicate = [&]( bool usePrefix ) { + bool localFound = false; + fs::file_time_type localNewestWrite = fs::file_time_type::min(); + fs::path localNewestPath = ""; + + for( const auto& entry : fs::directory_iterator( directory ) ) { + if( !entry.is_regular_file() || entry.path().extension() != ".csv" ) { + continue; + } + const std::string fileName = entry.path().filename().string(); + if( usePrefix && fileName.rfind( prefix, 0 ) != 0 ) { + continue; + } + if( !localFound || entry.last_write_time() > localNewestWrite ) { + localFound = true; + localNewestWrite = entry.last_write_time(); + localNewestPath = entry.path(); + } + } + if( localFound ) { + found = true; + newestWrite = localNewestWrite; + newestPath = localNewestPath; + } + }; + + selectNewestMatchingFileWithPredicate( enforcePrefix ); + if( !found ) { + selectNewestMatchingFileWithPredicate( false ); // Fallback to newest CSV file in directory + } + + return found ? newestPath.string() : ""; +} + +void compareHistoryCSV( const std::string& referenceFile, const std::string& testFile, double absTol = 1e-10, double relTol = 1e-8 ) { + std::ifstream referenceStream( referenceFile ); + std::ifstream testStream( testFile ); + REQUIRE( referenceStream.is_open() ); + REQUIRE( testStream.is_open() ); + + std::string lineRef, lineTest; + const char delim = ','; + + // Header row (skip token 0 because it contains runtime timestamp) + REQUIRE( std::getline( referenceStream, lineRef ) ); + REQUIRE( std::getline( testStream, lineTest ) ); + + std::vector headerRef; + std::vector headerTest; + tokenize( lineRef, delim, headerRef ); + tokenize( lineTest, delim, headerTest ); + + REQUIRE( headerRef.size() == headerTest.size() ); + REQUIRE( headerRef.size() > 1 ); + const unsigned expectedTokens = static_cast( headerRef.size() ); + + std::vector skipColumn( headerRef.size(), false ); + for( unsigned idx_token = 1; idx_token < headerRef.size(); idx_token++ ) { + REQUIRE( headerRef[idx_token] == headerTest[idx_token] ); + if( headerRef[idx_token] == "Wall_time_[s]" ) { + skipColumn[idx_token] = true; // wall time is machine/runtime dependent + } + } + + unsigned lineNumber = 1; + while( true ) { + const bool hasRefLine = static_cast( std::getline( referenceStream, lineRef ) ); + const bool hasTestLine = static_cast( std::getline( testStream, lineTest ) ); + + if( !hasRefLine || !hasTestLine ) { + REQUIRE( hasRefLine == hasTestLine ); + break; + } + + lineNumber++; + + std::vector tokensRef; + std::vector tokensTest; + tokenize( lineRef, delim, tokensRef ); + tokenize( lineTest, delim, tokensTest ); + + INFO( "Line " << lineNumber ); + REQUIRE( tokensRef.size() == tokensTest.size() ); + REQUIRE( tokensRef.size() > 1 ); + REQUIRE( tokensRef.size() == expectedTokens ); + + for( unsigned idx_token = 1; idx_token < tokensRef.size(); idx_token++ ) { // Skip date/time prefix in first column + if( skipColumn[idx_token] ) { + continue; + } + + double valueRef = 0.0; + double valueTest = 0.0; + bool refNumeric = parseTokenAsDouble( tokensRef[idx_token], valueRef ); + bool tstNumeric = parseTokenAsDouble( tokensTest[idx_token], valueTest ); + + INFO( "Line " << lineNumber << ", token " << idx_token ); + if( refNumeric && tstNumeric ) { + double scale = std::max( 1.0, std::max( std::fabs( valueRef ), std::fabs( valueTest ) ) ); + REQUIRE( std::fabs( valueRef - valueTest ) <= absTol + relTol * scale ); + } + else { + REQUIRE( tokensRef[idx_token] == tokensTest[idx_token] ); + } + } + } +} + +TEST_CASE( "SN_SOLVER_HPC_CSV_QOI_VALIDATION", "[validation_tests][hpc]" ) { + std::string hpcFileDir = "input/validation_tests/SN_solver_hpc/"; + resetTestLoggers(); // Ensure deterministic logger files for this test case + + SECTION( "lattice_200_cells_all_qois" ) { + std::string configFileName = std::string( TESTS_PATH ) + hpcFileDir + "lattice_hpc_200.cfg"; + std::string referenceCSV = std::string( TESTS_PATH ) + hpcFileDir + "lattice_hpc_200_csv_reference"; + std::string resultDir = std::string( TESTS_PATH ) + "result"; + std::string logDir = std::string( TESTS_PATH ) + "result/logs"; + std::filesystem::remove_all( resultDir ); + + Config* config = new Config( configFileName ); + SNSolverHPC* solver = new SNSolverHPC( config ); + solver->Solve(); + + requireAndFlushTestLoggers(); + + std::string outputCSV = findNewestCSVFile( logDir, "lattice_hpc_200_output" ); + REQUIRE( !outputCSV.empty() ); + compareHistoryCSV( referenceCSV, outputCSV ); + + delete solver; + delete config; + } + + SECTION( "symmetric_hohlraum_200_cells_all_qois" ) { + std::string configFileName = std::string( TESTS_PATH ) + hpcFileDir + "symmetric_hohlraum_hpc_200.cfg"; + std::string referenceCSV = std::string( TESTS_PATH ) + hpcFileDir + "symmetric_hohlraum_hpc_200_csv_reference"; + std::string resultDir = std::string( TESTS_PATH ) + "result"; + std::string logDir = std::string( TESTS_PATH ) + "result/logs"; + std::filesystem::remove_all( resultDir ); + + Config* config = new Config( configFileName ); + SNSolverHPC* solver = new SNSolverHPC( config ); + solver->Solve(); + + requireAndFlushTestLoggers(); + + std::string outputCSV = findNewestCSVFile( logDir, "symmetric_hohlraum_hpc_200_output" ); + REQUIRE( !outputCSV.empty() ); + compareHistoryCSV( referenceCSV, outputCSV ); + + delete solver; + delete config; + } +} + TEST_CASE( "screen_output", "[output]" ) { std::string out_fileDir = "input/validation_tests/output/"; - spdlog::drop_all(); // Make sure to write in own logging file + resetTestLoggers(); // Make sure to write in own logging file std::string config_file_name = std::string( TESTS_PATH ) + out_fileDir + "validate_logger.cfg"; std::string screenLoggerReference = std::string( TESTS_PATH ) + out_fileDir + "validate_logger_reference"; @@ -388,10 +592,7 @@ TEST_CASE( "screen_output", "[output]" ) { solver->Solve(); // Force Logger to flush - auto log = spdlog::get( "event" ); - auto logCSV = spdlog::get( "tabular" ); - log->flush(); - logCSV->flush(); + requireAndFlushTestLoggers(); // --- Read and validate logger --- std::ifstream screenLoggerReferenceStream( screenLoggerReference ); std::ifstream screenLoggerStream( screenLogger ); @@ -458,11 +659,14 @@ TEST_CASE( "screen_output", "[output]" ) { std::cout << "Files of unequal length!\n"; } REQUIRE( eqLen ); // Files must be of same length + + delete solver; + delete config; } TEST_CASE( "Data Generator Regression", "[dataGen]" ) { std::string out_fileDir = "input/validation_tests/dataGenerator/"; - spdlog::drop_all(); // Make sure to write in own logging file + resetTestLoggers(); // Make sure to write in own logging file std::string config_file_name = std::string( TESTS_PATH ) + out_fileDir + "validate_dataGen_regression.cfg"; std::string historyLoggerReference = std::string( TESTS_PATH ) + out_fileDir + "validate_dataGen_regression_csv_reference"; @@ -477,10 +681,7 @@ TEST_CASE( "Data Generator Regression", "[dataGen]" ) { datagen->ComputeTrainingData(); // --- Force Logger to flush - auto log = spdlog::get( "event" ); - auto logCSV = spdlog::get( "tabular" ); - log->flush(); - logCSV->flush(); + requireAndFlushTestLoggers(); // --- Read and validate logger --- std::ifstream historyLoggerReferenceStream( historyLoggerReference ); @@ -521,11 +722,12 @@ TEST_CASE( "Data Generator Regression", "[dataGen]" ) { REQUIRE( testPassed ); delete datagen; + delete config; } TEST_CASE( "Data Generator Classification", "[dataGen]" ) { std::string out_fileDir = "input/validation_tests/dataGenerator/"; - spdlog::drop_all(); // Make sure to write in own logging file + resetTestLoggers(); // Make sure to write in own logging file std::string config_file_name = std::string( TESTS_PATH ) + out_fileDir + "validate_dataGen_classification.cfg"; std::string historyLoggerReference = std::string( TESTS_PATH ) + out_fileDir + "validate_dataGen_csv_classification_reference"; @@ -540,10 +742,7 @@ TEST_CASE( "Data Generator Classification", "[dataGen]" ) { datagen->ComputeTrainingData(); // --- Force Logger to flush - auto log = spdlog::get( "event" ); - auto logCSV = spdlog::get( "tabular" ); - log->flush(); - logCSV->flush(); + requireAndFlushTestLoggers(); // --- Read and validate logger --- std::ifstream historyLoggerReferenceStream( historyLoggerReference ); @@ -589,11 +788,12 @@ TEST_CASE( "Data Generator Classification", "[dataGen]" ) { REQUIRE( testPassed ); delete datagen; + delete config; } TEST_CASE( "Data Generator Regularized Regression", "[dataGen]" ) { std::string out_fileDir = "input/validation_tests/dataGenerator/"; - spdlog::drop_all(); // Make sure to write in own logging file + resetTestLoggers(); // Make sure to write in own logging file std::string config_file_name = std::string( TESTS_PATH ) + out_fileDir + "validate_dataGen_regularized_regression.cfg"; std::string historyLoggerReference = std::string( TESTS_PATH ) + out_fileDir + "validate_dataGen_csv_regularized_reference"; @@ -608,10 +808,7 @@ TEST_CASE( "Data Generator Regularized Regression", "[dataGen]" ) { datagen->ComputeTrainingData(); // --- Force Logger to flush - auto log = spdlog::get( "event" ); - auto logCSV = spdlog::get( "tabular" ); - log->flush(); - logCSV->flush(); + requireAndFlushTestLoggers(); // --- Read and validate logger --- std::ifstream historyLoggerReferenceStream( historyLoggerReference ); @@ -653,4 +850,5 @@ TEST_CASE( "Data Generator Regularized Regression", "[dataGen]" ) { } REQUIRE( testPassed ); delete datagen; + delete config; } diff --git a/tests/test_config.cpp b/tests/test_config.cpp index ec1e7b0a..27d7d8bb 100644 --- a/tests/test_config.cpp +++ b/tests/test_config.cpp @@ -1,11 +1,12 @@ #include "catch.hpp" #include "common/config.hpp" +#include TEST_CASE( "Read in Config Template" ) { std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/common/unit_config.cfg"; // Load Settings from File - Config* config = new Config( filename ); + auto config = std::make_unique( filename ); // Test all set configurations bool allSatisfied = true; diff --git a/tests/test_config_extended.cpp b/tests/test_config_extended.cpp new file mode 100644 index 00000000..63b067df --- /dev/null +++ b/tests/test_config_extended.cpp @@ -0,0 +1,115 @@ +#include "catch.hpp" +#include "common/config.hpp" + +TEST_CASE( "extended config parameter tests", "[config_extended]" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/common/unit_config.cfg"; + Config* config = new Config( filename ); + + SECTION( "solver parameters" ) { + // Default solver is SN_SOLVER + REQUIRE( config->GetSolverName() == SN_SOLVER ); + // Default entropy + REQUIRE( config->GetEntropyName() == QUADRATIC ); + // SigmaS default for linesource + REQUIRE( config->GetSigmaS() >= 0.0 ); + // Max moment degree default + REQUIRE( config->GetMaxMomentDegree() >= 0 ); + } + + SECTION( "output paths are non-empty" ) { + REQUIRE( !config->GetOutputFile().empty() ); + REQUIRE( !config->GetLogDir().empty() ); + } + + SECTION( "mesh file path" ) { REQUIRE( !config->GetMeshFile().empty() ); } + + SECTION( "optimizer defaults" ) { + REQUIRE( config->GetNewtonOptimizerEpsilon() > 0.0 ); + REQUIRE( config->GetNewtonIter() > 0 ); + REQUIRE( config->GetNewtonStepSize() > 0.0 ); + REQUIRE( config->GetNewtonMaxLineSearches() > 0 ); + } + + SECTION( "CFL and time" ) { + REQUIRE( config->GetCFL() == Approx( 0.4 ) ); + REQUIRE( config->GetTEnd() == Approx( 0.3 ) ); + } + + SECTION( "quadrature" ) { + REQUIRE( config->GetQuadName() == QUAD_MonteCarlo ); + REQUIRE( config->GetQuadOrder() == 5000 ); + } + + SECTION( "problem and kernel" ) { + REQUIRE( config->GetProblemName() == PROBLEM_Linesource ); + REQUIRE( config->GetKernelName() == KERNEL_Isotropic ); + } + + SECTION( "boundary conditions" ) { + REQUIRE( config->GetBoundaryType( "DirichletTestMarker1" ) == DIRICHLET ); + REQUIRE( config->GetBoundaryType( "DirichletTestMarker2" ) == DIRICHLET ); + REQUIRE( config->GetBoundaryType( "NeumannTestMarker1" ) == NEUMANN ); + REQUIRE( config->GetBoundaryType( "NeumannTestMarker2" ) == NEUMANN ); + // Unknown marker should return INVALID + REQUIRE( config->GetBoundaryType( "UnknownMarker" ) == INVALID ); + } + + SECTION( "boolean defaults" ) { + // HPC default should be false + REQUIRE( config->GetHPC() == false ); + // CSD default + REQUIRE( config->GetIsCSD() == false ); + } + + delete config; +} + +TEST_CASE( "optimizer config parsing", "[config_extended]" ) { + + SECTION( "Newton optimizer config" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerNewton.cfg"; + Config* config = new Config( filename ); + REQUIRE( config->GetOptimizerName() == NEWTON ); + REQUIRE( config->GetEntropyName() == QUADRATIC ); + REQUIRE( config->GetMaxMomentDegree() == 1 ); + REQUIRE( config->GetNewtonFastMode() == false ); + REQUIRE( config->GetNewtonIter() == 100 ); + REQUIRE( config->GetNewtonOptimizerEpsilon() == Approx( 0.01 ) ); + REQUIRE( config->GetNewtonStepSize() == Approx( 0.5 ) ); + REQUIRE( config->GetNewtonMaxLineSearches() == 100 ); + REQUIRE( config->GetQuadName() == QUAD_GaussLegendreTensorized ); + REQUIRE( config->GetQuadOrder() == 4 ); + REQUIRE( config->GetSNAllGaussPts() == true ); + delete config; + } + + SECTION( "Regularized Newton optimizer config" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerRegularizedNewton.cfg"; + Config* config = new Config( filename ); + REQUIRE( config->GetOptimizerName() == REGULARIZED_NEWTON ); + REQUIRE( config->GetEntropyName() == MAXWELL_BOLTZMANN ); + REQUIRE( config->GetMaxMomentDegree() == 3 ); + REQUIRE( config->GetRegularizerGamma() == Approx( 0.001 ) ); + REQUIRE( config->GetSphericalBasisName() == SPHERICAL_MONOMIALS ); + REQUIRE( config->GetQuadName() == QUAD_GaussLegendreTensorized2D ); + REQUIRE( config->GetQuadOrder() == 20 ); + delete config; + } + + SECTION( "Part Regularized Newton optimizer config" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerPartRegularizedNewton.cfg"; + Config* config = new Config( filename ); + REQUIRE( config->GetOptimizerName() == PART_REGULARIZED_NEWTON ); + REQUIRE( config->GetEntropyName() == MAXWELL_BOLTZMANN ); + REQUIRE( config->GetMaxMomentDegree() == 3 ); + REQUIRE( config->GetRegularizerGamma() == Approx( 0.001 ) ); + delete config; + } + + SECTION( "MB Newton optimizer config" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerNewtonMB.cfg"; + Config* config = new Config( filename ); + REQUIRE( config->GetEntropyName() == MAXWELL_BOLTZMANN ); + delete config; + } +} diff --git a/tests/test_entropy.cpp b/tests/test_entropy.cpp new file mode 100644 index 00000000..ed223d82 --- /dev/null +++ b/tests/test_entropy.cpp @@ -0,0 +1,93 @@ +#include "catch.hpp" +#include "entropies/maxwellboltzmannentropy.hpp" +#include "entropies/quadraticentropy.hpp" + +TEST_CASE( "quadratic entropy", "[entropy]" ) { + QuadraticEntropy eta; + + SECTION( "known values" ) { + REQUIRE( eta.Entropy( 3.0 ) == Approx( 4.5 ) ); + REQUIRE( eta.Entropy( 0.0 ) == Approx( 0.0 ) ); + REQUIRE( eta.Entropy( -2.0 ) == Approx( 2.0 ) ); + + REQUIRE( eta.EntropyPrime( 3.0 ) == Approx( 3.0 ) ); + REQUIRE( eta.EntropyPrime( 0.0 ) == Approx( 0.0 ) ); + + REQUIRE( eta.EntropyDual( 3.0 ) == Approx( 4.5 ) ); + REQUIRE( eta.EntropyDual( 0.0 ) == Approx( 0.0 ) ); + + REQUIRE( eta.EntropyPrimeDual( 3.0 ) == Approx( 3.0 ) ); + + // Hessian of dual is constant 1 for quadratic entropy + REQUIRE( eta.EntropyHessianDual( 0.0 ) == Approx( 1.0 ) ); + REQUIRE( eta.EntropyHessianDual( 100.0 ) == Approx( 1.0 ) ); + REQUIRE( eta.EntropyHessianDual( -5.0 ) == Approx( 1.0 ) ); + } + + SECTION( "identity: eta_prime and eta_prime_dual are inverses" ) { + // For quadratic entropy, eta'(z) = z and eta*'(y) = y (identity) + for( double z = -3.0; z <= 3.0; z += 0.5 ) { + REQUIRE( eta.EntropyPrimeDual( eta.EntropyPrime( z ) ) == Approx( z ) ); + } + } +} + +TEST_CASE( "Maxwell-Boltzmann entropy", "[entropy]" ) { + MaxwellBoltzmannEntropy eta; + + SECTION( "known values" ) { + // eta(1) = 1*log(1) - 1 = -1 + REQUIRE( eta.Entropy( 1.0 ) == Approx( -1.0 ) ); + // eta'(1) = log(1) = 0 + REQUIRE( eta.EntropyPrime( 1.0 ) == Approx( 0.0 ) ); + // eta*(0) = exp(0) = 1 + REQUIRE( eta.EntropyDual( 0.0 ) == Approx( 1.0 ) ); + // eta*'(0) = exp(0) = 1 + REQUIRE( eta.EntropyPrimeDual( 0.0 ) == Approx( 1.0 ) ); + // Hessian of dual at y: exp(y) + REQUIRE( eta.EntropyHessianDual( 0.0 ) == Approx( 1.0 ) ); + REQUIRE( eta.EntropyHessianDual( 1.0 ) == Approx( std::exp( 1.0 ) ) ); + } + + SECTION( "Legendre-Fenchel duality" ) { + // eta(z) + eta*(eta'(z)) = z * eta'(z) + for( double z = 0.5; z <= 5.0; z += 0.5 ) { + double lhs = eta.Entropy( z ) + eta.EntropyDual( eta.EntropyPrime( z ) ); + double rhs = z * eta.EntropyPrime( z ); + REQUIRE( lhs == Approx( rhs ).epsilon( 1e-10 ) ); + } + } + + SECTION( "derivative consistency" ) { + // Numerical derivative of Entropy should match EntropyPrime + double h = 1e-7; + for( double z = 0.5; z <= 5.0; z += 0.5 ) { + double numDeriv = ( eta.Entropy( z + h ) - eta.Entropy( z - h ) ) / ( 2.0 * h ); + REQUIRE( numDeriv == Approx( eta.EntropyPrime( z ) ).epsilon( 1e-5 ) ); + } + } + + SECTION( "convexity of dual" ) { + // Hessian of dual should be positive for all y + for( double y = -5.0; y <= 5.0; y += 0.5 ) { + REQUIRE( eta.EntropyHessianDual( y ) > 0.0 ); + } + } + + SECTION( "eta_prime and eta_prime_dual are inverses" ) { + // eta'(z) = log(z), eta*'(y) = exp(y) + // eta*'(eta'(z)) = exp(log(z)) = z + for( double z = 0.5; z <= 5.0; z += 0.5 ) { + REQUIRE( eta.EntropyPrimeDual( eta.EntropyPrime( z ) ) == Approx( z ).epsilon( 1e-10 ) ); + } + } + + SECTION( "dual values" ) { + // eta*(y) = exp(y) for Maxwell-Boltzmann + for( double y = -3.0; y <= 3.0; y += 0.5 ) { + REQUIRE( eta.EntropyDual( y ) == Approx( std::exp( y ) ) ); + REQUIRE( eta.EntropyPrimeDual( y ) == Approx( std::exp( y ) ) ); + REQUIRE( eta.EntropyHessianDual( y ) == Approx( std::exp( y ) ) ); + } + } +} diff --git a/tests/test_io.cpp b/tests/test_io.cpp new file mode 100644 index 00000000..c83a8a63 --- /dev/null +++ b/tests/test_io.cpp @@ -0,0 +1,151 @@ +#include +#include + +#include "catch.hpp" +#include "common/io.hpp" +#include "common/typedef.hpp" + +TEST_CASE( "IO connectivity round-trip", "[io]" ) { + std::string testDir = std::string( TESTS_PATH ) + "result/"; + std::string testFile = testDir + "test_connectivity.csv"; + + // Ensure output directory exists + std::filesystem::create_directories( testDir ); + + unsigned nCells = 3; + unsigned nNodesPerCell = 3; + unsigned nDim = 2; + + // Create test data + std::vector> cellNeighbors( nCells ); + std::vector> cellInterfaceMidPoints( nCells ); + std::vector> cellNormals( nCells ); + std::vector cellBoundaryTypes( nCells ); + + // Fill with known values + for( unsigned i = 0; i < nCells; ++i ) { + cellNeighbors[i].resize( nNodesPerCell ); + cellInterfaceMidPoints[i].resize( nNodesPerCell ); + cellNormals[i].resize( nNodesPerCell ); + for( unsigned j = 0; j < nNodesPerCell; ++j ) { + cellNeighbors[i][j] = i * nNodesPerCell + j; + cellInterfaceMidPoints[i][j] = Vector( nDim, 0.0 ); + cellNormals[i][j] = Vector( nDim, 0.0 ); + for( unsigned k = 0; k < nDim; ++k ) { + cellInterfaceMidPoints[i][j][k] = ( i + 1 ) * 0.1 + ( j + 1 ) * 0.01 + ( k + 1 ) * 0.001; + cellNormals[i][j][k] = ( i + 1 ) * 1.0 + ( j + 1 ) * 0.1 + ( k + 1 ) * 0.01; + } + } + } + cellBoundaryTypes[0] = DIRICHLET; + cellBoundaryTypes[1] = NEUMANN; + cellBoundaryTypes[2] = NONE; + + // Write + WriteConnecitivityToFile( testFile, cellNeighbors, cellInterfaceMidPoints, cellNormals, cellBoundaryTypes, nCells, nDim ); + + // Read back + std::vector> readNeighbors( nCells ); + std::vector> readMidPoints( nCells ); + std::vector> readNormals( nCells ); + std::vector readBoundaryTypes( nCells ); + + LoadConnectivityFromFile( testFile, readNeighbors, readMidPoints, readNormals, readBoundaryTypes, nCells, nNodesPerCell, nDim ); + + SECTION( "neighbors match" ) { + for( unsigned i = 0; i < nCells; ++i ) { + REQUIRE( readNeighbors[i].size() == cellNeighbors[i].size() ); + for( unsigned j = 0; j < nNodesPerCell; ++j ) { + REQUIRE( readNeighbors[i][j] == cellNeighbors[i][j] ); + } + } + } + + SECTION( "midpoints match" ) { + for( unsigned i = 0; i < nCells; ++i ) { + REQUIRE( readMidPoints[i].size() == cellInterfaceMidPoints[i].size() ); + for( unsigned j = 0; j < nNodesPerCell; ++j ) { + for( unsigned k = 0; k < nDim; ++k ) { + REQUIRE( readMidPoints[i][j][k] == Approx( cellInterfaceMidPoints[i][j][k] ).epsilon( 1e-12 ) ); + } + } + } + } + + SECTION( "normals match" ) { + for( unsigned i = 0; i < nCells; ++i ) { + REQUIRE( readNormals[i].size() == cellNormals[i].size() ); + for( unsigned j = 0; j < nNodesPerCell; ++j ) { + for( unsigned k = 0; k < nDim; ++k ) { + REQUIRE( readNormals[i][j][k] == Approx( cellNormals[i][j][k] ).epsilon( 1e-12 ) ); + } + } + } + } + + SECTION( "boundary types match" ) { + REQUIRE( readBoundaryTypes[0] == DIRICHLET ); + REQUIRE( readBoundaryTypes[1] == NEUMANN ); + REQUIRE( readBoundaryTypes[2] == NONE ); + } + + // Cleanup + std::filesystem::remove( testFile ); +} + +TEST_CASE( "IO restart solution round-trip", "[io]" ) { + std::string testDir = std::string( TESTS_PATH ) + "result/"; + std::string baseFile = testDir + "test_restart"; + + // Ensure output directory exists + std::filesystem::create_directories( testDir ); + + // Create test data + int rank = 0; + int idx_iter = 42; + double totalAbsorptionCenter = 1.23456789; + double totalAbsorptionVertical = 2.34567890; + double totalAbsorptionHorizontal = 3.45678901; + double totalAbsorption = 4.56789012; + unsigned long nCells = 5; + std::vector solution = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 }; + std::vector scalarFlux = { 10.1, 10.2, 10.3, 10.4, 10.5 }; + + // Write + WriteRestartSolution( baseFile, solution, scalarFlux, rank, idx_iter, totalAbsorptionCenter, totalAbsorptionVertical, + totalAbsorptionHorizontal, totalAbsorption ); + + // Read back + std::vector readSolution; + std::vector readScalarFlux; + double readCenter = 0, readVertical = 0, readHorizontal = 0, readTotal = 0; + + int readIter = LoadRestartSolution( baseFile, readSolution, readScalarFlux, rank, nCells, readCenter, readVertical, readHorizontal, readTotal ); + + SECTION( "iteration number matches" ) { REQUIRE( readIter == idx_iter ); } + + SECTION( "absorption values match" ) { + REQUIRE( readCenter == Approx( totalAbsorptionCenter ) ); + REQUIRE( readVertical == Approx( totalAbsorptionVertical ) ); + REQUIRE( readHorizontal == Approx( totalAbsorptionHorizontal ) ); + REQUIRE( readTotal == Approx( totalAbsorption ) ); + } + + SECTION( "solution data matches" ) { + REQUIRE( readSolution.size() == solution.size() ); + for( unsigned i = 0; i < solution.size(); ++i ) { + REQUIRE( readSolution[i] == Approx( solution[i] ) ); + } + } + + SECTION( "scalar flux data matches" ) { + REQUIRE( readScalarFlux.size() == nCells ); + for( unsigned i = 0; i < nCells; ++i ) { + REQUIRE( readScalarFlux[i] == Approx( scalarFlux[i] ) ); + } + } + + // Cleanup + std::string restartFile = baseFile + "_restart_rank_0"; + std::filesystem::remove( restartFile ); +} diff --git a/tests/test_kernel.cpp b/tests/test_kernel.cpp index 0785aa59..c78e9280 100644 --- a/tests/test_kernel.cpp +++ b/tests/test_kernel.cpp @@ -1,4 +1,5 @@ #include +#include #include "catch.hpp" #include "common/config.hpp" @@ -9,14 +10,13 @@ TEST_CASE( "test all scattering kernels", "[kernel]" ) { std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/kernels/unit_kernel.cfg"; // Load Settings from File - Config* config = new Config( filename ); - - QuadratureBase* quad = QuadratureBase::Create( config ); //@TODO: swap out for different quadrature rule + auto config = std::make_unique( filename ); + auto quad = std::unique_ptr( QuadratureBase::Create( config.get() ) ); //@TODO: swap out for different quadrature rule SECTION( "isotropic scattering kernel" ) { auto weights = quad->GetWeights(); - Isotropic kernel( quad ); + Isotropic kernel( quad.get() ); Matrix scatteringMatrix = kernel.GetScatteringKernel(); bool errorWithinBounds = true; for( unsigned i = 0; i < scatteringMatrix.rows(); ++i ) { diff --git a/tests/test_mesh.cpp b/tests/test_mesh.cpp index b561bfc1..e6c834f7 100644 --- a/tests/test_mesh.cpp +++ b/tests/test_mesh.cpp @@ -7,13 +7,14 @@ #include "common/mesh.hpp" #include #include +#include TEST_CASE( "unit mesh tests", "[mesh]" ) { std::string config_file_name = std::string( TESTS_PATH ) + "input/unit_tests/common/unit_mesh.cfg"; - Config* config = new Config( config_file_name ); + auto config = std::make_unique( config_file_name ); config->SetForcedConnectivity( true ); - Mesh* mesh = LoadSU2MeshFromFile( config ); + auto mesh = std::unique_ptr( LoadSU2MeshFromFile( config.get() ) ); SECTION( "sum of all cell areas is equal to total domain volume" ) { double domainArea = 1.0; @@ -112,7 +113,7 @@ TEST_CASE( "unit mesh tests", "[mesh]" ) { } else { config->SetForcedConnectivity( false ); - Mesh* mesh2 = LoadSU2MeshFromFile( config ); + Mesh* mesh2 = LoadSU2MeshFromFile( config.get() ); // Check cell number REQUIRE( mesh2->GetNumCells() == mesh->GetNumCells() ); @@ -165,8 +166,8 @@ TEST_CASE( "reconstruction tests", "[mesh]" ) { SECTION( "ensure correct Gauss theorem" ) { std::string config_file_name = std::string( TESTS_PATH ) + "input/unit_tests/common/unit_mesh.cfg"; - Config* config = new Config( config_file_name ); - Mesh* mesh = LoadSU2MeshFromFile( config ); + auto config = std::make_unique( config_file_name ); + auto mesh = std::unique_ptr( LoadSU2MeshFromFile( config.get() ) ); int numCells = mesh->GetNumCells(); int nq = 5; diff --git a/tests/test_numericalflux_extended.cpp b/tests/test_numericalflux_extended.cpp new file mode 100644 index 00000000..c0ad816c --- /dev/null +++ b/tests/test_numericalflux_extended.cpp @@ -0,0 +1,288 @@ +#include + +#include "catch.hpp" +#include "common/typedef.hpp" +#include "fluxes/upwindflux.hpp" + +TEST_CASE( "upwind flux extended tests", "[numericalflux_extended]" ) { + UpwindFlux g; + + SECTION( "FluxXZ scalar - upwind selection" ) { + // FluxXZ uses Omega[0]*n[0] + Omega[2]*n[1] + Vector omega( 3, 0.0 ); + omega[0] = 1.0; + omega[1] = 0.0; // y component ignored + omega[2] = 0.5; + + Vector n( 2, 0.0 ); + n[0] = 1.0; + n[1] = 1.0; + + double psiL = 2.0; + double psiR = 5.0; + + // inner = 1.0*1.0 + 0.5*1.0 = 1.5 > 0 => use psiL + double result = g.FluxXZ( omega, psiL, psiR, n ); + REQUIRE( result == Approx( 1.5 * psiL ) ); + + // Flip normal to make inner negative + Vector nFlip( 2, 0.0 ); + nFlip[0] = -1.0; + nFlip[1] = -1.0; + double resultFlip = g.FluxXZ( omega, psiL, psiR, nFlip ); + REQUIRE( resultFlip == Approx( -1.5 * psiR ) ); + } + + SECTION( "FluxXZ scalar - symmetry" ) { + // F(omega, psiL, psiR, n) = -F(omega, psiR, psiL, -n) + Vector omega( 3, 0.0 ); + omega[0] = 0.7; + omega[1] = 0.3; // ignored + omega[2] = -0.4; + + Vector n( 2, 0.0 ); + n[0] = 0.5; + n[1] = -0.8; + + double psiL = 3.0; + double psiR = 7.0; + + double fPlus = g.FluxXZ( omega, psiL, psiR, n ); + double fMinus = g.FluxXZ( omega, psiR, psiL, -n ); + REQUIRE( std::fabs( fPlus + fMinus ) < 1e-12 ); + } + + SECTION( "FluxXZ scalar - consistency" ) { + // When psiL == psiR, flux = inner * psi + Vector omega( 3, 0.0 ); + omega[0] = 0.6; + omega[1] = 99.0; // ignored + omega[2] = -0.3; + + Vector n( 2, 0.0 ); + n[0] = 1.2; + n[1] = -0.7; + + double psi = 4.5; + double inner = omega[0] * n[0] + omega[2] * n[1]; + double result = g.FluxXZ( omega, psi, psi, n ); + REQUIRE( result == Approx( inner * psi ) ); + } + + SECTION( "Flux1D scalar - upwind selection" ) { + // Flux1D uses only Omega[0]*n[0] + Vector omega( 1, 0.0 ); + omega[0] = 2.0; + + Vector n( 1, 0.0 ); + n[0] = 1.0; + + double psiL = 3.0; + double psiR = 7.0; + + // inner = 2.0 > 0 => use psiL + double result = g.Flux1D( omega, psiL, psiR, n ); + REQUIRE( result == Approx( 2.0 * psiL ) ); + + // Negative direction + n[0] = -1.0; + result = g.Flux1D( omega, psiL, psiR, n ); + // inner = -2.0 < 0 => use psiR + REQUIRE( result == Approx( -2.0 * psiR ) ); + } + + SECTION( "Flux1D scalar - symmetry" ) { + Vector omega( 1, 0.0 ); + omega[0] = -1.5; + Vector n( 1, 0.0 ); + n[0] = 0.8; + + double psiL = 2.0; + double psiR = 5.0; + + double fPlus = g.Flux1D( omega, psiL, psiR, n ); + double fMinus = g.Flux1D( omega, psiR, psiL, -n ); + REQUIRE( std::fabs( fPlus + fMinus ) < 1e-12 ); + } + + SECTION( "Flux1D matrix - Steger-Warming" ) { + // 2x2 system: Flux1D uses only x-direction + Matrix AxP( 2, 2, 0.0 ); + Matrix AxM( 2, 2, 0.0 ); + AxP( 0, 0 ) = 1.0; + AxM( 1, 1 ) = -1.0; + + Vector psiL{ 3.0, 4.0 }; + Vector psiR{ 1.0, 2.0 }; + + // Positive normal + Vector nPos( 1, 0.0 ); + nPos[0] = 1.0; + Vector resultPos = g.Flux1D( AxP, AxM, psiL, psiR, nPos ); + // n[0]>0: result = n[0]*(AxP*psiL + AxM*psiR) + Vector expectedPos = 1.0 * ( AxP * psiL + AxM * psiR ); + REQUIRE( blaze::norm( resultPos - expectedPos ) < 1e-12 ); + + // Negative normal + Vector nNeg( 1, 0.0 ); + nNeg[0] = -1.0; + Vector resultNeg = g.Flux1D( AxP, AxM, psiL, psiR, nNeg ); + // n[0]<0: result = n[0]*(AxP*psiR + AxM*psiL) + Vector expectedNeg = -1.0 * ( AxP * psiR + AxM * psiL ); + REQUIRE( blaze::norm( resultNeg - expectedNeg ) < 1e-12 ); + } + + SECTION( "FluxXZ matrix - Steger-Warming" ) { + // FluxXZ matrix uses AxPlus/AxMinus for x and AzPlus/AzMinus for y-normal + Matrix AxP( 2, 2, 0.0 ); + Matrix AxM( 2, 2, 0.0 ); + Matrix AyP( 2, 2, 0.0 ); // ignored + Matrix AyM( 2, 2, 0.0 ); // ignored + Matrix AzP( 2, 2, 0.0 ); + Matrix AzM( 2, 2, 0.0 ); + + AxP( 0, 0 ) = 1.0; + AxM( 1, 1 ) = -0.5; + AzP( 0, 1 ) = 0.3; + AzM( 1, 0 ) = -0.2; + + Vector psiL{ 2.0, 3.0 }; + Vector psiR{ 5.0, 6.0 }; + + // Both normals positive + Vector n( 2, 0.0 ); + n[0] = 0.5; + n[1] = 0.8; + + Vector result = g.FluxXZ( AxP, AxM, AyP, AyM, AzP, AzM, psiL, psiR, n ); + + // Expected: n[0]*(AxP*psiL + AxM*psiR) + n[1]*(AzP*psiL + AzM*psiR) + Vector expected = n[0] * ( AxP * psiL + AxM * psiR ) + n[1] * ( AzP * psiL + AzM * psiR ); + REQUIRE( blaze::norm( result - expected ) < 1e-12 ); + } + + SECTION( "FluxXZ matrix - symmetry" ) { + Matrix AxP( 2, 2, 0.0 ); + Matrix AxM( 2, 2, 0.0 ); + Matrix AyP( 2, 2, 0.0 ); + Matrix AyM( 2, 2, 0.0 ); + Matrix AzP( 2, 2, 0.0 ); + Matrix AzM( 2, 2, 0.0 ); + + AxP( 0, 0 ) = 1.0; + AxM( 1, 1 ) = -1.0; + AzP( 0, 0 ) = 0.5; + AzM( 1, 1 ) = -0.5; + + Vector psiL{ 1.0, 2.0 }; + Vector psiR{ 3.0, 4.0 }; + Vector n( 2, 0.0 ); + n[0] = 0.7; + n[1] = -0.3; + + Vector fPlus = g.FluxXZ( AxP, AxM, AyP, AyM, AzP, AzM, psiL, psiR, n ); + Vector fMinus = g.FluxXZ( AxP, AxM, AyP, AyM, AzP, AzM, psiR, psiL, -n ); + REQUIRE( blaze::norm( fPlus + fMinus ) < 1e-12 ); + } + + SECTION( "FluxVanLeer" ) { + Matrix Ax( 2, 2, 0.0 ); + Matrix AxAbs( 2, 2, 0.0 ); + Matrix Ay( 2, 2, 0.0 ); // unused + Matrix AyAbs( 2, 2, 0.0 ); // unused + Matrix Az( 2, 2, 0.0 ); + Matrix AzAbs( 2, 2, 0.0 ); + + Ax( 0, 0 ) = 1.0; + Ax( 1, 1 ) = -0.5; + AxAbs( 0, 0 ) = 1.0; + AxAbs( 1, 1 ) = 0.5; + Az( 0, 0 ) = 0.3; + Az( 1, 1 ) = -0.2; + AzAbs( 0, 0 ) = 0.3; + AzAbs( 1, 1 ) = 0.2; + + Vector psiL{ 2.0, 3.0 }; + Vector psiR{ 5.0, 6.0 }; + Vector n( 2, 0.0 ); + n[0] = 0.7; + n[1] = 0.4; + + Vector resultFlux( 2, 0.0 ); + g.FluxVanLeer( Ax, AxAbs, Ay, AyAbs, Az, AzAbs, psiL, psiR, n, resultFlux ); + + // Expected: 0.5 * ( n[0] * (Ax*(psiL+psiR) - AxAbs*(psiR-psiL)) + n[1] * (Az*(psiL+psiR) - AzAbs*(psiR-psiL)) ) + Vector expected = 0.5 * ( n[0] * ( Ax * ( psiL + psiR ) - AxAbs * ( psiR - psiL ) ) + + n[1] * ( Az * ( psiL + psiR ) - AzAbs * ( psiR - psiL ) ) ); + REQUIRE( blaze::norm( resultFlux - expected ) < 1e-12 ); + } + + SECTION( "FluxVanLeer - consistency" ) { + // When psiL == psiR, difference term vanishes: flux = 0.5 * (n[0]*Ax + n[1]*Az) * 2*psi = (n[0]*Ax + n[1]*Az)*psi + Matrix Ax( 2, 2, 0.0 ); + Matrix AxAbs( 2, 2, 0.0 ); + Matrix Ay( 2, 2, 0.0 ); + Matrix AyAbs( 2, 2, 0.0 ); + Matrix Az( 2, 2, 0.0 ); + Matrix AzAbs( 2, 2, 0.0 ); + + Ax( 0, 0 ) = 2.0; + AxAbs( 0, 0 ) = 2.0; + Az( 1, 1 ) = 1.5; + AzAbs( 1, 1 ) = 1.5; + + Vector psi{ 3.0, 4.0 }; + Vector n( 2, 0.0 ); + n[0] = 0.6; + n[1] = 0.9; + + Vector resultFlux( 2, 0.0 ); + g.FluxVanLeer( Ax, AxAbs, Ay, AyAbs, Az, AzAbs, psi, psi, n, resultFlux ); + + Vector expected = ( n[0] * Ax + n[1] * Az ) * psi; + REQUIRE( blaze::norm( resultFlux - expected ) < 1e-12 ); + } + + SECTION( "vectorized Flux" ) { + // Test the vectorized Flux that accumulates over quadrature points + unsigned nq = 3; + VectorVector quadPts( nq ); + quadPts[0] = Vector{ 1.0, 0.0 }; + quadPts[1] = Vector{ 0.0, 1.0 }; + quadPts[2] = Vector{ -1.0, 0.0 }; + + Vector psiL{ 2.0, 3.0, 4.0 }; + Vector psiR{ 5.0, 6.0, 7.0 }; + Vector flux( nq, 0.0 ); + Vector n{ 1.0, 0.0 }; // normal in x direction + + g.Flux( quadPts, psiL, psiR, flux, n, nq ); + + // Point 0: inner = 1*1 + 0*0 = 1 > 0 => flux[0] += 1*psiL[0] = 2 + // Point 1: inner = 0*1 + 1*0 = 0 <= 0 => flux[1] += 0*psiR[1] = 0 + // Point 2: inner = -1*1 + 0*0 = -1 < 0 => flux[2] += -1*psiR[2] = -7 + REQUIRE( flux[0] == Approx( 2.0 ) ); + REQUIRE( flux[1] == Approx( 0.0 ) ); + REQUIRE( flux[2] == Approx( -7.0 ) ); + } + + SECTION( "vectorized Flux - accumulation" ) { + // Verify flux accumulates (+=) rather than assigns + unsigned nq = 2; + VectorVector quadPts( nq ); + quadPts[0] = Vector{ 1.0, 0.0 }; + quadPts[1] = Vector{ -1.0, 0.0 }; + + Vector psiL{ 1.0, 1.0 }; + Vector psiR{ 1.0, 1.0 }; + Vector flux{ 10.0, 20.0 }; // pre-existing values + Vector n{ 1.0, 0.0 }; + + g.Flux( quadPts, psiL, psiR, flux, n, nq ); + + // Point 0: inner=1 > 0, flux[0] += 1*1 = 11 + // Point 1: inner=-1 < 0, flux[1] += -1*1 = 19 + REQUIRE( flux[0] == Approx( 11.0 ) ); + REQUIRE( flux[1] == Approx( 19.0 ) ); + } +} diff --git a/tests/test_optimizer_extended.cpp b/tests/test_optimizer_extended.cpp new file mode 100644 index 00000000..8df07b64 --- /dev/null +++ b/tests/test_optimizer_extended.cpp @@ -0,0 +1,156 @@ +#include + +#include "catch.hpp" +#include "common/config.hpp" +#include "optimizers/optimizerbase.hpp" +#include "quadratures/quadraturebase.hpp" +#include "velocitybasis/sphericalbase.hpp" +#include "velocitybasis/sphericalharmonics.hpp" + +// Helper to set up moment basis from config +static VectorVector SetupMomentBasis( Config* config, SphericalBase* basis, QuadratureBase* quad ) { + VectorVector momentBasis( quad->GetNq() ); + VectorVector quadPointsSphere = quad->GetPointsSphere(); + for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) { + double my = quadPointsSphere[idx_quad][0]; + double phi = quadPointsSphere[idx_quad][1]; + momentBasis[idx_quad] = basis->ComputeSphericalBasis( my, phi ); + } + return momentBasis; +} + +TEST_CASE( "Newton optimizer with quadratic entropy - zero multipliers", "[optimizers_extended]" ) { + // Test that solving for the isotropic distribution gives near-zero multipliers + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerNewton.cfg"; + Config* config = new Config( filename ); + SphericalHarmonics basis( config->GetMaxMomentDegree() ); + QuadratureBase* quad = QuadratureBase::Create( config ); + OptimizerBase* optimizer = OptimizerBase::Create( config ); + + unsigned nTotalEntries = basis.GetGlobalIndexBasis( config->GetMaxMomentDegree(), config->GetMaxMomentDegree() ) + 1; + + // For quadratic entropy, alpha = u at optimum + // Test with u close to zero (except first component) + Vector u( nTotalEntries, 0.0 ); + u[0] = 1.0; // first moment = 1 + Vector alpha( nTotalEntries, 0.0 ); + + VectorVector moments = VectorVector( quad->GetNq() ); + VectorVector quadPointsSphere = quad->GetPointsSphere(); + for( unsigned idx_quad = 0; idx_quad < quad->GetNq(); idx_quad++ ) { + double my = quadPointsSphere[idx_quad][0]; + double phi = quadPointsSphere[idx_quad][1]; + moments[idx_quad] = basis.ComputeSphericalBasis( my, phi ); + } + + optimizer->Solve( alpha, u, moments ); + + REQUIRE( std::fabs( norm( alpha - u ) ) < config->GetNewtonOptimizerEpsilon() ); + + // Reconstruct and verify + Vector uRecons( nTotalEntries, 0.0 ); + optimizer->ReconstructMoments( uRecons, alpha, moments ); + REQUIRE( std::fabs( norm( uRecons - u ) ) < config->GetNewtonOptimizerEpsilon() ); + + delete optimizer; + delete quad; + delete config; +} + +TEST_CASE( "Newton optimizer with MB entropy - different initial guesses", "[optimizers_extended]" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerNewtonMB.cfg"; + Config* config = new Config( filename ); + SphericalBase* basis = SphericalBase::Create( config ); + QuadratureBase* quad = QuadratureBase::Create( config ); + OptimizerBase* optimizer = OptimizerBase::Create( config ); + + unsigned nTotalEntries = basis->GetBasisSize(); + Vector u( nTotalEntries, 0.5 ); + u[0] = 1.0; + + VectorVector momentBasis = SetupMomentBasis( config, basis, quad ); + + SECTION( "zero initial guess" ) { + Vector alpha( nTotalEntries, 0.0 ); + optimizer->Solve( alpha, u, momentBasis ); + Vector uRecons( nTotalEntries, 0.0 ); + optimizer->ReconstructMoments( uRecons, alpha, momentBasis ); + REQUIRE( std::fabs( norm( uRecons - u ) ) < config->GetNewtonOptimizerEpsilon() ); + } + + SECTION( "non-zero initial guess" ) { + Vector alpha( nTotalEntries, 0.1 ); + optimizer->Solve( alpha, u, momentBasis ); + Vector uRecons( nTotalEntries, 0.0 ); + optimizer->ReconstructMoments( uRecons, alpha, momentBasis ); + REQUIRE( std::fabs( norm( uRecons - u ) ) < config->GetNewtonOptimizerEpsilon() ); + } + + delete optimizer; + delete quad; + delete basis; + delete config; +} + +TEST_CASE( "Regularized Newton - moment reconstruction consistency", "[optimizers_extended]" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerRegularizedNewton.cfg"; + Config* config = new Config( filename ); + SphericalBase* basis = SphericalBase::Create( config ); + QuadratureBase* quad = QuadratureBase::Create( config ); + OptimizerBase* optimizer = OptimizerBase::Create( config ); + + unsigned nTotalEntries = basis->GetBasisSize(); + VectorVector momentBasis = SetupMomentBasis( config, basis, quad ); + + SECTION( "small perturbation from isotropic" ) { + Vector u( nTotalEntries, 0.1 ); + u[0] = 1.0; + Vector alpha( nTotalEntries, 0.0 ); + optimizer->Solve( alpha, u, momentBasis ); + Vector uRecons( nTotalEntries, 0.0 ); + optimizer->ReconstructMoments( uRecons, alpha, momentBasis ); + REQUIRE( std::fabs( norm( uRecons - u ) ) < config->GetNewtonOptimizerEpsilon() ); + } + + SECTION( "isotropic moment vector" ) { + // u = (1, 0, 0, ...) should yield alpha close to 0 + Vector u( nTotalEntries, 0.0 ); + u[0] = 1.0; + Vector alpha( nTotalEntries, 0.0 ); + optimizer->Solve( alpha, u, momentBasis ); + Vector uRecons( nTotalEntries, 0.0 ); + optimizer->ReconstructMoments( uRecons, alpha, momentBasis ); + REQUIRE( std::fabs( norm( uRecons - u ) ) < config->GetNewtonOptimizerEpsilon() ); + } + + delete optimizer; + delete quad; + delete basis; + delete config; +} + +TEST_CASE( "Part Regularized Newton - moment reconstruction", "[optimizers_extended]" ) { + std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/optimizers/unit_optimizerPartRegularizedNewton.cfg"; + Config* config = new Config( filename ); + SphericalBase* basis = SphericalBase::Create( config ); + QuadratureBase* quad = QuadratureBase::Create( config ); + OptimizerBase* optimizer = OptimizerBase::Create( config ); + + unsigned nTotalEntries = basis->GetBasisSize(); + VectorVector momentBasis = SetupMomentBasis( config, basis, quad ); + + // Test with different moment vectors + Vector u( nTotalEntries, 0.3 ); + u[0] = 1.0; + Vector alpha( nTotalEntries, 0.0 ); + optimizer->Solve( alpha, u, momentBasis ); + + Vector uRecons( nTotalEntries, 0.0 ); + optimizer->ReconstructMoments( uRecons, alpha, momentBasis ); + REQUIRE( std::fabs( norm( uRecons - u ) ) < config->GetNewtonOptimizerEpsilon() ); + + delete optimizer; + delete quad; + delete basis; + delete config; +} diff --git a/tests/test_physics.cpp b/tests/test_physics.cpp index 1d4d07fa..1f342f95 100644 --- a/tests/test_physics.cpp +++ b/tests/test_physics.cpp @@ -4,6 +4,7 @@ #include "quadratures/qgausslegendre1D.hpp" #include "toolboxes/epics.hpp" #include "toolboxes/icru.hpp" +#include TEST_CASE( "ENDL: checking angular integral of scattering cross sections to be unity", "[physics]" ) { @@ -92,8 +93,8 @@ TEST_CASE( "Test ICRU database" ) { Vector mu( quadPoints.size() ); for( unsigned i = 0; i < quadPoints.size(); ++i ) mu[i] = quadPoints[i][0]; Vector e( 1, 1 ); // 1 MeV - Config* settings = new Config( std::string( TESTS_PATH ) + "input/unit_tests/problems/icru.cfg" ); - ICRU db( mu, e, settings ); + auto settings = std::make_unique( std::string( TESTS_PATH ) + "input/unit_tests/problems/icru.cfg" ); + ICRU db( mu, e, settings.get() ); SECTION( "angular xs" ) { double totalXS_ref = 3044.818373; @@ -143,6 +144,4 @@ TEST_CASE( "Test ICRU database" ) { REQUIRE( errorWithinBounds ); } - // TODO: causes segmentation fault for unknown reasons - // delete settings; } diff --git a/tests/test_problemgeometry.cpp b/tests/test_problemgeometry.cpp new file mode 100644 index 00000000..3470f1b5 --- /dev/null +++ b/tests/test_problemgeometry.cpp @@ -0,0 +1,219 @@ +#include "catch.hpp" +#include "common/config.hpp" +#include "common/io.hpp" +#include "common/mesh.hpp" +#include "problems/checkerboard.hpp" +#include "problems/lattice.hpp" +#include "problems/linesource.hpp" +#include "quadratures/quadraturebase.hpp" + +TEST_CASE( "checkerboard problem cross sections", "[problem_geometry]" ) { + std::string config_file = std::string( TESTS_PATH ) + "input/validation_tests/SN_solver/checkerboard_SN.cfg"; + Config* config = new Config( config_file ); + Mesh* mesh = LoadSU2MeshFromFile( config ); + QuadratureBase* quad = QuadratureBase::Create( config ); + + Checkerboard_SN problem( config, mesh, quad ); + Vector energies( 1, 1.0 ); + + VectorVector sigmaS = problem.GetScatteringXS( energies ); + VectorVector sigmaT = problem.GetTotalXS( energies ); + + SECTION( "cross section dimensions" ) { + REQUIRE( sigmaS.size() == 1u ); + REQUIRE( sigmaT.size() == 1u ); + REQUIRE( sigmaS[0].size() == mesh->GetNumCells() ); + REQUIRE( sigmaT[0].size() == mesh->GetNumCells() ); + } + + SECTION( "sigma_s <= sigma_t everywhere" ) { + bool valid = true; + for( unsigned i = 0; i < sigmaS[0].size(); ++i ) { + if( sigmaS[0][i] > sigmaT[0][i] + 1e-12 ) valid = false; + } + REQUIRE( valid ); + } + + SECTION( "cross section values are non-negative" ) { + bool valid = true; + for( unsigned i = 0; i < sigmaS[0].size(); ++i ) { + if( sigmaS[0][i] < -1e-12 || sigmaT[0][i] < -1e-12 ) valid = false; + } + REQUIRE( valid ); + } + + SECTION( "source is zero outside source region" ) { + auto externalSource = problem.GetExternalSource( energies ); + auto cellMids = mesh->GetCellMidPoints(); + bool valid = true; + for( unsigned j = 0; j < cellMids.size(); ++j ) { + double x = cellMids[j][0]; + double y = cellMids[j][1]; + // Source should only be non-zero in [3,4] x [3,4] + if( ( x < 3.0 || x > 4.0 || y < 3.0 || y > 4.0 ) && externalSource[0][j][0] > 1e-12 ) { + valid = false; + } + } + REQUIRE( valid ); + } + + SECTION( "initial condition is uniform small value" ) { + VectorVector ic = problem.SetupIC(); + REQUIRE( ic.size() == mesh->GetNumCells() ); + bool valid = true; + for( unsigned j = 0; j < ic.size(); ++j ) { + // All initial values should be 1e-10 + for( unsigned k = 0; k < ic[j].size(); ++k ) { + if( std::fabs( ic[j][k] - 1e-10 ) > 1e-12 ) valid = false; + } + } + REQUIRE( valid ); + } + + delete quad; + delete mesh; + delete config; +} + +TEST_CASE( "lattice problem cross sections", "[problem_geometry]" ) { + // Use the checkerboard SN config as a starting point and modify for lattice + // We need a config with PROBLEM = LATTICE to make ProblemBase::Create work, + // but we can construct Lattice_SN directly if we have a mesh in [-3.5, 3.5]^2 + // For now, use the checkerboard config but construct Lattice directly + std::string config_file = std::string( TESTS_PATH ) + "input/validation_tests/SN_solver/checkerboard_SN.cfg"; + Config* config = new Config( config_file ); + Mesh* mesh = LoadSU2MeshFromFile( config ); + QuadratureBase* quad = QuadratureBase::Create( config ); + + // Even though the mesh doesn't match Lattice domain exactly, + // we can test the interface and cross-section properties + Lattice_SN problem( config, mesh, quad ); + Vector energies( 1, 1.0 ); + + VectorVector sigmaS = problem.GetScatteringXS( energies ); + VectorVector sigmaT = problem.GetTotalXS( energies ); + + SECTION( "cross section dimensions" ) { + REQUIRE( sigmaS.size() == 1u ); + REQUIRE( sigmaT.size() == 1u ); + REQUIRE( sigmaS[0].size() == mesh->GetNumCells() ); + REQUIRE( sigmaT[0].size() == mesh->GetNumCells() ); + } + + SECTION( "sigma_s <= sigma_t everywhere" ) { + bool valid = true; + for( unsigned i = 0; i < sigmaS[0].size(); ++i ) { + if( sigmaS[0][i] > sigmaT[0][i] + 1e-12 ) valid = false; + } + REQUIRE( valid ); + } + + SECTION( "cross sections are non-negative" ) { + bool valid = true; + for( unsigned i = 0; i < sigmaS[0].size(); ++i ) { + if( sigmaS[0][i] < -1e-12 || sigmaT[0][i] < -1e-12 ) valid = false; + } + REQUIRE( valid ); + } + + SECTION( "QOI initialization" ) { + REQUIRE( problem.GetCurAbsorptionLattice() == Approx( 0.0 ) ); + REQUIRE( problem.GetTotalAbsorptionLattice() == Approx( 0.0 ) ); + REQUIRE( problem.GetMaxAbsorptionLattice() == Approx( 0.0 ) ); + } + + delete quad; + delete mesh; + delete config; +} + +TEST_CASE( "line source analytical solution", "[problem_geometry]" ) { + std::string config_file = std::string( TESTS_PATH ) + "input/validation_tests/SN_solver/linesource_SN.cfg"; + Config* config = new Config( config_file ); + Mesh* mesh = LoadSU2MeshFromFile( config ); + QuadratureBase* quad = QuadratureBase::Create( config ); + + LineSource_SN problem( config, mesh, quad ); + Vector energies( 1, 1.0 ); + + SECTION( "cross sections are uniform" ) { + VectorVector sigmaS = problem.GetScatteringXS( energies ); + VectorVector sigmaT = problem.GetTotalXS( energies ); + REQUIRE( sigmaS.size() == 1u ); + REQUIRE( sigmaT.size() == 1u ); + + // All cells should have same scattering XS + bool uniformS = true; + bool uniformT = true; + for( unsigned i = 1; i < sigmaS[0].size(); ++i ) { + if( std::fabs( sigmaS[0][i] - sigmaS[0][0] ) > 1e-12 ) uniformS = false; + if( std::fabs( sigmaT[0][i] - sigmaT[0][0] ) > 1e-12 ) uniformT = false; + } + REQUIRE( uniformS ); + REQUIRE( uniformT ); + } + + SECTION( "external source is zero" ) { + auto Q = problem.GetExternalSource( energies ); + bool allZero = true; + for( unsigned j = 0; j < Q[0].size(); ++j ) { + if( std::fabs( Q[0][j][0] ) > 1e-12 ) allZero = false; + } + REQUIRE( allZero ); + } + + SECTION( "analytical solution - causality" ) { + // For sigma_s = 1: solution should be 0 when t < R (before wavefront arrives) + double t = 0.5; + double R = 0.8; // R > t, so solution should be 0 + double x = R; + double y = 0.0; + double sol = problem.GetAnalyticalSolution( x, y, t, 1.0 ); + REQUIRE( sol == Approx( 0.0 ).margin( 1e-12 ) ); + } + + SECTION( "analytical solution - zero at t=0" ) { + double sol = problem.GetAnalyticalSolution( 0.5, 0.5, 0.0, 1.0 ); + REQUIRE( sol == Approx( 0.0 ) ); + } + + SECTION( "analytical solution - positive inside wavefront" ) { + // For sigma_s = 1: solution should be positive when t > R + double t = 1.0; + double x = 0.1; + double y = 0.0; + double sol = problem.GetAnalyticalSolution( x, y, t, 1.0 ); + REQUIRE( sol > 0.0 ); + } + + SECTION( "analytical solution - radial symmetry" ) { + // Solution at (x,y) should equal solution at (y,x) due to rotational symmetry + double t = 1.0; + double x = 0.3; + double y = 0.4; + double sol1 = problem.GetAnalyticalSolution( x, y, t, 1.0 ); + double sol2 = problem.GetAnalyticalSolution( y, x, t, 1.0 ); + REQUIRE( sol1 == Approx( sol2 ).epsilon( 1e-6 ) ); + + // Also test negation symmetry + double sol3 = problem.GetAnalyticalSolution( -x, y, t, 1.0 ); + REQUIRE( sol1 == Approx( sol3 ).epsilon( 1e-6 ) ); + } + + SECTION( "initial condition is smooth Gaussian" ) { + VectorVector ic = problem.SetupIC(); + REQUIRE( ic.size() == mesh->GetNumCells() ); + // Values should be positive + bool allPositive = true; + for( unsigned j = 0; j < ic.size(); ++j ) { + for( unsigned k = 0; k < ic[j].size(); ++k ) { + if( ic[j][k] < 0 ) allPositive = false; + } + } + REQUIRE( allPositive ); + } + + delete quad; + delete mesh; + delete config; +} diff --git a/tests/test_quadrature.cpp b/tests/test_quadrature.cpp index 330e2cb9..6bc7798e 100644 --- a/tests/test_quadrature.cpp +++ b/tests/test_quadrature.cpp @@ -3,6 +3,7 @@ #include "common/globalconstants.hpp" #include "quadratures/quadraturebase.hpp" +#include #include std::vector quadraturenames = { QUAD_MonteCarlo, @@ -95,7 +96,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Set quadOrder config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); if( quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_Midpoint1D || quadraturename == QUAD_Rectangular1D ) { if( !approxequal( Q->SumUpWeights(), 2, lowAccuracyTesting ) ) { @@ -136,16 +137,15 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Special case for Gauss Legendre with half weights if( quadraturename == QUAD_GaussLegendreTensorized ) { config->SetSNAllGaussPts( false ); - QuadratureBase* Q = QuadratureBase::Create( config ); - if( !approxequal( Q->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) ) { + std::unique_ptr QHalf( QuadratureBase::Create( config ) ); + if( !approxequal( QHalf->SumUpWeights(), 4 * M_PI, lowAccuracyTesting ) ) { testPassed = false; - PrintErrorMsg( config, std::abs( Q->SumUpWeights() - 4 * M_PI ), Q->SumUpWeights(), lowAccuracyTesting ); + PrintErrorMsg( config, std::abs( QHalf->SumUpWeights() - 4 * M_PI ), QHalf->SumUpWeights(), lowAccuracyTesting ); printf( "Reduced number of quadrature was points used. \n" ); } config->SetSNAllGaussPts( true ); } - delete Q; } } REQUIRE( testPassed ); @@ -174,7 +174,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Set quadOrder config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); VectorVector points = Q->GetPoints(); for( unsigned i = 0; i < Q->GetNq(); i++ ) { if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) { @@ -187,10 +187,10 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Special case for Gauss Legendre with half weights if( quadraturename == QUAD_GaussLegendreTensorized ) { config->SetSNAllGaussPts( false ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr QHalf( QuadratureBase::Create( config ) ); - VectorVector points = Q->GetPoints(); - for( unsigned i = 0; i < Q->GetNq(); i++ ) { + VectorVector points = QHalf->GetPoints(); + for( unsigned i = 0; i < QHalf->GetNq(); i++ ) { if( !approxequal( 1.0, norm( points[i] ), lowAccuracyTesting ) ) { testPassed = false; PrintErrorMsg( config, std::abs( norm( points[i] ) - 1.0 ), norm( points[i] ), lowAccuracyTesting ); @@ -201,7 +201,6 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { config->SetSNAllGaussPts( true ); } - delete Q; } } REQUIRE( testPassed ); @@ -224,7 +223,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Set quadOrder config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); VectorVector points = Q->GetPoints(); //(v_x,v_y,v_z) VectorVector pointsSphere = Q->GetPointsSphere(); // (mu, phi, r) @@ -255,7 +254,6 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { printf( "Errorous index: %d\n", i ); } } - delete Q; } } REQUIRE( testPassed ); @@ -271,18 +269,17 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Set quadOrder config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); REQUIRE( Q->GetNq() == size( Q->GetWeights() ) ); // Special case for Gauss Legendre with half weights if( quadraturename == QUAD_GaussLegendreTensorized ) { config->SetSNAllGaussPts( false ); - QuadratureBase* Q = QuadratureBase::Create( config ); - REQUIRE( Q->GetNq() == size( Q->GetWeights() ) ); + std::unique_ptr QHalf( QuadratureBase::Create( config ) ); + REQUIRE( QHalf->GetNq() == size( QHalf->GetWeights() ) ); config->SetSNAllGaussPts( true ); } - delete Q; } } } @@ -297,17 +294,16 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Set quadOrder config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); REQUIRE( Q->GetNq() == size( Q->GetPoints() ) ); // Special case for Gauss Legendre with half weights if( quadraturename == QUAD_GaussLegendreTensorized ) { config->SetSNAllGaussPts( false ); - QuadratureBase* Q = QuadratureBase::Create( config ); - REQUIRE( Q->GetNq() == size( Q->GetPoints() ) ); + std::unique_ptr QHalf( QuadratureBase::Create( config ) ); + REQUIRE( QHalf->GetNq() == size( QHalf->GetPoints() ) ); config->SetSNAllGaussPts( true ); } - delete Q; } } } @@ -329,7 +325,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Set quadOrder config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); if( quadraturename == QUAD_Rectangular1D || quadraturename == QUAD_GaussLegendre1D || quadraturename == QUAD_Midpoint1D ) { result = Q->Integrate( sin ); @@ -374,9 +370,9 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Special case for Gauss Legendre with half weights if( quadraturename == QUAD_GaussLegendreTensorized ) { config->SetSNAllGaussPts( false ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr QHalf( QuadratureBase::Create( config ) ); - result = Q->Integrate( f ); + result = QHalf->Integrate( f ); if( !approxequal( result, 4.0 * M_PI, lowAccuracyTesting ) ) { testPassed = false; PrintErrorMsg( config, std::abs( result - 4.0 * M_PI ), result, lowAccuracyTesting ); @@ -404,7 +400,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { config->SetQuadName( quadraturename ); for( auto quadratureorder : quadratureorders[quadraturename] ) { config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); points = Q->GetPoints(); pointsSphere = Q->GetPointsSphere(); for( unsigned idx_nq = 0; idx_nq < Q->GetNq(); idx_nq++ ) { @@ -434,7 +430,6 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { } } } - delete Q; } } } @@ -465,7 +460,7 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { // Set quadOrder config->SetQuadOrder( quadratureorder ); - QuadratureBase* Q = QuadratureBase::Create( config ); + std::unique_ptr Q( QuadratureBase::Create( config ) ); // Note: Leaving out Quad_GaussLegendreTensorized with half weights... (to be added) if( quadraturename != QUAD_SphericalTessalation2D && quadraturename != QUAD_GaussLegendreTensorized2D && @@ -609,8 +604,6 @@ TEST_CASE( "Quadrature Tests", "[quadrature]" ) { printf( "Error at integrating Polygon.\n" ); } } - - delete Q; } } REQUIRE( testPassed ); diff --git a/tests/test_sphericalBasis.cpp b/tests/test_sphericalBasis.cpp index 84aa8664..0a5976c8 100644 --- a/tests/test_sphericalBasis.cpp +++ b/tests/test_sphericalBasis.cpp @@ -6,6 +6,7 @@ #include #include +#include #include double Y0_0( double, double ) { return sqrt( 1 / ( 4 * M_PI ) ); } @@ -33,7 +34,7 @@ TEST_CASE( "test spherical harmonics basis ", "[spherical_harmonics]" ) { std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/solvers/unit_harmonics.cfg"; // Load Settings from File - Config* config = new Config( filename ); + auto config = std::make_unique( filename ); unsigned maxMomentDegree = 2; @@ -127,7 +128,7 @@ TEST_CASE( "test spherical harmonics basis ", "[spherical_harmonics]" ) { SECTION( "test orthonormality - spherical coordinates" ) { // Caution: Integration only works with spherical coordinates! - QGaussLegendreTensorized quad( config ); + QGaussLegendreTensorized quad( config.get() ); double my, phi, w; Vector moment = testBase.ComputeSphericalBasisKarthesian( 0, 1, 0 ); @@ -173,7 +174,7 @@ TEST_CASE( "test spherical harmonics basis ", "[spherical_harmonics]" ) { Vector moment2 = testBase.ComputeSphericalBasis( 0, 0 ); // Parity in carthesian coordinates - QGaussLegendreTensorized quad( config ); + QGaussLegendreTensorized quad( config.get() ); double x, y, z; diff --git a/tests/test_utilities.cpp b/tests/test_utilities.cpp new file mode 100644 index 00000000..9788f480 --- /dev/null +++ b/tests/test_utilities.cpp @@ -0,0 +1,121 @@ +#include "catch.hpp" +#include "toolboxes/textprocessingtoolbox.hpp" + +// Include reconstructor functions (declared in header after the class) +#include "toolboxes/reconstructor.hpp" + +TEST_CASE( "text processing utilities", "[utilities]" ) { + + SECTION( "Split" ) { + auto tokens = TextProcessingToolbox::Split( "hello,world,test", ',' ); + REQUIRE( tokens.size() == 3 ); + REQUIRE( tokens[0] == "hello" ); + REQUIRE( tokens[1] == "world" ); + REQUIRE( tokens[2] == "test" ); + + // Single element (no delimiter found) + auto single = TextProcessingToolbox::Split( "hello", ',' ); + REQUIRE( single.size() == 1 ); + REQUIRE( single[0] == "hello" ); + } + + SECTION( "StringEndsWith" ) { + REQUIRE( TextProcessingToolbox::StringEndsWith( "hello.vtk", ".vtk" ) ); + REQUIRE( TextProcessingToolbox::StringEndsWith( "test", "test" ) ); + REQUIRE_FALSE( TextProcessingToolbox::StringEndsWith( "hello.vtk", ".su2" ) ); + REQUIRE_FALSE( TextProcessingToolbox::StringEndsWith( "hi", "longer" ) ); + REQUIRE( TextProcessingToolbox::StringEndsWith( "anything", "" ) ); + } + + SECTION( "GetTrailingNumber" ) { + REQUIRE( TextProcessingToolbox::GetTrailingNumber( "NDIME= 2" ) == 2 ); + REQUIRE( TextProcessingToolbox::GetTrailingNumber( "NPOIN= 100" ) == 100 ); + REQUIRE( TextProcessingToolbox::GetTrailingNumber( "test42" ) == 42 ); + } + + SECTION( "StringToUpperCase" ) { + std::string s = "hello World"; + TextProcessingToolbox::StringToUpperCase( s ); + REQUIRE( s == "HELLO WORLD" ); + + std::string already = "UPPER"; + TextProcessingToolbox::StringToUpperCase( already ); + REQUIRE( already == "UPPER" ); + } + + SECTION( "DoubleToScientificNotation" ) { + // Should produce uppercase 'E' with 6 decimal places + std::string result = TextProcessingToolbox::DoubleToScientificNotation( 1.5e3 ); + REQUIRE( result.find( 'E' ) != std::string::npos ); + + // Check known value + std::string r2 = TextProcessingToolbox::DoubleToScientificNotation( 0.0 ); + REQUIRE( r2.find( 'E' ) != std::string::npos ); + } + + SECTION( "DoubleToScientificNotation2" ) { + // Should produce scientific notation with 4 decimal precision + std::string result = TextProcessingToolbox::DoubleToScientificNotation2( 1.0 ); + REQUIRE( result.find( 'e' ) != std::string::npos ); + } +} + +TEST_CASE( "slope limiter functions", "[utilities]" ) { + + SECTION( "FortSign" ) { + // FortSign(a, b) returns |a| if b > 0, -|a| if b < 0, 0 if b == 0 + REQUIRE( FortSign( 3.0, 1.0 ) == Approx( 3.0 ) ); + REQUIRE( FortSign( 3.0, -1.0 ) == Approx( -3.0 ) ); + REQUIRE( FortSign( -3.0, 1.0 ) == Approx( 3.0 ) ); + REQUIRE( FortSign( -3.0, -1.0 ) == Approx( -3.0 ) ); + REQUIRE( FortSign( 5.0, 0.0 ) == Approx( 0.0 ) ); + REQUIRE( FortSign( -5.0, 0.0 ) == Approx( 0.0 ) ); + } + + SECTION( "LMinMod" ) { + // Same sign: returns sign * min(|sL|, |sR|) + REQUIRE( LMinMod( 2.0, 3.0 ) == Approx( 2.0 ) ); + REQUIRE( LMinMod( 3.0, 2.0 ) == Approx( 2.0 ) ); + // Different signs: returns 0 + REQUIRE( LMinMod( 2.0, -3.0 ) == Approx( 0.0 ) ); + REQUIRE( LMinMod( -2.0, 3.0 ) == Approx( 0.0 ) ); + // Both negative: returns -min(|sL|, |sR|) + REQUIRE( LMinMod( -2.0, -3.0 ) == Approx( -2.0 ) ); + } + + SECTION( "LVanLeer" ) { + // Same sign: returns harmonic mean type + double result = LVanLeer( 2.0, 4.0 ); + // Expected: (1+1) * 2 * 4 / (2 + 4 + eps) ~= 2*8/6 ~= 2.667 + REQUIRE( result == Approx( 2.0 * 2.0 * 4.0 / ( 2.0 + 4.0 + 0.0000001 ) ).epsilon( 1e-5 ) ); + // Different signs: approximately 0 + REQUIRE( std::fabs( LVanLeer( 2.0, -4.0 ) ) < 1e-4 ); + } + + SECTION( "LVanAlbaba" ) { + // Symmetric: LVanAlbaba(a, b) should equal LVanAlbaba(b, a) for positive values + double r1 = LVanAlbaba( 2.0, 3.0 ); + double r2 = LVanAlbaba( 3.0, 2.0 ); + REQUIRE( r1 == Approx( r2 ) ); + // Known value: (sL^2*sR + sL*sR^2) / (sL^2 + sR^2 + eps) + double expected = ( 4.0 * 3.0 + 2.0 * 9.0 ) / ( 4.0 + 9.0 + 0.0000001 ); + REQUIRE( r1 == Approx( expected ).epsilon( 1e-5 ) ); + } + + SECTION( "LSuperBee" ) { + // When sR is between 0.5*sL and 2*sL, returns max + double r1 = LSuperBee( 2.0, 3.0 ); + // 3.0 >= 0.5*2.0 = 1.0 and 3.0 <= 2.0*2.0 = 4.0, so returns max(2,3) = 3 + REQUIRE( r1 == Approx( 3.0 ) ); + + // Different signs: returns 0 + REQUIRE( LSuperBee( 2.0, -1.0 ) == Approx( 0.0 ) ); + } + + SECTION( "LWENOJS" ) { + // Currently returns 0 for all inputs + REQUIRE( LWENOJS( 1.0 ) == Approx( 0.0 ) ); + REQUIRE( LWENOJS( -5.0 ) == Approx( 0.0 ) ); + REQUIRE( LWENOJS( 0.0 ) == Approx( 0.0 ) ); + } +}