Sensitivity analysis for tallies

CMakeLists.txt

@@ -440,6 +440,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_distribcell.cpp
   src/tallies/filter_energy.cpp
   src/tallies/filter_energyfunc.cpp
+  src/tallies/filter_importance.cpp
   src/tallies/filter_legendre.cpp
   src/tallies/filter_material.cpp
   src/tallies/filter_materialfrom.cpp
@@ -463,6 +464,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_zernike.cpp
   src/tallies/tally.cpp
   src/tallies/tally_scoring.cpp
+  src/tallies/sensitivity.cpp
   src/tallies/trigger.cpp
   src/thermal.cpp
   src/timer.cpp

docs/source/pythonapi/base.rst

@@ -128,6 +128,7 @@ Constructing Tallies
    openmc.SurfaceFilter
    openmc.MeshFilter
    openmc.MeshBornFilter
+   openmc.ImportanceFilter
    openmc.MeshMaterialFilter
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
@@ -153,6 +154,8 @@ Constructing Tallies
    openmc.MeshMaterialVolumes
    openmc.Trigger
    openmc.TallyDerivative
+   openmc.SensitivityTally
+   openmc.Sensitivity
    openmc.Tally
    openmc.Tallies
 

docs/source/pythonapi/capi.rst

@@ -77,6 +77,7 @@ Classes
    Mesh
    MeshFilter
    MeshBornFilter
+   ImportanceFilter
    MeshSurfaceFilter
    MuFilter
    Nuclide
@@ -93,6 +94,7 @@ Classes
    SolidRayTracePlot
    SurfaceFilter
    Tally
+   SensitivityTally
    TemporarySession
    UniverseFilter
    UnstructuredMesh

include/openmc/constants.h

@@ -296,6 +296,8 @@ enum class MgxsType {
 
 enum class TallyResult { VALUE, SUM, SUM_SQ, SUM_THIRD, SUM_FOURTH };
 
+enum class SensitivityTallyResult { VALUE, SUM, SUM_SQ, PREVIOUS_VALUE };
+
 enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT };
 
 enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };

include/openmc/filter_importance.h

@@ -0,0 +1,58 @@
+#ifndef OPENMC_TALLIES_FILTER_IMPORTANCE_H
+#define OPENMC_TALLIES_FILTER_IMPORTANCE_H
+
+#include <cstdint>
+
+#include <gsl/gsl-lite.hpp>
+#include "openmc/vector.h"
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Multiplies each tally by the importance corresponding to the mesh index 
+//==============================================================================
+
+class ImportanceFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~ImportanceFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override {return "importance";}
+  FilterType type() const override { return FilterType::IMPORTANCE; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const std::vector<double>& importance() const { return importance_; }
+  void set_importance(gsl::span<const double> importance);
+
+  virtual int32_t mesh() const {return mesh_;}
+
+  virtual void set_mesh(int32_t mesh);
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  std::vector<double> importance_;
+
+  int32_t mesh_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_IMPORTANCE_H

include/openmc/particle_data.h

@@ -42,9 +42,11 @@ struct SourceSite {
   Position r;
   Direction u;
   double E;
+  double E_parent;
   double time {0.0};
   double wgt {1.0};
   int delayed_group {0};
+  int fission_nuclide;
   int surf_id {SURFACE_NONE};
   ParticleType particle;
 
@@ -498,8 +500,15 @@ class ParticleData : public GeometryState {
 
   double E_;
   double E_last_;
+  double E_parent_; //!< energy of parent neutron in eV
+  double E_pprod_;
   int g_ {0};
   int g_last_;
+
+  // Other birth data
+  int fission_nuclide_;  //!< this particle was born as a result of this nuclide fissioning
+  int pprod_nuclide_;
+  // a double for fission cross section at birth? if so, I need to also add it to the bank...
 
   double wgt_ {1.0};
   double wgt_born_ {1.0};
@@ -543,6 +552,11 @@ class ParticleData : public GeometryState {
   int64_t current_work_;
 
   vector<double> flux_derivs_;
+
+  std::vector<std::vector<double>> cumulative_sensitivities_;  // for sensitivities for this particle
+  std::vector<std::vector<double>> cum_sens_;
+  std::vector<std::vector<double>> pprod_sens_;
+  double neutron_flux_last_ = 0.0;
 
   vector<FilterMatch> filter_matches_;
 
@@ -605,6 +619,10 @@ class ParticleData : public GeometryState {
   const double& E() const { return E_; }
   double& E_last() { return E_last_; }
   const double& E_last() const { return E_last_; }
+  double& E_parent() { return E_parent_; }              // for sensitivity analysis
+  const double& E_parent() const { return E_parent_; }  // for SA
+  double& E_pprod() { return E_pprod_; }
+  const double& E_pprod() const { return E_pprod_; }
   int& g() { return g_; }
   const int& g() const { return g_; }
   int& g_last() { return g_last_; }
@@ -716,6 +734,21 @@ class ParticleData : public GeometryState {
   // Used in tally derivatives
   double& flux_derivs(int i) { return flux_derivs_[i]; }
   const double& flux_derivs(int i) const { return flux_derivs_[i]; }
+
+  // Used in sensitivity analysis
+  std::vector<double>& cumulative_sensitivities(int i) { return cumulative_sensitivities_[i]; }
+  const std::vector<double>& cumulative_sensitivities(int i) const { return cumulative_sensitivities_[i]; }
+  std::vector<double>& cum_sens(int i) { return cum_sens_[i]; }
+  const std::vector<double>& cum_sens(int i) const { return cum_sens_[i]; }
+  std::vector<double>& pprod_sens(int i) { return pprod_sens_[i]; }
+  const std::vector<double>& pprod_sens(int i) const { return pprod_sens_[i]; }
+  double neutron_flux_last() const { return neutron_flux_last_; }
+
+  int& fission_nuclide() { return fission_nuclide_; }
+  const int& fission_nuclide() const { return fission_nuclide_; }
+
+  int& pprod_nuclide() { return pprod_nuclide_; }
+  const int& pprod_nuclide() const { return pprod_nuclide_; }
 
   // Matches of tallies
   decltype(filter_matches_)& filter_matches() { return filter_matches_; }
@@ -784,6 +817,79 @@ class ParticleData : public GeometryState {
       d = 0;
     }
   }
+
+  void resize_flux_derivs(int newSize)
+  {
+    flux_derivs_.resize(newSize, 0.0);
+  }
+
+  // Resize and initialize sensitivity vectors
+  void resize_init_cumulative_sensitivities(int newSize)
+  {
+    cumulative_sensitivities_.resize(newSize, {0.0});
+  }
+
+  void resize_init_cumulative_sensitivities_vec(int indx, int newSize)
+  {
+    cumulative_sensitivities_[indx].resize(newSize, 0.0);
+  }
+
+  void resize_init_cum_sens(int newSize)
+  {
+    cum_sens_.resize(newSize, {0.0});
+  }
+
+  void resize_init_cum_sens_vec(int indx, int newSize)
+  {
+    cum_sens_[indx].resize(newSize, 0.0);
+  }
+
+  void resize_init_pprod_sens(int newSize)
+  {
+    pprod_sens_.resize(newSize, {0.0});
+  }
+
+  void resize_init_pprod_sens_vec(int indx, int newSize)
+  {
+    pprod_sens_[indx].resize(newSize, 0.0);
+  }
+
+  void resize_alloc_filter_matches(int newSize)
+  {
+    filter_matches_.resize(newSize);
+    for (auto &match: filter_matches_){
+      match.bins_.resize(1, 0.0);
+      match.weights_.resize(1, 0.0);
+      match.i_bin_ = 0;
+      // bins_present_ default is false from header
+    }
+  }
+
+  void initialize_cumulative_sensitivities()
+  {
+    for (auto& it : cumulative_sensitivities_){
+      std::fill(it.begin(), it.end(), 0.0);
+    }
+  }
+
+  void initialize_cum_sens()
+  {
+    for (auto& it : cum_sens_){
+      std::fill(it.begin(), it.end(), 0.0);
+    }
+  }
+
+  void initialize_pprod_sens()
+  {
+    for (auto& it : pprod_sens_){
+      std::fill(it.begin(), it.end(), 0.0);
+    }
+  }
+
+  void set_neutron_flux_last(double v) 
+  { 
+    neutron_flux_last_ = v;
+  }
 };
 
 } // namespace openmc

include/openmc/simulation.h

@@ -73,6 +73,11 @@ void initialize_generation();
 //! Full initialization of a particle history
 void initialize_history(Particle& p, int64_t index_source);
 
+//! Helper function for initialize_history() that is called independently elsewhere
+void initialize_history_partial(Particle& p);
+
+void initialize_history_in_part(Particle& p);
+
 //! Finalize a batch
 //!
 //! Handles synchronization and accumulation of tallies, calculation of Shannon

include/openmc/tallies/filter.h

@@ -33,6 +33,7 @@ enum class FilterType {
   MATERIALFROM,
   MESH,
   MESHBORN,
+  IMPORTANCE,
   MESH_MATERIAL,
   MESH_SURFACE,
   MU,

include/openmc/tallies/filter_importance.h

@@ -0,0 +1,58 @@
+#ifndef OPENMC_TALLIES_FILTER_IMPORTANCE_H
+#define OPENMC_TALLIES_FILTER_IMPORTANCE_H
+
+#include <cstdint>
+
+#include <gsl/gsl-lite.hpp>
+#include "openmc/vector.h"
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Multiplies each tally by the importance corresponding to the mesh index 
+//==============================================================================
+
+class ImportanceFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~ImportanceFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override {return "importance";}
+  FilterType type() const override { return FilterType::IMPORTANCE; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const std::vector<double>& importance() const { return importance_; }
+  void set_importance(gsl::span<const double> importance);
+
+  virtual int32_t mesh() const {return mesh_;}
+
+  virtual void set_mesh(int32_t mesh);
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  std::vector<double> importance_;
+
+  int32_t mesh_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_IMPORTANCE_H