Merge pull request #1162 from smharper/cpp_tallies

Move tally scoring to C++

CMakeLists.txt

@@ -349,16 +349,10 @@ add_library(libopenmc SHARED
   src/tallies/tally_derivative_header.F90
   src/tallies/tally_filter.F90
   src/tallies/tally_filter_header.F90
-  src/tallies/tally_filter_delayedgroup.F90
   src/tallies/tally_filter_distribcell.F90
-  src/tallies/tally_filter_energy.F90
-  src/tallies/tally_filter_mesh.F90
-  src/tallies/tally_filter_meshsurface.F90
   src/tallies/tally_filter_particle.F90
   src/tallies/tally_filter_sph_harm.F90
   src/tallies/tally_header.F90
-  src/tallies/trigger.F90
-  src/tallies/trigger_header.F90
   src/bank.cpp
   src/bremsstrahlung.cpp
   src/dagmc.cpp
@@ -412,6 +406,7 @@ add_library(libopenmc SHARED
   src/string_utils.cpp
   src/summary.cpp
   src/surface.cpp
+  src/tallies/derivative.cpp
   src/tallies/filter.cpp
   src/tallies/filter_azimuthal.cpp
   src/tallies/filter_cellborn.cpp
@@ -434,6 +429,8 @@ add_library(libopenmc SHARED
   src/tallies/filter_universe.cpp
   src/tallies/filter_zernike.cpp
   src/tallies/tally.cpp
+  src/tallies/tally_scoring.cpp
+  src/tallies/trigger.cpp
   src/timer.cpp
   src/thermal.cpp
   src/wmp.cpp

include/openmc/capi.h

@@ -105,7 +105,7 @@ extern "C" {
   int openmc_tally_get_active(int32_t index, bool* active);
   int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
   int openmc_tally_get_id(int32_t index, int32_t* id);
-  int openmc_tally_get_filters(int32_t index, int32_t** indices, int* n);
+  int openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n);
   int openmc_tally_get_n_realizations(int32_t index, int32_t* n);
   int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n);
   int openmc_tally_get_scores(int32_t index, int** scores, int* n);

include/openmc/constants.h

@@ -368,24 +368,6 @@ constexpr int NO_BIN_FOUND {-1};
 constexpr int FILTER_UNIVERSE       {1};
 constexpr int FILTER_MATERIAL       {2};
 constexpr int FILTER_CELL           {3};
-constexpr int FILTER_CELLBORN       {4};
-constexpr int FILTER_SURFACE        {5};
-constexpr int FILTER_MESH           {6};
-constexpr int FILTER_ENERGYIN       {7};
-constexpr int FILTER_ENERGYOUT      {8};
-constexpr int FILTER_DISTRIBCELL    {9};
-constexpr int FILTER_MU             {10};
-constexpr int FILTER_POLAR          {11};
-constexpr int FILTER_AZIMUTHAL      {12};
-constexpr int FILTER_DELAYEDGROUP   {13};
-constexpr int FILTER_ENERGYFUNCTION {14};
-constexpr int FILTER_CELLFROM       {15};
-constexpr int FILTER_MESHSURFACE    {16};
-constexpr int FILTER_LEGENDRE       {17};
-constexpr int FILTER_SPH_HARMONICS  {18};
-constexpr int FILTER_SPTL_LEGENDRE  {19};
-constexpr int FILTER_ZERNIKE        {20};
-constexpr int FILTER_PARTICLE       {21};
 
 // Mesh types
 constexpr int MESH_REGULAR {1};
@@ -404,11 +386,6 @@ constexpr int IN_BOTTOM  {10}; // z min
 constexpr int OUT_TOP    {11}; // z max
 constexpr int IN_TOP     {12}; // z max
 
-// Tally trigger types and threshold
-constexpr int VARIANCE           {1};
-constexpr int RELATIVE_ERROR     {2};
-constexpr int STANDARD_DEVIATION {3};
-
 // Global tally parameters
 constexpr int N_GLOBAL_TALLIES {4};
 constexpr int K_COLLISION   {0};
@@ -416,11 +393,6 @@ constexpr int K_ABSORPTION  {1};
 constexpr int K_TRACKLENGTH {2};
 constexpr int LEAKAGE       {3};
 
-// Differential tally independent variables
-constexpr int DIFF_DENSITY {1};
-constexpr int DIFF_NUCLIDE_DENSITY {2};
-constexpr int DIFF_TEMPERATURE {3};
-
 // Miscellaneous
 constexpr int C_NONE {-1};
 constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE

include/openmc/mgxs.h

@@ -134,7 +134,8 @@ class Mgxs {
     //! @param dg delayed group index; use nullptr if irrelevant.
     //! @return Requested cross section value.
     double
-    get_xs(int xstype, int gin, int* gout, double* mu, int* dg);
+    get_xs(int xstype, int gin, const int* gout, const double* mu,
+      const int* dg);
 
     //! \brief Samples the fission neutron energy and if prompt or delayed.
     //!

include/openmc/mgxs_interface.h

@@ -51,20 +51,21 @@ extern "C" void
 calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
      double& total_xs, double& abs_xs, double& nu_fiss_xs);
 
-extern "C" double
-get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
+double
+get_nuclide_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg);
 
-extern "C" double
-get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
+inline double
+get_nuclide_xs(int index, int xstype, int gin)
+{return get_nuclide_xs(index, xstype, gin, nullptr, nullptr, nullptr);}
 
-extern "C" void
-set_nuclide_angle_index_c(int index, const double uvw[3]);
+double
+get_macro_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg);
 
-extern "C" void
-set_macro_angle_index_c(int index, const double uvw[3]);
-
-extern "C" void
-set_nuclide_temperature_index_c(int index, double sqrtkT);
+inline double
+get_macro_xs(int index, int xstype, int gin)
+{return get_macro_xs(index, xstype, gin, nullptr, nullptr, nullptr);}
 
 //==============================================================================
 // General Mgxs methods

include/openmc/nuclide.h

@@ -26,10 +26,10 @@ namespace openmc {
 
 constexpr double CACHE_INVALID {-1.0};
 
-//===============================================================================
+//==============================================================================
 //! Cached microscopic cross sections for a particular nuclide at the current
 //! energy
-//===============================================================================
+//==============================================================================
 
 struct NuclideMicroXS {
   // Microscopic cross sections in barns
@@ -64,10 +64,10 @@ struct NuclideMicroXS {
                             //!< * temperature (eV))
 };
 
-//===============================================================================
+//==============================================================================
 // MATERIALMACROXS contains cached macroscopic cross sections for the material a
 // particle is traveling through
-//===============================================================================
+//==============================================================================
 
 struct MaterialMacroXS {
   double total;         //!< macroscopic total xs
@@ -83,9 +83,9 @@ struct MaterialMacroXS {
   double pair_production; //!< macroscopic pair production xs
 };
 
-//===============================================================================
+//==============================================================================
 // Data for a nuclide
-//===============================================================================
+//==============================================================================
 
 class Nuclide {
 public:
@@ -176,6 +176,8 @@ private:
 //! Checks for the right version of nuclear data within HDF5 files
 void check_data_version(hid_t file_id);
 
+extern "C" bool multipole_in_range(const Nuclide* nuc, double E);
+
 //==============================================================================
 // Global variables
 //==============================================================================

include/openmc/tallies/derivative.h

@@ -0,0 +1,64 @@
+#ifndef OPENMC_TALLIES_DERIVATIVE_H
+#define OPENMC_TALLIES_DERIVATIVE_H
+
+#include "openmc/particle.h"
+
+#include <unordered_map>
+#include <vector>
+
+#include "pugixml.hpp"
+
+//==============================================================================
+//! Describes a first-order derivative that can be applied to tallies.
+//==============================================================================
+
+namespace openmc {
+
+struct TallyDerivative {
+  int id;  //!< User-defined identifier
+  int variable;  //!< Independent variable (like temperature)
+  int diff_material;  //!< Material this derivative is applied to
+  int diff_nuclide;  //!< Nuclide this material is applied to
+  double flux_deriv;  //!< Derivative of the current particle's weight
+
+  TallyDerivative() {}
+  explicit TallyDerivative(pugi::xml_node node);
+};
+
+//==============================================================================
+// Non-method functions
+//==============================================================================
+
+//! Scale the given score by its logarithmic derivative
+
+void
+apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
+  double atom_density, int score_bin, double& score);
+
+} // namespace openmc
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+// Explicit vector template specialization declaration of threadprivate variable
+// outside of the openmc namespace for the picky Intel compiler.
+extern template class std::vector<openmc::TallyDerivative>;
+
+namespace openmc {
+
+namespace model {
+extern std::vector<TallyDerivative> tally_derivs;
+#pragma omp threadprivate(tally_derivs)
+extern std::unordered_map<int, int> tally_deriv_map;
+} // namespace model
+
+// Independent variables
+//TODO: convert to enum
+constexpr int DIFF_DENSITY {1};
+constexpr int DIFF_NUCLIDE_DENSITY {2};
+constexpr int DIFF_TEMPERATURE {3};
+
+} // namespace openmc
+
+#endif // OPENMC_TALLIES_DERIVATIVE_H

include/openmc/tallies/filter.h

@@ -22,6 +22,8 @@ class FilterMatch
 public:
   std::vector<int> bins_;
   std::vector<double> weights_;
+  int i_bin_;
+  bool bins_present_ {false};
 };
 
 } // namespace openmc
@@ -70,6 +72,8 @@ public:
 
   virtual void initialize() {}
 
+  int32_t id_;
+
   int n_bins_;
 };
 
@@ -93,10 +97,6 @@ extern std::vector<std::unique_ptr<Filter>> tally_filters;
 
 //==============================================================================
 
-extern "C" void free_memory_tally_c();
-
-//==============================================================================
-
 // Filter-related Fortran functions that will be called from C++
 extern "C" int verify_filter(int32_t index);
 extern "C" Filter* filter_from_f(int32_t index);

include/openmc/tallies/tally.h

@@ -2,17 +2,120 @@
 #define OPENMC_TALLIES_TALLY_H
 
 #include "openmc/constants.h"
+#include "openmc/tallies/trigger.h"
 
+#include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
 
+#include <memory> // for unique_ptr
+#include <string>
+#include <vector>
+
 namespace openmc {
 
+//==============================================================================
+//! A user-specified flux-weighted (or current) measurement.
+//==============================================================================
+
+class Tally {
+public:
+  Tally() {}
+
+  void init_from_xml(pugi::xml_node node);
+
+  void set_scores(pugi::xml_node node);
+
+  void set_scores(std::vector<std::string> scores);
+
+  void set_nuclides(pugi::xml_node node);
+
+  //----------------------------------------------------------------------------
+  // Methods for getting and setting filter/stride data.
+
+  const std::vector<int32_t>& filters() const {return filters_;}
+
+  int32_t filters(int i) const {return filters_[i];}
+
+  void set_filters(const int32_t filter_indices[], int n);
+
+  int32_t strides(int i) const {return strides_[i];}
+
+  int32_t n_filter_bins() const {return n_filter_bins_;}
+
+  //----------------------------------------------------------------------------
+  // Other methods.
+
+  void init_triggers(pugi::xml_node node, int i_tally);
+
+  //----------------------------------------------------------------------------
+  // Major public data members.
+
+  int id_; //!< User-defined identifier
+
+  std::string name_; //!< User-defined name
+
+  int type_ {TALLY_VOLUME}; //!< e.g. volume, surface current
+
+  //! Event type that contributes to this tally
+  int estimator_ {ESTIMATOR_TRACKLENGTH};
+
+  //! Whether this tally is currently being updated
+  bool active_ {false};
+
+  std::vector<int> scores_; //!< Filter integrands (e.g. flux, fission)
+
+  //! Index of each nuclide to be tallied.  -1 indicates total material.
+  std::vector<int> nuclides_ {-1};
+
+  //! True if this tally has a bin for every nuclide in the problem
+  bool all_nuclides_ {false};
+
+  //----------------------------------------------------------------------------
+  // Miscellaneous public members.
+
+  // We need to have quick access to some filters.  The following gives indices
+  // for various filters that could be in the tally or C_NONE if they are not
+  // present.
+  int energyout_filter_ {C_NONE};
+  int delayedgroup_filter_ {C_NONE};
+
+  bool depletion_rx_ {false}; //!< Has depletion reactions (e.g. (n,2n))
+
+  std::vector<Trigger> triggers_;
+
+  int deriv_ {C_NONE}; //!< Index of a TallyDerivative object for diff tallies.
+
+private:
+  //----------------------------------------------------------------------------
+  // Private data.
+
+  std::vector<int32_t> filters_; //!< Filter indices in global filters array
+
+  //! Index strides assigned to each filter to support 1D indexing.
+  std::vector<int32_t> strides_;
+
+  int32_t n_filter_bins_ {0};
+};
+
 //==============================================================================
 // Global variable declarations
 //==============================================================================
 
 extern "C" double total_weight;
 
+namespace model {
+
+extern std::vector<std::unique_ptr<Tally>> tallies;
+
+extern std::vector<int> active_tallies;
+extern std::vector<int> active_analog_tallies;
+extern std::vector<int> active_tracklength_tallies;
+extern std::vector<int> active_collision_tallies;
+extern std::vector<int> active_meshsurf_tallies;
+extern std::vector<int> active_surface_tallies;
+
+} // namespace model
+
 // Threadprivate variables
 extern "C" double global_tally_absorption;
 extern "C" double global_tally_collision;
@@ -44,6 +147,8 @@ adaptor_type<3> tally_results(int idx);
 extern "C" void reduce_tally_results();
 #endif
 
+extern "C" void free_memory_tally_c();
+
 } // namespace openmc
 
 #endif // OPENMC_TALLIES_TALLY_H

include/openmc/tallies/tally_scoring.h

@@ -0,0 +1,55 @@
+#ifndef OPENMC_TALLIES_TALLY_SCORING_H
+#define OPENMC_TALLIES_TALLY_SCORING_H
+
+#include "openmc/particle.h"
+#include "openmc/tallies/tally.h"
+
+namespace openmc {
+
+//==============================================================================
+//! An iterator over all combinations of a tally's matching filter bins.
+//
+//! This iterator handles two distinct tasks.  First, it maps the N-dimensional
+//! space created by the indices of N filters onto a 1D sequence.  In other
+//! words, it provides a single number that uniquely identifies a combination of
+//! bins for many filters.  Second, it handles the task of finding each valid
+//! combination of filter bins given that each filter can have 1 or 2 or many
+//! bins that are valid for the current tally event.
+//==============================================================================
+
+class FilterBinIter
+{
+public:
+
+  //! Construct an iterator over bins that match a given particle's state.
+  FilterBinIter(const Tally& tally, const Particle* p);
+
+  //! Construct an iterator over all filter bin combinations.
+  //
+  //! \param end if true, the returned iterator indicates the end of a loop.
+  FilterBinIter(const Tally& tally, bool end);
+
+  bool operator==(const FilterBinIter& other) const
+  {return index_ == other.index_;}
+
+  bool operator!=(const FilterBinIter& other) const
+  {return !(*this == other);}
+
+  FilterBinIter& operator++();
+
+  int index_ {1};
+  double weight_ {1.};
+
+private:
+  void compute_index_weight();
+
+  const Tally& tally_;
+};
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+} // namespace openmc
+
+#endif // OPENMC_TALLIES_TALLY_SCORING_H

include/openmc/tallies/trigger.h

@@ -0,0 +1,50 @@
+#ifndef OPENMC_TALLIES_TRIGGER_H
+#define OPENMC_TALLIES_TRIGGER_H
+
+#include <string>
+
+#include "pugixml.hpp"
+
+namespace openmc {
+
+//==============================================================================
+// Type definitions
+//==============================================================================
+
+enum class TriggerMetric {
+  variance, relative_error, standard_deviation, not_active
+};
+
+//! Stops the simulation early if a desired tally uncertainty is reached.
+
+struct Trigger {
+  TriggerMetric metric;  //!< The type of uncertainty (e.g. std dev) measured
+  double threshold;  //!< Uncertainty value below which trigger is satisfied
+  int score_index;  //!< Index of the relevant score in the tally's arrays
+};
+
+//! Stops the simulation early if a desired k-effective uncertainty is reached.
+
+struct KTrigger
+{
+  TriggerMetric metric {TriggerMetric::not_active};
+  double threshold {0.};
+};
+
+//==============================================================================
+// Global variable declarations
+//==============================================================================
+
+//TODO: consider a different namespace
+namespace settings {
+  extern KTrigger keff_trigger;
+}
+
+//==============================================================================
+// Non-memeber functions
+//==============================================================================
+
+void check_triggers();
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_TRIGGER_H

openmc/capi/tally.py

@@ -254,13 +254,13 @@ class Tally(_FortranObjectWithID):
         filt_idx = POINTER(c_int32)()
         n = c_int()
         _dll.openmc_tally_get_filters(self._index, filt_idx, n)
-        return [_get_filter(filt_idx[i]) for i in range(n.value)]
+        return [_get_filter(filt_idx[i]+1) for i in range(n.value)]
 
     @filters.setter
     def filters(self, filters):
         # Get filter indices as int32_t[]
         n = len(filters)
-        indices = (c_int32*n)(*(f._index for f in filters))
+        indices = (c_int32*n)(*(f._index-1 for f in filters))
 
         _dll.openmc_tally_set_filters(self._index, n, indices)
 
@@ -278,7 +278,7 @@ class Tally(_FortranObjectWithID):
         nucs = POINTER(c_int)()
         n = c_int()
         _dll.openmc_tally_get_nuclides(self._index, nucs, n)
-        return [Nuclide(nucs[i]).name if nucs[i] > 0 else 'total'
+        return [Nuclide(nucs[i]+1).name if nucs[i] >= 0 else 'total'
                 for i in range(n.value)]
 
     @nuclides.setter

src/api.F90

@@ -65,10 +65,8 @@ contains
     use settings
     use simulation_header
     use surface_header
-    use tally_derivative_header
     use tally_filter_header
     use tally_header
-    use trigger_header
     use volume_header
 
     interface
@@ -107,7 +105,6 @@ contains
     call free_memory_mesh()
     call free_memory_tally()
     call free_memory_tally_filter()
-    call free_memory_tally_derivative()
     call free_memory_bank()
     call free_memory_cmfd()
 #ifdef DAGMC

src/constants.F90

@@ -201,30 +201,10 @@ module constants
   integer, parameter :: NO_BIN_FOUND = -1
 
   ! Tally filter and map types
-  integer, parameter :: N_FILTER_TYPES = 22
   integer, parameter :: &
        FILTER_UNIVERSE       = 1,  &
        FILTER_MATERIAL       = 2,  &
-       FILTER_CELL           = 3,  &
-       FILTER_CELLBORN       = 4,  &
-       FILTER_SURFACE        = 5,  &
-       FILTER_MESH           = 6,  &
-       FILTER_ENERGYIN       = 7,  &
-       FILTER_ENERGYOUT      = 8,  &
-       FILTER_DISTRIBCELL    = 9,  &
-       FILTER_MU             = 10, &
-       FILTER_POLAR          = 11, &
-       FILTER_AZIMUTHAL      = 12, &
-       FILTER_DELAYEDGROUP   = 13, &
-       FILTER_ENERGYFUNCTION = 14, &
-       FILTER_CELLFROM       = 15, &
-       FILTER_MESHSURFACE    = 16, &
-       FILTER_LEGENDRE       = 17, &
-       FILTER_SPH_HARMONICS  = 18, &
-       FILTER_SPTL_LEGENDRE  = 19, &
-       FILTER_ZERNIKE        = 20, &
-       FILTER_ZERNIKE_RADIAL = 21, &
-       FILTER_PARTICLE       = 22
+       FILTER_CELL           = 3
 
   ! Tally surface current directions
   integer, parameter :: &
@@ -241,12 +221,6 @@ module constants
        OUT_TOP    = 11,  &   ! z max
        IN_TOP     = 12       ! z max
 
-  ! Tally trigger types and threshold
-  integer, parameter :: &
-       VARIANCE           = 1, &
-       RELATIVE_ERROR     = 2, &
-       STANDARD_DEVIATION = 3
-
   ! Global tally parameters
   integer, parameter :: N_GLOBAL_TALLIES = 4
   integer, parameter :: &

src/input_xml.F90

@@ -29,7 +29,6 @@ module input_xml
   use tally_derivative_header
   use tally_filter_header
   use tally_filter
-  use trigger_header
   use volume_header
   use xml_interface
 
@@ -446,53 +445,67 @@ contains
 
     integer :: i             ! loop over user-specified tallies
     integer :: j             ! loop over words
-    integer :: k             ! another loop index
     integer :: l             ! loop over bins
     integer :: filter_id     ! user-specified identifier for filter
     integer :: tally_id      ! user-specified identifier for filter
+    integer :: deriv_id
     integer :: i_filt        ! index in filters array
-    integer :: i_elem        ! index of entry in dictionary
     integer :: n             ! size of arrays in mesh specification
     integer :: n_words       ! number of words read
     integer :: n_filter      ! number of filters
-    integer :: n_scores      ! number of scores
-    integer :: n_user_trig   ! number of user-specified tally triggers
-    integer :: trig_ind      ! index of triggers array for each tally
-    integer :: user_trig_ind ! index of user-specified triggers for each tally
     integer :: i_start, i_end
     integer(C_INT) :: err
-    real(8) :: threshold     ! trigger convergence threshold
-    integer :: MT            ! user-specified MT for score
     logical :: file_exists   ! does tallies.xml file exist?
     integer, allocatable :: temp_filter(:) ! temporary filter indices
+    logical :: has_energyout
+    integer :: particle_filter_index
     character(MAX_LINE_LEN) :: filename
-    character(MAX_WORD_LEN) :: word
-    character(MAX_WORD_LEN) :: score_name
     character(MAX_WORD_LEN) :: temp_str
-    character(MAX_WORD_LEN), allocatable :: sarray(:)
-    type(DictCharInt) :: trigger_scores
     type(TallyFilterContainer), pointer :: f
     type(XMLDocument) :: doc
     type(XMLNode) :: root
     type(XMLNode) :: node_tal
     type(XMLNode) :: node_filt
-    type(XMLNode) :: node_trigger
     type(XMLNode), allocatable :: node_tal_list(:)
     type(XMLNode), allocatable :: node_filt_list(:)
-    type(XMLNode), allocatable :: node_trigger_list(:)
-    type(XMLNode), allocatable :: node_deriv_list(:)
-    type(DictEntryCI) :: elem
+    type(TallyDerivative), pointer :: deriv
+
+    interface
+      subroutine tally_init_from_xml(tally_ptr, xml_node) bind(C)
+        import C_PTR
+        type(C_PTR), value :: tally_ptr
+        type(C_PTR) :: xml_node
+      end subroutine
+
+      subroutine tally_set_scores(tally_ptr, xml_node) bind(C)
+        import C_PTR
+        type(C_PTR), value :: tally_ptr
+        type(C_PTR) :: xml_node
+      end subroutine
+
+      subroutine tally_set_nuclides(tally_ptr, xml_node) bind(C)
+        import C_PTR
+        type(C_PTR), value :: tally_ptr
+        type(C_PTR) :: xml_node
+      end subroutine
+
+      subroutine tally_init_triggers(tally_ptr, i_tally, xml_node) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: tally_ptr
+        integer(C_INT), value :: i_tally
+        type(C_PTR) :: xml_node
+      end subroutine
+
+      subroutine read_tally_derivatives(node_ptr) bind(C)
+        import C_PTR
+        type(C_PTR) :: node_ptr
+      end subroutine
+    end interface
 
     ! Check if tallies.xml exists
     filename = trim(path_input) // "tallies.xml"
     inquire(FILE=filename, EXIST=file_exists)
     if (.not. file_exists) then
-      ! We need to allocate tally_derivs to avoid segfaults.  Also needs to be
-      ! done in parallel because tally derivs are threadprivate.
-!$omp parallel
-      allocate(tally_derivs(0))
-!$omp end parallel
-
       ! Since a tallies.xml file is optional, no error is issued here
       return
     end if
@@ -533,29 +546,7 @@ contains
     ! ==========================================================================
     ! READ DATA FOR DERIVATIVES
 
-    ! Get pointer list to XML <derivative> nodes and allocate global array.
-    ! The array is threadprivate so it must be allocated in parallel.
-    call get_node_list(root, "derivative", node_deriv_list)
-!$omp parallel
-    allocate(tally_derivs(size(node_deriv_list)))
-!$omp end parallel
-
-    ! Make sure this is not an MG run.
-    if (.not. run_CE .and. size(node_deriv_list) > 0) then
-      call fatal_error("Differential tallies not supported in multi-group mode")
-    end if
-
-    ! Read derivative attributes.
-    do i = 1, size(node_deriv_list)
-!$omp parallel
-!$omp critical (ReadTallyDeriv)
-      call tally_derivs(i) % from_xml(node_deriv_list(i))
-!$omp end critical (ReadTallyDeriv)
-!$omp end parallel
-
-      ! Update tally derivative dictionary
-      call tally_deriv_dict % set(tally_derivs(i) % id, i)
-    end do
+    call read_tally_derivatives(root % ptr)
 
     ! ==========================================================================
     ! READ FILTER DATA
@@ -643,6 +634,8 @@ contains
       ! Get pointer to tally xml node
       node_tal = node_tal_list(i)
 
+      call tally_init_from_xml(t % ptr, node_tal % ptr)
+
       ! Copy and set tally id
       if (check_for_node(node_tal, "id")) then
         call get_node_value(node_tal, "id", tally_id)
@@ -688,11 +681,11 @@ contains
           else
             call fatal_error("Could not find filter " &
                  // trim(to_str(temp_filter(j))) // " specified on tally " &
-                 // trim(to_str(t % id)))
+                 // trim(to_str(t % id())))
           end if
 
           ! Store the index of the filter
-          temp_filter(j) = i_filt
+          temp_filter(j) = i_filt - 1
         end do
 
         ! Set the filters
@@ -700,368 +693,53 @@ contains
       end if
       deallocate(temp_filter)
 
+      ! Check for the presence of certain filter types
+      has_energyout = (t % energyout_filter() > 0)
+      particle_filter_index = 0
+      do j = 1, t % n_filters()
+        select type (filt => filters(t % filter(j) + 1) % obj)
+        type is (ParticleFilter)
+          particle_filter_index = j
+        end select
+      end do
+
+      ! Change the tally estimator if a filter demands it
+      do j = 1, t % n_filters()
+        select type (filt => filters(t % filter(j) + 1) % obj)
+        type is (EnergyoutFilter)
+          call t % set_estimator(ESTIMATOR_ANALOG)
+        type is (LegendreFilter)
+          call t % set_estimator(ESTIMATOR_ANALOG)
+        type is (SphericalHarmonicsFilter)
+          if (filt % cosine() == COSINE_SCATTER) then
+            call t % set_estimator(ESTIMATOR_ANALOG)
+          end if
+        type is (SpatialLegendreFilter)
+          call t % set_estimator(ESTIMATOR_COLLISION)
+        type is (ZernikeFilter)
+          call t % set_estimator(ESTIMATOR_COLLISION)
+        type is (ZernikeRadialFilter)
+          call t % set_estimator(ESTIMATOR_COLLISION)
+        end select
+      end do
+
       ! =======================================================================
       ! READ DATA FOR NUCLIDES
 
-      if (check_for_node(node_tal, "nuclides")) then
-
-        ! Allocate a temporary string array for nuclides and copy values over
-        allocate(sarray(node_word_count(node_tal, "nuclides")))
-        call get_node_array(node_tal, "nuclides", sarray)
-
-        if (trim(sarray(1)) == 'all') then
-          ! Handle special case <nuclides>all</nuclides>
-          allocate(t % nuclide_bins(n_nuclides + 1))
-
-          ! Set bins to 1, 2, 3, ..., n_nuclides, -1
-          t % nuclide_bins(1:n_nuclides) = &
-               (/ (j, j=1, n_nuclides) /)
-          t % nuclide_bins(n_nuclides + 1) = -1
-
-          ! Set number of nuclide bins
-          t % n_nuclide_bins = n_nuclides + 1
-
-          ! Set flag so we can treat this case specially
-          t % all_nuclides = .true.
-        else
-          ! Any other case, e.g. <nuclides>U-235 Pu-239</nuclides>
-          n_words = node_word_count(node_tal, "nuclides")
-          allocate(t % nuclide_bins(n_words))
-          do j = 1, n_words
-
-            ! Check if total material was specified
-            if (trim(sarray(j)) == 'total') then
-              t % nuclide_bins(j) = -1
-              cycle
-            end if
-
-            ! If a specific nuclide was specified
-            word = sarray(j)
-
-            ! Search through nuclides
-            k = nuclide_map_get(to_c_string(word))
-            if (k == -1) then
-              call fatal_error("Could not find the nuclide " &
-                   // trim(word) // " specified in tally " &
-                   // trim(to_str(t % id)) // " in any material.")
-            end if
-
-            ! Set bin to index in nuclides array
-            t % nuclide_bins(j) = k
-          end do
-
-          ! Set number of nuclide bins
-          t % n_nuclide_bins = n_words
-        end if
-
-        ! Deallocate temporary string array
-        deallocate(sarray)
-
-      else
-        ! No <nuclides> were specified -- create only one bin will be added
-        ! for the total material.
-        allocate(t % nuclide_bins(1))
-        t % nuclide_bins(1) = -1
-        t % n_nuclide_bins = 1
-      end if
+      call tally_set_nuclides(t % ptr, node_tal % ptr)
 
       ! =======================================================================
       ! READ DATA FOR SCORES
 
+      call tally_set_scores(t % ptr, node_tal % ptr)
+
       if (check_for_node(node_tal, "scores")) then
         n_words = node_word_count(node_tal, "scores")
-        allocate(sarray(n_words))
-        call get_node_array(node_tal, "scores", sarray)
-
-        ! Append the score to the list of possible trigger scores
-        do j = 1, n_words
-          sarray(j) = to_lower(sarray(j))
-          score_name = trim(sarray(j))
-
-          if (trigger_on) call trigger_scores % set(trim(score_name), j)
-
-        end do
-        n_scores = n_words
-
-        ! Allocate score storage accordingly
-        allocate(t % score_bins(n_scores))
-
-        ! Check the validity of the scores and their filters
-        do j = 1, n_scores
-          score_name = sarray(j)
-
-          ! Check if delayed group filter is used with any score besides
-          ! delayed-nu-fission or decay-rate
-          if ((score_name /= 'delayed-nu-fission' .and. &
-               score_name /= 'decay-rate') .and. &
-               t % find_filter(FILTER_DELAYEDGROUP) > 0) then
-            call fatal_error("Cannot tally " // trim(score_name) // " with a &
-                 &delayedgroup filter.")
-          end if
-
-          select case (trim(score_name))
-          case ('flux')
-            ! Prohibit user from tallying flux for an individual nuclide
-            if (.not. (t % n_nuclide_bins == 1 .and. &
-                 t % nuclide_bins(1) == -1)) then
-              call fatal_error("Cannot tally flux for an individual nuclide.")
-            end if
-
-            t % score_bins(j) = SCORE_FLUX
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              call fatal_error("Cannot tally flux with an outgoing energy &
-                   &filter.")
-            end if
-
-          case ('total', '(n,total)')
-            t % score_bins(j) = SCORE_TOTAL
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              call fatal_error("Cannot tally total reaction rate with an &
-                   &outgoing energy filter.")
-            end if
-
-          case ('scatter')
-            t % score_bins(j) = SCORE_SCATTER
-            if (t % find_filter(FILTER_ENERGYOUT) > 0 .or. &
-                 t % find_filter(FILTER_LEGENDRE) > 0) then
-              ! Set tally estimator to analog
-              t % estimator = ESTIMATOR_ANALOG
-            end if
-
-          case ('nu-scatter')
-            t % score_bins(j) = SCORE_NU_SCATTER
-
-            ! Set tally estimator to analog for CE mode
-            ! (MG mode has all data available without a collision being
-            ! necessary)
-            if (run_CE) then
-              t % estimator = ESTIMATOR_ANALOG
-            else
-              if (t % find_filter(FILTER_ENERGYOUT) > 0 .or. &
-                   t % find_filter(FILTER_LEGENDRE) > 0) then
-                ! Set tally estimator to analog
-                t % estimator = ESTIMATOR_ANALOG
-              end if
-            end if
-
-          case ('n2n', '(n,2n)')
-            t % score_bins(j) = N_2N
-            t % depletion_rx = .true.
-
-          case ('n3n', '(n,3n)')
-            t % score_bins(j) = N_3N
-            t % depletion_rx = .true.
-
-          case ('n4n', '(n,4n)')
-            t % score_bins(j) = N_4N
-            t % depletion_rx = .true.
-
-          case ('absorption')
-            t % score_bins(j) = SCORE_ABSORPTION
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              call fatal_error("Cannot tally absorption rate with an outgoing &
-                   &energy filter.")
-            end if
-          case ('fission', '18')
-            t % score_bins(j) = SCORE_FISSION
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              call fatal_error("Cannot tally fission rate with an outgoing &
-                   &energy filter.")
-            end if
-          case ('nu-fission')
-            t % score_bins(j) = SCORE_NU_FISSION
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              ! Set tally estimator to analog
-              t % estimator = ESTIMATOR_ANALOG
-            end if
-          case ('decay-rate')
-            t % score_bins(j) = SCORE_DECAY_RATE
-          case ('delayed-nu-fission')
-            t % score_bins(j) = SCORE_DELAYED_NU_FISSION
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              ! Set tally estimator to analog
-              t % estimator = ESTIMATOR_ANALOG
-            end if
-          case ('prompt-nu-fission')
-            t % score_bins(j) = SCORE_PROMPT_NU_FISSION
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              ! Set tally estimator to analog
-              t % estimator = ESTIMATOR_ANALOG
-            end if
-          case ('kappa-fission')
-            t % score_bins(j) = SCORE_KAPPA_FISSION
-          case ('inverse-velocity')
-            t % score_bins(j) = SCORE_INVERSE_VELOCITY
-          case ('fission-q-prompt')
-            t % score_bins(j) = SCORE_FISS_Q_PROMPT
-          case ('fission-q-recoverable')
-            t % score_bins(j) = SCORE_FISS_Q_RECOV
-          case ('current')
-
-            ! Check which type of current is desired: mesh currents or
-            ! surface currents
-            if (t % find_filter(FILTER_SURFACE) > 0 .or. &
-                 &t % find_filter(FILTER_CELL) > 0 .or. &
-                 &t % find_filter(FILTER_CELLFROM) > 0) then
-
-              ! Check to make sure that mesh surface currents are not desired as well
-              if (t % find_filter(FILTER_MESHSURFACE) > 0) then
-                call fatal_error("Cannot tally mesh surface currents &
-                     &in the same tally as normal surface currents")
-              end if
-
-              t % type = TALLY_SURFACE
-              t % score_bins(j) = SCORE_CURRENT
-
-            else if (t % find_filter(FILTER_MESHSURFACE) > 0) then
-              t % score_bins(j) = SCORE_CURRENT
-              t % type = TALLY_MESH_SURFACE
-
-              ! Check to make sure that current is the only desired response
-              ! for this tally
-              if (n_words > 1) then
-                call fatal_error("Cannot tally other scores in the &
-                     &same tally as surface currents")
-              end if
-            else
-                call fatal_error("Cannot tally currents without surface &
-                     &type filters")
-            end if
-
-          case ('events')
-            t % score_bins(j) = SCORE_EVENTS
-          case ('elastic', '(n,elastic)')
-            t % score_bins(j) = ELASTIC
-          case ('(n,2nd)')
-            t % score_bins(j) = N_2ND
-          case ('(n,na)')
-            t % score_bins(j) = N_2NA
-          case ('(n,n3a)')
-            t % score_bins(j) = N_N3A
-          case ('(n,2na)')
-            t % score_bins(j) = N_2NA
-          case ('(n,3na)')
-            t % score_bins(j) = N_3NA
-          case ('(n,np)')
-            t % score_bins(j) = N_NP
-          case ('(n,n2a)')
-            t % score_bins(j) = N_N2A
-          case ('(n,2n2a)')
-            t % score_bins(j) = N_2N2A
-          case ('(n,nd)')
-            t % score_bins(j) = N_ND
-          case ('(n,nt)')
-            t % score_bins(j) = N_NT
-          case ('(n,nHe-3)')
-            t % score_bins(j) = N_N3HE
-          case ('(n,nd2a)')
-            t % score_bins(j) = N_ND2A
-          case ('(n,nt2a)')
-            t % score_bins(j) = N_NT2A
-          case ('(n,3nf)')
-            t % score_bins(j) = N_3NF
-          case ('(n,2np)')
-            t % score_bins(j) = N_2NP
-          case ('(n,3np)')
-            t % score_bins(j) = N_3NP
-          case ('(n,n2p)')
-            t % score_bins(j) = N_N2P
-          case ('(n,npa)')
-            t % score_bins(j) = N_NPA
-          case ('(n,n1)')
-            t % score_bins(j) = N_N1
-          case ('(n,nc)')
-            t % score_bins(j) = N_NC
-          case ('(n,gamma)')
-            t % score_bins(j) = N_GAMMA
-            t % depletion_rx = .true.
-          case ('(n,p)')
-            t % score_bins(j) = N_P
-            t % depletion_rx = .true.
-          case ('(n,d)')
-            t % score_bins(j) = N_D
-          case ('(n,t)')
-            t % score_bins(j) = N_T
-          case ('(n,3He)')
-            t % score_bins(j) = N_3HE
-          case ('(n,a)')
-            t % score_bins(j) = N_A
-            t % depletion_rx = .true.
-          case ('(n,2a)')
-            t % score_bins(j) = N_2A
-          case ('(n,3a)')
-            t % score_bins(j) = N_3A
-          case ('(n,2p)')
-            t % score_bins(j) = N_2P
-          case ('(n,pa)')
-            t % score_bins(j) = N_PA
-          case ('(n,t2a)')
-            t % score_bins(j) = N_T2A
-          case ('(n,d2a)')
-            t % score_bins(j) = N_D2A
-          case ('(n,pd)')
-            t % score_bins(j) = N_PD
-          case ('(n,pt)')
-            t % score_bins(j) = N_PT
-          case ('(n,da)')
-            t % score_bins(j) = N_DA
-
-          case default
-            ! First look for deprecated scores
-            if (starts_with(trim(score_name), 'scatter-') .or. &
-                 starts_with(trim(score_name), 'nu-scatter-') .or. &
-                 starts_with(trim(score_name), 'total-y') .or. &
-                 starts_with(trim(score_name), 'flux-y')) then
-              call fatal_error(trim(score_name) // " is no longer available.")
-            end if
-
-            ! Assume that user has specified an MT number
-            MT = int(str_to_int(score_name))
-
-            if (MT /= ERROR_INT) then
-              ! Specified score was an integer
-              if (MT > 1) then
-                t % score_bins(j) = MT
-              else
-                call fatal_error("Invalid MT on <scores>: " // trim(score_name))
-              end if
-
-            else
-              ! Specified score was not an integer
-              call fatal_error("Unknown scoring function: " // trim(score_name))
-            end if
-
-          end select
-
-          ! Do a check at the end (instead of for every case) to make sure
-          ! the tallies are compatible with MG mode where we have less detailed
-          ! nuclear data
-          if (.not. run_CE .and. t % score_bins(j) > 0) then
-            call fatal_error("Cannot tally " // trim(score_name) // &
-                             " reaction rate in multi-group mode")
-          end if
-        end do
-
-        t % n_score_bins = n_scores
-
-        ! Deallocate temporary string array of scores
-        deallocate(sarray)
-
-        ! Check that no duplicate scores exist
-        do j = 1, n_scores - 1
-          do k = j + 1, n_scores
-            if (t % score_bins(j) == t % score_bins(k)) then
-              call fatal_error("Duplicate score of type '" // trim(&
-                   reaction_name(t % score_bins(j))) // "' found in tally " &
-                   // trim(to_str(t % id)))
-            end if
-          end do
-        end do
 
         ! Check if tally is compatible with particle type
         if (photon_transport) then
-          if (t % find_filter(FILTER_PARTICLE) == 0) then
-            do j = 1, n_scores
+          if (particle_filter_index == 0) then
+            do j = 1, t % n_score_bins()
               select case (t % score_bins(j))
               case (SCORE_INVERSE_VELOCITY)
                 call fatal_error("Particle filter must be used with photon &
@@ -1075,18 +753,18 @@ contains
               end select
             end do
           else
-            select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
+            select type(filt => filters(particle_filter_index) % obj)
             type is (ParticleFilter)
               do l = 1, filt % n_bins
                 if (filt % particles(l) == ELECTRON .or. filt % particles(l) == POSITRON) then
-                  t % estimator = ESTIMATOR_ANALOG
+                  call t % set_estimator(ESTIMATOR_ANALOG)
                 end if
               end do
             end select
           end if
         else
-          if (t % find_filter(FILTER_PARTICLE) > 0) then
-            select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
+          if (particle_filter_index > 0) then
+            select type(filt => filters(particle_filter_index) % obj)
             type is (ParticleFilter)
               do l = 1, filt % n_bins
                 if (filt % particles(l) /= NEUTRON) then
@@ -1100,38 +778,39 @@ contains
         end if
       else
         call fatal_error("No <scores> specified on tally " &
-             // trim(to_str(t % id)) // ".")
+             // trim(to_str(t % id())) // ".")
       end if
 
       ! Check for a tally derivative.
       if (check_for_node(node_tal, "derivative")) then
-        ! Temporarily store the derivative id.
-        call get_node_value(node_tal, "derivative", t % deriv)
+        call get_node_value(node_tal, "derivative", deriv_id)
 
         ! Find the derivative with the given id, and store it's index.
-        do j = 1, size(tally_derivs)
-          if (tally_derivs(j) % id == t % deriv) then
-            t % deriv = j
+        do j = 0, n_tally_derivs() - 1
+          deriv => tally_deriv_c(j)
+          if (deriv % id == deriv_id) then
+            call t % set_deriv(j)
             ! Only analog or collision estimators are supported for differential
             ! tallies.
-            if (t % estimator == ESTIMATOR_TRACKLENGTH) then
-              t % estimator = ESTIMATOR_COLLISION
+            if (t % estimator() == ESTIMATOR_TRACKLENGTH) then
+              call t % set_estimator(ESTIMATOR_COLLISION)
             end if
             ! We found the derivative we were looking for; exit the do loop.
             exit
           end if
-          if (j == size(tally_derivs)) then
+          if (j == n_tally_derivs()) then
             call fatal_error("Could not find derivative " &
-                 // trim(to_str(t % deriv)) // " specified on tally " &
-                 // trim(to_str(t % id)))
+                 // trim(to_str(deriv_id)) // " specified on tally " &
+                 // trim(to_str(t % id())))
           end if
         end do
 
-        if (tally_derivs(t % deriv) % variable == DIFF_NUCLIDE_DENSITY &
-             .or. tally_derivs(t % deriv) % variable == DIFF_TEMPERATURE) then
-          if (any(t % nuclide_bins == -1)) then
-            if (t % find_filter(FILTER_ENERGYOUT) > 0) then
-              call fatal_error("Error on tally " // trim(to_str(t % id)) &
+        deriv => tally_deriv_c(t % deriv())
+        if (deriv % variable == DIFF_NUCLIDE_DENSITY &
+             .or. deriv % variable == DIFF_TEMPERATURE) then
+          do j = 1, t % n_nuclide_bins()
+            if (has_energyout .and. t % nuclide_bins(j) == -1) then
+              call fatal_error("Error on tally " // trim(to_str(t % id())) &
                    // ": Cannot use a 'nuclide_density' or 'temperature' &
                    &derivative on a tally with an outgoing energy filter and &
                    &'total' nuclide rate. Instead, tally each nuclide in the &
@@ -1141,176 +820,14 @@ contains
               ! (e.g. pertrubing moderator but only tallying fuel), but this
               ! case would be hard to check for by only reading inputs.
             end if
-          end if
+          end do
         end if
       end if
 
       ! If settings.xml trigger is turned on, create tally triggers
       if (trigger_on) then
-
-        ! Get list of trigger nodes for this tally
-        call get_node_list(node_tal, "trigger", node_trigger_list)
-
-        ! Initialize the number of triggers
-        n_user_trig = size(node_trigger_list)
-
-        ! Count the number of triggers needed for all scores including "all"
-        t % n_triggers = 0
-        COUNT_TRIGGERS: do user_trig_ind = 1, n_user_trig
-
-          ! Get pointer to trigger node
-          node_trigger = node_trigger_list(user_trig_ind)
-
-          ! Get scores for this trigger
-          if (check_for_node(node_trigger, "scores")) then
-            n_words = node_word_count(node_trigger, "scores")
-            allocate(sarray(n_words))
-            call get_node_array(node_trigger, "scores", sarray)
-          else
-            n_words = 1
-            allocate(sarray(n_words))
-            sarray(1) = "all"
-          end if
-
-          ! Count the number of scores for this trigger
-          do j = 1, n_words
-            score_name = trim(to_lower(sarray(j)))
-
-            if (score_name == "all") then
-              t % n_triggers = t % n_triggers + trigger_scores % size()
-            else
-              t % n_triggers = t % n_triggers + 1
-            end if
-
-          end do
-
-          deallocate(sarray)
-
-        end do COUNT_TRIGGERS
-
-        ! Allocate array of triggers for this tally
-        if (t % n_triggers > 0) then
-          allocate(t % triggers(t % n_triggers))
-        end if
-
-        ! Initialize overall trigger index for this tally to zero
-        trig_ind = 1
-
-        ! Create triggers for all scores specified on each trigger
-        TRIGGER_LOOP: do user_trig_ind = 1, n_user_trig
-
-          ! Get pointer to trigger node
-          node_trigger = node_trigger_list(user_trig_ind)
-
-          ! Get the trigger type - "variance", "std_dev" or "rel_err"
-          if (check_for_node(node_trigger, "type")) then
-            call get_node_value(node_trigger, "type", temp_str)
-            temp_str = to_lower(temp_str)
-          else
-            call fatal_error("Must specify trigger type for tally " // &
-                 trim(to_str(t % id)) // " in tally XML file.")
-          end if
-
-          ! Get the convergence threshold for the trigger
-          if (check_for_node(node_trigger, "threshold")) then
-            call get_node_value(node_trigger, "threshold", threshold)
-          else
-            call fatal_error("Must specify trigger threshold for tally " // &
-                 trim(to_str(t % id)) // " in tally XML file.")
-          end if
-
-          ! Get list scores for this trigger
-          if (check_for_node(node_trigger, "scores")) then
-            n_words = node_word_count(node_trigger, "scores")
-            allocate(sarray(n_words))
-            call get_node_array(node_trigger, "scores", sarray)
-          else
-            n_words = 1
-            allocate(sarray(n_words))
-            sarray(1) = "all"
-          end if
-
-          ! Create a trigger for each score
-          SCORE_LOOP: do j = 1, n_words
-            score_name = trim(to_lower(sarray(j)))
-
-            ! Expand "all" to include TriggerObjects for each score in tally
-            if (score_name == "all") then
-
-              ! Loop over all tally scores
-              i_elem = 0
-              do
-                ! Move to next score
-                call trigger_scores % next_entry(elem, i_elem)
-                if (i_elem == 0) exit
-
-                score_name = trim(elem % key)
-
-                ! Store the score name and index in the trigger
-                t % triggers(trig_ind) % score_name = score_name
-                t % triggers(trig_ind) % score_index = elem % value
-
-                ! Set the trigger convergence threshold type
-                select case (temp_str)
-                case ('std_dev')
-                  t % triggers(trig_ind) % type = STANDARD_DEVIATION
-                case ('variance')
-                  t % triggers(trig_ind) % type = VARIANCE
-                case ('rel_err')
-                  t % triggers(trig_ind) % type = RELATIVE_ERROR
-                case default
-                  call fatal_error("Unknown trigger type " // &
-                       trim(temp_str) // " in tally " // trim(to_str(t % id)))
-                end select
-
-                ! Store the trigger convergence threshold
-                t % triggers(trig_ind) % threshold = threshold
-
-                ! Increment the overall trigger index
-                trig_ind = trig_ind + 1
-              end do
-
-            ! Scores other than the "all" placeholder
-            else
-
-              ! Store the score name and index
-              t % triggers(trig_ind) % score_name = trim(score_name)
-              t % triggers(trig_ind) % score_index = &
-                   trigger_scores % get(trim(score_name))
-
-              ! Check if an invalid score was set for the trigger
-              if (t % triggers(trig_ind) % score_index == 0) then
-                call fatal_error("The trigger score " // trim(score_name) // &
-                     " is not set for tally " // trim(to_str(t % id)))
-              end if
-
-              ! Store the trigger convergence threshold
-              t % triggers(trig_ind) % threshold = threshold
-
-              ! Set the trigger convergence threshold type
-              select case (temp_str)
-              case ('std_dev')
-                t % triggers(trig_ind) % type = STANDARD_DEVIATION
-              case ('variance')
-                t % triggers(trig_ind) % type = VARIANCE
-              case ('rel_err')
-                t % triggers(trig_ind) % type = RELATIVE_ERROR
-              case default
-                call fatal_error("Unknown trigger type " // trim(temp_str) // &
-                     " in tally " // trim(to_str(t % id)))
-              end select
-
-              ! Increment the overall trigger index
-              trig_ind = trig_ind + 1
-            end if
-          end do SCORE_LOOP
-
-          ! Deallocate the list of tally scores used to create triggers
-          deallocate(sarray)
-        end do TRIGGER_LOOP
-
-        ! Deallocate dictionary of scores/indices used to populate triggers
-        call trigger_scores % clear()
+        !TODO: off-by-one
+        call tally_init_triggers(t % ptr, i_start + i - 1 - 1, node_tal % ptr)
       end if
 
       ! =======================================================================
@@ -1322,33 +839,33 @@ contains
         call get_node_value(node_tal, "estimator", temp_str)
         select case(trim(temp_str))
         case ('analog')
-          t % estimator = ESTIMATOR_ANALOG
+          call t % set_estimator(ESTIMATOR_ANALOG)
 
         case ('tracklength', 'track-length', 'pathlength', 'path-length')
           ! If the estimator was set to an analog estimator, this means the
           ! tally needs post-collision information
-          if (t % estimator == ESTIMATOR_ANALOG) then
+          if (t % estimator() == ESTIMATOR_ANALOG) then
             call fatal_error("Cannot use track-length estimator for tally " &
-                 // to_str(t % id))
+                 // to_str(t % id()))
           end if
 
           ! Set estimator to track-length estimator
-          t % estimator = ESTIMATOR_TRACKLENGTH
+          call t % set_estimator(ESTIMATOR_TRACKLENGTH)
 
         case ('collision')
           ! If the estimator was set to an analog estimator, this means the
           ! tally needs post-collision information
-          if (t % estimator == ESTIMATOR_ANALOG) then
+          if (t % estimator() == ESTIMATOR_ANALOG) then
             call fatal_error("Cannot use collision estimator for tally " &
-                 // to_str(t % id))
+                 // to_str(t % id()))
           end if
 
           ! Set estimator to collision estimator
-          t % estimator = ESTIMATOR_COLLISION
+          call t % set_estimator(ESTIMATOR_COLLISION)
 
         case default
           call fatal_error("Invalid estimator '" // trim(temp_str) &
-               // "' on tally " // to_str(t % id))
+               // "' on tally " // to_str(t % id()))
         end select
       end if
 

src/material.cpp

@@ -636,7 +636,7 @@ void Material::init_nuclide_index()
   int n = settings::run_CE ?
     data::nuclides.size() : data::nuclides_MG.size();
   mat_nuclide_index_.resize(n);
-  std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), -1);
+  std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), C_NONE);
   for (int i = 0; i < nuclide_.size(); ++i) {
     mat_nuclide_index_[nuclide_[i]] = i;
   }
@@ -1135,11 +1135,6 @@ extern "C" {
     return model::materials[i_mat - 1]->density_gpcc_;
   }
 
-  int material_nuclide_index(int32_t i_mat, int i_nuc)
-  {
-    return model::materials[i_mat - 1]->mat_nuclide_index_[i_nuc - 1] + 1;
-  }
-
   void material_calculate_xs(const Particle* p)
   {
     model::materials[p->material - 1]->calculate_xs(*p);

src/material_header.F90

@@ -11,7 +11,6 @@ module material_header
   public :: material_id
   public :: material_nuclide
   public :: material_nuclide_size
-  public :: material_nuclide_index
   public :: material_atom_density
   public :: material_density_gpcc
 
@@ -40,13 +39,6 @@ module material_header
       integer(C_INT) :: n
     end function
 
-    function material_nuclide_index(i_mat, i_nuc) bind(C) result(idx)
-      import C_INT32_T, C_INT
-      integer(C_INT32_T), value :: i_mat
-      integer(C_INT), value :: i_nuc
-      integer(C_INT) :: idx
-    end function
-
     function material_atom_density(i_mat, idx) bind(C) result(density)
       import C_INT32_T, C_INT, C_DOUBLE
       integer(C_INT32_T), value :: i_mat

src/mgxs.cpp

@@ -424,7 +424,8 @@ Mgxs::combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scala
 //==============================================================================
 
 double
-Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
+Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
+  const int* dg)
 {
   // This method assumes that the temperature and angle indices are set
 #ifdef _OPENMP

src/mgxs_interface.F90

@@ -35,53 +35,6 @@ module mgxs_interface
       real(C_DOUBLE)                    :: awr
     end function get_awr_c
 
-    function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
-         bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT),  value,   intent(in) :: index
-      integer(C_INT),  value,   intent(in) :: xstype
-      integer(C_INT),  value,   intent(in) :: gin
-      integer(C_INT), optional, intent(in) :: gout
-      real(C_DOUBLE), optional, intent(in) :: mu
-      integer(C_INT), optional, intent(in) :: dg
-      real(C_DOUBLE)                       :: val
-    end function get_nuclide_xs_c
-
-    function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
-         bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT),  value,   intent(in) :: index
-      integer(C_INT),  value,   intent(in) :: xstype
-      integer(C_INT),  value,   intent(in) :: gin
-      integer(C_INT), optional, intent(in) :: gout
-      real(C_DOUBLE), optional, intent(in) :: mu
-      integer(C_INT), optional, intent(in) :: dg
-      real(C_DOUBLE)                       :: val
-    end function get_macro_xs_c
-
-    subroutine set_nuclide_angle_index_c(index, uvw) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: index
-      real(C_DOUBLE),        intent(in) :: uvw(1:3)
-    end subroutine set_nuclide_angle_index_c
-
-    subroutine set_macro_angle_index_c(index, uvw) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: index
-      real(C_DOUBLE),        intent(in) :: uvw(1:3)
-    end subroutine set_macro_angle_index_c
-
-    subroutine set_nuclide_temperature_index_c(index, sqrtkT) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: index
-      real(C_DOUBLE), value, intent(in) :: sqrtkT
-    end subroutine set_nuclide_temperature_index_c
-
   end interface
 
   ! Number of energy groups

src/mgxs_interface.cpp

@@ -231,12 +231,13 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
 //==============================================================================
 
 double
-get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
+get_nuclide_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg)
 {
   int gout_c;
-  int* gout_c_p;
+  const int* gout_c_p;
   int dg_c;
-  int* dg_c_p;
+  const int* dg_c_p;
   if (gout != nullptr) {
     gout_c = *gout - 1;
     gout_c_p = &gout_c;
@@ -255,12 +256,13 @@ get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
 //==============================================================================
 
 double
-get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
+get_macro_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg)
 {
   int gout_c;
-  int* gout_c_p;
+  const int* gout_c_p;
   int dg_c;
-  int* dg_c_p;
+  const int* dg_c_p;
   if (gout != nullptr) {
     gout_c = *gout - 1;
     gout_c_p = &gout_c;
@@ -276,33 +278,6 @@ get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
   return data::macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
 }
 
-//==============================================================================
-
-void
-set_nuclide_angle_index_c(int index, const double uvw[3])
-{
-  // Update the values
-  data::nuclides_MG[index - 1].set_angle_index(uvw);
-}
-
-//==============================================================================
-
-void
-set_macro_angle_index_c(int index, const double uvw[3])
-{
-  // Update the values
-  data::macro_xs[index - 1].set_angle_index(uvw);
-}
-
-//==============================================================================
-
-void
-set_nuclide_temperature_index_c(int index, double sqrtkT)
-{
-  // Update the values
-  data::nuclides_MG[index - 1].set_temperature_index(sqrtkT);
-}
-
 //==============================================================================
 // General Mgxs methods
 //==============================================================================

src/nuclide.cpp

@@ -274,7 +274,7 @@ void Nuclide::create_derived()
     xs_.emplace_back(shape, 0.0);
   }
 
-  reaction_index_.fill(-1);
+  reaction_index_.fill(C_NONE);
   for (int i = 0; i < reactions_.size(); ++i) {
     const auto& rx {reactions_[i]};
 
@@ -997,7 +997,7 @@ extern "C" void multipole_deriv_eval(Nuclide* nuc, double E, double sqrtkT,
   std::tie(*sig_s, *sig_a, *sig_f) = nuc->multipole_->evaluate_deriv(E, sqrtkT);
 }
 
-extern "C" bool multipole_in_range(Nuclide* nuc, double E)
+extern "C" bool multipole_in_range(const Nuclide* nuc, double E)
 {
   return nuc->multipole_ && E >= nuc->multipole_->E_min_&&
     E <= nuc->multipole_->E_max_;

src/output.F90

@@ -20,8 +20,6 @@ module output
   use tally_header
   use tally_derivative_header
   use tally_filter
-  use tally_filter_mesh, only: MeshFilter
-  use tally_filter_header, only: TallyFilterMatch
 
   implicit none
 
@@ -34,6 +32,9 @@ module output
       import Particle
       type(Particle), intent(in) :: p
     end subroutine
+
+    subroutine write_tallies() bind(C)
+    end subroutine
   end interface
 
 contains
@@ -92,249 +93,4 @@ contains
 
   end subroutine header
 
-!===============================================================================
-! WRITE_TALLIES creates an output file and writes out the mean values of all
-! tallies and their standard deviations
-!===============================================================================
-
-  subroutine write_tallies() bind(C)
-
-    integer :: i            ! index in tallies array
-    integer :: j            ! level in tally hierarchy
-    integer :: k            ! loop index for scoring bins
-    integer :: n            ! loop index for nuclides
-    integer :: h            ! loop index for tally filters
-    integer :: indent       ! number of spaces to preceed output
-    integer :: filter_index ! index in results array for filters
-    integer :: score_index  ! scoring bin index
-    integer :: i_nuclide    ! index in nuclides array
-    integer :: unit_tally   ! tallies.out file unit
-    integer :: nr           ! number of realizations
-    real(8) :: t_value      ! t-values for confidence intervals
-    real(8) :: alpha        ! significance level for CI
-    real(8) :: x(2)         ! mean and standard deviation
-    character(MAX_FILE_LEN) :: filename                    ! name of output file
-    character(36)           :: score_names(N_SCORE_TYPES)  ! names of scoring function
-    character(36)           :: score_name                  ! names of scoring function
-                                                           ! to be applied at write-time
-    integer, allocatable :: filter_bins(:)
-    character(MAX_WORD_LEN) :: temp_name
-
-    ! Skip if there are no tallies
-    if (n_tallies == 0) return
-
-    allocate(filter_bins(n_filters))
-
-    ! Initialize names for scores
-    score_names(abs(SCORE_FLUX))               = "Flux"
-    score_names(abs(SCORE_TOTAL))              = "Total Reaction Rate"
-    score_names(abs(SCORE_SCATTER))            = "Scattering Rate"
-    score_names(abs(SCORE_NU_SCATTER))         = "Scattering Production Rate"
-    score_names(abs(SCORE_ABSORPTION))         = "Absorption Rate"
-    score_names(abs(SCORE_FISSION))            = "Fission Rate"
-    score_names(abs(SCORE_NU_FISSION))         = "Nu-Fission Rate"
-    score_names(abs(SCORE_KAPPA_FISSION))      = "Kappa-Fission Rate"
-    score_names(abs(SCORE_EVENTS))             = "Events"
-    score_names(abs(SCORE_DECAY_RATE))         = "Decay Rate"
-    score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
-    score_names(abs(SCORE_PROMPT_NU_FISSION))  = "Prompt-Nu-Fission Rate"
-    score_names(abs(SCORE_INVERSE_VELOCITY))   = "Flux-Weighted Inverse Velocity"
-    score_names(abs(SCORE_FISS_Q_PROMPT))      = "Prompt fission power"
-    score_names(abs(SCORE_FISS_Q_RECOV))       = "Recoverable fission power"
-    score_names(abs(SCORE_CURRENT))            = "Current"
-
-    ! Create filename for tally output
-    filename = trim(path_output) // "tallies.out"
-
-    ! Open tally file for writing
-    open(FILE=filename, NEWUNIT=unit_tally, STATUS='replace', ACTION='write')
-
-    ! Calculate t-value for confidence intervals
-    if (confidence_intervals) then
-      alpha = ONE - CONFIDENCE_LEVEL
-      t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1)
-    else
-      t_value = ONE
-    end if
-
-    TALLY_LOOP: do i = 1, n_tallies
-      associate (t => tallies(i) % obj)
-      nr = t % n_realizations
-
-      if (confidence_intervals) then
-        ! Calculate t-value for confidence intervals
-        alpha = ONE - CONFIDENCE_LEVEL
-        t_value = t_percentile(ONE - alpha/TWO, nr - 1)
-      else
-        t_value = ONE
-      end if
-
-      ! Write header block
-      if (t % name == "") then
-        call header("TALLY " // trim(to_str(t % id)), 1, unit=unit_tally)
-      else
-        call header("TALLY " // trim(to_str(t % id)) // ": " &
-             // trim(t % name), 1, unit=unit_tally)
-      endif
-
-      ! Write derivative information.
-      if (t % deriv /= NONE) then
-        associate(deriv => tally_derivs(t % deriv))
-          select case (deriv % variable)
-          case (DIFF_DENSITY)
-            write(unit=unit_tally, fmt="(' Density derivative  Material ',A)") &
-                 to_str(deriv % diff_material)
-          case (DIFF_NUCLIDE_DENSITY)
-            write(unit=unit_tally, fmt="(' Nuclide density derivative  &
-                 &Material ',A,'  Nuclide ',A)") &
-                 trim(to_str(deriv % diff_material)), &
-                 trim(nuclides(deriv % diff_nuclide) % name)
-          case (DIFF_TEMPERATURE)
-            write(unit=unit_tally, fmt="(' Temperature derivative  Material ',&
-                 &A)") to_str(deriv % diff_material)
-          case default
-            call fatal_error("Differential tally dependent variable for tally "&
-                 // trim(to_str(t % id)) // " not defined in output.F90.")
-          end select
-        end associate
-      end if
-
-      ! WARNING: Admittedly, the logic for moving for printing results is
-      ! extremely confusing and took quite a bit of time to get correct. The
-      ! logic is structured this way since it is not practical to have a do
-      ! loop for each filter variable (given that only a few filters are likely
-      ! to be used for a given tally.
-
-      ! Initialize bins, filter level, and indentation
-      do h = 1, size(t % filter)
-        filter_bins(t % filter(h)) = 0
-      end do
-      j = 1
-      indent = 0
-
-      print_bin: do
-        find_bin: do
-          ! Check for no filters
-          if (size(t % filter) == 0) exit find_bin
-
-          ! Increment bin combination
-          filter_bins(t % filter(j)) = filter_bins(t % filter(j)) + 1
-
-          ! =================================================================
-          ! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
-
-          if (filter_bins(t % filter(j)) > &
-               filters(t % filter(j)) % obj % n_bins) then
-            ! If this is the first filter, then exit
-            if (j == 1) exit print_bin
-
-            filter_bins(t % filter(j)) = 0
-            j = j - 1
-            indent = indent - 2
-
-            ! =================================================================
-            ! VALID BIN -- WRITE FILTER INFORMATION OR EXIT TO WRITE RESULTS
-
-          else
-            ! Check if this is last filter
-            if (j == size(t % filter)) exit find_bin
-
-            ! Print current filter information
-            write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
-                 trim(filters(t % filter(j)) % obj % &
-                 text_label(filter_bins(t % filter(j))))
-            indent = indent + 2
-            j = j + 1
-          end if
-
-        end do find_bin
-
-        ! Print filter information
-        if (size(t % filter) > 0) then
-          write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
-               trim(filters(t % filter(j)) % obj % &
-               text_label(filter_bins(t % filter(j))))
-        end if
-
-        ! Determine scoring index for this bin combination -- note that unlike
-        ! in the score_tally subroutine, we have to use max(bins,1) since all
-        ! bins below the lowest filter level will be zeros
-
-        filter_index = 1
-        do h = 1, size(t % filter)
-          filter_index = filter_index &
-               + (max(filter_bins(t % filter(h)) ,1) - 1) * t % stride(h)
-        end do
-
-        ! Write results for this filter bin combination
-        score_index = 0
-        if (size(t % filter) > 0) indent = indent + 2
-        do n = 1, t % n_nuclide_bins
-          ! Write label for nuclide
-          i_nuclide = t % nuclide_bins(n)
-          if (i_nuclide == -1) then
-            write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
-                 "Total Material"
-          else
-            if (run_CE) then
-              write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
-                   trim(nuclides(i_nuclide) % name)
-            else
-              call get_name_c(i_nuclide, len(temp_name), temp_name)
-              write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
-                   trim(temp_name)
-            end if
-          end if
-
-          indent = indent + 2
-          do k = 1, t % n_score_bins
-            score_index = score_index + 1
-
-            associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :))
-
-            if (t % score_bins(k) > 0) then
-              score_name = reaction_name(t % score_bins(k))
-            else
-              score_name = score_names(abs(t % score_bins(k)))
-            end if
-            x(:) = mean_stdev(r(:, score_index, filter_index), nr)
-            write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
-                 repeat(" ", indent), score_name, &
-                 to_str(x(1)), trim(to_str(t_value * x(2)))
-            end associate
-          end do
-          indent = indent - 2
-
-        end do
-        indent = indent - 2
-
-        if (size(t % filter) == 0) exit print_bin
-
-      end do print_bin
-
-      end associate
-    end do TALLY_LOOP
-
-    close(UNIT=unit_tally)
-
-  end subroutine write_tallies
-
-!===============================================================================
-! MEAN_STDEV computes the sample mean and standard deviation of the mean of a
-! single tally score
-!===============================================================================
-
-  pure function mean_stdev(result_, n) result(x)
-    real(8), intent(in) :: result_(2) ! sum and sum-of-squares
-    integer, intent(in) :: n          ! number of realizations
-    real(8)  :: x(2)                  ! mean and standard deviation
-
-    x(1) = result_(1) / n
-    if (n > 1) then
-      x(2) = sqrt((result_(2) / n - x(1)*x(1))/(n - 1))
-    else
-      x(2) = ZERO
-    end if
-  end function mean_stdev
-
 end module output

src/output.cpp

@@ -6,7 +6,9 @@
 #include <iomanip>  // for setw, setprecision, put_time
 #include <ios> // for fixed, scientific, left
 #include <iostream>
+#include <fstream>
 #include <sstream>
+#include <unordered_map>
 #include <utility> // for pair
 
 #include <omp.h>
@@ -21,12 +23,18 @@
 #include "openmc/lattice.h"
 #include "openmc/math_functions.h"
 #include "openmc/message_passing.h"
+#include "openmc/mgxs_interface.h"
+#include "openmc/nuclide.h"
 #include "openmc/plot.h"
+#include "openmc/reaction.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/surface.h"
-#include "openmc/timer.h"
+#include "openmc/tallies/derivative.h"
+#include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/tallies/tally_scoring.h"
+#include "openmc/timer.h"
 
 namespace openmc {
 
@@ -87,8 +95,8 @@ void title()
 
 //==============================================================================
 
-void
-header(const char* msg, int level) {
+std::string
+header(const char* msg) {
   // Determine how many times to repeat the '=' character.
   int n_prefix = (63 - strlen(msg)) / 2;
   int n_suffix = n_prefix;
@@ -105,10 +113,18 @@ header(const char* msg, int level) {
   out << ">     " << upper << "     <";
   for (int i = 0; i < n_suffix; i++) out << '=';
 
+  return out.str();
+}
+
+std::string header(const std::string& msg) {return header(msg.c_str());}
+
+void
+header(const char* msg, int level) {
+  auto out = header(msg);
+
   // Print header based on verbosity level.
-  if (settings::verbosity >= level) {
-    std::cout << '\n' << out.str() << "\n\n";
-  }
+  if (settings::verbosity >= level)
+    std::cout << '\n' << out << "\n\n";
 }
 
 //==============================================================================
@@ -502,7 +518,8 @@ void print_runtime()
 
 //==============================================================================
 
-std::pair<double, double> mean_stdev(const double* x, int n)
+std::pair<double, double>
+mean_stdev(const double* x, int n)
 {
   double mean = x[RESULT_SUM] / n;
   double stdev = n > 1 ? std::sqrt((x[RESULT_SUM_SQ]/n
@@ -580,4 +597,138 @@ void print_results()
   std::cout.flags(f);
 }
 
+//==============================================================================
+
+const std::unordered_map<int, const char*> score_names = {
+  {SCORE_FLUX,               "Flux"},
+  {SCORE_TOTAL,              "Total Reaction Rate"},
+  {SCORE_SCATTER,            "Scattering Rate"},
+  {SCORE_NU_SCATTER,         "Scattering Production Rate"},
+  {SCORE_ABSORPTION,         "Absorption Rate"},
+  {SCORE_FISSION,            "Fission Rate"},
+  {SCORE_NU_FISSION,         "Nu-Fission Rate"},
+  {SCORE_KAPPA_FISSION,      "Kappa-Fission Rate"},
+  {SCORE_EVENTS,             "Events"},
+  {SCORE_DECAY_RATE,         "Decay Rate"},
+  {SCORE_DELAYED_NU_FISSION, "Delayed-Nu-Fission Rate"},
+  {SCORE_PROMPT_NU_FISSION,  "Prompt-Nu-Fission Rate"},
+  {SCORE_INVERSE_VELOCITY,   "Flux-Weighted Inverse Velocity"},
+  {SCORE_FISS_Q_PROMPT,      "Prompt fission power"},
+  {SCORE_FISS_Q_RECOV,       "Recoverable fission power"},
+  {SCORE_CURRENT,            "Current"},
+};
+
+//! Create an ASCII output file showing all tally results.
+
+extern "C" void
+write_tallies()
+{
+  if (model::tallies.empty()) return;
+
+  // Open the tallies.out file.
+  std::ofstream tallies_out;
+  tallies_out.open("tallies.out", std::ios::out | std::ios::trunc);
+  tallies_out << std::setprecision(6);
+
+  // Loop over each tally.
+  for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
+    const auto& tally {*model::tallies[i_tally]};
+    auto results = tally_results(i_tally);
+    // TODO: get this directly from the tally object when it's been translated
+    int32_t n_realizations;
+    auto err = openmc_tally_get_n_realizations(i_tally+1, &n_realizations);
+
+    // Calculate t-value for confidence intervals
+    double t_value = 1;
+    if (settings::confidence_intervals) {
+      auto alpha = 1 - CONFIDENCE_LEVEL;
+      t_value = t_percentile(1 - alpha*0.5, n_realizations - 1);
+    }
+
+    // Write header block.
+    std::string tally_header("TALLY " + std::to_string(tally.id_));
+    if (!tally.name_.empty()) tally_header += ": " + tally.name_;
+    tallies_out << header(tally_header) << "\n\n";
+
+    // Write derivative information.
+    if (tally.deriv_ != C_NONE) {
+      const auto& deriv {model::tally_derivs[tally.deriv_]};
+      switch (deriv.variable) {
+      case DIFF_DENSITY:
+        tallies_out << " Density derivative  Material "
+          << std::to_string(deriv.diff_material) << "\n";
+        break;
+      case DIFF_NUCLIDE_DENSITY:
+        tallies_out << " Nuclide density derivative  Material "
+          << std::to_string(deriv.diff_material) << "  Nuclide "
+          // TODO: off-by-one
+          << data::nuclides[deriv.diff_nuclide-1]->name_ << "\n";
+        break;
+      case DIFF_TEMPERATURE:
+        tallies_out << " Temperature derivative  Material "
+          << std::to_string(deriv.diff_material) << "\n";
+        break;
+      default:
+        fatal_error("Differential tally dependent variable for tally "
+          + std::to_string(tally.id_) + " not defined in output.cpp");
+      }
+    }
+
+    // Loop over all filter bin combinations.
+    auto filter_iter = FilterBinIter(tally, false);
+    auto end = FilterBinIter(tally, true);
+    for (; filter_iter != end; ++filter_iter) {
+      auto filter_index = filter_iter.index_;
+
+      // Print info about this combination of filter bins.  The stride check
+      // prevents redundant output.
+      int indent = 0;
+      for (auto i = 0; i < tally.filters().size(); ++i) {
+        if ((filter_index-1) % tally.strides(i) == 0) {
+          auto i_filt = tally.filters(i);
+          const auto& filt {*model::tally_filters[i_filt]};
+          auto& match {simulation::filter_matches[i_filt]};
+          tallies_out << std::string(indent+1, ' ')
+            // TODO: off-by-one
+            << filt.text_label(match.i_bin_+1) << "\n";
+        }
+        indent += 2;
+      }
+
+      // Loop over all nuclide and score combinations.
+      int score_index = 0;
+      for (auto i_nuclide : tally.nuclides_) {
+        // Write label for this nuclide bin.
+        if (i_nuclide == -1) {
+          tallies_out << std::string(indent+1, ' ') << "Total Material\n";
+        } else {
+          if (settings::run_CE) {
+            tallies_out << std::string(indent+1, ' ')
+              << data::nuclides[i_nuclide]->name_ << "\n";
+          } else {
+            tallies_out << std::string(indent+1, ' ')
+              << data::nuclides_MG[i_nuclide].name << "\n";
+          }
+        }
+
+        // Write the score, mean, and uncertainty.
+        indent += 2;
+        for (auto score : tally.scores_) {
+          std::string score_name = score > 0 ? reaction_name(score)
+            : score_names.at(score);
+          double mean, stdev;
+          //TODO: off-by-one
+          std::tie(mean, stdev) = mean_stdev(
+            &results(filter_index-1, score_index, 0), n_realizations);
+          tallies_out << std::string(indent+1, ' ')  << std::left
+            << std::setw(36) << score_name << " " << mean << " +/- "
+            << t_value * stdev << "\n";
+          score_index += 1;
+        }
+        indent -= 2;
+      }
+    }
+  }
+}
+
 } // namespace openmc

src/settings.F90

@@ -51,8 +51,6 @@ module settings
   integer(C_INT32_T), bind(C) :: gen_per_batch ! # of generations per batch
 
   integer(C_INT), bind(C) :: n_max_batches             ! max # of batches
-  integer(C_INT), bind(C, name='trigger_batch_interval') :: n_batch_interval      ! batch interval for triggers
-  logical(C_BOOL), bind(C, name='trigger_predict') :: pred_batches   ! predict batches for triggers
   logical(C_BOOL), bind(C) :: trigger_on      ! flag for turning triggers on/off
 
   logical(C_BOOL), bind(C) :: entropy_on

src/settings.cpp

@@ -25,6 +25,7 @@
 #include "openmc/simulation.h"
 #include "openmc/source.h"
 #include "openmc/string_utils.h"
+#include "openmc/tallies/trigger.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -105,13 +106,6 @@ int verbosity {7};
 double weight_cutoff {0.25};
 double weight_survive {1.0};
 
-// TODO: Move to separate file
-struct KTrigger {
-  int type;
-  double threshold;
-};
-extern "C" KTrigger keff_trigger;
-
 } // namespace settings
 
 //==============================================================================
@@ -160,11 +154,11 @@ void get_run_parameters(pugi::xml_node node_base)
       if (check_for_node(node_keff_trigger, "type")) {
         auto temp = get_node_value(node_keff_trigger, "type", true, true);
         if (temp == "std_dev") {
-          keff_trigger.type = STANDARD_DEVIATION;
+          keff_trigger.metric = TriggerMetric::standard_deviation;
         } else if (temp == "variance") {
-          keff_trigger.type = VARIANCE;
+          keff_trigger.metric = TriggerMetric::variance;
         } else if (temp == "rel_err") {
-          keff_trigger.type = RELATIVE_ERROR;
+          keff_trigger.metric = TriggerMetric::relative_error;
         } else {
           fatal_error("Unrecognized keff trigger type " + temp);
         }

src/simulation.F90

@@ -4,7 +4,6 @@ module simulation
 
   use nuclide_header,  only: micro_xs, n_nuclides
   use photon_header,   only: micro_photon_xs, n_elements
-  use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer
 
   implicit none
   private
@@ -17,18 +16,10 @@ contains
 
   subroutine simulation_init_f() bind(C)
 
-    integer :: i
-
 !$omp parallel
     ! Allocate array for microscopic cross section cache
     allocate(micro_xs(n_nuclides))
     allocate(micro_photon_xs(n_elements))
-
-    ! Allocate array for matching filter bins
-    allocate(filter_matches(n_filters))
-    do i = 1, n_filters
-      filter_matches(i) % ptr = filter_match_pointer(i - 1)
-    end do
 !$omp end parallel
 
   end subroutine
@@ -39,12 +30,10 @@ contains
 !===============================================================================
 
   subroutine simulation_finalize_f() bind(C)
-
     ! Free up simulation-specific memory
 !$omp parallel
-    deallocate(micro_xs, micro_photon_xs, filter_matches)
+    deallocate(micro_xs, micro_photon_xs)
 !$omp end parallel
-
   end subroutine
 
 end module simulation

src/simulation.cpp

@@ -17,6 +17,7 @@
 #include "openmc/timer.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/tallies/trigger.h"
 
 #include <omp.h>
 
@@ -28,7 +29,6 @@ namespace openmc {
 // data/functions from Fortran side
 extern "C" void accumulate_tallies();
 extern "C" void allocate_tally_results();
-extern "C" void check_triggers();
 extern "C" void init_tally_routines();
 extern "C" void setup_active_tallies();
 extern "C" void simulation_init_f();
@@ -138,11 +138,6 @@ int openmc_simulation_finalize()
   simulation::time_active.stop();
   simulation::time_finalize.start();
 
-#pragma omp parallel
-  {
-    simulation::filter_matches.clear();
-  }
-
   // Deallocate Fortran variables, set tallies to inactive
   for (auto& mat : model::materials) {
     mat->mat_nuclide_index_.clear();
@@ -159,6 +154,11 @@ int openmc_simulation_finalize()
   // Write tally results to tallies.out
   if (settings::output_tallies && mpi::master) write_tallies();
 
+#pragma omp parallel
+  {
+    simulation::filter_matches.clear();
+  }
+
   // Deactivate all tallies
   for (int i = 1; i <= n_tallies; ++i) {
     openmc_tally_set_active(i, false);
@@ -540,7 +540,7 @@ void calculate_work()
 #ifdef OPENMC_MPI
 void broadcast_results() {
   // Broadcast tally results so that each process has access to results
-  for (int i = 1; i <= n_tallies; ++i) {
+  for (int i = 0; i < n_tallies; ++i) {
     // Create a new datatype that consists of all values for a given filter
     // bin and then use that to broadcast. This is done to minimize the
     // chance of the 'count' argument of MPI_BCAST exceeding 2**31

src/simulation_header.F90

@@ -23,11 +23,6 @@ module simulation_header
   integer(C_INT), bind(C) :: total_gen         ! total number of generations simulated
   logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx?
 
-  ! ============================================================================
-  ! TALLY PRECISION TRIGGER VARIABLES
-
-  logical(C_BOOL), bind(C) :: satisfy_triggers  ! whether triggers are satisfied
-
   integer(C_INT64_T), bind(C) :: work         ! number of particles per processor
 
   ! ============================================================================

src/state_point.F90

@@ -13,19 +13,10 @@ module state_point
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
-  use endf,               only: reaction_name
-  use error,              only: fatal_error, warning, write_message
   use hdf5_interface
-  use message_passing
-  use mgxs_interface
-  use nuclide_header,     only: nuclides
   use settings
-  use simulation_header
   use string,             only: to_str
   use tally_header
-  use tally_filter_header
-  use tally_derivative_header, only: tally_derivs
 
   implicit none
 
@@ -38,174 +29,18 @@ contains
   subroutine statepoint_write_f(file_id) bind(C)
     integer(HID_T), value :: file_id
 
-    integer :: i, j
-    integer :: i_xs
-    integer, allocatable :: id_array(:)
-    integer(HID_T) :: tallies_group, tally_group, &
-                      filters_group, filter_group, derivs_group, &
-                      deriv_group
-    character(MAX_WORD_LEN), allocatable :: str_array(:)
-    character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
+    integer :: i
+    integer(HID_T) :: tallies_group, tally_group
 
     ! Open tallies group
     tallies_group = open_group(file_id, "tallies")
 
-    ! Write information for derivatives.
-    if (size(tally_derivs) > 0) then
-      derivs_group = create_group(tallies_group, "derivatives")
-      do i = 1, size(tally_derivs)
-        associate(deriv => tally_derivs(i))
-          deriv_group = create_group(derivs_group, "derivative " &
-               // trim(to_str(deriv % id)))
-          select case (deriv % variable)
-          case (DIFF_DENSITY)
-            call write_dataset(deriv_group, "independent variable", "density")
-            call write_dataset(deriv_group, "material", deriv % diff_material)
-          case (DIFF_NUCLIDE_DENSITY)
-            call write_dataset(deriv_group, "independent variable", &
-                 "nuclide_density")
-            call write_dataset(deriv_group, "material", deriv % diff_material)
-            call write_dataset(deriv_group, "nuclide", &
-                 nuclides(deriv % diff_nuclide) % name)
-          case (DIFF_TEMPERATURE)
-            call write_dataset(deriv_group, "independent variable", &
-                 "temperature")
-            call write_dataset(deriv_group, "material", deriv % diff_material)
-          case default
-            call fatal_error("Independent variable for derivative " &
-                 // trim(to_str(deriv % id)) // " not defined in &
-                 &state_point.F90.")
-          end select
-          call close_group(deriv_group)
-        end associate
-      end do
-      call close_group(derivs_group)
-    end if
-
-    ! Write number of filters
-    filters_group = create_group(tallies_group, "filters")
-    call write_attribute(filters_group, "n_filters", n_filters)
-
-    if (n_filters > 0) then
-      ! Write IDs of filters
-      allocate(id_array(n_filters))
-      do i = 1, n_filters
-        id_array(i) = filters(i) % obj % id
-      end do
-      call write_attribute(filters_group, "ids", id_array)
-      deallocate(id_array)
-
-      ! Write filter information
-      FILTER_LOOP: do i = 1, n_filters
-        filter_group = create_group(filters_group, "filter " // &
-             trim(to_str(filters(i) % obj % id)))
-        call filters(i) % obj % to_statepoint(filter_group)
-        call close_group(filter_group)
-      end do FILTER_LOOP
-    end if
-
-    call close_group(filters_group)
-
-      ! Write number of tallies
-    call write_attribute(tallies_group, "n_tallies", n_tallies)
-
-    if (n_tallies > 0) then
-      ! Write array of tally IDs
-      allocate(id_array(n_tallies))
-      do i = 1, n_tallies
-        id_array(i) = tallies(i) % obj % id
-      end do
-      call write_attribute(tallies_group, "ids", id_array)
-      deallocate(id_array)
-
-      ! Write all tally information except results
-      TALLY_METADATA: do i = 1, n_tallies
-
-        ! Get pointer to tally
-        associate (tally => tallies(i) % obj)
-        tally_group = create_group(tallies_group, "tally " // &
-             trim(to_str(tally % id)))
-
-        ! Write the name for this tally
-        call write_dataset(tally_group, "name", tally % name)
-
-        select case(tally % estimator)
-        case (ESTIMATOR_ANALOG)
-          call write_dataset(tally_group, "estimator", "analog")
-        case (ESTIMATOR_TRACKLENGTH)
-          call write_dataset(tally_group, "estimator", "tracklength")
-        case (ESTIMATOR_COLLISION)
-          call write_dataset(tally_group, "estimator", "collision")
-        end select
-        call write_dataset(tally_group, "n_realizations", &
-             tally % n_realizations)
-
-        call write_dataset(tally_group, "n_filters", size(tally % filter))
-        if (size(tally % filter) > 0) then
-          ! Write IDs of filters
-          allocate(id_array(size(tally % filter)))
-          do j = 1, size(tally % filter)
-            id_array(j) = filters(tally % filter(j)) % obj % id
-          end do
-          call write_dataset(tally_group, "filters", id_array)
-          deallocate(id_array)
-        end if
-
-        ! Set up nuclide bin array and then write
-        allocate(str_array(tally % n_nuclide_bins))
-        NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
-          if (tally % nuclide_bins(j) > 0) then
-            if (run_CE) then
-              i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
-              if (i_xs > 0) then
-                str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
-              else
-                str_array(j) = nuclides(tally % nuclide_bins(j)) % name
-              end if
-            else
-              call get_name_c(tally % nuclide_bins(j), len(temp_name), &
-                              temp_name)
-              i_xs = index(temp_name, '.')
-              if (i_xs > 0) then
-                str_array(j) = trim(temp_name(1 : i_xs-1))
-              else
-                str_array(j) = trim(temp_name)
-              end if
-            end if
-          else
-            str_array(j) = 'total'
-          end if
-        end do NUCLIDE_LOOP
-        call write_dataset(tally_group, "nuclides", str_array)
-        deallocate(str_array)
-
-        ! Write derivative information.
-        if (tally % deriv /= NONE) then
-          call write_dataset(tally_group, "derivative", &
-               tally_derivs(tally % deriv) % id)
-        end if
-
-        ! Write scores.
-        call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
-        allocate(str_array(size(tally % score_bins)))
-        do j = 1, size(tally % score_bins)
-          str_array(j) = reaction_name(tally % score_bins(j))
-        end do
-        call write_dataset(tally_group, "score_bins", str_array)
-
-        deallocate(str_array)
-
-        call close_group(tally_group)
-        end associate
-      end do TALLY_METADATA
-    end if
-
     if (reduce_tallies) then
       ! Write global tallies
       call write_dataset(file_id, "global_tallies", global_tallies)
 
       ! Write tallies
-      if (active_tallies % size() > 0) then
+      if (active_tallies_size() > 0) then
         ! Indicate that tallies are on
         call write_attribute(file_id, "tallies_present", 1)
 
@@ -214,7 +49,7 @@ contains
           associate (tally => tallies(i) % obj)
             ! Write sum and sum_sq for each bin
             tally_group = open_group(tallies_group, "tally " &
-                 // to_str(tally % id))
+                 // to_str(tally % id()))
             call tally % write_results_hdf5(tally_group)
             call close_group(tally_group)
           end associate
@@ -254,7 +89,7 @@ contains
         associate (t => tallies(i) % obj)
           ! Read sum, sum_sq, and N for each bin
           tally_group = open_group(tallies_group, "tally " // &
-                trim(to_str(t % id)))
+                trim(to_str(t % id())))
           call t % read_results_hdf5(tally_group)
           call read_dataset(t % n_realizations, tally_group, &
                 "n_realizations")

src/state_point.cpp

@@ -17,9 +17,12 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/mesh.h"
 #include "openmc/message_passing.h"
+#include "openmc/mgxs_interface.h"
+#include "openmc/nuclide.h"
 #include "openmc/output.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+#include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/timer.h"
@@ -101,13 +104,135 @@ openmc_statepoint_write(const char* filename, bool* write_source)
     write_attribute(file_id, "source_present", write_source_);
 
     // Write out information for eigenvalue run
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) write_eigenvalue_hdf5(file_id);
+    if (settings::run_mode == RUN_MODE_EIGENVALUE)
+      write_eigenvalue_hdf5(file_id);
 
-    // Write tally-related data
     hid_t tallies_group = create_group(file_id, "tallies");
+
+    // Write meshes
     meshes_to_hdf5(tallies_group);
-    statepoint_write_f(file_id);
+
+    // Write information for derivatives
+    if (!model::tally_derivs.empty()) {
+      hid_t derivs_group = create_group(tallies_group, "derivatives");
+      for (const auto& deriv : model::tally_derivs) {
+        hid_t deriv_group = create_group(derivs_group,
+          "derivative " + std::to_string(deriv.id));
+        write_dataset(deriv_group, "material", deriv.diff_material);
+        if (deriv.variable == DIFF_DENSITY) {
+          write_dataset(deriv_group, "independent variable", "density");
+        } else if (deriv.variable == DIFF_NUCLIDE_DENSITY) {
+          write_dataset(deriv_group, "independent variable", "nuclide_density");
+          //TODO: off-by-one
+          write_dataset(deriv_group, "nuclide",
+            data::nuclides[deriv.diff_nuclide-1]->name_);
+        } else if (deriv.variable == DIFF_TEMPERATURE) {
+          write_dataset(deriv_group, "independent variable", "temperature");
+        } else {
+          fatal_error("Independent variable for derivative "
+            + std::to_string(deriv.id) + " not defined in state_point.cpp");
+        }
+        close_group(deriv_group);
+      }
+      close_group(derivs_group);
+    }
+
+    // Write information for filters
+    hid_t filters_group = create_group(tallies_group, "filters");
+    write_attribute(filters_group, "n_filters", model::tally_filters.size());
+    if (!model::tally_filters.empty()) {
+      // Write filter IDs
+      std::vector<int32_t> filter_ids;
+      filter_ids.reserve(model::tally_filters.size());
+      for (const auto& filt : model::tally_filters)
+        filter_ids.push_back(filt->id_);
+      write_attribute(filters_group, "ids", filter_ids);
+
+      // Write info for each filter
+      for (const auto& filt : model::tally_filters) {
+        hid_t filter_group = create_group(filters_group,
+          "filter " + std::to_string(filt->id_));
+        filt->to_statepoint(filter_group);
+        close_group(filter_group);
+      }
+    }
+    close_group(filters_group);
+
+    // Write information for tallies
+    write_attribute(tallies_group, "n_tallies", model::tallies.size());
+    if (!model::tallies.empty()) {
+      // Write tally IDs
+      std::vector<int32_t> tally_ids;
+      tally_ids.reserve(model::tallies.size());
+      for (const auto& tally : model::tallies)
+        tally_ids.push_back(tally->id_);
+      write_attribute(tallies_group, "ids", tally_ids);
+
+      // Write all tally information except results
+      //TODO: use these two lines when openmc_get_n_realizations isn't needed
+      //for (const auto& tally_ptr : model::tallies) {
+      //  const auto& tally {*tally_ptr};
+      for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
+        const auto& tally {*model::tallies[i_tally]};
+        hid_t tally_group = create_group(tallies_group,
+          "tally " + std::to_string(tally.id_));
+
+        write_dataset(tally_group, "name",  tally.name_);
+
+        if (tally.estimator_ == ESTIMATOR_ANALOG) {
+          write_dataset(tally_group, "estimator", "analog");
+        } else if (tally.estimator_ == ESTIMATOR_TRACKLENGTH) {
+          write_dataset(tally_group, "estimator", "tracklength");
+        } else if (tally.estimator_ == ESTIMATOR_COLLISION) {
+          write_dataset(tally_group, "estimator", "collision");
+        }
+
+        // TODO: get this directly from the tally object when it's moved to C++
+        int32_t n_realizations;
+        auto err = openmc_tally_get_n_realizations(i_tally+1, &n_realizations);
+        write_dataset(tally_group, "n_realizations", n_realizations);
+
+        // Write the ID of each filter attached to this tally
+        write_dataset(tally_group, "n_filters", tally.filters().size());
+        if (!tally.filters().empty()) {
+          std::vector<int32_t> filter_ids;
+          filter_ids.reserve(tally.filters().size());
+          for (auto i_filt : tally.filters())
+            filter_ids.push_back(model::tally_filters[i_filt]->id_);
+          write_dataset(tally_group, "filters", filter_ids);
+        }
+
+        // Write the nuclides this tally scores
+        std::vector<std::string> nuclides;
+        for (auto i_nuclide : tally.nuclides_) {
+          if (i_nuclide == -1) {
+            nuclides.push_back("total");
+          } else {
+            if (settings::run_CE) {
+              nuclides.push_back(data::nuclides[i_nuclide]->name_);
+            } else {
+              nuclides.push_back(data::nuclides_MG[i_nuclide].name);
+            }
+          }
+        }
+        write_dataset(tally_group, "nuclides", nuclides);
+
+        if (tally.deriv_ != C_NONE) write_dataset(tally_group, "derivative",
+          model::tally_derivs[tally.deriv_].id);
+
+        // Write the tally score bins
+        std::vector<std::string> scores;
+        for (auto sc : tally.scores_) scores.push_back(reaction_name(sc));
+        write_dataset(tally_group, "n_score_bins", scores.size());
+        write_dataset(tally_group, "score_bins", scores);
+
+        close_group(tally_group);
+      }
+
+    }
+
     close_group(tallies_group);
+    statepoint_write_f(file_id);
   }
 
   // Check for the no-tally-reduction method
@@ -529,10 +654,11 @@ void write_tally_results_nr(hid_t file_id)
     write_dataset(file_id, "global_tallies", gt);
   }
 
-  for (int i = 1; i <= n_tallies; ++i) {
+  for (int i = 0; i < n_tallies; ++i) {
     // Skip any tallies that are not active
     bool active;
-    openmc_tally_get_active(i, &active);
+    // TODO: off-by-one
+    openmc_tally_get_active(i+1, &active);
     if (!active) continue;
 
     if (mpi::master && !object_exists(file_id, "tallies_present")) {
@@ -550,7 +676,8 @@ void write_tally_results_nr(hid_t file_id)
     if (mpi::master) {
       // Open group for tally
       int id;
-      openmc_tally_get_id(i, &id);
+      // TODO: off-by-one
+      openmc_tally_get_id(i+1, &id);
       std::string groupname {"tally " + std::to_string(id)};
       hid_t tally_group = open_group(tallies_group, groupname.c_str());
 

src/tallies/derivative.cpp

@@ -0,0 +1,720 @@
+#include "openmc/tallies/derivative.h"
+
+#include "openmc/error.h"
+#include "openmc/material.h"
+#include "openmc/nuclide.h"
+#include "openmc/settings.h"
+#include "openmc/tallies/tally.h"
+#include "openmc/xml_interface.h"
+
+#include <sstream>
+
+template class std::vector<openmc::TallyDerivative>;
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+namespace model {
+  std::vector<TallyDerivative> tally_derivs;
+  std::unordered_map<int, int> tally_deriv_map;
+}
+
+//==============================================================================
+// TallyDerivative implementation
+//==============================================================================
+
+TallyDerivative::TallyDerivative(pugi::xml_node node)
+{
+  if (check_for_node(node, "id")) {
+    id = std::stoi(get_node_value(node, "id"));
+  } else {
+    fatal_error("Must specify an ID for <derivative> elements in the tally "
+                "XML file");
+  }
+
+  if (id <= 0)
+    fatal_error("<derivative> IDs must be an integer greater than zero");
+
+  std::string variable_str = get_node_value(node, "variable");
+
+  if (variable_str == "density") {
+    variable = DIFF_DENSITY;
+
+  } else if (variable_str == "nuclide_density") {
+    variable = DIFF_NUCLIDE_DENSITY;
+
+    std::string nuclide_name = get_node_value(node, "nuclide");
+    bool found = false;
+    for (auto i = 0; i < data::nuclides.size(); ++i) {
+      if (data::nuclides[i]->name_ == nuclide_name) {
+        found = true;
+        //TODO: off-by-one
+        diff_nuclide = i + 1;
+      }
+    }
+    if (!found) {
+      std::stringstream out;
+      out << "Could not find the nuclide \"" << nuclide_name
+          << "\" specified in derivative " << id << " in any material.";
+      fatal_error(out);
+    }
+
+  } else if (variable_str == "temperature") {
+    variable = DIFF_TEMPERATURE;
+
+  } else {
+    std::stringstream out;
+    out << "Unrecognized variable \"" << variable_str
+        << "\" on derivative " << id;
+    fatal_error(out);
+  }
+
+  diff_material = std::stoi(get_node_value(node, "material"));
+}
+
+//==============================================================================
+// Non-method functions
+//==============================================================================
+
+extern "C" void
+read_tally_derivatives(pugi::xml_node* node)
+{
+  // Populate the derivatives array.  This must be done in parallel because
+  // the derivatives are threadprivate.
+  #pragma omp parallel
+  {
+    for (auto deriv_node : node->children("derivative"))
+      model::tally_derivs.emplace_back(deriv_node);
+  }
+
+  // Fill the derivative map.
+  for (auto i = 0; i < model::tally_derivs.size(); ++i) {
+    auto id = model::tally_derivs[i].id;
+    auto search = model::tally_deriv_map.find(id);
+    if (search == model::tally_deriv_map.end()) {
+      model::tally_deriv_map[id] = i;
+    } else {
+      fatal_error("Two or more derivatives use the same unique ID: "
+                  + std::to_string(id));
+    }
+  }
+
+  // Make sure derivatives were not requested for an MG run.
+  if (!settings::run_CE && !model::tally_derivs.empty())
+    fatal_error("Differential tallies not supported in multi-group mode");
+}
+
+void
+apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
+  double atom_density, int score_bin, double& score)
+{
+  const Tally& tally {*model::tallies[i_tally]};
+
+  if (score == 0.0) return;
+
+  // If our score was previously c then the new score is
+  // c * (1/f * d_f/d_p + 1/c * d_c/d_p)
+  // where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the
+  // perturbated variable.
+
+  const auto& deriv {model::tally_derivs[tally.deriv_]};
+  auto flux_deriv = deriv.flux_deriv;
+
+  // Handle special cases where we know that d_c/d_p must be zero.
+  if (score_bin == SCORE_FLUX) {
+    score *= flux_deriv;
+    return;
+  } else if (p->material == MATERIAL_VOID) {
+    score *= flux_deriv;
+    return;
+  }
+  //TODO: off-by-one
+  const Material& material {*model::materials[p->material-1]};
+  if (material.id_ != deriv.diff_material) {
+    score *= flux_deriv;
+    return;
+  }
+
+  switch (deriv.variable) {
+
+  //============================================================================
+  // Density derivative:
+  // c = Sigma_MT
+  // c = sigma_MT * N
+  // c = sigma_MT * rho * const
+  // d_c / d_rho = sigma_MT * const
+  // (1 / c) * (d_c / d_rho) = 1 / rho
+
+  case DIFF_DENSITY:
+    switch (tally.estimator_) {
+
+    case ESTIMATOR_ANALOG:
+    case ESTIMATOR_COLLISION:
+      switch (score_bin) {
+
+      case SCORE_TOTAL:
+      case SCORE_SCATTER:
+      case SCORE_ABSORPTION:
+      case SCORE_FISSION:
+      case SCORE_NU_FISSION:
+        score *= flux_deriv + 1. / material.density_gpcc_;
+        break;
+
+      default:
+        fatal_error("Tally derivative not defined for a score on tally "
+          + std::to_string(tally.id_));
+      }
+      break;
+
+    default:
+      fatal_error("Differential tallies are only implemented for analog and "
+        "collision estimators.");
+    }
+    break;
+
+  //============================================================================
+  // Nuclide density derivative:
+  // If we are scoring a reaction rate for a single nuclide then
+  // c = Sigma_MT_i
+  // c = sigma_MT_i * N_i
+  // d_c / d_N_i = sigma_MT_i
+  // (1 / c) * (d_c / d_N_i) = 1 / N_i
+  // If the score is for the total material (i_nuclide = -1)
+  // c = Sum_i(Sigma_MT_i)
+  // d_c / d_N_i = sigma_MT_i
+  // (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT
+  // where i is the perturbed nuclide.
+
+  case DIFF_NUCLIDE_DENSITY:
+    //TODO: off-by-one throughout on diff_nuclide
+    switch (tally.estimator_) {
+
+    case ESTIMATOR_ANALOG:
+      if (p->event_nuclide != deriv.diff_nuclide) {
+        score *= flux_deriv;
+        return;
+      }
+
+      switch (score_bin) {
+
+      case SCORE_TOTAL:
+      case SCORE_SCATTER:
+      case SCORE_ABSORPTION:
+      case SCORE_FISSION:
+      case SCORE_NU_FISSION:
+        {
+          // Find the index of the perturbed nuclide.
+          int i;
+          for (i = 0; i < material.nuclide_.size(); ++i)
+            if (material.nuclide_[i] == deriv.diff_nuclide - 1) break;
+          score *= flux_deriv + 1. / material.atom_density_(i);
+        }
+        break;
+
+      default:
+        fatal_error("Tally derivative not defined for a score on tally "
+          + std::to_string(tally.id_));
+      }
+      break;
+
+    case ESTIMATOR_COLLISION:
+      switch (score_bin) {
+
+      case SCORE_TOTAL:
+        if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
+          score *= flux_deriv
+            + simulation::micro_xs[deriv.diff_nuclide-1].total
+            / simulation::material_xs.total;
+        } else if (i_nuclide == deriv.diff_nuclide-1
+                   && simulation::micro_xs[i_nuclide].total) {
+          score *= flux_deriv + 1. / atom_density;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_SCATTER:
+        if (i_nuclide == -1 && (simulation::material_xs.total
+                                - simulation::material_xs.absorption) > 0.0) {
+          score *= flux_deriv
+            + (simulation::micro_xs[deriv.diff_nuclide-1].total
+            - simulation::micro_xs[deriv.diff_nuclide-1].absorption)
+            / (simulation::material_xs.total
+            - simulation::material_xs.absorption);
+        } else if (i_nuclide == deriv.diff_nuclide-1) {
+          score *= flux_deriv + 1. / atom_density;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_ABSORPTION:
+        if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
+          score *= flux_deriv
+            + simulation::micro_xs[deriv.diff_nuclide-1].absorption
+            / simulation::material_xs.absorption;
+        } else if (i_nuclide == deriv.diff_nuclide-1
+                   && simulation::micro_xs[i_nuclide].absorption) {
+          score *= flux_deriv + 1. / atom_density;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_FISSION:
+        if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
+          score *= flux_deriv
+            + simulation::micro_xs[deriv.diff_nuclide-1].fission
+            / simulation::material_xs.fission;
+        } else if (i_nuclide == deriv.diff_nuclide-1
+                   && simulation::micro_xs[i_nuclide].fission) {
+          score *= flux_deriv + 1. / atom_density;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_NU_FISSION:
+        if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
+          score *= flux_deriv
+            + simulation::micro_xs[deriv.diff_nuclide-1].nu_fission
+            / simulation::material_xs.nu_fission;
+        } else if (i_nuclide == deriv.diff_nuclide-1
+                   && simulation::micro_xs[i_nuclide].nu_fission) {
+          score *= flux_deriv + 1. / atom_density;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      default:
+        fatal_error("Tally derivative not defined for a score on tally "
+          + std::to_string(tally.id_));
+      }
+      break;
+
+    default:
+      fatal_error("Differential tallies are only implemented for analog and "
+        "collision estimators.");
+    }
+    break;
+
+  //============================================================================
+  // Temperature derivative:
+  // If we are scoring a reaction rate for a single nuclide then
+  // c = Sigma_MT_i
+  // c = sigma_MT_i * N_i
+  // d_c / d_T = (d_sigma_Mt_i / d_T) * N_i
+  // (1 / c) * (d_c / d_T) = (d_sigma_MT_i / d_T) / sigma_MT_i
+  // If the score is for the total material (i_nuclide = -1)
+  // (1 / c) * (d_c / d_T) = Sum_i((d_sigma_MT_i / d_T) * N_i) / Sigma_MT_i
+  // where i is the perturbed nuclide.  The d_sigma_MT_i / d_T term is
+  // computed by multipole_deriv_eval.  It only works for the resolved
+  // resonance range and requires multipole data.
+
+  case DIFF_TEMPERATURE:
+    switch (tally.estimator_) {
+
+    case ESTIMATOR_ANALOG:
+      {
+        // Find the index of the event nuclide.
+        int i;
+        for (i = 0; i < material.nuclide_.size(); ++i)
+          if (material.nuclide_[i] == p->event_nuclide-1) break;
+
+        const auto& nuc {*data::nuclides[p->event_nuclide-1]};
+        if (!multipole_in_range(&nuc, p->last_E)) {
+          score *= flux_deriv;
+          break;
+        }
+
+        switch (score_bin) {
+
+        case SCORE_TOTAL:
+          if (simulation::micro_xs[p->event_nuclide-1].total) {
+            double dsig_s, dsig_a, dsig_f;
+            std::tie(dsig_s, dsig_a, dsig_f)
+              = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+            score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i)
+              / simulation::material_xs.total;
+          } else {
+            score *= flux_deriv;
+          }
+          break;
+
+        case SCORE_SCATTER:
+          if (simulation::micro_xs[p->event_nuclide-1].total
+              - simulation::micro_xs[p->event_nuclide-1].absorption) {
+            double dsig_s, dsig_a, dsig_f;
+            std::tie(dsig_s, dsig_a, dsig_f)
+              = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+            score *= flux_deriv + dsig_s * material.atom_density_(i)
+              / (simulation::material_xs.total
+              - simulation::material_xs.absorption);
+          } else {
+            score *= flux_deriv;
+          }
+          break;
+
+        case SCORE_ABSORPTION:
+          if (simulation::micro_xs[p->event_nuclide-1].absorption) {
+            double dsig_s, dsig_a, dsig_f;
+            std::tie(dsig_s, dsig_a, dsig_f)
+              = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+            score *= flux_deriv + dsig_a * material.atom_density_(i)
+              / simulation::material_xs.absorption;
+          } else {
+            score *= flux_deriv;
+          }
+          break;
+
+        case SCORE_FISSION:
+          if (simulation::micro_xs[p->event_nuclide-1].fission) {
+            double dsig_s, dsig_a, dsig_f;
+            std::tie(dsig_s, dsig_a, dsig_f)
+              = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+            score *= flux_deriv + dsig_f * material.atom_density_(i)
+              / simulation::material_xs.fission;
+          } else {
+            score *= flux_deriv;
+          }
+          break;
+
+        case SCORE_NU_FISSION:
+          if (simulation::micro_xs[p->event_nuclide-1].fission) {
+            double nu = simulation::micro_xs[p->event_nuclide-1].nu_fission
+              / simulation::micro_xs[p->event_nuclide-1].fission;
+            double dsig_s, dsig_a, dsig_f;
+            std::tie(dsig_s, dsig_a, dsig_f)
+              = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+            score *= flux_deriv + nu * dsig_f * material.atom_density_(i)
+              / simulation::material_xs.nu_fission;
+          } else {
+            score *= flux_deriv;
+          }
+          break;
+
+        default:
+          fatal_error("Tally derivative not defined for a score on tally "
+            + std::to_string(tally.id_));
+        }
+      }
+      break;
+
+    case ESTIMATOR_COLLISION:
+      if (i_nuclide != -1) {
+        const auto& nuc {data::nuclides[i_nuclide]};
+        if (!multipole_in_range(nuc.get(), p->last_E)) {
+          score *= flux_deriv;
+          return;
+        }
+      }
+
+      switch (score_bin) {
+
+      case SCORE_TOTAL:
+        if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
+          double cum_dsig = 0;
+          for (auto i = 0; i < material.nuclide_.size(); ++i) {
+            auto i_nuc = material.nuclide_[i];
+            const auto& nuc {*data::nuclides[i_nuc]};
+            if (multipole_in_range(&nuc, p->last_E)
+                && simulation::micro_xs[i_nuc].total) {
+              double dsig_s, dsig_a, dsig_f;
+              std::tie(dsig_s, dsig_a, dsig_f)
+                = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+              cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i);
+            }
+          }
+          score *= flux_deriv + cum_dsig / simulation::material_xs.total;
+        } else if (simulation::micro_xs[i_nuclide].total) {
+          const auto& nuc {*data::nuclides[i_nuclide]};
+          double dsig_s, dsig_a, dsig_f;
+          std::tie(dsig_s, dsig_a, dsig_f)
+            = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+          score *= flux_deriv
+            + (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_SCATTER:
+        if (i_nuclide == -1 && (simulation::material_xs.total
+            - simulation::material_xs.absorption)) {
+          double cum_dsig = 0;
+          for (auto i = 0; i < material.nuclide_.size(); ++i) {
+            auto i_nuc = material.nuclide_[i];
+            const auto& nuc {*data::nuclides[i_nuc]};
+            if (multipole_in_range(&nuc, p->last_E)
+                && (simulation::micro_xs[i_nuc].total
+                - simulation::micro_xs[i_nuc].absorption)) {
+              double dsig_s, dsig_a, dsig_f;
+              std::tie(dsig_s, dsig_a, dsig_f)
+                = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+              cum_dsig += dsig_s * material.atom_density_(i);
+            }
+          }
+          score *= flux_deriv + cum_dsig / (simulation::material_xs.total
+            - simulation::material_xs.absorption);
+        } else if (simulation::micro_xs[i_nuclide].total
+                   - simulation::micro_xs[i_nuclide].absorption) {
+          const auto& nuc {*data::nuclides[i_nuclide]};
+          double dsig_s, dsig_a, dsig_f;
+          std::tie(dsig_s, dsig_a, dsig_f)
+            = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+          score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total
+            - simulation::micro_xs[i_nuclide].absorption);
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_ABSORPTION:
+        if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
+          double cum_dsig = 0;
+          for (auto i = 0; i < material.nuclide_.size(); ++i) {
+            auto i_nuc = material.nuclide_[i];
+            const auto& nuc {*data::nuclides[i_nuc]};
+            if (multipole_in_range(&nuc, p->last_E)
+                && simulation::micro_xs[i_nuc].absorption) {
+              double dsig_s, dsig_a, dsig_f;
+              std::tie(dsig_s, dsig_a, dsig_f)
+                = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+              cum_dsig += dsig_a * material.atom_density_(i);
+            }
+          }
+          score *= flux_deriv + cum_dsig / simulation::material_xs.absorption;
+        } else if (simulation::micro_xs[i_nuclide].absorption) {
+          const auto& nuc {*data::nuclides[i_nuclide]};
+          double dsig_s, dsig_a, dsig_f;
+          std::tie(dsig_s, dsig_a, dsig_f)
+            = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+          score *= flux_deriv
+            + dsig_a / simulation::micro_xs[i_nuclide].absorption;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_FISSION:
+        if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
+          double cum_dsig = 0;
+          for (auto i = 0; i < material.nuclide_.size(); ++i) {
+            auto i_nuc = material.nuclide_[i];
+            const auto& nuc {*data::nuclides[i_nuc]};
+            if (multipole_in_range(&nuc, p->last_E)
+                && simulation::micro_xs[i_nuc].fission) {
+              double dsig_s, dsig_a, dsig_f;
+              std::tie(dsig_s, dsig_a, dsig_f)
+                = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+              cum_dsig += dsig_f * material.atom_density_(i);
+            }
+          }
+          score *= flux_deriv + cum_dsig / simulation::material_xs.fission;
+        } else if (simulation::micro_xs[i_nuclide].fission) {
+          const auto& nuc {*data::nuclides[i_nuclide]};
+          double dsig_s, dsig_a, dsig_f;
+          std::tie(dsig_s, dsig_a, dsig_f)
+            = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+          score *= flux_deriv
+            + dsig_f / simulation::micro_xs[i_nuclide].fission;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      case SCORE_NU_FISSION:
+        if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
+          double cum_dsig = 0;
+          for (auto i = 0; i < material.nuclide_.size(); ++i) {
+            auto i_nuc = material.nuclide_[i];
+            const auto& nuc {*data::nuclides[i_nuc]};
+            if (multipole_in_range(&nuc, p->last_E)
+                && simulation::micro_xs[i_nuc].fission) {
+              double nu = simulation::micro_xs[i_nuc].nu_fission
+                / simulation::micro_xs[i_nuc].fission;
+              double dsig_s, dsig_a, dsig_f;
+              std::tie(dsig_s, dsig_a, dsig_f)
+                = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+              cum_dsig += nu * dsig_f * material.atom_density_(i);
+            }
+          }
+          score *= flux_deriv + cum_dsig / simulation::material_xs.nu_fission;
+        } else if (simulation::micro_xs[i_nuclide].fission) {
+          const auto& nuc {*data::nuclides[i_nuclide]};
+          double dsig_s, dsig_a, dsig_f;
+          std::tie(dsig_s, dsig_a, dsig_f)
+            = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+          score *= flux_deriv
+            + dsig_f / simulation::micro_xs[i_nuclide].fission;
+        } else {
+          score *= flux_deriv;
+        }
+        break;
+
+      default:
+        break;
+      }
+      break;
+
+    default:
+      fatal_error("Differential tallies are only implemented for analog and "
+        "collision estimators.");
+    }
+    break;
+  }
+}
+
+//! Adjust diff tally flux derivatives for a particle tracking event.
+
+extern "C" void
+score_track_derivative(const Particle* p, double distance)
+{
+  // A void material cannot be perturbed so it will not affect flux derivatives.
+  if (p->material == MATERIAL_VOID) return;
+  //TODO: off-by-one
+  const Material& material {*model::materials[p->material-1]};
+
+  for (auto& deriv : model::tally_derivs) {
+    if (deriv.diff_material != material.id_) continue;
+
+    switch (deriv.variable) {
+
+    case DIFF_DENSITY:
+      // phi is proportional to e^(-Sigma_tot * dist)
+      // (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
+      // (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
+      deriv.flux_deriv -= distance * simulation::material_xs.total
+        / material.density_gpcc_;
+      break;
+
+    case DIFF_NUCLIDE_DENSITY:
+      // phi is proportional to e^(-Sigma_tot * dist)
+      // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
+      // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
+      //TODO: off-by-one
+      deriv.flux_deriv -= distance
+        * simulation::micro_xs[deriv.diff_nuclide-1].total;
+      break;
+
+    case DIFF_TEMPERATURE:
+      for (auto i = 0; i < material.nuclide_.size(); ++i) {
+        const auto& nuc {*data::nuclides[material.nuclide_[i]]};
+        if (multipole_in_range(&nuc, p->last_E)) {
+          // phi is proportional to e^(-Sigma_tot * dist)
+          // (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
+          // (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
+          double dsig_s, dsig_a, dsig_f;
+          std::tie(dsig_s, dsig_a, dsig_f)
+            = nuc.multipole_->evaluate_deriv(p->E, p->sqrtkT);
+          deriv.flux_deriv -= distance * (dsig_s + dsig_a)
+            * material.atom_density_(i);
+        }
+      }
+      break;
+    }
+  }
+}
+
+//! Adjust diff tally flux derivatives for a particle scattering event.
+//
+//! Note that this subroutine will be called after absorption events in
+//! addition to scattering events, but any flux derivatives scored after an
+//! absorption will never be tallied.  The paricle will be killed before any
+//! further tallies are scored.
+
+extern "C" void
+score_collision_derivative(const Particle* p)
+{
+  // A void material cannot be perturbed so it will not affect flux derivatives.
+  if (p->material == MATERIAL_VOID) return;
+  //TODO: off-by-one
+  const Material& material {*model::materials[p->material-1]};
+
+  for (auto& deriv : model::tally_derivs) {
+    if (deriv.diff_material != material.id_) continue;
+
+    switch (deriv.variable) {
+
+    case DIFF_DENSITY:
+      // phi is proportional to Sigma_s
+      // (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
+      // (1 / phi) * (d_phi / d_rho) = 1 / rho
+      deriv.flux_deriv += 1. / material.density_gpcc_;
+      break;
+
+    case DIFF_NUCLIDE_DENSITY:
+      //TODO: off-by-one throughout on diff_nuclide
+      if (p->event_nuclide != deriv.diff_nuclide) continue;
+      // Find the index in this material for the diff_nuclide.
+      int i;
+      for (i = 0; i < material.nuclide_.size(); ++i)
+        if (material.nuclide_[i] == deriv.diff_nuclide - 1) break;
+      // Make sure we found the nuclide.
+      if (material.nuclide_[i] != deriv.diff_nuclide - 1) {
+        std::stringstream err_msg;
+        err_msg << "Could not find nuclide "
+          << data::nuclides[deriv.diff_nuclide-1]->name_ << " in material "
+          << material.id_ << " for tally derivative " << deriv.id;
+        fatal_error(err_msg);
+      }
+      // phi is proportional to Sigma_s
+      // (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s
+      // (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s
+      // (1 / phi) * (d_phi / d_N) = 1 / N
+      deriv.flux_deriv += 1. / material.atom_density_(i);
+      break;
+
+    case DIFF_TEMPERATURE:
+      // Loop over the material's nuclides until we find the event nuclide.
+      for (auto i_nuc : material.nuclide_) {
+        const auto& nuc {*data::nuclides[i_nuc]};
+        //TODO: off-by-one
+        if (i_nuc == p->event_nuclide - 1
+            && multipole_in_range(&nuc, p->last_E)) {
+          // phi is proportional to Sigma_s
+          // (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
+          // (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
+          const auto& micro_xs {simulation::micro_xs[i_nuc]};
+          double dsig_s, dsig_a, dsig_f;
+          std::tie(dsig_s, dsig_a, dsig_f)
+            = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
+          deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
+          // Note that this is an approximation!  The real scattering cross
+          // section is
+          // Sigma_s(E'->E, uvw'->uvw) = Sigma_s(E') * P(E'->E, uvw'->uvw).
+          // We are assuming that d_P(E'->E, uvw'->uvw) / d_T = 0 and only
+          // computing d_S(E') / d_T.  Using this approximation in the vicinity
+          // of low-energy resonances causes errors (~2-5% for PWR pincell
+          // eigenvalue derivatives).
+        }
+      }
+      break;
+    }
+  }
+}
+
+//! Set the flux derivatives on differential tallies to zero.
+extern "C" void
+zero_flux_derivs()
+{
+  for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.;
+}
+
+//==============================================================================
+// Fortran interop
+//==============================================================================
+
+extern "C" int n_tally_derivs() {return model::tally_derivs.size();}
+
+extern "C" TallyDerivative*
+tally_deriv_c(int i)
+{
+  return &model::tally_derivs[i];
+}
+
+}// namespace openmc

src/tallies/filter.cpp

@@ -43,21 +43,6 @@ namespace model {
 std::vector<std::unique_ptr<Filter>> tally_filters;
 } // namespace model
 
-//==============================================================================
-// Non-member functions
-//==============================================================================
-
-void
-free_memory_tally_c()
-{
-  #pragma omp parallel
-  {
-    simulation::filter_matches.clear();
-  }
-
-  model::tally_filters.clear();
-}
-
 //==============================================================================
 // Fortran compatibility functions
 //==============================================================================
@@ -65,38 +50,6 @@ free_memory_tally_c()
 extern "C" {
   // filter_match_point moved to simulation.cpp
 
-  void
-  filter_match_bins_push_back(FilterMatch* match, int val)
-  {match->bins_.push_back(val);}
-
-  void
-  filter_match_weights_push_back(FilterMatch* match, double val)
-  {match->weights_.push_back(val);}
-
-  void
-  filter_match_bins_clear(FilterMatch* match)
-  {match->bins_.clear();}
-
-  void
-  filter_match_weights_clear(FilterMatch* match)
-  {match->weights_.clear();}
-
-  int
-  filter_match_bins_size(FilterMatch* match)
-  {return match->bins_.size();}
-
-  int
-  filter_match_bins_data(FilterMatch* match, int indx)
-  {return match->bins_[indx-1];}
-
-  double
-  filter_match_weights_data(FilterMatch* match, int indx)
-  {return match->weights_[indx-1];}
-
-  void
-  filter_match_bins_set_data(FilterMatch* match, int indx, int val)
-  {match->bins_[indx-1] = val;}
-
   Filter*
   allocate_filter(const char* type)
   {
@@ -151,6 +104,10 @@ extern "C" {
     return model::tally_filters.back().get();
   }
 
+  int32_t filter_get_id(Filter* filt) {return filt->id_;}
+
+  void filter_set_id(Filter* filt, int32_t id) {filt->id_ = id;}
+
   void filter_from_xml(Filter* filt, pugi::xml_node* node)
   {filt->from_xml(*node);}
 

src/tallies/filter_delayedgroup.cpp

@@ -48,11 +48,4 @@ DelayedGroupFilter::text_label(int bin) const
   return "Delayed Group " + std::to_string(groups_[bin-1]);
 }
 
-//==============================================================================
-// Fortran interoperability
-//==============================================================================
-
-extern "C" int delayedgroup_filter_groups(DelayedGroupFilter* filt, int i)
-{return filt->groups_[i-1];}
-
 } // namespace openmc

src/tallies/filter_energy.cpp

@@ -112,24 +112,6 @@ EnergyoutFilter::text_label(int bin) const
   return out.str();
 }
 
-//==============================================================================
-// Fortran interoperability
-//==============================================================================
-
-extern "C" bool energy_filter_matches_transport_groups(EnergyFilter* filt)
-{return filt->matches_transport_groups_;}
-
-extern "C" int
-energy_filter_search(EnergyFilter* filt, double val)
-{
-  if (val < filt->bins_.front() || val > filt->bins_.back()) {
-    return -1;
-  } else {
-    //TODO: off-by-one
-    return lower_bound_index(filt->bins_.begin(), filt->bins_.end(), val) + 1;
-  }
-}
-
 //==============================================================================
 // C-API functions
 //==============================================================================

src/tallies/tally.cpp

@@ -2,7 +2,24 @@
 
 #include "openmc/capi.h"
 #include "openmc/constants.h"
+#include "openmc/error.h"
 #include "openmc/message_passing.h"
+#include "openmc/mgxs_interface.h"
+#include "openmc/nuclide.h"
+#include "openmc/reaction_product.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+#include "openmc/tallies/derivative.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/filter_cell.h"
+#include "openmc/tallies/filter_cellfrom.h"
+#include "openmc/tallies/filter_delayedgroup.h"
+#include "openmc/tallies/filter_energy.h"
+#include "openmc/tallies/filter_legendre.h"
+#include "openmc/tallies/filter_mesh.h"
+#include "openmc/tallies/filter_meshsurface.h"
+#include "openmc/tallies/filter_surface.h"
+#include "openmc/xml_interface.h"
 
 #include "xtensor/xadapt.hpp"
 #include "xtensor/xbuilder.hpp" // for empty_like
@@ -10,6 +27,8 @@
 
 #include <array>
 #include <cstddef>
+#include <sstream>
+#include <string>
 
 namespace openmc {
 
@@ -17,11 +36,473 @@ namespace openmc {
 // Global variable definitions
 //==============================================================================
 
+namespace model {
+  std::vector<std::unique_ptr<Tally>> tallies;
+
+  std::vector<int> active_tallies;
+  std::vector<int> active_analog_tallies;
+  std::vector<int> active_tracklength_tallies;
+  std::vector<int> active_collision_tallies;
+  std::vector<int> active_meshsurf_tallies;
+  std::vector<int> active_surface_tallies;
+}
+
 double global_tally_absorption;
 double global_tally_collision;
 double global_tally_tracklength;
 double global_tally_leakage;
 
+int
+score_str_to_int(std::string score_str)
+{
+  if (score_str == "flux")
+    return SCORE_FLUX;
+
+  if (score_str == "total" || score_str == "(n,total)")
+    return SCORE_TOTAL;
+
+  if (score_str == "scatter")
+    return SCORE_SCATTER;
+
+  if (score_str == "nu-scatter")
+    return SCORE_NU_SCATTER;
+
+  if (score_str == "absorption")
+    return SCORE_ABSORPTION;
+
+  if (score_str == "fission" || score_str == "18")
+    return SCORE_FISSION;
+
+  if (score_str == "nu-fission")
+    return SCORE_NU_FISSION;
+
+  if (score_str == "decay-rate")
+    return SCORE_DECAY_RATE;
+
+  if (score_str == "delayed-nu-fission")
+    return SCORE_DELAYED_NU_FISSION;
+
+  if (score_str == "prompt-nu-fission")
+    return SCORE_PROMPT_NU_FISSION;
+
+  if (score_str == "kappa-fission")
+    return SCORE_KAPPA_FISSION;
+
+  if (score_str == "inverse-velocity")
+    return SCORE_INVERSE_VELOCITY;
+
+  if (score_str == "fission-q-prompt")
+    return SCORE_FISS_Q_PROMPT;
+
+  if (score_str == "fission-q-recoverable")
+    return SCORE_FISS_Q_RECOV;
+
+  if (score_str == "current")
+    return SCORE_CURRENT;
+
+  if (score_str == "events")
+    return SCORE_EVENTS;
+
+  if (score_str == "elastic" || score_str == "(n,elastic)")
+    return ELASTIC;
+
+  if (score_str == "n2n" || score_str == "(n,2n)")
+    return N_2N;
+
+  if (score_str == "n3n" || score_str == "(n,3n)")
+    return N_3N;
+
+  if (score_str == "n4n" || score_str == "(n,4n)")
+    return N_4N;
+
+  if (score_str == "(n,2nd)")
+    return N_2ND;
+  if (score_str == "(n,na)")
+    return N_2NA;
+  if (score_str == "(n,n3a)")
+    return N_N3A;
+  if (score_str == "(n,2na)")
+    return N_2NA;
+  if (score_str == "(n,3na)")
+    return N_3NA;
+  if (score_str == "(n,np)")
+    return N_NP;
+  if (score_str == "(n,n2a)")
+    return N_N2A;
+  if (score_str == "(n,2n2a)")
+    return N_2N2A;
+  if (score_str == "(n,nd)")
+    return N_ND;
+  if (score_str == "(n,nt)")
+    return N_NT;
+  if (score_str == "(n,nHe-3)")
+    return N_N3HE;
+  if (score_str == "(n,nd2a)")
+    return N_ND2A;
+  if (score_str == "(n,nt2a)")
+    return N_NT2A;
+  if (score_str == "(n,3nf)")
+    return N_3NF;
+  if (score_str == "(n,2np)")
+    return N_2NP;
+  if (score_str == "(n,3np)")
+    return N_3NP;
+  if (score_str == "(n,n2p)")
+    return N_N2P;
+  if (score_str == "(n,npa)")
+    return N_NPA;
+  if (score_str == "(n,n1)")
+    return N_N1;
+  if (score_str == "(n,nc)")
+    return N_NC;
+  if (score_str == "(n,gamma)")
+    return N_GAMMA;
+  if (score_str == "(n,p)")
+    return N_P;
+  if (score_str == "(n,d)")
+    return N_D;
+  if (score_str == "(n,t)")
+    return N_T;
+  if (score_str == "(n,3He)")
+    return N_3HE;
+  if (score_str == "(n,a)")
+    return N_A;
+  if (score_str == "(n,2a)")
+    return N_2A;
+  if (score_str == "(n,3a)")
+    return N_3A;
+  if (score_str == "(n,2p)")
+    return N_2P;
+  if (score_str == "(n,pa)")
+    return N_PA;
+  if (score_str == "(n,t2a)")
+    return N_T2A;
+  if (score_str == "(n,d2a)")
+    return N_D2A;
+  if (score_str == "(n,pd)")
+    return N_PD;
+  if (score_str == "(n,pt)")
+    return N_PT;
+  if (score_str == "(n,da)")
+    return N_DA;
+
+  // So far we have not identified this score string.  Check to see if it is a
+  // deprecated score.
+  if (score_str.rfind("scatter-", 0) == 0
+      || score_str.rfind("nu-scatter-", 0) == 0
+      || score_str.rfind("total-y", 0) == 0
+      || score_str.rfind("flux-y", 0) == 0)
+    fatal_error(score_str + " is no longer an available score");
+
+
+  // Assume the given string is a reaction MT number.  Make sure it's a natural
+  // number then return.
+  int MT;
+  try {
+    MT = std::stoi(score_str);
+  } catch (const std::invalid_argument& ex) {
+    throw std::invalid_argument("Invalid tally score \"" + score_str + "\"");
+  }
+  if (MT < 1)
+    throw std::invalid_argument("Invalid tally score \"" + score_str + "\"");
+  return MT;
+}
+
+//==============================================================================
+// Tally object implementation
+//==============================================================================
+
+void
+Tally::init_from_xml(pugi::xml_node node)
+{
+  if (check_for_node(node, "name")) name_ = get_node_value(node, "name");
+}
+
+void
+Tally::set_filters(const int32_t filter_indices[], int n)
+{
+  // Clear old data.
+  filters_.clear();
+  strides_.clear();
+
+  // Copy in the given filter indices.
+  filters_.assign(filter_indices, filter_indices + n);
+
+  for (int i = 0; i < n; ++i) {
+    auto i_filt = filters_[i];
+    if (i_filt < 0 || i_filt >= model::tally_filters.size())
+      throw std::out_of_range("Index in tally filter array out of bounds.");
+
+    // Keep track of indices for special filters.
+    const auto* filt = model::tally_filters[i_filt].get();
+    if (dynamic_cast<const EnergyoutFilter*>(filt)) {
+      energyout_filter_ = i;
+    } else if (dynamic_cast<const DelayedGroupFilter*>(filt)) {
+      delayedgroup_filter_ = i;
+    }
+  }
+
+  // Set the strides.  Filters are traversed in reverse so that the last filter
+  // has the shortest stride in memory and the first filter has the longest
+  // stride.
+  strides_.resize(n, 0);
+  int stride = 1;
+  for (int i = n-1; i >= 0; --i) {
+    strides_[i] = stride;
+    stride *= model::tally_filters[filters_[i]]->n_bins_;
+  }
+  n_filter_bins_ = stride;
+}
+
+void
+Tally::set_scores(pugi::xml_node node)
+{
+  if (!check_for_node(node, "scores"))
+    fatal_error("No scores specified on tally " + std::to_string(id_));
+
+  auto scores = get_node_array<std::string>(node, "scores");
+  set_scores(scores);
+}
+
+void
+Tally::set_scores(std::vector<std::string> scores)
+{
+  // Reset state and prepare for the new scores.
+  scores_.clear();
+  depletion_rx_ = false;
+  scores_.reserve(scores.size());
+
+  // Check for the presence of certain restrictive filters.
+  bool energyout_present = energyout_filter_ != C_NONE;
+  bool legendre_present = false;
+  bool cell_present = false;
+  bool cellfrom_present = false;
+  bool surface_present = false;
+  bool meshsurface_present = false;
+  for (auto i_filt : filters_) {
+    const auto* filt {model::tally_filters[i_filt].get()};
+    if (dynamic_cast<const LegendreFilter*>(filt)) {
+      legendre_present = true;
+    } else if (dynamic_cast<const CellFromFilter*>(filt)) {
+      cellfrom_present = true;
+    } else if (dynamic_cast<const CellFilter*>(filt)) {
+      cell_present = true;
+    } else if (dynamic_cast<const SurfaceFilter*>(filt)) {
+      surface_present = true;
+    } else if (dynamic_cast<const MeshSurfaceFilter*>(filt)) {
+      meshsurface_present = true;
+    }
+  }
+
+  // Iterate over the given scores.
+  for (auto score_str : scores) {
+    // Make sure a delayed group filter wasn't used with an incompatible score.
+    bool has_delayedgroup = delayedgroup_filter_ != C_NONE;
+    if (delayedgroup_filter_ != C_NONE) {
+      if (score_str != "delayed-nu-fission" && score_str != "decay-rate")
+        fatal_error("Cannot tally " + score_str + "with a delayedgroup filter");
+    }
+
+    auto score = score_str_to_int(score_str);
+
+    switch (score) {
+    case SCORE_FLUX:
+      if (!nuclides_.empty())
+        if (!(nuclides_.size() == 1 && nuclides_[0] == -1))
+          fatal_error("Cannot tally flux for an individual nuclide.");
+      if (energyout_present)
+        fatal_error("Cannot tally flux with an outgoing energy filter.");
+      break;
+
+    case SCORE_TOTAL:
+    case SCORE_ABSORPTION:
+    case SCORE_FISSION:
+      if (energyout_present)
+        fatal_error("Cannot tally " + score_str + " reaction rate with an "
+          "outgoing energy filter");
+      break;
+
+    case SCORE_SCATTER:
+      if (legendre_present)
+        estimator_ = ESTIMATOR_ANALOG;
+    case SCORE_NU_FISSION:
+    case SCORE_DELAYED_NU_FISSION:
+    case SCORE_PROMPT_NU_FISSION:
+      if (energyout_present)
+        estimator_ = ESTIMATOR_ANALOG;
+      break;
+
+    case SCORE_NU_SCATTER:
+      if (settings::run_CE) {
+        estimator_ = ESTIMATOR_ANALOG;
+      } else {
+        if (energyout_present || legendre_present)
+          estimator_ = ESTIMATOR_ANALOG;
+      }
+      break;
+
+    case N_2N:
+    case N_3N:
+    case N_4N:
+    case N_GAMMA:
+    case N_P:
+    case N_A:
+      depletion_rx_ = true;
+      break;
+
+    case SCORE_CURRENT:
+      // Check which type of current is desired: mesh or surface currents.
+      if (surface_present || cell_present || cellfrom_present) {
+        if (meshsurface_present)
+          fatal_error("Cannot tally mesh surface currents in the same tally as "
+            "normal surface currents");
+        type_ = TALLY_SURFACE;
+      } else if (meshsurface_present) {
+        type_ = TALLY_MESH_SURFACE;
+      } else {
+        fatal_error("Cannot tally currents without surface type filters");
+      }
+      break;
+    }
+
+    scores_.push_back(score);
+  }
+
+  // Make sure that no duplicate scores exist.
+  for (auto it1 = scores_.begin(); it1 != scores_.end(); ++it1) {
+    for (auto it2 = it1 + 1; it2 != scores_.end(); ++it2) {
+      if (*it1 == *it2)
+        fatal_error("Duplicate score of type \"" + reaction_name(*it1)
+          + "\" found in tally " + std::to_string(id_));
+    }
+  }
+
+  // Make sure all scores are compatible with multigroup mode.
+  if (!settings::run_CE) {
+    for (auto sc : scores_)
+      if (sc > 0)
+        fatal_error("Cannot tally " + reaction_name(sc) + " reaction rate "
+          "in multi-group mode");
+  }
+
+  // Make sure current scores are not mixed in with volumetric scores.
+  if (type_ == TALLY_SURFACE || type_ == TALLY_MESH_SURFACE) {
+    if (scores_.size() != 1)
+      fatal_error("Cannot tally other scores in the same tally as surface "
+        "currents");
+  }
+}
+
+void
+Tally::set_nuclides(pugi::xml_node node)
+{
+  nuclides_.clear();
+
+  // By default, we tally just the total material rates.
+  if (!check_for_node(node, "nuclides")) {
+    nuclides_.push_back(-1);
+    return;
+  }
+
+  if (get_node_value(node, "nuclides") == "all") {
+    // This tally should bin every nuclide in the problem.  It should also bin
+    // the total material rates.  To achieve this, set the nuclides_ vector to
+    // 0, 1, 2, ..., -1.
+    nuclides_.reserve(data::nuclides.size() + 1);
+    for (auto i = 0; i < data::nuclides.size(); ++i)
+      nuclides_.push_back(i);
+    nuclides_.push_back(-1);
+    all_nuclides_ = true;
+
+  } else {
+    // The user provided specifics nuclides.  Parse it as an array with either
+    // "total" or a nuclide name like "U-235" in each position.
+    auto words = get_node_array<std::string>(node, "nuclides");
+    for (auto word : words) {
+      if (word == "total") {
+        nuclides_.push_back(-1);
+      } else {
+        auto search = data::nuclide_map.find(word);
+        if (search == data::nuclide_map.end())
+          fatal_error("Could not find the nuclide " + word
+            + " specified in tally " + std::to_string(id_)
+            + " in any material");
+        nuclides_.push_back(search->second);
+      }
+    }
+  }
+}
+
+void
+Tally::init_triggers(pugi::xml_node node, int i_tally)
+{
+  for (auto trigger_node: node.children("trigger")) {
+    // Read the trigger type.
+    TriggerMetric metric;
+    if (check_for_node(trigger_node, "type")) {
+      auto type_str = get_node_value(trigger_node, "type");
+      if (type_str == "std_dev") {
+        metric = TriggerMetric::standard_deviation;
+      } else if (type_str == "variance") {
+        metric = TriggerMetric::variance;
+      } else if (type_str == "rel_err") {
+        metric = TriggerMetric::relative_error;
+      } else {
+        std::stringstream msg;
+        msg << "Unknown trigger type \"" << type_str << "\" in tally "  << id_;
+        fatal_error(msg);
+      }
+    } else {
+      std::stringstream msg;
+      msg << "Must specify trigger type for tally " << id_
+          << " in tally XML file";
+      fatal_error(msg);
+    }
+
+    // Read the trigger threshold.
+    double threshold;
+    if (check_for_node(trigger_node, "threshold")) {
+      threshold = std::stod(get_node_value(trigger_node, "threshold"));
+    } else {
+      std::stringstream msg;
+      msg << "Must specify trigger threshold for tally " << id_
+          << " in tally XML file";
+      fatal_error(msg);
+    }
+
+    // Read the trigger scores.
+    std::vector<std::string> trigger_scores;
+    if (check_for_node(trigger_node, "scores")) {
+      trigger_scores = get_node_array<std::string>(trigger_node, "scores");
+    } else {
+      trigger_scores.push_back("all");
+    }
+
+    // Parse the trigger scores and populate the triggers_ vector.
+    const auto& tally {*model::tallies[i_tally]};
+    for (auto score_str : trigger_scores) {
+      if (score_str == "all") {
+        triggers_.reserve(triggers_.size() + tally.scores_.size());
+        for (auto i_score = 0; i_score < tally.scores_.size(); ++i_score) {
+          triggers_.push_back({metric, threshold, i_score});
+        }
+      } else {
+        int i_score = 0;
+        for (; i_score < tally.scores_.size(); ++i_score) {
+          if (reaction_name(tally.scores_[i_score]) == score_str) break;
+        }
+        if (i_score == tally.scores_.size()) {
+          std::stringstream msg;
+          msg << "Could not find the score \"" << score_str << "\" in tally "
+              << id_ << " but it was listed in a trigger on that tally";
+          fatal_error(msg);
+        }
+        triggers_.push_back({metric, threshold, i_score});
+      }
+    }
+  }
+}
+
 //==============================================================================
 // Non-member functions
 //==============================================================================
@@ -44,7 +525,8 @@ adaptor_type<3> tally_results(int idx)
   // Get pointer to tally results
   double* results;
   std::array<std::size_t, 3> shape;
-  openmc_tally_results(idx, &results, shape.data());
+  // TODO: off-by-one
+  openmc_tally_results(idx+1, &results, shape.data());
 
   // Adapt array into xtensor with no ownership
   std::size_t size {shape[0] * shape[1] * shape[2]};
@@ -54,10 +536,11 @@ adaptor_type<3> tally_results(int idx)
 #ifdef OPENMC_MPI
 void reduce_tally_results()
 {
-  for (int i = 1; i <= n_tallies; ++i) {
+  for (int i = 0; i < n_tallies; ++i) {
     // Skip any tallies that are not active
     bool active;
-    openmc_tally_get_active(i, &active);
+    // TODO: off-by-one
+    openmc_tally_get_active(i+1, &active);
     if (!active) continue;
 
     // Get view of accumulated tally values
@@ -108,4 +591,335 @@ void reduce_tally_results()
 }
 #endif
 
+extern "C" void
+setup_active_tallies_c()
+{
+  model::active_tallies.clear();
+  model::active_analog_tallies.clear();
+  model::active_tracklength_tallies.clear();
+  model::active_collision_tallies.clear();
+  model::active_meshsurf_tallies.clear();
+  model::active_surface_tallies.clear();
+
+  for (auto i = 0; i < model::tallies.size(); ++i) {
+    const auto& tally {*model::tallies[i]};
+
+    if (tally.active_) {
+      model::active_tallies.push_back(i);
+      switch (tally.type_) {
+
+      case TALLY_VOLUME:
+        switch (tally.estimator_) {
+          case ESTIMATOR_ANALOG:
+            model::active_analog_tallies.push_back(i);
+            break;
+          case ESTIMATOR_TRACKLENGTH:
+            model::active_tracklength_tallies.push_back(i);
+            break;
+          case ESTIMATOR_COLLISION:
+            model::active_collision_tallies.push_back(i);
+        }
+        break;
+
+      case TALLY_MESH_SURFACE:
+        model::active_meshsurf_tallies.push_back(i);
+        break;
+
+      case TALLY_SURFACE:
+        model::active_surface_tallies.push_back(i);
+      }
+    }
+  }
+}
+
+extern "C" void
+free_memory_tally_c()
+{
+  #pragma omp parallel
+  {
+    model::tally_derivs.clear();
+  }
+
+  model::tally_filters.clear();
+
+  model::tallies.clear();
+
+  model::active_tallies.clear();
+  model::active_analog_tallies.clear();
+  model::active_tracklength_tallies.clear();
+  model::active_collision_tallies.clear();
+  model::active_meshsurf_tallies.clear();
+  model::active_surface_tallies.clear();
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern "C" int
+openmc_tally_get_type(int32_t index, int32_t* type)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  //TODO: off-by-one
+  *type = model::tallies[index-1]->type_;
+
+  return 0;
+}
+
+extern "C" int
+openmc_tally_set_type(int32_t index, const char* type)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  if (strcmp(type, "volume") == 0) {
+    model::tallies[index-1]->type_ = TALLY_VOLUME;
+  } else if (strcmp(type, "mesh-surface") == 0) {
+    model::tallies[index-1]->type_ = TALLY_MESH_SURFACE;
+  } else if (strcmp(type, "surface") == 0) {
+    model::tallies[index-1]->type_ = TALLY_SURFACE;
+  } else {
+    std::stringstream errmsg;
+    errmsg << "Unknown tally type: " << type;
+    set_errmsg(errmsg);
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+
+  return 0;
+}
+
+extern "C" int
+openmc_tally_get_active(int32_t index, bool* active)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  //TODO: off-by-one
+  *active = model::tallies[index-1]->active_;
+
+  return 0;
+}
+
+extern "C" int
+openmc_tally_set_active(int32_t index, bool active)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  //TODO: off-by-one
+  model::tallies[index-1]->active_ = active;
+
+  return 0;
+}
+
+extern "C" int
+openmc_tally_get_scores(int32_t index, int** scores, int* n)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  //TODO: off-by-one
+  *scores = model::tallies[index-1]->scores_.data();
+  *n = model::tallies[index-1]->scores_.size();
+  return 0;
+}
+
+extern "C" int
+openmc_tally_set_scores(int32_t index, int n, const char** scores)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  std::vector<std::string> scores_str(scores, scores+n);
+  try {
+    //TODO: off-by-one
+    model::tallies[index-1]->set_scores(scores_str);
+  } catch (const std::invalid_argument& ex) {
+    set_errmsg(ex.what());
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+
+  return 0;
+}
+
+extern "C" int
+openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n)
+{
+  // Make sure the index fits in the array bounds.
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  //TODO: off-by-one
+  *n = model::tallies[index-1]->nuclides_.size();
+  *nuclides = model::tallies[index-1]->nuclides_.data();
+
+  return 0;
+}
+
+extern "C" int
+openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides)
+{
+  // Make sure the index fits in the array bounds.
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  std::vector<std::string> words(nuclides, nuclides+n);
+  std::vector<int> nucs;
+  for (auto word : words){
+    if (word == "total") {
+      nucs.push_back(-1);
+    } else {
+      auto search = data::nuclide_map.find(word);
+      if (search == data::nuclide_map.end()) {
+        set_errmsg("Nuclide \"" + word  + "\" has not been loaded yet");
+        return OPENMC_E_DATA;
+      }
+      nucs.push_back(search->second);
+    }
+  }
+
+  //TODO: off-by-one
+  model::tallies[index-1]->nuclides_ = nucs;
+
+  return 0;
+}
+
+extern "C" int
+openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  //TODO: off-by-one
+  *indices = model::tallies[index-1]->filters().data();
+  *n = model::tallies[index-1]->filters().size();
+  return 0;
+}
+
+extern "C" int
+openmc_tally_set_filters(int32_t index, int n, const int32_t* indices)
+{
+  // Make sure the index fits in the array bounds.
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  // Set the filters.
+  try {
+    //TODO: off-by-one
+    model::tallies[index-1]->set_filters(indices, n);
+  } catch (const std::out_of_range& ex) {
+    set_errmsg(ex.what());
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  return 0;
+}
+
+//==============================================================================
+// Fortran compatibility functions
+//==============================================================================
+
+extern "C" {
+  Tally* tally_pointer(int indx) {return model::tallies[indx].get();}
+
+  void
+  extend_tallies_c(int n)
+  {
+    for (int i = 0; i < n; ++i)
+      model::tallies.push_back(std::make_unique<Tally>());
+  }
+
+  int active_tallies_data(int i)
+  {return model::active_tallies[i-1];}
+
+  int active_tallies_size()
+  {return model::active_tallies.size();}
+
+  int active_analog_tallies_size()
+  {return model::active_analog_tallies.size();}
+
+  int active_tracklength_tallies_size()
+  {return model::active_tracklength_tallies.size();}
+
+  int active_collision_tallies_size()
+  {return model::active_collision_tallies.size();}
+
+  int active_meshsurf_tallies_size()
+  {return model::active_meshsurf_tallies.size();}
+
+  int active_surface_tallies_size()
+  {return model::active_surface_tallies.size();}
+
+  void tally_init_from_xml(Tally* tally, pugi::xml_node* node)
+  {tally->init_from_xml(*node);}
+
+  int tally_get_id_c(Tally* tally) {return tally->id_;}
+
+  void tally_set_id_c(Tally* tally, int id) {tally->id_ = id;}
+
+  int tally_get_type_c(Tally* tally) {return tally->type_;}
+
+  void tally_set_type_c(Tally* tally, int type) {tally->type_ = type;}
+
+  int tally_get_estimator_c(Tally* tally) {return tally->estimator_;}
+
+  void tally_set_estimator_c(Tally* tally, int e) {tally->estimator_ = e;}
+
+  bool tally_get_depletion_rx_c(Tally* tally) {return tally->depletion_rx_;}
+
+  int tally_get_n_scores_c(Tally* tally) {return tally->scores_.size();}
+
+  int tally_get_score_c(Tally* tally, int i) {return tally->scores_[i];}
+
+  void tally_set_filters_c(Tally* tally, int n, int32_t filter_indices[])
+  {tally->set_filters(filter_indices, n);}
+
+  int tally_get_n_filters_c(Tally* tally) {return tally->filters().size();}
+
+  int32_t tally_get_filter_c(Tally* tally, int i) {return tally->filters(i);}
+
+  int32_t tally_get_n_filter_bins_c(Tally* tally)
+  {return tally->n_filter_bins();}
+
+  int tally_get_n_nuclide_bins_c(Tally* tally)
+  {return tally->nuclides_.size();}
+
+  int tally_get_nuclide_bins_c(Tally* tally, int i)
+  {return tally->nuclides_[i-1];}
+
+  int tally_get_energyout_filter_c(Tally* tally)
+  {return tally->energyout_filter_;}
+
+  void tally_set_scores(Tally* tally, pugi::xml_node* node)
+  {tally->set_scores(*node);}
+
+  void tally_set_nuclides(Tally* tally, pugi::xml_node* node)
+  {tally->set_nuclides(*node);}
+
+  void tally_init_triggers(Tally* tally, int i_tally, pugi::xml_node* node)
+  {tally->init_triggers(*node, i_tally);}
+
+  int tally_get_deriv_c(Tally* tally) {return tally->deriv_;}
+
+  int tally_set_deriv_c(Tally* tally, int deriv) {tally->deriv_ = deriv;}
+}
+
 } // namespace openmc

src/tallies/tally_derivative_header.F90

@@ -1,103 +1,33 @@
 module tally_derivative_header
 
-  use constants
-  use dict_header, only: DictIntInt
-  use error, only: fatal_error
-  use nuclide_header, only: nuclide_map_get
-  use string, only: to_str, to_lower, to_c_string
-  use xml_interface
+  use, intrinsic :: ISO_C_BINDING
 
   implicit none
-  private
-  public :: free_memory_tally_derivative
 
 !===============================================================================
 ! TALLYDERIVATIVE describes a first-order derivative that can be applied to
 ! tallies.
 !===============================================================================
 
-  type, public :: TallyDerivative
-    integer :: id
-    integer :: variable
-    integer :: diff_material
-    integer :: diff_nuclide
-    real(8) :: flux_deriv
-  contains
-    procedure :: from_xml
+  type, bind(C), public :: TallyDerivative
+    integer(C_INT) :: id
+    integer(C_INT) :: variable
+    integer(C_INT) :: diff_material
+    integer(C_INT) :: diff_nuclide
+    real(C_DOUBLE) :: flux_deriv
   end type TallyDerivative
 
-  type(TallyDerivative), public, allocatable :: tally_derivs(:)
-!$omp threadprivate(tally_derivs)
+  interface
+    function n_tally_derivs() result(n) bind(C)
+      import C_INT
+      integer(C_INT) :: n
+    end function
 
-  ! Dictionary that maps user IDs to indices in 'tally_derivs'
-  type(DictIntInt), public :: tally_deriv_dict
-
-contains
-
-  subroutine from_xml(this, node)
-    class(TallyDerivative), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    character(MAX_WORD_LEN) :: temp_str
-    character(MAX_WORD_LEN) :: word
-    integer :: val
-
-    ! Copy the derivative id.
-    if (check_for_node(node, "id")) then
-      call get_node_value(node, "id", this % id)
-    else
-      call fatal_error("Must specify an ID for <derivative> elements in the&
-           & tally XML file")
-    end if
-
-    ! Make sure the id is > 0.
-    if (this % id <= 0) then
-      call fatal_error("<derivative> IDs must be an integer greater than &
-           &zero")
-    end if
-
-    ! Make sure this id has not already been used.
-    if (tally_deriv_dict % has(this % id)) then
-      call fatal_error("Two or more <derivative>'s use the same unique &
-           &ID: " // trim(to_str(this % id)))
-    end if
-
-    ! Read the independent variable name.
-    call get_node_value(node, "variable", temp_str)
-    temp_str = to_lower(temp_str)
-
-    select case(temp_str)
-    case("density")
-      this % variable = DIFF_DENSITY
-
-    case("nuclide_density")
-      this % variable = DIFF_NUCLIDE_DENSITY
-
-      call get_node_value(node, "nuclide", word)
-      val = nuclide_map_get(to_c_string(word))
-      if (val == -1) then
-        call fatal_error("Could not find the nuclide " &
-             // trim(word) // " specified in derivative " &
-             // trim(to_str(this % id)) // " in any material.")
-      end if
-      this % diff_nuclide = val
-
-    case("temperature")
-      this % variable = DIFF_TEMPERATURE
-    end select
-
-    call get_node_value(node, "material", this % diff_material)
-
-  end subroutine from_xml
-
-!===============================================================================
-! FREE_MEMORY_TALLY_DERIVATIVE deallocates global arrays defined in this module
-!===============================================================================
-
-  subroutine free_memory_tally_derivative()
-!$omp parallel
-    if (allocated(tally_derivs)) deallocate(tally_derivs)
-!$omp end parallel
-  end subroutine free_memory_tally_derivative
+    function tally_deriv_c(i) result(deriv) bind(C)
+      import C_INT, TallyDerivative
+      integer(C_INT), value :: i
+      type(TallyDerivative), pointer :: deriv
+    end function
+  end interface
 
 end module tally_derivative_header

src/tallies/tally_filter.F90

@@ -7,11 +7,7 @@ module tally_filter
   use tally_filter_header
 
   ! Inherit other filters
-  use tally_filter_delayedgroup
   use tally_filter_distribcell
-  use tally_filter_energy
-  use tally_filter_mesh
-  use tally_filter_meshsurface
   use tally_filter_particle
   use tally_filter_sph_harm
 

src/tallies/tally_filter_delayedgroup.F90

@@ -1,38 +0,0 @@
-module tally_filter_delayedgroup
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use tally_filter_header
-
-  implicit none
-  private
-
-!===============================================================================
-! DELAYEDGROUPFILTER bins outgoing fission neutrons in their delayed groups.
-! The get_all_bins functionality is not actually used.  The bins are manually
-! iterated over in the scoring subroutines.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: DelayedGroupFilter
-  contains
-    procedure :: groups
-  end type DelayedGroupFilter
-
-contains
-
-  function groups(this, i) result(group)
-    class(DelayedGroupFilter), intent(in) :: this
-    integer,                   intent(in) :: i
-    integer                               :: group
-    interface
-      function delayedgroup_filter_groups(filt, i) result(group) bind(C)
-        import C_PTR, C_INT
-        type(C_PTR),    value :: filt
-        integer(C_INT), value :: i
-        integer(C_INT)        :: group
-      end function
-    end interface
-    group = delayedgroup_filter_groups(this % ptr, i)
-  end function groups
-
-end module tally_filter_delayedgroup

src/tallies/tally_filter_energy.F90

@@ -1,91 +0,0 @@
-module tally_filter_energy
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-  public :: openmc_energy_filter_get_bins
-  public :: openmc_energy_filter_set_bins
-
-  interface
-    function openmc_energy_filter_get_bins(index, energies, n) result(err) bind(C)
-      import C_INT32_T, C_PTR, C_INT
-      integer(C_INT32_T), value :: index
-      type(C_PTR), intent(out) :: energies
-      integer(C_INT32_T), intent(out) :: n
-      integer(C_INT) :: err
-    end function
-    function openmc_energy_filter_set_bins(index, n, energies) result(err) bind(C)
-      import C_INT32_T, C_DOUBLE, C_INT
-      integer(C_INT32_T), value, intent(in) :: index
-      integer(C_INT32_T), value, intent(in) :: n
-      real(C_DOUBLE), intent(in) :: energies(n)
-      integer(C_INT) :: err
-    end function
-  end interface
-
-!===============================================================================
-! ENERGYFILTER bins the incident neutron energy.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: EnergyFilter
-    ! True if transport group number can be used directly to get bin number
-    logical :: matches_transport_groups = .false.
-  contains
-    procedure :: from_xml => from_xml_energy
-    procedure :: search
-  end type EnergyFilter
-
-!===============================================================================
-! ENERGYOUTFILTER bins the outgoing neutron energy.  Only scattering events use
-! the get_all_bins functionality.  Nu-fission tallies manually iterate over the
-! filter bins.
-!===============================================================================
-
-  type, public, extends(EnergyFilter) :: EnergyoutFilter
-  end type EnergyoutFilter
-
-contains
-
-!===============================================================================
-! EnergyFilter methods
-!===============================================================================
-
-  subroutine from_xml_energy(this, node)
-    class(EnergyFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    interface
-      function energy_filter_matches_transport_groups(filt) result(matches) &
-           bind(C)
-        import C_PTR, C_BOOL
-        type(C_PTR), value :: filt
-        logical(C_BOOL)    :: matches
-      end function
-    end interface
-
-    call this % from_xml_cpp(node)
-    this % n_bins = this % n_bins_cpp()
-    this % matches_transport_groups = &
-         energy_filter_matches_transport_groups(this % ptr)
-  end subroutine from_xml_energy
-
-  function search(this, val) result(bin)
-    class(EnergyFilter), intent(in) :: this
-    real(8),             intent(in) :: val
-    integer                         :: bin
-    interface
-      function energy_filter_search(filt, val) result(bin) bind(C)
-        import C_PTR, C_DOUBLE, C_INT
-        type(C_PTR), value    :: filt
-        real(C_DOUBLE), value :: val
-        integer(C_INT)        :: bin
-      end function
-    end interface
-    bin = energy_filter_search(this % ptr, val)
-  end function search
-
-end module tally_filter_energy

src/tallies/tally_filter_header.F90

@@ -19,39 +19,6 @@ module tally_filter_header
   public :: openmc_filter_set_id
   public :: openmc_get_filter_index
   public :: openmc_get_filter_next_id
-  public :: filter_match_pointer
-
-  interface
-    function filter_match_pointer(indx) bind(C) result(ptr)
-      import C_PTR, C_INT
-      integer(C_INT), intent(in), value :: indx
-      type(C_PTR)                       :: ptr
-    end function filter_match_pointer
-  end interface
-
-!===============================================================================
-! TALLYFILTERMATCH stores every valid bin and weight for a filter
-!===============================================================================
-
-  type, public :: TallyFilterMatch
-    type(C_PTR) :: ptr
-
-    ! Index of the bin and weight being used in the current filter combination
-    integer          :: i_bin
-
-    ! Indicates whether all valid bins for this filter have been found
-    logical          :: bins_present = .false.
-
-  contains
-    procedure :: bins_push_back
-    procedure :: weights_push_back
-    procedure :: bins_clear
-    procedure :: weights_clear
-    procedure :: bins_size
-    procedure :: bins_data
-    procedure :: weights_data
-    procedure :: bins_set_data
-  end type TallyFilterMatch
 
 !===============================================================================
 ! TALLYFILTER describes a filter that limits what events score to a tally. For
@@ -60,14 +27,12 @@ module tally_filter_header
 !===============================================================================
 
   type, public, abstract :: TallyFilter
-    integer :: id
     integer :: n_bins = 0
     type(C_PTR) :: ptr
   contains
+    procedure :: id => filter_get_id_f90
     procedure :: from_xml
-    procedure :: get_all_bins
     procedure :: to_statepoint
-    procedure :: text_label
     procedure :: initialize
     procedure :: n_bins_cpp
     procedure :: from_xml_cpp
@@ -90,6 +55,15 @@ module tally_filter_header
   type, public, extends(CellFilter) :: CellFromFilter
   end type
 
+  type, public, extends(TallyFilter) :: DelayedGroupFilter
+  end type
+
+  type, public, extends(TallyFilter) :: EnergyFilter
+  end type
+
+  type, public, extends(EnergyFilter) :: EnergyoutFilter
+  end type
+
   type, public, extends(TallyFilter) :: EnergyFunctionFilter
   end type
 
@@ -99,6 +73,12 @@ module tally_filter_header
   type, public, extends(TallyFilter) :: MaterialFilter
   end type
 
+  type, public, extends(TallyFilter) :: MeshFilter
+  end type
+
+  type, public, extends(TallyFilter) :: MeshSurfaceFilter
+  end type
+
   type, public, extends(TallyFilter) :: MuFilter
   end type
 
@@ -134,8 +114,6 @@ module tally_filter_header
   integer(C_INT32_T), public, bind(C) :: n_filters = 0 ! # of filters
 
   type(TallyFilterContainer), public, allocatable, target :: filters(:)
-  type(TallyFilterMatch), public, allocatable :: filter_matches(:)
-!$omp threadprivate(filter_matches)
 
   ! Dictionary that maps user IDs to indices in 'filters'
   type(DictIntInt), public :: filter_dict
@@ -146,120 +124,23 @@ module tally_filter_header
 
 contains
 
-!===============================================================================
-! TallyFilterMatch implementation
-!===============================================================================
-
-  subroutine bins_push_back(this, val)
-    class(TallyFilterMatch), intent(inout) :: this
-    integer,                 intent(in)    :: val
-    interface
-      subroutine filter_match_bins_push_back(ptr, val) bind(C)
-        import C_PTR, C_INT
-        type(C_PTR),                value :: ptr
-        integer(C_INT), intent(in), value :: val
-      end subroutine
-    end interface
-    call filter_match_bins_push_back(this % ptr, val)
-  end subroutine bins_push_back
-
-  subroutine weights_push_back(this, val)
-    class(TallyFilterMatch), intent(inout) :: this
-    real(8),                 intent(in)    :: val
-    interface
-      subroutine filter_match_weights_push_back(ptr, val) bind(C)
-        import C_PTR, C_DOUBLE
-        type(C_PTR),                value :: ptr
-        real(C_DOUBLE), intent(in), value :: val
-      end subroutine
-    end interface
-    call filter_match_weights_push_back(this % ptr, val)
-  end subroutine weights_push_back
-
-  subroutine bins_clear(this)
-    class(TallyFilterMatch), intent(inout) :: this
-    interface
-      subroutine filter_match_bins_clear(ptr) bind(C)
-        import C_PTR
-        type(C_PTR), value :: ptr
-      end subroutine
-    end interface
-    call filter_match_bins_clear(this % ptr)
-  end subroutine bins_clear
-
-  subroutine weights_clear(this)
-    class(TallyFilterMatch), intent(inout) :: this
-    interface
-      subroutine filter_match_weights_clear(ptr) bind(C)
-        import C_PTR
-        type(C_PTR), value :: ptr
-      end subroutine
-    end interface
-    call filter_match_weights_clear(this % ptr)
-  end subroutine weights_clear
-
-  function bins_size(this) result(len)
-    class(TallyFilterMatch), intent(inout) :: this
-    integer                                :: len
-    interface
-      function filter_match_bins_size(ptr) bind(C) result(len)
-        import C_PTR, C_INT
-        type(C_PTR), value :: ptr
-        integer(C_INT)     :: len
-      end function
-    end interface
-    len = filter_match_bins_size(this % ptr)
-  end function bins_size
-
-  function bins_data(this, indx) result(val)
-    class(TallyFilterMatch), intent(inout) :: this
-    integer,                 intent(in)    :: indx
-    integer                                :: val
-    interface
-      function filter_match_bins_data(ptr, indx) bind(C) result(val)
-        import C_PTR, C_INT
-        type(C_PTR),                value :: ptr
-        integer(C_INT), intent(in), value :: indx
-        integer(C_INT)                    :: val
-      end function
-    end interface
-    val = filter_match_bins_data(this % ptr, indx)
-  end function bins_data
-
-  function weights_data(this, indx) result(val)
-    class(TallyFilterMatch), intent(inout) :: this
-    integer,                 intent(in)    :: indx
-    real(8)                                :: val
-    interface
-      function filter_match_weights_data(ptr, indx) bind(C) result(val)
-        import C_PTR, C_INT, C_DOUBLE
-        type(C_PTR),                value :: ptr
-        integer(C_INT), intent(in), value :: indx
-        real(C_DOUBLE)                    :: val
-      end function
-    end interface
-    val = filter_match_weights_data(this % ptr, indx)
-  end function weights_data
-
-  subroutine bins_set_data(this, indx, val)
-    class(TallyFilterMatch), intent(inout) :: this
-    integer,                 intent(in)    :: indx
-    integer,                 intent(in)    :: val
-    interface
-      subroutine filter_match_bins_set_data(ptr, indx, val) bind(C)
-        import C_PTR, C_INT
-        type(C_PTR),    value             :: ptr
-        integer(C_INT), value, intent(in) :: indx
-        integer(C_INT), value, intent(in) :: val
-      end subroutine
-    end interface
-    call filter_match_bins_set_data(this % ptr, indx, val)
-  end subroutine bins_set_data
-
 !===============================================================================
 ! TallyFilter implementation
 !===============================================================================
 
+  function filter_get_id_f90(this) result(id)
+    class(TallyFilter) :: this
+    integer(C_INT32_T) :: id
+    interface
+      function filter_get_id(filt) result(id) bind(C)
+        import C_PTR, C_INT32_T
+        type(C_PTR), value :: filt
+        integer(C_INT32_T) :: id
+      end function
+    end interface
+    id = filter_get_id(this % ptr)
+  end function
+
   subroutine from_xml(this, node)
     class(TallyFilter), intent(inout) :: this
     type(XMLNode),      intent(in)    :: node
@@ -267,23 +148,6 @@ contains
     this % n_bins = this % n_bins_cpp()
   end subroutine from_xml
 
-  subroutine get_all_bins(this, p, estimator, match)
-    class(TallyFilter),     intent(in)    :: this
-    type(Particle),         intent(in)    :: p
-    integer,                intent(in)    :: estimator
-    type(TallyFilterMatch), intent(inout) :: match
-    interface
-      subroutine filter_get_all_bins(filt, p, estimator, match) bind(C)
-        import C_PTR, Particle, C_INT
-        type(C_PTR),                value :: filt
-        type(Particle), intent(in)        :: p
-        integer(C_INT), intent(in), value :: estimator
-        type(C_PTR),                value :: match
-      end subroutine filter_get_all_bins
-    end interface
-    call filter_get_all_bins(this % ptr, p, estimator, match % ptr)
-  end subroutine get_all_bins
-
   subroutine to_statepoint(this, filter_group)
     class(TallyFilter), intent(in) :: this
     integer(HID_T),     intent(in) :: filter_group
@@ -297,28 +161,6 @@ contains
     call filter_to_statepoint(this % ptr, filter_group)
   end subroutine to_statepoint
 
-  function text_label(this, bin) result(label)
-    class(TallyFilter), intent(in) :: this
-    integer,            intent(in) :: bin
-    character(MAX_LINE_LEN)        :: label
-    character(kind=C_CHAR)         :: label_(MAX_LINE_LEN+1)
-    integer :: i
-    interface
-      subroutine filter_text_label(filt, bin, label) bind(C)
-        import C_PTR, C_INT, C_CHAR
-        type(C_PTR), value     :: filt
-        integer(C_INT), value  :: bin
-        character(kind=C_CHAR) :: label(*)
-      end subroutine filter_text_label
-    end interface
-    call filter_text_label(this % ptr, bin, label_)
-    label = " "
-    do i = 1, MAX_LINE_LEN
-      if (label_(i) == C_NULL_CHAR) exit
-      label(i:i) = label_(i)
-    end do
-  end function text_label
-
   subroutine initialize(this)
     class(TallyFilter), intent(inout) :: this
     call this % initialize_cpp()
@@ -456,7 +298,7 @@ contains
     integer(C_INT) :: err
 
     if (index >= 1 .and. index <= n_filters) then
-      id = filters(index) % obj % id
+      id = filters(index) % obj % id()
       err = 0
     else
       err = E_OUT_OF_BOUNDS
@@ -471,9 +313,17 @@ contains
     integer(C_INT32_T), value, intent(in) :: id
     integer(C_INT) :: err
 
+    interface
+      subroutine filter_set_id(filt, id) bind(C)
+        import C_PTR, C_INT32_T
+        type(C_PTR), value :: filt
+        integer(C_INT32_T), value :: id
+      end subroutine
+    end interface
+
     if (index >= 1 .and. index <= n_filters) then
       if (allocated(filters(index) % obj)) then
-        filters(index) % obj % id = id
+        call filter_set_id(filters(index) % obj % ptr, id)
         call filter_dict % set(id, index)
         if (id > largest_filter_id) largest_filter_id = id
 

src/tallies/tally_filter_mesh.F90

@@ -1,38 +0,0 @@
-module tally_filter_mesh
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use tally_filter_header
-
-  implicit none
-
-  interface
-    function openmc_mesh_filter_set_mesh(index, index_mesh) result(err) bind(C)
-      import C_INT32_T, C_INT
-      integer(C_INT32_T), value, intent(in) :: index
-      integer(C_INT32_T), value, intent(in) :: index_mesh
-      integer(C_INT) :: err
-    end function
-  end interface
-
-  type, public, extends(TallyFilter) :: MeshFilter
-  contains
-    procedure :: mesh => get_mesh
-  end type MeshFilter
-
-contains
-
-  function get_mesh(this) result(mesh)
-    class(MeshFilter) :: this
-    integer :: mesh
-    interface
-      function mesh_filter_get_mesh(filt) result(index_mesh) bind(C)
-        import C_PTR, C_INT
-        type(C_PTR), value :: filt
-        integer(C_INT)     :: index_mesh
-      end function mesh_filter_get_mesh
-    end interface
-    mesh = mesh_filter_get_mesh(this % ptr)
-  end function get_mesh
-
-end module tally_filter_mesh

src/tallies/tally_filter_meshsurface.F90

@@ -1,22 +0,0 @@
-module tally_filter_meshsurface
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use tally_filter_header
-
-  implicit none
-
-  interface
-    function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) &
-         bind(C)
-      import C_INT32_T, C_INT
-      integer(C_INT32_T), value, intent(in) :: index
-      integer(C_INT32_T), value, intent(in) :: index_mesh
-      integer(C_INT) :: err
-    end function
-  end interface
-
-  type, extends(TallyFilter) :: MeshSurfaceFilter
-  end type MeshSurfaceFilter
-
-end module tally_filter_meshsurface

src/tallies/trigger.F90

@@ -1,443 +0,0 @@
-module trigger
-
-  use, intrinsic :: ISO_C_BINDING
-
-#ifdef OPENMC_MPI
-  use message_passing
-#endif
-
-  use constants
-  use eigenvalue,     only: openmc_get_keff
-  use error,          only: warning, write_message
-  use string,         only: to_str
-  use mesh_header,    only: RegularMesh, meshes
-  use message_passing, only: master
-  use settings
-  use simulation_header
-  use trigger_header
-  use tally,          only: TallyObject
-  use tally_filter_mesh, only: MeshFilter
-  use tally_filter_header, only: filter_matches, filters
-  use tally_header, only: tallies, n_tallies
-
-  implicit none
-
-contains
-
-!===============================================================================
-! CHECK_TRIGGERS checks any user-specified precision triggers' for convergence
-! and predicts the number of remainining batches to convergence.
-!===============================================================================
-
-  subroutine check_triggers() bind(C)
-
-    implicit none
-
-    ! Variables to reflect distance to trigger convergence criteria
-    real(8)            :: max_ratio       ! max uncertainty/thresh ratio
-    integer            :: tally_id        ! id for tally with max ratio
-    character(len=52)  :: name            ! "eigenvalue" or tally score
-
-    integer    :: n_pred_batches  ! predicted # batches to satisfy all triggers
-
-    ! Checks if current_batch is one for which the triggers must be checked
-    if (current_batch < n_batches .or. (.not. trigger_on)) return
-    if (mod((current_batch - n_batches), n_batch_interval) /= 0 .and. &
-         current_batch /= n_max_batches) return
-
-    ! Check the trigger and output the result
-    call check_tally_triggers(max_ratio, tally_id, name)
-
-    ! When trigger threshold is reached, write information
-    if (satisfy_triggers) then
-      call write_message("Triggers satisfied for batch " // &
-           trim(to_str(current_batch)), 7)
-
-    ! When trigger is not reached write convergence info for user
-    elseif (name == "eigenvalue") then
-      call write_message("Triggers unsatisfied, max unc./thresh. is " // &
-           trim(to_str(max_ratio)) //  " for " // trim(name), 7)
-    else
-      call write_message("Triggers unsatisfied, max unc./thresh. is " // &
-           trim(to_str(max_ratio)) // " for " // trim(name) // &
-           " in tally " // trim(to_str(tally_id)), 7)
-    end if
-
-    ! If batch_interval is not set, estimate batches till triggers are satisfied
-    if (pred_batches .and. .not. satisfy_triggers) then
-
-      ! Estimate the number of remaining batches to convergence
-      ! The prediction uses the fact that tally variances are proportional
-      ! to 1/N where N is the number of the batches/particles
-      n_batch_interval = int((current_batch-n_inactive) * &
-           (max_ratio ** 2)) + n_inactive-n_batches + 1
-      n_pred_batches = n_batch_interval + n_batches
-
-      ! Write the predicted number of batches for the user
-      if (n_pred_batches > n_max_batches) then
-        call warning("The estimated number of batches is " // &
-             trim(to_str(n_pred_batches)) // &
-             " --  greater than max batches. ")
-      else
-        call write_message("The estimated number of batches is " // &
-             trim(to_str(n_pred_batches)), 7)
-      end if
-    end if
-  end subroutine check_triggers
-
-
-!===============================================================================
-! CHECK_TALLY_TRIGGERS checks whether uncertainties are below the threshold,
-! and finds the maximum  uncertainty/threshold ratio for all triggers
-!===============================================================================
-
-  subroutine check_tally_triggers(max_ratio, tally_id, name)
-
-    ! Variables to reflect distance to trigger convergence criteria
-    real(8), intent(inout) :: max_ratio       ! max uncertainty/thresh ratio
-    integer, intent(inout) :: tally_id        ! id for tally with max ratio
-    character(len=52), intent(inout) :: name  ! "eigenvalue" or tally score
-
-    integer :: i              ! index in tallies array
-    integer :: j              ! index in tally filters
-    integer :: n              ! loop index for nuclides
-    integer :: s              ! loop index for triggers
-    integer :: filter_index   ! index in results array for filters
-    integer :: score_index    ! scoring bin index
-    integer(C_INT) :: err
-    real(8) :: uncertainty    ! trigger uncertainty
-    real(8) :: std_dev = ZERO ! trigger standard deviation
-    real(8) :: rel_err = ZERO ! trigger relative error
-    real(8) :: ratio          ! ratio of the uncertainty/trigger threshold
-    real(C_DOUBLE) :: k_combined(2)
-
-    ! Initialize tally trigger maximum uncertainty ratio to zero
-    max_ratio = 0
-
-    if (master) then
-
-      ! By default, assume all triggers are satisfied
-      satisfy_triggers = .true.
-
-      ! Check eigenvalue trigger
-      if (run_mode == MODE_EIGENVALUE) then
-        if (keff_trigger % trigger_type /= 0) then
-          err = openmc_get_keff(k_combined)
-          select case (keff_trigger % trigger_type)
-          case(VARIANCE)
-            uncertainty = k_combined(2) ** 2
-          case(STANDARD_DEVIATION)
-            uncertainty = k_combined(2)
-          case default
-            uncertainty = k_combined(2) / k_combined(1)
-          end select
-
-          ! If uncertainty is above threshold, store uncertainty ratio
-          if (uncertainty > keff_trigger % threshold) then
-            satisfy_triggers = .false.
-            if (keff_trigger % trigger_type == VARIANCE) then
-              ratio = sqrt(uncertainty / keff_trigger % threshold)
-            else
-              ratio = uncertainty / keff_trigger % threshold
-            end if
-            if (max_ratio < ratio) then
-              max_ratio = ratio
-              name = "eigenvalue"
-            end if
-          end if
-        end if
-      end if
-
-      ! Compute uncertainties for all tallies, scores with triggers
-      TALLY_LOOP: do i = 1, n_tallies
-        associate (t => tallies(i) % obj)
-
-        ! Cycle through if only one batch has been simumlate
-        if (t % n_realizations == 1) then
-          cycle TALLY_LOOP
-        end if
-
-        TRIGGER_LOOP: do s = 1, t % n_triggers
-          associate (trigger => t % triggers(s))
-
-          ! Initialize trigger uncertainties to zero
-          trigger % std_dev = ZERO
-          trigger % rel_err = ZERO
-          trigger % variance = ZERO
-
-          ! Mesh current tally triggers require special treatment
-          if (t % type == TALLY_MESH_SURFACE) then
-            call compute_tally_current(t, trigger)
-
-          else
-
-            ! Initialize bins, filter level
-            do j = 1, size(t % filter)
-              call filter_matches(t % filter(j)) % bins_clear()
-              call filter_matches(t % filter(j)) % bins_push_back(0)
-            end do
-
-            FILTER_LOOP: do filter_index = 1, t % n_filter_bins
-
-              ! Initialize score index
-              score_index = trigger % score_index
-
-              ! Initialize score bin index
-              NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins
-
-                call get_trigger_uncertainty(std_dev, rel_err, &
-                     score_index, filter_index, t)
-
-                if (trigger % variance < variance) then
-                  trigger % variance = std_dev ** 2
-                end if
-                if (trigger % std_dev < std_dev) then
-                  trigger % std_dev = std_dev
-                end if
-                if (trigger % rel_err < rel_err) then
-                  trigger % rel_err = rel_err
-                end if
-
-                select case (t % triggers(s) % type)
-                case(VARIANCE)
-                  uncertainty = trigger % variance
-                case(STANDARD_DEVIATION)
-                  uncertainty = trigger % std_dev
-                case default
-                  uncertainty = trigger % rel_err
-                end select
-
-                if (uncertainty > t % triggers(s) % threshold) then
-                  satisfy_triggers = .false.
-
-                  if (t % triggers(s) % type == VARIANCE) then
-                    ratio = sqrt(uncertainty / t % triggers(s) % threshold)
-                  else
-                    ratio = uncertainty / t % triggers(s) % threshold
-                  end if
-
-                  if (max_ratio < ratio) then
-                    max_ratio = ratio
-                    name  = t % triggers(s) % score_name
-                    tally_id = t % id
-                  end if
-                end if
-              end do NUCLIDE_LOOP
-              if (size(t % filter) == 0) exit FILTER_LOOP
-            end do FILTER_LOOP
-          end if
-          end associate
-        end do TRIGGER_LOOP
-        end associate
-      end do TALLY_LOOP
-    end if
-  end subroutine check_tally_triggers
-
-
-!===============================================================================
-! COMPUTE_TALLY_CURRENT computes the current for a mesh current tally with
-! precision trigger(s).
-!===============================================================================
-
-  subroutine compute_tally_current(t, trigger)
-    type(TallyObject),   intent(in)    :: t        ! mesh current tally
-    type(TriggerObject), intent(inout) :: trigger  ! mesh current tally trigger
-
-    integer :: i                    ! mesh index
-    integer :: j                    ! loop index for tally filters
-    integer :: ijk(3)               ! indices of mesh cells
-    integer :: n_dim                ! number of mesh dimensions
-    integer :: n_cells              ! number of mesh cells
-    integer :: l                    ! index for energy
-    integer :: i_filter_mesh        ! index for mesh filter
-    integer :: i_filter_ein         ! index for incoming energy filter
-    integer :: i_filter_surf        ! index for surface filter
-    integer :: n                    ! number of incoming energy bins
-    integer :: filter_index         ! index in results array for filters
-    logical :: print_ebin           ! should incoming energy bin be displayed?
-    real(8) :: rel_err  = ZERO      ! temporary relative error of result
-    real(8) :: std_dev  = ZERO      ! temporary standard deviration of result
-    type(RegularMesh) :: m        ! surface current mesh
-
-    ! Get pointer to mesh
-    i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
-    i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
-    select type(filt => filters(i_filter_mesh) % obj)
-    type is (MeshFilter)
-      m = meshes(filt % mesh())
-    end select
-
-    ! initialize bins array
-    do j = 1, size(t % filter)
-      call filter_matches(t % filter(j)) % bins_clear()
-      call filter_matches(t % filter(j)) % bins_push_back(1)
-    end do
-
-    ! determine how many energyin bins there are
-    i_filter_ein = t % find_filter(FILTER_ENERGYIN)
-    if (i_filter_ein > 0) then
-      print_ebin = .true.
-      n = filters(t % filter(i_filter_ein)) % obj % n_bins
-      i_filter_ein = t % filter(i_filter_ein)
-    else
-      print_ebin = .false.
-      n = 1
-    end if
-
-    ! Get the dimensions and number of cells in the mesh
-    n_dim = m % n_dimension()
-    n_cells = 1
-    do j = 1, n_dim
-      n_cells = n_cells * m % dimension(j)
-    end do
-
-    ! Loop over all the mesh cells
-    do i = 1, n_cells
-
-      ! Get the indices for this cell
-      call m % get_indices_from_bin(i, ijk)
-      call filter_matches(i_filter_mesh) % bins_set_data(1, i)
-
-      do l = 1, n
-
-        if (print_ebin) then
-          call filter_matches(i_filter_ein) % bins_set_data(1, l)
-        end if
-
-        ! Left Surface
-        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_LEFT)
-        filter_index = 1
-        do j = 1, size(t % filter)
-          filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins_data(1) - 1) * t % stride(j)
-        end do
-        call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-        if (trigger % std_dev < std_dev) then
-          trigger % std_dev = std_dev
-        end if
-        if (trigger % rel_err < rel_err) then
-          trigger % rel_err = rel_err
-        end if
-        trigger % variance = std_dev**2
-
-        ! Right Surface
-        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_RIGHT)
-        filter_index = 1
-        do j = 1, size(t % filter)
-          filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins_data(1) - 1) * t % stride(j)
-        end do
-        call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-        if (trigger % std_dev < std_dev) then
-          trigger % std_dev = std_dev
-        end if
-        if (trigger % rel_err < rel_err) then
-          trigger % rel_err = rel_err
-        end if
-        trigger % variance = trigger % std_dev**2
-
-        ! Back Surface
-        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BACK)
-        filter_index = 1
-        do j = 1, size(t % filter)
-          filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins_data(1) - 1) * t % stride(j)
-        end do
-        call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-        if (trigger % std_dev < std_dev) then
-          trigger % std_dev = std_dev
-        end if
-        if (trigger % rel_err < rel_err) then
-          trigger % rel_err = rel_err
-        end if
-        trigger % variance = trigger % std_dev**2
-
-        ! Front Surface
-        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_FRONT)
-        filter_index = 1
-        do j = 1, size(t % filter)
-          filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins_data(1) - 1) * t % stride(j)
-        end do
-        call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-        if (trigger % std_dev < std_dev) then
-          trigger % std_dev = std_dev
-        end if
-        if (trigger % rel_err < rel_err) then
-          trigger % rel_err = rel_err
-        end if
-        trigger % variance = trigger % std_dev**2
-
-        ! Bottom Surface
-        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BOTTOM)
-        filter_index = 1
-        do j = 1, size(t % filter)
-          filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins_data(1) - 1) * t % stride(j)
-        end do
-        call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-        if (trigger % std_dev < std_dev) then
-          trigger % std_dev = std_dev
-        end if
-        if (trigger % rel_err < rel_err) then
-          trigger % rel_err = rel_err
-        end if
-        trigger % variance = trigger % std_dev**2
-
-        ! Top Surface
-        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_TOP)
-        filter_index = 1
-        do j = 1, size(t % filter)
-          filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins_data(1) - 1) * t % stride(j)
-        end do
-        call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
-        if (trigger % std_dev < std_dev) then
-          trigger % std_dev = std_dev
-        end if
-        if (trigger % rel_err < rel_err) then
-          trigger % rel_err = rel_err
-        end if
-        trigger % variance = trigger % std_dev**2
-
-      end do
-    end do
-
-  end subroutine compute_tally_current
-
-!===============================================================================
-! GET_TRIGGER_UNCERTAINTY computes the standard deviation and relative error
-! for a single tally bin for CHECK_TALLY_TRIGGERS.
-!===============================================================================
-
-  subroutine get_trigger_uncertainty(std_dev, rel_err, score_index, &
-       filter_index, t)
-
-    real(8), intent(inout)     :: std_dev         ! tally standard deviation
-    real(8), intent(inout)     :: rel_err         ! tally relative error
-    integer, intent(in)        :: score_index     ! tally results score index
-    integer, intent(in)        :: filter_index    ! tally results filter index
-    type(TallyObject), intent(in) :: t            ! tally
-
-    integer :: n               ! number of realizations
-    real(8) :: mean            ! tally mean
-    real(8) :: tally_sum, tally_sum_sq     ! results for a single tally bin
-
-    n = t % n_realizations
-    tally_sum = t % results(RESULT_SUM, score_index, filter_index)
-    tally_sum_sq = t % results(RESULT_SUM_SQ, score_index, filter_index)
-
-    ! Compute the tally mean and standard deviation
-    mean    = tally_sum / n
-    std_dev = sqrt((tally_sum_sq / n - mean * mean) / (n - 1))
-
-    ! Compute the relative error if the mean is non-zero
-    if (mean == ZERO) then
-      rel_err = ZERO
-    else
-      rel_err = std_dev / mean
-    end if
-
-  end subroutine get_trigger_uncertainty
-
-end module trigger

src/tallies/trigger.cpp

@@ -0,0 +1,201 @@
+#include "openmc/tallies/trigger.h"
+
+#include <cmath>
+#include <sstream>
+#include <utility> // for std::pair
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/reaction.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+#include "openmc/tallies/tally.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variable definitions
+//==============================================================================
+
+namespace settings {
+  KTrigger keff_trigger;
+}
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+std::pair<double, double>
+get_tally_uncertainty(int i_tally, int score_index, int filter_index)
+{
+  int n;
+  //TODO: off-by-one
+  int err = openmc_tally_get_n_realizations(i_tally+1, &n);
+
+  auto results = tally_results(i_tally);
+  auto sum = results(filter_index, score_index, RESULT_SUM);
+  auto sum_sq = results(filter_index, score_index, RESULT_SUM_SQ);
+
+  auto mean = sum / n;
+  double std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1));
+  double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;
+
+  return {std_dev, rel_err};
+}
+
+//! Find the limiting limiting tally trigger.
+//
+//! param[out] ratio The uncertainty/threshold ratio for the most limiting
+//!   tally trigger
+//! param[out] tally_id The ID number of the most limiting tally
+//! param[out] score The most limiting tally score bin
+
+void
+check_tally_triggers(double& ratio, int& tally_id, int& score)
+{
+  ratio = 0.;
+  for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
+    const Tally& t {*model::tallies[i_tally]};
+
+    // Ignore tallies with less than two realizations.
+    int n_reals;
+    //TODO: off-by-one
+    int err = openmc_tally_get_n_realizations(i_tally+1, &n_reals);
+    if (n_reals < 2) continue;
+
+    for (const auto& trigger : t.triggers_) {
+      const auto& results = tally_results(i_tally);
+      for (auto filter_index = 0; filter_index < results.shape()[0];
+           ++filter_index) {
+        for (auto score_index = 0; score_index < results.shape()[1];
+             ++score_index) {
+          // Compute the tally uncertainty metrics.
+          auto uncert_pair = get_tally_uncertainty(i_tally, score_index,
+            filter_index);
+          double std_dev = uncert_pair.first;
+          double rel_err = uncert_pair.second;
+
+          // Pick out the relevant uncertainty metric for this trigger.
+          double uncertainty;
+          switch (trigger.metric) {
+            case TriggerMetric::variance:
+              uncertainty = std_dev * std_dev;
+              break;
+            case TriggerMetric::standard_deviation:
+              uncertainty = std_dev;
+              break;
+            case TriggerMetric::relative_error:
+              uncertainty = rel_err;
+              break;
+          }
+
+          // Compute the uncertainty / threshold ratio.
+          double this_ratio = uncertainty / trigger.threshold;
+          if (trigger.metric == TriggerMetric::variance) {
+            this_ratio = std::sqrt(ratio);
+          }
+
+          // If this is the most uncertain value, set the output variables.
+          if (this_ratio > ratio) {
+            ratio = this_ratio;
+            score = t.scores_[trigger.score_index];
+            tally_id = t.id_;
+          }
+        }
+      }
+    }
+  }
+}
+
+//! Compute the uncertainty/threshold ratio for the eigenvalue trigger.
+
+double
+check_keff_trigger()
+{
+  if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.;
+  if (settings::keff_trigger.metric == TriggerMetric::not_active) return 0.;
+
+  double k_combined[2];
+  int err = openmc_get_keff(k_combined);
+
+  double uncertainty = 0.;
+  switch (settings::keff_trigger.metric) {
+  case TriggerMetric::variance:
+    uncertainty = k_combined[1] * k_combined[1];
+    break;
+  case TriggerMetric::standard_deviation:
+    uncertainty = k_combined[1];
+    break;
+  case TriggerMetric::relative_error:
+    uncertainty = k_combined[1] / k_combined[0];
+  }
+
+  double ratio = uncertainty / settings::keff_trigger.threshold;
+  if (settings::keff_trigger.metric == TriggerMetric::variance)
+    ratio = std::sqrt(ratio);
+  return ratio;
+}
+
+//! See if tally and eigenvalue uncertainties are under trigger thresholds.
+
+void
+check_triggers()
+{
+  // Make some aliases.
+  const auto current_batch {simulation::current_batch};
+  const auto n_batches {settings::n_batches};
+  const auto interval {settings::trigger_batch_interval};
+
+  // See if the current batch is one for which the triggers must be checked.
+  if (!settings::trigger_on) return;
+  if (current_batch < n_batches) return;
+  if (((current_batch - n_batches) % interval) != 0) return;
+
+  // Check the eigenvalue and tally triggers.
+  double keff_ratio = check_keff_trigger();
+  double tally_ratio;
+  int tally_id, score;
+  check_tally_triggers(tally_ratio, tally_id, score);
+
+  // If all the triggers are satisfied, alert the user and return.
+  if (std::max(keff_ratio, tally_ratio) <= 1.) {
+    simulation::satisfy_triggers = true;
+    write_message("Triggers satisfied for batch "
+      + std::to_string(current_batch), 7);
+    return;
+  }
+
+  // At least one trigger is unsatisfied.  Let the user know which one.
+  simulation::satisfy_triggers = false;
+  std::stringstream msg;
+  msg << "Triggers unsatisfied, max unc./thresh. is ";
+  if (keff_ratio >= tally_ratio) {
+    msg << keff_ratio << " for eigenvalue";
+  } else {
+    msg << tally_ratio << " for " << reaction_name(score) << " in tally "
+        << tally_id;
+  }
+  write_message(msg, 7);
+
+  // Estimate batches til triggers are satisfied.
+  if (settings::trigger_predict) {
+    // This calculation assumes tally variance is proportional to 1/N where N is
+    // the number of batches.
+    auto max_ratio = std::max(keff_ratio, tally_ratio);
+    auto n_active = current_batch - settings::n_inactive;
+    auto n_pred_batches = static_cast<int>(n_active * max_ratio * max_ratio)
+      + settings::n_inactive + 1;
+
+    std::stringstream msg;
+    msg << "The estimated number of batches is " << n_pred_batches;
+    if (n_pred_batches > settings::n_max_batches) {
+      msg << " --- greater than max batches";
+      warning(msg);
+    } else {
+      write_message(msg, 7);
+    }
+  }
+}
+
+} // namespace openmc

src/tallies/trigger_header.F90

@@ -1,34 +0,0 @@
-module trigger_header
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants, only: NONE, N_FILTER_TYPES, ZERO
-
-  implicit none
-  private
-
-!===============================================================================
-! TRIGGEROBJECT stores the variance, relative error and standard deviation
-! for some user-specified trigger.
-!===============================================================================
-  type, public :: TriggerObject
-    integer            :: type          ! "variance", "std_dev" or "rel_err"
-    real(8)            :: threshold     ! a convergence threshold
-    character(len=52)  :: score_name    ! the name of the score
-    integer            :: score_index   ! the index of the score
-    real(8)            :: variance = ZERO ! temp variance container
-    real(8)            :: std_dev  = ZERO ! temp std. dev. container
-    real(8)            :: rel_err  = ZERO ! temp rel. err. container
-  end type TriggerObject
-
-!===============================================================================
-! KTRIGGER describes a user-specified precision trigger for k-effective
-!===============================================================================
-  type, public, bind(C) :: KTrigger
-    integer(C_INT)    :: trigger_type = 0
-    real(C_DOUBLE)    :: threshold    = ZERO
-  end type KTrigger
-
-  type(KTrigger), public, bind(C) :: keff_trigger  ! trigger for k-effective
-
-end module trigger_header

src/tracking.F90

@@ -24,6 +24,7 @@ module tracking
   use tally_header
   use tally,              only: score_analog_tally, score_tracklength_tally, &
                                 score_collision_tally, score_surface_tally, &
+                                score_meshsurface_tally, &
                                 score_track_derivative, zero_flux_derivs, &
                                 score_collision_derivative
   use track_output,       only: initialize_particle_track, write_particle_track, &
@@ -86,7 +87,7 @@ contains
     endif
 
     ! Every particle starts with no accumulated flux derivative.
-    if (active_tallies % size() > 0) call zero_flux_derivs()
+    if (active_tallies_size() > 0) call zero_flux_derivs()
 
     EVENT_LOOP: do
       ! Set the random number stream
@@ -172,7 +173,7 @@ contains
       end do
 
       ! Score track-length tallies
-      if (active_tracklength_tallies % size() > 0) then
+      if (active_tracklength_tallies_size() > 0) then
         call score_tracklength_tally(p, distance)
       end if
 
@@ -183,7 +184,7 @@ contains
       end if
 
       ! Score flux derivative accumulators for differential tallies.
-      if (active_tallies % size() > 0) call score_track_derivative(p, distance)
+      if (active_tallies_size() > 0) call score_track_derivative(p, distance)
 
       if (d_collision > d_boundary) then
         ! ====================================================================
@@ -210,8 +211,7 @@ contains
           p % event = EVENT_SURFACE
         end if
         ! Score cell to cell partial currents
-        if(active_surface_tallies % size() > 0) &
-             call score_surface_tally(p, active_surface_tallies)
+        if (active_surface_tallies_size() > 0) call score_surface_tally(p)
       else
         ! ====================================================================
         ! PARTICLE HAS COLLISION
@@ -226,8 +226,7 @@ contains
         ! since the direction of the particle will change and we need to use the
         ! pre-collision direction to figure out what mesh surfaces were crossed
 
-        if (active_meshsurf_tallies % size() > 0) &
-             call score_surface_tally(p, active_meshsurf_tallies)
+        if (active_meshsurf_tallies_size() > 0) call score_meshsurface_tally(p)
 
         ! Clear surface component
         p % surface = ERROR_INT
@@ -241,8 +240,8 @@ contains
         ! Score collision estimator tallies -- this is done after a collision
         ! has occurred rather than before because we need information on the
         ! outgoing energy for any tallies with an outgoing energy filter
-        if (active_collision_tallies % size() > 0) call score_collision_tally(p)
-        if (active_analog_tallies % size() > 0) call score_analog_tally(p)
+        if (active_collision_tallies_size() > 0) call score_collision_tally(p)
+        if (active_analog_tallies_size() > 0) call score_analog_tally(p)
 
         ! Reset banked weight during collision
         p % n_bank   = 0
@@ -273,7 +272,7 @@ contains
         end do
 
         ! Score flux derivative accumulators for differential tallies.
-        if (active_tallies % size() > 0) call score_collision_derivative(p)
+        if (active_tallies_size() > 0) call score_collision_derivative(p)
       end if
 
       ! If particle has too many events, display warning and kill it
@@ -346,12 +345,12 @@ contains
       ! forward slightly so that if the mesh boundary is on the surface, it is
       ! still processed
 
-      if (active_meshsurf_tallies % size() > 0) then
+      if (active_meshsurf_tallies_size() > 0) then
         ! TODO: Find a better solution to score surface currents than
         ! physically moving the particle forward slightly
 
         p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
-        call score_surface_tally(p, active_meshsurf_tallies)
+        call score_meshsurface_tally(p)
       end if
 
       ! Score to global leakage tally
@@ -382,14 +381,13 @@ contains
       ! the particle slightly back in case the surface crossing is coincident
       ! with a mesh boundary
 
-      if(active_surface_tallies % size() > 0) &
-           call score_surface_tally(p, active_surface_tallies)
+      if (active_surface_tallies_size() > 0) call score_surface_tally(p)
 
 
-      if (active_meshsurf_tallies % size() > 0) then
+      if (active_meshsurf_tallies_size() > 0) then
         xyz = p % coord(1) % xyz
         p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
-        call score_surface_tally(p, active_meshsurf_tallies)
+        call score_meshsurface_tally(p)
         p % coord(1) % xyz = xyz
       end if
 
@@ -443,10 +441,10 @@ contains
       ! Score surface currents since reflection causes the direction of the
       ! particle to change -- artificially move the particle slightly back in
       ! case the surface crossing is coincident with a mesh boundary
-      if (active_meshsurf_tallies % size() > 0) then
+      if (active_meshsurf_tallies_size() > 0) then
         xyz = p % coord(1) % xyz
         p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
-        call score_surface_tally(p, active_meshsurf_tallies)
+        call score_meshsurface_tally(p)
         p % coord(1) % xyz = xyz
       end if
 

tests/unit_tests/test_capi.py

@@ -156,6 +156,7 @@ def test_tally_mapping(capi_init):
 
 def test_tally(capi_init):
     t = openmc.capi.tallies[1]
+    assert t.type == 'volume'
     assert len(t.filters) == 2
     assert isinstance(t.filters[0], openmc.capi.MaterialFilter)
     assert isinstance(t.filters[1], openmc.capi.EnergyFilter)