Merge pull request #2116 from ebknudsen/mcpl_input

Add support for MCPL files

CMakeLists.txt

@@ -36,6 +36,7 @@ option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags"
 option(OPENMC_USE_DAGMC       "Enable support for DAGMC (CAD) geometry"              OFF)
 option(OPENMC_USE_LIBMESH     "Enable support for libMesh unstructured mesh tallies" OFF)
 option(OPENMC_USE_MPI         "Enable MPI"                                           OFF)
+option(OPENMC_USE_MCPL        "Enable MCPL"                                          OFF)
 
 # Warnings for deprecated options
 foreach(OLD_OPT IN ITEMS "openmp" "profile" "coverage" "dagmc" "libmesh")
@@ -158,6 +159,15 @@ if(${HDF5_VERSION} VERSION_GREATER_EQUAL 1.12.0)
   list(APPEND cxxflags -DH5Oget_info_by_idx_vers=1 -DH5O_info_t_vers=1)
 endif()
 
+#===============================================================================
+# MCPL
+#===============================================================================
+
+if (OPENMC_USE_MCPL)
+  find_package(MCPL REQUIRED)
+  message(STATUS "Found MCPL: ${MCPL_DIR} (found version \"${MCPL_VERSION}\")")
+endif()
+
 #===============================================================================
 # Set compile/link flags based on which compiler is being used
 #===============================================================================
@@ -323,6 +333,7 @@ list(APPEND libopenmc_SOURCES
   src/lattice.cpp
   src/material.cpp
   src/math_functions.cpp
+  src/mcpl_interface.cpp
   src/mesh.cpp
   src/message_passing.cpp
   src/mgxs.cpp
@@ -482,6 +493,11 @@ if (OPENMC_USE_MPI)
   target_link_libraries(libopenmc MPI::MPI_CXX)
 endif()
 
+if (OPENMC_USE_MCPL)
+  target_compile_definitions(libopenmc PUBLIC OPENMC_MCPL)
+  target_link_libraries(libopenmc MCPL::mcpl)
+endif()
+
 #===============================================================================
 # Log build info that this executable can report later
 #===============================================================================

cmake/OpenMCConfig.cmake.in

@@ -25,3 +25,7 @@ endif()
 if(@OPENMC_USE_MPI@)
   find_package(MPI REQUIRED)
 endif()
+
+if(@OPENMC_USE_MCPL@)
+  find_package(MCPL REQUIRED)
+endif()

docs/source/io_formats/settings.rst

@@ -732,6 +732,14 @@ attributes/sub-elements:
 
     *Default*: false
 
+  :mcpl:
+    If this element is set to "true", the source point file containing the
+    source bank will be written as an MCPL_ file name ``source.mcpl`` instead of
+    an HDF5 file. This option is only applicable if the ``<separate>`` element
+    is set to true.
+
+    *Default*: false
+
 ------------------------------
 ``<surf_source_read>`` Element
 ------------------------------
@@ -767,6 +775,16 @@ certain surfaces and write out the source bank in a separate file called
 
     *Default*: None
 
+  :mcpl:
+    An optional boolean which indicates if the banked particles should be
+    written to a file in the MCPL_-format instead of the native HDF5-based
+    format. If activated the output file name is changed to
+    ``surface_source.mcpl``.
+
+    *Default*: false
+
+    .. _MCPL: https://mctools.github.io/mcpl/mcpl.pdf
+
 ------------------------------
 ``<survival_biasing>`` Element
 ------------------------------

include/openmc/mcpl_interface.h

@@ -0,0 +1,36 @@
+#ifndef OPENMC_MCPL_INTERFACE_H
+#define OPENMC_MCPL_INTERFACE_H
+
+#include "openmc/particle_data.h"
+#include "openmc/vector.h"
+
+#include <string>
+
+namespace openmc {
+
+//==============================================================================
+// Constants
+//==============================================================================
+
+extern "C" const bool MCPL_ENABLED;
+
+//==============================================================================
+// Functions
+//==============================================================================
+
+//! Get a vector of source sites from an MCPL file
+//
+//! \param[in] path  Path to MCPL file
+//! \return  Vector of source sites
+vector<SourceSite> mcpl_source_sites(std::string path);
+
+//! Write an MCPL source file
+//
+//! \param[in] filename  Path to MCPL file
+//! \param[in] surf_source_bank  Whether to use the surface source bank
+void write_mcpl_source_point(
+  const char* filename, bool surf_source_bank = false);
+
+} // namespace openmc
+
+#endif // OPENMC_MCPL_INTERFACE_H

include/openmc/settings.h

@@ -47,7 +47,9 @@ extern "C" bool run_CE;            //!< run with continuous-energy data?
 extern bool source_latest;         //!< write latest source at each batch?
 extern bool source_separate;       //!< write source to separate file?
 extern bool source_write;          //!< write source in HDF5 files?
+extern bool source_mcpl_write;     //!< write source in mcpl files?
 extern bool surf_source_write;     //!< write surface source file?
+extern bool surf_mcpl_write;       //!< write surface mcpl file?
 extern bool surf_source_read;      //!< read surface source file?
 extern bool survival_biasing;      //!< use survival biasing?
 extern bool temperature_multipole; //!< use multipole data?

include/openmc/source.h

@@ -102,6 +102,7 @@ class FileSource : public Source {
 public:
   // Constructors
   explicit FileSource(std::string path);
+  explicit FileSource(const vector<SourceSite>& sites) : sites_ {sites} {}
 
   // Methods
   SourceSite sample(uint64_t* seed) const override;

openmc/lib/__init__.py

@@ -48,6 +48,9 @@ def _coord_levels():
 def _libmesh_enabled():
     return c_bool.in_dll(_dll, "LIBMESH_ENABLED").value
 
+def _mcpl_enabled():
+    return c_bool.in_dll(_dll, "MCPL_ENABLED").value
+
 from .error import *
 from .core import *
 from .nuclide import *

openmc/settings.py

@@ -157,6 +157,7 @@ class Settings:
         :separate: bool indicating whether the source should be written as a
                    separate file
         :write: bool indicating whether or not to write the source
+        :mcpl: bool indicating whether to write the source as an MCPL file
     statepoint : dict
         Options for writing state points. Acceptable keys are:
 
@@ -172,6 +173,7 @@ class Settings:
                    banked (int)
         :max_particles: Maximum number of particles to be banked on
                    surfaces per process (int)
+        :mcpl: Output in the form of an MCPL-file (bool)
     survival_biasing : bool
         Indicate whether survival biasing is to be used
     tabular_legendre : dict
@@ -622,6 +624,8 @@ class Settings:
                 cv.check_type('sourcepoint write', value, bool)
             elif key == 'overwrite':
                 cv.check_type('sourcepoint overwrite', value, bool)
+            elif key == 'mcpl':
+                cv.check_type('sourcepoint mcpl', value, bool)
             else:
                 raise ValueError(f"Unknown key '{key}' encountered when "
                                  "setting sourcepoint options.")
@@ -655,7 +659,7 @@ class Settings:
         cv.check_type('surface source writing options', surf_source_write, Mapping)
         for key, value in surf_source_write.items():
             cv.check_value('surface source writing key', key,
-                           ('surface_ids', 'max_particles'))
+                           ('surface_ids', 'max_particles', 'mcpl'))
             if key == 'surface_ids':
                 cv.check_type('surface ids for source banking', value,
                               Iterable, Integral)
@@ -667,6 +671,9 @@ class Settings:
                               value, Integral)
                 cv.check_greater_than('maximum particle banks on surfaces per process',
                                       value, 0)
+            elif key == 'mcpl':
+                cv.check_type('write to an MCPL-format file', value, bool)
+
         self._surf_source_write = surf_source_write
 
     @confidence_intervals.setter
@@ -1018,6 +1025,10 @@ class Settings:
                 subelement = ET.SubElement(element, "overwrite_latest")
                 subelement.text = str(self._sourcepoint['overwrite']).lower()
 
+            if 'mcpl' in self._sourcepoint:
+                subelement = ET.SubElement(element, "mcpl")
+                subelement.text = str(self._sourcepoint['mcpl']).lower()
+
     def _create_surf_source_read_subelement(self, root):
         if self._surf_source_read:
             element = ET.SubElement(root, "surf_source_read")
@@ -1035,6 +1046,9 @@ class Settings:
             if 'max_particles' in self._surf_source_write:
                 subelement = ET.SubElement(element, "max_particles")
                 subelement.text = str(self._surf_source_write['max_particles'])
+            if 'mcpl' in self._surf_source_write:
+                subelement = ET.SubElement(element, "mcpl")
+                subelement.text = str(self._surf_source_write['mcpl']).lower()
 
     def _create_confidence_intervals(self, root):
         if self._confidence_intervals is not None:
@@ -1315,10 +1329,10 @@ class Settings:
     def _sourcepoint_from_xml_element(self, root):
         elem = root.find('source_point')
         if elem is not None:
-            for key in ('separate', 'write', 'overwrite_latest', 'batches'):
+            for key in ('separate', 'write', 'overwrite_latest', 'batches', 'mcpl'):
                 value = get_text(elem, key)
                 if value is not None:
-                    if key in ('separate', 'write'):
+                    if key in ('separate', 'write', 'mcpl'):
                         value = value in ('true', '1')
                     elif key == 'overwrite_latest':
                         value = value in ('true', '1')
@@ -1337,13 +1351,15 @@ class Settings:
     def _surf_source_write_from_xml_element(self, root):
         elem = root.find('surf_source_write')
         if elem is not None:
-            for key in ('surface_ids', 'max_particles'):
+            for key in ('surface_ids', 'max_particles','mcpl'):
                 value = get_text(elem, key)
                 if value is not None:
                     if key == 'surface_ids':
                         value = [int(x) for x in value.split()]
                     elif key in ('max_particles'):
                         value = int(value)
+                    elif key == 'mcpl':
+                        value = value in ('true', '1')
                     self.surf_source_write[key] = value
 
     def _confidence_intervals_from_xml_element(self, root):

src/mcpl_interface.cpp

@@ -0,0 +1,247 @@
+#include "openmc/mcpl_interface.h"
+
+#include "openmc/bank.h"
+#include "openmc/error.h"
+#include "openmc/message_passing.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+#include "openmc/state_point.h"
+#include "openmc/vector.h"
+
+#include <fmt/core.h>
+
+#ifdef OPENMC_MCPL
+#include <mcpl.h>
+#endif
+
+namespace openmc {
+
+//==============================================================================
+// Constants
+//==============================================================================
+
+#ifdef OPENMC_MCPL
+const bool MCPL_ENABLED = true;
+#else
+const bool MCPL_ENABLED = false;
+#endif
+
+//==============================================================================
+// Functions
+//==============================================================================
+
+#ifdef OPENMC_MCPL
+SourceSite mcpl_particle_to_site(const mcpl_particle_t* particle)
+{
+  SourceSite site;
+
+  switch (particle->pdgcode) {
+  case 2112:
+    site.particle = ParticleType::neutron;
+    break;
+  case 22:
+    site.particle = ParticleType::photon;
+    break;
+  case 11:
+    site.particle = ParticleType::electron;
+    break;
+  case -11:
+    site.particle = ParticleType::positron;
+    break;
+  }
+
+  // Copy position and direction
+  site.r.x = particle->position[0];
+  site.r.y = particle->position[1];
+  site.r.z = particle->position[2];
+  site.u.x = particle->direction[0];
+  site.u.y = particle->direction[1];
+  site.u.z = particle->direction[2];
+
+  // MCPL stores kinetic energy in [MeV], time in [ms]
+  site.E = particle->ekin * 1e6;
+  site.time = particle->time * 1e-3;
+  site.wgt = particle->weight;
+
+  return site;
+}
+#endif
+
+//==============================================================================
+
+vector<SourceSite> mcpl_source_sites(std::string path)
+{
+  vector<SourceSite> sites;
+
+#ifdef OPENMC_MCPL
+  // Open MCPL file and determine number of particles
+  auto mcpl_file = mcpl_open_file(path.c_str());
+  size_t n_sites = mcpl_hdr_nparticles(mcpl_file);
+
+  for (int i = 0; i < n_sites; i++) {
+    // Extract particle from mcpl-file, checking if it is a neutron, photon,
+    // electron, or positron. Otherwise skip.
+    const mcpl_particle_t* particle;
+    int pdg = 0;
+    while (pdg != 2112 && pdg != 22 && pdg != 11 && pdg != -11) {
+      particle = mcpl_read(mcpl_file);
+      pdg = particle->pdgcode;
+    }
+
+    // Convert to source site and add to vector
+    sites.push_back(mcpl_particle_to_site(particle));
+  }
+
+  // Check that some sites were read
+  if (sites.empty()) {
+    fatal_error("MCPL file contained no neutron, photon, electron, or positron "
+                "source particles.");
+  }
+
+  mcpl_close_file(mcpl_file);
+#else
+  fatal_error(
+    "Your build of OpenMC does not support reading MCPL source files.");
+#endif
+
+  return sites;
+}
+
+//==============================================================================
+
+#ifdef OPENMC_MCPL
+void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
+{
+  int64_t dims_size = settings::n_particles;
+  int64_t count_size = simulation::work_per_rank;
+
+  // Set vectors for source bank and starting bank index of each process
+  vector<int64_t>* bank_index = &simulation::work_index;
+  vector<SourceSite>* source_bank = &simulation::source_bank;
+  vector<int64_t> surf_source_index_vector;
+  vector<SourceSite> surf_source_bank_vector;
+
+  if (surf_source_bank) {
+    surf_source_index_vector = calculate_surf_source_size();
+    dims_size = surf_source_index_vector[mpi::n_procs];
+    count_size = simulation::surf_source_bank.size();
+
+    bank_index = &surf_source_index_vector;
+
+    // Copy data in a SharedArray into a vector.
+    surf_source_bank_vector.resize(count_size);
+    surf_source_bank_vector.assign(simulation::surf_source_bank.data(),
+      simulation::surf_source_bank.data() + count_size);
+    source_bank = &surf_source_bank_vector;
+  }
+
+  if (mpi::master) {
+    // Particles are writeen to disk from the master node only
+
+    // Save source bank sites since the array is overwritten below
+#ifdef OPENMC_MPI
+    vector<SourceSite> temp_source {source_bank->begin(), source_bank->end()};
+#endif
+
+    // loop over the other nodes and receive data - then write those.
+    for (int i = 0; i < mpi::n_procs; ++i) {
+      // number of particles for node node i
+      size_t count[] {
+        static_cast<size_t>((*bank_index)[i + 1] - (*bank_index)[i])};
+
+#ifdef OPENMC_MPI
+      if (i > 0)
+        MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i,
+          mpi::intracomm, MPI_STATUS_IGNORE);
+#endif
+      // now write the source_bank data again.
+      for (const auto& site : *source_bank) {
+        // particle is now at the iterator
+        // write it to the mcpl-file
+        mcpl_particle_t p;
+        p.position[0] = site.r.x;
+        p.position[1] = site.r.y;
+        p.position[2] = site.r.z;
+
+        // mcpl requires that the direction vector is unit length
+        // which is also the case in openmc
+        p.direction[0] = site.u.x;
+        p.direction[1] = site.u.y;
+        p.direction[2] = site.u.z;
+
+        // MCPL stores kinetic energy in [MeV], time in [ms]
+        p.ekin = site.E * 1e-6;
+        p.time = site.time * 1e3;
+        p.weight = site.wgt;
+
+        switch (site.particle) {
+        case ParticleType::neutron:
+          p.pdgcode = 2112;
+          break;
+        case ParticleType::photon:
+          p.pdgcode = 22;
+          break;
+        case ParticleType::electron:
+          p.pdgcode = 11;
+          break;
+        case ParticleType::positron:
+          p.pdgcode = -11;
+          break;
+        }
+
+        mcpl_add_particle(file_id, &p);
+      }
+    }
+#ifdef OPENMC_MPI
+    // Restore state of source bank
+    std::copy(temp_source.begin(), temp_source.end(), source_bank->begin());
+#endif
+  } else {
+#ifdef OPENMC_MPI
+    MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank,
+      mpi::intracomm);
+#endif
+  }
+}
+#endif
+
+//==============================================================================
+
+void write_mcpl_source_point(const char* filename, bool surf_source_bank)
+{
+  std::string filename_;
+  if (filename) {
+    filename_ = filename;
+  } else {
+    // Determine width for zero padding
+    int w = std::to_string(settings::n_max_batches).size();
+
+    filename_ = fmt::format("{0}source.{1:0{2}}.mcpl", settings::path_output,
+      simulation::current_batch, w);
+  }
+
+#ifdef OPENMC_MCPL
+  mcpl_outfile_t file_id;
+
+  std::string line;
+  if (mpi::master) {
+    file_id = mcpl_create_outfile(filename_.c_str());
+    if (VERSION_DEV) {
+      line = fmt::format("OpenMC {0}.{1}.{2}-development", VERSION_MAJOR,
+        VERSION_MINOR, VERSION_RELEASE);
+    } else {
+      line = fmt::format(
+        "OpenMC {0}.{1}.{2}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
+    }
+    mcpl_hdr_set_srcname(file_id, line.c_str());
+  }
+
+  write_mcpl_source_bank(file_id, surf_source_bank);
+
+  if (mpi::master) {
+    mcpl_close_outfile(file_id);
+  }
+#endif
+}
+
+} // namespace openmc

src/settings.cpp

@@ -18,6 +18,7 @@
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
+#include "openmc/mcpl_interface.h"
 #include "openmc/mesh.h"
 #include "openmc/message_passing.h"
 #include "openmc/output.h"
@@ -60,7 +61,9 @@ bool run_CE {true};
 bool source_latest {false};
 bool source_separate {false};
 bool source_write {true};
+bool source_mcpl_write {false};
 bool surf_source_write {false};
+bool surf_mcpl_write {false};
 bool surf_source_read {false};
 bool survival_biasing {false};
 bool temperature_multipole {false};
@@ -429,7 +432,12 @@ void read_settings_xml()
   for (pugi::xml_node node : root.children("source")) {
     if (check_for_node(node, "file")) {
       auto path = get_node_value(node, "file", false, true);
-      model::external_sources.push_back(make_unique<FileSource>(path));
+      if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) {
+        auto sites = mcpl_source_sites(path);
+        model::external_sources.push_back(make_unique<FileSource>(sites));
+      } else {
+        model::external_sources.push_back(make_unique<FileSource>(path));
+      }
     } else if (check_for_node(node, "library")) {
       // Get shared library path and parameters
       auto path = get_node_value(node, "library", false, true);
@@ -647,6 +655,15 @@ void read_settings_xml()
     if (check_for_node(node_sp, "write")) {
       source_write = get_node_value_bool(node_sp, "write");
     }
+    if (check_for_node(node_sp, "mcpl")) {
+      source_mcpl_write = get_node_value_bool(node_sp, "mcpl");
+
+      // Make sure MCPL support is enabled
+      if (source_mcpl_write && !MCPL_ENABLED) {
+        fatal_error(
+          "Your build of OpenMC does not support writing MCPL source files.");
+      }
+    }
     if (check_for_node(node_sp, "overwrite_latest")) {
       source_latest = get_node_value_bool(node_sp, "overwrite_latest");
       source_separate = source_latest;
@@ -677,9 +694,18 @@ void read_settings_xml()
       max_surface_particles =
         std::stoll(get_node_value(node_ssw, "max_particles"));
     }
+    if (check_for_node(node_ssw, "mcpl")) {
+      surf_mcpl_write = get_node_value_bool(node_ssw, "mcpl");
+
+      // Make sure MCPL support is enabled
+      if (surf_mcpl_write && !MCPL_ENABLED) {
+        fatal_error("Your build of OpenMC does not support writing MCPL "
+                    "surface source files.");
+      }
+    }
   }
 
-  // If source is not seperate and is to be written out in the statepoint file,
+  // If source is not separate and is to be written out in the statepoint file,
   // make sure that the sourcepoint batch numbers are contained in the
   // statepoint list
   if (!source_separate) {

src/simulation.cpp

@@ -8,6 +8,7 @@
 #include "openmc/event.h"
 #include "openmc/geometry_aux.h"
 #include "openmc/material.h"
+#include "openmc/mcpl_interface.h"
 #include "openmc/message_passing.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
@@ -391,21 +392,35 @@ void finalize_batch()
     // Write out a separate source point if it's been specified for this batch
     if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
         settings::source_write && settings::source_separate) {
-      write_source_point(nullptr);
+      if (settings::source_mcpl_write) {
+        write_mcpl_source_point(nullptr);
+      } else {
+        write_source_point(nullptr);
+      }
     }
 
     // Write a continously-overwritten source point if requested.
     if (settings::source_latest) {
-      auto filename = settings::path_output + "source.h5";
-      write_source_point(filename.c_str());
+      if (settings::source_mcpl_write) {
+        auto filename = settings::path_output + "source.mcpl";
+        write_mcpl_source_point(filename.c_str());
+      } else {
+        auto filename = settings::path_output + "source.h5";
+        write_source_point(filename.c_str());
+      }
     }
   }
 
   // Write out surface source if requested.
   if (settings::surf_source_write &&
       simulation::current_batch == settings::n_batches) {
-    auto filename = settings::path_output + "surface_source.h5";
-    write_source_point(filename.c_str(), true);
+    if (settings::surf_mcpl_write) {
+      auto filename = settings::path_output + "surface_source.mcpl";
+      write_mcpl_source_point(filename.c_str(), true);
+    } else {
+      auto filename = settings::path_output + "surface_source.h5";
+      write_source_point(filename.c_str(), true);
+    }
   }
 }
 

tests/regression_tests/source_mcpl_file/geometry.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <!-- Sphere with radius 10 -->
+  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="vacuum"/>
+  <cell id="1" material="1" region="-1" />
+    
+</geometry>

tests/regression_tests/source_mcpl_file/materials.xml

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+</materials>

tests/regression_tests/source_mcpl_file/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+3.009416E-01 3.229999E-03

tests/regression_tests/source_mcpl_file/settings.xml

@@ -0,0 +1,15 @@
+<?xml version="1.0"?>
+<settings>
+  <state_point batches="10" />
+  <source_point mcpl="true" separate="true" />
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+</settings>

tests/regression_tests/source_mcpl_file/test.py

@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+import openmc.lib
+import pytest
+import glob
+import os
+
+from tests.testing_harness import *
+pytestmark = pytest.mark.skipif(
+    not openmc.lib._mcpl_enabled(),
+    reason="MCPL is not enabled.")
+
+settings1="""<?xml version="1.0"?>
+<settings>
+  <state_point batches="10" />
+  <source_point mcpl="true" separate="true" />
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+  <source>
+    <space type="box">
+      <parameters>-4 -4 -4  4  4  4</parameters>
+    </space>
+  </source>
+</settings>
+"""
+
+settings2 = """<?xml version="1.0"?>
+<settings>
+  <eigenvalue>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <particles>1000</particles>
+  </eigenvalue>
+  <source>
+    <file>source.10.{}</file>
+  </source>
+</settings>
+"""
+
+
+class SourceFileTestHarness(TestHarness):
+    def execute_test(self):
+        """Run OpenMC with the appropriate arguments and check the outputs."""
+        try:
+            self._run_openmc()
+            self._test_output_created()
+            self._run_openmc_restart()
+            results = self._get_results()
+            self._write_results(results)
+            self._compare_results()
+        finally:
+            self._cleanup()
+
+    def update_results(self):
+        """Update the results_true using the current version of OpenMC."""
+        try:
+            self._run_openmc()
+            self._test_output_created()
+            self._run_openmc_restart()
+            results = self._get_results()
+            self._write_results(results)
+            self._overwrite_results()
+        finally:
+            self._cleanup()
+
+    def _test_output_created(self):
+        """Make sure statepoint and source files have been created."""
+        statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))
+        assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \
+             'exist.'
+        assert statepoint[0].endswith('h5'), \
+             'Statepoint file is not a HDF5 file.'
+
+        source = glob.glob(os.path.join(os.getcwd(), 'source.10.mcpl*'))
+        assert len(source) == 1, 'Either multiple or no source files exist.'
+        assert source[0].endswith('mcpl') or source[0].endswith('mcpl.gz'), \
+             'Source file is not a MCPL file.'
+
+    def _run_openmc_restart(self):
+        # Get the name of the source file.
+        source = glob.glob(os.path.join(os.getcwd(), 'source.10.*'))
+
+        # Write the new settings.xml file.
+        with open('settings.xml','w') as fh:
+            fh.write(settings2.format(source[0].split('.')[-1]))
+
+        # Run OpenMC.
+        self._run_openmc()
+
+    def _cleanup(self):
+        TestHarness._cleanup(self)
+        output = glob.glob(os.path.join(os.getcwd(), 'source.*'))
+        with open('settings.xml','w') as fh:
+            fh.write(settings1)
+
+
+def test_source_file():
+    harness = SourceFileTestHarness('statepoint.10.h5')
+    harness.main()

tests/unit_tests/test_settings.py

@@ -17,7 +17,7 @@ def test_export_to_xml(run_in_tmpdir):
     s.output = {'summary': True, 'tallies': False, 'path': 'here'}
     s.verbosity = 7
     s.sourcepoint = {'batches': [50, 150, 500, 1000], 'separate': True,
-                     'write': True, 'overwrite': True}
+                     'write': True, 'overwrite': True, 'mcpl': True}
     s.statepoint = {'batches': [50, 150, 500, 1000]}
     s.surf_source_read = {'path': 'surface_source_1.h5'}
     s.surf_source_write = {'surface_ids': [2], 'max_particles': 200}
@@ -75,7 +75,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.output == {'summary': True, 'tallies': False, 'path': 'here'}
     assert s.verbosity == 7
     assert s.sourcepoint == {'batches': [50, 150, 500, 1000], 'separate': True,
-                             'write': True, 'overwrite': True}
+                             'write': True, 'overwrite': True, 'mcpl': True}
     assert s.statepoint == {'batches': [50, 150, 500, 1000]}
     assert s.surf_source_read == {'path': 'surface_source_1.h5'}
     assert s.surf_source_write == {'surface_ids': [2], 'max_particles': 200}

tools/ci/gha-install-mcpl.sh

@@ -0,0 +1,7 @@
+#!/bin/bash
+set -ex
+cd $HOME
+git clone https://github.com/mctools/mcpl
+cd mcpl
+mkdir build && cd build
+cmake .. && make 2>/dev/null && sudo make install

tools/ci/gha-install.py

@@ -25,8 +25,8 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False):
     os.mkdir('build')
     os.chdir('build')
 
-    # Build in debug mode by default
-    cmake_cmd = ['cmake', '-DCMAKE_BUILD_TYPE=Debug']
+    # Build in debug mode by default with support for MCPL
+    cmake_cmd = ['cmake', '-DCMAKE_BUILD_TYPE=Debug', '-DOPENMC_USE_MCPL=on']
 
     # Turn off OpenMP if specified
     if not omp:

tools/ci/gha-install.sh

@@ -27,6 +27,9 @@ if [[ $LIBMESH = 'y' ]]; then
     ./tools/ci/gha-install-libmesh.sh
 fi
 
+# Install MCPL
+./tools/ci/gha-install-mcpl.sh
+
 # For MPI configurations, make sure mpi4py and h5py are built against the
 # correct version of MPI
 if [[ $MPI == 'y' ]]; then