From df0ec91a2352aa2efaab3c1cc75ea04af4d86a59 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 22 Sep 2025 15:31:49 +0300 Subject: [PATCH 01/18] wip --- include/openmc/constants.h | 10 +- include/openmc/particle.h | 8 ++ include/openmc/physics.h | 2 +- include/openmc/position.h | 11 ++ include/openmc/simulation.h | 3 + include/openmc/tallies/filter.h | 1 + include/openmc/tallies/filter_point.h | 55 +++++++++ include/openmc/tallies/tally.h | 4 + include/openmc/tallies/tally_scoring.h | 12 ++ src/particle.cpp | 121 +++++++++++++++--- src/physics.cpp | 30 ++++- src/simulation.cpp | 24 ++++ src/source.cpp | 5 + src/tallies/filter_point.cpp | 74 +++++++++++ src/tallies/tally.cpp | 22 ++++ src/tallies/tally_scoring.cpp | 162 +++++++++++++++++++++++++ 16 files changed, 522 insertions(+), 22 deletions(-) create mode 100644 include/openmc/tallies/filter_point.h create mode 100644 src/tallies/filter_point.cpp diff --git a/include/openmc/constants.h b/include/openmc/constants.h index df13da3707d..a6dfe32685e 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -288,9 +288,15 @@ enum class MgxsType { enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE }; -enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT }; +enum class TallyType { + VOLUME, + MESH_SURFACE, + SURFACE, + PULSE_HEIGHT, + NEXT_EVENT +}; -enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION }; +enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION, POINT }; enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB }; diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 0f37719b94f..14face2d059 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -37,6 +37,8 @@ class Particle : public ParticleData { //========================================================================== // Methods + double mass() const; + double speed() const; //! create a secondary particle @@ -64,6 +66,8 @@ class Particle : public ParticleData { //! \param src Source site data void from_source(const SourceSite* src); + void initialize_pseudoparticle(Particle& p, Direction u, double E); + // Coarse-grained particle events void event_calculate_xs(); void event_advance(); @@ -72,6 +76,10 @@ class Particle : public ParticleData { void event_revive_from_secondary(); void event_death(); + // Pseudo-particles events + void event_advance_pseudo(double& total_distance, double& mfp); + void event_cross_surface_pseudo(); + //! pulse-height recording void pht_collision_energy(); void pht_secondary_particles(); diff --git a/include/openmc/physics.h b/include/openmc/physics.h index f62f43a02f2..eb81bb66a31 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -92,7 +92,7 @@ void sample_fission_neutron( //! handles all reactions with a single secondary neutron (other than fission), //! i.e. level scattering, (n,np), (n,na), etc. -void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p); +void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle& p); void sample_secondary_photons(Particle& p, int i_nuclide); diff --git a/include/openmc/position.h b/include/openmc/position.h index 5d291d26b95..92346e27ad4 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -204,6 +204,17 @@ inline Position operator/(double a, Position b) return b /= a; } +inline bool operator<(Position a, Position b) +{ + if (a[0] != b[0]) + return (a[0] < b[0]); + if (a[1] != b[1]) + return (a[1] < b[1]); + if (a[2] != b[2]) + return (a[2] < b[2]); + return false; +} + inline Position Position::reflect(Position n) const { const double projection = n.dot(*this); diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 3e4e24e1d06..0e99ea5b16a 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -91,6 +91,9 @@ void broadcast_results(); void free_memory_simulation(); +//! Simulate a single pseudoparticle history +void transport_pseudoparticle(Particle& p, double total_distance, double& mfp); + //! Simulate a single particle history (and all generated secondary particles, //! if enabled), from birth to death void transport_history_based_single_particle(Particle& p); diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 65098597a5d..e900d35af9d 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -39,6 +39,7 @@ enum class FilterType { MUSURFACE, PARENT_NUCLIDE, PARTICLE, + POINT, POLAR, SPHERICAL_HARMONICS, SPATIAL_LEGENDRE, diff --git a/include/openmc/tallies/filter_point.h b/include/openmc/tallies/filter_point.h new file mode 100644 index 00000000000..909f0baa27d --- /dev/null +++ b/include/openmc/tallies/filter_point.h @@ -0,0 +1,55 @@ +#ifndef OPENMC_TALLIES_FILTER_POINT_H +#define OPENMC_TALLIES_FILTER_POINT_H + +#include "openmc/position.h" +#include "openmc/span.h" +#include "openmc/tallies/filter.h" +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +//! Bins tally by point detectors +//============================================================================== + +class PointFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~PointFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "point"; } + FilterType type() const override { return FilterType::POINT; } + + 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 vector>& detectors() const + { + return detectors_; + } + + void set_detectors(span> detectors); + +private: + //---------------------------------------------------------------------------- + // Data members + + vector> detectors_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_PARTICLE_H diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 3beeb9d5ac5..e5ce8b98948 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -3,6 +3,7 @@ #include "openmc/constants.h" #include "openmc/memory.h" // for unique_ptr +#include "openmc/position.h" #include "openmc/span.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/trigger.h" @@ -12,6 +13,7 @@ #include "xtensor/xfixed.hpp" #include "xtensor/xtensor.hpp" +#include #include #include @@ -205,11 +207,13 @@ extern vector active_analog_tallies; extern vector active_tracklength_tallies; extern vector active_timed_tracklength_tallies; extern vector active_collision_tallies; +extern vector active_point_tallies; extern vector active_meshsurf_tallies; extern vector active_surface_tallies; extern vector active_pulse_height_tallies; extern vector pulse_height_cells; extern vector time_grid; +extern std::set active_point_detectors; } // namespace model diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index c3ab779e6a1..2ebce2ec35d 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -1,9 +1,11 @@ #ifndef OPENMC_TALLIES_TALLY_SCORING_H #define OPENMC_TALLIES_TALLY_SCORING_H +#include "openmc/nuclide.h" #include "openmc/particle.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" +#include "openmc/thermal.h" namespace openmc { @@ -114,6 +116,16 @@ void score_surface_tally(Particle& p, const vector& tallies); //! \param tallies A vector of the indices of the tallies to score to void score_pulse_height_tally(Particle& p, const vector& tallies); +void score_point_tally(Particle& p, int i_nuclide, const Reaction& rx, + int i_product, Direction* v_t); + +void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, + const NuclideMicroXS& micro); + +void score_point_tally(SourceSite& site, int source_index); + +void score_pseudoparticle_tally( + Particle& p, double distance, double mfp, double pdf); } // namespace openmc #endif // OPENMC_TALLIES_TALLY_SCORING_H diff --git a/src/particle.cpp b/src/particle.cpp index f5ad45d80fe..7cc71eaab5d 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -42,23 +42,25 @@ namespace openmc { // Particle implementation //============================================================================== +double Particle::mass() const +{ + // Determine mass in eV/c^2 + switch (this->type()) { + case ParticleType::neutron: + return MASS_NEUTRON_EV; + case ParticleType::photon: + return 0.0; + case ParticleType::electron: + case ParticleType::positron: + return MASS_ELECTRON_EV; + } + UNREACHABLE(); +} + double Particle::speed() const { if (settings::run_CE) { - // Determine mass in eV/c^2 - double mass; - switch (this->type()) { - case ParticleType::neutron: - mass = MASS_NEUTRON_EV; - break; - case ParticleType::photon: - mass = 0.0; - break; - case ParticleType::electron: - case ParticleType::positron: - mass = MASS_ELECTRON_EV; - break; - } + double mass = this->mass(); // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<E() * (this->E() + 2 * mass)) / (this->E() + mass); @@ -152,6 +154,42 @@ void Particle::from_source(const SourceSite* src) } } +void Particle::initialize_pseudoparticle( + Particle& p, Direction u_new, double E_new) +{ + // Reset some attributes + clear(); + surface() = SURFACE_NONE; + cell_born() = C_NONE; + material() = C_NONE; + n_collision() = 0; + fission() = false; + zero_flux_derivs(); + lifetime() = 0.0; + + // Copy attributes from particle + type() = p.type(); + wgt() = p.wgt(); + wgt_last() = p.wgt(); + r() = p.r(); + u() = u_new; + r_born() = p.r(); + r_last_current() = p.r(); + r_last() = p.r(); + u_last() = p.u(); + if (settings::run_CE) { + E() = E_new; + g() = 0; + } else { + g() = static_cast(E_new); + g_last() = static_cast(E_new); + E() = data::mg.energy_bin_avg_[g()]; + } + E_last() = E(); + time() = p.time(); + time_last() = p.time(); +} + void Particle::event_calculate_xs() { // Set the random number stream @@ -501,6 +539,61 @@ void Particle::event_death() } } +void Particle::event_advance_pseudo(double& total_distance, double& mfp) +{ + // Find the distance to the nearest boundary + boundary() = distance_to_boundary(*this); + + double speed = this->speed(); + double time_cutoff = settings::time_cutoff[static_cast(type())]; + double distance_cutoff = + (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; + + if (distance_cutoff < boundary().distance()) { + wgt() = 0.0; + return; + } + + // Select smaller of the two distances + double distance = std::min({boundary().distance(), total_distance}); + + // Advance particle in space and time + this->move_distance(distance); + double dt = distance / speed; + this->time() += dt; + this->lifetime() += dt; + mfp += distance * macro_xs().total; + total_distance -= distance; +} + +void Particle::event_cross_surface_pseudo() +{ + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + + // Set surface that particle is on and adjust coordinate levels + surface() = boundary().surface(); + n_coord() = boundary().coord_level(); + + if (boundary().lattice_translation()[0] != 0 || + boundary().lattice_translation()[1] != 0 || + boundary().lattice_translation()[2] != 0) { + // Particle crosses lattice boundary + + bool verbose = settings::verbosity >= 10 || trace(); + cross_lattice(*this, boundary(), verbose); + event() = TallyEvent::LATTICE; + } else { + // Particle crosses surface + const auto& surf {model::surfaces[surface_index()].get()}; + this->cross_surface(*surf); + event() = TallyEvent::SURFACE; + } +} + void Particle::pht_collision_energy() { // Adds the energy particles lose in a collision to the pulse-height diff --git a/src/physics.cpp b/src/physics.cpp index f667fd586d7..d24cfee0984 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -25,6 +25,7 @@ #include "openmc/simulation.h" #include "openmc/string_utils.h" #include "openmc/tallies/tally.h" +#include "openmc/tallies/tally_scoring.h" #include "openmc/thermal.h" #include "openmc/weight_windows.h" @@ -732,7 +733,7 @@ void scatter(Particle& p, int i_nuclide) // Perform collision physics for inelastic scattering const auto& rx {nuc->reactions_[i]}; - inelastic_scatter(*nuc, *rx, p); + inelastic_scatter(i_nuclide, *rx, p); p.event_mt() = rx->mt_; } @@ -778,6 +779,10 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p) // Find speed of neutron in CM vel = v_n.norm(); + if (!model::active_point_tallies.empty()) { + score_point_tally(p, i_nuclide, rx, 0, &v_t); + } + // Sample scattering angle, checking if angle distribution is present (assume // isotropic otherwise) double mu_cm; @@ -825,8 +830,12 @@ void sab_scatter(int i_nuclide, int i_sab, Particle& p) // Sample energy and angle double E_out; - data::thermal_scatt[i_sab]->data_[i_temp].sample( - micro, p.E(), &E_out, &p.mu(), p.current_seed()); + auto& sab = data::thermal_scatt[i_sab]->data_[i_temp]; + sab.sample(micro, p.E(), &E_out, &p.mu(), p.current_seed()); + + if (!model::active_point_tallies.empty()) { + score_point_tally(p, i_nuclide, sab, micro); + } // Set energy to outgoing, change direction of particle p.E() = E_out; @@ -1091,6 +1100,10 @@ void sample_fission_neutron( site->delayed_group = 0; } + if (!model::active_point_tallies.empty()) { + score_point_tally(p, i_nuclide, rx, site->delayed_group, nullptr); + } + // sample from prompt neutron energy distribution int n_sample = 0; double mu; @@ -1116,8 +1129,11 @@ void sample_fission_neutron( site->u = rotate_angle(p.u(), mu, nullptr, seed); } -void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) +void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle& p) { + // Get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide]}; + // copy energy of neutron double E_in = p.E(); @@ -1126,13 +1142,17 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) double mu; rx.products_[0].sample(E_in, E, mu, p.current_seed()); + if (!model::active_point_tallies.empty()) { + score_point_tally(p, i_nuclide, rx, 0, nullptr); + } + // if scattering system is in center-of-mass, transfer cosine of scattering // angle and outgoing energy from CM to LAB if (rx.scatter_in_cm_) { double E_cm = E; // determine outgoing energy in lab - double A = nuc.awr_; + double A = nuc->awr_; E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) / ((A + 1.0) * (A + 1.0)); diff --git a/src/simulation.cpp b/src/simulation.cpp index b9986f47375..a0fd11a69a0 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -22,6 +22,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" @@ -801,6 +802,29 @@ void transport_history_based_single_particle(Particle& p) p.event_death(); } +void transport_pseudoparticle(Particle& p, double total_distance, double& mfp) +{ + if (!p.alive()) + return; + + double distance = total_distance; + bool cross_surface = false; + while (distance > 0.0) { + p.event_calculate_xs(); + if (p.alive()) { + cross_surface = (distance > p.boundary().distance()); + p.event_advance_pseudo(distance, mfp); + } else { + return; + } + if (p.alive() && cross_surface) { + p.event_cross_surface_pseudo(); + } + if (!p.alive()) + return; + } +} + void transport_history_based() { #pragma omp parallel for schedule(runtime) diff --git a/src/source.cpp b/src/source.cpp index 12323f7bd7e..fb2388704be 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -33,6 +33,7 @@ #include "openmc/simulation.h" #include "openmc/state_point.h" #include "openmc/string_utils.h" +#include "openmc/tallies/tally_scoring.h" #include "openmc/xml_interface.h" namespace openmc { @@ -650,6 +651,10 @@ SourceSite sample_external_source(uint64_t* seed) site.E = data::mg.num_energy_groups_ - site.E - 1.; } + if (!model::active_point_tallies.empty()) { + score_point_tally(site, i); + } + return site; } diff --git a/src/tallies/filter_point.cpp b/src/tallies/filter_point.cpp new file mode 100644 index 00000000000..8d09da6f3d0 --- /dev/null +++ b/src/tallies/filter_point.cpp @@ -0,0 +1,74 @@ +#include "openmc/tallies/filter_point.h" + +#include + +#include "openmc/math_functions.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void PointFilter::from_xml(pugi::xml_node node) +{ + auto bins = get_node_array(node, "bins"); + + // Convert to vector of detectors + vector> detectors; + size_t n = bins.size() / 4; + for (int i = 0; i < n, ++i) { + Position pos {bins[4 * i], bins[4 * i + 1], bins[4 * i + 2]}; + detectors.push_back(std::make_pair(pos, bins[4 * i + 3])); + } + this->set_detectors(detectors); +} + +void PointFilter::set_detectors(span> detectors) +{ + // Clear existing detectors + detectors_.clear(); + detectors_.reserve(detectors.size()); + + // Set detectors and number of bins + for (auto d : detectors) { + detectors_.push_back(d); + } + n_bins_ = detectors_.size(); +} + +void PointFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + for (auto i = 0; i < detectors_.size(); i++) { + auto [pos, r] = detectors_[i]; + if ((p.r() - pos).norm() < FP_COINCIDENT) { + match.bins_.push_back(i); + double weight; + if ((p.r_last() - pos).norm() < r) { + weight = 3.0 * exprel(-r * p.macro_xs().total); + } else { + weight = p.wgt() / p.wgt_last(); + } + match.weights_.push_back(weight); + } + } +} + +void PointFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + vector detectors; + for (auto [pos, r] : detectors_) { + detectors.push_back(pos[0]); + detectors.push_back(pos[1]); + detectors.push_back(pos[2]); + detectors.push_back(r); + } + write_dataset(filter_group, "bins", particles); +} + +std::string PointFilter::text_label(int bin) const +{ + auto [pos, r] = detectors_.at(bin); + return fmt::format("Point: {} {}", pos, r); +} + +} // namespace openmc diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ae0bffe6eb0..2e8cf53f47a 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -11,6 +11,7 @@ #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" #include "openmc/particle.h" +#include "openmc/position.h" #include "openmc/reaction.h" #include "openmc/reaction_product.h" #include "openmc/settings.h" @@ -61,11 +62,13 @@ vector active_analog_tallies; vector active_tracklength_tallies; vector active_timed_tracklength_tallies; vector active_collision_tallies; +vector active_point_tallies; vector active_meshsurf_tallies; vector active_surface_tallies; vector active_pulse_height_tallies; vector pulse_height_cells; vector time_grid; +std::set active_point_detectors; } // namespace model namespace simulation { @@ -1104,6 +1107,14 @@ void add_to_time_grid(vector grid) model::time_grid.swap(merged); } +//! Add new points to the global time grid +// +//! \param detector Position of new point detector to add +void add_point_detector(Position& detector) +{ + model::active_point_detectors.insert(detector); +} + void setup_active_tallies() { model::active_tallies.clear(); @@ -1111,10 +1122,12 @@ void setup_active_tallies() model::active_tracklength_tallies.clear(); model::active_timed_tracklength_tallies.clear(); model::active_collision_tallies.clear(); + model::active_point_tallies.clear(); model::active_meshsurf_tallies.clear(); model::active_surface_tallies.clear(); model::active_pulse_height_tallies.clear(); model::time_grid.clear(); + model::active_point_detectors.clear(); for (auto i = 0; i < model::tallies.size(); ++i) { const auto& tally {*model::tallies[i]}; @@ -1155,6 +1168,13 @@ void setup_active_tallies() case TallyType::PULSE_HEIGHT: model::active_pulse_height_tallies.push_back(i); break; + + case TallyType::NEXT_EVENT: + switch (tally.estimator_) { + case TallyEstimator::POINT: + model::active_point_tallies.push_back(i); + break; + } } } } @@ -1175,10 +1195,12 @@ void free_memory_tally() model::active_tracklength_tallies.clear(); model::active_timed_tracklength_tallies.clear(); model::active_collision_tallies.clear(); + model::active_point_tallies.clear(); model::active_meshsurf_tallies.clear(); model::active_surface_tallies.clear(); model::active_pulse_height_tallies.clear(); model::time_grid.clear(); + model::active_point_detectors.clear(); model::tally_map.clear(); } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 0df80a2398f..14945689e5a 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -13,6 +13,7 @@ #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/source.h" #include "openmc/string_utils.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" @@ -24,6 +25,10 @@ namespace openmc { +namespace simulation { +thread_local Particle pseudoparticle; +} + //============================================================================== // FilterBinIter implementation //============================================================================== @@ -2729,4 +2734,161 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) p.E_last() = orig_E_last; } } + +void score_pseudoparticle_tally( + Particle& p, double distance, double mfp, double pdf) +{ + double attenuation = std::exp(-mfp); + double flux = pdf * p.wgt() / (2 * PI * distance * distance); + + // Save the attenuation for point filter handling + p.wgt_last() = p.wgt(); + p.wgt() *= attenuation; + + for (auto i_tally : model::active_point_tallies) { + const Tally& tally {*model::tallies[i_tally]}; + + // 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 nuclide bins. + for (auto i = 0; i < tally.nuclides_.size(); ++i) { + auto i_nuclide = tally.nuclides_[i]; + + // Tally this event in the present nuclide bin if that bin represents + // the event nuclide or the total material. Note that the atomic + // density argument for score_general is not used for analog tallies. + if (i_nuclide == p.event_nuclide() || i_nuclide == -1) + score_general_ce_analog(p, i_tally, i * tally.scores_.size(), + filter_index, filter_weight, i_nuclide, -1.0, flux); + } + } + + // If the user has specified that we can assume all tallies are spatially + // separate, this implies that once a tally has been scored to, we needn't + // check the others. This cuts down on overhead when there are many + // tallies specified + if (settings::assume_separate) + break; + } + + // Reset all the filter matches for the next tally event. + for (auto& match : p.filter_matches()) + match.bins_present_ = false; +} + +void score_point_tally(Particle& p, const Nuclide* nuc, const Reaction& rx, + int i_product, Direction* v_t) +{ + double awr = nuc->awr_; + double E_in = p.E(); + double vel = std::sqrt(E_in); + Direction v_n = vel * p.u(); + Direction v_cm = (v_n + awr * (v_t != nullptr ? *v_t : 0)) / (awr + 1.0); + Direction u_cm = (v_n - v_cm); + u_cm /= u_cm.norm(); + double mfp; + + for (auto& det : model::active_point_detectors) { + + auto u = (det - p.r()); + double distance = u.norm(); + u /= distance; + + double mu; + if (rx.scatter_in_cm_) { + mu = u.dot(u_cm); + } else { + mu = u.dot(p.u()); + } + double E_out; + double pdf = rx.products_[i_product].sample_energy_and_pdf( + p.E(), mu, E_out, p.current_seed()); + if (rx.scatter_in_cm_) { + double E_cm = E_out; + Direction v_out = std::sqrt(E_cm) * u_cm; + v_out += v_cm; + E_out = std::pow(v_out.norm(), 2); + pdf *= std::sqrt(E_out / E_cm) / + (1 - mu / (awr + 1) * std::sqrt(E_in / E_out)); + } + simulation::pseudoparticle.initialize_pseudoparticle(p, u, E_out); + mfp = 0.0; + transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); + if (!simulation::pseudoparticle.alive()) + continue; + score_pseudoparticle_tally(simulation::pseudoparticle, distance, mfp, pdf); + } +} + +void score_point_tally(Particle& p, const Nuclide* nuc, const ThermalData& sab, + const NuclideMicroXS& micro) +{ + double awr = nuc->awr_; + double E_in = p.E(); + double vel = std::sqrt(E_in); + Direction v_n = vel * p.u(); + Direction v_cm = v_n / (awr + 1.0); + Direction u_cm = (v_n - v_cm); + u_cm /= u_cm.norm(); + double mfp; + + for (auto& det : model::active_point_detectors) { + + auto u = (det - p.r()); + double distance = u.norm(); + u /= distance; + + double mu = u.dot(p.u()); + double E_out; + double pdf = + sab.sample_energy_and_pdf(micro, p.E(), mu, E_out, p.current_seed()); + simulation::pseudoparticle.initialize_pseudoparticle(p, u, E_out); + mfp = 0.0; + transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); + if (!simulation::pseudoparticle.alive()) + continue; + score_pseudoparticle_tally(simulation::pseudoparticle, distance, mfp, pdf); + } +} + +void score_point_tally(SourceSite& site, int source_index) +{ + double E_out = site.E; + Position r = site.r; + auto src_ = model::external_sources[source_index].get(); + auto src = dynamic_cast(src_); + if (!src) + fatal_error("Only independent source is valid for point detectors."); + auto angle = src->angle(); + + double mfp; + + for (auto& det : model::active_point_detectors) { + + auto u = (det - r); + double distance = u.norm(); + u /= distance; + + double pdf = angle->evaluate(u); + simulation::pseudoparticle.from_source(&site); + simulation::pseudoparticle.u() = u; + mfp = 0.0; + transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); + if (!simulation::pseudoparticle.alive()) + continue; + score_pseudoparticle_tally(simulation::pseudoparticle, distance, mfp, pdf); + } +} + } // namespace openmc From 312f327a63b63286403bbdf8c0a3f534b0646e8d Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 11:32:31 +0200 Subject: [PATCH 02/18] merge get_pdf updates --- include/openmc/angle_energy.h | 2 + include/openmc/chain.h | 2 + include/openmc/distribution_angle.h | 1 + include/openmc/distribution_multi.h | 10 ++ include/openmc/reaction_product.h | 5 + include/openmc/secondary_correlated.h | 3 + include/openmc/secondary_kalbach.h | 4 + include/openmc/secondary_nbody.h | 3 + include/openmc/secondary_thermal.h | 63 +++++++++ include/openmc/secondary_uncorrelated.h | 2 + include/openmc/thermal.h | 6 +- src/chain.cpp | 7 + src/distribution_angle.cpp | 15 ++ src/distribution_multi.cpp | 14 ++ src/reaction_product.cpp | 19 ++- src/secondary_correlated.cpp | 21 ++- src/secondary_kalbach.cpp | 22 ++- src/secondary_nbody.cpp | 28 ++-- src/secondary_thermal.cpp | 178 +++++++++++++++++++++--- src/secondary_uncorrelated.cpp | 13 ++ src/thermal.cpp | 20 ++- 21 files changed, 396 insertions(+), 42 deletions(-) diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h index ac931b1b533..35e5b68b64b 100644 --- a/include/openmc/angle_energy.h +++ b/include/openmc/angle_energy.h @@ -16,6 +16,8 @@ class AngleEnergy { public: virtual void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const = 0; + virtual double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const = 0; virtual ~AngleEnergy() = default; }; diff --git a/include/openmc/chain.h b/include/openmc/chain.h index a3bc6f3a364..c2a45f75adc 100644 --- a/include/openmc/chain.h +++ b/include/openmc/chain.h @@ -70,6 +70,8 @@ class DecayPhotonAngleEnergy : public AngleEnergy { //! \param[inout] seed Pseudorandom seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: const Distribution* photon_energy_; diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h index efd4e58425b..a7e021085f0 100644 --- a/include/openmc/distribution_angle.h +++ b/include/openmc/distribution_angle.h @@ -25,6 +25,7 @@ class AngleDistribution { //! \param[inout] seed pseudorandom number seed pointer //! \return Cosine of the angle in the range [-1,1] double sample(double E, uint64_t* seed) const; + double evaluate(double E, double mu) const; //! Determine whether angle distribution is empty //! \return Whether distribution is empty diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h index 7b9c2abf8ce..d53e5ed560e 100644 --- a/include/openmc/distribution_multi.h +++ b/include/openmc/distribution_multi.h @@ -6,6 +6,7 @@ #include "pugixml.hpp" #include "openmc/distribution.h" +#include "openmc/error.h" #include "openmc/position.h" namespace openmc { @@ -29,6 +30,11 @@ class UnitSphereDistribution { //! \return (sampled Direction, sample weight) virtual std::pair sample(uint64_t* seed) const = 0; + virtual double evaluate(Direction u) const + { + fatal_error("evaluate not available for this UnitSphereDistribution type"); + } + Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction }; @@ -52,6 +58,8 @@ class PolarAzimuthal : public UnitSphereDistribution { //! \return (sampled Direction, value of the PDF at this Direction) std::pair sample_as_bias(uint64_t* seed) const; + double evaluate(Direction u) const override; + // Observing pointers Distribution* mu() const { return mu_.get(); } Distribution* phi() const { return phi_.get(); } @@ -87,6 +95,8 @@ class Isotropic : public UnitSphereDistribution { //! \return (sampled direction, sample weight) std::pair sample(uint64_t* seed) const override; + double evaluate(Direction u) const override; + // Set or get bias distribution void set_bias(std::unique_ptr bias) { diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index 4fbbc1b626a..98e5910db9b 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -49,6 +49,11 @@ class ReactionProduct { //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const; + AngleEnergy& sample_dist(double E_in, uint64_t* seed) const; + + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const; + ParticleType particle_; //!< Particle type EmissionMode emission_mode_; //!< Emission mode double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s] diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h index 6905c38e369..e7c933bb614 100644 --- a/include/openmc/secondary_correlated.h +++ b/include/openmc/secondary_correlated.h @@ -40,6 +40,9 @@ class CorrelatedAngleEnergy : public AngleEnergy { //! \param[inout] seed Pseudorandom seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + Distribution& sample_dist(double E_in, double& E_out, uint64_t* seed) const; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; // energy property vector& energy() { return energy_; } diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index 83806d35248..0168ee13ced 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -31,6 +31,10 @@ class KalbachMann : public AngleEnergy { //! \param[inout] seed Pseudorandom seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + void sample_params(double E_in, double& E_out, double& km_a, double& km_r, + uint64_t* seed) const; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: //! Outgoing energy/angle at a single incoming energy diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h index efb4fd75ba1..dad896cd740 100644 --- a/include/openmc/secondary_nbody.h +++ b/include/openmc/secondary_nbody.h @@ -27,6 +27,9 @@ class NBodyPhaseSpace : public AngleEnergy { //! \param[inout] seed Pseudorandom seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + double sample_energy(double E_in, uint64_t* seed) const; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: int n_bodies_; //!< Number of particles distributed diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 5b18902afbb..b7de6ab4e61 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -6,6 +6,7 @@ #include "openmc/angle_energy.h" #include "openmc/endf.h" +#include "openmc/search.h" #include "openmc/secondary_correlated.h" #include "openmc/vector.h" @@ -32,6 +33,8 @@ class CoherentElasticAE : public AngleEnergy { //! \param[inout] seed Pseudorandom seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section @@ -55,6 +58,8 @@ class IncoherentElasticAE : public AngleEnergy { //! \param[inout] seed Pseudorandom number seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: double debye_waller_; @@ -80,6 +85,8 @@ class IncoherentElasticAEDiscrete : public AngleEnergy { //! \param[inout] seed Pseudorandom number seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: const vector& energy_; //!< Energies at which cosines are tabulated @@ -106,6 +113,9 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy { //! \param[inout] seed Pseudorandom number seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + void sample_params(double E_in, double& E_out, int& j, uint64_t* seed) const; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: const vector& energy_; //!< Incident energies @@ -134,6 +144,10 @@ class IncoherentInelasticAE : public AngleEnergy { //! \param[inout] seed Pseudorandom number seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + void sample_params(double E_in, double& E_out, double& f, int& l, int& j, + uint64_t* seed) const; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: //! Secondary energy/angle distribution @@ -169,6 +183,9 @@ class MixedElasticAE : public AngleEnergy { //! \param[inout] seed Pseudorandom number seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + const AngleEnergy& sample_dist(double E_in, uint64_t* seed) const; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; private: CoherentElasticAE coherent_dist_; //!< Coherent distribution @@ -178,6 +195,52 @@ class MixedElasticAE : public AngleEnergy { const Function1D& incoherent_xs_; //!< Polymorphic ref. to incoherent XS }; +struct DoubleVector { + double data; + const double& operator[](size_t index) const { return data; } +}; + +template +double get_pdf_discrete(const vector& mu, const T& w, double mu_0) +{ + // Make sure mu is in range [-1,1] + if (std::abs(mu_0) > 1.0) + mu_0 = std::copysign(1.0, mu_0); + double a0; + double a1; + double b0; + double b1; + int32_t ai = -1; + int32_t bi = -1; + if (mu_0 > mu[0]) { + ai = lower_bound_index(mu.begin(), mu.end(), mu_0); + a0 = mu[ai]; + a1 = (ai > 1) ? mu[ai - 1] : -1.0; + } else { + a0 = -1.0; + a1 = -1.0; + } + if (mu_0 < mu[mu.size() - 1]) { + bi = upper_bound_index(mu.begin(), mu.end(), mu_0); + b0 = mu[bi]; + b1 = (bi < mu.size() - 1) ? mu[bi + 1] : 1.0; + } else { + b0 = 1.0; + b1 = 1.0; + } + + // Calculate Delta_a and Delta_b + double delta_a = 0.5 * std::min(b0 - a0, a0 - a1); + double delta_b = 0.5 * std::min(b1 - b0, b0 - a0); + + if (mu_0 < a0 + delta_a) + return w[ai] / (2.0 * delta_a); + else if (mu_0 + delta_b < b0) + return w[bi] / (2.0 * delta_b); + else + return 0.0; +} + } // namespace openmc #endif // OPENMC_SECONDARY_THERMAL_H diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h index 3afa3d9ceb7..d971231b277 100644 --- a/include/openmc/secondary_uncorrelated.h +++ b/include/openmc/secondary_uncorrelated.h @@ -31,6 +31,8 @@ class UncorrelatedAngleEnergy : public AngleEnergy { //! \param[inout] seed Pseudorandom seed pointer void sample( double E_in, double& E_out, double& mu, uint64_t* seed) const override; + double sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const override; // Accessors AngleDistribution& angle() { return angle_; } diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index de0767d0af0..73c343ad6a7 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -51,7 +51,11 @@ class ThermalData { //! \param[out] mu Outgoing scattering angle cosine //! \param[inout] seed Pseudorandom seed pointer void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out, - double* mu, uint64_t* seed); + double* mu, uint64_t* seed) const; + AngleEnergy& sample_dist( + const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const; + double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in, + double mu, double& E_out, uint64_t* seed) const; private: struct Reaction { diff --git a/src/chain.cpp b/src/chain.cpp index a60ef12cd6c..85d9c119791 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -74,6 +74,13 @@ void DecayPhotonAngleEnergy::sample( mu = Uniform(-1., 1.).sample(seed).first; } +double DecayPhotonAngleEnergy::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + E_out = photon_energy_->sample(seed).first; + return 0.5; +} + //============================================================================== // Global variables //============================================================================== diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp index 50f1aca112b..66accac25ae 100644 --- a/src/distribution_angle.cpp +++ b/src/distribution_angle.cpp @@ -83,4 +83,19 @@ double AngleDistribution::sample(double E, uint64_t* seed) const return mu; } +double AngleDistribution::evaluate(double E, double mu) const +{ + // Find energy bin and calculate interpolation factor + int i; + double r; + get_energy_index(energy_, E, i, r); + + double pdf = 0.0; + if (r > 0.0) + pdf += r * distribution_[i + 1]->evaluate(mu); + if (r < 1.0) + pdf += (1.0 - r) * distribution_[i]->evaluate(mu); + return pdf; +} + } // namespace openmc diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp index 857e1c30b40..8ceb712ac08 100644 --- a/src/distribution_multi.cpp +++ b/src/distribution_multi.cpp @@ -44,6 +44,7 @@ UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node) fatal_error("Angular distribution reference direction must have " "three parameters specified."); u_ref_ = Direction(u_ref.data()); + u_ref_ /= u_ref_.norm(); } } @@ -65,6 +66,7 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node) fatal_error("Angular distribution reference v direction must have " "three parameters specified."); v_ref_ = Direction(v_ref.data()); + v_ref_ /= v_ref_.norm(); } w_ref_ = u_ref_.cross(v_ref_); if (check_for_node(node, "mu")) { @@ -116,6 +118,13 @@ std::pair PolarAzimuthal::sample_impl( weight}; } +double PolarAzimuthal::evaluate(Direction u) const +{ + double mu = u.dot(u_ref_); + double phi = std::acos(u.dot(v_ref_) / std::sqrt(1 - mu * mu)); + return mu_->evaluate(mu) * phi_->evaluate(phi); +} + //============================================================================== // Isotropic implementation //============================================================================== @@ -157,6 +166,11 @@ std::pair Isotropic::sample(uint64_t* seed) const } } +double Isotropic::evaluate(Direction u) const +{ + return 1.0 / (4.0 * PI); +} + //============================================================================== // Monodirectional implementation //============================================================================== diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 3ba2c0cfb05..398f72e7967 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -106,8 +106,7 @@ ReactionProduct::ReactionProduct(const ChainNuclide::Product& product) make_unique(chain_nuc->photon_energy())); } -void ReactionProduct::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const +AngleEnergy& ReactionProduct::sample_dist(double E_in, uint64_t* seed) const { auto n = applicability_.size(); if (n > 1) { @@ -119,14 +118,26 @@ void ReactionProduct::sample( // If i-th distribution is sampled, sample energy from the distribution if (c <= prob) { - distribution_[i]->sample(E_in, E_out, mu, seed); + return *distribution_[i]; break; } } } else { // If only one distribution is present, go ahead and sample it - distribution_[0]->sample(E_in, E_out, mu, seed); + return *distribution_[0]; } } +void ReactionProduct::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + sample_dist(E_in, seed).sample(E_in, E_out, mu, seed); +} + +double ReactionProduct::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed); +} + } // namespace openmc diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index e820419b484..dfde3f1e03b 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -153,9 +153,8 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) distribution_.push_back(std::move(d)); } // incoming energies } - -void CorrelatedAngleEnergy::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const +Distribution& CorrelatedAngleEnergy::sample_dist( + double E_in, double& E_out, uint64_t* seed) const { // Find energy bin and calculate interpolation factor int i; @@ -247,10 +246,22 @@ void CorrelatedAngleEnergy::sample( // Find correlated angular distribution for closest outgoing energy bin if (r1 - c_k < c_k1 - r1 || distribution_[l].interpolation == Interpolation::histogram) { - mu = distribution_[l].angle[k]->sample(seed).first; + return *distribution_[l].angle[k]; } else { - mu = distribution_[l].angle[k + 1]->sample(seed).first; + return *distribution_[l].angle[k + 1]; } } +void CorrelatedAngleEnergy::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + mu = sample_dist(E_in, E_out, seed).sample(seed).first; +} + +double CorrelatedAngleEnergy::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + return sample_dist(E_in, E_out, seed).evaluate(mu); +} + } // namespace openmc diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 6ac91e665b9..43a40c0cd4e 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -115,8 +115,8 @@ KalbachMann::KalbachMann(hid_t group) } // incoming energies } -void KalbachMann::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const +void KalbachMann::sample_params( + double E_in, double& E_out, double& km_a, double& km_r, uint64_t* seed) const { // Find energy bin and calculate interpolation factor int i; @@ -171,7 +171,6 @@ void KalbachMann::sample( double E_l_k = distribution_[l].e_out[k]; double p_l_k = distribution_[l].p[k]; - double km_r, km_a; if (distribution_[l].interpolation == Interpolation::histogram) { // Histogram interpolation if (p_l_k > 0.0 && k >= n_discrete) { @@ -217,6 +216,13 @@ void KalbachMann::sample( E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); } } +} + +void KalbachMann::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + double km_r, km_a; + sample_params(E_in, E_out, km_a, km_r, seed); // Sampled correlated angle from Kalbach-Mann parameters if (prn(seed) > km_r) { @@ -227,5 +233,15 @@ void KalbachMann::sample( mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a; } } +double KalbachMann::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + double km_r, km_a; + sample_params(E_in, E_out, km_a, km_r, seed); + + // https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle + return km_a / (2 * std::sinh(km_a)) * + (std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu)); +} } // namespace openmc diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp index da0bb81c471..72f0b0b92d0 100644 --- a/src/secondary_nbody.cpp +++ b/src/secondary_nbody.cpp @@ -22,13 +22,8 @@ NBodyPhaseSpace::NBodyPhaseSpace(hid_t group) read_attribute(group, "q_value", Q_); } -void NBodyPhaseSpace::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const +double NBodyPhaseSpace::sample_energy(double E_in, uint64_t* seed) const { - // By definition, the distribution of the angle is isotropic for an N-body - // phase space distribution - mu = uniform_distribution(-1., 1., seed); - // Determine E_max parameter double Ap = mass_ratio_; double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_); @@ -59,12 +54,29 @@ void NBodyPhaseSpace::sample( std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2); break; default: - throw std::runtime_error {"N-body phase space with >5 bodies."}; + fatal_error("N-body phase space with >5 bodies."); } // Now determine v and E_out double v = x / (x + y); - E_out = E_max * v; + return E_max * v; +} + +void NBodyPhaseSpace::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + // By definition, the distribution of the angle is isotropic for an N-body + // phase space distribution + mu = uniform_distribution(-1., 1., seed); + + E_out = sample_energy(E_in, seed); +} + +double NBodyPhaseSpace::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + E_out = sample_energy(E_in, seed); + return 0.5; } } // namespace openmc diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 030d398aabe..240bbf8cf36 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -4,6 +4,7 @@ #include "openmc/math_functions.h" #include "openmc/random_lcg.h" #include "openmc/search.h" +#include "openmc/vector.h" #include "xtensor/xview.hpp" @@ -12,6 +13,12 @@ namespace openmc { +double get_pdf_discrete(const vector& mu, double mu_0) +{ + DoubleVector w {1.0 / mu.size()}; + return get_pdf_discrete(mu, w, mu_0); +} + //============================================================================== // CoherentElasticAE implementation //============================================================================== @@ -42,6 +49,39 @@ void CoherentElasticAE::sample( mu = 1.0 - 2.0 * energies[k] / E_in; } +double CoherentElasticAE::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1) + + double pdf; + E_out = E_in; + const auto& energies {xs_.bragg_edges()}; + const auto& factors = xs_.factors(); + + if (E_in < energies.front() || E_in > energies.back()) { + return 0; + } + + const int n = upper_bound_index(energies.begin(), energies.end(), E_in); + + vector mu_vector; + mu_vector.reserve(n); + + std::transform(energies.rbegin() + n - 1, energies.rend(), + std::back_inserter(mu_vector), + [E_in](double Ei) { return 1 - 2 * Ei / E_in; }); + + vector weights; + weights.reserve(n); + + weights.emplace_back(factors[0] / factors[n]); + for (int i = 1; i <= n; ++i) { + weights.emplace_back((factors[i] - factors[i - 1]) / factors[n]); + } + return get_pdf_discrete(mu_vector, weights, mu); +} + //============================================================================== // IncoherentElasticAE implementation //============================================================================== @@ -54,12 +94,21 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group) void IncoherentElasticAE::sample( double E_in, double& E_out, double& mu, uint64_t* seed) const { + E_out = E_in; + // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4 double c = 2 * E_in * debye_waller_; mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0; - - // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4) +} +double IncoherentElasticAE::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ E_out = E_in; + + // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4 + double c = 2 * E_in * debye_waller_; + double A = c / (1 - std::exp(-2.0 * c)); // normalization factor + return A * std::exp(-c * (1 - mu)); } //============================================================================== @@ -116,6 +165,28 @@ void IncoherentElasticAEDiscrete::sample( E_out = E_in; } +double IncoherentElasticAEDiscrete::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + // Get index and interpolation factor for elastic grid + int i; + double f; + get_energy_index(energy_, E_in, i, f); + // Energy doesn't change in elastic scattering + E_out = E_in; + int n_mu = mu_out_.shape()[1]; + + std::vector mu_vector; + mu_vector.reserve(n_mu); + + for (int k = 0; k < n_mu; ++k) { + double mu_k = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k)); + mu_vector.push_back(mu_k); + } + + return get_pdf_discrete(mu_vector, mu); +} + //============================================================================== // IncoherentInelasticAEDiscrete implementation //============================================================================== @@ -129,8 +200,8 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete( read_dataset(group, "skewed", skewed_); } -void IncoherentInelasticAEDiscrete::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const +void IncoherentInelasticAEDiscrete::sample_params( + double E_in, double& E_out, int& j, uint64_t* seed) const { // Get index and interpolation factor for inelastic grid int i; @@ -144,7 +215,6 @@ void IncoherentInelasticAEDiscrete::sample( // for the second and second to last bins, relative to a normal bin // probability of 1). Otherwise, each bin is equally probable. - int j; int n = energy_out_.shape()[1]; if (!skewed_) { // All bins equally likely @@ -176,6 +246,18 @@ void IncoherentInelasticAEDiscrete::sample( // Outgoing energy E_out = (1 - f) * E_ij + f * E_i1j; +} + +void IncoherentInelasticAEDiscrete::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + // Get index and interpolation factor for inelastic grid + int i; + double f; + get_energy_index(energy_, E_in, i, f); + + int j; + sample_params(E_in, E_out, j, seed); // Sample outgoing cosine bin int m = mu_out_.shape()[2]; @@ -189,6 +271,30 @@ void IncoherentInelasticAEDiscrete::sample( mu = (1 - f) * mu_ijk + f * mu_i1jk; } +double IncoherentInelasticAEDiscrete::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + // Get index and interpolation factor for inelastic grid + int i; + double f; + get_energy_index(energy_, E_in, i, f); + int j; + sample_params(E_in, E_out, j, seed); + + int m = mu_out_.shape()[2]; + std::vector mu_vector; + mu_vector.reserve(m); + + for (int k = 0; k < m; ++k) { + double mu_ijk = mu_out_(i, j, k); + double mu_i1jk = mu_out_(i + 1, j, k); + double mu_k = (1 - f) * mu_ijk + f * mu_i1jk; + mu_vector.push_back(mu_k); + } + + return get_pdf_discrete(mu_vector, mu); +} + //============================================================================== // IncoherentInelasticAE implementation //============================================================================== @@ -231,24 +337,23 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group) } } -void IncoherentInelasticAE::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const +void IncoherentInelasticAE::sample_params( + double E_in, double& E_out, double& f, int& l, int& j, uint64_t* seed) const { // Get index and interpolation factor for inelastic grid int i; - double f; - get_energy_index(energy_, E_in, i, f); + double f0; + get_energy_index(energy_, E_in, i, f0); // Pick closer energy based on interpolation factor - int l = f > 0.5 ? i + 1 : i; + l = f0 > 0.5 ? i + 1 : i; // Determine outgoing energy bin // (First reset n_energy_out to the right value) - auto n = distribution_[l].n_e_out; + int n = distribution_[l].n_e_out; double r1 = prn(seed); double c_j = distribution_[l].e_out_cdf[0]; double c_j1; - std::size_t j; for (j = 0; j < n - 1; ++j) { c_j1 = distribution_[l].e_out_cdf[j + 1]; if (r1 < c_j1) @@ -286,6 +391,15 @@ void IncoherentInelasticAE::sample( E_out += E_in - E_l; } + f = (r1 - c_j) / (c_j1 - c_j); +} +void IncoherentInelasticAE::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + double f; + int l, j; + sample_params(E_in, E_out, f, l, j, seed); + // Sample outgoing cosine bin int n_mu = distribution_[l].mu.shape()[1]; std::size_t k = prn(seed) * n_mu; @@ -294,7 +408,6 @@ void IncoherentInelasticAE::sample( // a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the // discrete mu value itself. const auto& mu_l = distribution_[l].mu; - f = (r1 - c_j) / (c_j1 - c_j); // Interpolate kth mu value between distributions at energies j and j+1 mu = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k)); @@ -318,6 +431,25 @@ void IncoherentInelasticAE::sample( mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5); } +double IncoherentInelasticAE::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + double f; + int l, j; + sample_params(E_in, E_out, f, l, j, seed); + + int n_mu = distribution_[l].mu.shape()[1]; + const auto& mu_l = distribution_[l].mu; + std::vector mu_vector; + + for (int k = 0; k < n_mu; ++k) { + double mu_k = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k)); + mu_vector.push_back(mu_k); + } + + return get_pdf_discrete(mu_vector, mu); +} + //============================================================================== // MixedElasticAE implementation //============================================================================== @@ -340,18 +472,30 @@ MixedElasticAE::MixedElasticAE( close_group(incoherent_group); } -void MixedElasticAE::sample( - double E_in, double& E_out, double& mu, uint64_t* seed) const +const AngleEnergy& MixedElasticAE::sample_dist( + double E_in, uint64_t* seed) const { // Evaluate coherent and incoherent elastic cross sections double xs_coh = coherent_xs_(E_in); double xs_incoh = incoherent_xs_(E_in); if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) { - coherent_dist_.sample(E_in, E_out, mu, seed); + return coherent_dist_; } else { - incoherent_dist_->sample(E_in, E_out, mu, seed); + return *incoherent_dist_; } } +void MixedElasticAE::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + sample_dist(E_in, seed).sample(E_in, E_out, mu, seed); +} + +double MixedElasticAE::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed); +} + } // namespace openmc diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp index 5cbb76fb9b9..5746f824d17 100644 --- a/src/secondary_uncorrelated.cpp +++ b/src/secondary_uncorrelated.cpp @@ -65,4 +65,17 @@ void UncorrelatedAngleEnergy::sample( E_out = energy_->sample(E_in, seed); } +double UncorrelatedAngleEnergy::sample_energy_and_pdf( + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + // Sample outgoing energy + E_out = energy_->sample(E_in, seed); + if (!angle_.empty()) { + return angle_.evaluate(E_in, mu); + } else { + // no angle distribution given => assume isotropic for all energies + return 0.5; + } +} + } // namespace openmc diff --git a/src/thermal.cpp b/src/thermal.cpp index cbe0983ed65..1a06d1d19a3 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -296,16 +296,21 @@ void ThermalData::calculate_xs( *inelastic = (*inelastic_.xs)(E); } -void ThermalData::sample(const NuclideMicroXS& micro_xs, double E, - double* E_out, double* mu, uint64_t* seed) +AngleEnergy& ThermalData::sample_dist( + const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const { // Determine whether inelastic or elastic scattering will occur if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) { - elastic_.distribution->sample(E, *E_out, *mu, seed); + return *elastic_.distribution; } else { - inelastic_.distribution->sample(E, *E_out, *mu, seed); + return *inelastic_.distribution; } +} +void ThermalData::sample(const NuclideMicroXS& micro_xs, double E, + double* E_out, double* mu, uint64_t* seed) const +{ + sample_dist(micro_xs, E, seed).sample(E, *E_out, *mu, seed); // Because of floating-point roundoff, it may be possible for mu to be // outside of the range [-1,1). In these cases, we just set mu to exactly // -1 or 1 @@ -313,6 +318,13 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E, *mu = std::copysign(1.0, *mu); } +double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs, + double E_in, double mu, double& E_out, uint64_t* seed) const +{ + return sample_dist(micro_xs, E_in, seed) + .sample_energy_and_pdf(E_in, mu, E_out, seed); +} + void free_memory_thermal() { data::thermal_scatt.clear(); From 93d905c27b2579cbce334e1f68b7ad3725edf25a Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 15:03:40 +0200 Subject: [PATCH 03/18] wip --- CMakeLists.txt | 1 + include/openmc/math_functions.h | 6 ++++ openmc/filter.py | 49 +++++++++++++++++++++++++++++---- openmc/lib/filter.py | 4 ++- src/math_functions.cpp | 9 ++++++ src/tallies/filter.cpp | 3 ++ src/tallies/filter_point.cpp | 4 +-- src/tallies/tally_scoring.cpp | 8 ++++-- 8 files changed, 73 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 87b8789d101..bd9e93e7d7c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -425,6 +425,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_musurface.cpp src/tallies/filter_parent_nuclide.cpp src/tallies/filter_particle.cpp + src/tallies/filter_point.cpp src/tallies/filter_polar.cpp src/tallies/filter_sph_harm.cpp src/tallies/filter_sptl_legendre.cpp diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index 0d960c33db7..c5823635ebc 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -201,6 +201,12 @@ std::complex faddeeva(std::complex z); //! \return Derivative of Faddeeva function evaluated at z std::complex w_derivative(std::complex z, int order); +//! Evaluate relative exponential function +//! +//! \param x Real argument +//! \return (exp(x)-1)/x without loss of precision near 0 +double exprel(double x); + //! Helper function to get index and interpolation function on an incident //! energy grid //! diff --git a/openmc/filter.py b/openmc/filter.py index e7c54f2f0ee..b36c3c4ec91 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -22,11 +22,11 @@ _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', - 'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', - 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance', - 'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface', - 'meshmaterial', + 'energyout', 'mu', 'musurface', 'point', 'polar', 'azimuthal', 'distribcell', + 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', + 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', + 'cellinstance', 'collision', 'time', 'parentnuclide', 'weight', 'meshborn', + 'meshsurface', 'meshmaterial', ) _CURRENT_NAMES = ( @@ -787,6 +787,45 @@ def from_xml_element(cls, elem, **kwargs): return cls(bins, filter_id=filter_id) +class PointFilter(Filter): + """Bins tally events based on point detectors. + + Parameters + ---------- + bins : sequence of tuple[tuple[Real, Real, Real], Real] + Point detectors positions and exclusion radii. + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + bins : sequence of tuple[tuple[Real, Real, Real], Real] + Point detectors positions and exclusion radii. + id : int + Unique identifier for the filter + num_bins : Integral + The number of filter bins + + """ + def __eq__(self, other): + if type(self) is not type(other): + return False + elif len(self.bins) != len(other.bins): + return False + else: + return all(b1==b2 for b1,b2 in zip(self.bins,other.bins)) + + @Filter.bins.setter + def bins(self, bins): + cv.check_type('bins', bins, Sequence, tuple) + for i, item in enumerate(bins): + cv.check_type(f'bins[{i}]', item, tuple) + cv.check_length(f'bins[{i}]', item, 2, 2) + cv.check_type(f'bins[{i}][0]', item[0], tuple, Real) + cv.check_length(f'bins[{i}][0]', item[0], 3, 3) + cv.check_type(f'bins[{i}][1]', item[1], Real) + self._bins = bins + class ParentNuclideFilter(ParticleFilter): """Bins tally events based on the parent nuclide diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 5cf496107b3..e3cb96d65af 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -22,7 +22,7 @@ 'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter', 'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', - 'ParentNuclideFilter', 'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter', + 'ParentNuclideFilter', 'ParticleFilter', 'PointFilter', 'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' ] @@ -603,6 +603,8 @@ def bins(self): self._index, particle_i.ctypes.data_as(POINTER(c_int))) return [ParticleType(i) for i in particle_i] +class PointFilter(Filter): + filter_type = 'point' class PolarFilter(Filter): filter_type = 'polar' diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 9473f928b2c..20ca0fb52c5 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -919,6 +919,15 @@ std::complex w_derivative(std::complex z, int order) } } +double exprel(double x) +{ + if (std::abs(x) < 1e-16) + return 1.0; + else { + return std::expm1(x) / x; + } +} + // Helper function to get index and interpolation function on an incident energy // grid void get_energy_index( diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 79817981db0..29737b7b23a 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -31,6 +31,7 @@ #include "openmc/tallies/filter_musurface.h" #include "openmc/tallies/filter_parent_nuclide.h" #include "openmc/tallies/filter_particle.h" +#include "openmc/tallies/filter_point.h" #include "openmc/tallies/filter_polar.h" #include "openmc/tallies/filter_sph_harm.h" #include "openmc/tallies/filter_sptl_legendre.h" @@ -146,6 +147,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "particle") { return Filter::create(id); + } else if (type == "point") { + return Filter::create(id); } else if (type == "polar") { return Filter::create(id); } else if (type == "surface") { diff --git a/src/tallies/filter_point.cpp b/src/tallies/filter_point.cpp index 8d09da6f3d0..c044a33d12a 100644 --- a/src/tallies/filter_point.cpp +++ b/src/tallies/filter_point.cpp @@ -14,7 +14,7 @@ void PointFilter::from_xml(pugi::xml_node node) // Convert to vector of detectors vector> detectors; size_t n = bins.size() / 4; - for (int i = 0; i < n, ++i) { + for (int i = 0; i < n; ++i) { Position pos {bins[4 * i], bins[4 * i + 1], bins[4 * i + 2]}; detectors.push_back(std::make_pair(pos, bins[4 * i + 3])); } @@ -62,7 +62,7 @@ void PointFilter::to_statepoint(hid_t filter_group) const detectors.push_back(pos[2]); detectors.push_back(r); } - write_dataset(filter_group, "bins", particles); + write_dataset(filter_group, "bins", detectors); } std::string PointFilter::text_label(int bin) const diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 6d6a8632730..57f45cbf853 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2796,9 +2796,10 @@ void score_pseudoparticle_tally( match.bins_present_ = false; } -void score_point_tally(Particle& p, const Nuclide* nuc, const Reaction& rx, - int i_product, Direction* v_t) +void score_point_tally( + Particle& p, int i_nuclide, const Reaction& rx, int i_product, Direction* v_t) { + const auto& nuc {data::nuclides[i_nuclide]}; double awr = nuc->awr_; double E_in = p.E(); double vel = std::sqrt(E_in); @@ -2840,9 +2841,10 @@ void score_point_tally(Particle& p, const Nuclide* nuc, const Reaction& rx, } } -void score_point_tally(Particle& p, const Nuclide* nuc, const ThermalData& sab, +void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, const NuclideMicroXS& micro) { + const auto& nuc {data::nuclides[i_nuclide]}; double awr = nuc->awr_; double E_in = p.E(); double vel = std::sqrt(E_in); From 04c979721bece39e21091df07ff64260fbe7a468 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 15:20:15 +0200 Subject: [PATCH 04/18] wip --- openmc/filter.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openmc/filter.py b/openmc/filter.py index b36c3c4ec91..dd1793cd781 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -807,6 +807,9 @@ class PointFilter(Filter): The number of filter bins """ + + __hash__ = Filter.__hash__ + def __eq__(self, other): if type(self) is not type(other): return False From 8afed32bb7ffe61cd6c954582103c10771c7a279 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 15:39:40 +0200 Subject: [PATCH 05/18] wip --- openmc/filter.py | 18 ++++++++++++++++++ src/tallies/tally.cpp | 33 ++++++++++++++++++++++++++++++++- src/tallies/tally_scoring.cpp | 11 +++++++---- 3 files changed, 57 insertions(+), 5 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index dd1793cd781..572b22270f8 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -828,6 +828,24 @@ def bins(self, bins): cv.check_length(f'bins[{i}][0]', item[0], 3, 3) cv.check_type(f'bins[{i}][1]', item[1], Real) self._bins = bins + + def to_xml_element(self): + """Return XML Element representing the Filter. + + Returns + ------- + element : lxml.etree._Element + XML element containing filter data + + """ + element = ET.Element('filter') + element.set('id', str(self.id)) + element.set('type', self.short_name.lower()) + + subelement = ET.SubElement(element, 'bins') + subelement.text = ' '.join(str(b) for item in self.bins + for b in list(item[0])+[item[1]]) + return element class ParentNuclideFilter(ParticleFilter): """Bins tally events based on the parent nuclide diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 93b4f598c61..4c22d716cf1 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -31,6 +31,7 @@ #include "openmc/tallies/filter_meshmaterial.h" #include "openmc/tallies/filter_meshsurface.h" #include "openmc/tallies/filter_particle.h" +#include "openmc/tallies/filter_point.h" #include "openmc/tallies/filter_sph_harm.h" #include "openmc/tallies/filter_surface.h" #include "openmc/tallies/filter_time.h" @@ -538,6 +539,7 @@ void Tally::set_scores(const vector& scores) bool legendre_present = false; bool cell_present = false; bool cellfrom_present = false; + bool point_present = false; bool surface_present = false; bool meshsurface_present = false; bool non_cell_energy_present = false; @@ -554,6 +556,10 @@ void Tally::set_scores(const vector& scores) cellfrom_present = true; } else if (filt->type() == FilterType::CELL) { cell_present = true; + } else if (filt->type() == FilterType::POINT) { + point_present = true; + type_ = TallyType::NEXT_EVENT; + estimator_ = TallyEstimator::POINT; } else if (filt->type() == FilterType::SURFACE) { surface_present = true; } else if (filt->type() == FilterType::MESH_SURFACE) { @@ -561,6 +567,16 @@ void Tally::set_scores(const vector& scores) } } + if (point_present) { + if (legendre_present) + fatal_error("Cannot use LegendreFilter with PointFilter."); + if (energyout_present) + fatal_error("Cannot use EnergyoutFilter with PointFilter."); + if (surface_present || meshsurface_present) + fatal_error( + "Cannot use surface or mesh-surface filters with PointFilter."); + } + // Iterate over the given scores. for (auto score_str : scores) { // Make sure a delayed group filter wasn't used with an incompatible @@ -604,6 +620,9 @@ void Tally::set_scores(const vector& scores) case SCORE_NU_SCATTER: if (settings::run_CE) { + if (point_present) + fatal_error("Cannot use nu-scatter score with PointFilter in " + "continuous energy mode."); estimator_ = TallyEstimator::ANALOG; } else { if (energyout_present || legendre_present) @@ -612,6 +631,8 @@ void Tally::set_scores(const vector& scores) break; case SCORE_CURRENT: + if (point_present) + fatal_error("Cannot use current score with PointFilter."); // Check which type of current is desired: mesh or surface currents. if (surface_present || cell_present || cellfrom_present) { if (meshsurface_present) @@ -627,11 +648,15 @@ void Tally::set_scores(const vector& scores) break; case HEATING: + if (point_present) + fatal_error("Cannot use heating score with PointFilter."); if (settings::photon_transport) estimator_ = TallyEstimator::COLLISION; break; case SCORE_PULSE_HEIGHT: + if (point_present) + fatal_error("Cannot use pulse-height score with PointFilter."); if (non_cell_energy_present) { fatal_error("Pulse-height tallies are not compatible with filters " "other than CellFilter and EnergyFilter"); @@ -658,6 +683,8 @@ void Tally::set_scores(const vector& scores) case SCORE_IFP_TIME_NUM: case SCORE_IFP_BETA_NUM: case SCORE_IFP_DENOM: + if (point_present) + fatal_error("Cannot use ifp scores with PointFilter."); estimator_ = TallyEstimator::COLLISION; break; } @@ -1139,7 +1166,7 @@ void add_to_time_grid(vector grid) model::time_grid.swap(merged); } -//! Add new points to the global time grid +//! Add new points to global point_detectors // //! \param detector Position of new point detector to add void add_point_detector(Position& detector) @@ -1205,6 +1232,10 @@ void setup_active_tallies() switch (tally.estimator_) { case TallyEstimator::POINT: model::active_point_tallies.push_back(i); + auto pf = tally.get_filter(); + for (auto [det, r] : pf->detectors()) { + add_point_detector(det); + } break; } } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 57f45cbf853..e8c4e0050ed 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2748,11 +2748,13 @@ void score_pseudoparticle_tally( Particle& p, double distance, double mfp, double pdf) { double attenuation = std::exp(-mfp); - double flux = pdf * p.wgt() / (2 * PI * distance * distance); + double flux = pdf / (distance * distance); + + write_message(2, "wgt {} wgt_last {} flux {}", p.wgt(), p.wgt_last(), flux); // Save the attenuation for point filter handling - p.wgt_last() = p.wgt(); - p.wgt() *= attenuation; + //p.wgt_last() = p.wgt(); + //p.wgt() *= attenuation; for (auto i_tally : model::active_point_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2808,6 +2810,7 @@ void score_point_tally( Direction u_cm = (v_n - v_cm); u_cm /= u_cm.norm(); double mfp; + fatal_error("moo"); for (auto& det : model::active_point_detectors) { @@ -2844,6 +2847,7 @@ void score_point_tally( void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, const NuclideMicroXS& micro) { + fatal_error("moo"); const auto& nuc {data::nuclides[i_nuclide]}; double awr = nuc->awr_; double E_in = p.E(); @@ -2886,7 +2890,6 @@ void score_point_tally(SourceSite& site, int source_index) double mfp; for (auto& det : model::active_point_detectors) { - auto u = (det - r); double distance = u.norm(); u /= distance; From a27738759123cd73686f4c5119cd6dff091b2758 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 17:17:51 +0200 Subject: [PATCH 06/18] wip --- src/tallies/filter_point.cpp | 8 +++--- src/tallies/tally_scoring.cpp | 48 +++-------------------------------- 2 files changed, 9 insertions(+), 47 deletions(-) diff --git a/src/tallies/filter_point.cpp b/src/tallies/filter_point.cpp index c044a33d12a..2dd637ee728 100644 --- a/src/tallies/filter_point.cpp +++ b/src/tallies/filter_point.cpp @@ -37,15 +37,17 @@ void PointFilter::set_detectors(span> detectors) void PointFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { + double attenuation = p.wgt() / p.wgt_last(); for (auto i = 0; i < detectors_.size(); i++) { auto [pos, r] = detectors_[i]; if ((p.r() - pos).norm() < FP_COINCIDENT) { match.bins_.push_back(i); double weight; - if ((p.r_last() - pos).norm() < r) { - weight = 3.0 * exprel(-r * p.macro_xs().total); + double distance = (p.r_last() - pos).norm(); + if (distance > r) { + weight = attenuation / (distance * distance); } else { - weight = p.wgt() / p.wgt_last(); + weight = 3.0 * -exprel(-p.macro_xs().total * r) / (r * r); } match.weights_.push_back(weight); } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e8c4e0050ed..1278b68c04f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2748,54 +2748,14 @@ void score_pseudoparticle_tally( Particle& p, double distance, double mfp, double pdf) { double attenuation = std::exp(-mfp); - double flux = pdf / (distance * distance); - - write_message(2, "wgt {} wgt_last {} flux {}", p.wgt(), p.wgt_last(), flux); // Save the attenuation for point filter handling - //p.wgt_last() = p.wgt(); - //p.wgt() *= attenuation; + p.wgt_last() = p.wgt(); + p.wgt() *= attenuation; - for (auto i_tally : model::active_point_tallies) { - const Tally& tally {*model::tallies[i_tally]}; - - // 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 nuclide bins. - for (auto i = 0; i < tally.nuclides_.size(); ++i) { - auto i_nuclide = tally.nuclides_[i]; - - // Tally this event in the present nuclide bin if that bin represents - // the event nuclide or the total material. Note that the atomic - // density argument for score_general is not used for analog tallies. - if (i_nuclide == p.event_nuclide() || i_nuclide == -1) - score_general_ce_analog(p, i_tally, i * tally.scores_.size(), - filter_index, filter_weight, i_nuclide, -1.0, flux); - } - } + double flux = p.wgt_last() * pdf; - // If the user has specified that we can assume all tallies are spatially - // separate, this implies that once a tally has been scored to, we needn't - // check the others. This cuts down on overhead when there are many - // tallies specified - if (settings::assume_separate) - break; - } - - // Reset all the filter matches for the next tally event. - for (auto& match : p.filter_matches()) - match.bins_present_ = false; + score_tracklength_tally_general(p, flux, model::active_point_tallies); } void score_point_tally( From dd49099492cd9733f096b746af47a355fc86055c Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 17:33:59 +0200 Subject: [PATCH 07/18] wip --- include/openmc/simulation.h | 14 ++++++++------ src/simulation.cpp | 1 + src/surface.cpp | 3 +++ src/tallies/tally.cpp | 3 +++ 4 files changed, 15 insertions(+), 6 deletions(-) diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index c2e04fb628e..64902eba235 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -36,12 +36,14 @@ extern "C" double extern double log_spacing; //!< lethargy spacing for energy grid searches extern "C" int n_lost_particles; //!< cumulative number of lost particles extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? -extern "C" int restart_batch; //!< batch at which a restart job resumed -extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? -extern int ssw_current_file; //!< current surface source file -extern "C" int total_gen; //!< total number of generations simulated -extern double total_weight; //!< Total source weight in a batch -extern int64_t work_per_rank; //!< number of particles per MPI rank +extern bool + nonvacuum_boundary_present; //!< Does the geometry contain non-vacuum b.c. +extern "C" int restart_batch; //!< batch at which a restart job resumed +extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? +extern int ssw_current_file; //!< current surface source file +extern "C" int total_gen; //!< total number of generations simulated +extern double total_weight; //!< Total source weight in a batch +extern int64_t work_per_rank; //!< number of particles per MPI rank extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; diff --git a/src/simulation.cpp b/src/simulation.cpp index df4c789ad20..16f78cb95cd 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -311,6 +311,7 @@ double k_abs_tra {0.0}; double log_spacing; int n_lost_particles {0}; bool need_depletion_rx {false}; +bool nonvacuum_boundary_present {false}; int restart_batch; bool satisfy_triggers {false}; int ssw_current_file; diff --git a/src/surface.cpp b/src/surface.cpp index 81b756deae7..d52a20d5f81 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1175,6 +1175,7 @@ void read_surfaces(pugi::xml_node node, std::unordered_map& albedo_map, std::unordered_map& periodic_sense_map) { + simulation::nonvacuum_boundary_present = false; // Count the number of surfaces int n_surfaces = 0; for (pugi::xml_node surf_node : node.children("surface")) { @@ -1245,6 +1246,8 @@ void read_surfaces(pugi::xml_node node, // Check for a periodic surface if (check_for_node(surf_node, "boundary")) { std::string surf_bc = get_node_value(surf_node, "boundary", true, true); + if (surf_bc != "vacuum") + simulation::nonvacuum_boundary_present = true; if (surf_bc == "periodic") { periodic_sense_map[model::surfaces.back()->id_] = 0; // Check for surface albedo. Skip sanity check as it is already done diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 4c22d716cf1..23694ec85a8 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -568,6 +568,9 @@ void Tally::set_scores(const vector& scores) } if (point_present) { + if (simulation::nonvacuum_boundary_present) + fatal_error( + "Cannot use point detectors with non-vacuum boundary conditions."); if (legendre_present) fatal_error("Cannot use LegendreFilter with PointFilter."); if (energyout_present) From 9710d593155666d5732bda11b657df6b2d29a3be Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 17:36:54 +0200 Subject: [PATCH 08/18] wip --- src/tallies/tally_scoring.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 1278b68c04f..75af18075ca 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2770,7 +2770,6 @@ void score_point_tally( Direction u_cm = (v_n - v_cm); u_cm /= u_cm.norm(); double mfp; - fatal_error("moo"); for (auto& det : model::active_point_detectors) { @@ -2807,7 +2806,6 @@ void score_point_tally( void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, const NuclideMicroXS& micro) { - fatal_error("moo"); const auto& nuc {data::nuclides[i_nuclide]}; double awr = nuc->awr_; double E_in = p.E(); From e800c4712a3f1c55059a7d0fdc3cfc1c1ed98c61 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 18:19:30 +0200 Subject: [PATCH 09/18] wip --- include/openmc/tallies/filter_point.h | 2 ++ include/openmc/tallies/tally_scoring.h | 5 +++++ src/surface.cpp | 1 + src/tallies/filter_point.cpp | 18 +++++++++++++++++- src/tallies/tally.cpp | 7 +++++++ src/tallies/tally_scoring.cpp | 10 +++++++++- 6 files changed, 41 insertions(+), 2 deletions(-) diff --git a/include/openmc/tallies/filter_point.h b/include/openmc/tallies/filter_point.h index 909f0baa27d..e2f3908d746 100644 --- a/include/openmc/tallies/filter_point.h +++ b/include/openmc/tallies/filter_point.h @@ -43,12 +43,14 @@ class PointFilter : public Filter { } void set_detectors(span> detectors); + void reset_indices(); private: //---------------------------------------------------------------------------- // Data members vector> detectors_; + vector> indices_; }; } // namespace openmc diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 2ebce2ec35d..0e132a8a48d 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -9,6 +9,11 @@ namespace openmc { +namespace simulation { +extern thread_local Particle pseudoparticle; +extern thread_local int i_det; +} // namespace simulation + //============================================================================== //! An iterator over all combinations of a tally's matching filter bins. // diff --git a/src/surface.cpp b/src/surface.cpp index d52a20d5f81..08cdc494807 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -17,6 +17,7 @@ #include "openmc/math_functions.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/string_utils.h" #include "openmc/xml_interface.h" diff --git a/src/tallies/filter_point.cpp b/src/tallies/filter_point.cpp index 2dd637ee728..cfc5ee0ff40 100644 --- a/src/tallies/filter_point.cpp +++ b/src/tallies/filter_point.cpp @@ -1,8 +1,10 @@ #include "openmc/tallies/filter_point.h" +#include "openmc/tallies/tally_scoring.h" #include #include "openmc/math_functions.h" +#include "openmc/simulation.h" #include "openmc/xml_interface.h" namespace openmc { @@ -34,11 +36,25 @@ void PointFilter::set_detectors(span> detectors) n_bins_ = detectors_.size(); } +void PointFilter::reset_indices() +{ + indices_.clear(); + for (auto det : model::active_point_detectors) { + vector temp; + for (auto i = 0; i < detectors_.size(); i++) { + auto [d, r] = detectors_[i]; + if (d == det) + temp.push_back(i); + } + indices_.push_back(temp); + } +} + void PointFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { double attenuation = p.wgt() / p.wgt_last(); - for (auto i = 0; i < detectors_.size(); i++) { + for (auto i : indices_[simulation::i_det]) { auto [pos, r] = detectors_[i]; if ((p.r() - pos).norm() < FP_COINCIDENT) { match.bins_.push_back(i); diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 23694ec85a8..1bb5f74c4a5 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1191,6 +1191,8 @@ void setup_active_tallies() model::time_grid.clear(); model::active_point_detectors.clear(); + std::set particle_filter_ids; + for (auto i = 0; i < model::tallies.size(); ++i) { const auto& tally {*model::tallies[i]}; @@ -1239,11 +1241,16 @@ void setup_active_tallies() for (auto [det, r] : pf->detectors()) { add_point_detector(det); } + particle_filter_ids.insert(pf->id()); break; } } } } + for (auto id : particle_filter_ids) { + auto pf = dynamic_cast(model::tally_filters.at(id).get()); + pf->reset_indices(); + } } void free_memory_tally() diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 75af18075ca..26ccc26bdc0 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -27,7 +27,8 @@ namespace openmc { namespace simulation { thread_local Particle pseudoparticle; -} +thread_local int i_det; +} // namespace simulation //============================================================================== // FilterBinIter implementation @@ -2771,7 +2772,9 @@ void score_point_tally( u_cm /= u_cm.norm(); double mfp; + simulation::i_det = -1; for (auto& det : model::active_point_detectors) { + ++simulation::i_det; auto u = (det - p.r()); double distance = u.norm(); @@ -2816,7 +2819,9 @@ void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, u_cm /= u_cm.norm(); double mfp; + simulation::i_det = -1; for (auto& det : model::active_point_detectors) { + ++simulation::i_det; auto u = (det - p.r()); double distance = u.norm(); @@ -2847,7 +2852,10 @@ void score_point_tally(SourceSite& site, int source_index) double mfp; + simulation::i_det = -1; for (auto& det : model::active_point_detectors) { + ++simulation::i_det; + auto u = (det - r); double distance = u.norm(); u /= distance; From a5f480996cae950d51c0d6a6f780fc6cbc651412 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 18:25:11 +0200 Subject: [PATCH 10/18] wip --- src/tallies/tally.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 1bb5f74c4a5..c2f81dd469b 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1241,7 +1241,7 @@ void setup_active_tallies() for (auto [det, r] : pf->detectors()) { add_point_detector(det); } - particle_filter_ids.insert(pf->id()); + particle_filter_ids.insert(model::tally_map[pf->id()]); break; } } From 753b8496f509fd315a46a17db722e9d7857ce968 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 18:42:13 +0200 Subject: [PATCH 11/18] wip --- include/openmc/tallies/tally_scoring.h | 3 +-- src/tallies/tally_scoring.cpp | 9 ++++----- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 0e132a8a48d..b3b653383d8 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -129,8 +129,7 @@ void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, void score_point_tally(SourceSite& site, int source_index); -void score_pseudoparticle_tally( - Particle& p, double distance, double mfp, double pdf); +void score_pseudoparticle_tally(Particle& p, double mfp, double pdf); } // namespace openmc #endif // OPENMC_TALLIES_TALLY_SCORING_H diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 26ccc26bdc0..44d816630a1 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2745,8 +2745,7 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) } } -void score_pseudoparticle_tally( - Particle& p, double distance, double mfp, double pdf) +void score_pseudoparticle_tally(Particle& p, double mfp, double pdf) { double attenuation = std::exp(-mfp); @@ -2802,7 +2801,7 @@ void score_point_tally( transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); if (!simulation::pseudoparticle.alive()) continue; - score_pseudoparticle_tally(simulation::pseudoparticle, distance, mfp, pdf); + score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); } } @@ -2836,7 +2835,7 @@ void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); if (!simulation::pseudoparticle.alive()) continue; - score_pseudoparticle_tally(simulation::pseudoparticle, distance, mfp, pdf); + score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); } } @@ -2867,7 +2866,7 @@ void score_point_tally(SourceSite& site, int source_index) transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); if (!simulation::pseudoparticle.alive()) continue; - score_pseudoparticle_tally(simulation::pseudoparticle, distance, mfp, pdf); + score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); } } From 9d16542dd8c88bf6cc51595e733c8cd0b7f82c9f Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 19:37:39 +0200 Subject: [PATCH 12/18] wip --- src/secondary_uncorrelated.cpp | 7 ++++++- src/tallies/tally_scoring.cpp | 4 +++- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp index 5746f824d17..ec2af710258 100644 --- a/src/secondary_uncorrelated.cpp +++ b/src/secondary_uncorrelated.cpp @@ -69,7 +69,12 @@ double UncorrelatedAngleEnergy::sample_energy_and_pdf( double E_in, double mu, double& E_out, uint64_t* seed) const { // Sample outgoing energy - E_out = energy_->sample(E_in, seed); + if (energy_ != nullptr) { + E_out = energy_->sample(E_in, seed); + } else { + E_out = E_in; + } + if (!angle_.empty()) { return angle_.evaluate(E_in, mu); } else { diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 44d816630a1..0c19803a509 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2766,7 +2766,9 @@ void score_point_tally( double E_in = p.E(); double vel = std::sqrt(E_in); Direction v_n = vel * p.u(); - Direction v_cm = (v_n + awr * (v_t != nullptr ? *v_t : 0)) / (awr + 1.0); + Direction v_cm = v_n / (awr + 1.0); + if (v_t != nullptr) + v_cm += awr / (awr + 1.0) * (*v_t); Direction u_cm = (v_n - v_cm); u_cm /= u_cm.norm(); double mfp; From 19827823dce1e1aa40447128701d7398aeacee87 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 20:20:03 +0200 Subject: [PATCH 13/18] try to use another stream for next event tracking --- include/openmc/random_lcg.h | 6 +++++- src/tallies/tally_scoring.cpp | 8 ++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h index 5aecafed3cb..88a3fe23739 100644 --- a/include/openmc/random_lcg.h +++ b/include/openmc/random_lcg.h @@ -9,11 +9,15 @@ namespace openmc { // Module constants. //============================================================================== -constexpr int N_STREAMS {4}; +constexpr int N_STREAMS {5}; constexpr int STREAM_TRACKING {0}; constexpr int STREAM_SOURCE {1}; constexpr int STREAM_URR_PTABLE {2}; constexpr int STREAM_VOLUME {3}; +constexpr int STREAM_NEXT_EVENT +{ + 4 +} constexpr int64_t DEFAULT_SEED {1}; constexpr uint64_t DEFAULT_STRIDE {152917ULL}; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 0c19803a509..4ff3dc0d312 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2773,6 +2773,9 @@ void score_point_tally( u_cm /= u_cm.norm(); double mfp; + auto old_stream = p.stream(); + p.stream() = STREAM_NEXT_EVENT; + simulation::i_det = -1; for (auto& det : model::active_point_detectors) { ++simulation::i_det; @@ -2805,6 +2808,7 @@ void score_point_tally( continue; score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); } + p.stream() = old_stream; } void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, @@ -2820,6 +2824,9 @@ void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, u_cm /= u_cm.norm(); double mfp; + auto old_stream = p.stream(); + p.stream() = STREAM_NEXT_EVENT; + simulation::i_det = -1; for (auto& det : model::active_point_detectors) { ++simulation::i_det; @@ -2839,6 +2846,7 @@ void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, continue; score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); } + p.stream() = old_stream; } void score_point_tally(SourceSite& site, int source_index) From d08626053556d2fd204c09884d66029e5efff7a6 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 20:21:02 +0200 Subject: [PATCH 14/18] try to use another stream for next event tracking --- include/openmc/random_lcg.h | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h index 88a3fe23739..cc6ccfac47c 100644 --- a/include/openmc/random_lcg.h +++ b/include/openmc/random_lcg.h @@ -14,10 +14,7 @@ constexpr int STREAM_TRACKING {0}; constexpr int STREAM_SOURCE {1}; constexpr int STREAM_URR_PTABLE {2}; constexpr int STREAM_VOLUME {3}; -constexpr int STREAM_NEXT_EVENT -{ - 4 -} +constexpr int STREAM_NEXT_EVENT {4}; constexpr int64_t DEFAULT_SEED {1}; constexpr uint64_t DEFAULT_STRIDE {152917ULL}; From 1524e64c4f54fa1aae7dc9fa1f1102d0e008ddd8 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 30 Jan 2026 00:36:23 +0200 Subject: [PATCH 15/18] wip edit --- include/openmc/particle.h | 4 --- src/particle.cpp | 55 ----------------------------------- src/reaction_product.cpp | 4 +-- src/simulation.cpp | 36 ++++++++++++----------- src/tallies/tally_scoring.cpp | 6 ++-- 5 files changed, 23 insertions(+), 82 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 14face2d059..362c5956dc7 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -76,10 +76,6 @@ class Particle : public ParticleData { void event_revive_from_secondary(); void event_death(); - // Pseudo-particles events - void event_advance_pseudo(double& total_distance, double& mfp); - void event_cross_surface_pseudo(); - //! pulse-height recording void pht_collision_energy(); void pht_secondary_particles(); diff --git a/src/particle.cpp b/src/particle.cpp index 1843f28bb5f..922e6245700 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -549,61 +549,6 @@ void Particle::event_death() } } -void Particle::event_advance_pseudo(double& total_distance, double& mfp) -{ - // Find the distance to the nearest boundary - boundary() = distance_to_boundary(*this); - - double speed = this->speed(); - double time_cutoff = settings::time_cutoff[static_cast(type())]; - double distance_cutoff = - (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; - - if (distance_cutoff < boundary().distance()) { - wgt() = 0.0; - return; - } - - // Select smaller of the two distances - double distance = std::min({boundary().distance(), total_distance}); - - // Advance particle in space and time - this->move_distance(distance); - double dt = distance / speed; - this->time() += dt; - this->lifetime() += dt; - mfp += distance * macro_xs().total; - total_distance -= distance; -} - -void Particle::event_cross_surface_pseudo() -{ - // Saving previous cell data - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); - - // Set surface that particle is on and adjust coordinate levels - surface() = boundary().surface(); - n_coord() = boundary().coord_level(); - - if (boundary().lattice_translation()[0] != 0 || - boundary().lattice_translation()[1] != 0 || - boundary().lattice_translation()[2] != 0) { - // Particle crosses lattice boundary - - bool verbose = settings::verbosity >= 10 || trace(); - cross_lattice(*this, boundary(), verbose); - event() = TallyEvent::LATTICE; - } else { - // Particle crosses surface - const auto& surf {model::surfaces[surface_index()].get()}; - this->cross_surface(*surf); - event() = TallyEvent::SURFACE; - } -} - void Particle::pht_collision_energy() { // Adds the energy particles lose in a collision to the pulse-height diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 398f72e7967..6cfd164a01d 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -117,10 +117,8 @@ AngleEnergy& ReactionProduct::sample_dist(double E_in, uint64_t* seed) const prob += applicability_[i](E_in); // If i-th distribution is sampled, sample energy from the distribution - if (c <= prob) { + if (c <= prob) return *distribution_[i]; - break; - } } } else { // If only one distribution is present, go ahead and sample it diff --git a/src/simulation.cpp b/src/simulation.cpp index 16f78cb95cd..3262ad446cb 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -7,10 +7,12 @@ #include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/event.h" +#include "openmc/geometry.h" #include "openmc/geometry_aux.h" #include "openmc/ifp.h" #include "openmc/material.h" #include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" #include "openmc/output.h" #include "openmc/particle.h" @@ -824,25 +826,25 @@ void transport_history_based_single_particle(Particle& p) void transport_pseudoparticle(Particle& p, double total_distance, double& mfp) { - if (!p.alive()) - return; - - double distance = total_distance; - bool cross_surface = false; - while (distance > 0.0) { + double remaining_distance = total_distance; + p.event_calculate_xs(); + p.boundary() = distance_to_boundary(p); + double advance_distance = p.boundary().distance(); + + while (advance_distance < remaining_distance) { + mfp += advance_distance * p.macro_xs().total; + // Advance particle in space and time + p.move_distance(advance_distance); + remaining_distance -= advance_distance; + p.time() += advance_distance / p.speed(); + p.event_cross_surface(); p.event_calculate_xs(); - if (p.alive()) { - cross_surface = (distance > p.boundary().distance()); - p.event_advance_pseudo(distance, mfp); - } else { - return; - } - if (p.alive() && cross_surface) { - p.event_cross_surface_pseudo(); - } - if (!p.alive()) - return; + p.boundary() = distance_to_boundary(p); + advance_distance = p.boundary().distance(); } + mfp += remaining_distance * p.macro_xs().total; + p.move_distance(remaining_distance); + p.time() += remaining_distance / p.speed(); } void transport_history_based() diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 4ff3dc0d312..4040d2320bb 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2771,7 +2771,7 @@ void score_point_tally( v_cm += awr / (awr + 1.0) * (*v_t); Direction u_cm = (v_n - v_cm); u_cm /= u_cm.norm(); - double mfp; + double mfp = 0; auto old_stream = p.stream(); p.stream() = STREAM_NEXT_EVENT; @@ -2822,7 +2822,7 @@ void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, Direction v_cm = v_n / (awr + 1.0); Direction u_cm = (v_n - v_cm); u_cm /= u_cm.norm(); - double mfp; + double mfp = 0; auto old_stream = p.stream(); p.stream() = STREAM_NEXT_EVENT; @@ -2859,7 +2859,7 @@ void score_point_tally(SourceSite& site, int source_index) fatal_error("Only independent source is valid for point detectors."); auto angle = src->angle(); - double mfp; + double mfp = 0; simulation::i_det = -1; for (auto& det : model::active_point_detectors) { From 84ebc85c056d6cb14f8fb9579d4c57c86453fce1 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 30 Jan 2026 00:55:54 +0200 Subject: [PATCH 16/18] still works --- include/openmc/particle.h | 4 ++-- src/particle.cpp | 16 ++++++++-------- src/simulation.cpp | 6 +++--- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 362c5956dc7..6b92b12f190 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -69,9 +69,9 @@ class Particle : public ParticleData { void initialize_pseudoparticle(Particle& p, Direction u, double E); // Coarse-grained particle events - void event_calculate_xs(); + void event_calculate_xs(bool is_pseudo = false); void event_advance(); - void event_cross_surface(); + void event_cross_surface(bool is_pseudo = false); void event_collide(); void event_revive_from_secondary(); void event_death(); diff --git a/src/particle.cpp b/src/particle.cpp index 922e6245700..337914de138 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -195,10 +195,10 @@ void Particle::initialize_pseudoparticle( time_last() = p.time(); } -void Particle::event_calculate_xs() +void Particle::event_calculate_xs(bool is_pseudo) { // Set the random number stream - stream() = STREAM_TRACKING; + stream() = is_pseudo ? STREAM_NEXT_EVENT : STREAM_TRACKING; // Store pre-collision particle properties wgt_last() = wgt(); @@ -234,7 +234,7 @@ void Particle::event_calculate_xs() } // Write particle track. - if (write_track()) + if (write_track() && !is_pseudo) write_particle_track(*this); if (settings::check_overlaps) @@ -323,7 +323,7 @@ void Particle::event_advance() } } -void Particle::event_cross_surface() +void Particle::event_cross_surface(bool is_pseudo) { // Saving previous cell data for (int j = 0; j < n_coord(); ++j) { @@ -347,21 +347,21 @@ void Particle::event_cross_surface() // Particle crosses surface const auto& surf {model::surfaces[surface_index()].get()}; // If BC, add particle to surface source before crossing surface - if (surf->surf_source_ && surf->bc_) { + if (surf->surf_source_ && surf->bc_ && !is_pseudo) { add_surf_source_to_bank(*this, *surf); } this->cross_surface(*surf); // If no BC, add particle to surface source after crossing surface - if (surf->surf_source_ && !surf->bc_) { + if (surf->surf_source_ && !surf->bc_ && !is_pseudo) { add_surf_source_to_bank(*this, *surf); } - if (settings::weight_window_checkpoint_surface) { + if (settings::weight_window_checkpoint_surface && !is_pseudo) { apply_weight_windows(*this); } event() = TallyEvent::SURFACE; } // Score cell to cell partial currents - if (!model::active_surface_tallies.empty()) { + if (!model::active_surface_tallies.empty() && !is_pseudo) { score_surface_tally(*this, model::active_surface_tallies); } } diff --git a/src/simulation.cpp b/src/simulation.cpp index 3262ad446cb..7f27cde8a98 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -827,7 +827,7 @@ void transport_history_based_single_particle(Particle& p) void transport_pseudoparticle(Particle& p, double total_distance, double& mfp) { double remaining_distance = total_distance; - p.event_calculate_xs(); + p.event_calculate_xs(true); p.boundary() = distance_to_boundary(p); double advance_distance = p.boundary().distance(); @@ -837,8 +837,8 @@ void transport_pseudoparticle(Particle& p, double total_distance, double& mfp) p.move_distance(advance_distance); remaining_distance -= advance_distance; p.time() += advance_distance / p.speed(); - p.event_cross_surface(); - p.event_calculate_xs(); + p.event_cross_surface(true); + p.event_calculate_xs(true); p.boundary() = distance_to_boundary(p); advance_distance = p.boundary().distance(); } From 3e2a26bee4bba09344d41ee04b93c60bc2666622 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 30 Jan 2026 01:06:57 +0200 Subject: [PATCH 17/18] still works --- include/openmc/simulation.h | 2 +- src/simulation.cpp | 38 ++++++++++++++++++++--------------- src/tallies/tally_scoring.cpp | 9 +++------ 3 files changed, 26 insertions(+), 23 deletions(-) diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 64902eba235..42ab43c1fca 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -95,7 +95,7 @@ void broadcast_results(); void free_memory_simulation(); //! Simulate a single pseudoparticle history -void transport_pseudoparticle(Particle& p, double total_distance, double& mfp); +double transport_pseudoparticle(Particle& p, double total_distance); //! Simulate a single particle history (and all generated secondary particles, //! if enabled), from birth to death diff --git a/src/simulation.cpp b/src/simulation.cpp index 7f27cde8a98..d3f94165fca 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -824,27 +824,33 @@ void transport_history_based_single_particle(Particle& p) p.event_death(); } -void transport_pseudoparticle(Particle& p, double total_distance, double& mfp) +double transport_pseudoparticle(Particle& p, double total_distance) { - double remaining_distance = total_distance; - p.event_calculate_xs(true); - p.boundary() = distance_to_boundary(p); - double advance_distance = p.boundary().distance(); + double mfp = 0.0; + while (total_distance > 0.0) { + p.event_calculate_xs(true); + p.boundary() = distance_to_boundary(p); + double speed = p.speed(); + double time_cutoff = settings::time_cutoff[static_cast(p.type())]; + double distance_cutoff = + (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; + + double distance = + std::min({p.boundary().distance(), distance_cutoff, total_distance}); + if (distance == distance_cutoff) { + p.wgt() = 0.0; + return; + } - while (advance_distance < remaining_distance) { - mfp += advance_distance * p.macro_xs().total; // Advance particle in space and time p.move_distance(advance_distance); - remaining_distance -= advance_distance; - p.time() += advance_distance / p.speed(); - p.event_cross_surface(true); - p.event_calculate_xs(true); - p.boundary() = distance_to_boundary(p); - advance_distance = p.boundary().distance(); + p.time() += distance / speed; + p.lifetime() += distance / speed; + mfp += distance * p.macro_xs().total; + if (distance == p.boundary().distance()) + p.event_cross_surface(true); } - mfp += remaining_distance * p.macro_xs().total; - p.move_distance(remaining_distance); - p.time() += remaining_distance / p.speed(); + return mfp; } void transport_history_based() diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 4040d2320bb..4db3bec70f8 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2802,8 +2802,7 @@ void score_point_tally( (1 - mu / (awr + 1) * std::sqrt(E_in / E_out)); } simulation::pseudoparticle.initialize_pseudoparticle(p, u, E_out); - mfp = 0.0; - transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); + double mfp = transport_pseudoparticle(simulation::pseudoparticle, distance); if (!simulation::pseudoparticle.alive()) continue; score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); @@ -2840,8 +2839,7 @@ void score_point_tally(Particle& p, int i_nuclide, const ThermalData& sab, double pdf = sab.sample_energy_and_pdf(micro, p.E(), mu, E_out, p.current_seed()); simulation::pseudoparticle.initialize_pseudoparticle(p, u, E_out); - mfp = 0.0; - transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); + double mfp = transport_pseudoparticle(simulation::pseudoparticle, distance); if (!simulation::pseudoparticle.alive()) continue; score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); @@ -2872,8 +2870,7 @@ void score_point_tally(SourceSite& site, int source_index) double pdf = angle->evaluate(u); simulation::pseudoparticle.from_source(&site); simulation::pseudoparticle.u() = u; - mfp = 0.0; - transport_pseudoparticle(simulation::pseudoparticle, distance, mfp); + double mfp = transport_pseudoparticle(simulation::pseudoparticle, distance); if (!simulation::pseudoparticle.alive()) continue; score_pseudoparticle_tally(simulation::pseudoparticle, mfp, pdf); From 00995505f730b54e10d2f9bd8e2d1333f7d83786 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 30 Jan 2026 01:22:39 +0200 Subject: [PATCH 18/18] fix typo --- src/simulation.cpp | 22 +++++++++++++++------- src/tallies/filter_point.cpp | 2 +- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/simulation.cpp b/src/simulation.cpp index d3f94165fca..453bb6f20dc 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -826,29 +826,37 @@ void transport_history_based_single_particle(Particle& p) double transport_pseudoparticle(Particle& p, double total_distance) { + double time_cutoff = settings::time_cutoff[static_cast(p.type())]; + double speed = p.speed(); + p.event_calculate_xs(true); + double mfp = 0.0; while (total_distance > 0.0) { - p.event_calculate_xs(true); p.boundary() = distance_to_boundary(p); - double speed = p.speed(); - double time_cutoff = settings::time_cutoff[static_cast(p.type())]; + double distance_cutoff = - (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; + (time_cutoff < INFTY) ? (time_cutoff - p.time()) * speed : INFTY; double distance = std::min({p.boundary().distance(), distance_cutoff, total_distance}); if (distance == distance_cutoff) { p.wgt() = 0.0; - return; + return mfp; } // Advance particle in space and time - p.move_distance(advance_distance); + p.move_distance(distance); p.time() += distance / speed; p.lifetime() += distance / speed; + total_distance -= distance; mfp += distance * p.macro_xs().total; - if (distance == p.boundary().distance()) + + if (distance == p.boundary().distance()) { p.event_cross_surface(true); + p.event_calculate_xs(true); + if (!settings::run_CE) + speed = p.speed(); + } } return mfp; } diff --git a/src/tallies/filter_point.cpp b/src/tallies/filter_point.cpp index cfc5ee0ff40..0a9587f6650 100644 --- a/src/tallies/filter_point.cpp +++ b/src/tallies/filter_point.cpp @@ -63,7 +63,7 @@ void PointFilter::get_all_bins( if (distance > r) { weight = attenuation / (distance * distance); } else { - weight = 3.0 * -exprel(-p.macro_xs().total * r) / (r * r); + weight = 3.0 * exprel(-p.macro_xs().total * r) / (r * r); } match.weights_.push_back(weight); }