Fixing treatment of cosmology and mean density calculation

CMakeLists.txt

@@ -175,6 +175,7 @@ set(SOURCE_FILES
   src/partition.cpp
   src/params.cpp
   src/simulation.cpp
+  src/debugging_utils.cpp
 )
 
 # ========== Target ==========

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:

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

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<Cell> &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];

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<int>(cid))
+            continue;
+
           // Attach the neighbour to the cell
           cell->neighbours.push_back(&cells[cjd]);
         }

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) {