Add stat:sum field to MCPL files for proper weight normalization (#3522)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: bpolania

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.

src/mcpl_interface.cpp

@@ -67,6 +67,8 @@ 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_stat_sum_fpt = void (*)(
+  mcpl_outfile_t* outfile_handle, const char* key, double value);
 
 namespace openmc {
 
@@ -110,6 +112,7 @@ 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_stat_sum_fpt hdr_add_stat_sum;
 
   explicit McplApi(LibraryHandleType lib_handle)
   {
@@ -147,6 +150,15 @@ 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_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 +510,36 @@ 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 (not the number of particles in the file)
+      if (g_mcpl_api->hdr_add_stat_sum) {
+        // Calculate total source particles from active batches
+        // Per issue #3514: this should be the original number of source
+        // particles, not the number written to the file
+        int64_t total_source_particles =
+          static_cast<int64_t>(settings::n_batches - settings::n_inactive) *
+          settings::gen_per_batch * settings::n_particles;
+        // 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_source_particles));
+      }
+
       g_mcpl_api->close_outfile(file_id);
     }
   }

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
 )
 

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());
+  }
+}

tests/unit_tests/test_mcpl_stat_sum.py

@@ -0,0 +1,69 @@
+"""Test for MCPL stat:sum functionality"""
+
+from pathlib import Path
+import shutil
+
+import pytest
+import openmc
+
+
+@pytest.mark.skipif(shutil.which("mcpl-config") is None, reason="MCPL is not available.")
+def test_mcpl_stat_sum_field(run_in_tmpdir):
+    """Test that MCPL files contain proper stat:sum field with particle count.
+
+    This test verifies that when OpenMC creates MCPL source files, they contain
+    the stat:sum field. Since MCPL functions are not exposed in the Python API,
+    this test creates an actual OpenMC simulation to generate MCPL files and
+    then checks their content.
+    """
+
+    mcpl = pytest.importorskip("mcpl")
+
+    # Create a minimal working model that will generate MCPL files
+    model = openmc.examples.pwr_pin_cell()
+    model.settings.batches = 5
+    model.settings.inactive = 2
+    model.settings.particles = 1000
+    model.settings.sourcepoint = {'mcpl': True, 'separate': True}
+
+    # Run a short simulation to generate MCPL files
+    model.run(output=False)
+
+    # Find the generated MCPL file (from the last batch)
+    mcpl_file = Path('source.5.mcpl')
+    assert mcpl_file.exists(), "No MCPL files were generated"
+
+    # Open and verify the stat:sum field exists
+    with mcpl.MCPLFile(mcpl_file) as f:
+        # Check if stat:sum field exists using convenience property
+        if hasattr(f, 'stat_sum'):
+            # Use the convenience .stat_sum property directly
+            stat_sum_dict = f.stat_sum
+            assert 'openmc_np1' in stat_sum_dict, "openmc_np1 key not found in stat_sum"
+            stat_sum_value = int(stat_sum_dict['openmc_np1'])
+        else:
+            # Fallback to checking comments for older MCPL versions
+            comments = f.comments
+
+            # Check for stat:sum in comments (MCPL stores these as comments)
+            stat_sum_value = None
+
+            for comment in comments:
+                if 'stat:sum:openmc_np1' in comment:
+                    # Extract the value
+                    parts = comment.split(':')
+                    if len(parts) >= 4:
+                        stat_sum_value = int(parts[3].strip())
+                    break
+            else:
+                pytest.skip("stat:sum field not found - may be running with MCPL < 2.1.0")
+
+    # Verify the stat:sum value is reasonable
+    assert stat_sum_value != -1, "stat:sum was not updated from initial -1 value"
+
+    # In eigenvalue mode, active batches generate source particles
+    active_batches = model.settings.batches - model.settings.inactive  # 3 active batches
+    expected_particles = active_batches * model.settings.particles     # 3000 total
+
+    assert stat_sum_value == expected_particles, \
+        f"stat:sum value {stat_sum_value} doesn't match expected {expected_particles}"