Create C++ TallyDerivative struct

CMakeLists.txt

@@ -419,6 +419,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

include/openmc/constants.h

@@ -367,24 +367,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};
@@ -415,11 +397,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/nuclide.h

@@ -25,10 +25,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
@@ -63,10 +63,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
@@ -82,9 +82,9 @@ struct MaterialMacroXS {
   double pair_production; //!< macroscopic pair production xs
 };
 
-//===============================================================================
+//==============================================================================
 // Data for a nuclide
-//===============================================================================
+//==============================================================================
 
 class Nuclide {
 public:

include/openmc/tallies/derivative.h

@@ -0,0 +1,51 @@
+#ifndef OPENMC_TALLIES_DERIVATIVE_H
+#define OPENMC_TALLIES_DERIVATIVE_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() {}
+  TallyDerivative(pugi::xml_node node);
+};
+
+} // 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;
+  }
+
+  // Independent variables
+  //TODO: convert to enum
+  constexpr int DIFF_DENSITY {1};
+  constexpr int DIFF_NUCLIDE_DENSITY {2};
+  constexpr int DIFF_TEMPERATURE {3};
+}
+
+#endif // OPENMC_TALLIES_DERIVATIVE_H

src/constants.F90

@@ -252,30 +252,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 :: &

src/input_xml.F90

@@ -949,6 +949,13 @@ contains
     type(XMLNode), allocatable :: node_deriv_list(:)
     type(DictEntryCI) :: elem
 
+    interface
+      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)
@@ -999,6 +1006,8 @@ contains
     ! ==========================================================================
     ! READ DATA FOR DERIVATIVES
 
+    call read_tally_derivatives(root % ptr)
+
     ! 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)

src/tallies/derivative.cpp

@@ -0,0 +1,103 @@
+#include "openmc/error.h"
+#include "openmc/nuclide.h"
+#include "openmc/settings.h"
+#include "openmc/tallies/derivative.h"
+#include "openmc/xml_interface.h"
+
+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 = get_node_value(node, "variable");
+
+  if (variable == "density") {
+    variable = DIFF_DENSITY;
+
+  } else if (variable == "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;
+        diff_nuclide = i;
+      }
+    }
+    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 == "temperature") {
+    variable = DIFF_TEMPERATURE;
+
+  } else {
+    std::stringstream out;
+    out << "Unrecognized variable \"" << variable << "\" 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.push_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");
+}
+
+}// namespace openmc

src/tallies/trigger_header.F90

@@ -2,7 +2,7 @@ module trigger_header
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants, only: NONE, N_FILTER_TYPES, ZERO
+  use constants, only: NONE, ZERO
 
   implicit none
   private