From e3cd406a9b8ac6acf3b3f04419f132cead52459b Mon Sep 17 00:00:00 2001 From: Christopher Fichtlscherer <29277544+cfichtlscherer@users.noreply.github.com> Date: Mon, 10 Jul 2023 05:50:08 +0200 Subject: [PATCH] Pulse height tally for photons (#2452) Co-authored-by: Paul Romano --- docs/source/usersguide/tallies.rst | 5 + include/openmc/constants.h | 5 +- include/openmc/particle.h | 4 + include/openmc/particle_data.h | 3 + include/openmc/tallies/tally.h | 4 + include/openmc/tallies/tally_scoring.h | 9 +- openmc/lib/tally.py | 4 +- src/output.cpp | 1 + src/particle.cpp | 64 ++++++ src/particle_data.cpp | 5 + src/reaction.cpp | 1 + src/simulation.cpp | 3 + src/tallies/tally.cpp | 72 ++++++- src/tallies/tally_scoring.cpp | 76 +++++++ .../regression_tests/pulse_height/__init__.py | 0 .../pulse_height/inputs_true.dat | 42 ++++ .../pulse_height/results_true.dat | 201 ++++++++++++++++++ tests/regression_tests/pulse_height/test.py | 51 +++++ 18 files changed, 538 insertions(+), 12 deletions(-) create mode 100644 tests/regression_tests/pulse_height/__init__.py create mode 100644 tests/regression_tests/pulse_height/inputs_true.dat create mode 100644 tests/regression_tests/pulse_height/results_true.dat create mode 100644 tests/regression_tests/pulse_height/test.py diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 97ca319a2..28ce6795f 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -317,6 +317,11 @@ The following tables show all valid scores: | |particle. This corresponds to MT=444 produced by | | |NJOY's HEATR module. | +----------------------+---------------------------------------------------+ + |pulse-height |The energy deposited by an entire photon's history | + | |(including its progeny). Units are eV per source | + | |particle. Note that this score can only be combined| + | |with a cell filter and an energy filter. | + +----------------------+---------------------------------------------------+ .. _usersguide_tally_normalization: diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 6a77396b7..2b9eaca29 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -279,7 +279,7 @@ enum class MgxsType { enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE }; -enum class TallyType { VOLUME, MESH_SURFACE, SURFACE }; +enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT }; enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION }; @@ -307,7 +307,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 }; // Global tally parameters diff --git a/include/openmc/particle.h b/include/openmc/particle.h index e3f23fd5b..376c9abf3 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -66,6 +66,10 @@ public: void event_revive_from_secondary(); void event_death(); + //! pulse-height recording + void pht_collision_energy(); + void pht_secondary_particles(); + //! Cross a surface and handle boundary conditions void cross_surface(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index c4f959666..6d86d528e 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -313,6 +313,8 @@ private: vector nu_bank_; // bank of most recently fissioned particles + vector pht_storage_; // interim pulse-height results + // Global tally accumulators double keff_tally_absorption_ {0.0}; double keff_tally_collision_ {0.0}; @@ -444,6 +446,7 @@ public: decltype(tracks_)& tracks() { return tracks_; } decltype(nu_bank_)& nu_bank() { return nu_bank_; } NuBank& nu_bank(int i) { return nu_bank_[i]; } + vector& pht_storage() { return pht_storage_; } double& keff_tally_absorption() { return keff_tally_absorption_; } double& keff_tally_collision() { return keff_tally_collision_; } diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 303b1ebb1..2f0cddcf0 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -169,8 +169,10 @@ public: // 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 energy_filter_ {C_NONE}; int energyout_filter_ {C_NONE}; int delayedgroup_filter_ {C_NONE}; + int cell_filter_ {C_NONE}; vector triggers_; @@ -206,6 +208,8 @@ extern vector active_tracklength_tallies; extern vector active_collision_tallies; extern vector active_meshsurf_tallies; extern vector active_surface_tallies; +extern vector active_pulse_height_tallies; +extern vector pulse_height_cells; } // namespace model namespace simulation { diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 85d0fc58f..28f1f1622 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -94,9 +94,16 @@ void score_tracklength_tally(Particle& p, double distance); //! Score surface or mesh-surface tallies for particle currents. // //! \param p The particle being tracked -//! \param tallies A vector of tallies to score to +//! \param tallies A vector of the indices of the tallies to score to void score_surface_tally(Particle& p, const vector& tallies); +//! Score the pulse-height tally +//! This is triggered at the end of every particle history +// +//! \param p The particle being tracked +//! \param tallies A vector of the indices of the tallies to score to +void score_pulse_height_tally(Particle& p, const vector& tallies); + } // namespace openmc #endif // OPENMC_TALLIES_TALLY_SCORING_H diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index f327692e7..417d084a6 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -104,13 +104,13 @@ _SCORES = { -5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission', -9: 'current', -10: 'events', -11: 'delayed-nu-fission', -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt', - -15: 'fission-q-recoverable', -16: 'decay-rate' + -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'pulse-height' } _ESTIMATORS = { 0: 'analog', 1: 'tracklength', 2: 'collision' } _TALLY_TYPES = { - 0: 'volume', 1: 'mesh-surface', 2: 'surface' + 0: 'volume', 1: 'mesh-surface', 2: 'surface', 3: 'pulse-height' } diff --git a/src/output.cpp b/src/output.cpp index 2b0a4024c..17e158d3f 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -594,6 +594,7 @@ const std::unordered_map 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 index f455ab834..8dd4a12e5 100644 --- 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" @@ -290,6 +291,11 @@ void Particle::event_collide() } } + if (!model::active_pulse_height_tallies.empty() && + type() == ParticleType::photon) { + pht_collision_energy(); + } + // Reset banked weight during collision n_bank() = 0; n_bank_second() = 0; @@ -354,6 +360,25 @@ void Particle::event_revive_from_secondary() secondary_bank().pop_back(); n_event() = 0; + // Subtract secondary particle energy from interim pulse-height results + if (!model::active_pulse_height_tallies.empty() && + this->type() == ParticleType::photon) { + // Since the birth cell of the particle has not been set we + // have to determine it before the energy of the secondary particle can be + // removed from the pulse-height of this cell. + 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; + } + pht_secondary_particles(); + } + // Enter new particle in particle track file if (write_track()) add_particle_track(*this); @@ -387,6 +412,10 @@ void Particle::event_death() keff_tally_tracklength() = 0.0; keff_tally_leakage() = 0.0; + if (!model::active_pulse_height_tallies.empty()) { + score_pulse_height_tally(*this, model::active_pulse_height_tallies); + } + // Record the number of progeny created by this particle. // This data will be used to efficiently sort the fission bank. if (settings::run_mode == RunMode::EIGENVALUE) { @@ -395,6 +424,41 @@ void Particle::event_death() } } +void Particle::pht_collision_energy() +{ + // Adds the energy particles lose in a collision to the pulse-height + + // determine index of cell in pulse_height_cells + auto it = std::find(model::pulse_height_cells.begin(), + model::pulse_height_cells.end(), coord(n_coord() - 1).cell); + + if (it != model::pulse_height_cells.end()) { + int index = std::distance(model::pulse_height_cells.begin(), it); + pht_storage()[index] += E_last() - E(); + + // If the energy of the particle is below the cutoff, it will not be sampled + // so its energy is added to the pulse-height in the cell + int photon = static_cast(ParticleType::photon); + if (E() < settings::energy_cutoff[photon]) { + pht_storage()[index] += E(); + } + } +} + +void Particle::pht_secondary_particles() +{ + // Removes the energy of secondary produced particles from the pulse-height + + // determine index of cell in pulse_height_cells + auto it = std::find(model::pulse_height_cells.begin(), + model::pulse_height_cells.end(), cell_born()); + + if (it != model::pulse_height_cells.end()) { + int index = std::distance(model::pulse_height_cells.begin(), it); + pht_storage()[index] -= E(); + } +} + void Particle::cross_surface() { int i_surface = std::abs(surface()); diff --git a/src/particle_data.cpp b/src/particle_data.cpp index a2be72084..0fc7202c5 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -53,6 +53,11 @@ ParticleData::ParticleData() // Create microscopic cross section caches neutron_xs_.resize(data::nuclides.size()); photon_xs_.resize(data::elements.size()); + + // Creates the pulse-height storage for the particle + if (!model::pulse_height_cells.empty()) { + pht_storage_.resize(model::pulse_height_cells.size(), 0.0); + } } TrackState ParticleData::get_track_state() const diff --git a/src/reaction.cpp b/src/reaction.cpp index ec9d32f9a..0643cb08b 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -171,6 +171,7 @@ std::unordered_map 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 index 310629936..a7aea6942 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -550,6 +550,9 @@ void initialize_history(Particle& p, int64_t index_source) // Reset weight window ratio p.ww_factor() = 0.0; + // Reset pulse_height_storage + std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0); + // set random number seed int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) * settings::n_particles + diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index c9c977ea9..7fd444a58 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -3,6 +3,7 @@ #include "openmc/array.h" #include "openmc/capi.h" #include "openmc/constants.h" +#include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/mesh.h" @@ -18,6 +19,7 @@ #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_cell.h" +#include "openmc/tallies/filter_cellborn.h" #include "openmc/tallies/filter_cellfrom.h" #include "openmc/tallies/filter_collision.h" #include "openmc/tallies/filter_delayedgroup.h" @@ -55,6 +57,8 @@ vector active_tracklength_tallies; vector active_collision_tallies; vector active_meshsurf_tallies; vector active_surface_tallies; +vector active_pulse_height_tallies; +vector pulse_height_cells; } // namespace model namespace simulation { @@ -175,6 +179,16 @@ Tally::Tally(pugi::xml_node node) } // Check if tally is compatible with particle type + if (!settings::photon_transport) { + for (int score : scores_) { + switch (score) { + case SCORE_PULSE_HEIGHT: + fatal_error( + "For pulse-height tallies, photon transport needs to be activated."); + break; + } + } + } if (settings::photon_transport) { if (particle_filter_index == C_NONE) { for (int score : scores_) { @@ -409,10 +423,14 @@ void Tally::add_filter(Filter* filter) return; // Keep track of indices for special filters - if (dynamic_cast(filter)) { + if (filter->type() == FilterType::ENERGY_OUT) { energyout_filter_ = filters_.size(); - } else if (dynamic_cast(filter)) { + } else if (filter->type() == FilterType::DELAYED_GROUP) { delayedgroup_filter_ = filters_.size(); + } else if (filter->type() == FilterType::CELL) { + cell_filter_ = filters_.size(); + } else if (filter->type() == FilterType::ENERGY) { + energy_filter_ = filters_.size(); } filters_.push_back(filter_idx); } @@ -454,17 +472,23 @@ void Tally::set_scores(const vector& scores) bool cellfrom_present = false; bool surface_present = false; bool meshsurface_present = false; + bool non_cell_energy_present = false; for (auto i_filt : filters_) { const auto* filt {model::tally_filters[i_filt].get()}; - if (dynamic_cast(filt)) { + // Checking for only cell and energy filters for pulse-height tally + if (!(filt->type() == FilterType::CELL || + filt->type() == FilterType::ENERGY)) { + non_cell_energy_present = true; + } + if (filt->type() == FilterType::LEGENDRE) { legendre_present = true; - } else if (dynamic_cast(filt)) { + } else if (filt->type() == FilterType::CELLFROM) { cellfrom_present = true; - } else if (dynamic_cast(filt)) { + } else if (filt->type() == FilterType::CELL) { cell_present = true; - } else if (dynamic_cast(filt)) { + } else if (filt->type() == FilterType::SURFACE) { surface_present = true; - } else if (dynamic_cast(filt)) { + } else if (filt->type() == FilterType::MESH_SURFACE) { meshsurface_present = true; } } @@ -536,6 +560,31 @@ void Tally::set_scores(const vector& scores) if (settings::photon_transport) estimator_ = TallyEstimator::COLLISION; break; + + case SCORE_PULSE_HEIGHT: + if (non_cell_energy_present) { + fatal_error("Pulse-height tallies are not compatible with filters " + "other than CellFilter and EnergyFilter"); + } + type_ = TallyType::PULSE_HEIGHT; + + // Collecting indices of all cells covered by the filters in the pulse + // height tally in global variable pulse_height_cells + for (const auto& i_filt : filters_) { + auto cell_filter = + dynamic_cast(model::tally_filters[i_filt].get()); + if (cell_filter) { + const auto& cells = cell_filter->cells(); + for (int i = 0; i < cell_filter->n_bins(); i++) { + int cell_index = cells[i]; + if (!contains(model::pulse_height_cells, cell_index)) { + model::pulse_height_cells.push_back(cell_index); + } + } + } + } + + break; } scores_.push_back(score); @@ -953,6 +1002,7 @@ void setup_active_tallies() model::active_collision_tallies.clear(); model::active_meshsurf_tallies.clear(); model::active_surface_tallies.clear(); + model::active_pulse_height_tallies.clear(); for (auto i = 0; i < model::tallies.size(); ++i) { const auto& tally {*model::tallies[i]}; @@ -980,6 +1030,11 @@ void setup_active_tallies() case TallyType::SURFACE: model::active_surface_tallies.push_back(i); + break; + + case TallyType::PULSE_HEIGHT: + model::active_pulse_height_tallies.push_back(i); + break; } } } @@ -1001,6 +1056,7 @@ void free_memory_tally() model::active_collision_tallies.clear(); model::active_meshsurf_tallies.clear(); model::active_surface_tallies.clear(); + model::active_pulse_height_tallies.clear(); model::tally_map.clear(); } @@ -1122,6 +1178,8 @@ extern "C" int openmc_tally_set_type(int32_t index, const char* type) model::tallies[index]->type_ = TallyType::MESH_SURFACE; } else if (strcmp(type, "surface") == 0) { model::tallies[index]->type_ = TallyType::SURFACE; + } else if (strcmp(type, "pulse-height") == 0) { + model::tallies[index]->type_ = TallyType::PULSE_HEIGHT; } else { set_errmsg(fmt::format("Unknown tally type: {}", type)); return OPENMC_E_INVALID_ARGUMENT; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 3e58fec05..47fe4d1be 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -15,6 +15,7 @@ #include "openmc/string_utils.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_cell.h" #include "openmc/tallies/filter_delayedgroup.h" #include "openmc/tallies/filter_energy.h" @@ -2498,4 +2499,79 @@ void score_surface_tally(Particle& p, const vector& tallies) match.bins_present_ = false; } +void score_pulse_height_tally(Particle& p, const vector& tallies) +{ + // The pulse height tally in OpenMC hijacks the logic of CellFilter and + // EnergyFilter to score specific quantities related to particle pulse height. + // This is achieved by setting the pulse-height cell of the tally to the cell + // of the particle being scored, and the energy to the particle's last + // recorded energy (E_last()). After the tally is scored, the values are reset + // to ensure proper accounting and avoid interference with subsequent + // calculations or tallies. + + // Save original cell/energy information + int orig_n_coord = p.n_coord(); + int orig_cell = p.coord(0).cell; + double orig_E_last = p.E_last(); + + for (auto i_tally : tallies) { + auto& tally {*model::tallies[i_tally]}; + + // Determine all CellFilter in the tally + for (const auto& filter : tally.filters()) { + auto cell_filter = + dynamic_cast(model::tally_filters[filter].get()); + if (cell_filter != nullptr) { + + const auto& cells = cell_filter->cells(); + // Loop over all cells in the CellFilter + for (auto cell_index = 0; cell_index < cells.size(); ++cell_index) { + int cell_id = cells[cell_index]; + + // Temporarily change cell of particle + p.n_coord() = 1; + p.coord(0).cell = cell_id; + + // Determine index of cell in model::pulse_height_cells + auto it = std::find(model::pulse_height_cells.begin(), + model::pulse_height_cells.end(), cell_id); + int index = std::distance(model::pulse_height_cells.begin(), it); + + // Temporarily change energy of particle to pulse-height value + p.E_last() = p.pht_storage()[index]; + + // Initialize an iterator over valid filter bin combinations. If + // there are no valid combinations, use a continue statement to ensure + // we skip the assume_separate break below. + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true, &p.filter_matches()); + if (filter_iter == end) + continue; + + // Loop over filter bins. + for (; filter_iter != end; ++filter_iter) { + auto filter_index = filter_iter.index_; + auto filter_weight = filter_iter.weight_; + + // Loop over scores. + for (auto score_index = 0; score_index < tally.scores_.size(); + ++score_index) { +#pragma omp atomic + tally.results_(filter_index, score_index, TallyResult::VALUE) += + filter_weight; + } + } + + // Reset all the filter matches for the next tally event. + for (auto& match : p.filter_matches()) + match.bins_present_ = false; + } + } + } + // Restore cell/energy + p.n_coord() = orig_n_coord; + p.coord(0).cell = orig_cell; + p.E_last() = orig_E_last; + } +} } // namespace openmc diff --git a/tests/regression_tests/pulse_height/__init__.py b/tests/regression_tests/pulse_height/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/pulse_height/inputs_true.dat b/tests/regression_tests/pulse_height/inputs_true.dat new file mode 100644 index 000000000..544ddb097 --- /dev/null +++ b/tests/regression_tests/pulse_height/inputs_true.dat @@ -0,0 +1,42 @@ + + + + + + + + + + + + + + + + + fixed source + 100 + 5 + + + 0.0 0.0 0.0 + + + 1000000.0 1.0 + + + true + + + + 1 + + + 0.0 10000.0 20000.0 30000.0 40000.0 50000.0 60000.0 70000.0 80000.0 90000.0 100000.0 110000.0 120000.0 130000.0 140000.0 150000.0 160000.0 170000.0 180000.0 190000.0 200000.0 210000.0 220000.0 230000.0 240000.0 250000.0 260000.0 270000.0 280000.0 290000.0 300000.0 310000.0 320000.0 330000.0 340000.0 350000.0 360000.0 370000.0 380000.0 390000.0 400000.0 410000.0 420000.0 430000.0 440000.0 450000.0 460000.0 470000.0 480000.0 490000.0 500000.0 510000.0 520000.0 530000.0 540000.0 550000.0 560000.0 570000.0 580000.0 590000.0 600000.0 610000.0 620000.0 630000.0 640000.0 650000.0 660000.0 670000.0 680000.0 690000.0 700000.0 710000.0 720000.0 730000.0 740000.0 750000.0 760000.0 770000.0 780000.0 790000.0 800000.0 810000.0 820000.0 830000.0 840000.0 850000.0 860000.0 870000.0 880000.0 890000.0 900000.0 910000.0 920000.0 930000.0 940000.0 950000.0 960000.0 970000.0 980000.0 990000.0 1000000.0 + + + 1 2 + pulse-height + + + diff --git a/tests/regression_tests/pulse_height/results_true.dat b/tests/regression_tests/pulse_height/results_true.dat new file mode 100644 index 000000000..c57e8ff1c --- /dev/null +++ b/tests/regression_tests/pulse_height/results_true.dat @@ -0,0 +1,201 @@ +tally 1: +4.140000E+00 +3.443000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +2.000000E-02 +4.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.000000E-02 +3.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.000000E-02 +6.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +4.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +2.000000E-02 +2.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +3.000000E-02 +5.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-02 +2.000000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.000000E-01 +8.600000E-03 diff --git a/tests/regression_tests/pulse_height/test.py b/tests/regression_tests/pulse_height/test.py new file mode 100644 index 000000000..c6459a866 --- /dev/null +++ b/tests/regression_tests/pulse_height/test.py @@ -0,0 +1,51 @@ +import numpy as np +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def sphere_model(): + + model = openmc.model.Model() + + # Define materials + NaI = openmc.Material() + NaI.set_density('g/cc', 3.7) + NaI.add_element('Na', 1.0) + NaI.add_element('I', 1.0) + + model.materials = openmc.Materials([NaI]) + + # Define geometry: two spheres in each other + s1 = openmc.Sphere(r=1) + s2 = openmc.Sphere(r=2, boundary_type='vacuum') + inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1) + outer_sphere = openmc.Cell(name='outer sphere', region=+s1 & -s2) + model.geometry = openmc.Geometry([inner_sphere, outer_sphere]) + + # Define settings + model.settings.run_mode = 'fixed source' + model.settings.batches = 5 + model.settings.particles = 100 + model.settings.photon_transport = True + model.settings.source = openmc.Source(space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1e6],[1]), + particle='photon') + + # Define tallies + tally = openmc.Tally(name="pht tally") + tally.scores = ['pulse-height'] + cell_filter = openmc.CellFilter(inner_sphere) + energy_filter = openmc.EnergyFilter(np.linspace(0, 1_000_000, 101)) + tally.filters = [cell_filter, energy_filter] + model.tallies = [tally] + + return model + + + +def test_pulse_height(sphere_model): + harness = PyAPITestHarness('statepoint.5.h5', sphere_model) + harness.main()