Merge branch 'develop' into ccfe-spp-2

examples/jupyter/pincell.ipynb

@@ -1373,7 +1373,7 @@
    "source": [
     "## Geometry plotting\n",
     "\n",
-    "We saw before that we could call the `Universe.plot()` method to show a universe while we were creating our geometry. There is also a built-in plotter in the Fortran codebase that is much faster than the Python plotter and has more options. The interface looks somewhat similar to the `Universe.plot()` method. Instead though, we create `Plot` instances, assign them to a `Plots` collection, export it to XML, and then run OpenMC in geometry plotting mode. As an example, let's specify that we want the plot to be colored by material (rather than by cell) and we assign yellow to fuel and blue to water."
+    "We saw before that we could call the `Universe.plot()` method to show a universe while we were creating our geometry. There is also a built-in plotter in the codebase that is much faster than the Python plotter and has more options. The interface looks somewhat similar to the `Universe.plot()` method. Instead though, we create `Plot` instances, assign them to a `Plots` collection, export it to XML, and then run OpenMC in geometry plotting mode. As an example, let's specify that we want the plot to be colored by material (rather than by cell) and we assign yellow to fuel and blue to water."
    ]
   },
   {

examples/jupyter/post-processing.ipynb

@@ -911,7 +911,7 @@
     "print(sum(probability*np.diff(energy_bins)))\n",
     "\n",
     "# Plot source energy PDF\n",
-    "plt.semilogx(energy_bins[:-1], probability*np.diff(energy_bins), linestyle='steps')\n",
+    "plt.semilogx(energy_bins[:-1], probability*np.diff(energy_bins), drawstyle='steps')\n",
     "plt.xlabel('Energy (eV)')\n",
     "plt.ylabel('Probability/eV')"
    ]

include/openmc/error.h

@@ -5,7 +5,10 @@
 #include <string>
 #include <sstream>
 
+#include <fmt/format.h>
+
 #include "openmc/capi.h"
+#include "openmc/settings.h"
 
 #ifdef __GNUC__
 #define UNREACHABLE() __builtin_unreachable()
@@ -63,6 +66,21 @@ void write_message(const std::stringstream& message, int level)
   write_message(message.str(), level);
 }
 
+template<typename... Params>
+void write_message(
+  int level, const std::string& message, const Params&... fmt_args)
+{
+  if (settings::verbosity >= level) {
+    write_message(fmt::format(message, fmt_args...));
+  }
+}
+
+template<typename... Params>
+void write_message(const std::string& message, const Params&... fmt_args)
+{
+  write_message(fmt::format(message, fmt_args...));
+}
+
 #ifdef OPENMC_MPI
 extern "C" void abort_mpi(int code);
 #endif

openmc/deplete/abc.py

@@ -879,16 +879,29 @@ class Integrator(ABC):
         return (self.operator.prev_res[-1].time[-1],
                 len(self.operator.prev_res) - 1)
 
-    def integrate(self):
-        """Perform the entire depletion process across all steps"""
+    def integrate(self, final_step=True):
+        """Perform the entire depletion process across all steps
+
+        Parameters
+        ----------
+        final_step : bool, optional
+            Indicate whether or not a transport solve should be run at the end
+            of the last timestep.
+
+            .. versionadded:: 0.12.1
+
+        """
         with self.operator as conc:
             t, self._i_res = self._get_start_data()
 
             for i, (dt, source_rate) in enumerate(self):
+                # Solve transport equation (or obtain result from restart)
                 if i > 0 or self.operator.prev_res is None:
                     conc, res = self._get_bos_data_from_operator(i, source_rate, conc)
                 else:
                     conc, res = self._get_bos_data_from_restart(i, source_rate, conc)
+
+                # Solve Bateman equations over time interval
                 proc_time, conc_list, res_list = self(conc, res.rates, dt, source_rate, i)
 
                 # Insert BOS concentration, transport results
@@ -903,8 +916,11 @@ class Integrator(ABC):
 
                 t += dt
 
-            # Final simulation
-            res_list = [self.operator(conc, source_rate)]
+            # Final simulation -- in the case that final_step is False, a zero
+            # source rate is passed to the transport operator (which knows to
+            # just return zero reaction rates without actually doing a transport
+            # solve)
+            res_list = [self.operator(conc, source_rate if final_step else 0.0)]
             Results.save(self.operator, [conc], res_list, [t, t],
                          source_rate, self._i_res + len(self), proc_time)
             self.operator.write_bos_data(len(self) + self._i_res)

openmc/deplete/operator.py

@@ -277,6 +277,16 @@ class Operator(TransportOperator):
             Eigenvalue and reaction rates resulting from transport operator
 
         """
+        # Reset results in OpenMC
+        openmc.lib.reset()
+
+        # If the source rate is zero, return zero reaction rates without running
+        # a transport solve
+        if source_rate == 0.0:
+            rates = self.reaction_rates.copy()
+            rates.fill(0.0)
+            return OperatorResult(ufloat(0.0, 0.0), rates)
+
         # Prevent OpenMC from complaining about re-creating tallies
         openmc.reset_auto_ids()
 
@@ -291,7 +301,6 @@ class Operator(TransportOperator):
         self._yield_helper.update_tally_nuclides(nuclides)
 
         # Run OpenMC
-        openmc.lib.reset()
         openmc.lib.run()
         openmc.lib.reset_timers()
 

src/cross_sections.cpp

@@ -228,7 +228,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
         int idx = data::library_map[key];
         std::string& filename = data::libraries[idx].path_;
 
-        write_message("Reading " + name + " from " + filename, 6);
+        write_message(6, "Reading {} from {}", name, filename);
 
         // Open file and make sure version matches
         hid_t file_id = file_open(filename, 'r');
@@ -256,10 +256,8 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
   }
 
   // Show minimum/maximum temperature
-  write_message("Minimum neutron data temperature: " +
-    std::to_string(data::temperature_min) + " K", 4);
-  write_message("Maximum neutron data temperature: " +
-    std::to_string(data::temperature_max) + " K", 4);
+  write_message(4, "Minimum neutron data temperature: {} K", data::temperature_min);
+  write_message(4, "Maximum neutron data temperature: {} K", data::temperature_max);
 
   // If the user wants multipole, make sure we found a multipole library.
   if (settings::temperature_multipole) {

src/mesh.cpp

@@ -2051,7 +2051,7 @@ void
 UnstructuredMesh::write(std::string base_filename) const {
   // add extension to the base name
   auto filename = base_filename + ".vtk";
-  write_message("Writing unstructured mesh " + filename + "...", 5);
+  write_message(5, "Writing unstructured mesh {}...", filename);
   filename = settings::path_output + filename;
 
   // write the tetrahedral elements of the mesh only

src/mgxs_interface.cpp

@@ -3,6 +3,8 @@
 #include <string>
 #include <unordered_set>
 
+#include <fmt/format.h>
+
 #include "openmc/cell.h"
 #include "openmc/cross_sections.h"
 #include "openmc/container_util.h"
@@ -57,8 +59,7 @@ void MgxsInterface::init()
   // Check if MGXS Library exists
   if (!file_exists(cross_sections_path_)) {
     // Could not find MGXS Library file
-    fatal_error("Cross sections HDF5 file '" + cross_sections_path_ +
-      "' does not exist.");
+    fatal_error(fmt::format("Cross sections HDF5 file '{}' does not exist!", cross_sections_path_));
   }
 
   write_message("Loading cross section data...", 5);
@@ -98,15 +99,14 @@ void
 MgxsInterface::add_mgxs(hid_t file_id, const std::string& name,
   const std::vector<double>& temperature)
 {
-  write_message("Loading " + std::string(name) + " data...", 6);
+  write_message(5, "Loading {} data...", name);
 
   // Check to make sure cross section set exists in the library
   hid_t xs_grp;
   if (object_exists(file_id, name.c_str())) {
     xs_grp = open_group(file_id, name.c_str());
   } else {
-    fatal_error("Data for " + std::string(name) + " does not exist in "
-                + "provided MGXS Library");
+    fatal_error(fmt::format("Data for {} does not exist in provided MGXS Library", name));
   }
 
   nuclides_.emplace_back(xs_grp, temperature, num_energy_groups_,
@@ -187,8 +187,7 @@ void MgxsInterface::read_header(const std::string& path_cross_sections)
   // Check if MGXS Library exists
   if (!file_exists(cross_sections_path_)) {
     // Could not find MGXS Library file
-    fatal_error("Cross sections HDF5 file '" + cross_sections_path_ +
-      "' does not exist.");
+    fatal_error(fmt::format("Cross section HDF5 file '{}' does not exist", cross_sections_path_));
   }
   write_message("Reading cross sections HDF5 file...", 5);
 

src/nuclide.cpp

@@ -944,7 +944,7 @@ extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n)
     // Get filename for library containing nuclide
     int idx = it->second;
     const auto& filename = data::libraries[idx].path_;
-    write_message("Reading " + std::string{name} + " from " + filename, 6);
+    write_message(6, "Reading {} from {}", name, filename);
 
     // Open file and make sure version is sufficient
     hid_t file_id = file_open(filename, 'r');
@@ -977,7 +977,7 @@ extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n)
 
         int idx = it->second;
         const auto& filename = data::libraries[idx].path_;
-        write_message("Reading " + element + " from " + filename, 6);
+        write_message(6, "Reading {} from {} ", element, filename);
 
         // Open file and make sure version is sufficient
         hid_t file_id = file_open(filename, 'r');

src/particle.cpp

@@ -410,7 +410,7 @@ Particle::cross_surface()
   // TODO: off-by-one
   const auto& surf {model::surfaces[i_surface - 1].get()};
   if (settings::verbosity >= 10 || trace_) {
-    write_message("    Crossing surface " + std::to_string(surf->id_));
+    write_message("    Crossing surface {}", surf->id_);
   }
 
   if (surf->bc_ == Surface::BoundaryType::VACUUM && (settings::run_mode != RunMode::PLOTTING)) {
@@ -437,7 +437,7 @@ Particle::cross_surface()
 
     // Display message
     if (settings::verbosity >= 10 || trace_) {
-      write_message("    Leaked out of surface " + std::to_string(surf->id_));
+      write_message("    Leaked out of surface {}", surf->id_);
     }
     return;
 
@@ -502,7 +502,7 @@ Particle::cross_surface()
 
     // Diagnostic message
     if (settings::verbosity >= 10 || trace_) {
-      write_message("    Reflected from surface " + std::to_string(surf->id_));
+      write_message("    Reflected from surface {}", surf->id_);
     }
     return;
 
@@ -556,8 +556,7 @@ Particle::cross_surface()
 
     // Diagnostic message
     if (settings::verbosity >= 10 || trace_) {
-      write_message("    Hit periodic boundary on surface " +
-        std::to_string(surf->id_));
+      write_message("    Hit periodic boundary on surface {}", surf->id_);
     }
     return;
   }

src/particle_restart.cpp

@@ -24,8 +24,7 @@ namespace openmc {
 void read_particle_restart(Particle& p, RunMode& previous_run_mode)
 {
   // Write meessage
-  write_message("Loading particle restart file " +
-    settings::path_particle_restart + "...", 5);
+  write_message(5, "Loading particle restart file {}", settings::path_particle_restart);
 
   // Open file
   hid_t file_id = file_open(settings::path_particle_restart, 'r');

src/plot.cpp

@@ -94,8 +94,7 @@ extern "C"
 int openmc_plot_geometry()
 {
   for (auto pl : model::plots) {
-    write_message(fmt::format("Processing plot {}: {}...",
-      pl.id_, pl.path_plot_), 5);
+    write_message(5, "Processing plot {}: {}...", pl.id_, pl.path_plot_);
 
     if (PlotType::slice == pl.type_) {
       // create 2D image
@@ -114,7 +113,7 @@ void read_plots_xml()
   // Check if plots.xml exists
   std::string filename = settings::path_input + "plots.xml";
   if (!file_exists(filename)) {
-    fatal_error("Plots XML file '" + filename + "' does not exist!");
+    fatal_error(fmt::format("Plots XML file '{}' does not exist!", filename));
   }
 
   write_message("Reading plot XML file...", 5);

src/reaction_product.cpp

@@ -3,7 +3,10 @@
 #include <memory> // for unique_ptr
 #include <string> // for string
 
+#include <fmt/core.h>
+
 #include "openmc/endf.h"
+#include "openmc/error.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/particle.h"
 #include "openmc/random_lcg.h"
@@ -36,8 +39,14 @@ ReactionProduct::ReactionProduct(hid_t group)
   }
 
   // Read decay rate for delayed emission
-  if (emission_mode_ == EmissionMode::delayed)
-    read_attribute(group, "decay_rate", decay_rate_);
+  if (emission_mode_ == EmissionMode::delayed) {
+    if (attribute_exists(group, "decay_rate")) {
+      read_attribute(group, "decay_rate", decay_rate_);
+    } else if (particle_ == Particle::Type::neutron) {
+      warning(fmt::format("Decay rate doesn't exist for delayed neutron "
+        "emission ({}).", object_name(group)));
+    }
+  }
 
   // Read secondary particle yield
   yield_ = read_function(group, "yield");

src/simulation.cpp

@@ -33,6 +33,8 @@
 #include <mpi.h>
 #endif
 
+#include <fmt/format.h>
+
 #include <algorithm>
 #include <cmath>
 #include <string>
@@ -305,8 +307,7 @@ void initialize_batch()
   ++simulation::current_batch;
 
   if (settings::run_mode == RunMode::FIXED_SOURCE) {
-    int b = simulation::current_batch;
-    write_message("Simulating batch " + std::to_string(b), 6);
+    write_message(6, "Simulating batch {}", simulation::current_batch);
   }
 
   // Reset total starting particle weight used for normalizing tallies
@@ -500,7 +501,7 @@ void initialize_history(Particle& p, int64_t index_source)
 
   // Display message if high verbosity or trace is on
   if (settings::verbosity >= 9 || p.trace_) {
-    write_message("Simulating Particle " + std::to_string(p.id_));
+    write_message("Simulating Particle {}", p.id_);
   }
 
   // Add paricle's starting weight to count for normalizing tallies later
@@ -610,9 +611,8 @@ void initialize_data()
       double max_E = nuc->grid_[0].energy.back();
       int neutron = static_cast<int>(Particle::Type::neutron);
       if (max_E == data::energy_max[neutron]) {
-        write_message("Maximum neutron transport energy: " +
-          std::to_string(data::energy_max[neutron]) + " eV for " +
-          nuc->name_, 7);
+        write_message(7, "Maximum neutron transport energy: {} eV for {}",
+          data::energy_max[neutron], nuc->name_);
         if (mpi::master && data::energy_max[neutron] < 20.0e6) {
           warning("Maximum neutron energy is below 20 MeV. This may bias "
             "the results.");

src/source.cpp

@@ -263,8 +263,7 @@ void initialize_source()
     // Read the source from a binary file instead of sampling from some
     // assumed source distribution
 
-    write_message(fmt::format("Reading source file from {}...",
-      settings::path_source), 6);
+    write_message(6, "Reading source file from {}...", settings::path_source);
 
     // Open the binary file
     hid_t file_id = file_open(settings::path_source, 'r', true);
@@ -285,13 +284,13 @@ void initialize_source()
     file_close(file_id);
   } else if (!settings::path_source_library.empty()) {
 
-    write_message(fmt::format("Sampling from library source {}...",
-      settings::path_source), 6);
+    write_message(6, "Sampling library source {}...", settings::path_source);
 
     fill_source_bank_custom_source();
 
   } else {
     // Generation source sites from specified distribution in user input
+    #pragma omp parallel for
     for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
       // initialize random number seed
       int64_t id = simulation::total_gen*settings::n_particles +

src/tallies/tally.cpp

@@ -239,10 +239,12 @@ score_str_to_int(std::string score_str)
   try {
     MT = std::stoi(score_str);
   } catch (const std::invalid_argument& ex) {
-    throw std::invalid_argument("Invalid tally score \"" + score_str + "\"");
+    throw std::invalid_argument("Invalid tally score \"" + score_str + "\". See the docs "
+      "for details: https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");
   }
   if (MT < 1)
-    throw std::invalid_argument("Invalid tally score \"" + score_str + "\"");
+    throw std::invalid_argument("Invalid tally score \"" + score_str + "\". See the docs "
+      "for details: https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");
   return MT;
 }
 

src/tallies/trigger.cpp

@@ -164,8 +164,7 @@ check_triggers()
   // If all the triggers are satisfied, alert the user and return.
   if (std::max(keff_ratio, tally_ratio) <= 1.) {
     simulation::satisfy_triggers = true;
-    write_message("Triggers satisfied for batch "
-      + std::to_string(current_batch), 7);
+    write_message(7, "Triggers satisfied for batch {}", current_batch);
     return;
   }
 

src/volume_calc.cpp

@@ -468,9 +468,7 @@ int openmc_calculate_volumes() {
   time_volume.start();
 
   for (int i = 0; i < model::volume_calcs.size(); ++i) {
-    if (mpi::master) {
-      write_message(fmt::format("Running volume calculation {}...", i + 1), 4);
-    }
+    write_message(4, "Running volume calculation {}", i+1);
 
     // Run volume calculation
     const auto& vol_calc {model::volume_calcs[i]};
@@ -488,8 +486,8 @@ int openmc_calculate_volumes() {
 
       // Display domain volumes
       for (int j = 0; j < vol_calc.domain_ids_.size(); j++) {
-        write_message(fmt::format("{}{}: {} +/- {} cm^3", domain_type,
-          vol_calc.domain_ids_[j], results[j].volume[0], results[j].volume[1]), 4);
+        write_message(4, "{}{}: {} +/- {} cm^3", domain_type,
+          vol_calc.domain_ids_[j], results[j].volume[0], results[j].volume[1]);
       }
 
       // Write volumes to HDF5 file
@@ -502,9 +500,7 @@ int openmc_calculate_volumes() {
 
   // Show elapsed time
   time_volume.stop();
-  if (mpi::master) {
-    write_message(fmt::format("Elapsed time: {} s", time_volume.elapsed()), 6);
-  }
+  write_message(6, "Elapsed time: {} s", time_volume.elapsed());
 
   return 0;
 }

src/wmp.cpp

@@ -5,6 +5,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"
 #include "openmc/nuclide.h"
+#include "openmc/error.h"  // for writing messages
 
 #include <fmt/core.h>
 
@@ -227,7 +228,7 @@ void read_multipole_data(int i_nuclide)
   std::string& filename = data::libraries[idx].path_;
 
   // Display message
-  write_message("Reading " + nuc->name_ + " WMP data from " + filename, 6);
+  write_message(6, "Reading {} WMP data from {}", nuc->name_, filename);
 
   // Open file and make sure version is sufficient
   hid_t file = file_open(filename, 'r');