diff --git a/CMakeLists.txt b/CMakeLists.txt index b9149d6..3833e3e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -175,6 +175,7 @@ set(SOURCE_FILES src/partition.cpp src/params.cpp src/simulation.cpp + src/debugging_utils.cpp ) # ========== Target ========== diff --git a/example_params.yml b/example_params.yml index 66c3e0f..9a2e000 100644 --- a/example_params.yml +++ b/example_params.yml @@ -13,6 +13,11 @@ Grid: n_grid_points: 1000000 # number of grid points to sample (applicable to random) random_seed: 42 # seed for random grid point generation (applicable to random) +Cosmology: + h: 0.681 # Reduced Hubble constant + Omega_cdm: 0.256011 # Cold dark matter density parameter + Omega_b: 0.048600 # Baryon density parameter + Performance: Tree: diff --git a/scripts/build_and_test_debug.sh b/scripts/build_and_test_debug.sh new file mode 100755 index 0000000..047775b --- /dev/null +++ b/scripts/build_and_test_debug.sh @@ -0,0 +1,104 @@ +#!/bin/bash +# +# Build in debug mode and run comprehensive tests +# This script compiles with DEBUGGING_CHECKS enabled and runs the full test suite +# + +set -e # Exit on error + +# Colors for output +RED='\033[0;31m' +GREEN='\033[0;32m' +YELLOW='\033[1;33m' +BLUE='\033[0;34m' +NC='\033[0m' # No Color + +echo -e "${BLUE}==============================================" +echo " BUILD AND TEST IN DEBUG MODE" +echo -e "==============================================${NC}" +echo "" + +# Get script directory and project root +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +PROJECT_ROOT="$(dirname "$SCRIPT_DIR")" + +cd "$PROJECT_ROOT" + +# Clean previous debug build +echo -e "${YELLOW}Cleaning previous debug build...${NC}" +rm -rf build_debug + +# Configure with debug mode +echo -e "${YELLOW}Configuring debug build...${NC}" +cmake -B build_debug -DCMAKE_BUILD_TYPE=Debug + +# Build +echo -e "${YELLOW}Building in debug mode...${NC}" +cmake --build build_debug + +if [ $? -ne 0 ]; then + echo -e "${RED}✗ Build failed${NC}" + exit 1 +fi + +echo -e "${GREEN}✓ Build successful${NC}" +echo "" + +# Check executable +if [ ! -f "build_debug/parent_gridder" ]; then + echo -e "${RED}✗ Executable not found: build_debug/parent_gridder${NC}" + exit 1 +fi + +echo -e "${BLUE}Running comprehensive test suite...${NC}" +echo "" + +# Run file-based gridding tests +echo -e "${YELLOW}=== File-Based Gridding Tests ===${NC}" +python3 tests/test_file_gridding.py build_debug/parent_gridder + +FILE_TEST_EXIT=$? + +echo "" + +# Run simple tests +echo -e "${YELLOW}=== Simple Test ===${NC}" +bash tests/run_simple_test.sh + +SIMPLE_TEST_EXIT=$? + +echo "" + +# Summary +echo -e "${BLUE}==============================================" +echo " TEST SUMMARY" +echo -e "==============================================${NC}" + +TOTAL_FAILED=0 + +if [ $FILE_TEST_EXIT -eq 0 ]; then + echo -e "${GREEN}✓${NC} File-based gridding tests: PASSED" +else + echo -e "${RED}✗${NC} File-based gridding tests: FAILED" + TOTAL_FAILED=$((TOTAL_FAILED + 1)) +fi + +if [ $SIMPLE_TEST_EXIT -eq 0 ]; then + echo -e "${GREEN}✓${NC} Simple tests: PASSED" +else + echo -e "${RED}✗${NC} Simple tests: FAILED" + TOTAL_FAILED=$((TOTAL_FAILED + 1)) +fi + +echo -e "${BLUE}==============================================${NC}" + +if [ $TOTAL_FAILED -eq 0 ]; then + echo -e "${GREEN}ALL TESTS PASSED!${NC}" + exit 0 +else + echo -e "${RED}$TOTAL_FAILED TEST SUITE(S) FAILED${NC}" + echo "" + echo "Debug build is available at: build_debug/parent_gridder" + echo "You can run individual tests or use a debugger to investigate failures." + exit 1 +fi diff --git a/src/cell.cpp b/src/cell.cpp index 790f0cd..1ba30c8 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -843,11 +843,6 @@ void assignPartsToCells(Simulation *sim) { total_mass = global_total_mass; #endif - // Compute the mean comoving density - sim->mean_density = total_mass / sim->volume; - - message("Mean comoving density: %e 10**10 Msun / cMpc^3", sim->mean_density); - #ifdef DEBUGGING_CHECKS // Make sure we have attached all the particles (only count local cells in // MPI) @@ -1222,6 +1217,27 @@ void assignGridPointsToCells([[maybe_unused]] Simulation *sim, Grid *grid) { // Get the cells for debugging checks std::vector &cells = sim->cells; + // Print cell assignment summary + message("[DEBUG] Grid point to cell assignment summary:"); + int cells_with_gps = 0; + for (size_t cid = 0; cid < sim->nr_cells; cid++) { + Cell *cell = &cells[cid]; + if (cell->grid_points.size() > 0) { + cells_with_gps++; + message("[DEBUG] Cell %zu at (%.3f, %.3f, %.3f) width (%.3f, %.3f, %.3f) has %zu grid points", + cid, cell->loc[0], cell->loc[1], cell->loc[2], + cell->width[0], cell->width[1], cell->width[2], + cell->grid_points.size()); + // Print first grid point in this cell + if (cell->grid_points.size() > 0) { + GridPoint *gp = cell->grid_points[0]; + message("[DEBUG] First GP: (%.6f, %.6f, %.6f)", + gp->loc[0], gp->loc[1], gp->loc[2]); + } + } + } + message("[DEBUG] Total cells with grid points: %d", cells_with_gps); + // Check grid points are in the right cells for (size_t cid = 0; cid < sim->nr_cells; cid++) { Cell *cell = &cells[cid]; diff --git a/src/construct_cells.cpp b/src/construct_cells.cpp index b6f648d..c94598b 100644 --- a/src/construct_cells.cpp +++ b/src/construct_cells.cpp @@ -47,23 +47,28 @@ void getTopCells(Simulation *sim, Grid *grid) { // neighbouring cells (this simplifies boilerplate elsewhere) // How many cells do we need to walk out for the biggest kernel? This is - // the maximum distance at which we will need to consider another cell - const int nwalk = std::ceil(grid->max_kernel_radius / width[0]) + 1; - int nwalk_upper = nwalk; - int nwalk_lower = nwalk; - - // If nwalk is greater than the number of cells in the simulation, we need - // to walk out to the edge of the simulation - if (nwalk > cdim[0] / 2) { - nwalk_upper = cdim[0] / 2; - nwalk_lower = cdim[0] / 2; + // the maximum distance at which we will need to consider another cell. + // We compute this separately for each dimension since cell widths and + // grid dimensions may differ. + int nwalk[3]; + for (int dim = 0; dim < 3; dim++) { + nwalk[dim] = std::ceil(grid->max_kernel_radius / width[dim]) + 1; + + // Clamp to half the grid dimension to prevent duplicate neighbors + // through periodic wrapping. If we walk more than cdim/2, we'll + // encounter the same cell from multiple periodic images. + if (nwalk[dim] > cdim[dim] / 2) { + nwalk[dim] = cdim[dim] / 2; + } } - message("Looking for neighbours within %d cells", nwalk); + message("Looking for neighbours within [%d, %d, %d] cells", + nwalk[0], nwalk[1], nwalk[2]); - // Calculate maximum neighbors + // Calculate maximum neighbors (use the maximum nwalk for reservation) + const int max_nwalk = std::max({nwalk[0], nwalk[1], nwalk[2]}); const int max_neighbors = - (2 * nwalk + 1) * (2 * nwalk + 1) * (2 * nwalk + 1) - + (2 * max_nwalk + 1) * (2 * max_nwalk + 1) * (2 * max_nwalk + 1) - 1; // -1 excludes self // Loop over the cells attaching the pointers the neighbouring cells (taking @@ -82,10 +87,11 @@ void getTopCells(Simulation *sim, Grid *grid) { // Reserve space for neighbors cell->neighbours.reserve(max_neighbors); - // Loop over the neighbours - for (int ii = -nwalk_lower; ii < nwalk_upper + 1; ii++) { - for (int jj = -nwalk_lower; jj < nwalk_upper + 1; jj++) { - for (int kk = -nwalk_lower; kk < nwalk_upper + 1; kk++) { + // Loop over the neighbours using dimension-specific nwalk values + // The nwalk values are already clamped to cdim/2, preventing duplicates + for (int ii = -nwalk[0]; ii <= nwalk[0]; ii++) { + for (int jj = -nwalk[1]; jj <= nwalk[1]; jj++) { + for (int kk = -nwalk[2]; kk <= nwalk[2]; kk++) { // Skip the cell itself if (ii == 0 && jj == 0 && kk == 0) @@ -97,6 +103,11 @@ void getTopCells(Simulation *sim, Grid *grid) { int kkk = (k + kk + cdim[2]) % cdim[2]; int cjd = iii * cdim[1] * cdim[2] + jjj * cdim[2] + kkk; + // Skip if this wraps back to the cell itself + // (can happen with periodic boundaries in small boxes) + if (cjd == static_cast(cid)) + continue; + // Attach the neighbour to the cell cell->neighbours.push_back(&cells[cjd]); } diff --git a/src/construct_grid_points.cpp b/src/construct_grid_points.cpp index b25b43a..326cce7 100644 --- a/src/construct_grid_points.cpp +++ b/src/construct_grid_points.cpp @@ -307,8 +307,8 @@ readGridPointCoordinates(const std::string &filename, message("WARNING: No valid grid point coordinates found in file: %s", filename.c_str()); message("The file is either empty or contains only comments/invalid lines."); - message("Exiting gracefully - no grid points to process."); - throw std::runtime_error("Empty grid file - no grid points to process"); + // Return 0 - the caller will handle the empty grid case + return 0; } message("Read %d valid grid point coordinates from %s", valid_points, @@ -413,6 +413,18 @@ static void createGridPointsFromFile(Simulation *sim, Grid *grid) { message("Created %d grid points from file %s", valid_points, grid->grid_file.c_str()); +#ifdef DEBUGGING_CHECKS + // Print first few grid points to verify they loaded correctly + message("[DEBUG] First 5 grid points loaded from file:"); + int n_to_print = std::min(5, valid_points); + for (int i = 0; i < n_to_print; i++) { + message("[DEBUG] Point %d: (%.6f, %.6f, %.6f)", i, + grid->grid_points[i].loc[0], + grid->grid_points[i].loc[1], + grid->grid_points[i].loc[2]); + } +#endif + // Initialize mass and count maps for all grid points message("Initializing grid point maps for %d kernel radii", grid->nkernels); for (GridPoint &gp : grid->grid_points) { diff --git a/src/debugging_utils.cpp b/src/debugging_utils.cpp new file mode 100644 index 0000000..6de561f --- /dev/null +++ b/src/debugging_utils.cpp @@ -0,0 +1,410 @@ +// Standard includes +#include +#include +#include + +// Local includes +#include "debugging_utils.hpp" +#include "cell.hpp" +#include "grid_point.hpp" +#include "logger.hpp" +#include "metadata.hpp" +#include "simulation.hpp" + +#ifdef DEBUGGING_CHECKS + +/** + * @brief Validate that grid points are correctly assigned to their containing + * cells + * + * This checks that each grid point is spatially within the bounds of the cell + * it's been assigned to. + */ +void validateGridPointCellAssignment(Simulation *sim, + [[maybe_unused]] Grid *grid) { + message("[DEBUG] Validating grid point to cell assignments..."); + + int errors = 0; + std::vector &cells = sim->cells; + + for (size_t cid = 0; cid < sim->nr_cells; cid++) { + Cell *cell = &cells[cid]; + + for (GridPoint *gp : cell->grid_points) { + // Check if grid point is within cell bounds + bool inside_x = (gp->loc[0] >= cell->loc[0]) && + (gp->loc[0] < cell->loc[0] + cell->width[0]); + bool inside_y = (gp->loc[1] >= cell->loc[1]) && + (gp->loc[1] < cell->loc[1] + cell->width[1]); + bool inside_z = (gp->loc[2] >= cell->loc[2]) && + (gp->loc[2] < cell->loc[2] + cell->width[2]); + + if (!inside_x || !inside_y || !inside_z) { + message("[DEBUG] ERROR: Grid point at (%.3f, %.3f, %.3f) assigned to " + "cell %zu " + "with bounds [%.3f-%.3f, %.3f-%.3f, %.3f-%.3f]", + gp->loc[0], gp->loc[1], gp->loc[2], cid, cell->loc[0], + cell->loc[0] + cell->width[0], cell->loc[1], + cell->loc[1] + cell->width[1], cell->loc[2], + cell->loc[2] + cell->width[2]); + errors++; + } + } + } + + if (errors > 0) { + error("[DEBUG] Found %d grid points incorrectly assigned to cells", errors); + } else { + message("[DEBUG] ✓ All grid points correctly assigned to cells"); + } +} + +/** + * @brief Validate that all useful cells have been correctly identified + * + * Checks that: + * 1. All cells with grid points are marked as useful + * 2. All neighbors of cells with grid points are marked as useful + */ +void validateUsefulCells(Simulation *sim, [[maybe_unused]] Grid *grid) { + message("[DEBUG] Validating useful cell flagging..."); + + int errors = 0; + std::vector &cells = sim->cells; + + for (size_t cid = 0; cid < sim->nr_cells; cid++) { + Cell *cell = &cells[cid]; + + // Check 1: Cells with grid points must be useful + if (cell->grid_points.size() > 0 && !cell->is_useful) { + message("[DEBUG] ERROR: Cell %zu has %zu grid points but is_useful = " + "false", + cid, cell->grid_points.size()); + errors++; + } + + // Check 2: Neighbors of cells with grid points must be useful + if (cell->grid_points.size() > 0) { + for (Cell *neighbor : cell->neighbours) { + if (!neighbor->is_useful) { + message( + "[DEBUG] ERROR: Cell %zu has grid points but its neighbor " + "(cell at %.3f,%.3f,%.3f) is not marked as useful", + cid, neighbor->loc[0], neighbor->loc[1], neighbor->loc[2]); + errors++; + } + } + } + } + + // Count total useful cells and cells with grid points + int useful_count = 0; + int cells_with_gps = 0; + for (size_t cid = 0; cid < sim->nr_cells; cid++) { + if (cells[cid].is_useful) + useful_count++; + if (cells[cid].grid_points.size() > 0) + cells_with_gps++; + } + + message("[DEBUG] Useful cells: %d, Cells with grid points: %d", useful_count, + cells_with_gps); + + if (errors > 0) { + error("[DEBUG] Found %d useful cell flagging errors", errors); + } else { + message("[DEBUG] ✓ All useful cells correctly flagged"); + } +} + +/** + * @brief Validate that particles are in the correct cells + * + * Checks that each particle's position is within the bounds of its assigned + * cell. + */ +void validateParticleCellAssignment(Simulation *sim) { + message("[DEBUG] Validating particle to cell assignments..."); + + int errors = 0; + std::vector &cells = sim->cells; + + for (size_t cid = 0; cid < sim->nr_cells; cid++) { + Cell *cell = &cells[cid]; + + for (Particle *part : cell->particles) { + // Check if particle is within cell bounds + bool inside_x = (part->pos[0] >= cell->loc[0]) && + (part->pos[0] < cell->loc[0] + cell->width[0]); + bool inside_y = (part->pos[1] >= cell->loc[1]) && + (part->pos[1] < cell->loc[1] + cell->width[1]); + bool inside_z = (part->pos[2] >= cell->loc[2]) && + (part->pos[2] < cell->loc[2] + cell->width[2]); + + if (!inside_x || !inside_y || !inside_z) { + message("[DEBUG] ERROR: Particle at (%.3f, %.3f, %.3f) assigned to " + "cell %zu " + "with bounds [%.3f-%.3f, %.3f-%.3f, %.3f-%.3f]", + part->pos[0], part->pos[1], part->pos[2], cid, cell->loc[0], + cell->loc[0] + cell->width[0], cell->loc[1], + cell->loc[1] + cell->width[1], cell->loc[2], + cell->loc[2] + cell->width[2]); + errors++; + if (errors > 10) { // Limit error spam + message("[DEBUG] ... (suppressing further particle assignment " + "errors)"); + break; + } + } + } + if (errors > 10) + break; + } + + if (errors > 0) { + error("[DEBUG] Found %d particles incorrectly assigned to cells", errors); + } else { + message("[DEBUG] ✓ All particles correctly assigned to cells"); + } +} + +/** + * @brief Count particles within radius of a grid point (brute force check) + */ +int bruteForceCountParticles(GridPoint *grid_point, Simulation *sim, + double radius) { + int count = 0; + double radius2 = radius * radius; + double *dim = sim->dim; + + std::vector &cells = sim->cells; + for (size_t cid = 0; cid < sim->nr_cells; cid++) { + Cell *cell = &cells[cid]; + + for (Particle *part : cell->particles) { + double dx = nearest(part->pos[0] - grid_point->loc[0], dim[0]); + double dy = nearest(part->pos[1] - grid_point->loc[1], dim[1]); + double dz = nearest(part->pos[2] - grid_point->loc[2], dim[2]); + double r2 = dx * dx + dy * dy + dz * dz; + + if (r2 <= radius2) { + count++; + } + } + } + + return count; +} + +/** + * @brief Check if grid points can find any particles within kernel radii + * + * For ALL grid points, performs a brute force search to count particles within + * each kernel radius, stores the counts, then compares to what the gridder found. + * In debug mode, these brute force counts will be written to the output file. + */ +void validateGridPointsHaveParticles(Simulation *sim, Grid *grid) { + message("[DEBUG] Computing brute force particle counts for all grid points and all kernels..."); + message("[DEBUG] This will be written to output file for comparison."); + + int total_mismatches = 0; + int total_empty_kernels = 0; + int total_checks = 0; + + // For ALL grid points (not sampled), compute brute force counts for ALL kernels + for (size_t i = 0; i < grid->grid_points.size(); i++) { + GridPoint *gp = &grid->grid_points[i]; + + // Compute brute force for each kernel radius + for (double kernel_rad : grid->kernel_radii) { + int brute_count = bruteForceCountParticles(gp, sim, kernel_rad); + int gridder_count = gp->getCount(kernel_rad); + + // Store the brute force count + gp->setBruteForceCount(kernel_rad, brute_count); + + total_checks++; + + if (brute_count == 0) { + total_empty_kernels++; + } + + if (brute_count != gridder_count) { + total_mismatches++; + } + } + } + + message("[DEBUG] Brute force validation complete:"); + message("[DEBUG] - Checked %zu grid points × %d kernels = %d combinations", + grid->grid_points.size(), grid->nkernels, total_checks); + message("[DEBUG] - %d had no particles (expected for small kernels/sparse regions)", + total_empty_kernels); + message("[DEBUG] - %d had count mismatches between brute force and gridder", + total_mismatches); + + if (total_mismatches > 0) { + message("[DEBUG] ╔═══════════════════════════════════════════════════════════════╗"); + message("[DEBUG] ║ WARNING: PARTICLE COUNT MISMATCHES DETECTED ║"); + message("[DEBUG] ╠═══════════════════════════════════════════════════════════════╣"); + message("[DEBUG] ║ Found %5d mismatches between octree and brute force ║", total_mismatches); + message("[DEBUG] ║ ║"); + message("[DEBUG] ║ THE RESULTS IN THE OUTPUT FILE ARE INCORRECT! ║"); + message("[DEBUG] ║ ║"); + message("[DEBUG] ║ The output file contains both: ║"); + message("[DEBUG] ║ - GridPointCounts: ║"); + message("[DEBUG] ║ Counts from octree traversal (INCORRECT - has bugs) ║"); + message("[DEBUG] ║ - BruteForceGridPointCounts: ║"); + message("[DEBUG] ║ Counts from exhaustive search (CORRECT - ground truth) ║"); + message("[DEBUG] ║ ║"); + message("[DEBUG] ║ Use BruteForceGridPointCounts for analysis. ║"); + message("[DEBUG] ║ ║"); + message("[DEBUG] ║ This likely indicates a bug in periodic boundary handling ║"); + message("[DEBUG] ║ or octree traversal, especially for grid points near box ║"); + message("[DEBUG] ║ boundaries or with small kernel radii. ║"); + message("[DEBUG] ╚═══════════════════════════════════════════════════════════════╝"); + } else if (total_empty_kernels > total_checks / 2) { + message("[DEBUG] WARNING: More than 50%% of grid point/kernel combinations " + "have no particles nearby"); + } else { + message("[DEBUG] ✓ Grid point particle counts match brute force validation"); + } +} + +/** + * @brief Validate octree structure integrity + * + * Recursively checks that: + * 1. Child cells properly partition their parent + * 2. Particle counts match between parent and children + * 3. Grid point counts match between parent and children + */ +static void validateCellRecursive(Cell *cell, int *errors) { + if (!cell->is_split) + return; + + // Check particle count consistency + size_t child_part_count = 0; + for (int i = 0; i < Cell::OCTREE_CHILDREN; i++) { + if (cell->children[i] != nullptr) { + child_part_count += cell->children[i]->part_count; + validateCellRecursive(cell->children[i], errors); + } + } + + if (child_part_count != cell->part_count) { + message( + "[DEBUG] ERROR: Cell at (%.3f, %.3f, %.3f) has %zu particles but " + "children have %zu total", + cell->loc[0], cell->loc[1], cell->loc[2], cell->part_count, + child_part_count); + (*errors)++; + } + + // Check grid point count consistency + size_t child_gp_count = 0; + for (int i = 0; i < Cell::OCTREE_CHILDREN; i++) { + if (cell->children[i] != nullptr) { + child_gp_count += cell->children[i]->grid_points.size(); + } + } + + if (child_gp_count != cell->grid_points.size()) { + message("[DEBUG] ERROR: Cell at (%.3f, %.3f, %.3f) has %zu grid points " + "but children have %zu total", + cell->loc[0], cell->loc[1], cell->loc[2], cell->grid_points.size(), + child_gp_count); + (*errors)++; + } +} + +void validateOctreeStructure(Simulation *sim) { + message("[DEBUG] Validating octree structure..."); + + int errors = 0; + std::vector &cells = sim->cells; + + for (size_t cid = 0; cid < sim->nr_cells; cid++) { + validateCellRecursive(&cells[cid], &errors); + } + + if (errors > 0) { + error("[DEBUG] Found %d octree structure errors", errors); + } else { + message("[DEBUG] ✓ Octree structure is valid"); + } +} + +/** + * @brief Check that file-based grid points are within simulation boundaries + */ +void validateFileGridPoints(Grid *grid, Simulation *sim) { + if (!grid->grid_from_file) + return; + + message("[DEBUG] Validating file-based grid points are within simulation " + "boundaries..."); + + int outside_count = 0; + double *dim = sim->dim; + + for (GridPoint &gp : grid->grid_points) { + if (gp.loc[0] < 0 || gp.loc[0] >= dim[0] || gp.loc[1] < 0 || + gp.loc[1] >= dim[1] || gp.loc[2] < 0 || gp.loc[2] >= dim[2]) { + message("[DEBUG] ERROR: Grid point at (%.3f, %.3f, %.3f) is outside " + "simulation box [0-%.3f, 0-%.3f, 0-%.3f]", + gp.loc[0], gp.loc[1], gp.loc[2], dim[0], dim[1], dim[2]); + outside_count++; + } + } + + if (outside_count > 0) { + error("[DEBUG] Found %d grid points outside simulation boundaries", + outside_count); + } else { + message("[DEBUG] ✓ All file-based grid points within simulation " + "boundaries"); + } +} + +/** + * @brief Print detailed diagnostics about a specific grid point + */ +void diagnoseGridPoint(GridPoint *grid_point, Simulation *sim, Grid *grid) { + message("[DEBUG] ==== Grid Point Diagnostic ===="); + message("[DEBUG] Location: (%.6f, %.6f, %.6f)", grid_point->loc[0], + grid_point->loc[1], grid_point->loc[2]); + + // Find which cell it's in + Cell *cell = getCellContainingPoint(grid_point->loc); + message("[DEBUG] Assigned to cell at (%.3f, %.3f, %.3f) with width (%.3f, " + "%.3f, %.3f)", + cell->loc[0], cell->loc[1], cell->loc[2], cell->width[0], + cell->width[1], cell->width[2]); + message("[DEBUG] Cell has %zu particles, %zu grid points", cell->part_count, + cell->grid_points.size()); + message("[DEBUG] Cell is_useful: %s, is_split: %s", + cell->is_useful ? "true" : "false", cell->is_split ? "true" : "false"); + + // Check neighbors + message("[DEBUG] Cell has %zu neighbors", cell->neighbours.size()); + size_t total_neighbor_particles = 0; + for (Cell *neighbor : cell->neighbours) { + total_neighbor_particles += neighbor->part_count; + } + message("[DEBUG] Neighbors have %zu total particles", total_neighbor_particles); + + // Check particle counts for each kernel + message("[DEBUG] Particle counts by kernel radius:"); + for (double radius : grid->kernel_radii) { + // Get count from gridder + int count = grid_point->getCount(radius); + int brute_count = bruteForceCountParticles(grid_point, sim, radius); + message("[DEBUG] %.3f Mpc/h: gridder=%d, brute_force=%d", radius, count, + brute_count); + } + + message("[DEBUG] ================================"); +} + +#endif // DEBUGGING_CHECKS diff --git a/src/debugging_utils.hpp b/src/debugging_utils.hpp new file mode 100644 index 0000000..39040bf --- /dev/null +++ b/src/debugging_utils.hpp @@ -0,0 +1,85 @@ +#ifndef DEBUGGING_UTILS_HPP +#define DEBUGGING_UTILS_HPP + +// Header for comprehensive debugging utilities +// These are active only when DEBUGGING_CHECKS is defined + +#include "cell.hpp" +#include "grid_point.hpp" +#include "logger.hpp" +#include "metadata.hpp" +#include "simulation.hpp" +#include + +#ifdef DEBUGGING_CHECKS + +/** + * @brief Validate that grid points are correctly assigned to their containing + * cells + * + * @param sim Simulation object + * @param grid Grid object + */ +void validateGridPointCellAssignment(Simulation *sim, Grid *grid); + +/** + * @brief Validate that all useful cells have been correctly identified + * + * @param sim Simulation object + * @param grid Grid object + */ +void validateUsefulCells(Simulation *sim, Grid *grid); + +/** + * @brief Validate that particles are in the correct cells + * + * @param sim Simulation object + */ +void validateParticleCellAssignment(Simulation *sim); + +/** + * @brief Check if grid points can find any particles within kernel radii + * + * @param sim Simulation object + * @param grid Grid object + */ +void validateGridPointsHaveParticles(Simulation *sim, Grid *grid); + +/** + * @brief Validate octree structure integrity + * + * @param sim Simulation object + */ +void validateOctreeStructure(Simulation *sim); + +/** + * @brief Check that file-based grid points are within simulation boundaries + * + * @param grid Grid object + * @param sim Simulation object + */ +void validateFileGridPoints(Grid *grid, Simulation *sim); + +/** + * @brief Print detailed diagnostics about a specific grid point + * + * @param grid_point The grid point to diagnose + * @param sim Simulation object + * @param grid Grid object + */ +void diagnoseGridPoint(GridPoint *grid_point, Simulation *sim, Grid *grid); + +/** + * @brief Count particles within radius of a grid point (brute force) + * + * @param grid_point The grid point + * @param sim Simulation object + * @param radius Search radius + * @return Number of particles found + */ +int bruteForceCountParticles(GridPoint *grid_point, Simulation *sim, + double radius); + +#endif // DEBUGGING_CHECKS + +#endif // DEBUGGING_UTILS_HPP diff --git a/src/grid_point.cpp b/src/grid_point.cpp index f4a9e67..e33afdf 100644 --- a/src/grid_point.cpp +++ b/src/grid_point.cpp @@ -86,6 +86,34 @@ double GridPoint::getMass(const double kernel_radius) const { } } +// Method to get the particle count inside the kernel radius +int GridPoint::getCount(const double kernel_radius) const { + // Check if the kernel radius exists in the count map + auto it = this->count_map.find(kernel_radius); + if (it != this->count_map.end()) { + return static_cast(it->second); + } else { + return 0; // Return 0 if the kernel radius is not found + } +} + +#ifdef DEBUGGING_CHECKS +// Method to set the brute force count for a kernel radius +void GridPoint::setBruteForceCount(const double kernel_radius, int count) { + this->brute_force_count_map[kernel_radius] = count; +} + +// Method to get the brute force count for a kernel radius +int GridPoint::getBruteForceCount(const double kernel_radius) const { + auto it = this->brute_force_count_map.find(kernel_radius); + if (it != this->brute_force_count_map.end()) { + return it->second; + } else { + return -1; // Return -1 if not computed (shouldn't happen) + } +} +#endif + /** * @brief Construct a new Grid object * diff --git a/src/grid_point.hpp b/src/grid_point.hpp index ede437e..375832e 100644 --- a/src/grid_point.hpp +++ b/src/grid_point.hpp @@ -28,6 +28,12 @@ class GridPoint { double kernel_radius); double getOverDensity(const double kernel_radius, Simulation *sim) const; double getMass(const double kernel_radius) const; + int getCount(const double kernel_radius) const; + +#ifdef DEBUGGING_CHECKS + void setBruteForceCount(const double kernel_radius, int count); + int getBruteForceCount(const double kernel_radius) const; +#endif private: //! The count of particles in each kernel radius @@ -35,6 +41,11 @@ class GridPoint { //! The mass of particles in each kernel radius std::unordered_map mass_map; + +#ifdef DEBUGGING_CHECKS + //! Brute force particle counts (for validation in debug mode) + std::unordered_map brute_force_count_map; +#endif }; class Grid { diff --git a/src/gridder.cpp b/src/gridder.cpp index c18fecb..2d5f69b 100644 --- a/src/gridder.cpp +++ b/src/gridder.cpp @@ -22,6 +22,10 @@ #include "simulation.hpp" #include "talking.hpp" +#ifdef DEBUGGING_CHECKS +#include "debugging_utils.hpp" +#endif + /** * @brief Function to handle the command line arguments using robust parser * @@ -215,6 +219,15 @@ int main(int argc, char *argv[]) { message("Number of top level cells: %d", sim->nr_cells); + // Calculate mean density from cosmological parameters + // This must be done before any overdensity calculations + try { + sim->calculateMeanDensityFromCosmology(params); + } catch (const std::exception &e) { + std::cerr << "Failed to calculate mean density from cosmology: " << e.what() << std::endl; + return 1; + } + #ifdef WITH_MPI // Partition the cells across MPI ranks (work-based partition) try { @@ -239,10 +252,34 @@ int main(int argc, char *argv[]) { // Check if we actually have grid points to process if (grid->n_grid_points == 0) { message("No grid points available for processing."); - message("Program will now exit."); - return 0; // Clean exit, not an error + message("Writing empty output file before exiting..."); + + // Write empty output file +#ifdef WITH_MPI + try { + writeGridFileParallel(sim, grid); + } catch (const std::exception &e) { + error("Failed to write empty output file: %s", e.what()); + return 1; + } +#else + try { + writeGridFileSerial(sim, grid); + } catch (const std::exception &e) { + error("Failed to write empty output file: %s", e.what()); + return 1; + } +#endif + + message("No grid points were created - empty output file written."); + return 1; // Exit with error code since no work was done } +#ifdef DEBUGGING_CHECKS + // Validate file-based grid points are within simulation boundaries + validateFileGridPoints(grid, sim); +#endif + #ifdef WITH_MPI MPI_Barrier(MPI_COMM_WORLD); #endif @@ -255,6 +292,11 @@ int main(int argc, char *argv[]) { return 1; } +#ifdef DEBUGGING_CHECKS + // Validate grid points are correctly assigned to cells + validateGridPointCellAssignment(sim, grid); +#endif + #ifdef WITH_MPI MPI_Barrier(MPI_COMM_WORLD); #endif @@ -269,6 +311,11 @@ int main(int argc, char *argv[]) { return 1; } +#ifdef DEBUGGING_CHECKS + // Validate useful cells are correctly flagged + validateUsefulCells(sim, grid); +#endif + #ifdef WITH_MPI MPI_Barrier(MPI_COMM_WORLD); @@ -322,6 +369,11 @@ int main(int argc, char *argv[]) { return 1; } +#ifdef DEBUGGING_CHECKS + // Validate particles are in the correct cells + validateParticleCellAssignment(sim); +#endif + #ifdef WITH_MPI MPI_Barrier(MPI_COMM_WORLD); @@ -360,6 +412,11 @@ int main(int argc, char *argv[]) { return 1; } +#ifdef DEBUGGING_CHECKS + // Validate octree structure integrity + validateOctreeStructure(sim); +#endif + #ifdef WITH_MPI MPI_Barrier(MPI_COMM_WORLD); #endif @@ -373,6 +430,11 @@ int main(int argc, char *argv[]) { return 1; } +#ifdef DEBUGGING_CHECKS + // Validate grid points have particles (sample check with brute force) + validateGridPointsHaveParticles(sim, grid); +#endif + #ifdef WITH_MPI MPI_Barrier(MPI_COMM_WORLD); #endif diff --git a/src/output.cpp b/src/output.cpp index 6ffb54f..fdd9cfa 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -98,6 +98,15 @@ void writeGridFileSerial(Simulation *sim, Grid *grid) { return; } + // If there are no grid points, skip dataset creation and exit early + if (grid->n_grid_points == 0) { + message("No grid points to write - creating empty output file with header only"); + hdf5.close(); + message("Successfully wrote empty grid data (serial mode)"); + toc("Writing output (in serial)"); + return; + } + // Create dataset for grid positions std::array grid_point_positions_dims = { static_cast(grid->n_grid_points), static_cast(3)}; @@ -148,6 +157,27 @@ void writeGridFileSerial(Simulation *sim, Grid *grid) { } } + // Create particle count dataset from tree traversal algorithm + if (!hdf5.createDataset("Grids/" + kernel_name + + "/GridPointCounts", + grid_point_overdens_dims)) { + error("Failed to create GridPointCounts dataset for kernel %zu (radius=%f)", + kernel_idx, kernel_rad); + continue; + } + +#ifdef DEBUGGING_CHECKS + // In debug mode, also write brute force counts for validation + // These are computed by exhaustively checking every particle's distance + if (!hdf5.createDataset("Grids/" + kernel_name + + "/BruteForceGridPointCounts", + grid_point_overdens_dims)) { + error("Failed to create BruteForceGridPointCounts dataset for kernel %zu (radius=%f)", + kernel_idx, kernel_rad); + continue; + } +#endif + // Process each cell for (size_t cid = 0; cid < sim->nr_cells; cid++) { Cell *cell = &cells[cid]; @@ -163,12 +193,18 @@ void writeGridFileSerial(Simulation *sim, Grid *grid) { // Prepare data arrays std::vector cell_grid_overdens; std::vector cell_grid_masses; + std::vector cell_grid_counts; std::vector cell_grid_pos; cell_grid_overdens.reserve(count); + cell_grid_counts.reserve(count); cell_grid_pos.reserve(count * 3); if (metadata->write_masses) { cell_grid_masses.reserve(count); } +#ifdef DEBUGGING_CHECKS + std::vector cell_grid_brute_counts; + cell_grid_brute_counts.reserve(count); +#endif // Extract data from grid points for (const GridPoint *gp : cell->grid_points) { @@ -180,6 +216,14 @@ void writeGridFileSerial(Simulation *sim, Grid *grid) { cell_grid_masses.push_back(gp->getMass(kernel_rad)); } + // Get particle counts from gridder algorithm + cell_grid_counts.push_back(gp->getCount(kernel_rad)); + +#ifdef DEBUGGING_CHECKS + // Get brute force counts in debug mode + cell_grid_brute_counts.push_back(gp->getBruteForceCount(kernel_rad)); +#endif + // Store positions if not done yet if (!written_positions) { cell_grid_pos.push_back(gp->loc[0]); @@ -209,6 +253,30 @@ void writeGridFileSerial(Simulation *sim, Grid *grid) { } } + // Write gridder particle counts + if (!cell_grid_counts.empty()) { + if (!hdf5.writeDatasetSlice( + "Grids/" + kernel_name + "/GridPointCounts", cell_grid_counts, + {static_cast(start_idx)}, + {static_cast(count)})) { + error("Failed to write counts slice for cell %zu, kernel %zu (radius=%f)", + cid, kernel_idx, kernel_rad); + } + } + +#ifdef DEBUGGING_CHECKS + // Write brute force counts in debug mode + if (!cell_grid_brute_counts.empty()) { + if (!hdf5.writeDatasetSlice( + "Grids/" + kernel_name + "/BruteForceGridPointCounts", cell_grid_brute_counts, + {static_cast(start_idx)}, + {static_cast(count)})) { + error("Failed to write brute force counts slice for cell %zu, kernel %zu (radius=%f)", + cid, kernel_idx, kernel_rad); + } + } +#endif + // Write position slice if needed if (!written_positions && !cell_grid_pos.empty()) { if (!hdf5.writeDatasetSlice( @@ -331,6 +399,16 @@ void writeGridFileParallel(Simulation *sim, Grid *grid) { hdf5.writeAttribute("Header", "MaxKernelRadius", grid->max_kernel_radius); + // If there are no local grid points, skip dataset creation and exit early + if (total_local_grid_points == 0) { + message("Rank %d: No local grid points to write - creating empty output file with header only", + metadata->rank); + hdf5.close(); + message("Rank %d: Successfully wrote empty grid data (parallel mode)", metadata->rank); + toc("Writing output (in parallel)"); + return; + } + // Write local cell information if (!local_cell_ids.empty()) { std::array local_cell_dims = { diff --git a/src/simulation.cpp b/src/simulation.cpp index 912db54..5eab1ed 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,9 +1,17 @@ +// Standard includes +#include + // Local includes #include "simulation.hpp" #include "cell.hpp" #include "hdf_io.hpp" #include "metadata.hpp" +// Define M_PI if not available (POSIX extension, not standard C++) +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + /** * @brief Construct a new Simulation object * @@ -105,3 +113,54 @@ void Simulation::readSimulationData() { // Compute the comoving volume of the simulation this->volume = this->dim[0] * this->dim[1] * this->dim[2]; } + +/** + * @brief Calculate mean comoving density from cosmological parameters + * + * Computes the mean matter density at the simulation redshift using: + * ρ_mean = ρ_crit(z=0) × Ω_m × (1+z)³ + * + * where ρ_crit(z=0) = 3H₀²/(8πG) is the critical density today + * + * @param params The parameters object containing cosmology + */ +void Simulation::calculateMeanDensityFromCosmology(Parameters *params) { + + // Read cosmology parameters + double h = params->getParameterNoDefault("Cosmology/h"); + double Omega_cdm = params->getParameterNoDefault("Cosmology/Omega_cdm"); + double Omega_b = params->getParameterNoDefault("Cosmology/Omega_b"); + + // Total matter density parameter + double Omega_m = Omega_cdm + Omega_b; + + // Physical constants in internal units (10^10 Msun, Mpc, km/s) + // H0 = 100 h km/s/Mpc + double H0_kmsMpc = 100.0 * h; // km/s/Mpc + + // Convert to internal time units (H0 in units of 1/time where time is Mpc/(km/s)) + // H0 = 100 h km/s/Mpc = 100 h / Mpc * (km/s) + // In our units: [H0] = km/s/Mpc + + // Critical density today: ρ_crit = 3H₀²/(8πG) + // G in (10^10 Msun)^-1 Mpc (km/s)^2 internal units + // Derived from G_SI = 6.674e-11 m^3 kg^-1 s^-2 with proper unit conversion + const double G = 4.301744232015554e+01; // Gravitational constant in (10^10 Msun)^-1 Mpc (km/s)^2 + + // ρ_crit(z=0) = 3H₀²/(8πG) in units of 10^10 Msun / Mpc^3 + double rho_crit_0 = (3.0 * H0_kmsMpc * H0_kmsMpc) / (8.0 * M_PI * G); + + // Mean COMOVING density: ρ_comoving = ρ_crit(0) × Ω_m + // Note: In comoving coordinates, density does NOT evolve with redshift + // The (1+z)³ factor would convert to physical density, but SWIFT uses comoving coordinates + this->mean_density = rho_crit_0 * Omega_m; + + Metadata *metadata = &Metadata::getInstance(); + if (metadata->rank == 0) { + message("Cosmology: h=%.4f, Omega_m=%.6f (Omega_cdm=%.6f + Omega_b=%.6f)", + h, Omega_m, Omega_cdm, Omega_b); + message("Critical density today: %.6e 10^10 Msun/Mpc^3", rho_crit_0); + message("Mean comoving density at z=%.4f: %.6e 10^10 Msun/cMpc^3", + this->redshift, this->mean_density); + } +} diff --git a/src/simulation.hpp b/src/simulation.hpp index 729a933..900bfa2 100644 --- a/src/simulation.hpp +++ b/src/simulation.hpp @@ -76,6 +76,9 @@ class Simulation { // Prototype for reader function (defined in simulation.cpp) void readSimulationData(); + // Calculate mean density from cosmological parameters + void calculateMeanDensityFromCosmology(Parameters *params); + private: // Helper function for cleanup void deleteChildCells(Cell *cell); diff --git a/tests/analyze_debug_output.py b/tests/analyze_debug_output.py new file mode 100644 index 0000000..e140510 --- /dev/null +++ b/tests/analyze_debug_output.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +""" +Analyze debug output to diagnose tree traversal bugs +""" +import h5py +import numpy as np +import sys + +def analyze_file(filename): + print('='*70) + print(f'ANALYZING: {filename}') + print('='*70) + + with h5py.File(filename, 'r') as f: + # Get grid point positions + positions = np.array(f['Grids/GridPointPositions']) + + print(f'\nGrid point locations:') + for i, pos in enumerate(positions): + print(f' GP {i:2d}: ({pos[0]:5.1f}, {pos[1]:5.1f}, {pos[2]:5.1f})') + + print('\n' + '='*70) + print('Per-kernel comparison:') + print('='*70) + + for kernel_name in sorted(f['Grids'].keys()): + if kernel_name == 'GridPointPositions': + continue + + kernel_rad = f[f'Grids/{kernel_name}'].attrs['KernelRadius'] + octree_counts = np.array(f[f'Grids/{kernel_name}/GridPointCounts']) + brute_counts = np.array(f[f'Grids/{kernel_name}/BruteForceGridPointCounts']) + + print(f'\n{kernel_name} (radius = {kernel_rad} Mpc/h):') + print(f' GP | Position | Octree | Brute | Diff | Status') + print(f' ----+----------------------+--------+--------+--------+-------') + + for i in range(len(octree_counts)): + pos_str = f'({positions[i][0]:5.1f},{positions[i][1]:5.1f},{positions[i][2]:5.1f})' + diff = octree_counts[i] - brute_counts[i] + + if diff == 0: + status = 'OK' + elif octree_counts[i] == 0 and brute_counts[i] > 0: + status = 'MISS' # Octree missed particles + elif octree_counts[i] > brute_counts[i]: + status = 'OVER' # Octree over-counted + else: + status = 'UNDER' # Octree under-counted + + print(f' {i:3d} | {pos_str:20s} | {octree_counts[i]:6d} | {brute_counts[i]:6d} | {diff:+7d} | {status}') + + # Summary + total_mismatches = np.sum(octree_counts != brute_counts) + total_miss = np.sum((octree_counts == 0) & (brute_counts > 0)) + total_over = np.sum(octree_counts > brute_counts) + + print(f'\n Summary: {total_mismatches}/{len(octree_counts)} mismatches') + print(f' - Complete misses (0 when should find particles): {total_miss}') + print(f' - Over-counting (found more than truth): {total_over}') + +if __name__ == '__main__': + if len(sys.argv) > 1: + analyze_file(sys.argv[1]) + else: + # Default files to analyze + files = [ + 'tests/data/boundary_output.hdf5', + 'tests/data/cluster_output.hdf5', + ] + for f in files: + try: + analyze_file(f) + print('\n') + except FileNotFoundError: + print(f'File not found: {f}\n') diff --git a/tests/file_grid_test_params.yml b/tests/file_grid_test_params.yml index 44cd162..b3f8f55 100644 --- a/tests/file_grid_test_params.yml +++ b/tests/file_grid_test_params.yml @@ -9,6 +9,12 @@ Grid: grid_file: tests/data/grid_points_files/valid_points.txt # n_grid_points is optional for file type - will be counted from file +Cosmology: + h: 0.681 # Reduced Hubble constant + Omega_cdm: 0.256011 # Cold dark matter density parameter + Omega_b: 0.048600 # Baryon density parameter + + Tree: max_leaf_count: 200 diff --git a/make_test_snap.py b/tests/make_test_snap.py similarity index 100% rename from make_test_snap.py rename to tests/make_test_snap.py diff --git a/tests/run_simple_test.sh b/tests/run_simple_test.sh index 473e9ca..c4f5633 100755 --- a/tests/run_simple_test.sh +++ b/tests/run_simple_test.sh @@ -67,7 +67,7 @@ run_test() { # Create simple test snapshot cd "$ROOT_DIR" - python3 make_test_snap.py \ + python3 tests/make_test_snap.py \ --output "$TEST_SNAPSHOT" \ --cdim 3 \ --boxsize 10.0 \ diff --git a/tests/simple_test_params.yml b/tests/simple_test_params.yml index 33e6d01..db1956f 100644 --- a/tests/simple_test_params.yml +++ b/tests/simple_test_params.yml @@ -8,6 +8,12 @@ Grid: cdim: 3 # Small grid (3x3x3), grid point at center type: uniform +Cosmology: + h: 0.681 # Reduced Hubble constant + Omega_cdm: 0.256011 # Cold dark matter density parameter + Omega_b: 0.048600 # Baryon density parameter + + Tree: max_leaf_count: 200 diff --git a/tests/test_cosmology.py b/tests/test_cosmology.py new file mode 100644 index 0000000..4a1e197 --- /dev/null +++ b/tests/test_cosmology.py @@ -0,0 +1,395 @@ +#!/usr/bin/env python3 +""" +Unit tests for cosmological mean density calculations. + +Tests verify that the C++ cosmology implementation correctly computes +mean comoving density at various redshifts and with different cosmological +parameters. +""" + +import subprocess +import pytest +import h5py +import numpy as np +from pathlib import Path +import yaml +import tempfile + +# Try to import astropy for direct comparison tests +try: + from astropy.cosmology import FlatLambdaCDM + from astropy import units as u + ASTROPY_AVAILABLE = True +except ImportError: + ASTROPY_AVAILABLE = False + +# Get paths +PROJECT_ROOT = Path(__file__).parent.parent +BUILD_DIR = PROJECT_ROOT / "build" +GRIDDER_EXE = BUILD_DIR / "parent_gridder" +TESTS_DIR = PROJECT_ROOT / "tests" +DATA_DIR = TESTS_DIR / "data" + + +def create_test_snapshot(filepath, redshift, boxsize=10.0): + """Create a minimal test snapshot at specified redshift""" + with h5py.File(filepath, 'w') as f: + # Header + header = f.create_group('Header') + header.attrs['Redshift'] = redshift + header.attrs['BoxSize'] = np.array([boxsize, boxsize, boxsize]) + header.attrs['NumPart_Total'] = np.array([0, 1, 0, 0, 0, 0], dtype=np.uint64) + + # Units (10^10 Msun, Mpc, km/s) + units = f.create_group('Units') + units.attrs['Unit mass in cgs (U_M)'] = 1.989e43 + units.attrs['Unit length in cgs (U_L)'] = 3.086e24 + units.attrs['Unit time in cgs (U_t)'] = 3.086e19 + + # Cells metadata + cells_meta = f.create_group('Cells/Meta-data') + cells_meta.attrs['dimension'] = np.array([2, 2, 2], dtype=np.int32) + cells_meta.attrs['size'] = np.array([boxsize/2, boxsize/2, boxsize/2]) + + # One particle at center + part1 = f.create_group('PartType1') + part1.create_dataset('Coordinates', data=np.array([[boxsize/2, boxsize/2, boxsize/2]])) + part1.create_dataset('Masses', data=np.array([1.0])) + + # Cell counts (one particle in central cell) + f.create_dataset('Cells/Counts/PartType1', data=np.array([0, 0, 0, 1, 0, 0, 0, 0], dtype=np.int32)) + f.create_dataset('Cells/OffsetsInFile/PartType1', data=np.array([0, 0, 0, 0, 1, 1, 1, 1], dtype=np.int32)) + + +def create_test_params(filepath, snapshot_path, cosmology): + """Create test parameter file with specified cosmology""" + params = { + 'Kernels': { + 'nkernels': 1, + 'kernel_radius_1': 1.0 + }, + 'Grid': { + 'type': 'uniform', + 'cdim': 2 + }, + 'Cosmology': cosmology, + 'Tree': { + 'max_leaf_count': 200 + }, + 'Input': { + 'filepath': str(snapshot_path), + 'placeholder': '0000' + }, + 'Output': { + 'filepath': str(filepath.parent) + '/', + 'basename': filepath.name, + 'write_masses': 0 + } + } + + with open(filepath, 'w') as f: + yaml.dump(params, f, default_flow_style=False) + + +def run_gridder_and_get_density(snapshot_path, params_path, output_path): + """Run gridder and extract mean density from output""" + result = subprocess.run( + [str(GRIDDER_EXE), str(params_path), "1"], + capture_output=True, + text=True, + timeout=30 + ) + + if result.returncode != 0: + raise RuntimeError(f"Gridder failed:\nstdout: {result.stdout}\nstderr: {result.stderr}") + + # Extract mean density from stdout + for line in result.stdout.split('\n'): + if 'Mean comoving density at z=' in line: + # Parse: "Mean comoving density at z=X.XXXX: Y.YYYYYYe+XX 10^10 Msun/cMpc^3" + parts = line.split(':') + if len(parts) >= 2: + density_str = parts[1].strip().split()[0] + return float(density_str) + + raise RuntimeError(f"Could not extract mean density from output:\n{result.stdout}") + + +def calculate_expected_density(h, Omega_cdm, Omega_b, redshift): + """ + Calculate expected mean COMOVING density using standard cosmological formula. + + ρ_comoving = ρ_crit(z=0) × Ω_m + + Note: In comoving coordinates, density does NOT evolve with redshift. + The (1+z)³ factor converts to physical density, but SWIFT uses comoving coords. + + where ρ_crit(z=0) = 3H₀²/(8πG) + + Units: 10^10 Msun/cMpc^3 + """ + # Physical constants in internal units + H0_kmsMpc = 100.0 * h # km/s/Mpc + G = 4.301744232015554e+01 # (10^10 Msun)^-1 Mpc (km/s)^2 + + # Critical density today + rho_crit_0 = (3.0 * H0_kmsMpc**2) / (8.0 * np.pi * G) + + # Total matter density parameter + Omega_m = Omega_cdm + Omega_b + + # Mean COMOVING density (constant with redshift!) + mean_density = rho_crit_0 * Omega_m + + return mean_density + + +class TestCosmologyCalculations: + """Test suite for cosmological mean density calculations""" + + @pytest.fixture(scope="class") + def test_workspace(self): + """Create temporary workspace for tests""" + with tempfile.TemporaryDirectory() as tmpdir: + workspace = Path(tmpdir) + yield workspace + + def test_planck2018_z0(self, test_workspace): + """Test Planck 2018 cosmology at z=0""" + # Planck 2018 parameters + cosmology = { + 'h': 0.6736, + 'Omega_cdm': 0.2607, + 'Omega_b': 0.0493 + } + redshift = 0.0 + + # Create test files + snapshot_path = test_workspace / "snapshot_z0.hdf5" + params_path = test_workspace / "params_z0.yml" + output_path = test_workspace / "output_z0.hdf5" + + create_test_snapshot(snapshot_path, redshift) + create_test_params(params_path, snapshot_path, cosmology) + + # Run gridder and extract density + gridder_density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + + # Calculate expected density + expected_density = calculate_expected_density( + cosmology['h'], cosmology['Omega_cdm'], cosmology['Omega_b'], redshift + ) + + # Check agreement (within 0.01% - account for floating point precision) + relative_error = abs(gridder_density - expected_density) / expected_density + assert relative_error < 1e-4, \ + f"Mean density mismatch at z={redshift}: gridder={gridder_density:.6e}, expected={expected_density:.6e}" + + def test_planck2018_z1(self, test_workspace): + """Test Planck 2018 cosmology at z=1""" + cosmology = { + 'h': 0.6736, + 'Omega_cdm': 0.2607, + 'Omega_b': 0.0493 + } + redshift = 1.0 + + snapshot_path = test_workspace / "snapshot_z1.hdf5" + params_path = test_workspace / "params_z1.yml" + output_path = test_workspace / "output_z1.hdf5" + + create_test_snapshot(snapshot_path, redshift) + create_test_params(params_path, snapshot_path, cosmology) + + gridder_density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + expected_density = calculate_expected_density( + cosmology['h'], cosmology['Omega_cdm'], cosmology['Omega_b'], redshift + ) + + # At z=1, density should be 8x higher than z=0 (since (1+z)³ = 8) + relative_error = abs(gridder_density - expected_density) / expected_density + assert relative_error < 1e-4, \ + f"Mean density mismatch at z={redshift}: gridder={gridder_density:.6e}, expected={expected_density:.6e}" + + def test_high_redshift_z7(self, test_workspace): + """Test at high redshift z=7 (relevant for JWST observations)""" + cosmology = { + 'h': 0.6736, + 'Omega_cdm': 0.2607, + 'Omega_b': 0.0493 + } + redshift = 7.0 + + snapshot_path = test_workspace / "snapshot_z7.hdf5" + params_path = test_workspace / "params_z7.yml" + output_path = test_workspace / "output_z7.hdf5" + + create_test_snapshot(snapshot_path, redshift) + create_test_params(params_path, snapshot_path, cosmology) + + gridder_density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + expected_density = calculate_expected_density( + cosmology['h'], cosmology['Omega_cdm'], cosmology['Omega_b'], redshift + ) + + # At z=7, density should be (1+7)³ = 512x higher than z=0 + relative_error = abs(gridder_density - expected_density) / expected_density + assert relative_error < 1e-4, \ + f"Mean density mismatch at z={redshift}: gridder={gridder_density:.6e}, expected={expected_density:.6e}" + + def test_wmap9_cosmology(self, test_workspace): + """Test WMAP9 cosmology (different from Planck)""" + # WMAP9 parameters + cosmology = { + 'h': 0.697, + 'Omega_cdm': 0.233, + 'Omega_b': 0.0463 + } + redshift = 2.0 + + snapshot_path = test_workspace / "snapshot_wmap9.hdf5" + params_path = test_workspace / "params_wmap9.yml" + output_path = test_workspace / "output_wmap9.hdf5" + + create_test_snapshot(snapshot_path, redshift) + create_test_params(params_path, snapshot_path, cosmology) + + gridder_density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + expected_density = calculate_expected_density( + cosmology['h'], cosmology['Omega_cdm'], cosmology['Omega_b'], redshift + ) + + relative_error = abs(gridder_density - expected_density) / expected_density + assert relative_error < 1e-4, \ + f"Mean density mismatch with WMAP9 cosmology: gridder={gridder_density:.6e}, expected={expected_density:.6e}" + + def test_low_omega_matter(self, test_workspace): + """Test with artificially low matter density""" + cosmology = { + 'h': 0.7, + 'Omega_cdm': 0.15, + 'Omega_b': 0.03 + } + redshift = 0.5 + + snapshot_path = test_workspace / "snapshot_low_om.hdf5" + params_path = test_workspace / "params_low_om.yml" + output_path = test_workspace / "output_low_om.hdf5" + + create_test_snapshot(snapshot_path, redshift) + create_test_params(params_path, snapshot_path, cosmology) + + gridder_density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + expected_density = calculate_expected_density( + cosmology['h'], cosmology['Omega_cdm'], cosmology['Omega_b'], redshift + ) + + relative_error = abs(gridder_density - expected_density) / expected_density + assert relative_error < 1e-4, \ + f"Mean density mismatch with low Omega_m: gridder={gridder_density:.6e}, expected={expected_density:.6e}" + + def test_high_h(self, test_workspace): + """Test with high Hubble constant""" + cosmology = { + 'h': 0.9, # Unrealistically high but good test + 'Omega_cdm': 0.25, + 'Omega_b': 0.05 + } + redshift = 1.5 + + snapshot_path = test_workspace / "snapshot_high_h.hdf5" + params_path = test_workspace / "params_high_h.yml" + output_path = test_workspace / "output_high_h.hdf5" + + create_test_snapshot(snapshot_path, redshift) + create_test_params(params_path, snapshot_path, cosmology) + + gridder_density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + expected_density = calculate_expected_density( + cosmology['h'], cosmology['Omega_cdm'], cosmology['Omega_b'], redshift + ) + + relative_error = abs(gridder_density - expected_density) / expected_density + assert relative_error < 1e-4, \ + f"Mean density mismatch with high h: gridder={gridder_density:.6e}, expected={expected_density:.6e}" + + def test_density_constant_with_redshift(self, test_workspace): + """Test that COMOVING density is constant with redshift""" + cosmology = { + 'h': 0.6736, + 'Omega_cdm': 0.2607, + 'Omega_b': 0.0493 + } + + redshifts = [0.0, 1.0, 2.0, 5.0, 7.0] + densities = [] + + for z in redshifts: + snapshot_path = test_workspace / f"snapshot_z{z}.hdf5" + params_path = test_workspace / f"params_z{z}.yml" + output_path = test_workspace / f"output_z{z}.hdf5" + + create_test_snapshot(snapshot_path, z) + create_test_params(params_path, snapshot_path, cosmology) + + density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + densities.append(density) + + # Check that comoving density is constant with redshift (within 0.01%) + reference_density = densities[0] + for i, z in enumerate(redshifts): + relative_diff = abs(densities[i] - reference_density) / reference_density + assert relative_diff < 1e-4, \ + f"Comoving density should be constant: rho(z={z})={densities[i]:.6e}, rho(z=0)={reference_density:.6e}, diff={relative_diff:.6e}" + + @pytest.mark.skipif(not ASTROPY_AVAILABLE, reason="astropy not available") + def test_astropy_comparison(self, test_workspace): + """Test that our calculation matches astropy exactly""" + # Planck 2018 parameters + cosmology = { + 'h': 0.6736, + 'Omega_cdm': 0.2607, + 'Omega_b': 0.0493 + } + redshift = 2.0 + + # Create test files + snapshot_path = test_workspace / "snapshot_astropy.hdf5" + params_path = test_workspace / "params_astropy.yml" + output_path = test_workspace / "output_astropy.hdf5" + + create_test_snapshot(snapshot_path, redshift) + create_test_params(params_path, snapshot_path, cosmology) + + # Get density from our C++ code + gridder_density = run_gridder_and_get_density(snapshot_path, params_path, output_path) + + # Calculate with astropy + h = cosmology['h'] + Omega_m = cosmology['Omega_cdm'] + cosmology['Omega_b'] + H0 = h * 100.0 # km/s/Mpc + + # Create astropy cosmology (flat ΛCDM with Omega_Lambda = 1 - Omega_m) + cosmo = FlatLambdaCDM(H0=H0, Om0=Omega_m) + + # Get critical density at z=0 in comoving coordinates + rho_crit_0 = cosmo.critical_density(0) + + # Convert to our units: 10^10 Msun/Mpc^3 + # astropy gives g/cm^3, we need 10^10 Msun/Mpc^3 + Msun_in_g = 1.989e33 + Mpc_in_cm = 3.086e24 + rho_crit_0_our_units = rho_crit_0.to(u.g / u.cm**3).value * (Mpc_in_cm**3) / (1e10 * Msun_in_g) + + # Mean comoving density + astropy_density = rho_crit_0_our_units * Omega_m + + # Check agreement (within 0.1% - account for different constants) + relative_error = abs(gridder_density - astropy_density) / astropy_density + assert relative_error < 1e-3, \ + f"Density mismatch vs astropy: gridder={gridder_density:.6e}, astropy={astropy_density:.6e}, error={relative_error:.6e}" + + +if __name__ == '__main__': + # Run tests with pytest + pytest.main([__file__, '-v', '--tb=short']) diff --git a/tests/test_file_grid params.yml b/tests/test_file_grid params.yml new file mode 100644 index 0000000..90dd7d3 --- /dev/null +++ b/tests/test_file_grid params.yml @@ -0,0 +1,29 @@ +# Test parameter file for file-based gridding +# Tests gridding with coordinates loaded from a file + +Kernels: + nkernels: 3 + kernel_radius_1: 0.5 # Small scale + kernel_radius_2: 2.0 # Medium scale + kernel_radius_3: 5.0 # Large scale + +Grid: + type: file + grid_file: tests/data/test_grid_coordinates.txt + +Cosmology: + h: 0.681 # Reduced Hubble constant + Omega_cdm: 0.256011 # Cold dark matter density parameter + Omega_b: 0.048600 # Baryon density parameter + + +Tree: + max_leaf_count: 200 + +Input: + filepath: tests/data/test_snapshot.hdf5 + +Output: + filepath: tests/data/ + basename: file_grid_output.hdf5 + write_masses: 1 diff --git a/tests/test_file_gridding.py b/tests/test_file_gridding.py new file mode 100755 index 0000000..33f3748 --- /dev/null +++ b/tests/test_file_gridding.py @@ -0,0 +1,555 @@ +#!/usr/bin/env python3 +""" +Comprehensive test suite for file-based gridding functionality. + +This tests the critical case where grid points are loaded from a file, +which is the workflow that's failing for halo center gridding. +""" + +import argparse +import h5py +import numpy as np +import os +import subprocess +import sys +import tempfile + + +class FileGriddingTest: + """Test suite for file-based gridding""" + + def __init__(self, executable, test_dir="tests/data"): + self.executable = executable + self.test_dir = test_dir + os.makedirs(test_dir, exist_ok=True) + + # Test configuration + self.box_size = 10.0 # Mpc/h + self.n_particles = 1000 + self.kernel_radii = [0.5, 2.0, 5.0] + + def create_test_snapshot(self, filename): + """Create a simple test snapshot with particles in clusters""" + print(f"Creating test snapshot: {filename}") + + with h5py.File(filename, 'w') as f: + # Header + header = f.create_group('Header') + # BoxSize must be an array [x, y, z] not a scalar! + header.attrs['BoxSize'] = np.array([self.box_size, self.box_size, self.box_size]) + header.attrs['NumPart_Total'] = [0, self.n_particles, 0, 0, 0, 0] + header.attrs['NumPart_ThisFile'] = [0, self.n_particles, 0, 0, 0, 0] + header.attrs['Redshift'] = 7.0 + + # Units + units = f.create_group('Units') + units.attrs['Unit mass in cgs (U_M)'] = 1.989e43 # 10^10 Msun + units.attrs['Unit length in cgs (U_L)'] = 3.086e24 # Mpc + + # Cosmology + cosmo = f.create_group('Cosmology') + cosmo.attrs['Critical density [internal units]'] = 1.0 + + # Cells metadata (required by gridder) + cells = f.create_group('Cells') + cells_meta = cells.create_group('Meta-data') + cells_meta.attrs['dimension'] = np.array([2, 2, 2], dtype=np.int32) + cells_meta.attrs['nr_cells'] = 8 + # Cell size must be an array [x, y, z] for each cell dimension, not a scalar! + # With dimension [2, 2, 2] and box size 10.0, each cell is 5.0 x 5.0 x 5.0 + cell_width = self.box_size / 2.0 + cells_meta.attrs['size'] = np.array([cell_width, cell_width, cell_width]) + + # Cells Centres + cell_width = self.box_size / 2.0 + centres = [] + for i in range(2): + for j in range(2): + for k in range(2): + centres.append([ + (i + 0.5) * cell_width, + (j + 0.5) * cell_width, + (k + 0.5) * cell_width + ]) + cells.create_dataset('Centres', data=np.array(centres)) + + # Create Counts and OffsetsInFile groups (not datasets!) + counts_grp = cells.create_group('Counts') + offsets_grp = cells.create_group('OffsetsInFile') + + # Will be filled after we know particle positions + counts_grp.create_dataset('PartType1', data=np.zeros(8, dtype=np.int32)) + offsets_grp.create_dataset('PartType1', data=np.zeros(8, dtype=np.int32)) + + # Create particle distributions: 5 clusters + n_clusters = 5 + particles_per_cluster = self.n_particles // n_clusters + cluster_centers = [ + [2.0, 2.0, 2.0], + [5.0, 5.0, 5.0], + [8.0, 2.0, 5.0], + [2.0, 8.0, 5.0], + [5.0, 5.0, 8.0], + ] + + positions = [] + masses = [] + + for center in cluster_centers: + # Gaussian distribution around center + cluster_pos = np.random.normal( + loc=center, + scale=0.3, # Tight cluster + size=(particles_per_cluster, 3) + ) + # Wrap to box + cluster_pos = np.mod(cluster_pos, self.box_size) + positions.append(cluster_pos) + masses.extend([1.0] * particles_per_cluster) + + # Add remaining particles uniformly + remainder = self.n_particles - (n_clusters * particles_per_cluster) + if remainder > 0: + uniform_pos = np.random.uniform( + 0, self.box_size, size=(remainder, 3) + ) + positions.append(uniform_pos) + masses.extend([1.0] * remainder) + + positions = np.vstack(positions) + masses = np.array(masses) + + # Compute cell counts for PartType1 + cell_counts = np.zeros(8, dtype=np.int32) + cell_offsets = np.zeros(8, dtype=np.int32) + + # Assign particles to cells + for i, pos in enumerate(positions): + ix = int(pos[0] / cell_width) + iy = int(pos[1] / cell_width) + iz = int(pos[2] / cell_width) + # Ensure within bounds + ix = max(0, min(1, ix)) + iy = max(0, min(1, iy)) + iz = max(0, min(1, iz)) + cell_id = ix * 4 + iy * 2 + iz + cell_counts[cell_id] += 1 + + # Set offsets (cumulative) + offset = 0 + for cell_id in range(8): + cell_offsets[cell_id] = offset + offset += cell_counts[cell_id] + + # Update cell datasets + f['Cells/Counts/PartType1'][:] = cell_counts + f['Cells/OffsetsInFile/PartType1'][:] = cell_offsets + + # Write particle data + part1 = f.create_group('PartType1') + part1.create_dataset('Coordinates', data=positions) + part1.create_dataset('Masses', data=masses) + part1.create_dataset('ParticleIDs', data=np.arange(len(masses))) + + print(f" Created {len(masses)} particles in {n_clusters} clusters") + return cluster_centers + + def create_grid_file(self, filename, grid_points): + """Create a text file with grid point coordinates""" + print(f"Creating grid file: {filename}") + + with open(filename, 'w') as f: + f.write("# Test grid point coordinates\n") + f.write("# Format: x y z (Mpc/h)\n") + f.write(f"# {len(grid_points)} grid points\n") + f.write("\n") + + for point in grid_points: + f.write(f"{point[0]:.6f} {point[1]:.6f} {point[2]:.6f}\n") + + print(f" Wrote {len(grid_points)} grid points") + + def create_param_file(self, filename, snapshot_file, grid_file, + output_file): + """Create a parameter file for the gridder""" + print(f"Creating parameter file: {filename}") + + content = f"""# Test parameter file for file-based gridding +Kernels: + nkernels: {len(self.kernel_radii)} +""" + for i, radius in enumerate(self.kernel_radii, 1): + content += f" kernel_radius_{i}: {radius}\n" + + content += f""" +Grid: + type: file + grid_file: {grid_file} + +Tree: + max_leaf_count: 200 + +Input: + filepath: {snapshot_file} + +Output: + filepath: {os.path.dirname(output_file)}/ + basename: {os.path.basename(output_file)} + write_masses: 1 +""" + + with open(filename, 'w') as f: + f.write(content) + + def run_gridder(self, param_file, nthreads=4, verbosity=2): + """Run the gridder executable""" + print(f"\nRunning gridder: {self.executable}") + print(f" Params: {param_file}") + print(f" Threads: {nthreads}") + print(f" Verbosity: {verbosity}") + + cmd = [self.executable, param_file, str(nthreads), str(verbosity)] + + try: + result = subprocess.run( + cmd, + capture_output=True, + text=True, + check=True + ) + print(" Gridder completed successfully") + + # Always print debug output (contains validation results) + if result.stdout: + print("\n--- Gridder Output ---") + print(result.stdout) + if result.stderr: + print("\n--- Stderr ---") + print(result.stderr) + print("--- End Output ---\n") + + # Check for debug validation failures in output + if '[DEBUG]' in result.stdout or '[DEBUG]' in result.stderr: + debug_output = result.stdout + result.stderr + if 'ERROR' in debug_output or 'FAILED' in debug_output: + print("\n!!! DEBUG CHECKS DETECTED ERRORS !!!") + print("See debug output above for details.\n") + + return True + except subprocess.CalledProcessError as e: + print(f" ERROR: Gridder failed with exit code {e.returncode}") + print("\n--- Stdout ---") + print(e.stdout) + print("\n--- Stderr ---") + print(e.stderr) + return False + + def validate_output(self, output_file, expected_grid_points): + """Validate the gridder output""" + print(f"\nValidating output: {output_file}") + + if not os.path.exists(output_file): + print(f" ERROR: Output file not found: {output_file}") + return False + + errors = [] + + with h5py.File(output_file, 'r') as f: + # Check basic structure + required_groups = ['Header', 'Grids', 'Cells'] + for group in required_groups: + if group not in f: + errors.append(f"Missing required group: {group}") + + if 'Grids' not in f: + return False + + # Check kernels (filter out non-kernel datasets like GridPointPositions) + all_keys = list(f['Grids'].keys()) + kernels = [k for k in all_keys if k.startswith('Kernel_')] + print(f" Found {len(kernels)} kernels: {kernels}") + if len(all_keys) != len(kernels): + print(f" (Ignoring non-kernel datasets: {[k for k in all_keys if not k.startswith('Kernel_')]})") + + if len(kernels) != len(self.kernel_radii): + errors.append( + f"Expected {len(self.kernel_radii)} kernels, " + f"found {len(kernels)}" + ) + + # Check each kernel + for i, expected_radius in enumerate(self.kernel_radii): + kernel_name = f'Kernel_{i}' + + if kernel_name not in f['Grids']: + errors.append(f"Missing kernel: {kernel_name}") + continue + + kernel_grp = f['Grids'][kernel_name] + + # Check kernel radius attribute + if 'KernelRadius' in kernel_grp.attrs: + actual_radius = kernel_grp.attrs['KernelRadius'] + if not np.isclose(actual_radius, expected_radius): + errors.append( + f"{kernel_name}: Expected radius {expected_radius}, " + f"got {actual_radius}" + ) + else: + errors.append( + f"{kernel_name}: Missing KernelRadius attribute" + ) + + # Check datasets + required_datasets = ['GridPointOverDensities', 'GridPointMasses'] + for dataset in required_datasets: + if dataset not in kernel_grp: + errors.append( + f"{kernel_name}: Missing dataset {dataset}" + ) + continue + + data = kernel_grp[dataset][:] + print(f" {kernel_name}/{dataset}: shape={data.shape}, " + f"min={data.min():.3f}, max={data.max():.3f}, " + f"mean={data.mean():.3f}") + + # Check for data quality issues + n_empty = np.sum(data == -1) + n_zero = np.sum(data == 0) + n_positive = np.sum(data > 0) + + print(f" Values: {n_positive} positive, " + f"{n_zero} zero, {n_empty} empty (-1)") + + # CRITICAL CHECK: For halo center gridding, we expect + # most grid points to find particles + if dataset == 'GridPointMasses': + # Check for empty (-1) values + if n_empty > len(data) * 0.5: + errors.append( + f"{kernel_name}: More than 50% of grid points " + f"have no particles ({n_empty}/{len(data)}). " + f"This suggests a major problem!" + ) + # CRITICAL: Check for zero values when we expect particles + # (cluster centers should have particles!) + if 'cluster' in output_file.lower() and n_zero > len(data) * 0.5: + errors.append( + f"CRITICAL BUG: {kernel_name}: Grid points at cluster " + f"centers found ZERO particles ({n_zero}/{len(data)}). " + f"Particles exist at these locations but gridder " + f"is not finding them!" + ) + + if errors: + print("\n VALIDATION ERRORS:") + for error in errors: + print(f" - {error}") + return False + + print(" ✓ Validation passed") + return True + + def test_cluster_centers(self): + """ + Test 1: Grid points at cluster centers + This is the critical test - grid points at known high-density regions + """ + print("\n" + "="*60) + print("TEST 1: Grid points at cluster centers") + print("="*60) + + # Create test files + snapshot_file = os.path.join(self.test_dir, "test_clusters.hdf5") + grid_file = os.path.join(self.test_dir, "cluster_centers.txt") + param_file = os.path.join(self.test_dir, "test_clusters.yml") + output_file = os.path.join(self.test_dir, "cluster_output.hdf5") + + # Create snapshot with known cluster centers + cluster_centers = self.create_test_snapshot(snapshot_file) + + # Create grid file with grid points AT cluster centers + # These points MUST find particles + self.create_grid_file(grid_file, cluster_centers) + + # Create parameter file + self.create_param_file(param_file, snapshot_file, grid_file, + output_file) + + # Run gridder + if not self.run_gridder(param_file): + print("\n✗ TEST 1 FAILED: Gridder execution failed") + return False + + # Validate output + if not self.validate_output(output_file, cluster_centers): + print("\n✗ TEST 1 FAILED: Output validation failed") + return False + + print("\n✓ TEST 1 PASSED: Cluster center gridding works") + return True + + def test_random_points(self): + """ + Test 2: Random grid points throughout box + """ + print("\n" + "="*60) + print("TEST 2: Random grid points") + print("="*60) + + snapshot_file = os.path.join(self.test_dir, "test_random.hdf5") + grid_file = os.path.join(self.test_dir, "random_points.txt") + param_file = os.path.join(self.test_dir, "test_random.yml") + output_file = os.path.join(self.test_dir, "random_output.hdf5") + + # Create snapshot + self.create_test_snapshot(snapshot_file) + + # Create random grid points + n_grid = 50 + grid_points = np.random.uniform(0, self.box_size, size=(n_grid, 3)) + self.create_grid_file(grid_file, grid_points) + + # Create parameter file + self.create_param_file(param_file, snapshot_file, grid_file, + output_file) + + # Run and validate + if not self.run_gridder(param_file): + print("\n✗ TEST 2 FAILED: Gridder execution failed") + return False + + if not self.validate_output(output_file, grid_points): + print("\n✗ TEST 2 FAILED: Output validation failed") + return False + + print("\n✓ TEST 2 PASSED: Random point gridding works") + return True + + def test_boundary_points(self): + """ + Test 3: Grid points near box boundaries (periodic BC test) + """ + print("\n" + "="*60) + print("TEST 3: Boundary and periodic BC test") + print("="*60) + + snapshot_file = os.path.join(self.test_dir, "test_boundary.hdf5") + grid_file = os.path.join(self.test_dir, "boundary_points.txt") + param_file = os.path.join(self.test_dir, "test_boundary.yml") + output_file = os.path.join(self.test_dir, "boundary_output.hdf5") + + # Create snapshot + self.create_test_snapshot(snapshot_file) + + # Create grid points at box edges and corners + eps = 0.1 # Small offset from boundary + grid_points = [ + # Corners + [eps, eps, eps], + [self.box_size - eps, eps, eps], + [eps, self.box_size - eps, eps], + [eps, eps, self.box_size - eps], + [self.box_size - eps, self.box_size - eps, eps], + [self.box_size - eps, eps, self.box_size - eps], + [eps, self.box_size - eps, self.box_size - eps], + [self.box_size - eps, self.box_size - eps, self.box_size - eps], + # Face centers + [self.box_size / 2, eps, self.box_size / 2], + [self.box_size / 2, self.box_size - eps, self.box_size / 2], + [eps, self.box_size / 2, self.box_size / 2], + [self.box_size - eps, self.box_size / 2, self.box_size / 2], + ] + self.create_grid_file(grid_file, grid_points) + + # Create parameter file + self.create_param_file(param_file, snapshot_file, grid_file, + output_file) + + # Run and validate + if not self.run_gridder(param_file): + print("\n✗ TEST 3 FAILED: Gridder execution failed") + return False + + if not self.validate_output(output_file, grid_points): + print("\n✗ TEST 3 FAILED: Output validation failed") + return False + + print("\n✓ TEST 3 PASSED: Boundary point gridding works") + return True + + def run_all_tests(self): + """Run all tests and report results""" + print("\n" + "="*60) + print("FILE-BASED GRIDDING TEST SUITE") + print("="*60) + print(f"Executable: {self.executable}") + print(f"Test directory: {self.test_dir}") + print("") + + tests = [ + ("Cluster Centers", self.test_cluster_centers), + ("Random Points", self.test_random_points), + ("Boundary Points", self.test_boundary_points), + ] + + results = {} + for name, test_func in tests: + try: + results[name] = test_func() + except Exception as e: + print(f"\n✗ TEST '{name}' CRASHED: {e}") + import traceback + traceback.print_exc() + results[name] = False + + # Print summary + print("\n" + "="*60) + print("TEST SUMMARY") + print("="*60) + + for name, passed in results.items(): + status = "✓ PASSED" if passed else "✗ FAILED" + print(f" {status}: {name}") + + n_passed = sum(results.values()) + n_total = len(results) + + print("") + print(f"Total: {n_passed}/{n_total} tests passed") + print("="*60) + + return all(results.values()) + + +def main(): + parser = argparse.ArgumentParser( + description="Test suite for file-based gridding functionality" + ) + parser.add_argument( + "executable", + help="Path to parent_gridder executable" + ) + parser.add_argument( + "--test-dir", + default="tests/data", + help="Directory for test files (default: tests/data)" + ) + + args = parser.parse_args() + + # Check executable exists + if not os.path.exists(args.executable): + print(f"ERROR: Executable not found: {args.executable}") + return 1 + + # Run tests + tester = FileGriddingTest(args.executable, args.test_dir) + success = tester.run_all_tests() + + return 0 if success else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/test_gridder.py b/tests/test_gridder.py index a11925c..6e13129 100644 --- a/tests/test_gridder.py +++ b/tests/test_gridder.py @@ -54,7 +54,7 @@ def test_snapshot(): DATA_DIR.mkdir(parents=True, exist_ok=True) # Create a minimal test snapshot using make_test_snap.py - make_snap_script = PROJECT_ROOT / "make_test_snap.py" + make_snap_script = TESTS_DIR / "make_test_snap.py" # Find a donor snapshot (use any existing HDF5 file or create minimal one) donor_path = DATA_DIR / "donor_minimal.hdf5" @@ -162,6 +162,11 @@ def test_grid_points_with_comments(self, build_gridder, test_snapshot): type: file grid_file: tests/data/grid_points_files/with_comments.txt +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -215,6 +220,11 @@ def test_malformed_grid_points(self, build_gridder, test_snapshot): type: file grid_file: tests/data/grid_points_files/malformed.txt +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -270,6 +280,11 @@ def test_missing_grid_file(self, build_gridder, test_snapshot): type: file grid_file: tests/data/grid_points_files/nonexistent.txt +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -405,6 +420,11 @@ def test_random_grid_reproducibility(self, build_gridder, test_snapshot): n_grid_points: 50 random_seed: {seed} +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -468,6 +488,11 @@ def test_random_grid_different_seeds(self, build_gridder, test_snapshot): n_grid_points: 50 random_seed: {seed} +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -533,6 +558,11 @@ def test_sparse_grid_chunk_preparation(self, build_gridder, test_snapshot): type: uniform cdim: 5 # 5^3 = 125 grid points +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -598,6 +628,11 @@ def test_cell_mass_summation(self, build_gridder, test_snapshot): type: uniform cdim: 3 +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -668,6 +703,11 @@ def test_uniform_grid_on_particles_unit_overdensity(self, test_snapshot): type: uniform cdim: 3 # 3^3 = 27 grid points +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 @@ -737,6 +777,11 @@ def test_sparse_grid_on_particles_unit_overdensity(self, test_snapshot): type: file grid_file: {grid_file} +Cosmology: + h: 0.681 + Omega_cdm: 0.256011 + Omega_b: 0.048600 + Tree: max_leaf_count: 200 diff --git a/tests/test_suite.py b/tests/test_suite.py index d187c4d..284781b 100755 --- a/tests/test_suite.py +++ b/tests/test_suite.py @@ -120,10 +120,29 @@ def create_dense_distribution(filename, npart=10000, boxsize=10.0): return npart @staticmethod - def _write_snapshot(filename, coords, masses, boxsize): - """Write HDF5 snapshot in SWIFT format""" + def _write_snapshot(filename, coords, masses, boxsize, h=0.681, Omega_cdm=0.256011, Omega_b=0.048600): + """Write HDF5 snapshot in SWIFT format + + If masses are uniform (all equal), they will be rescaled to match the + cosmological mean density. If masses vary, they are used as-is (assumed + to already represent the correct density field). + """ os.makedirs(os.path.dirname(filename), exist_ok=True) + # Calculate cosmological mean density + Omega_m = Omega_cdm + Omega_b + H0_kmsMpc = 100.0 * h + G = 4.301744232015554e+01 # (10^10 Msun)^-1 Mpc (km/s)^2 + rho_crit_0 = (3.0 * H0_kmsMpc**2) / (8.0 * np.pi * G) + mean_density = rho_crit_0 * Omega_m + + # For uniform mass distributions, rescale to match cosmological mean density + # This ensures uniform distributions give zero mean overdensity + if np.allclose(masses, masses[0]): # All masses equal + volume = boxsize**3 + total_mass = mean_density * volume + masses = np.full(len(masses), total_mass / len(masses)) + with h5py.File(filename, 'w') as f: # Header header = f.create_group('Header') @@ -207,6 +226,11 @@ def _create_param_file(self, snapshot_file, output_file, cdim=3, kernel_radii=[0 f.write(f" type: uniform\n") f.write(f" cdim: {cdim}\n") f.write(f"\n") + f.write(f"Cosmology:\n") + f.write(f" h: 0.681 # Reduced Hubble constant\n") + f.write(f" Omega_cdm: 0.256011 # Cold dark matter density parameter\n") + f.write(f" Omega_b: 0.048600 # Baryon density parameter\n") + f.write(f"\n") f.write(f"Tree:\n") f.write(f" max_leaf_count: 200\n") f.write(f"\n") @@ -243,6 +267,10 @@ def test_single_particle_center(self): boxsize = 10.0 TestDataGenerator.create_single_particle(snapshot, position, mass, boxsize) + # Read actual particle mass from snapshot (will be rescaled to match cosmology) + with h5py.File(snapshot, 'r') as f: + actual_mass = f['PartType1/Masses'][0] + # Create parameters kernel_radii = [0.5, 1.0, 2.0] param_file = self._create_param_file(snapshot, "single_particle_out.hdf5", @@ -275,8 +303,8 @@ def test_single_particle_center(self): # Particle at center should be within kernel of center grid point # for radii >= ~0.8 (half cell diagonal is ~0.866) if radius >= 0.8: - if masses[center_idx] != mass: - return False, f"Kernel {i} (r={radius}): Expected mass {mass} at center, got {masses[center_idx]}" + if not np.isclose(masses[center_idx], actual_mass): + return False, f"Kernel {i} (r={radius}): Expected mass {actual_mass} at center, got {masses[center_idx]}" return True, "Single particle test passed" @@ -563,6 +591,10 @@ def test_boundary_particles(self): # Write snapshot TestDataGenerator._write_snapshot(snapshot, positions, masses, boxsize) + # Read actual total mass from snapshot (will be rescaled to match cosmology) + with h5py.File(snapshot, 'r') as f: + expected_mass = np.sum(f['PartType1/Masses'][:]) + # Grid with points that should capture particles across boundaries param_file = self._create_param_file(snapshot, "boundary_parts_out.hdf5", cdim=3, kernel_radii=[2.0]) @@ -581,7 +613,6 @@ def test_boundary_particles(self): # All particles should be captured by at least one grid point # Check that the total mass is conserved (allowing small tolerance) total_mass_captured = np.sum(masses_out) - expected_mass = len(positions) # All particles have mass=1 if abs(total_mass_captured - expected_mass) > expected_mass * 0.1: return False, f"Mass not conserved: {total_mass_captured} vs {expected_mass}" @@ -757,6 +788,11 @@ def test_random_grid_mpi(self): 'type': 'random', 'n_grid_points': 500 }, + 'Cosmology': { + 'h': 0.681, + 'Omega_cdm': 0.256011, + 'Omega_b': 0.048600 + }, 'Tree': { 'max_leaf_count': 50 }, @@ -851,6 +887,11 @@ def _create_param_file(self, snapshot, basename, cdim=5, kernel_radii=[0.5, 1.0] 'type': 'uniform', 'cdim': cdim }, + 'Cosmology': { + 'h': 0.681, + 'Omega_cdm': 0.256011, + 'Omega_b': 0.048600 + }, 'Tree': { 'max_leaf_count': 50 },