paulromano · Aug 8, 2025 · Aug 6, 2025 · Aug 6, 2025 · Aug 6, 2025 · Aug 6, 2025
diff --git a/include/openmc/mcpl_interface.h b/include/openmc/mcpl_interface.h
@@ -19,8 +19,17 @@ namespace openmc {
 //! \return  Vector of source sites
 vector<SourceSite> mcpl_source_sites(std::string path);
 
-//! Write an MCPL source file
-//
+//! Write an MCPL source file with stat:sum metadata
+//!
+//! This function writes particle data to an MCPL file. For MCPL >= 2.1.0,
+//! it includes a stat:sum field (key: "openmc_np1") containing the total
+//! number of source particles, which is essential for proper file merging
+//! and weight normalization when using MCPL files with McStas/McXtrace.
+//!
+//! The stat:sum field follows the crash-safety pattern:
+//! - Initially set to -1 when opening (indicates incomplete file)
+//! - Updated with actual particle count before closing
+//!
 //! \param[in] filename     Path to MCPL file
 //! \param[in] source_bank  Vector of SourceSites to write to file for this
 //!                         MPI rank.

diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp
@@ -58,7 +58,8 @@ struct mcpl_outfile_t {
 // Function pointer types for the dynamically loaded MCPL library
 using mcpl_open_file_fpt = mcpl_file_t* (*)(const char* filename);
 using mcpl_hdr_nparticles_fpt = uint64_t (*)(mcpl_file_t* file_handle);
-using mcpl_read_fpt = const mcpl_particle_repr_t* (*)(mcpl_file_t* file_handle);
+using mcpl_read_fpt = const mcpl_particle_repr_t* (*)(mcpl_file_t *
+                                                      file_handle);
 using mcpl_close_file_fpt = void (*)(mcpl_file_t* file_handle);
 
 using mcpl_create_outfile_fpt = mcpl_outfile_t* (*)(const char* filename);
@@ -67,6 +68,10 @@ using mcpl_hdr_set_srcname_fpt = void (*)(
 using mcpl_add_particle_fpt = void (*)(
   mcpl_outfile_t* outfile_handle, const mcpl_particle_repr_t* particle);
 using mcpl_close_outfile_fpt = void (*)(mcpl_outfile_t* outfile_handle);
+using mcpl_hdr_add_data_fpt = void (*)(mcpl_outfile_t* outfile_handle,
+  const char* key, uint32_t datalength, const char* data);
+using mcpl_hdr_add_stat_sum_fpt = void (*)(
+  mcpl_outfile_t* outfile_handle, const char* key, double value);
 
 namespace openmc {
 
@@ -110,6 +115,8 @@ struct McplApi {
   mcpl_hdr_set_srcname_fpt hdr_set_srcname;
   mcpl_add_particle_fpt add_particle;
   mcpl_close_outfile_fpt close_outfile;
+  mcpl_hdr_add_data_fpt hdr_add_data;
+  mcpl_hdr_add_stat_sum_fpt hdr_add_stat_sum;
 
   explicit McplApi(LibraryHandleType lib_handle)
   {
@@ -147,6 +154,24 @@ struct McplApi {
       load_symbol_platform("mcpl_add_particle"));
     close_outfile = reinterpret_cast<mcpl_close_outfile_fpt>(
       load_symbol_platform("mcpl_close_outfile"));
+
+    // Try to load mcpl_hdr_add_data (available in MCPL >= 2.1.0)
+    // Set to nullptr if not available for graceful fallback
+    try {
+      hdr_add_data = reinterpret_cast<mcpl_hdr_add_data_fpt>(
+        load_symbol_platform("mcpl_hdr_add_data"));
+    } catch (const std::runtime_error&) {
+      hdr_add_data = nullptr;
+    }
+
+    // Try to load mcpl_hdr_add_stat_sum (available in MCPL >= 2.1.0)
+    // Set to nullptr if not available for graceful fallback
+    try {
+      hdr_add_stat_sum = reinterpret_cast<mcpl_hdr_add_stat_sum_fpt>(
+        load_symbol_platform("mcpl_hdr_add_stat_sum"));
+    } catch (const std::runtime_error&) {
+      hdr_add_stat_sum = nullptr;
+    }
   }
 };
 
@@ -498,12 +523,31 @@ void write_mcpl_source_point(const char* filename, span<SourceSite> source_bank,
         "OpenMC {}.{}.{}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
     }
     g_mcpl_api->hdr_set_srcname(file_id, src_line.c_str());
+
+    // Initialize stat:sum with -1 to indicate incomplete file (issue #3514)
+    // This follows MCPL >= 2.1.0 convention for tracking simulation statistics
+    // The -1 value indicates "not available" if file creation is interrupted
+    if (g_mcpl_api->hdr_add_stat_sum) {
+      // Using key "openmc_np1" following tkittel's recommendation
+      // Initial value of -1 prevents misleading values in case of crashes
+      g_mcpl_api->hdr_add_stat_sum(file_id, "openmc_np1", -1.0);
+    }
   }
 
   write_mcpl_source_bank_internal(file_id, source_bank, bank_index);
 
   if (mpi::master) {
     if (file_id) {
+      // Update stat:sum with actual particle count before closing (issue #3514)
+      // This represents the original number of source particles in the
+      // simulation
+      if (g_mcpl_api->hdr_add_stat_sum) {
+        int64_t total_particles = bank_index.empty() ? 0 : bank_index.back();
+        // Update with actual count - this overwrites the initial -1 value
+        g_mcpl_api->hdr_add_stat_sum(
+          file_id, "openmc_np1", static_cast<double>(total_particles));
+      }
+
       g_mcpl_api->close_outfile(file_id);
     }
   }

diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt
@@ -4,6 +4,7 @@ set(TEST_NAMES
   test_tally
   test_interpolate
   test_math
+  test_mcpl_stat_sum
   # Add additional unit test files here
 )
 

diff --git a/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp b/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp
@@ -0,0 +1,108 @@
+#include <catch2/catch_test_macros.hpp>
+#include <cstdio>
+#include <string>
+#include <vector>
+
+#include "openmc/bank.h"
+#include "openmc/mcpl_interface.h"
+
+// Test the MCPL stat:sum functionality (issue #3514)
+TEST_CASE("MCPL stat:sum field")
+{
+  // Check if MCPL interface is available
+  if (!openmc::is_mcpl_interface_available()) {
+    SKIP("MCPL library not available");
+  }
+
+  SECTION("stat:sum field is written to MCPL files")
+  {
+    // Create a temporary filename
+    std::string filename = "test_stat_sum.mcpl";
+
+    // Create some test particles
+    std::vector<openmc::SourceSite> source_bank(100);
+    std::vector<int64_t> bank_index = {0, 100}; // 100 particles total
+
+    // Initialize test particles
+    for (int i = 0; i < 100; ++i) {
+      source_bank[i].particle = openmc::ParticleType::neutron;
+      source_bank[i].r = {i * 0.1, i * 0.2, i * 0.3};
+      source_bank[i].u = {0.0, 0.0, 1.0};
+      source_bank[i].E = 2.0e6; // 2 MeV
+      source_bank[i].time = 0.0;
+      source_bank[i].wgt = 1.0;
+    }
+
+    // Write the MCPL file
+    openmc::write_mcpl_source_point(filename.c_str(), source_bank, bank_index);
+
+    // Verify the file was created
+    FILE* f = std::fopen(filename.c_str(), "r");
+    REQUIRE(f != nullptr);
+    std::fclose(f);
+
+    // Read the file back to check stat:sum
+    // Note: This would require mcpl_open_file and checking the header
+    // Since we can't easily read MCPL headers in C++ without the full MCPL API,
+    // we rely on the Python test to verify the actual content
+
+    // Clean up
+    std::remove(filename.c_str());
+  }
+
+  SECTION("stat:sum uses correct particle count")
+  {
+    std::string filename = "test_count.mcpl";
+
+    // Test with different particle counts
+    std::vector<int> test_counts = {1, 10, 100, 1000};
+
+    for (int count : test_counts) {
+      std::vector<openmc::SourceSite> source_bank(count);
+      std::vector<int64_t> bank_index = {0, count};
+
+      // Initialize particles
+      for (int i = 0; i < count; ++i) {
+        source_bank[i].particle = openmc::ParticleType::neutron;
+        source_bank[i].r = {0.0, 0.0, 0.0};
+        source_bank[i].u = {0.0, 0.0, 1.0};
+        source_bank[i].E = 1.0e6;
+        source_bank[i].time = 0.0;
+        source_bank[i].wgt = 1.0;
+      }
+
+      // Write MCPL file
+      openmc::write_mcpl_source_point(
+        filename.c_str(), source_bank, bank_index);
+
+      // The stat:sum should equal count (verified by Python test)
+      // Here we just verify the file was created successfully
+      FILE* f = std::fopen(filename.c_str(), "r");
+      REQUIRE(f != nullptr);
+      std::fclose(f);
+
+      // Clean up
+      std::remove(filename.c_str());
+    }
+  }
+
+  SECTION("stat:sum handles empty particle bank")
+  {
+    std::string filename = "test_empty.mcpl";
+
+    // Create empty particle bank
+    std::vector<openmc::SourceSite> source_bank;
+    std::vector<int64_t> bank_index = {0};
+
+    // This should still create a valid MCPL file with stat:sum = 0
+    openmc::write_mcpl_source_point(filename.c_str(), source_bank, bank_index);
+
+    // Verify file was created
+    FILE* f = std::fopen(filename.c_str(), "r");
+    REQUIRE(f != nullptr);
+    std::fclose(f);
+
+    // Clean up
+    std::remove(filename.c_str());
+  }
+}
diff --git a/tests/unit_tests/test_mcpl_stat_sum.py b/tests/unit_tests/test_mcpl_stat_sum.py
@@ -0,0 +1,115 @@
+#!/usr/bin/env python
+"""Test for MCPL stat:sum functionality (issue #3514)"""
+
+import os
+import shutil
+import tempfile
+import pytest
+import openmc
+
+# Skip test if MCPL is not available
+pytestmark = pytest.mark.skipif(
+    shutil.which("mcpl-config") is None,
+    reason="mcpl-config command not found in PATH; MCPL is likely not available."
+)
+
+
+def test_mcpl_stat_sum_field():
+    """Test that MCPL files contain proper stat:sum field with particle count."""
+
+    # Only run if mcpl module is available for verification
+    mcpl = pytest.importorskip("mcpl")
+
+    with tempfile.TemporaryDirectory() as tmpdir:
+        # Create a simple model
+        model = openmc.Model()
+
+        # Create a simple material
+        mat = openmc.Material()
+        mat.add_nuclide('U235', 1.0)
+        mat.set_density('g/cm3', 10.0)
+        model.materials = [mat]
+
+        # Create a simple geometry (sphere)
+        sphere = openmc.Sphere(r=10.0, boundary_type='vacuum')
+        cell = openmc.Cell(fill=mat, region=-sphere)
+        model.geometry = openmc.Geometry([cell])
+
+        # Configure settings for MCPL output
+        model.settings = openmc.Settings()
+        model.settings.batches = 5
+        model.settings.inactive = 2
+        model.settings.particles = 1000
+        model.settings.source = openmc.Source(space=openmc.stats.Point())
+
+        # Enable MCPL source file writing
+        model.settings.source_point = {'mcpl': True, 'separate': True}
+
+        # Run the simulation
+        cwd = os.getcwd()
+        try:
+            os.chdir(tmpdir)
+            model.run(output=False)
+
+            # Find the MCPL file
+            import glob
+            mcpl_files = glob.glob('source.*.mcpl*')
+            assert len(mcpl_files) > 0, "No MCPL source files were created"
+
+            # Open and check the MCPL file
+            mcpl_file = mcpl_files[0]
+            with mcpl.MCPLFile(mcpl_file) as f:
+                # Check if stat:sum field exists
+                # The stat:sum should be accessible through the header
+                # Look for the "openmc_np1" key
+                comments = f.comments if hasattr(f, 'comments') else []
+
+                # Check for stat:sum in comments (MCPL stores these as comments)
+                stat_sum_found = False
+                stat_sum_value = None
+
+                for comment in comments:
+                    if 'stat:sum:openmc_np1' in comment:
+                        stat_sum_found = True
+                        # Extract the value
+                        parts = comment.split(':')
+                        if len(parts) >= 4:
+                            stat_sum_value = float(parts[3].strip())
+                        break
+
+                assert stat_sum_found, "stat:sum:openmc_np1 field not found in MCPL file"
+
+                # The value should be the total number of source particles
+                # For 5 batches with 1000 particles each = 5000 total
+                expected_particles = model.settings.batches * model.settings.particles
+
+                # If stat:sum was properly updated, it should equal expected_particles
+                # If it's still -1, the update before closing didn't work
+                assert stat_sum_value != -1, "stat:sum was not updated from initial -1 value"
+                assert stat_sum_value == expected_particles, \
+                    f"stat:sum value {stat_sum_value} doesn't match expected {expected_particles}"
+
+        finally:
+            os.chdir(cwd)
+
+
+def test_mcpl_stat_sum_crash_safety():
+    """Test that incomplete MCPL files have stat:sum = -1."""
+
+    # This test would verify that if file creation is interrupted,
+    # the stat:sum field remains at -1 to indicate incomplete file
+    # This is harder to test without mocking the MCPL library internals
+
+    # For now, we can at least document the expected behavior:
+    # 1. When mcpl_create_outfile is called, stat:sum should be set to -1
+    # 2. Only when mcpl_close_outfile is called should it be updated
+    # 3. If the program crashes between these calls, stat:sum remains -1
+
+    # This could be tested with a C++ unit test that directly uses the
+    # mcpl_interface functions and simulates a crash
+    pass
+
+
+if __name__ == "__main__":
+    # Allow running this test directly
+    test_mcpl_stat_sum_field()