Mark multigroup source files and keep their groups in FileSource

docs/source/io_formats/source.rst

@@ -17,6 +17,18 @@ following the same format.
 :Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
              - **version** (*int[2]*) -- Major and minor version of the source
                file format.
+             - **energy_mode** (*char[]*) -- Whether ``E`` holds energy group
+               indices ('multi-group') or energies in [eV]
+               ('continuous-energy'). A file without this attribute is read as
+               'continuous-energy'. When the file is read by
+               :class:`openmc.FileSource`, any other value is an error, and so
+               is a 'multi-group' file in a continuous-energy simulation.
+               OpenMC writes this attribute, with the value 'multi-group', only
+               to source point files written by a multigroup simulation.
+               Surface source files and files written by a continuous-energy
+               simulation, by :func:`openmc.write_source_file` or by
+               :meth:`openmc.ParticleList.export_to_hdf5` do not have this
+               attribute, so their ``E`` is read as an energy in [eV].
 
 :Datasets:
 
@@ -25,4 +37,6 @@ following the same format.
              ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``,
              which represent the position, direction, energy, time, weight,
              delayed group, surface ID, and particle type (PDG number),
-             respectively.
+             respectively. When the **energy_mode** attribute is
+             'multi-group', ``E`` is the energy group index of the particle
+             (0 for the highest-energy group) rather than its energy in [eV].

docs/source/io_formats/statepoint.rst

@@ -56,8 +56,10 @@ The current version of the statepoint file format is 18.2.
              ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and
              ``particle``, which represent the position, direction, energy,
              time, weight, delayed group, surface ID, and particle type
-             (PDG number), respectively. Only present when `run_mode` is
-             'eigenvalue'.
+             (PDG number), respectively. When `energy_mode` is 'multi-group',
+             ``E`` is the energy group index of the particle (0 for the
+             highest-energy group) rather than its energy in [eV]. Only present
+             when `run_mode` is 'eigenvalue'.
 
 **/tallies/**
 

docs/source/usersguide/settings.rst

@@ -352,6 +352,16 @@ source file can be manually generated with the :func:`openmc.write_source_file`
 function. This is particularly useful for coupling OpenMC with another program
 that generates a source to be used in OpenMC.
 
+In a multigroup simulation, the source point files and state points store the
+energy group index of each site rather than its energy, and source point files
+record this with an ``energy_mode`` attribute (see :ref:`io_source`). When a
+multigroup simulation reads such a file, each particle starts in its stored
+group; energy bounds given as source constraints are not applied to these sites,
+and a warning says so. A continuous-energy simulation cannot use such a file.
+Files without this record, such as surface source files and files written by
+:func:`openmc.write_source_file`, are read as energies in [eV], which a
+multigroup simulation converts to groups.
+
 Surface Sources
 +++++++++++++++
 

include/openmc/source.h

@@ -94,6 +94,10 @@ class Source {
 
   bool satisfies_spatial_constraints(Position r) const;
 
+  //! Whether the energy of a sampled site is a multigroup energy group index
+  //! rather than an energy in [eV]
+  virtual bool energy_is_group_index() const { return false; }
+
 protected:
   // Strategy used for rejecting sites when constraints are applied. KILL means
   // that sites are always accepted but if they don't satisfy constraints, they
@@ -183,11 +187,14 @@ class FileSource : public Source {
   void load_sites_from_file(
     const std::string& path); //!< Load source sites from file
 
+  bool energy_is_group_index() const override { return group_index_; }
+
 protected:
   SourceSite sample(uint64_t* seed) const override;
 
 private:
   vector<SourceSite> sites_; //!< Source sites
+  bool group_index_ {false}; //!< Whether site energies are group indices
 };
 
 //==============================================================================

include/openmc/state_point.h

@@ -36,11 +36,16 @@ void load_state_point();
 //
 // The source_bank variable is used as work space if MPI is used,
 // so it cannot be given as a const span.
+//
+// If group_indices is true, the energies of the sites are multigroup energy
+// group indices, which the file records with an "energy_mode" attribute equal
+// to "multi-group" so that a FileSource reads them as group indices. MCPL
+// files do not record it.
 void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
-  const vector<int64_t>& bank_index);
+  const vector<int64_t>& bank_index, bool group_indices = false);
 
 void write_source_point(std::string, span<SourceSite> source_bank,
-  const vector<int64_t>& bank_index, bool use_mcpl);
+  const vector<int64_t>& bank_index, bool use_mcpl, bool group_indices = false);
 
 // This appends a source bank specification to an HDF5 file
 // that's already open. It is used internally by write_source_point.

openmc/source.py

@@ -816,11 +816,13 @@ class FileSource(SourceBase):
         :class:`openmc.Universe` for which sampled sites must be within. For
         'time_bounds' and 'energy_bounds', the corresponding value is a sequence
         of floats giving the lower and upper bounds on time in [s] or energy in
-        [eV] that the sampled particle must be within. For 'fissionable', the
-        value is a bool indicating that only sites in fissionable material
-        should be accepted. The 'rejection_strategy' indicates what should
-        happen when a source particle is rejected: either 'resample' (pick a new
-        particle) or 'kill' (accept and terminate).
+        [eV] that the sampled particle must be within. Energy bounds are not
+        applied to the sites of a file that holds multigroup energy group
+        indices (see below). For 'fissionable', the value is a bool indicating
+        that only sites in fissionable material should be accepted. The
+        'rejection_strategy' indicates what should happen when a source
+        particle is rejected: either 'resample' (pick a new particle) or 'kill'
+        (accept and terminate).
 
     Attributes
     ----------
@@ -835,6 +837,15 @@ class FileSource(SourceBase):
         'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds',
         'fissionable', and 'rejection_strategy'.
 
+    Notes
+    -----
+    Source point files and state points written by a multigroup simulation
+    store the energy group index of each site rather than its energy (see
+    :ref:`io_source`). A multigroup simulation starts each particle from such a
+    file in its stored group, and a continuous-energy simulation cannot use
+    such a file. Other files, such as surface source files and files written by
+    :func:`openmc.write_source_file`, are read as energies in [eV].
+
     """
 
     def __init__(
@@ -1337,7 +1348,8 @@ class SourceParticle:
     u : iterable of float
         Directional cosines
     E : float
-        Energy of particle in [eV]
+        Energy of particle in [eV]. For a particle read from a source file
+        written by a multigroup simulation, the energy group index.
     time : float
         Time of particle in [s]
     wgt : float
@@ -1435,6 +1447,10 @@ class ParticleList(list):
     def from_hdf5(cls, filename: PathLike) -> ParticleList:
         """Create particle list from an HDF5 file.
 
+        Energies are returned as stored. For a source file written by a
+        multigroup simulation (``energy_mode`` attribute 'multi-group'), they
+        are energy group indices.
+
         Parameters
         ----------
         filename : path-like
@@ -1548,7 +1564,10 @@ def export_to_hdf5(self, filename: PathLike, **kwargs):
         """Export particle list to an HDF5 file.
 
         This method write out an .h5 file that can be used as a source file in
-        conjunction with the :class:`openmc.FileSource` class.
+        conjunction with the :class:`openmc.FileSource` class. The file has no
+        ``energy_mode`` attribute, so OpenMC reads its energies in [eV], also
+        when the particles were read from a source file written by a multigroup
+        simulation.
 
         Parameters
         ----------

src/simulation.cpp

@@ -573,15 +573,15 @@ void finalize_batch()
         settings::path_output, simulation::current_batch, w);
       span<SourceSite> bankspan(simulation::source_bank);
       write_source_point(source_point_filename, bankspan,
-        simulation::work_index, settings::source_mcpl_write);
+        simulation::work_index, settings::source_mcpl_write, !settings::run_CE);
     }
 
     // Write a continously-overwritten source point if requested.
     if (settings::source_latest) {
       auto filename = settings::path_output + "source";
       span<SourceSite> bankspan(simulation::source_bank);
       write_source_point(filename, bankspan, simulation::work_index,
-        settings::source_mcpl_write);
+        settings::source_mcpl_write, !settings::run_CE);
     }
   }
 

src/source.cpp

@@ -214,10 +214,12 @@ SourceSite Source::sample_with_constraints(uint64_t* seed) const
     if (constraints_applied()) {
       accepted = true;
     } else {
-      // Check whether sampled site satisfies constraints
-      accepted = satisfies_spatial_constraints(site.r) &&
-                 satisfies_energy_constraints(site.E) &&
-                 satisfies_time_constraints(site.time);
+      // Check whether sampled site satisfies constraints. Energy bounds are
+      // not applied to multigroup energy group indices.
+      accepted =
+        satisfies_spatial_constraints(site.r) &&
+        (energy_is_group_index() || satisfies_energy_constraints(site.E)) &&
+        satisfies_time_constraints(site.time);
       if (!accepted) {
         ++n_local_reject;
 
@@ -524,6 +526,44 @@ void FileSource::load_sites_from_file(const std::string& path)
       fatal_error("Specified starting source file not a source file type.");
     }
 
+    // Source files written by a multigroup simulation store an energy group
+    // index rather than an energy for each site, which an "energy_mode"
+    // attribute equal to "multi-group" indicates. A state point records the
+    // energy mode of its simulation as a dataset. Without either, the sites
+    // hold energies in [eV].
+    std::string energy_mode {"continuous-energy"};
+    if (filetype == "source") {
+      if (attribute_exists(file_id, "energy_mode")) {
+        read_attribute(file_id, "energy_mode", energy_mode);
+      }
+    } else if (object_exists(file_id, "energy_mode")) {
+      read_dataset(file_id, "energy_mode", energy_mode);
+    }
+    if (energy_mode == "multi-group") {
+      group_index_ = true;
+    } else if (energy_mode != "continuous-energy") {
+      fatal_error(
+        fmt::format("Source file '{}' has an unrecognized energy mode '{}'.",
+          path, energy_mode));
+    }
+
+    if (group_index_) {
+      if (settings::run_CE) {
+        fatal_error(fmt::format(
+          "Source file '{}' holds multigroup energy group indices and cannot "
+          "be used in a continuous-energy simulation.",
+          path));
+      }
+      if (mpi::master &&
+          (energy_bounds_.first != 0.0 ||
+            energy_bounds_.second != std::numeric_limits<double>::max())) {
+        warning(fmt::format(
+          "Energy bounds are not applied to the sites of source file '{}', "
+          "which holds multigroup energy group indices.",
+          path));
+      }
+    }
+
     // Read in the source particles
     read_source_bank(file_id, sites_, false);
 
@@ -1222,6 +1262,18 @@ SourceSite sample_external_source(uint64_t* seed)
       model::external_sources[i]->strength() * n_sources / total_strength;
   }
 
+  // Sites from a source file written by a multigroup simulation already hold
+  // an energy group index, which only needs to be checked
+  if (model::external_sources[i]->energy_is_group_index()) {
+    int n_groups = data::mg.num_energy_groups_;
+    if (!(site.E >= 0.0 && site.E < n_groups && site.E == std::floor(site.E))) {
+      fatal_error(fmt::format("Source site has energy group index {}, which is "
+                              "not valid for {} energy groups.",
+        site.E, n_groups));
+    }
+    return site;
+  }
+
   // If running in MG, convert site.E to group
   if (!settings::run_CE) {
     site.E = lower_bound_index(data::mg.rev_energy_bins_.begin(),

src/state_point.cpp

@@ -596,7 +596,7 @@ hid_t h5banktype(bool memory)
 }
 
 void write_source_point(std::string filename, span<SourceSite> source_bank,
-  const vector<int64_t>& bank_index, bool use_mcpl)
+  const vector<int64_t>& bank_index, bool use_mcpl, bool group_indices)
 {
   std::string ext = use_mcpl ? "mcpl" : "h5";
 
@@ -616,12 +616,13 @@ void write_source_point(std::string filename, span<SourceSite> source_bank,
     write_mcpl_source_point(filename.c_str(), source_bank, bank_index);
   } else {
     filename.append(".h5");
-    write_h5_source_point(filename.c_str(), source_bank, bank_index);
+    write_h5_source_point(
+      filename.c_str(), source_bank, bank_index, group_indices);
   }
 }
 
 void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
-  const vector<int64_t>& bank_index)
+  const vector<int64_t>& bank_index, bool group_indices)
 {
   // 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
@@ -648,6 +649,8 @@ void write_h5_source_point(const char* filename, span<SourceSite> source_bank,
     file_id = file_open(filename_.c_str(), 'w', true);
     write_attribute(file_id, "filetype", "source");
     write_attribute(file_id, "version", VERSION_STATEPOINT);
+    if (group_indices)
+      write_attribute(file_id, "energy_mode", "multi-group");
   }
 
   // Get pointer to source bank and write to file