Merge pull request #1637 from ypark234/surf_src

Add surface source read/write capability

docs/source/io_formats/index.rst

@@ -44,6 +44,7 @@ Output Files
 
    statepoint
    source
+   surface_source
    summary
    depletion_results
    particle_restart

docs/source/io_formats/settings.rst

@@ -723,6 +723,42 @@ attributes/sub-elements:
 
     *Default*: false
 
+---------------------------
+``<surf_src_read>`` Element
+---------------------------
+
+The ``<surf_src_read>`` element specifies a surface source file for OpenMC to
+read source bank for initializing histories.
+This element has the following attributes/sub-elements:
+
+  :path:
+    Absolute or relative path to a surface source file to read in source bank.
+
+    *Default*: ``surface_source.h5`` in current working directory
+
+----------------------------
+``<surf_src_write>`` Element
+----------------------------
+
+The ``<surf_src_write>`` element triggers OpenMC to bank particles crossing
+certain surfaces and write out the source bank in a separate file called
+``surface_source.h5``.
+This element has the following attributes/sub-elements:
+
+  :surface_ids:
+    A list of integers separated by spaces indicating the unique IDs of surfaces
+    for which crossing particles will be banked.
+
+    *Default*: None
+
+  :max_particles:
+    An integer indicating the maximum number of particles to be banked on
+    specified surfaces per processor. The size of source bank in
+    ``surface_source.h5`` is limited to this value times the number of
+    processors.
+
+    *Default*: None
+
 ------------------------------
 ``<survival_biasing>`` Element
 ------------------------------

docs/source/io_formats/source.rst

@@ -8,6 +8,10 @@ Normally, source data is stored in a state point file. However, it is possible
 to request that the source be written separately, in which case the format used
 is that documented here.
 
+When surface source writing is triggered, a source file named
+``surface_source.h5`` is written with only the sources on specified surfaces,
+following the same format.
+
 **/**
 
 :Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
@@ -16,6 +20,7 @@ is that documented here.
 
            - **source_bank** (Compound type) -- Source bank information for each
              particle. The compound type has fields ``r``, ``u``, ``E``,
-             ``wgt``, ``delayed_group``, and ``particle``, which represent the
-             position, direction, energy, weight, delayed group, and particle
-             type (0=neutron, 1=photon, 2=electron, 3=positron), respectively.
+             ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``,
+             which represent the position, direction, energy, weight,
+             delayed group, surface ID, and particle type (0=neutron, 1=photon,
+             2=electron, 3=positron), respectively.

docs/source/io_formats/statepoint.rst

@@ -52,9 +52,10 @@ The current version of the statepoint file format is 17.0.
              sum-of-squares for each global tally.
            - **source_bank** (Compound type) -- Source bank information for each
              particle. The compound type has fields ``r``, ``u``, ``E``,
-             ``wgt``, ``delayed_group``, and ``particle``, which represent the
-             position, direction, energy, weight, delayed group, and particle
-             type (0=neutron, 1=photon, 2=electron, 3=positron), respectively.
+             ``wgt``, ``delayed_group``, ``surf_id``, and ``particle``, which
+             represent the position, direction, energy, weight, delayed group,
+             surface ID, and particle type (0=neutron, 1=photon, 2=electron,
+             3=positron), respectively.
              Only present when `run_mode` is 'eigenvalue'.
 
 **/tallies/**

include/openmc/bank.h

@@ -18,6 +18,8 @@ namespace simulation {
 
 extern std::vector<Particle::Bank> source_bank;
 
+extern SharedArray<Particle::Bank> surf_source_bank;
+
 extern SharedArray<Particle::Bank> fission_bank;
 
 extern std::vector<int64_t> progeny_per_particle;

include/openmc/particle.h

@@ -172,6 +172,7 @@ public:
     double E;
     double wgt;
     int delayed_group;
+    int surf_id;
     Type particle;
     int64_t parent_id;
     int64_t progeny_id;

include/openmc/settings.h

@@ -45,6 +45,8 @@ 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 surf_source_write;        //!< write surface source file?
+extern bool surf_source_read;         //!< read surface source file?
 extern bool survival_biasing;         //!< use survival biasing?
 extern bool temperature_multipole;    //!< use multipole data?
 extern "C" bool trigger_on;           //!< tally triggers enabled?
@@ -85,6 +87,8 @@ extern std::vector<std::string> res_scat_nuclides;  //!< Nuclides using res. ups
 extern RunMode run_mode;                 //!< Run mode (eigenvalue, fixed src, etc.)
 extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
 extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written
+extern std::unordered_set<int> source_write_surf_id; //!< Surface ids where sources will be written
+extern int64_t max_particles;            //!< maximum number of particles to be banked on surfaces per process
 extern TemperatureMethod temperature_method;           //!< method for choosing temperatures
 extern double temperature_tolerance;     //!< Tolerance in [K] on choosing temperatures
 extern double temperature_default;       //!< Default T in [K]

include/openmc/state_point.h

@@ -12,8 +12,9 @@
 namespace openmc {
 
 void load_state_point();
-void write_source_point(const char* filename);
-void write_source_bank(hid_t group_id);
+std::vector<int64_t> calculate_surf_source_size();
+void write_source_point(const char* filename, bool surf_source_bank = false);
+void write_source_bank(hid_t group_id, bool surf_source_bank);
 void read_source_bank(hid_t group_id, std::vector<Particle::Bank>& sites, bool distribute);
 void write_tally_results_nr(hid_t file_id);
 void restart_set_keff();

include/openmc/surface.h

@@ -88,6 +88,7 @@ public:
   int id_; //!< Unique ID
   std::string name_; //!< User-defined name
   std::shared_ptr<BoundaryCondition> bc_ {nullptr}; //!< Boundary condition
+  bool surf_source_ {false};     //!< Activate source banking for the surface?
 
   explicit Surface(pugi::xml_node surf_node);
   Surface();

openmc/lib/core.py

@@ -17,6 +17,7 @@ class _Bank(Structure):
                 ('E', c_double),
                 ('wgt', c_double),
                 ('delayed_group', c_int),
+                ('surf_id', c_int),
                 ('particle', c_int),
                 ('parent_id', c_int64),
                 ('progeny_id', c_int64)]

openmc/settings.py

@@ -143,6 +143,17 @@ class Settings:
         Options for writing state points. Acceptable keys are:
 
         :batches: list of batches at which to write source
+    surf_source_read : dict
+        Options for reading surface source points. Acceptable keys are:
+
+        :path: Path to surface source file (str).
+    surf_source_write : dict
+        Options for writing surface source points. Acceptable keys are:
+
+        :surface_ids: List of surface ids at which crossing particles are to be
+                   banked (int)
+        :max_particles: Maximum number of particles to be banked on surfaces
+                         per process (int)
     survival_biasing : bool
         Indicate whether survival biasing is to be used
     tabular_legendre : dict
@@ -229,6 +240,9 @@ class Settings:
         self._statepoint = {}
         self._sourcepoint = {}
 
+        self._surf_source_read = {}
+        self._surf_source_write = {}
+
         self._no_reduce = None
 
         self._verbosity = None
@@ -356,6 +370,14 @@ class Settings:
     def statepoint(self):
         return self._statepoint
 
+    @property
+    def surf_source_read(self):
+        return self._surf_source_read
+
+    @property
+    def surf_source_write(self):
+        return self._surf_source_write
+
     @property
     def no_reduce(self):
         return self._no_reduce
@@ -567,6 +589,35 @@ class Settings:
                                  "statepoint options.".format(key))
         self._statepoint = statepoint
 
+    @surf_source_read.setter
+    def surf_source_read(self, surf_source_read):
+        cv.check_type('surface source reading options', surf_source_read, Mapping)
+        for key, value in surf_source_read.items():
+            cv.check_value('surface source reading key', key,
+                           ('path'))
+            if key == 'path':
+                cv.check_type('path to surface source file', value, str)
+        self._surf_source_read = surf_source_read
+
+    @surf_source_write.setter
+    def surf_source_write(self, surf_source_write):
+        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'))
+            if key == 'surface_ids':
+                cv.check_type('surface ids for source banking', value,
+                              Iterable, Integral)
+                for surf_id in value:
+                    cv.check_greater_than('surface id for source banking',
+                                          surf_id, 0)
+            elif key == 'max_particles':
+                cv.check_type('maximum particle banks on surfaces per process',
+                              value, Integral)
+                cv.check_greater_than('maximum particle banks on surfaces per process',
+                                      value, 0)
+        self._surf_source_write = surf_source_write
+
     @confidence_intervals.setter
     def confidence_intervals(self, confidence_intervals):
         cv.check_type('confidence interval', confidence_intervals, bool)
@@ -889,6 +940,24 @@ class Settings:
                 subelement = ET.SubElement(element, "overwrite_latest")
                 subelement.text = str(self._sourcepoint['overwrite']).lower()
 
+    def _create_surf_source_read_subelement(self, root):
+        if self._surf_source_read:
+            element = ET.SubElement(root, "surf_source_read")
+            if 'path' in self._surf_source_read:
+                subelement = ET.SubElement(element, "path")
+                subelement.text = self._surf_source_read['path']
+
+    def _create_surf_source_write_subelement(self, root):
+        if self._surf_source_write:
+            element = ET.SubElement(root, "surf_source_write")
+            if 'surface_ids' in self._surf_source_write:
+                subelement = ET.SubElement(element, "surface_ids")
+                subelement.text = ' '.join(
+                    str(x) for x in self._surf_source_write['surface_ids'])
+            if 'max_particles' in self._surf_source_write:
+                subelement = ET.SubElement(element, "max_particles")
+                subelement.text = str(self._surf_source_write['max_particles'])
+
     def _create_confidence_intervals(self, root):
         if self._confidence_intervals is not None:
             element = ET.SubElement(root, "confidence_intervals")
@@ -1151,6 +1220,25 @@ class Settings:
                         value = [int(x) for x in value.split()]
                     self.sourcepoint[key] = value
 
+    def _surf_source_read_from_xml_element(self, root):
+        elem = root.find('surf_source_read')
+        if elem is not None:
+            value = get_text(elem, 'path')
+            if value is not None:
+                self.surf_source_read['path'] = value
+
+    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'):
+                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)
+                    self.surf_source_write[key] = value
+
     def _confidence_intervals_from_xml_element(self, root):
         text = get_text(root, 'confidence_intervals')
         if text is not None:
@@ -1349,6 +1437,8 @@ class Settings:
         self._create_output_subelement(root_element)
         self._create_statepoint_subelement(root_element)
         self._create_sourcepoint_subelement(root_element)
+        self._create_surf_source_read_subelement(root_element)
+        self._create_surf_source_write_subelement(root_element)
         self._create_confidence_intervals(root_element)
         self._create_electron_treatment_subelement(root_element)
         self._create_energy_mode_subelement(root_element)
@@ -1422,6 +1512,8 @@ class Settings:
         settings._output_from_xml_element(root)
         settings._statepoint_from_xml_element(root)
         settings._sourcepoint_from_xml_element(root)
+        settings._surf_source_read_from_xml_element(root)
+        settings._surf_source_write_from_xml_element(root)
         settings._confidence_intervals_from_xml_element(root)
         settings._electron_treatment_from_xml_element(root)
         settings._energy_mode_from_xml_element(root)

openmc/source.py

@@ -257,17 +257,20 @@ class SourceParticle:
         Weight of the particle
     delayed_group : int
         Delayed group particle was created in (neutrons only)
+    surf_id : int
+        Surface ID where particle is at, if any.
     particle : ParticleType
         Type of the particle
 
     """
     def __init__(self, r=(0., 0., 0.), u=(0., 0., 1.), E=1.0e6, wgt=1.0,
-                 delayed_group=0, particle=ParticleType.NEUTRON):
+                 delayed_group=0, surf_id=0, particle=ParticleType.NEUTRON):
         self.r = tuple(r)
         self.u = tuple(u)
         self.E = float(E)
         self.wgt = float(wgt)
         self.delayed_group = delayed_group
+        self.surf_id = surf_id
         self.particle = particle
 
     def to_tuple(self):
@@ -280,7 +283,7 @@ class SourceParticle:
 
         """
         return (self.r, self.u, self.E, self.wgt,
-                self.delayed_group, self.particle.value)
+                self.delayed_group, self.surf_id, self.particle.value)
 
 
 def write_source_file(source_particles, filename, **kwargs):
@@ -304,6 +307,7 @@ def write_source_file(source_particles, filename, **kwargs):
         ('E', '<f8'),
         ('wgt', '<f8'),
         ('delayed_group', '<i4'),
+        ('surf_id', '<i4'),
         ('particle', '<i4'),
     ])
 

openmc/statepoint.py

@@ -94,9 +94,9 @@ class StatePoint:
         Pseudorandom number generator seed
     source : numpy.ndarray of compound datatype
         Array of source sites. The compound datatype has fields 'r', 'u',
-        'E', 'wgt', 'delayed_group', and 'particle', corresponding to the
-        position, direction, energy, weight, delayed group, and particle type
-        of the source site, respectively.
+        'E', 'wgt', 'delayed_group', 'surf_id', and 'particle', corresponding to
+        the position, direction, energy, weight, delayed group, surface ID and
+        particle type of the source site, respectively.
     source_present : bool
         Indicate whether source sites are present
     sparse : bool

src/bank.cpp

@@ -17,6 +17,8 @@ namespace simulation {
 
 std::vector<Particle::Bank> source_bank;
 
+SharedArray<Particle::Bank> surf_source_bank;
+
 // The fission bank is allocated as a SharedArray, rather than a vector, as it will
 // be shared by all threads in the simulation. It will be allocated to a fixed
 // maximum capacity in the init_fission_bank() function. Then, Elements will be
@@ -37,6 +39,7 @@ std::vector<int64_t> progeny_per_particle;
 void free_memory_bank()
 {
   simulation::source_bank.clear();
+  simulation::surf_source_bank.clear();
   simulation::fission_bank.clear();
   simulation::progeny_per_particle.clear();
 }

src/dagmc.cpp

@@ -2,6 +2,7 @@
 
 #include "openmc/cell.h"
 #include "openmc/constants.h"
+#include "openmc/container_util.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
 #include "openmc/geometry.h"
@@ -296,6 +297,10 @@ void load_dagmc_geometry()
     s->id_ = model::DAG->id_by_index(2, s->dag_index_);
     s->dagmc_ptr_ = model::DAG;
 
+    if (contains(settings::source_write_surf_id, s->id_)) {
+      s->surf_source_ = true;
+    }
+
     // set BCs
     std::string bc_value = DMD.get_surface_property("boundary", surf_handle);
     to_lower(bc_value);

src/initialize.cpp

@@ -102,22 +102,23 @@ void initialize_mpi(MPI_Comm intracomm)
 
   // Create bank datatype
   Particle::Bank b;
-  MPI_Aint disp[8];
+  MPI_Aint disp[9];
   MPI_Get_address(&b.r, &disp[0]);
   MPI_Get_address(&b.u, &disp[1]);
   MPI_Get_address(&b.E, &disp[2]);
   MPI_Get_address(&b.wgt, &disp[3]);
   MPI_Get_address(&b.delayed_group, &disp[4]);
-  MPI_Get_address(&b.particle, &disp[5]);
-  MPI_Get_address(&b.parent_id, &disp[6]);
-  MPI_Get_address(&b.progeny_id, &disp[7]);
-  for (int i = 7; i >= 0; --i) {
+  MPI_Get_address(&b.surf_id, &disp[5]);
+  MPI_Get_address(&b.particle, &disp[6]);
+  MPI_Get_address(&b.parent_id, &disp[7]);
+  MPI_Get_address(&b.progeny_id, &disp[8]);
+  for (int i = 8; i >= 0; --i) {
     disp[i] -= disp[0];
   }
 
-  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1};
-  MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
-  MPI_Type_create_struct(8, blocks, disp, types, &mpi::bank);
+  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1};
+  MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
+  MPI_Type_create_struct(9, blocks, disp, types, &mpi::bank);
   MPI_Type_commit(&mpi::bank);
 }
 #endif // OPENMC_MPI

src/particle.cpp

@@ -413,6 +413,20 @@ Particle::cross_surface()
     write_message(1, "    Crossing surface {}", surf->id_);
   }
 
+  if (surf->surf_source_ && simulation::current_batch == settings::n_batches) {
+    Particle::Bank site;
+    site.r = this->r();
+    site.u = this->u();
+    site.E = this->E_;
+    site.wgt = this->wgt_;
+    site.delayed_group = this->delayed_group_;
+    site.surf_id = surf->id_;
+    site.particle = this->type_;
+    site.parent_id = this->id_;
+    site.progeny_id = this->n_progeny_;
+    int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
+  }
+
   // Handle any applicable boundary conditions.
   if (surf->bc_ && settings::run_mode != RunMode::PLOTTING) {
     surf->bc_->handle_particle(*this, *surf);

src/physics.cpp

@@ -190,6 +190,7 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
     site.wgt = 1. / weight;
     site.parent_id = p.id_;
     site.progeny_id = p.n_progeny_++;
+    site.surf_id = 0;
 
     // Sample delayed group and angle/energy for fission reaction
     sample_fission_neutron(i_nuclide, rx, p.E_, &site, p.current_seed());

src/relaxng/settings.rnc

@@ -139,6 +139,17 @@ element settings {
       attribute overwrite_latest {xsd:boolean})?
   }? &
 
+  element surf_source_read {
+    (element path { xsd:string } | attribute path { xsd:string })
+  }? &
+
+  element surf_source_write {
+    (element surface_ids { list { xsd:positiveInteger+ } } |
+      attribute surface_ids { list { xsd:positiveInteger+ } }) &
+    (element max_particles { xsd:positiveInteger } |
+      attribute max_particles { xsd:positiveInteger })
+  }? &
+
   element survival_biasing { xsd:boolean }? &
 
   element temperature_default { xsd:double }? &

src/relaxng/settings.rng

@@ -638,6 +638,48 @@
         </interleave>
       </element>
     </optional>
+    <optional>
+      <element name="surf_source_read">
+        <choice>
+          <element name="path">
+            <data type="string"/>
+          </element>
+          <attribute name="path">
+            <data type="string"/>
+          </attribute>
+        </choice>
+      </element>
+    </optional>
+    <optional>
+      <element name="surf_source_write">
+        <interleave>
+          <choice>
+            <element name="surface_ids">
+              <list>
+                <oneOrMore>
+                  <data type="positiveInteger"/>
+                </oneOrMore>
+              </list>
+            </element>
+            <attribute name="surface_ids">
+              <list>
+                <oneOrMore>
+                  <data type="positiveInteger"/>
+                </oneOrMore>
+              </list>
+            </attribute>
+          </choice>
+          <choice>
+            <element name="max_particles">
+              <data type="positiveInteger"/>
+            </element>
+            <attribute name="max_particles">
+              <data type="positiveInteger"/>
+            </attribute>
+          </choice>
+        </interleave>
+      </element>
+    </optional>
     <optional>
       <element name="survival_biasing">
         <data type="boolean"/>

src/settings.cpp

@@ -60,6 +60,8 @@ bool run_CE                  {true};
 bool source_latest           {false};
 bool source_separate         {false};
 bool source_write            {true};
+bool surf_source_write       {false};
+bool surf_source_read        {false};
 bool survival_biasing        {false};
 bool temperature_multipole   {false};
 bool trigger_on              {false};
@@ -98,6 +100,8 @@ std::vector<std::string> res_scat_nuclides;
 RunMode run_mode {RunMode::UNSET};
 std::unordered_set<int> sourcepoint_batch;
 std::unordered_set<int> statepoint_batch;
+std::unordered_set<int> source_write_surf_id;
+int64_t max_particles;
 TemperatureMethod temperature_method {TemperatureMethod::NEAREST};
 double temperature_tolerance {10.0};
 double temperature_default {293.6};
@@ -437,6 +441,20 @@ void read_settings_xml()
     }
   }
 
+  // Check if the user has specified to read surface source
+  if (check_for_node(root, "surf_source_read")) {
+    surf_source_read = true;
+    // Get surface source read node
+    xml_node node_ssr = root.child("surf_source_read");
+
+    std::string path = "surface_source.h5";
+    // Check if the user has specified different file for surface source reading
+    if (check_for_node(node_ssr, "path")) {
+      path = get_node_value(node_ssr, "path", false, true);
+    }
+    model::external_sources.push_back(std::make_unique<FileSource>(path));
+  }
+
   // If no source specified, default to isotropic point source at origin with Watt spectrum
   if (model::external_sources.empty()) {
     model::external_sources.push_back(std::make_unique<IndependentSource>(
@@ -622,6 +640,26 @@ void read_settings_xml()
     sourcepoint_batch = statepoint_batch;
   }
 
+  // Check if the user has specified to write surface source
+  if (check_for_node(root, "surf_source_write")) {
+    surf_source_write = true;
+    // Get surface source write node
+    xml_node node_ssw = root.child("surf_source_write");
+
+    // Determine surface ids at which crossing particles are to be banked
+    if (check_for_node(node_ssw, "surface_ids")) {
+      auto temp = get_node_array<int>(node_ssw, "surface_ids");
+      for (const auto& b : temp) {
+        source_write_surf_id.insert(b);
+      }
+    }
+
+    // Get maximum number of particles to be banked per surface
+    if (check_for_node(node_ssw, "max_particles")) {
+      max_particles = std::stoi(get_node_value(node_ssw, "max_particles"));
+    }
+  }
+
   // If source is not seperate and is to be written out in the statepoint file,
   // make sure that the sourcepoint batch numbers are contained in the
   // statepoint list
@@ -794,6 +832,7 @@ void read_settings_xml()
 void free_memory_settings() {
   settings::statepoint_batch.clear();
   settings::sourcepoint_batch.clear();
+  settings::source_write_surf_id.clear();
   settings::res_scat_nuclides.clear();
 }
 

src/simulation.cpp

@@ -79,10 +79,8 @@ int openmc_simulation_init()
   // Determine how much work each process should do
   calculate_work();
 
-  // Allocate source and fission banks for eigenvalue simulations
-  if (settings::run_mode == RunMode::EIGENVALUE) {
-    allocate_banks();
-  }
+  // Allocate source, fission and surface source banks.
+  allocate_banks();
 
   // If doing an event-based simulation, intialize the particle buffer
   // and event queues
@@ -281,10 +279,19 @@ std::vector<int64_t> work_index;
 
 void allocate_banks()
 {
-  // Allocate source bank
-  simulation::source_bank.resize(simulation::work_per_rank);
-  // Allocate fission bank
-  init_fission_bank(3*simulation::work_per_rank);
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    // Allocate source bank
+    simulation::source_bank.resize(simulation::work_per_rank);
+
+    // Allocate fission bank
+    init_fission_bank(3*simulation::work_per_rank);
+  }
+
+  if (settings::surf_source_write) {
+    // Allocate surface source bank
+    simulation::surf_source_bank.reserve(settings::max_particles);
+  }
+
 }
 
 void initialize_batch()
@@ -375,6 +382,12 @@ void finalize_batch()
       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);
+  }
 }
 
 void initialize_generation()
@@ -439,7 +452,10 @@ void finalize_generation()
 void initialize_history(Particle& p, int64_t index_source)
 {
   // set defaults
-  if (settings::run_mode == RunMode::FIXED_SOURCE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
+    // set defaults for eigenvalue simulations from primary bank
+    p.from_source(&simulation::source_bank[index_source - 1]);
+  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
     // initialize random number seed
     int64_t id = (simulation::total_gen + overall_generation() - 1)*settings::n_particles +
       simulation::work_index[mpi::rank] + index_source;
@@ -447,9 +463,6 @@ void initialize_history(Particle& p, int64_t index_source)
     // sample from external source distribution or custom library then set
     auto site = sample_external_source(&seed);
     p.from_source(&site);
-  } else if (settings::run_mode == RunMode::EIGENVALUE) {
-    // set defaults for eigenvalue simulations from primary bank
-    p.from_source(&simulation::source_bank[index_source - 1]);
   }
   p.current_work_ = index_source;
 

src/source.cpp

@@ -225,6 +225,8 @@ Particle::Bank IndependentSource::sample(uint64_t* seed) const
 
   // Set delayed group
   site.delayed_group = 0;
+  // Set surface ID
+  site.surf_id = 0;
 
   return site;
 }
@@ -342,7 +344,7 @@ void initialize_source()
     write_message("Writing out initial source...", 5);
     std::string filename = settings::path_output + "initial_source.h5";
     hid_t file_id = file_open(filename, 'w', true);
-    write_source_bank(file_id);
+    write_source_bank(file_id, false);
     file_close(file_id);
   }
 }

src/state_point.cpp

@@ -311,7 +311,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
   // Write the source bank if desired
   if (write_source_) {
     if (mpi::master || parallel) file_id = file_open(filename_, 'a', true);
-    write_source_bank(file_id);
+    write_source_bank(file_id, false);
     if (mpi::master || parallel) file_close(file_id);
   }
 
@@ -518,14 +518,38 @@ hid_t h5banktype() {
   H5Tinsert(banktype, "E", HOFFSET(Particle::Bank, E), H5T_NATIVE_DOUBLE);
   H5Tinsert(banktype, "wgt", HOFFSET(Particle::Bank, wgt), H5T_NATIVE_DOUBLE);
   H5Tinsert(banktype, "delayed_group", HOFFSET(Particle::Bank, delayed_group), H5T_NATIVE_INT);
+  H5Tinsert(banktype, "surf_id", HOFFSET(Particle::Bank, surf_id), H5T_NATIVE_INT);
   H5Tinsert(banktype, "particle", HOFFSET(Particle::Bank, particle), H5T_NATIVE_INT);
 
   H5Tclose(postype);
   return banktype;
 }
 
+std::vector<int64_t> calculate_surf_source_size()
+{
+  std::vector<int64_t> surf_source_index;
+  surf_source_index.reserve(mpi::n_procs + 1);
+
+#ifdef OPENMC_MPI
+  surf_source_index.resize(mpi::n_procs);
+  std::vector<int64_t> bank_size(mpi::n_procs);
+
+  // Populate the surf_source_index with cumulative sum of the number of
+  // surface source banks per process
+  int64_t size = simulation::surf_source_bank.size();
+  MPI_Scan(&size, bank_size.data(), 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
+  MPI_Allgather(bank_size.data(), 1, MPI_INT64_T, surf_source_index.data(), 1, MPI_INT64_T, mpi::intracomm);
+  surf_source_index.insert(surf_source_index.begin(), 0);
+#else
+  surf_source_index.push_back(0);
+  surf_source_index.push_back(simulation::surf_source_bank.size());
+#endif
+
+  return surf_source_index;
+}
+
 void
-write_source_point(const char* filename)
+write_source_point(const char* filename, bool surf_source_bank)
 {
   // When using parallel HDF5, the file is written to collectively by all
   // processes. With MPI-only, the file is opened and written by the master
@@ -555,29 +579,54 @@ write_source_point(const char* filename)
   }
 
   // Get pointer to source bank and write to file
-  write_source_bank(file_id);
+  write_source_bank(file_id, surf_source_bank);
 
   if (mpi::master || parallel) file_close(file_id);
 }
 
 void
-write_source_bank(hid_t group_id)
+write_source_bank(hid_t group_id, bool surf_source_bank)
 {
   hid_t banktype = h5banktype();
 
+  // Set total and individual process dataspace sizes for 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
+  std::vector<int64_t>* bank_index = &simulation::work_index;
+  std::vector<Particle::Bank>* source_bank = &simulation::source_bank;
+  std::vector<int64_t> surf_source_index_vector;
+  std::vector<Particle::Bank> surf_source_bank_vector;
+
+  // Reset dataspace sizes and vectors for surface source bank
+  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;
+  }
+
 #ifdef PHDF5
   // Set size of total dataspace for all procs and rank
-  hsize_t dims[] {static_cast<hsize_t>(settings::n_particles)};
+  hsize_t dims[] {static_cast<hsize_t>(dims_size)};
   hid_t dspace = H5Screate_simple(1, dims, nullptr);
   hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace,
                          H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
   // Create another data space but for each proc individually
-  hsize_t count[] {static_cast<hsize_t>(simulation::work_per_rank)};
+  hsize_t count[] {static_cast<hsize_t>(count_size)};
   hid_t memspace = H5Screate_simple(1, count, nullptr);
 
   // Select hyperslab for this dataspace
-  hsize_t start[] {static_cast<hsize_t>(simulation::work_index[mpi::rank])};
+  hsize_t start[] {static_cast<hsize_t>((*bank_index)[mpi::rank])};
   H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
   // Set up the property list for parallel writing
@@ -585,7 +634,7 @@ write_source_bank(hid_t group_id)
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
 
   // Write data to file in parallel
-  H5Dwrite(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
+  H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank->data());
 
   // Free resources
   H5Sclose(dspace);
@@ -597,38 +646,35 @@ write_source_bank(hid_t group_id)
 
   if (mpi::master) {
     // Create dataset big enough to hold all source sites
-    hsize_t dims[] {static_cast<hsize_t>(settings::n_particles)};
+    hsize_t dims[] {static_cast<hsize_t>(dims_size)};
     hid_t dspace = H5Screate_simple(1, dims, nullptr);
     hid_t dset = H5Dcreate(group_id, "source_bank", banktype, dspace,
                            H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
-    // Save source bank sites since the souce_bank array is overwritten below
+    // Save source bank sites since the array is overwritten below
 #ifdef OPENMC_MPI
-    std::vector<Particle::Bank> temp_source {simulation::source_bank.begin(),
-      simulation::source_bank.begin() + simulation::work_per_rank};
+    std::vector<Particle::Bank> temp_source {source_bank->begin(), source_bank->end()};
 #endif
 
     for (int i = 0; i < mpi::n_procs; ++i) {
       // Create memory space
-      hsize_t count[] {static_cast<hsize_t>(simulation::work_index[i+1] -
-        simulation::work_index[i])};
+      hsize_t count[] {static_cast<hsize_t>((*bank_index)[i+1] - (*bank_index)[i])};
       hid_t memspace = H5Screate_simple(1, count, nullptr);
 
 #ifdef OPENMC_MPI
       // Receive source sites from other processes
       if (i > 0)
-        MPI_Recv(simulation::source_bank.data(), count[0], mpi::bank, i, i,
+        MPI_Recv(source_bank->data(), count[0], mpi::bank, i, i,
                  mpi::intracomm, MPI_STATUS_IGNORE);
 #endif
 
       // Select hyperslab for this dataspace
       dspace = H5Dget_space(dset);
-      hsize_t start[] {static_cast<hsize_t>(simulation::work_index[i])};
+      hsize_t start[] {static_cast<hsize_t>((*bank_index)[i])};
       H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
       // Write data to hyperslab
-      H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT,
-        simulation::source_bank.data());
+      H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, (*source_bank).data());
 
       H5Sclose(memspace);
       H5Sclose(dspace);
@@ -639,11 +685,11 @@ write_source_bank(hid_t group_id)
 
 #ifdef OPENMC_MPI
     // Restore state of source bank
-    std::copy(temp_source.begin(), temp_source.end(), simulation::source_bank.begin());
+    std::copy(temp_source.begin(), temp_source.end(), source_bank->begin());
 #endif
   } else {
 #ifdef OPENMC_MPI
-    MPI_Send(simulation::source_bank.data(), simulation::work_per_rank, mpi::bank,
+    MPI_Send(source_bank->data(), count_size, mpi::bank,
       0, mpi::rank, mpi::intracomm);
 #endif
   }

src/surface.cpp

@@ -8,6 +8,7 @@
 #include <fmt/core.h>
 #include <gsl/gsl>
 
+#include "openmc/container_util.h"
 #include "openmc/error.h"
 #include "openmc/dagmc.h"
 #include "openmc/hdf5_interface.h"
@@ -115,6 +116,9 @@ Surface::Surface(pugi::xml_node surf_node)
 {
   if (check_for_node(surf_node, "id")) {
     id_ = std::stoi(get_node_value(surf_node, "id"));
+    if (contains(settings::source_write_surf_id, id_)) {
+      surf_source_ = true;
+    }
   } else {
     fatal_error("Must specify id of surface in geometry XML file.");
   }

tests/regression_tests/surface_source/inputs_true_read.dat

@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="void" region="-1" universe="1" />
+  <cell id="2" material="void" region="1 -2" universe="1" />
+  <cell id="3" material="void" region="2 -3" universe="1" />
+  <cell id="4" material="void" region="3 -4" universe="1" />
+  <surface coeffs="0.0 0.0 0.0 1.0" id="1" type="sphere" />
+  <surface coeffs="0.0 0.0 0.0 2.0" id="2" type="sphere" />
+  <surface coeffs="0.0 0.0 0.0 2.5" id="3" type="sphere" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 0.0 4.0" id="4" type="sphere" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>1000</particles>
+  <batches>10</batches>
+  <surf_source_read>
+    <path>surface_source_true.h5</path>
+  </surf_source_read>
+  <seed>1</seed>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="cell">
+    <bins>3</bins>
+  </filter>
+  <tally id="1">
+    <filters>1</filters>
+    <scores>flux</scores>
+  </tally>
+</tallies>

tests/regression_tests/surface_source/inputs_true_write.dat

@@ -0,0 +1,40 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="void" region="-1" universe="1" />
+  <cell id="2" material="void" region="1 -2" universe="1" />
+  <cell id="3" material="void" region="2 -3" universe="1" />
+  <cell id="4" material="void" region="3 -4" universe="1" />
+  <surface coeffs="0.0 0.0 0.0 1.0" id="1" type="sphere" />
+  <surface coeffs="0.0 0.0 0.0 2.0" id="2" type="sphere" />
+  <surface coeffs="0.0 0.0 0.0 2.5" id="3" type="sphere" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 0.0 4.0" id="4" type="sphere" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>1000</particles>
+  <batches>10</batches>
+  <source strength="1.0">
+    <space type="point">
+      <parameters>0 0 0</parameters>
+    </space>
+  </source>
+  <surf_source_write>
+    <surface_ids>1</surface_ids>
+    <max_particles>1000</max_particles>
+  </surf_source_write>
+  <seed>1</seed>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="cell">
+    <bins>3</bins>
+  </filter>
+  <tally id="1">
+    <filters>1</filters>
+    <scores>flux</scores>
+  </tally>
+</tallies>

tests/regression_tests/surface_source/results_true.dat

@@ -0,0 +1,3 @@
+tally 1:
+5.000000E+00
+2.500000E+00

tests/regression_tests/surface_source/test.py

@@ -0,0 +1,125 @@
+import os
+
+import h5py
+import numpy as np
+import pytest
+import openmc
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model(request):
+    openmc.reset_auto_ids()
+    marker = request.node.get_closest_marker("surf_source_op")
+    surf_source_op = marker.args[0]
+
+    openmc_model = openmc.model.Model()
+
+    # Materials
+    # None
+
+    # Geometry
+    # Concentric void spheres
+    # - Innermost sphere to bank surface sources
+    # - Second shell to tally cell flux
+    # - Outermost sphere as vacuum boundary
+    sph_1 = openmc.Sphere(r=1.0)  # Surface to bank/write sources.
+    sph_2 = openmc.Sphere(r=2.0)
+    sph_3 = openmc.Sphere(r=2.5)
+    sph_4 = openmc.Sphere(r=4.0, boundary_type='vacuum')
+    cell_1 = openmc.Cell(region=-sph_1)
+    cell_2 = openmc.Cell(region=+sph_1&-sph_2)
+    cell_3 = openmc.Cell(region=+sph_2&-sph_3)  # Cell to tally flux.
+    cell_4 = openmc.Cell(region=+sph_3&-sph_4)
+    root = openmc.Universe(cells=[cell_1, cell_2, cell_3, cell_4])
+    openmc_model.geometry = openmc.Geometry(root)
+
+    # Settings
+    openmc_model.settings.run_mode = 'fixed source'
+    openmc_model.settings.particles = 1000
+    openmc_model.settings.batches = 10
+    openmc_model.settings.seed = 1
+
+    if surf_source_op == 'write':
+        point = openmc.stats.Point((0, 0, 0))
+        pt_src = openmc.Source(space=point)
+        openmc_model.settings.source = pt_src
+
+        openmc_model.settings.surf_source_write = {'surface_ids': [1],
+                                                   'max_particles': 1000}
+    elif surf_source_op == 'read':
+        openmc_model.settings.surf_source_read = {'path': 'surface_source_true.h5'}
+
+    # Tallies
+    tal = openmc.Tally()
+    cell_filter = openmc.CellFilter(cell_3)
+    tal.filters = [cell_filter]
+    tal.scores = ['flux']
+    openmc_model.tallies.append(tal)
+
+    return openmc_model
+
+
+class SurfaceSourceTestHarness(PyAPITestHarness):
+    def _test_output_created(self):
+        """Make sure surface_source.h5 has also been created."""
+        super()._test_output_created()
+        # Check if 'surface_source.h5' has been created.
+        if self._model.settings.surf_source_write:
+            assert os.path.exists('surface_source.h5'), \
+                'Surface source file does not exist.'
+
+    def _compare_output(self):
+        """Make sure the current surface_source.h5 agree with the reference."""
+        if self._model.settings.surf_source_write:
+            with h5py.File("surface_source_true.h5", 'r') as f:
+                source_true = f['source_bank'][()]
+                # Convert dtye from mixed to a float for comparison assertion
+                source_true.dtype = 'float64'
+            with h5py.File("surface_source.h5", 'r') as f:
+                source_test = f['source_bank'][()]
+                # Convert dtye from mixed to a float for comparison assertion
+                source_test.dtype = 'float64'
+            np.testing.assert_allclose(np.sort(source_true),
+                                       np.sort(source_test),
+                                       atol=1e-07)
+
+    def execute_test(self):
+        """Build input XMLs, run OpenMC, check output and results."""
+        try:
+            self._build_inputs()
+            inputs = self._get_inputs()
+            self._write_inputs(inputs)
+            self._compare_inputs()
+            self._run_openmc()
+            self._test_output_created()
+            self._compare_output()
+            results = self._get_results()
+            self._write_results(results)
+            self._compare_results()
+        finally:
+            self._cleanup()
+
+    def _cleanup(self):
+        """Delete statepoints, tally, and test files."""
+        super()._cleanup()
+        fs = 'surface_source.h5'
+        if os.path.exists(fs):
+            os.remove(fs)
+
+
+@pytest.mark.surf_source_op('write')
+def test_surface_source_write(model):
+    harness = SurfaceSourceTestHarness('statepoint.10.h5',
+                                       model,
+                                       'inputs_true_write.dat')
+    harness.main()
+
+
+@pytest.mark.surf_source_op('read')
+def test_surface_source_read(model):
+    harness = SurfaceSourceTestHarness('statepoint.10.h5',
+                                       model,
+                                       'inputs_true_read.dat')
+    harness.main()

tests/unit_tests/test_settings.py

@@ -20,6 +20,8 @@ def test_export_to_xml(run_in_tmpdir):
     s.sourcepoint = {'batches': [50, 150, 500, 1000], 'separate': True,
                      'write': True, 'overwrite': 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}
     s.confidence_intervals = True
     s.ptables = True
     s.seed = 17
@@ -76,6 +78,8 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.sourcepoint == {'batches': [50, 150, 500, 1000], 'separate': True,
                              'write': True, 'overwrite': 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}
     assert s.confidence_intervals
     assert s.ptables
     assert s.seed == 17