Merge pull request #1128 from paulromano/cpp-physics

Move source/fission banks and neutron physics to C++

CMakeLists.txt

@@ -381,6 +381,7 @@ add_library(libopenmc SHARED
   src/tallies/tally_header.F90
   src/tallies/trigger.F90
   src/tallies/trigger_header.F90
+  src/bank.cpp
   src/dagmc.cpp
   src/cell.cpp
   src/cmfd_execute.cpp
@@ -408,6 +409,8 @@ add_library(libopenmc SHARED
   src/nuclide.cpp
   src/output.cpp
   src/particle.cpp
+  src/photon.cpp
+  src/physics.cpp
   src/physics_common.cpp
   src/physics_mg.cpp
   src/plot.cpp

include/openmc/bank.h

@@ -0,0 +1,47 @@
+#ifndef OPENMC_BANK_H
+#define OPENMC_BANK_H
+
+#include <cstdint>
+#include <vector>
+
+namespace openmc {
+
+struct Bank {
+  double wgt;
+  double xyz[3];
+  double uvw[3];
+  double E;
+  int delayed_group;
+  int particle;
+};
+
+} // namespace openmc
+
+// Without an explicit instantiation of vector<Bank>, the Intel compiler
+// will complain about the threadprivate directive on filter_matches. Note that
+// this has to happen *outside* of the openmc namespace
+extern template class std::vector<openmc::Bank>;
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace simulation {
+
+extern "C" int64_t n_bank;
+
+extern std::vector<Bank> source_bank;
+extern std::vector<Bank> fission_bank;
+#ifdef _OPENMP
+extern std::vector<Bank> master_fission_bank;
+#endif
+
+#pragma omp threadprivate(fission_bank, n_bank)
+
+} // namespace simulation
+
+} // namespace openmc
+
+#endif // OPENMC_BANK_H

include/openmc/capi.h

@@ -6,9 +6,11 @@
 #include <stddef.h>
 
 #ifdef __cplusplus
+#include "openmc/bank.h"
 extern "C" {
-#endif
-
+  int openmc_fission_bank(openmc::Bank** ptr, int64_t* n);
+  int openmc_source_bank(openmc::Bank** ptr, int64_t* n);
+#else
   struct Bank {
     double wgt;
     double xyz[3];
@@ -18,6 +20,10 @@ extern "C" {
     int particle;
   };
 
+  int openmc_fission_bank(struct Bank** ptr, int64_t* n);
+  int openmc_source_bank(struct Bank** ptr, int64_t* n);
+#endif
+
   int openmc_calculate_volumes();
   int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n);
   int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
@@ -38,7 +44,6 @@ extern "C" {
   int openmc_filter_set_type(int32_t index, const char* type);
   int openmc_finalize();
   int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance);
-  int openmc_fission_bank(struct Bank** ptr, int64_t* n);
   int openmc_get_cell_index(int32_t id, int32_t* index);
   int openmc_get_filter_index(int32_t id, int32_t* index);
   void openmc_get_filter_next_id(int32_t* id);
@@ -85,7 +90,6 @@ extern "C" {
   void openmc_set_seed(int64_t new_seed);
   int openmc_simulation_finalize();
   int openmc_simulation_init();
-  int openmc_source_bank(struct Bank** ptr, int64_t* n);
   int openmc_spatial_legendre_filter_get_order(int32_t index, int* order);
   int openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, double* min, double* max);
   int openmc_spatial_legendre_filter_set_order(int32_t index, int order);

include/openmc/eigenvalue.h

@@ -25,9 +25,6 @@ extern std::array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
 extern std::vector<double> entropy; //!< Shannon entropy at each generation
 extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
 
-extern "C" int64_t n_bank;
-#pragma omp threadprivate(n_bank)
-
 } // namespace simulation
 
 //==============================================================================
@@ -35,15 +32,22 @@ extern "C" int64_t n_bank;
 //==============================================================================
 
 //! Collect/normalize the tracklength keff from each process
-extern "C" void calculate_generation_keff();
+void calculate_generation_keff();
 
 //! Calculate mean/standard deviation of keff during active generations
 //!
 //! This function sets the global variables keff and keff_std which represent
 //! the mean and standard deviation of the mean of k-effective over active
 //! generations. It also broadcasts the value from the master process.
-extern "C" void calculate_average_keff();
+void calculate_average_keff();
 
+#ifdef _OPENMP
+//! Join threadprivate fission banks into a single fission bank
+//!
+//! Note that this operation is necessarily sequential to preserve the order of
+//! the bank when using varying numbers of threads.
+void join_bank_from_threads();
+#endif
 
 //! Calculates a minimum variance estimate of k-effective
 //!
@@ -60,16 +64,16 @@ extern "C" void calculate_average_keff();
 extern "C" int openmc_get_keff(double* k_combined);
 
 //! Sample/redistribute source sites from accumulated fission sites
-extern "C" void synchronize_bank();
+void synchronize_bank();
 
 //! Calculates the Shannon entropy of the fission source distribution to assess
 //! source convergence
-extern "C" void shannon_entropy();
+void shannon_entropy();
 
 //! Determines the source fraction in each UFS mesh cell and reweights the
 //! source bank so that the sum of the weights is equal to n_particles. The
 //! 'source_frac' variable is used later to bias the production of fission sites
-extern "C" void ufs_count_sites();
+void ufs_count_sites();
 
 //! Get UFS weight corresponding to particle's location
 extern "C" double ufs_get_weight(const Particle* p);

include/openmc/endf.h

@@ -22,6 +22,17 @@ Interpolation int2interp(int i);
 //! \return Whether corresponding reaction is a fission reaction
 bool is_fission(int MT);
 
+//! Determine if a given MT number is that of a disappearance reaction, i.e., a
+//! reaction with no neutron in the exit channel
+//! \param[in] MT ENDF MT value
+//! \return Whether corresponding reaction is a disappearance reaction
+bool is_disappearance(int MT);
+
+//! Determine if a given MT number is that of an inelastic scattering reaction
+//! \param[in] MT ENDF MT value
+//! \return Whether corresponding reaction is an inelastic scattering reaction
+bool is_inelastic_scatter(int MT);
+
 //==============================================================================
 //! Abstract one-dimensional function
 //==============================================================================

include/openmc/material.h

@@ -43,5 +43,11 @@ public:
   explicit Material(pugi::xml_node material_node);
 };
 
+//==============================================================================
+// Fortran compatibility
+//==============================================================================
+
+extern "C" bool material_isotropic(int i_material, int i_nuc_mat);
+
 } // namespace openmc
 #endif // OPENMC_MATERIAL_H

include/openmc/nuclide.h

@@ -5,23 +5,75 @@
 #define OPENMC_NUCLIDE_H
 
 #include <array>
+#include <memory> // for unique_ptr
+#include <vector>
+
+#include <hdf5.h>
 
 #include "openmc/constants.h"
+#include "openmc/endf.h"
+#include "openmc/reaction.h"
+#include "openmc/reaction_product.h"
 
 namespace openmc {
 
 //==============================================================================
-// Global variables
+// Constants
 //==============================================================================
 
-namespace data {
+constexpr double CACHE_INVALID {-1.0};
 
-// Minimum/maximum transport energy for each particle type. Order corresponds to
-// that of the ParticleType enum
-extern std::array<double, 2> energy_min;
-extern std::array<double, 2> energy_max;
+//===============================================================================
+// Data for a nuclide
+//===============================================================================
 
-} // namespace data
+class Nuclide {
+public:
+  // Types, aliases
+  using EmissionMode = ReactionProduct::EmissionMode;
+  struct EnergyGrid {
+    std::vector<int> grid_index;
+    std::vector<double> energy;
+  };
+
+  // Constructors
+  Nuclide(hid_t group, const double* temperature, int n);
+
+  // Methods
+  double nu(double E, EmissionMode mode, int group=0) const;
+  void calculate_elastic_xs(int i_nuclide) const;
+
+  //! Determines the microscopic 0K elastic cross section at a trial relative
+  //! energy used in resonance scattering
+  double elastic_xs_0K(double E) const;
+
+  // Data members
+  std::string name_; //! Name of nuclide, e.g. "U235"
+  int Z_; //! Atomic number
+  int A_; //! Mass number
+  int metastable_; //! Metastable state
+  double awr_; //! Atomic weight ratio
+  std::vector<double> kTs_; //! temperatures in eV (k*T)
+  std::vector<EnergyGrid> grid_; //! Energy grid at each temperature
+
+  bool fissionable_ {false}; //! Whether nuclide is fissionable
+  bool has_partial_fission_ {false}; //! has partial fission reactions?
+  std::vector<Reaction*> fission_rx_; //! Fission reactions
+  int n_precursor_ {0}; //! Number of delayed neutron precursors
+  std::unique_ptr<Function1D> total_nu_; //! Total neutron yield
+
+  // Resonance scattering information
+  bool resonant_ {false};
+  std::vector<double> energy_0K_;
+  std::vector<double> elastic_0K_;
+  std::vector<double> xs_cdf_;
+
+  std::vector<std::unique_ptr<Reaction>> reactions_; //! Reactions
+  std::vector<int> index_inelastic_scatter_;
+
+private:
+  void create_derived();
+};
 
 //===============================================================================
 //! Cached microscopic cross sections for a particular nuclide at the current
@@ -66,19 +118,53 @@ struct NuclideMicroXS {
 // particle is traveling through
 //===============================================================================
 
-  struct MaterialMacroXS {
-    double total;         //!< macroscopic total xs
-    double absorption;    //!< macroscopic absorption xs
-    double fission;       //!< macroscopic fission xs
-    double nu_fission;    //!< macroscopic production xs
-    double photon_prod;   //!< macroscopic photon production xs
+struct MaterialMacroXS {
+  double total;         //!< macroscopic total xs
+  double absorption;    //!< macroscopic absorption xs
+  double fission;       //!< macroscopic fission xs
+  double nu_fission;    //!< macroscopic production xs
+  double photon_prod;   //!< macroscopic photon production xs
+
+  // Photon cross sections
+  double coherent;        //!< macroscopic coherent xs
+  double incoherent;      //!< macroscopic incoherent xs
+  double photoelectric;   //!< macroscopic photoelectric xs
+  double pair_production; //!< macroscopic pair production xs
+};
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+
+// Minimum/maximum transport energy for each particle type. Order corresponds to
+// that of the ParticleType enum
+extern std::array<double, 2> energy_min;
+extern std::array<double, 2> energy_max;
+
+extern std::vector<std::unique_ptr<Nuclide>> nuclides;
+
+} // namespace data
+
+namespace simulation {
+
+// Cross section caches
+extern NuclideMicroXS* micro_xs;
+extern "C" MaterialMacroXS material_xs;
+#pragma omp threadprivate(micro_xs, material_xs)
+
+} // namespace simulation
+
+//==============================================================================
+// Fortran compatibility
+//==============================================================================
+
+extern "C" void set_micro_xs();
+extern "C" bool nuclide_wmp_present(int i_nuclide);
+extern "C" double nuclide_wmp_emin(int i_nuclide);
+extern "C" double nuclide_wmp_emax(int i_nuclide);
 
-    // Photon cross sections
-    double coherent;        //!< macroscopic coherent xs
-    double incoherent;      //!< macroscopic incoherent xs
-    double photoelectric;   //!< macroscopic photoelectric xs
-    double pair_production; //!< macroscopic pair production xs
-  };
 
 } // namespace openmc
 

include/openmc/photon.h

@@ -0,0 +1,37 @@
+#ifndef OPENMC_PHOTON_H
+#define OPENMC_PHOTON_H
+
+#include <hdf5.h>
+
+#include <string>
+#include <vector>
+
+namespace openmc {
+
+//==============================================================================
+//! Photon interaction data for a single element
+//==============================================================================
+
+class PhotonInteraction {
+public:
+  // Constructors
+  PhotonInteraction(hid_t group);
+
+  // Data members
+  std::string name_; //! Name of element, e.g. "Zr"
+  int Z_; //! Atomic number
+};
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+
+extern std::vector<PhotonInteraction> elements;
+
+} // namespace data
+
+} // namespace openmc
+
+#endif // OPENMC_PHOTON_H

include/openmc/physics.h

@@ -0,0 +1,89 @@
+#ifndef OPENMC_PHYSICS_H
+#define OPENMC_PHYSICS_H
+
+#include "openmc/bank.h"
+#include "openmc/nuclide.h"
+#include "openmc/particle.h"
+#include "openmc/position.h"
+#include "openmc/reaction.h"
+
+namespace openmc {
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+//! Sample a nuclide and reaction and then calls the appropriate routine
+extern "C" void collision(Particle* p);
+
+//! Samples an incident neutron reaction
+void sample_neutron_reaction(Particle* p);
+
+//! Samples an element based on the macroscopic cross sections for each nuclide
+//! within a material and then samples a reaction for that element and calls the
+//! appropriate routine to process the physics.
+extern "C" void sample_photon_reaction(Particle* p);
+
+//! Terminates the particle and either deposits all energy locally
+//! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung
+//! photons from electron deflections with charged particles (electron_treatment
+//! = ELECTRON_TTB).
+void sample_electron_reaction(Particle* p);
+
+//! Terminates the particle and either deposits all energy locally
+//! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung
+//! photons from electron deflections with charged particles (electron_treatment
+//! = ELECTRON_TTB). Two annihilation photons of energy MASS_ELECTRON_EV (0.511
+//! MeV) are created and travel in opposite directions.
+void sample_positron_reaction(Particle* p);
+
+void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat);
+
+//! Determine the average total, prompt, and delayed neutrons produced from
+//! fission and creates appropriate bank sites.
+void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
+  Bank* bank_array, int64_t* bank_size, int64_t bank_capacity);
+
+// void sample_element(Particle* p);
+
+Reaction* sample_fission(int i_nuclide, double E);
+
+void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
+
+void absorption(Particle* p, int i_nuclide);
+
+void scatter(Particle*, int i_nuclide, int i_nuc_mat);
+
+//! Treats the elastic scattering of a neutron with a target.
+void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
+  double* uvw, double* mu_lab, double* wgt);
+
+extern "C" void sab_scatter(int i_nuclide, int i_sab, double* E,
+  double* uvw, double* mu);
+
+//! samples the target velocity. The constant cross section free gas model is
+//! the default method. Methods for correctly accounting for the energy
+//! dependence of cross sections in treating resonance elastic scattering such
+//! as the DBRC, WCM, and a new, accelerated scheme are also implemented here.
+Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
+   Direction v_neut, double xs_eff, double kT, double* wgt);
+
+//! samples a target velocity based on the free gas scattering formulation, used
+//! by most Monte Carlo codes, in which cross section is assumed to be constant
+//! in energy. Excellent documentation for this method can be found in
+//! FRA-TM-123.
+Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT);
+
+void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site);
+
+//! handles all reactions with a single secondary neutron (other than fission),
+//! i.e. level scattering, (n,np), (n,na), etc.
+void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p);
+
+void sample_secondary_photons(Particle* p, int i_nuclide);
+
+extern "C" void thick_target_bremsstrahlung(Particle* p, double* E_lost);
+
+} // namespace openmc
+
+#endif // OPENMC_PHYSICS_H

include/openmc/position.h

@@ -25,6 +25,8 @@ struct Position {
   Position& operator-=(double);
   Position& operator*=(Position);
   Position& operator*=(double);
+  Position& operator/=(Position);
+  Position& operator/=(double);
 
   const double& operator[](int i) const {
     switch (i) {
@@ -76,6 +78,10 @@ inline Position operator*(Position a, Position b) { return a *= b; }
 inline Position operator*(Position a, double b)   { return a *= b; }
 inline Position operator*(double a, Position b)   { return b *= a; }
 
+inline Position operator/(Position a, Position b) { return a /= b; }
+inline Position operator/(Position a, double b)   { return a /= b; }
+inline Position operator/(double a, Position b)   { return b /= a; }
+
 inline bool operator==(Position a, Position b)
 {return a.x == b.x && a.y == b.y && a.z == b.z;}
 

include/openmc/reaction.h

@@ -4,6 +4,7 @@
 #ifndef OPENMC_REACTION_H
 #define OPENMC_REACTION_H
 
+#include <string>
 #include <vector>
 
 #include "hdf5.h"
@@ -39,13 +40,17 @@ public:
   std::vector<ReactionProduct> products_; //!< Reaction products
 };
 
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+std::string reaction_name(int mt);
+
 //==============================================================================
 // Fortran compatibility functions
 //==============================================================================
 
 extern "C" {
-  Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n);
-  void reaction_delete(Reaction* rx);
   int reaction_mt(Reaction* rx);
   double reaction_q_value(Reaction* rx);
   bool reaction_scatter_in_cm(Reaction* rx);

include/openmc/simulation.h

@@ -56,6 +56,9 @@ extern "C" int thread_id;  //!< ID of a given thread
 // Functions
 //==============================================================================
 
+//! Allocate space for source and fission banks
+void allocate_banks();
+
 //! Determine number of particles to transport per process
 void calculate_work();
 

include/openmc/source.h

@@ -9,9 +9,9 @@
 
 #include "pugixml.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/distribution_multi.h"
 #include "openmc/distribution_spatial.h"
-#include "openmc/capi.h"
 #include "openmc/particle.h"
 
 namespace openmc {
@@ -62,7 +62,7 @@ extern "C" void initialize_source();
 //! Sample a site from all external source distributions in proportion to their
 //! source strength
 //! \return Sampled source site
-extern "C" Bank sample_external_source();
+Bank sample_external_source();
 
 //! Fill source bank at end of generation for fixed source simulations
 void fill_source_bank_fixedsource();

include/openmc/state_point.h

@@ -10,8 +10,8 @@
 namespace openmc {
 
 void write_source_point(const char* filename);
-extern "C" void write_source_bank(hid_t group_id, Bank* source_bank);
-extern "C" void read_source_bank(hid_t group_id, Bank* source_bank);
+extern "C" void write_source_bank(hid_t group_id);
+extern "C" void read_source_bank(hid_t group_id);
 extern "C" void write_tally_results_nr(hid_t file_id);
 extern "C" void restart_set_keff();
 

openmc/capi/__init__.py

@@ -38,8 +38,11 @@ else:
     from unittest.mock import Mock
     _dll = Mock()
 
-dagmc_enabled = bool(c_bool.in_dll(_dll, "dagmc_enabled"))
-    
+
+def _dagmc_enabled():
+    return c_bool.in_dll(_dll, "dagmc_enabled").value
+
+
 from .error import *
 from .core import *
 from .nuclide import *

src/api.F90

@@ -82,6 +82,9 @@ contains
 
       subroutine free_memory_settings() bind(C)
       end subroutine free_memory_settings
+
+      subroutine free_memory_bank() bind(C)
+      end subroutine free_memory_bank
     end interface
 
     call free_memory_geometry()

src/bank.cpp

@@ -0,0 +1,107 @@
+#include "openmc/bank.h"
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+
+#include <cstdint>
+
+// Explicit template instantiation definition
+template class std::vector<openmc::Bank>;
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace simulation {
+
+int64_t n_bank;
+
+std::vector<Bank> source_bank;
+std::vector<Bank> fission_bank;
+#ifdef _OPENMP
+std::vector<Bank> master_fission_bank;
+#endif
+
+} // namespace simulation
+
+//==============================================================================
+// C API
+//==============================================================================
+
+extern "C" int openmc_source_bank(Bank** ptr, int64_t* n)
+{
+  if (simulation::source_bank.size() == 0) {
+    set_errmsg("Source bank has not been allocated.");
+    return OPENMC_E_ALLOCATE;
+  } else {
+    *ptr = simulation::source_bank.data();
+    *n = simulation::source_bank.size();
+    return 0;
+  }
+}
+
+extern "C" int openmc_fission_bank(Bank** ptr, int64_t* n)
+{
+  if (simulation::fission_bank.size() == 0) {
+    set_errmsg("Fission bank has not been allocated.");
+    return OPENMC_E_ALLOCATE;
+  } else {
+    *ptr = simulation::fission_bank.data();
+    *n = simulation::fission_bank.size();
+    return 0;
+  }
+}
+
+//==============================================================================
+// Fortran compatibility
+//==============================================================================
+
+extern "C" void free_memory_bank()
+{
+  simulation::source_bank.clear();
+#pragma omp parallel
+  {
+    simulation::fission_bank.clear();
+  }
+#ifdef _OPENMP
+  simulation::master_fission_bank.clear();
+#endif
+}
+
+extern "C" int fission_bank_delayed_group(int64_t i) {
+  return simulation::fission_bank[i-1].delayed_group;
+}
+
+extern "C" double fission_bank_E(int64_t i) {
+  return simulation::fission_bank[i-1].E;
+}
+
+extern "C" double fission_bank_wgt(int64_t i) {
+  return simulation::fission_bank[i-1].wgt;
+}
+
+extern "C" void source_bank_xyz(int64_t i, double* xyz)
+{
+  xyz[0] = simulation::source_bank[i-1].xyz[0];
+  xyz[1] = simulation::source_bank[i-1].xyz[1];
+  xyz[2] = simulation::source_bank[i-1].xyz[2];
+}
+
+extern "C" double source_bank_E(int64_t i)
+{
+  return simulation::source_bank[i-1].E;
+}
+
+extern "C" double source_bank_wgt(int64_t i)
+{
+  return simulation::source_bank[i-1].wgt;
+}
+
+extern "C" void source_bank_set_wgt(int64_t i, double wgt)
+{
+  simulation::source_bank[i-1].wgt = wgt;
+}
+
+} // namespace openmc

src/bank_header.F90

@@ -2,8 +2,6 @@ module bank_header
 
   use, intrinsic :: ISO_C_BINDING
 
-  use error, only: E_ALLOCATE, set_errmsg
-
   implicit none
 
 !===============================================================================
@@ -21,70 +19,58 @@ module bank_header
     integer(C_INT) :: particle      ! particle type (neutron, photon, etc.)
   end type Bank
 
-  ! Source and fission bank
-  type(Bank), allocatable, target :: source_bank(:)
-  type(Bank), allocatable, target :: fission_bank(:)
-#ifdef _OPENMP
-  type(Bank), allocatable, target :: master_fission_bank(:)
-#endif
-
   integer(C_INT64_T), bind(C) :: n_bank       ! # of sites in fission bank
+!$omp threadprivate(n_bank)
 
-!$omp threadprivate(fission_bank, n_bank)
+  interface
+    function openmc_fission_bank(ptr, n) result(err) bind(C)
+      import C_PTR, C_INT64_T, C_INT
+      type(C_PTR), intent(out) :: ptr
+      integer(C_INT64_T), intent(out) :: n
+      integer(C_INT) :: err
+    end function
 
-contains
+    function fission_bank_delayed_group(i) result(g) bind(C)
+      import C_INT64_T, C_INT
+      integer(C_INT64_T), value :: i
+      integer(C_INT) :: g
+    end function
 
-!===============================================================================
-! FREE_MEMORY_BANK deallocates global arrays defined in this module
-!===============================================================================
+    function fission_bank_E(i) result(E) bind(C, name='fission_bank_E')
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: E
+    end function
 
-  subroutine free_memory_bank()
+    function fission_bank_wgt(i) result(wgt) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: wgt
+    end function
 
-    ! Deallocate fission and source bank and entropy
-!$omp parallel
-    if (allocated(fission_bank)) deallocate(fission_bank)
-!$omp end parallel
-#ifdef _OPENMP
-    if (allocated(master_fission_bank)) deallocate(master_fission_bank)
-#endif
-    if (allocated(source_bank)) deallocate(source_bank)
+    subroutine source_bank_xyz(i, xyz) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE), intent(in) :: xyz(*)
+    end subroutine
 
-  end subroutine free_memory_bank
+    function source_bank_E(i) result(E) bind(C, name='source_bank_E')
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: E
+    end function
 
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
+    function source_bank_wgt(i) result(wgt) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE) :: wgt
+    end function
 
-  function openmc_source_bank(ptr, n) result(err) bind(C)
-    ! Return a pointer to the source bank
-    type(C_PTR), intent(out) :: ptr
-    integer(C_INT64_T), intent(out) :: n
-    integer(C_INT) :: err
-
-    if (.not. allocated(source_bank)) then
-      err = E_ALLOCATE
-      call set_errmsg("Source bank has not been allocated.")
-    else
-      err = 0
-      ptr = C_LOC(source_bank)
-      n = size(source_bank)
-    end if
-  end function openmc_source_bank
-
-  function openmc_fission_bank(ptr, n) result(err) bind(C)
-    ! Return a pointer to the source bank
-    type(C_PTR), intent(out) :: ptr
-    integer(C_INT64_T), intent(out) :: n
-    integer(C_INT) :: err
-
-    if (.not. allocated(fission_bank)) then
-      err = E_ALLOCATE
-      call set_errmsg("Fission bank has not been allocated.")
-    else
-      err = 0
-      ptr = C_LOC(fission_bank)
-      n = size(fission_bank)
-    end if
-  end function openmc_fission_bank
+    subroutine source_bank_set_wgt(i, wgt) bind(C)
+      import C_INT64_T, C_DOUBLE
+      integer(C_INT64_T), value :: i
+      real(C_DOUBLE), value :: wgt
+    end subroutine
+  end interface
 
 end module bank_header

src/cmfd_execute.F90

@@ -218,7 +218,7 @@ contains
   subroutine cmfd_reweight(new_weights)
 
     use algorithm,   only: binary_search
-    use bank_header, only: source_bank
+    use bank_header
     use constants,   only: ZERO, ONE
     use error,       only: warning, fatal_error
     use message_passing
@@ -230,13 +230,14 @@ contains
     integer :: ny       ! maximum number of cells in y direction
     integer :: nz       ! maximum number of cells in z direction
     integer(C_INT) :: ng       ! maximum number of energy groups
-    integer :: i        ! iteration counter
+    integer(8) :: i        ! iteration counter
     integer :: g        ! index for group
     integer :: ijk(3)   ! spatial bin location
     integer :: e_bin    ! energy bin of source particle
     integer :: mesh_bin ! mesh bin of soruce particle
     integer :: n_groups ! number of energy groups
     real(8) :: norm     ! normalization factor
+    real(C_DOUBLE) :: xyz(3)
     logical(C_BOOL) :: outside  ! any source sites outside mesh
     logical :: in_mesh  ! source site is inside mesh
 
@@ -309,21 +310,22 @@ contains
     end if
 
     ! begin loop over source bank
-    do i = 1, int(work,4)
+    do i = 1, work
 
       ! Determine spatial bin
-      call cmfd_mesh % get_indices(source_bank(i) % xyz, ijk, in_mesh)
+      call source_bank_xyz(i, xyz)
+      call cmfd_mesh % get_indices(xyz, ijk, in_mesh)
 
       ! Determine energy bin
       n_groups = size(cmfd % egrid) - 1
-      if (source_bank(i) % E < cmfd % egrid(1)) then
+      if (source_bank_E(i) < cmfd % egrid(1)) then
         e_bin = 1
         if (master) call warning('Source pt below energy grid')
-      elseif (source_bank(i) % E > cmfd % egrid(n_groups + 1)) then
+      elseif (source_bank_E(i) > cmfd % egrid(n_groups + 1)) then
         e_bin = n_groups
         if (master) call warning('Source pt above energy grid')
       else
-        e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank(i) % E)
+        e_bin = binary_search(cmfd % egrid, n_groups + 1, source_bank_E(i))
       end if
 
       ! Reverese energy bin (lowest grp is highest energy bin)
@@ -335,8 +337,8 @@ contains
       end if
 
       ! Reweight particle
-      source_bank(i) % wgt = source_bank(i) % wgt * &
-           cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3))
+      call source_bank_set_wgt(i, source_bank_wgt(i) * &
+           cmfd % weightfactors(e_bin, ijk(1), ijk(2), ijk(3)))
     end do
 
   end subroutine cmfd_reweight

src/cmfd_execute.cpp

@@ -5,6 +5,7 @@
 #include "xtensor/xarray.hpp"
 #include "xtensor/xio.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/mesh.h"
 #include "openmc/message_passing.h"
@@ -18,14 +19,10 @@ extern "C" void
 cmfd_populate_sourcecounts(int n_energy, const double* energies,
   double* source_counts, bool* outside)
 {
-  // Get pointer to source bank
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-
   // Get source counts in each mesh bin / energy bin
   auto& m = model::meshes.at(settings::index_cmfd_mesh);
-  xt::xarray<double> counts = m->count_sites(simulation::work, source_bank, n_energy, energies, outside);
+  xt::xarray<double> counts = m->count_sites(simulation::work,
+    simulation::source_bank.data(), n_energy, energies, outside);
 
   // Copy data from the xarray into the source counts array
   std::copy(counts.begin(), counts.end(), source_counts);

src/eigenvalue.F90

@@ -2,9 +2,6 @@ module eigenvalue
 
   use, intrinsic :: ISO_C_BINDING
 
-  use bank_header
-  use simulation_header
-
   implicit none
 
   interface
@@ -15,49 +12,4 @@ module eigenvalue
     end function
   end interface
 
-contains
-
-#ifdef _OPENMP
-!===============================================================================
-! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission
-! bank that can be sampled. Note that this operation is necessarily sequential
-! to preserve the order of the bank when using varying numbers of threads.
-!===============================================================================
-
-  subroutine join_bank_from_threads() bind(C)
-
-    integer(8) :: total ! total number of fission bank sites
-    integer    :: i     ! loop index for threads
-
-    ! Initialize the total number of fission bank sites
-    total = 0
-
-!$omp parallel
-
-    ! Copy thread fission bank sites to one shared copy
-!$omp do ordered schedule(static)
-    do i = 1, n_threads
-!$omp ordered
-      master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
-      total = total + n_bank
-!$omp end ordered
-    end do
-!$omp end do
-
-    ! Make sure all threads have made it to this point
-!$omp barrier
-
-    ! Now copy the shared fission bank sites back to the master thread's copy.
-    if (thread_id == 0) then
-      n_bank = total
-      fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
-    else
-      n_bank = 0
-    end if
-
-!$omp end parallel
-
-  end subroutine join_bank_from_threads
-#endif
-
 end module eigenvalue

src/eigenvalue.cpp

@@ -5,6 +5,7 @@
 #include "xtensor/xtensor.hpp"
 #include "xtensor/xview.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/error.h"
@@ -69,13 +70,6 @@ void synchronize_bank()
 {
   simulation::time_bank.start();
 
-  // Get pointers to source/fission bank
-  Bank* source_bank;
-  Bank* fission_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-  openmc_fission_bank(&fission_bank, &n);
-
   // In order to properly understand the fission bank algorithm, you need to
   // think of the fission and source bank as being one global array divided
   // over multiple processors. At the start, each processor has a random amount
@@ -147,14 +141,14 @@ void synchronize_bank()
     // and the remaining 100 would be randomly sampled.
     if (total < settings::n_particles) {
       for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
-        temp_sites[index_temp] = fission_bank[i];
+        temp_sites[index_temp] = simulation::fission_bank[i];
         ++index_temp;
       }
     }
 
     // Randomly sample sites needed
     if (prn() < p_sample) {
-      temp_sites[index_temp] = fission_bank[i];
+      temp_sites[index_temp] = simulation::fission_bank[i];
       ++index_temp;
     }
   }
@@ -195,7 +189,8 @@ void synchronize_bank()
       // fission bank
       sites_needed = settings::n_particles - finish;
       for (int i = 0; i < sites_needed; ++i) {
-        temp_sites[index_temp] = fission_bank[simulation::n_bank - sites_needed + i];
+        int i_bank = simulation::n_bank - sites_needed + i;
+        temp_sites[index_temp] = simulation::fission_bank[i_bank];
         ++index_temp;
       }
     }
@@ -274,7 +269,7 @@ void synchronize_bank()
       // asynchronous receive for the source sites
 
       requests.emplace_back();
-      MPI_Irecv(&source_bank[index_local], static_cast<int>(n), mpi::bank,
+      MPI_Irecv(&simulation::source_bank[index_local], static_cast<int>(n), mpi::bank,
             neighbor, neighbor, mpi::intracomm, &requests.back());
 
     } else {
@@ -283,7 +278,7 @@ void synchronize_bank()
 
       index_temp = start - bank_position[mpi::rank];
       std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n],
-        &source_bank[index_local]);
+        &simulation::source_bank[index_local]);
     }
 
     // Increment all indices
@@ -300,7 +295,8 @@ void synchronize_bank()
   MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE);
 
 #else
-  std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, source_bank);
+  std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles,
+    simulation::source_bank.begin());
 #endif
 
   simulation::time_bank_sendrecv.stop();
@@ -347,6 +343,46 @@ void calculate_average_keff()
   }
 }
 
+#ifdef _OPENMP
+void join_bank_from_threads()
+{
+  // Initialize the total number of fission bank sites
+  int64_t total = 0;
+
+#pragma omp parallel
+  {
+    // Copy thread fission bank sites to one shared copy
+#pragma omp for ordered schedule(static)
+    for (int i = 0; i < simulation::n_threads; ++i) {
+#pragma omp ordered
+      {
+        std::copy(
+          &simulation::fission_bank[0],
+          &simulation::fission_bank[0] + simulation::n_bank,
+          &simulation::master_fission_bank[total]
+        );
+        total += simulation::n_bank;
+      }
+    }
+
+    // Make sure all threads have made it to this point
+#pragma omp barrier
+
+    // Now copy the shared fission bank sites back to the master thread's copy.
+    if (simulation::thread_id == 0) {
+      simulation::n_bank = total;
+      std::copy(
+        &simulation::master_fission_bank[0],
+        &simulation::master_fission_bank[0] + simulation::n_bank,
+        &simulation::fission_bank[0]
+      );
+    } else {
+      simulation::n_bank = 0;
+    }
+  }
+}
+#endif
+
 int openmc_get_keff(double* k_combined)
 {
   k_combined[0] = 0.0;
@@ -499,15 +535,10 @@ void shannon_entropy()
   // Get pointer to entropy mesh
   auto& m = model::meshes[settings::index_entropy_mesh];
 
-  // Get pointer to fission bank
-  Bank* fission_bank;
-  int64_t n;
-  openmc_fission_bank(&fission_bank, &n);
-
   // Get source weight in each mesh bin
   bool sites_outside;
-  xt::xtensor<double, 1> p = m->count_sites(
-    simulation::n_bank, fission_bank, 0, nullptr, &sites_outside);
+  xt::xtensor<double, 1> p = m->count_sites(simulation::n_bank,
+    simulation::fission_bank.data(), 0, nullptr, &sites_outside);
 
   // display warning message if there were sites outside entropy box
   if (sites_outside) {
@@ -544,15 +575,10 @@ void ufs_count_sites()
     s = m->volume_frac_;
 
   } else {
-    // Get pointer to source bank
-    Bank* source_bank;
-    int64_t n;
-    openmc_source_bank(&source_bank, &n);
-
     // count number of source sites in each ufs mesh cell
     bool sites_outside;
-    simulation::source_frac = m->count_sites(simulation::work, source_bank, 0, nullptr,
-      &sites_outside);
+    simulation::source_frac = m->count_sites(simulation::work,
+      simulation::source_bank.data(), 0, nullptr, &sites_outside);
 
     // Check for sites outside of the mesh
     if (mpi::master && sites_outside) {
@@ -572,7 +598,7 @@ void ufs_count_sites()
     // Since the total starting weight is not equal to n_particles, we need to
     // renormalize the weight of the source sites
     for (int i = 0; i < simulation::work; ++i) {
-      source_bank[i].wgt *= settings::n_particles / total;
+      simulation::source_bank[i].wgt *= settings::n_particles / total;
     }
   }
 }

src/endf.cpp

@@ -45,6 +45,37 @@ bool is_fission(int mt)
   return mt == 18 || mt == 19 || mt == 20 || mt == 21 || mt == 38;
 }
 
+bool is_disappearance(int mt)
+{
+  if (mt >= N_DISAPPEAR && mt <= N_DA) {
+    return true;
+  } else if (mt >= N_P0 && mt <= N_AC) {
+    return true;
+  } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE
+    || mt == N_3HEA || mt == N_3P) {
+    return true;
+  } else {
+    return false;
+  }
+}
+
+bool is_inelastic_scatter(int mt)
+{
+  if (mt < 100) {
+    if (is_fission(mt)) {
+      return false;
+    } else {
+      return mt >= MISC && mt != 27;
+    }
+  } else if (mt <= 200) {
+    return !is_disappearance(mt);
+  } else if (mt >= N_2N0 && mt <= N_2NC) {
+    return true;
+  } else {
+    return false;
+  }
+}
+
 //==============================================================================
 // Polynomial implementation
 //==============================================================================

src/input_xml.F90

@@ -2113,9 +2113,6 @@ contains
           call close_group(group_id)
           call file_close(file_id)
 
-          ! Assign resonant scattering data
-          if (res_scat_on) call nuclides(i_nuclide) % assign_0K_elastic_scattering()
-
           ! Determine if minimum/maximum energy for this nuclide is greater/less
           ! than the previous
           if (size(nuclides(i_nuclide) % grid) >= 1) then
@@ -2297,7 +2294,6 @@ contains
     logical :: file_exists                 ! Does multipole library exist?
     character(7) :: readable               ! Is multipole library readable?
     character(MAX_FILE_LEN) :: filename    ! Path to multipole xs library
-    character(kind=C_CHAR), pointer :: string(:)
     integer(HID_T) :: file_id
     integer(HID_T) :: group_id
 

src/material_header.F90

@@ -1041,4 +1041,21 @@ contains
 
   end subroutine bremsstrahlung_init
 
+!===============================================================================
+! Fortran compatibility
+!===============================================================================
+
+  function material_isotropic(i_material, i_nuc_mat) result(iso) bind(C)
+    integer(C_INT), value :: i_material
+    integer(C_INT), value :: i_nuc_mat
+    logical(C_BOOL) :: iso
+
+    iso = .false.
+    associate (mat => materials(i_material))
+      if (mat % has_isotropic_nuclides) then
+        iso = mat % p0(i_nuc_mat)
+      end if
+    end associate
+  end function
+
 end module material_header

src/nuclide.cpp

@@ -1,5 +1,17 @@
 #include "openmc/nuclide.h"
 
+#include "openmc/container_util.h"
+#include "openmc/endf.h"
+#include "openmc/error.h"
+#include "openmc/hdf5_interface.h"
+#include "openmc/message_passing.h"
+#include "openmc/search.h"
+#include "openmc/settings.h"
+#include "openmc/string_utils.h"
+
+#include <algorithm> // for sort
+#include <string> // for to_string, stoi
+
 namespace openmc {
 
 //==============================================================================
@@ -7,12 +19,348 @@ namespace openmc {
 //==============================================================================
 
 namespace data {
-
 std::array<double, 2> energy_min {0.0, 0.0};
 std::array<double, 2> energy_max {INFTY, INFTY};
-
+std::vector<std::unique_ptr<Nuclide>> nuclides;
 } // namespace data
 
+namespace simulation {
+NuclideMicroXS* micro_xs;
+MaterialMacroXS material_xs;
+} // namespace simulation
+
+//==============================================================================
+// Nuclide implementation
+//==============================================================================
+
+Nuclide::Nuclide(hid_t group, const double* temperature, int n)
+{
+  // Get name of nuclide from group, removing leading '/'
+  name_ = object_name(group).substr(1);
+
+  read_attribute(group, "Z", Z_);
+  read_attribute(group, "A", A_);
+  read_attribute(group, "metastable", metastable_);
+  read_attribute(group, "atomic_weight_ratio", awr_);
+
+  // Determine temperatures available
+  hid_t kT_group = open_group(group, "kTs");
+  auto dset_names = dataset_names(kT_group);
+  std::vector<double> temps_available;
+  for (const auto& name : dset_names) {
+    double T;
+    read_dataset(kT_group, name.c_str(), T);
+    temps_available.push_back(T / K_BOLTZMANN);
+  }
+  std::sort(temps_available.begin(), temps_available.end());
+
+  // If only one temperature is available, revert to nearest temperature
+  if (temps_available.size() == 1 && settings::temperature_method == TEMPERATURE_INTERPOLATION) {
+    if (mpi::master) {
+      warning("Cross sections for " + name_ + " are only available at one "
+        "temperature. Reverting to nearest temperature method.");
+    }
+    settings::temperature_method = TEMPERATURE_NEAREST;
+  }
+
+  // Determine actual temperatures to read -- start by checking whether a
+  // temperature range was given, in which case all temperatures in the range
+  // are loaded irrespective of what temperatures actually appear in the model
+  std::vector<int> temps_to_read;
+  double T_min = n > 0 ? settings::temperature_range[0] : 0.0;
+  double T_max = n > 0 ? settings::temperature_range[1] : INFTY;
+  if (T_max > 0.0) {
+    for (auto T : temps_available) {
+      if (T_min <= T && T <= T_max) {
+        temps_to_read.push_back(std::round(T));
+      }
+    }
+  }
+
+  switch (settings::temperature_method) {
+  case TEMPERATURE_NEAREST:
+    // Find nearest temperatures
+    for (int i = 0; i < n; ++i) {
+      double T_desired = temperature[i];
+
+      // Determine closest temperature
+      double min_delta_T = INFTY;
+      double T_actual;
+      for (auto T : temps_available) {
+        double delta_T = std::abs(T - T_desired);
+        if (delta_T < min_delta_T) {
+          T_actual = T;
+          min_delta_T = delta_T;
+        }
+      }
+
+      if (std::abs(T_actual - T_desired) < settings::temperature_tolerance) {
+        if (!contains(temps_to_read, std::round(T_actual))) {
+          temps_to_read.push_back(std::round(T_actual));
+
+          // Write warning for resonance scattering data if 0K is not available
+          if (std::abs(T_actual - T_desired) > 0 && T_desired == 0 && mpi::master) {
+            warning(name_ + " does not contain 0K data needed for resonance "
+              "scattering options selected. Using data at " + std::to_string(T_actual)
+              + " K instead.");
+          }
+        }
+      } else {
+        fatal_error("Nuclear data library does not contain cross sections for " +
+          name_ + " at or near " + std::to_string(T_desired) + " K.");
+      }
+    }
+    break;
+
+  case TEMPERATURE_INTERPOLATION:
+    // If temperature interpolation or multipole is selected, get a list of
+    // bounding temperatures for each actual temperature present in the model
+    for (int i = 0; i < n; ++i) {
+      double T_desired = temperature[i];
+
+      bool found_pair = false;
+      for (int j = 0; j < temps_available.size() - 1; ++j) {
+        if (temps_available[j] <= T_desired && T_desired < temps_available[j + 1]) {
+          int T_j = std::round(temps_available[j]);
+          int T_j1 = std::round(temps_available[j+1]);
+          if (!contains(temps_to_read, T_j)) {
+            temps_to_read.push_back(T_j);
+          }
+          if (!contains(temps_to_read, T_j1)) {
+            temps_to_read.push_back(T_j1);
+          }
+          found_pair = true;
+        }
+      }
+
+      if (!found_pair) {
+        fatal_error("Nuclear data library does not contain cross sections for " +
+          name_ +" at temperatures that bound " + std::to_string(T_desired) + " K.");
+      }
+    }
+    break;
+  }
+
+  // Sort temperatures to read
+  std::sort(temps_to_read.begin(), temps_to_read.end());
+
+  hid_t energy_group = open_group(group, "energy");
+  for (const auto& T : temps_to_read) {
+    std::string dset {std::to_string(T) + "K"};
+
+    // Determine exact kT values
+    double kT;
+    read_dataset(kT_group, dset.c_str(), kT);
+    kTs_.push_back(kT);
+
+    // Read energy grid
+    grid_.emplace_back();
+    read_dataset(energy_group, dset.c_str(), grid_.back().energy);
+  }
+  close_group(kT_group);
+
+  // Check for 0K energy grid
+  if (object_exists(energy_group, "0K")) {
+    read_dataset(energy_group, "0K", energy_0K_);
+  }
+  close_group(energy_group);
+
+  // Read reactions
+  hid_t rxs_group = open_group(group, "reactions");
+  for (auto name : group_names(rxs_group)) {
+    if (starts_with(name, "reaction_")) {
+      hid_t rx_group = open_group(rxs_group, name.c_str());
+      reactions_.push_back(std::make_unique<Reaction>(rx_group, temps_to_read));
+
+      // Check for 0K elastic scattering
+      const auto& rx = reactions_.back();
+      if (rx->mt_ == ELASTIC) {
+        if (object_exists(rx_group, "0K")) {
+          hid_t temp_group = open_group(rx_group, "0K");
+          read_dataset(temp_group, "xs", elastic_0K_);
+          close_group(temp_group);
+        }
+      }
+      close_group(rx_group);
+
+      // Determine reaction indices for inelastic scattering reactions
+      if (is_inelastic_scatter(rx->mt_) && !rx->redundant_) {
+        index_inelastic_scatter_.push_back(reactions_.size() - 1);
+      }
+    }
+  }
+  close_group(rxs_group);
+
+  // Check for nu-total
+  if (object_exists(group, "total_nu")) {
+    // Read total nu data
+    hid_t nu_group = open_group(group, "total_nu");
+    hid_t nu_dset = open_dataset(nu_group, "yield");
+    std::string func_type;
+    read_attribute(nu_dset, "type", func_type);
+    if (func_type == "Tabulated1D") {
+      total_nu_ = std::make_unique<Tabulated1D>(nu_dset);
+    } else if (func_type == "Polynomial") {
+      total_nu_ = std::make_unique<Polynomial>(nu_dset);
+    }
+    close_dataset(nu_dset);
+    close_group(nu_group);
+  }
+
+  this->create_derived();
+}
+
+void Nuclide::create_derived()
+{
+  for (int i = 0; i < reactions_.size(); ++i) {
+    const auto& rx {reactions_[i]};
+
+    for (int t = 0; t < kTs_.size(); ++t) {
+      // Skip redundant reactions
+      if (rx->redundant_) continue;
+
+      if (is_fission(rx->mt_)) {
+        fissionable_ = true;
+
+        // Keep track of fission reactions
+        if (t == 0) {
+          fission_rx_.push_back(rx.get());
+          if (rx->mt_ == N_F) has_partial_fission_ = true;
+        }
+      }
+    }
+  }
+
+  // Determine number of delayed neutron precursors
+  if (fissionable_) {
+    for (const auto& product : fission_rx_[0]->products_) {
+      if (product.emission_mode_ == EmissionMode::delayed) {
+        ++n_precursor_;
+      }
+    }
+  }
+
+  if (settings::res_scat_on) {
+    // Determine if this nuclide should be treated as a resonant scatterer
+    if (!settings::res_scat_nuclides.empty()) {
+      // If resonant nuclides were specified, check the list explicitly
+      for (const auto& name : settings::res_scat_nuclides) {
+        if (name_ == name) {
+          resonant_ = true;
+
+          // Make sure nuclide has 0K data
+          if (energy_0K_.empty()) {
+            fatal_error("Cannot treat " + name_ + " as a resonant scatterer "
+              "because 0 K elastic scattering data is not present.");
+          }
+          break;
+        }
+      }
+    } else {
+      // Otherwise, assume that any that have 0 K elastic scattering data are
+      // resonant
+      resonant_ = !energy_0K_.empty();
+    }
+
+    if (resonant_) {
+      // Build CDF for 0K elastic scattering
+      double xs_cdf_sum = 0.0;
+      xs_cdf_.resize(energy_0K_.size());
+      xs_cdf_[0] = 0.0;
+
+      const auto& E = energy_0K_;
+      auto& xs = elastic_0K_;
+      for (int i = 0; i < E.size() - 1; ++i) {
+        // Negative cross sections result in a CDF that is not monotonically
+        // increasing. Set all negative xs values to zero.
+        if (xs[i] < 0.0) xs[i] = 0.0;
+
+        // build xs cdf
+        xs_cdf_sum += (std::sqrt(E[i])*xs[i] + std::sqrt(E[i+1])*xs[i+1])
+              / 2.0 * (E[i+1] - E[i]);
+        xs_cdf_[i] = xs_cdf_sum;
+      }
+    }
+  }
+}
+
+double Nuclide::nu(double E, EmissionMode mode, int group) const
+{
+  if (!fissionable_) return 0.0;
+
+  switch (mode) {
+  case EmissionMode::prompt:
+    return (*fission_rx_[0]->products_[0].yield_)(E);
+  case EmissionMode::delayed:
+    if (n_precursor_ > 0) {
+      auto rx = fission_rx_[0];
+      if (group >= 1 && group < rx->products_.size()) {
+        // If delayed group specified, determine yield immediately
+        return (*rx->products_[group].yield_)(E);
+      } else {
+        double nu {0.0};
+
+        for (int i = 1; i < rx->products_.size(); ++i) {
+          // Skip any non-neutron products
+          const auto& product = rx->products_[i];
+          if (product.particle_ != ParticleType::neutron) continue;
+
+          // Evaluate yield
+          if (product.emission_mode_ == EmissionMode::delayed) {
+            nu += (*product.yield_)(E);
+          }
+        }
+        return nu;
+      }
+    } else {
+      return 0.0;
+    }
+  case EmissionMode::total:
+    if (total_nu_) {
+      return (*total_nu_)(E);
+    } else {
+      return (*fission_rx_[0]->products_[0].yield_)(E);
+    }
+  }
+}
+
+void Nuclide::calculate_elastic_xs(int i_nuclide) const
+{
+  // Get temperature index, grid index, and interpolation factor
+  auto& micro = simulation::micro_xs[i_nuclide];
+  int i_temp = micro.index_temp - 1;
+  int i_grid = micro.index_grid - 1;
+  double f = micro.interp_factor;
+
+  if (i_temp >= 0) {
+    const auto& xs = reactions_[0]->xs_[i_temp].value;
+    micro.elastic = (1.0 - f)*xs[i_grid] + f*xs[i_grid + 1];
+  }
+}
+
+double Nuclide::elastic_xs_0K(double E) const
+{
+  // Determine index on nuclide energy grid
+  int i_grid;
+  if (E < energy_0K_.front()) {
+    i_grid = 0;
+  } else if (E > energy_0K_.back()) {
+    i_grid = energy_0K_.size() - 2;
+  } else {
+    i_grid = lower_bound_index(energy_0K_.begin(), energy_0K_.end(), E);
+  }
+
+  // check for rare case where two energy points are the same
+  if (energy_0K_[i_grid] == energy_0K_[i_grid+1]) ++i_grid;
+
+  // calculate interpolation factor
+  double f = (E - energy_0K_[i_grid]) /
+    (energy_0K_[i_grid + 1] - energy_0K_[i_grid]);
+
+  // Calculate microscopic nuclide elastic cross section
+  return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1];
+}
+
 //==============================================================================
 // Fortran compatibility functions
 //==============================================================================
@@ -24,4 +372,27 @@ set_particle_energy_bounds(int particle, double E_min, double E_max)
   data::energy_max[particle - 1] = E_max;
 }
 
+extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature, int n)
+{
+  data::nuclides.push_back(std::make_unique<Nuclide>(group, temperature, n));
+  return data::nuclides.back().get();
+}
+
+extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx)
+{
+  return nuc->reactions_[i_rx-1].get();
+}
+
+extern "C" void nuclides_clear() { data::nuclides.clear(); }
+
+
+extern "C" NuclideMicroXS* micro_xs_ptr();
+void set_micro_xs()
+{
+#pragma omp parallel
+  {
+    simulation::micro_xs = micro_xs_ptr();
+  }
+}
+
 } // namespace openmc

src/nuclide_header.F90

@@ -71,15 +71,7 @@ module nuclide_header
     ! Microscopic cross sections
     type(SumXS), allocatable :: xs(:)
 
-    ! Resonance scattering info
-    logical              :: resonant = .false. ! resonant scatterer?
-    real(8), allocatable :: energy_0K(:)  ! energy grid for 0K xs
-    real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section
-    real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section
-
     ! Fission information
-    logical :: has_partial_fission = .false. ! nuclide has partial fission reactions?
-    integer :: n_fission = 0                 ! # of fission reactions
     integer :: n_precursor = 0               ! # of delayed neutron precursors
     integer, allocatable :: index_fission(:) ! indices in reactions
     class(Function1D), allocatable :: total_nu
@@ -95,7 +87,6 @@ module nuclide_header
 
     ! Reactions
     type(Reaction), allocatable :: reactions(:)
-    integer, allocatable :: index_inelastic_scatter(:)
 
     ! Array that maps MT values to index in reactions; used at tally-time. Note
     ! that ENDF-102 does not assign any MT values above 891.
@@ -105,8 +96,9 @@ module nuclide_header
     class(Function1D), allocatable :: fission_q_prompt ! fragments and prompt neutrons, gammas
     class(Function1D), allocatable :: fission_q_recov  ! fragments, neutrons, gammas, betas
 
+    type(C_PTR) :: ptr
+
   contains
-    procedure :: assign_0K_elastic_scattering
     procedure :: clear => nuclide_clear
     procedure :: from_hdf5 => nuclide_from_hdf5
     procedure :: init_grid => nuclide_init_grid
@@ -184,7 +176,7 @@ module nuclide_header
 
   ! Cross section caches
   type(NuclideMicroXS), allocatable, target :: micro_xs(:)  ! Cache for each nuclide
-  type(MaterialMacroXS)             :: material_xs  ! Cache for current material
+  type(MaterialMacroXS), bind(C)            :: material_xs  ! Cache for current material
 !$omp threadprivate(micro_xs, material_xs)
 
   ! Minimum/maximum energies
@@ -231,78 +223,10 @@ contains
     b = library_present_c(type, to_c_string(name))
   end function
 
-!===============================================================================
-! ASSIGN_0K_ELASTIC_SCATTERING
-!===============================================================================
-
-  subroutine assign_0K_elastic_scattering(this)
-    class(Nuclide), intent(inout) :: this
-
-    integer :: i
-    real(8) :: xs_cdf_sum
-
-    interface
-      function res_scat_nuclides_empty() result(empty) bind(C)
-        import C_BOOL
-        logical(C_BOOL) :: empty
-      end function
-
-      function res_scat_nuclides_size() result(n) bind(C)
-        import C_INT
-        integer(C_INT) :: n
-      end function
-
-      function res_scat_nuclides_cmp(i, name) result(b) bind(C)
-        import C_INT, C_CHAR, C_BOOL
-        integer(C_INT), value :: i
-        character(kind=C_CHAR), intent(in) :: name(*)
-        logical(C_BOOL) :: b
-      end function
-    end interface
-
-    this % resonant = .false.
-    if (.not. res_scat_nuclides_empty()) then
-      ! If resonant nuclides were specified, check the list explicitly
-      do i = 1, res_scat_nuclides_size()
-        if (res_scat_nuclides_cmp(i, to_c_string(this % name))) then
-          this % resonant = .true.
-
-          ! Make sure nuclide has 0K data
-          if (.not. allocated(this % energy_0K)) then
-            call fatal_error("Cannot treat " // trim(this % name) // " as a &
-                 &resonant scatterer because 0 K elastic scattering data is &
-                 &not present.")
-          end if
-
-          exit
-        end if
-      end do
-    else
-      ! Otherwise, assume that any that have 0 K elastic scattering data are
-      ! resonant
-      this % resonant = allocated(this % energy_0K)
-    end if
-
-    if (this % resonant) then
-      ! Build CDF for 0K elastic scattering
-      xs_cdf_sum = ZERO
-      allocate(this % xs_cdf(0:size(this % energy_0K)))
-      this % xs_cdf(0) = ZERO
-
-      associate (E => this % energy_0K, xs => this % elastic_0K)
-        do i = 1, size(E) - 1
-          ! Negative cross sections result in a CDF that is not monotonically
-          ! increasing. Set all negative xs values to zero.
-          if (xs(i) < ZERO) xs(i) = ZERO
-
-          ! build xs cdf
-          xs_cdf_sum = xs_cdf_sum + (sqrt(E(i))*xs(i) + sqrt(E(i+1))*xs(i+1))&
-               / TWO * (E(i+1) - E(i))
-          this % xs_cdf(i) = xs_cdf_sum
-        end do
-      end associate
-    end if
-  end subroutine assign_0K_elastic_scattering
+  function micro_xs_ptr() result(ptr) bind(C)
+    type(C_PTR) :: ptr
+    ptr = C_LOC(micro_xs(1))
+  end function
 
 !===============================================================================
 ! NUCLIDE_CLEAR resets and deallocates data in Nuclide
@@ -310,19 +234,6 @@ contains
 
   subroutine nuclide_clear(this)
     class(Nuclide), intent(inout) :: this ! The Nuclide object to clear
-    integer :: i
-
-    interface
-      subroutine reaction_delete(rx) bind(C)
-        import C_PTR
-        type(C_PTR), value :: rx
-      end subroutine reaction_delete
-    end interface
-
-    do i = 1, size(this % reactions)
-      call reaction_delete(this % reactions(i) % ptr)
-    end do
-    deallocate(this % reactions)
 
     if (associated(this % multipole)) deallocate(this % multipole)
 
@@ -332,7 +243,7 @@ contains
                                minmax, master, i_nuclide)
     class(Nuclide),   intent(inout) :: this
     integer(HID_T),   intent(in)    :: group_id
-    type(VectorReal), intent(in)    :: temperature ! list of desired temperatures
+    type(VectorReal), intent(in), target    :: temperature ! list of desired temperatures
     integer,          intent(inout) :: method
     real(8),          intent(in)    :: tolerance
     real(8),          intent(in)    :: minmax(2)  ! range of temperatures
@@ -347,7 +258,6 @@ contains
     integer(HID_T) :: kT_group
     integer(HID_T) :: rxs_group
     integer(HID_T) :: rx_group
-    integer(HID_T) :: xs, temp_group
     integer(HID_T) :: total_nu
     integer(HID_T) :: fer_group                 ! fission_energy_release group
     integer(HID_T) :: fer_dset
@@ -361,7 +271,24 @@ contains
     real(8) :: temp_actual
     type(VectorInt) :: MTs
     type(VectorInt) :: temps_to_read
-    type(VectorInt) :: index_inelastic_scatter
+
+    interface
+      function nuclide_from_hdf5_c(group, temperature, n) result(ptr) bind(C)
+        import HID_T, C_DOUBLE, C_INT, C_PTR
+        integer(HID_T), value :: group
+        type(C_PTR), value :: temperature
+        integer(C_INT), value :: n
+        type(C_PTR) :: ptr
+      end function
+    end interface
+
+    ! Read data on C++ side
+    if (temperature % size() > 0) then
+      this % ptr = nuclide_from_hdf5_c(group_id, C_LOC(temperature % data(1)), &
+           temperature % size())
+    else
+      this % ptr = nuclide_from_hdf5_c(group_id, C_NULL_PTR, 0)
+    end if
 
     ! Get name of nuclide from group
     this % name = get_name(group_id)
@@ -487,15 +414,6 @@ contains
       call read_dataset(this % grid(i) % energy, energy_dset)
       call close_dataset(energy_dset)
     end do
-
-    ! Check for 0K energy grid
-    if (object_exists(energy_group, '0K')) then
-      energy_dset = open_dataset(energy_group, '0K')
-      call get_shape(energy_dset, dims)
-      allocate(this % energy_0K(int(dims(1), 4)))
-      call read_dataset(this % energy_0K, energy_dset)
-      call close_dataset(energy_dset)
-    end if
     call close_group(energy_group)
 
     ! Get MT values based on group names
@@ -513,36 +431,13 @@ contains
       rx_group = open_group(rxs_group, 'reaction_' // trim(&
            zero_padded(MTs % data(i), 3)))
 
-      call this % reactions(i) % from_hdf5(rx_group, temps_to_read)
-
-      ! Check for 0K elastic scattering
-      if (this % reactions(i) % MT == 2) then
-        if (object_exists(rx_group, '0K')) then
-          temp_group = open_group(rx_group, '0K')
-          xs = open_dataset(temp_group, 'xs')
-          call get_shape(xs, dims)
-          allocate(this % elastic_0K(int(dims(1), 4)))
-          call read_dataset(this % elastic_0K, xs)
-          call close_dataset(xs)
-          call close_group(temp_group)
-        end if
-      end if
-
-      ! Add the reaction index to the scattering array if this is an inelastic
-      ! scatter reaction
-      if (is_inelastic_scatter(MTs % data(i))) then
-        call index_inelastic_scatter % push_back(i)
-      end if
+      ! Set pointer for each reaction
+      call this % reactions(i) % init(this % ptr, i)
 
       call close_group(rx_group)
     end do
     call close_group(rxs_group)
 
-    ! Recast to a regular array to save space
-    allocate(this % index_inelastic_scatter(index_inelastic_scatter % size()))
-    this % index_inelastic_scatter = &
-         index_inelastic_scatter % data(1: index_inelastic_scatter % size())
-
     ! Read unresolved resonance probability tables if present
     if (object_exists(group_id, 'urr')) then
       this % urr_present = .true.
@@ -718,7 +613,6 @@ contains
               allocate(this % index_fission(1))
             elseif (rx % MT == N_F) then
               allocate(this % index_fission(PARTIAL_FISSION_MAX))
-              this % has_partial_fission = .true.
             end if
           end if
 
@@ -738,7 +632,6 @@ contains
             if (t == 1) then
               i_fission = i_fission + 1
               this % index_fission(i_fission) = i
-              this % n_fission = this % n_fission + 1
             end if
           end if  ! fission
         end do  ! temperature
@@ -1353,45 +1246,6 @@ contains
     end if
   end subroutine multipole_deriv_eval
 
-!===============================================================================
-! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
-! for a given nuclide at the trial relative energy used in resonance scattering
-!===============================================================================
-
-  pure function elastic_xs_0K(E, nuc) result(xs_out)
-    real(8),       intent(in) :: E      ! trial energy
-    type(Nuclide), intent(in) :: nuc    ! target nuclide at temperature
-    real(8)                   :: xs_out ! 0K xs at trial energy
-
-    integer :: i_grid ! index on nuclide energy grid
-    integer :: n_grid
-    real(8) :: f      ! interp factor on nuclide energy grid
-
-    ! Determine index on nuclide energy grid
-    n_grid = size(nuc % energy_0K)
-    if (E < nuc % energy_0K(1)) then
-      i_grid = 1
-    elseif (E > nuc % energy_0K(n_grid)) then
-      i_grid = n_grid - 1
-    else
-      i_grid = binary_search(nuc % energy_0K, n_grid, E)
-    end if
-
-    ! check for rare case where two energy points are the same
-    if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then
-      i_grid = i_grid + 1
-    end if
-
-    ! calculate interpolation factor
-    f = (E - nuc % energy_0K(i_grid)) &
-         & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid))
-
-    ! Calculate microscopic nuclide elastic cross section
-    xs_out = (ONE - f) * nuc % elastic_0K(i_grid) &
-         & + f * nuc % elastic_0K(i_grid + 1)
-
-  end function elastic_xs_0K
-
 !===============================================================================
 ! CALCULATE_URR_XS determines cross sections in the unresolved resonance range
 ! from probability tables
@@ -1573,6 +1427,9 @@ contains
     interface
       subroutine library_clear() bind(C)
       end subroutine
+
+      subroutine nuclides_clear() bind(C)
+      end subroutine
     end interface
 
     ! Deallocate cross section data, listings, and cache
@@ -1582,6 +1439,7 @@ contains
         call nuclides(i) % clear()
       end do
       deallocate(nuclides)
+      call nuclides_clear()
     end if
     n_nuclides = 0
 
@@ -1666,9 +1524,6 @@ contains
         call nuclide_dict % set(to_lower(name_), n)
         n_nuclides = n
 
-        ! Assign resonant scattering data
-        if (res_scat_on) call nuclides(n) % assign_0K_elastic_scattering()
-
         ! Initialize nuclide grid
         call nuclides(n) % init_grid(energy_min(NEUTRON), &
              energy_max(NEUTRON), n_log_bins)
@@ -1706,4 +1561,23 @@ contains
     end if
   end function openmc_nuclide_name
 
+  function nuclide_wmp_present(i_nuclide) result(b) bind(C)
+    integer(C_INT), value :: i_nuclide
+    logical(C_BOOL) :: b
+    b = nuclides(i_nuclide + 1) % mp_present
+  end function
+
+  function nuclide_wmp_emin(i_nuclide) result(E) bind(C)
+    integer(C_INT), value :: i_nuclide
+    real(C_DOUBLE) :: E
+    E = nuclides(i_nuclide + 1) % multipole % E_min
+  end function
+
+  function nuclide_wmp_emax(i_nuclide) result(E) bind(C)
+    integer(C_INT), value :: i_nuclide
+    real(C_DOUBLE) :: E
+    E = nuclides(i_nuclide + 1) % multipole % E_max
+  end function
+
+
 end module nuclide_header

src/particle.cpp

@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <sstream>
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/error.h"
@@ -184,17 +185,12 @@ Particle::write_restart() const
     write_dataset(file_id, "id", id);
     write_dataset(file_id, "type", type);
 
-    // Get pointer to source bank
-    Bank* src;
-    int64_t n;
-    openmc_source_bank(&src, &n);
-
     int64_t i = simulation::current_work;
-    write_dataset(file_id, "weight", src[i-1].wgt);
-    write_dataset(file_id, "energy", src[i-1].E);
+    write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt);
+    write_dataset(file_id, "energy", simulation::source_bank[i-1].E);
     hsize_t dims[] {3};
-    write_double(file_id, 1, dims, "xyz", src[i-1].xyz, false);
-    write_double(file_id, 1, dims, "uvw", src[i-1].uvw, false);
+    write_double(file_id, 1, dims, "xyz", simulation::source_bank[i-1].xyz, false);
+    write_double(file_id, 1, dims, "uvw", simulation::source_bank[i-1].uvw, false);
 
     // Close file
     file_close(file_id);

src/particle_restart.F90

@@ -33,12 +33,20 @@ contains
     integer :: previous_run_mode
     type(Particle) :: p
 
+    interface
+      subroutine set_micro_xs() bind(C)
+      end subroutine
+    end interface
+
     err = 0
 
     ! Set verbosity high
     verbosity = 10
 
+  !$omp parallel
     allocate(micro_xs(n_nuclides))
+  !$omp end parallel
+    call set_micro_xs()
 
     ! Initialize the particle to be tracked
     call particle_initialize(p)

src/photon.cpp

@@ -0,0 +1,38 @@
+#include "openmc/photon.h"
+
+#include "openmc/hdf5_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace data {
+
+std::vector<PhotonInteraction> elements;
+
+} // namespace data
+
+//==============================================================================
+// PhotonInteraction implementation
+//==============================================================================
+
+PhotonInteraction::PhotonInteraction(hid_t group)
+{
+  // Get name of nuclide from group, removing leading '/'
+  name_ = object_name(group).substr(1);
+
+  read_attribute(group, "Z", Z_);
+}
+
+//==============================================================================
+// Fortran compatibility
+//==============================================================================
+
+extern "C" void photon_from_hdf5_c(hid_t group)
+{
+  data::elements.emplace_back(group);
+}
+
+} // namespace openmc

src/photon_header.F90

@@ -130,6 +130,16 @@ contains
     real(8), allocatable :: matrix(:,:)
     real(8), allocatable :: dcs(:,:)
 
+    interface
+      subroutine photon_from_hdf5_c(group) bind(C)
+        import HID_T
+        integer(HID_T), value :: group
+      end subroutine
+    end interface
+
+    ! Read element data on C++ side
+    call photon_from_hdf5_c(group_id)
+
     ! Get name of nuclide from group
     this % name = get_name(group_id)
 

src/photon_physics.F90

@@ -540,9 +540,9 @@ contains
 ! THICK_TARGET_BREMSSTRAHLUNG
 !===============================================================================
 
-  subroutine thick_target_bremsstrahlung(p, E_lost)
+  subroutine thick_target_bremsstrahlung(p, E_lost) bind(C)
     type(Particle), intent(inout) :: p
-    real(8),        intent(inout) :: E_lost
+    real(C_DOUBLE), intent(out) :: E_lost
 
     integer :: i, j
     integer :: i_e, i_w

src/physics_mg.cpp

@@ -5,6 +5,7 @@
 
 #include "xtensor/xarray.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
@@ -47,12 +48,8 @@ sample_reaction(Particle* p, const MaterialMacroXS* material_xs)
 
   if (model::materials[p->material - 1]->fissionable) {
     if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      Bank* result_bank;
-      int64_t result_bank_size;
-      // Get pointer to fission bank from Fortran side
-      openmc_fission_bank(&result_bank, &result_bank_size);
-      create_fission_sites(p, result_bank, &simulation::n_bank, result_bank_size,
-                           material_xs);
+      create_fission_sites(p, simulation::fission_bank.data(), &simulation::n_bank,
+                           simulation::fission_bank.size(), material_xs);
     } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
                (settings::create_fission_neutrons)) {
       create_fission_sites(p, p->secondary_bank, &(p->n_secondary),

src/position.cpp

@@ -60,4 +60,22 @@ Position::operator*=(double v)
   return *this;
 }
 
+Position&
+Position::operator/=(Position other)
+{
+  x /= other.x;
+  y /= other.y;
+  z /= other.z;
+  return *this;
+}
+
+Position&
+Position::operator/=(double v)
+{
+  x /= v;
+  y /= v;
+  z /= v;
+  return *this;
+}
+
 } // namespace openmc

src/reaction.cpp

@@ -3,6 +3,7 @@
 #include <string>
 #include <utility> // for move
 
+#include "openmc/constants.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/endf.h"
 #include "openmc/random_lcg.h"
@@ -10,6 +11,10 @@
 
 namespace openmc {
 
+//==============================================================================
+// Reaction implementation
+//==============================================================================
+
 Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
 {
   read_attribute(group, "Q_value", q_value_);
@@ -77,16 +82,198 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
 }
 
 //==============================================================================
-// Fortran compatibility functions
+// Non-member functions
 //==============================================================================
 
-Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n)
+std::string reaction_name(int mt)
 {
-  std::vector<int> temps {temperatures, temperatures + n};
-  return new Reaction{group, temps};
+  if (mt == SCORE_FLUX) {
+    return  "flux";
+  } else if (mt == SCORE_TOTAL) {
+    return  "total";
+  } else if (mt == SCORE_SCATTER) {
+    return  "scatter";
+  } else if (mt == SCORE_NU_SCATTER) {
+    return  "nu-scatter";
+  } else if (mt == SCORE_ABSORPTION) {
+    return  "absorption";
+  } else if (mt == SCORE_FISSION) {
+    return  "fission";
+  } else if (mt == SCORE_NU_FISSION) {
+    return  "nu-fission";
+  } else if (mt == SCORE_DECAY_RATE) {
+    return  "decay-rate";
+  } else if (mt == SCORE_DELAYED_NU_FISSION) {
+    return  "delayed-nu-fission";
+  } else if (mt == SCORE_PROMPT_NU_FISSION) {
+    return  "prompt-nu-fission";
+  } else if (mt == SCORE_KAPPA_FISSION) {
+    return  "kappa-fission";
+  } else if (mt == SCORE_CURRENT) {
+    return  "current";
+  } else if (mt == SCORE_EVENTS) {
+    return  "events";
+  } else if (mt == SCORE_INVERSE_VELOCITY) {
+    return  "inverse-velocity";
+  } else if (mt == SCORE_FISS_Q_PROMPT) {
+    return  "fission-q-prompt";
+  } else if (mt == SCORE_FISS_Q_RECOV) {
+    return  "fission-q-recoverable";
+
+  // Normal ENDF-based reactions
+  } else if (mt == TOTAL_XS) {
+    return  "(n,total)";
+  } else if (mt == ELASTIC) {
+    return  "(n,elastic)";
+  } else if (mt == N_LEVEL) {
+    return  "(n,level)";
+  } else if (mt == N_2ND) {
+    return  "(n,2nd)";
+  } else if (mt == N_2N) {
+    return  "(n,2n)";
+  } else if (mt == N_3N) {
+    return  "(n,3n)";
+  } else if (mt == N_FISSION) {
+    return  "(n,fission)";
+  } else if (mt == N_F) {
+    return  "(n,f)";
+  } else if (mt == N_NF) {
+    return  "(n,nf)";
+  } else if (mt == N_2NF) {
+    return  "(n,2nf)";
+  } else if (mt == N_NA) {
+    return  "(n,na)";
+  } else if (mt == N_N3A) {
+    return  "(n,n3a)";
+  } else if (mt == N_2NA) {
+    return  "(n,2na)";
+  } else if (mt == N_3NA) {
+    return  "(n,3na)";
+  } else if (mt == N_NP) {
+    return  "(n,np)";
+  } else if (mt == N_N2A) {
+    return  "(n,n2a)";
+  } else if (mt == N_2N2A) {
+    return  "(n,2n2a)";
+  } else if (mt == N_ND) {
+    return  "(n,nd)";
+  } else if (mt == N_NT) {
+    return  "(n,nt)";
+  } else if (mt == N_N3HE) {
+    return  "(n,nHe-3)";
+  } else if (mt == N_ND2A) {
+    return  "(n,nd2a)";
+  } else if (mt == N_NT2A) {
+    return  "(n,nt2a)";
+  } else if (mt == N_4N) {
+    return  "(n,4n)";
+  } else if (mt == N_3NF) {
+    return  "(n,3nf)";
+  } else if (mt == N_2NP) {
+    return  "(n,2np)";
+  } else if (mt == N_3NP) {
+    return  "(n,3np)";
+  } else if (mt == N_N2P) {
+    return  "(n,n2p)";
+  } else if (mt == N_NPA) {
+    return  "(n,npa)";
+  } else if (N_N1 <= mt && mt <= N_N40) {
+    return  "(n,n" + std::to_string(mt-50) + ")";
+  } else if (mt == N_NC) {
+    return  "(n,nc)";
+  } else if (mt == N_DISAPPEAR) {
+    return  "(n,disappear)";
+  } else if (mt == N_GAMMA) {
+    return  "(n,gamma)";
+  } else if (mt == N_P) {
+    return  "(n,p)";
+  } else if (mt == N_D) {
+    return  "(n,d)";
+  } else if (mt == N_T) {
+    return  "(n,t)";
+  } else if (mt == N_3HE) {
+    return  "(n,3He)";
+  } else if (mt == N_A) {
+    return  "(n,a)";
+  } else if (mt == N_2A) {
+    return  "(n,2a)";
+  } else if (mt == N_3A) {
+    return  "(n,3a)";
+  } else if (mt == N_2P) {
+    return  "(n,2p)";
+  } else if (mt == N_PA) {
+    return  "(n,pa)";
+  } else if (mt == N_T2A) {
+    return  "(n,t2a)";
+  } else if (mt == N_D2A) {
+    return  "(n,d2a)";
+  } else if (mt == N_PD) {
+    return  "(n,pd)";
+  } else if (mt == N_PT) {
+    return  "(n,pt)";
+  } else if (mt == N_DA) {
+    return  "(n,da)";
+  } else if (mt == 201) {
+    return  "(n,Xn)";
+  } else if (mt == 202) {
+    return  "(n,Xgamma)";
+  } else if (mt == 203) {
+    return  "(n,Xp)";
+  } else if (mt == 204) {
+    return  "(n,Xd)";
+  } else if (mt == 205) {
+    return  "(n,Xt)";
+  } else if (mt == 206) {
+    return  "(n,X3He)";
+  } else if (mt == 207) {
+    return  "(n,Xa)";
+  } else if (mt == 444) {
+    return  "(damage)";
+  } else if (mt == COHERENT) {
+    return  "coherent scatter";
+  } else if (mt == INCOHERENT) {
+    return  "incoherent scatter";
+  } else if (mt == PAIR_PROD_ELEC) {
+    return  "pair production, electron";
+  } else if (mt == PAIR_PROD) {
+    return  "pair production";
+  } else if (mt == PAIR_PROD_NUC) {
+    return  "pair production, nuclear";
+  } else if (mt == PHOTOELECTRIC) {
+    return  "photoelectric";
+  } else if (534 <= mt && mt <= 572) {
+    std::stringstream name;
+    name << "photoelectric, " << SUBSHELLS[mt - 534] << " subshell";
+    return name.str();
+  } else if (600 <= mt && mt <= 648) {
+    return  "(n,p" + std::to_string(mt-600) + ")";
+  } else if (mt == 649) {
+    return  "(n,pc)";
+  } else if (650 <= mt && mt <= 698) {
+    return  "(n,d" + std::to_string(mt-650) + ")";
+  } else if (mt == 699) {
+    return  "(n,dc)";
+  } else if (700 <= mt && mt <= 748) {
+    return  "(n,t" + std::to_string(mt-700) + ")";
+  } else if (mt == 749) {
+    return  "(n,tc)";
+  } else if (750 <= mt && mt <= 798) {
+    return  "(n,3He" + std::to_string(mt-750) + ")";
+  } else if (mt == 799) {
+    return  "(n,3Hec)";
+  } else if (800 <= mt && mt <= 848) {
+    return  "(n,a" + std::to_string(mt-800) + ")";
+  } else if (mt == 849) {
+    return  "(n,ac)";
+  } else {
+    return  "MT=" + std::to_string(mt);
+  }
 }
 
-void reaction_delete(Reaction* rx) { delete rx; }
+
+//==============================================================================
+// Fortran compatibility functions
+//==============================================================================
 
 int reaction_mt(Reaction* rx) { return rx->mt_; }
 

src/reaction_header.F90

@@ -22,7 +22,7 @@ module reaction_header
     logical(C_BOOL) :: scatter_in_cm           ! scattering system in center-of-mass?
     logical(C_BOOL) :: redundant               ! redundant reaction?
   contains
-    procedure :: from_hdf5
+    procedure :: init
     procedure :: mt_
     procedure :: q_value_
     procedure :: scatter_in_cm_
@@ -40,11 +40,10 @@ module reaction_header
   end type Reaction
 
   interface
-    function reaction_from_hdf5(group, temperatures, n) result(ptr) bind(C)
-      import C_PTR, HID_T, C_INT
-      integer(HID_T), value :: group
-      integer(C_INT), intent(in) :: temperatures
-      integer(C_INT), value :: n
+    function nuclide_reaction(nuc_ptr, i_rx) result(ptr) bind(C)
+      import C_PTR, C_INT
+      type(C_PTR), value :: nuc_ptr
+      integer(C_INT), value :: i_rx
       type(C_PTR) :: ptr
     end function
 
@@ -148,27 +147,17 @@ module reaction_header
 
 contains
 
-  subroutine from_hdf5(this, group_id, temperatures)
+  subroutine init(this, nuc_ptr, i_rx)
     class(Reaction), intent(inout) :: this
-    integer(HID_T),  intent(in)    :: group_id
-    type(VectorInt), intent(in)    :: temperatures
+    type(C_PTR), intent(in) :: nuc_ptr
+    integer(C_INT), intent(in) :: i_rx
 
-    integer(C_INT) :: dummy
-    integer(C_INT) :: n
-
-    n = temperatures % size()
-    if (n > 0) then
-      this % ptr = reaction_from_hdf5(group_id, temperatures % data(1), n)
-    else
-      ! In this case, temperatures % data(1) doesn't exist, so we just pass a
-      ! dummy value
-      this % ptr = reaction_from_hdf5(group_id, dummy, n)
-    end if
+    this % ptr = nuclide_reaction(nuc_ptr, i_rx)
     this % MT = reaction_mt(this % ptr)
     this % Q_value = reaction_q_value(this % ptr)
     this % scatter_in_cm = reaction_scatter_in_cm(this % ptr)
     this % redundant = reaction_redundant(this % ptr)
-  end subroutine from_hdf5
+  end subroutine
 
   function mt_(this) result(mt)
     class(Reaction), intent(in) :: this

src/settings.cpp

@@ -806,18 +806,6 @@ void read_settings_xml()
 //==============================================================================
 
 extern "C" {
-  bool res_scat_nuclides_empty() {
-    return settings::res_scat_nuclides.empty();
-  }
-
-  int res_scat_nuclides_size() {
-    return settings::res_scat_nuclides.size();
-  }
-
-  bool res_scat_nuclides_cmp(int i, const char* name) {
-    return settings::res_scat_nuclides[i - 1] == name;
-  }
-
   const char* path_cross_sections_c() {
     return settings::path_cross_sections.c_str();
   }

src/simulation.F90

@@ -2,20 +2,9 @@ module simulation
 
   use, intrinsic :: ISO_C_BINDING
 
-#ifdef _OPENMP
-  use omp_lib
-#endif
-
-  use bank_header,     only: source_bank
-  use constants,       only: ZERO
-  use error,           only: fatal_error
   use material_header, only: n_materials, materials
-  use message_passing
   use nuclide_header,  only: micro_xs, n_nuclides
   use photon_header,   only: micro_photon_xs, n_elements
-  use settings
-  use simulation_header
-  use tally_header
   use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer
 
   implicit none
@@ -69,56 +58,4 @@ contains
 
   end subroutine
 
-!===============================================================================
-! ALLOCATE_BANKS allocates memory for the fission and source banks
-!===============================================================================
-
-  subroutine allocate_banks() bind(C)
-
-    integer :: alloc_err  ! allocation error code
-
-    ! Allocate source bank
-    if (allocated(source_bank)) deallocate(source_bank)
-    allocate(source_bank(work), STAT=alloc_err)
-
-    ! Check for allocation errors
-    if (alloc_err /= 0) then
-      call fatal_error("Failed to allocate source bank.")
-    end if
-
-    if (run_mode == MODE_EIGENVALUE) then
-
-#ifdef _OPENMP
-      ! If OpenMP is being used, each thread needs its own private fission
-      ! bank. Since the private fission banks need to be combined at the end of
-      ! a generation, there is also a 'master_fission_bank' that is used to
-      ! collect the sites from each thread.
-
-      n_threads = omp_get_max_threads()
-
-!$omp parallel
-      thread_id = omp_get_thread_num()
-
-      if (allocated(fission_bank)) deallocate(fission_bank)
-      if (thread_id == 0) then
-        allocate(fission_bank(3*work))
-      else
-        allocate(fission_bank(3*work/n_threads))
-      end if
-!$omp end parallel
-      if (allocated(master_fission_bank)) deallocate(master_fission_bank)
-      allocate(master_fission_bank(3*work), STAT=alloc_err)
-#else
-      if (allocated(fission_bank)) deallocate(fission_bank)
-      allocate(fission_bank(3*work), STAT=alloc_err)
-#endif
-
-      ! Check for allocation errors
-      if (alloc_err /= 0) then
-        call fatal_error("Failed to allocate fission bank.")
-      end if
-    end if
-
-  end subroutine allocate_banks
-
 end module simulation

src/simulation.cpp

@@ -1,10 +1,12 @@
 #include "openmc/simulation.h"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/container_util.h"
 #include "openmc/eigenvalue.h"
 #include "openmc/error.h"
 #include "openmc/message_passing.h"
+#include "openmc/nuclide.h"
 #include "openmc/output.h"
 #include "openmc/particle.h"
 #include "openmc/random_lcg.h"
@@ -15,6 +17,8 @@
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
 
+#include <omp.h>
+
 #include <algorithm>
 #include <string>
 
@@ -24,14 +28,12 @@ namespace openmc {
 extern "C" bool cmfd_on;
 
 extern "C" void accumulate_tallies();
-extern "C" void allocate_banks();
 extern "C" void allocate_tally_results();
 extern "C" void check_triggers();
 extern "C" void cmfd_init_batch();
 extern "C" void cmfd_tally_init();
 extern "C" void execute_cmfd();
 extern "C" void init_tally_routines();
-extern "C" void join_bank_from_threads();
 extern "C" void load_state_point();
 extern "C" void print_batch_keff();
 extern "C" void print_columns();
@@ -99,6 +101,7 @@ int openmc_simulation_init()
 
   // Call Fortran initialization
   simulation_init_f();
+  set_micro_xs();
 
   // Reset global variables -- this is done before loading state point (as that
   // will potentially populate k_generation and entropy)
@@ -280,6 +283,36 @@ int thread_id;  //!< ID of a given thread
 // Non-member functions
 //==============================================================================
 
+void allocate_banks()
+{
+  // Allocate source bank
+  simulation::source_bank.resize(simulation::work);
+
+  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+#ifdef _OPENMP
+    // If OpenMP is being used, each thread needs its own private fission
+    // bank. Since the private fission banks need to be combined at the end of
+    // a generation, there is also a 'master_fission_bank' that is used to
+    // collect the sites from each thread.
+
+    simulation::n_threads = omp_get_max_threads();
+
+#pragma omp parallel
+    {
+      simulation::thread_id = omp_get_thread_num();
+      if (simulation::thread_id == 0) {
+        simulation::fission_bank.resize(3*simulation::work);
+      } else {
+        simulation::fission_bank.resize(3*simulation::work / simulation::n_threads);
+      }
+    }
+    simulation::master_fission_bank.resize(3*simulation::work);
+#else
+    simulation::fission_bank.resize(3*simulation::work);
+#endif
+  }
+}
+
 void initialize_batch()
 {
   // Increment current batch
@@ -453,13 +486,8 @@ void finalize_generation()
 
 void initialize_history(Particle* p, int64_t index_source)
 {
-  // Get pointer to source bank
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-
   // set defaults
-  p->from_source(&source_bank[index_source - 1]);
+  p->from_source(&simulation::source_bank[index_source - 1]);
 
   // set identifier for particle
   p->id = simulation::work_index[mpi::rank] + index_source;

src/source.cpp

@@ -5,6 +5,7 @@
 
 #include "xtensor/xadapt.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/cell.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
@@ -242,11 +243,6 @@ void initialize_source()
 {
   write_message("Initializing source particles...", 5);
 
-  // Get pointer to source bank
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-
   if (settings::path_source != "") {
     // Read the source from a binary file instead of sampling from some
     // assumed source distribution
@@ -268,7 +264,7 @@ void initialize_source()
     }
 
     // Read in the source bank
-    read_source_bank(file_id, source_bank);
+    read_source_bank(file_id);
 
     // Close file
     file_close(file_id);
@@ -282,7 +278,7 @@ void initialize_source()
       set_particle_seed(id);
 
       // sample external source distribution
-      source_bank[i] = sample_external_source();
+      simulation::source_bank[i] = sample_external_source();
     }
   }
 
@@ -291,7 +287,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, source_bank);
+    write_source_bank(file_id);
     file_close(file_id);
   }
 }
@@ -354,11 +350,6 @@ extern "C" double total_source_strength()
 void fill_source_bank_fixedsource()
 {
   if (settings::path_source.empty()) {
-    // Get pointer to source bank
-    Bank* source_bank;
-    int64_t n;
-    openmc_source_bank(&source_bank, &n);
-
     for (int64_t i = 0; i < simulation::work; ++i) {
       // initialize random number seed
       int64_t id = (simulation::total_gen + overall_generation()) *
@@ -366,7 +357,7 @@ void fill_source_bank_fixedsource()
       set_particle_seed(id);
 
       // sample external source distribution
-      source_bank[i] = sample_external_source();
+      simulation::source_bank[i] = sample_external_source();
     }
   }
 }

src/state_point.F90

@@ -35,16 +35,14 @@ module state_point
   implicit none
 
   interface
-    subroutine write_source_bank(group_id, bank_) bind(C)
-      import HID_T, C_INT64_T, Bank
+    subroutine write_source_bank(group_id) bind(C)
+      import HID_T
       integer(HID_T), value :: group_id
-      type(Bank), intent(in) :: bank_(*)
     end subroutine write_source_bank
 
-    subroutine read_source_bank(group_id, bank_) bind(C)
-      import HID_T, C_INT64_T, Bank
+    subroutine read_source_bank(group_id) bind(C)
+      import HID_T
       integer(HID_T), value :: group_id
-      type(Bank), intent(out) :: bank_(*)
     end subroutine read_source_bank
   end interface
 
@@ -437,7 +435,7 @@ contains
       if (master .or. parallel) then
         file_id = file_open(filename_, 'a', parallel=.true.)
       end if
-      call write_source_bank(file_id, source_bank)
+      call write_source_bank(file_id)
       if (master .or. parallel) call file_close(file_id)
     end if
   end function openmc_statepoint_write
@@ -634,8 +632,8 @@ contains
         file_id = file_open(path_source_point, 'r', parallel=.true.)
       end if
 
-      ! Write out source
-      call read_source_bank(file_id, source_bank)
+      ! Read source
+      call read_source_bank(file_id)
 
     end if
 

src/state_point.cpp

@@ -8,6 +8,7 @@
 #include "xtensor/xbuilder.hpp" // for empty_like
 #include "xtensor/xview.hpp"
 
+#include "openmc/bank.h"
 #include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
@@ -70,16 +71,13 @@ write_source_point(const char* filename)
   }
 
   // Get pointer to source bank and write to file
-  Bank* source_bank;
-  int64_t n;
-  openmc_source_bank(&source_bank, &n);
-  write_source_bank(file_id, source_bank);
+  write_source_bank(file_id);
 
   if (mpi::master || parallel) file_close(file_id);
 }
 
 void
-write_source_bank(hid_t group_id, Bank* source_bank)
+write_source_bank(hid_t group_id)
 {
   hid_t banktype = h5banktype();
 
@@ -103,7 +101,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
 
   // Write data to file in parallel
-  H5Dwrite(dset, banktype, memspace, dspace, plist, source_bank);
+  H5Dwrite(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
 
   // Free resources
   H5Sclose(dspace);
@@ -122,7 +120,8 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 
     // Save source bank sites since the souce_bank array is overwritten below
 #ifdef OPENMC_MPI
-    std::vector<Bank> temp_source {source_bank, source_bank + simulation::work};
+    std::vector<Bank> temp_source {simulation::source_bank.begin(),
+      simulation::source_bank.begin() + simulation::work};
 #endif
 
     for (int i = 0; i < mpi::n_procs; ++i) {
@@ -134,7 +133,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 #ifdef OPENMC_MPI
       // Receive source sites from other processes
       if (i > 0)
-        MPI_Recv(source_bank, count[0], mpi::bank, i, i,
+        MPI_Recv(simulation::source_bank.data(), count[0], mpi::bank, i, i,
                  mpi::intracomm, MPI_STATUS_IGNORE);
 #endif
 
@@ -144,7 +143,8 @@ write_source_bank(hid_t group_id, Bank* source_bank)
       H5Sselect_hyperslab(dspace, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
       // Write data to hyperslab
-      H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
+      H5Dwrite(dset, banktype, memspace, dspace, H5P_DEFAULT,
+        simulation::source_bank.data());
 
       H5Sclose(memspace);
       H5Sclose(dspace);
@@ -155,12 +155,12 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 
 #ifdef OPENMC_MPI
     // Restore state of source bank
-    std::copy(temp_source.begin(), temp_source.end(), source_bank);
+    std::copy(temp_source.begin(), temp_source.end(), simulation::source_bank.begin());
 #endif
   } else {
 #ifdef OPENMC_MPI
-    MPI_Send(source_bank, simulation::work, mpi::bank, 0, mpi::rank,
-             mpi::intracomm);
+    MPI_Send(simulation::source_bank.data(), simulation::work, mpi::bank,
+      0, mpi::rank, mpi::intracomm);
 #endif
   }
 #endif
@@ -169,7 +169,7 @@ write_source_bank(hid_t group_id, Bank* source_bank)
 }
 
 
-void read_source_bank(hid_t group_id, Bank* source_bank)
+void read_source_bank(hid_t group_id)
 {
   hid_t banktype = h5banktype();
 
@@ -197,10 +197,10 @@ void read_source_bank(hid_t group_id, Bank* source_bank)
     // Read data in parallel
   hid_t plist = H5Pcreate(H5P_DATASET_XFER);
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
-  H5Dread(dset, banktype, memspace, dspace, plist, source_bank);
+  H5Dread(dset, banktype, memspace, dspace, plist, simulation::source_bank.data());
   H5Pclose(plist);
 #else
-  H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, source_bank);
+  H5Dread(dset, banktype, memspace, dspace, H5P_DEFAULT, simulation::source_bank.data());
 #endif
 
   // Close all ids

src/tallies/tally.F90

@@ -3,6 +3,7 @@ module tally
   use, intrinsic :: ISO_C_BINDING
 
   use algorithm,        only: binary_search
+  use bank_header
   use constants
   use dict_header,      only: EMPTY
   use error,            only: fatal_error
@@ -687,7 +688,7 @@ contains
             do k = 1, p % n_bank
 
               ! get the delayed group
-              g = fission_bank(n_bank - p % n_bank + k) % delayed_group
+              g = fission_bank_delayed_group(n_bank - p % n_bank + k)
 
               ! Case for tallying delayed emissions
               if (g /= 0) then
@@ -697,8 +698,8 @@ contains
                      reactions(nuclides(p % event_nuclide) % index_fission(1)))
 
                   ! determine score based on bank site weight and keff.
-                  score = score + keff * fission_bank(n_bank - p % n_bank + k) &
-                       % wgt * rxn % product_decay_rate(1 + g) * flux
+                  score = score + keff * fission_bank_wgt(n_bank - p % n_bank + k) &
+                       * rxn % product_decay_rate(1 + g) * flux
                 end associate
 
                 ! if the delayed group filter is present, tally to corresponding
@@ -1843,7 +1844,7 @@ contains
             do k = 1, p % n_bank
 
               ! get the delayed group
-              g = fission_bank(n_bank - p % n_bank + k) % delayed_group
+              g = fission_bank_delayed_group(n_bank - p % n_bank + k)
 
               ! Case for tallying delayed emissions
               if (g /= 0) then
@@ -1851,13 +1852,13 @@ contains
                 ! determine score based on bank site weight and keff.
                 if (i_nuclide > 0) then
                   score = score + keff * atom_density * &
-                       fission_bank(n_bank - p % n_bank + k) % wgt * &
+                       fission_bank_wgt(n_bank - p % n_bank + k) * &
                        get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
                        get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
                        get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
                 else
                   score = score + keff * &
-                       fission_bank(n_bank - p % n_bank + k) % wgt * &
+                       fission_bank_wgt(n_bank - p % n_bank + k) * &
                        get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
                 end if
 
@@ -2375,10 +2376,10 @@ contains
       do k = 1, p % n_bank
 
         ! get the delayed group
-        g = fission_bank(n_bank - p % n_bank + k) % delayed_group
+        g = fission_bank_delayed_group(n_bank - p % n_bank + k)
 
         ! determine score based on bank site weight and keff
-        score = keff * fission_bank(n_bank - p % n_bank + k) % wgt
+        score = keff * fission_bank_wgt(n_bank - p % n_bank + k)
 
         ! Add derivative information for differential tallies.  Note that the
         ! i_nuclide and atom_density arguments do not matter since this is an
@@ -2390,7 +2391,7 @@ contains
         if (.not. run_CE .and. eo_filt % matches_transport_groups) then
 
           ! determine outgoing energy group from fission bank
-          g_out = int(fission_bank(n_bank - p % n_bank + k) % E)
+          g_out = int(fission_bank_E(n_bank - p % n_bank + k))
 
           ! modify the value so that g_out = 1 corresponds to the highest
           ! energy bin
@@ -2403,10 +2404,9 @@ contains
 
           ! determine outgoing energy from fission bank
           if (run_CE) then
-            E_out = fission_bank(n_bank - p % n_bank + k) % E
+            E_out = fission_bank_E(n_bank - p % n_bank + k)
           else
-            E_out = energy_bin_avg(int(fission_bank(n_bank - p % n_bank + k) &
-                 % E))
+            E_out = energy_bin_avg(int(fission_bank_E(n_bank - p % n_bank + k)))
           end if
 
           ! If this outgoing energy falls within the energyout filter's range,

tests/regression_tests/dagmc/test.py

@@ -6,7 +6,7 @@ import pytest
 from tests.testing_harness import PyAPITestHarness
 
 pytestmark = pytest.mark.skipif(
-    not openmc.capi.dagmc_enabled,
+    not openmc.capi._dagmc_enabled(),
     reason="DAGMC CAD geometry is not enabled.")
 
 def test_dagmc():
@@ -22,7 +22,7 @@ def test_dagmc():
     model.settings.source = source
 
     model.settings.dagmc = True
-    
+
     # tally
     tally = openmc.Tally()
     tally.scores = ['total']