openmc-dev · cfichtlscherer · Sep 14, 2021 · Sep 17, 2021 · Sep 21, 2021 · Sep 21, 2021
diff --git a/examples/custom_source/CMakeLists.txt b/examples/custom_source/CMakeLists.txt
@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
 project(openmc_sources CXX)
 add_library(source SHARED source_ring.cpp)
-find_package(OpenMC REQUIRED)
+find_package(OpenMC REQUIRED HINTS "/home/cpf/openmc/build/bin/openmc")
 if (OpenMC_FOUND)
   message(STATUS "Found OpenMC: ${OpenMC_DIR}")
 endif()

diff --git a/include/openmc/constants.h b/include/openmc/constants.h
@@ -306,7 +306,8 @@ enum TallyScore {
   SCORE_INVERSE_VELOCITY = -13,   // flux-weighted inverse velocity
   SCORE_FISS_Q_PROMPT = -14,      // prompt fission Q-value
   SCORE_FISS_Q_RECOV = -15,       // recoverable fission Q-value
-  SCORE_DECAY_RATE = -16          // delayed neutron precursor decay rate
+  SCORE_DECAY_RATE = -16,         // delayed neutron precursor decay rate
+  SCORE_PULSE_HEIGHT = -17        // pulse-height of photons in a cell
 };
 
 // Global tally parameters

diff --git a/include/openmc/particle.h b/include/openmc/particle.h
@@ -63,6 +63,16 @@ class Particle : public ParticleData {
   void event_revive_from_secondary();
   void event_death();
 
+  //! needed for determining the value of the pulse-height tally
+  //
+  //! sums up the energy a photon looses after in a cell
+  void energy_delivered_in_cell();
+  //! adds the last energy of the particle if it dies.
+  void killed_particle_energy_delivered();
+  //! removes the energy from created secondary particles to avoid
+  //! double counting
+  void remove_energy_of_secondary();
+
   //! Cross a surface and handle boundary conditions
   void cross_surface();
 

diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
@@ -302,6 +302,8 @@ class ParticleData {
   // Secondary particle bank
   vector<SourceSite> secondary_bank_;
 
+  std::vector<double> pht_storage_; // array for the pulse-height tally
+
   int64_t current_work_; // current work index
 
   vector<double> flux_derivs_; // for derivatives for this particle
@@ -435,6 +437,7 @@ class ParticleData {
   SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
   decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
   int64_t& current_work() { return current_work_; }
+  std::vector<double>& pht_storage() { return pht_storage_; }
   const int64_t& current_work() const { return current_work_; }
   double& flux_derivs(int i) { return flux_derivs_[i]; }
   const double& flux_derivs(int i) const { return flux_derivs_[i]; }

diff --git a/include/openmc/settings.h b/include/openmc/settings.h
@@ -123,6 +123,8 @@ extern int trigger_batch_interval; //!< Batch interval for triggers
 extern "C" int verbosity;          //!< How verbose to make output
 extern double weight_cutoff;       //!< Weight cutoff for Russian roulette
 extern double weight_survive;      //!< Survival weight after Russian roulette
+extern std::vector<double> bins_pht; //! Used for the pulse-height tally
+extern int cell_pht;                 //! Used for the pulse-height tally
 } // namespace settings
 
 //==============================================================================

diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
@@ -47,6 +47,8 @@ extern const RegularMesh* ufs_mesh;
 extern vector<double> k_generation;
 extern vector<int64_t> work_index;
 
+extern std::vector<double> bins_pht;
+extern int cell_pht;
 } // namespace simulation
 
 //==============================================================================

diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
@@ -97,6 +97,11 @@ void score_tracklength_tally(Particle& p, double distance);
 //! \param tallies A vector of tallies to score to
 void score_surface_tally(Particle& p, const vector<int>& tallies);
 
+//! Score the pulse-height tally
+//
+//! \param p The particle being tracked
+void score_pht_tally(Particle& p);
+
 } // namespace openmc
 
 #endif // OPENMC_TALLIES_TALLY_SCORING_H
diff --git a/src/output.cpp b/src/output.cpp
@@ -639,6 +639,7 @@ const std::unordered_map<int, const char*> score_names = {
   {SCORE_FISS_Q_PROMPT, "Prompt fission power"},
   {SCORE_FISS_Q_RECOV, "Recoverable fission power"},
   {SCORE_CURRENT, "Current"},
+  {SCORE_PULSE_HEIGHT, "pulse-height"},
 };
 
 //! Create an ASCII output file showing all tally results.

diff --git a/src/particle.cpp b/src/particle.cpp
@@ -17,6 +17,7 @@
 #include "openmc/message_passing.h"
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
+#include "openmc/particle_data.h"
 #include "openmc/photon.h"
 #include "openmc/physics.h"
 #include "openmc/physics_mg.h"
@@ -85,6 +86,7 @@ void Particle::create_secondary(
 void Particle::from_source(const SourceSite* src)
 {
   // Reset some attributes
+
   clear();
   surface() = 0;
   cell_born() = C_NONE;
@@ -286,6 +288,11 @@ void Particle::event_collide()
     }
   }
 
+  // TODO: The assumption behind the modelling of the pulse-height tally is that
+  // particles lose energy only in collisions, for a more comprehensive modeling
+  // of charged particles this has to be changed. 
+  if (type() == ParticleType::photon){energy_delivered_in_cell();} // score the pht only for photons
+
   // Reset banked weight during collision
   n_bank() = 0;
   n_bank_second() = 0;
@@ -315,7 +322,6 @@ void Particle::event_collide()
       coord(j + 1).u = coord(j).u;
     }
   }
-
   // Score flux derivative accumulators for differential tallies.
   if (!model::active_tallies.empty())
     score_collision_derivative(*this);
@@ -350,6 +356,9 @@ void Particle::event_revive_from_secondary()
     secondary_bank().pop_back();
     n_event() = 0;
 
+    // remove energy of created secondary particles from pht value
+    remove_energy_of_secondary();
+
     // Enter new particle in particle track file
     if (write_track())
       add_particle_track(*this);
@@ -389,6 +398,41 @@ void Particle::event_death()
     int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
     simulation::progeny_per_particle[offset] = n_progeny();
   }
+killed_particle_energy_delivered();
+}
+
+void Particle::energy_delivered_in_cell()
+{
+  // Adds the energy the particle loses during the collision.
+  pht_storage()[coord(n_coord() - 1).cell] += E_last() - E();
+}
+
+void Particle::killed_particle_energy_delivered()
+{
+  // Adds the energy of killed particles.
+  pht_storage()[coord(n_coord() - 1).cell] += E();
+}
+
+void Particle::remove_energy_of_secondary()
+{
+
+  // determine where the particle was born
+  if (coord(n_coord() - 1).cell == C_NONE) {
+    if (!exhaustive_find_cell(*this)) {
+      mark_as_lost(
+        "Could not find the cell containing particle " + std::to_string(id()));
+      return;
+    }
+
+    // Set birth cell attribute
+    if (cell_born() == C_NONE)
+      cell_born() = coord(n_coord() - 1).cell;
+  }
+
+  // settings::energy_cutoff[1] is the cutoff energy of photons
+  // Remove the energy of created secondary particles from the pht_storage array
+  if (type() == ParticleType::photon)
+    pht_storage()[cell_born()] -= E();
 }
 
 void Particle::cross_surface()

diff --git a/src/particle_data.cpp b/src/particle_data.cpp
@@ -53,6 +53,8 @@ ParticleData::ParticleData()
   // Create microscopic cross section caches
   neutron_xs_.resize(data::nuclides.size());
   photon_xs_.resize(data::elements.size());
+
+  pht_storage_.resize(model::cells.size(), 0.0);
 }
 
 TrackState ParticleData::get_track_state() const

diff --git a/src/physics.cpp b/src/physics.cpp
@@ -64,7 +64,6 @@ void collision(Particle& p)
     sample_positron_reaction(p);
     break;
   }
-
   // Kill particle if energy falls below cutoff
   int type = static_cast<int>(p.type());
   if (p.E() < settings::energy_cutoff[type]) {

diff --git a/src/reaction.cpp b/src/reaction.cpp
@@ -171,6 +171,7 @@ std::unordered_map<int, std::string> REACTION_NAME_MAP {
   {SCORE_INVERSE_VELOCITY, "inverse-velocity"},
   {SCORE_FISS_Q_PROMPT, "fission-q-prompt"},
   {SCORE_FISS_Q_RECOV, "fission-q-recoverable"},
+  {SCORE_PULSE_HEIGHT, "pulse-height"},
   // Normal ENDF-based reactions
   {TOTAL_XS, "(n,total)"},
   {ELASTIC, "(n,elastic)"},

diff --git a/src/simulation.cpp b/src/simulation.cpp
@@ -21,6 +21,7 @@
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/tallies/tally_scoring.h"
 #include "openmc/tallies/trigger.h"
 #include "openmc/timer.h"
 #include "openmc/track_output.h"
@@ -289,6 +290,9 @@ const RegularMesh* ufs_mesh {nullptr};
 vector<double> k_generation;
 vector<int64_t> work_index;
 
+std::vector<double> bins_pht;
+int cell_pht;
+
 } // namespace simulation
 
 //==============================================================================
@@ -709,8 +713,15 @@ void transport_history_based_single_particle(Particle& p)
       p.event_collide();
     }
     p.event_revive_from_secondary();
-    if (!p.alive())
+    if (!p.alive()) {
+      for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
+        const auto& tally {*model::tallies[i_tally]};
+        if (tally.scores_[0] == -17) {
+          score_pht_tally(p);
+        }
+      }
       break;
+    }
   }
   p.event_death();
 }

diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp
@@ -5,6 +5,7 @@
 #include "openmc/capi.h"
 #include "openmc/cell.h"
 #include "openmc/error.h"
+#include "openmc/simulation.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {

diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp
@@ -7,6 +7,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/search.h"
 #include "openmc/settings.h"
+#include "openmc/simulation.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -18,6 +19,7 @@ namespace openmc {
 void EnergyFilter::from_xml(pugi::xml_node node)
 {
   auto bins = get_node_array<double>(node, "bins");
+  simulation::bins_pht = bins; //! store them as a global variable for the pht
   this->set_bins(bins);
 }
 

diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
@@ -30,6 +30,10 @@
 #include "openmc/tallies/filter_surface.h"
 #include "openmc/xml_interface.h"
 
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 #include "xtensor/xadapt.hpp"
 #include "xtensor/xbuilder.hpp" // for empty_like
 #include "xtensor/xview.hpp"
@@ -253,6 +257,17 @@ Tally::Tally(pugi::xml_node node)
     }
   }
 
+  // Currently the pulse-height tally is only working in a selected type of
+  // settings Here these settings are checked again
+  bool exists = std::find(std::begin(scores_), std::end(scores_),
+                  TallyScore::SCORE_PULSE_HEIGHT) != std::end(scores_);
+  if (exists && scores_.size() > 1) {
+    fatal_error("Error: The Pulse-Height Tally can currently not be used in combination with other tallys.");
+  }
+  if (exists && !settings::photon_transport) {
+      fatal_error("Error: The Pulse-Height Tally works only with photon transport True.");
+  }
+
   // If settings.xml trigger is turned on, create tally triggers
   if (settings::trigger_on) {
     this->init_triggers(node);

diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
@@ -17,6 +17,7 @@
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/filter_delayedgroup.h"
 #include "openmc/tallies/filter_energy.h"
+#include "openmc/xml_interface.h"
 
 #include <string>
 
@@ -2468,4 +2469,26 @@ void score_surface_tally(Particle& p, const vector<int>& tallies)
     match.bins_present_ = false;
 }
 
+void score_pht_tally(Particle& p)
+{
+  for (auto i_tally : model::active_tracklength_tallies) {
+    auto& tally {*model::tallies[i_tally]};
+
+    auto filter_iter = FilterBinIter(tally, p);
+    auto score_index = 0;
+
+    auto score = p.pht_storage()[simulation::cell_pht];
+    //auto score = p.pht_storage()[7]; 
+    //auto score = p.pht_storage()[0]; 
+    // the -1 in the next line is for the index shift
+    auto scoring_it = std::upper_bound(simulation::bins_pht.begin(), simulation::bins_pht.end(), score) - simulation::bins_pht.begin() - 1;    
+#pragma omp atomic    
+    tally.results_(scoring_it, score_index, TallyResult::VALUE) += 1;
+  }
+
+  // Reset all the filter matches for the next tally event.
+  for (auto& match : p.filter_matches())
+    match.bins_present_ = false;
+}
+
 } // namespace openmc
diff --git a/tests/regression_tests/cpp_driver/CMakeLists.txt b/tests/regression_tests/cpp_driver/CMakeLists.txt
@@ -0,0 +1,6 @@
+
+cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
+project(openmc_cpp_driver CXX)
+add_executable(cpp_driver driver.cpp)
+find_package(OpenMC REQUIRED HINTS /home/cpf/openmc/build)
+target_link_libraries(cpp_driver OpenMC::libopenmc)
diff --git a/tests/regression_tests/pulse_height_tally/__init__.py b/tests/regression_tests/pulse_height_tally/__init__.py
diff --git a/tests/regression_tests/pulse_height_tally/inputs_true.dat b/tests/regression_tests/pulse_height_tally/inputs_true.dat
@@ -0,0 +1,55 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="1" name="detector" region="-1 2 -3" universe="1" />
+  <cell id="2" material="2" name="surrounding" region="~(-1 2 -3) -4" universe="1" />
+  <cell fill="1" id="3" name="root cell" region="-4" universe="2" />
+  <surface coeffs="0 0 2.54" id="1" name="cylinder" type="z-cylinder" />
+  <surface coeffs="0" id="2" name="bottom" type="z-plane" />
+  <surface coeffs="5" id="3" name="top" type="z-plane" />
+  <surface boundary="vacuum" coeffs="0 0 0 10" id="4" name="outer sphere" type="sphere" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material id="1" name="sodium iodide detector material">
+    <density units="g/cc" value="3.67" />
+    <nuclide ao="1.0" name="Na23" />
+    <nuclide ao="1.0" name="I127" />
+  </material>
+  <material id="2" name="air">
+    <density units="g/cc" value="0.001225" />
+    <nuclide ao="3.985348" name="N14" />
+    <nuclide ao="0.014652" name="N15" />
+    <nuclide ao="1" name="O16" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>1000</particles>
+  <batches>1</batches>
+  <inactive>0</inactive>
+  <source particle="photon" strength="1.0">
+    <space type="point">
+      <parameters>0.0 0.0 -1.0</parameters>
+    </space>
+    <angle reference_uvw="0.0 0.0 1.0" type="monodirectional" />
+    <energy type="discrete">
+      <parameters>662000.0 1.0</parameters>
+    </energy>
+  </source>
+  <photon_transport>true</photon_transport>
+  <seed>1</seed>
+</settings>
+<?xml version='1.0' encoding='utf-8'?>
+<tallies>
+  <filter id="1" type="cell">
+    <bins>1</bins>
+  </filter>
+  <filter id="2" type="energy">
+    <bins>0.0 20000.0 40000.0 60000.0 80000.0 100000.0 120000.0 140000.0 160000.0 180000.0 200000.0 220000.0 240000.0 260000.0 280000.0 300000.0 320000.0 340000.0 360000.0 380000.0 400000.0 420000.0 440000.0 460000.0 480000.0 500000.0 520000.0 540000.0 560000.0 580000.0 600000.0 620000.0 640000.0 660000.0 680000.0 700000.0 720000.0 740000.0 760000.0 780000.0 800000.0 820000.0 840000.0 860000.0 880000.0 900000.0 920000.0 940000.0 960000.0 980000.0 1000000.0</bins>
+  </filter>
+  <tally id="1" name="pulse-height tally">
+    <filters>1 2</filters>
+    <scores>pulse-height</scores>
+  </tally>
+</tallies>