From 42d06c045384a5c8b6aa9722234f19d614e159e6 Mon Sep 17 00:00:00 2001 From: Itay Date: Sat, 2 Aug 2025 09:17:51 +0000 Subject: [PATCH] Implement relativistic point detector tally --- .gitmodules | 3 + include/openmc/constants.h | 2 +- include/openmc/distribution.h | 17 +- include/openmc/distribution_angle.h | 3 +- include/openmc/hdf5_interface.h | 2 +- include/openmc/particle.h | 11 +- include/openmc/particle_data.h | 7 + include/openmc/physics.h | 5 + include/openmc/reaction_product.h | 8 +- include/openmc/secondary_correlated.h | 5 + include/openmc/secondary_kalbach.h | 5 + include/openmc/secondary_nbody.h | 5 + include/openmc/secondary_thermal.h | 6 + include/openmc/secondary_uncorrelated.h | 5 + include/openmc/tallies/tally.h | 8 + include/openmc/tallies/tally_scoring.h | 114 +++++- include/openmc/thermal.h | 2 + include/openmc/urr.h | 8 +- openmc/tallies.py | 104 +++-- prints | 36 ++ src/distribution.cpp | 77 ++++ src/distribution_angle.cpp | 33 ++ src/output.cpp | 47 +-- src/particle.cpp | 95 +++++ src/particle_data.cpp | 7 + src/physics.cpp | 137 ++++++- src/reaction_product.cpp | 68 ++++ src/secondary_correlated.cpp | 151 ++++++++ src/secondary_kalbach.cpp | 141 +++++++ src/secondary_nbody.cpp | 63 ++++ src/secondary_thermal.cpp | 480 ++++++++++++++++++++++++ src/secondary_uncorrelated.cpp | 59 +++ src/source.cpp | 3 +- src/state_point.cpp | 2 + src/tallies/tally.cpp | 48 +++ src/tallies/tally_scoring.cpp | 477 ++++++++++++++++++++++- src/thermal.cpp | 54 ++- vendor/gsl-lite | 1 + 38 files changed, 2211 insertions(+), 88 deletions(-) create mode 100644 prints create mode 160000 vendor/gsl-lite diff --git a/.gitmodules b/.gitmodules index f84d09bb1f9..dc9b2c84373 100644 --- a/.gitmodules +++ b/.gitmodules @@ -13,3 +13,6 @@ [submodule "vendor/Catch2"] path = vendor/Catch2 url = https://github.com/catchorg/Catch2.git +[submodule "vendor/gsl-lite"] + path = vendor/gsl-lite + url = https://github.com/gsl-lite/gsl-lite.git diff --git a/include/openmc/constants.h b/include/openmc/constants.h index ae705607958..63137a148a4 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -290,7 +290,7 @@ enum class TallyResult { VALUE, SUM, SUM_SQ, SIZE }; enum class TallyType { VOLUME, MESH_SURFACE, SURFACE, PULSE_HEIGHT }; -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/distribution.h b/include/openmc/distribution.h index 854cf7d7719..4424f0e7b9f 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -23,7 +23,7 @@ class Distribution { public: virtual ~Distribution() = default; virtual double sample(uint64_t* seed) const = 0; - + virtual double get_pdf(double x) const = 0; //! Return integral of distribution //! \return Integral of distribution virtual double integral() const { return 1.0; }; @@ -84,6 +84,10 @@ class Discrete : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + //! Calculate the probability density function (PDF) at a given value + //! \param x The value at which to evaluate the PDF + //! \return The value of the PDF at the given point + double get_pdf(double x) const; double integral() const override { return di_.integral(); }; @@ -111,6 +115,7 @@ class Uniform : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; double a() const { return a_; } double b() const { return b_; } @@ -135,6 +140,7 @@ class PowerLaw : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; double a() const { return std::pow(offset_, ninv_); } double b() const { return std::pow(offset_ + span_, ninv_); } @@ -160,6 +166,7 @@ class Maxwell : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; double theta() const { return theta_; } @@ -180,6 +187,7 @@ class Watt : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; double a() const { return a_; } double b() const { return b_; } @@ -204,6 +212,7 @@ class Normal : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; double mean_value() const { return mean_value_; } double std_dev() const { return std_dev_; } @@ -227,8 +236,10 @@ class Tabular : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; - // properties + // double get_pdf_value(double x,uint64_t* seed) const; + // properties vector& x() { return x_; } const vector& x() const { return x_; } const vector& p() const { return p_; } @@ -263,6 +274,7 @@ class Equiprobable : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; const vector& x() const { return x_; } @@ -282,6 +294,7 @@ class Mixture : public Distribution { //! \param seed Pseudorandom number seed pointer //! \return Sampled value double sample(uint64_t* seed) const override; + double get_pdf(double x) const; double integral() const override { return integral_; } diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h index efd4e58425b..b0547ac6678 100644 --- a/include/openmc/distribution_angle.h +++ b/include/openmc/distribution_angle.h @@ -25,10 +25,11 @@ 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 get_pdf(double E, double mu, uint64_t* seed) const; //! Determine whether angle distribution is empty //! \return Whether distribution is empty bool empty() const { return energy_.empty(); } + double get_energy(int num) { return energy_[num]; } private: vector energy_; diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 0092c08f8df..880c352041e 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -114,7 +114,7 @@ void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name, void write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const long long* buffer, bool indep); void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen, - const char* name, char const* buffer, bool indep); + const char* name, const char* buffer, bool indep); void write_tally_results( hid_t group_id, hsize_t n_filter, hsize_t n_score, const double* results); } // extern "C" diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 0f37719b94f..d03b28f5c78 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -34,13 +34,14 @@ class Particle : public ParticleData { Particle() = default; - //========================================================================== - // Methods - double speed() const; //! create a secondary particle // + // + double getMass() const; + // + // //! stores the current phase space attributes of the particle in the //! secondary bank and increments the number of sites in the secondary bank. //! \param wgt Weight of the secondary particle @@ -63,7 +64,9 @@ class Particle : public ParticleData { //! simply as a secondary particle. //! \param src Source site data void from_source(const SourceSite* src); - + void initialize_ghost_particle(Particle& p, Direction u, double E); + void initialize_ghost_particle_from_source( + const SourceSite* src, Direction u_new); // Coarse-grained particle events void event_calculate_xs(); void event_advance(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index ec393845f81..3e657da2647 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -54,6 +54,8 @@ struct SourceSite { int parent_nuclide {-1}; int64_t parent_id; int64_t progeny_id; + int source_index; + bool ext {true}; //!< Is this source site external? }; //! State of a particle used for particle track files @@ -488,7 +490,9 @@ class ParticleData : public GeometryState { int event_nuclide_; int event_mt_; int delayed_group_ {0}; + int event_index_mt_; int parent_nuclide_ {-1}; + Direction v_t_; int n_bank_ {0}; double bank_second_E_ {0.0}; @@ -622,7 +626,10 @@ class ParticleData : public GeometryState { const int& event_nuclide() const { return event_nuclide_; } int& event_mt() { return event_mt_; } // MT number of collision const int& event_mt() const { return event_mt_; } + int& event_index_mt() { return event_index_mt_; } int& delayed_group() { return delayed_group_; } // delayed group + Position& v_t() { return v_t_; } + const Position& v_t() const { return v_t_; } const int& parent_nuclide() const { return parent_nuclide_; } int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide diff --git a/include/openmc/physics.h b/include/openmc/physics.h index f62f43a02f2..4f99ee01999 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -90,6 +90,11 @@ Direction sample_cxs_target_velocity( void sample_fission_neutron( int i_nuclide, const Reaction& rx, SourceSite* site, Particle& p); +void score_fission_neutron(int i_tally, int i_nuclide, const Reaction& rx, + SourceSite* site, Particle& p, std::vector& mu_cm, + std::vector& Js, std::vector& ghost_particles, + std::vector& pdfs_lab); + //! 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); diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index 4fbbc1b626a..3d2262c7100 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -11,8 +11,9 @@ #include "openmc/endf.h" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" +#include "openmc/tallies/filter.h" +#include "openmc/tallies/tally.h" #include "openmc/vector.h" // for vector - namespace openmc { //============================================================================== @@ -48,7 +49,10 @@ class ReactionProduct { //! \param[out] mu Outgoing cosine with respect to current direction //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const; - + void get_pdf(int i_tally, double E_in, double& E_out, uint64_t* seed, + Particle& p, std::vector& mu_cm, std::vector& Js, + std::vector& ghost_particles, + std::vector& pdfs_lab) 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..991f8d61138 100644 --- a/include/openmc/secondary_correlated.h +++ b/include/openmc/secondary_correlated.h @@ -10,6 +10,7 @@ #include "openmc/angle_energy.h" #include "openmc/distribution.h" #include "openmc/endf.h" +#include "openmc/particle.h" #include "openmc/vector.h" namespace openmc { @@ -40,6 +41,10 @@ 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; + void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed, + Particle& p, std::vector& mu_cm, std::vector& Js, + std::vector& ghost_particles, + std::vector& pdfs_lab) const; // energy property vector& energy() { return energy_; } diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index 83806d35248..ff5f13bd7dc 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -10,6 +10,7 @@ #include "openmc/angle_energy.h" #include "openmc/constants.h" #include "openmc/endf.h" +#include "openmc/particle.h" #include "openmc/vector.h" namespace openmc { @@ -31,6 +32,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 get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed, + Particle& p, std::vector& mu_cm, std::vector& Js, + std::vector& ghost_particles, + std::vector& pdfs_lab) const; 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..8395c68290d 100644 --- a/include/openmc/secondary_nbody.h +++ b/include/openmc/secondary_nbody.h @@ -7,6 +7,7 @@ #include "hdf5.h" #include "openmc/angle_energy.h" +#include "openmc/particle.h" namespace openmc { @@ -27,6 +28,10 @@ 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; + void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed, + Particle& p, std::vector& mu_cm, std::vector& Js, + std::vector& ghost_particles, + std::vector& pdfs_lab) const; private: int n_bodies_; //!< Number of particles distributed diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 5b18902afbb..7fb758ec570 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -32,6 +32,7 @@ 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 get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const; private: const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section @@ -55,6 +56,7 @@ 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 get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const; private: double debye_waller_; @@ -80,6 +82,7 @@ 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 get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const; private: const vector& energy_; //!< Energies at which cosines are tabulated @@ -106,6 +109,8 @@ 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; + double get_pdf( + double E_in, double& E_out, double& mu, uint64_t* seed, int l = -1) const; private: const vector& energy_; //!< Incident energies @@ -134,6 +139,7 @@ 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; + double get_pdf(double E_in, double& E_out, double& mu, uint64_t* seed) const; private: //! Secondary energy/angle distribution diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h index 3afa3d9ceb7..30d9c4b10df 100644 --- a/include/openmc/secondary_uncorrelated.h +++ b/include/openmc/secondary_uncorrelated.h @@ -10,6 +10,7 @@ #include "openmc/distribution_angle.h" #include "openmc/distribution_energy.h" #include "openmc/memory.h" +#include "openmc/particle.h" #include "openmc/vector.h" namespace openmc { @@ -31,6 +32,10 @@ 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; + void get_pdf(double det_pos[4], double E_in, double& E_out, uint64_t* seed, + Particle& p, std::vector& mu_cm, std::vector& Js, + std::vector& ghost_particles, + std::vector& pdfs_lab) const; // Accessors AngleDistribution& angle() { return angle_; } diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 48f678ced0f..8420cee5472 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -55,6 +55,10 @@ class Tally { void set_nuclides(const vector& nuclides); + void set_positions(pugi::xml_node node); + + void set_positions(const vector& positions); + const xt::xtensor& results() const { return results_; } //! returns vector of indices corresponding to the tally this is called on @@ -151,6 +155,9 @@ class Tally { vector scores_; //!< Filter integrands (e.g. flux, fission) + //! Index of each pos to be tallied. + vector positions_; + //! Index of each nuclide to be tallied. -1 indicates total material. vector nuclides_ {-1}; @@ -206,6 +213,7 @@ extern vector active_tallies; extern vector active_analog_tallies; extern vector active_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; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 28f1f16223d..eec09ff9cbf 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -67,8 +67,114 @@ class FilterBinIter { //! \param p The particle being tracked void score_collision_tally(Particle& p); -//! Score tallies based on a simple count of events (for continuous energy). -// +//! Score tallies using the point estimator. +//! This is triggered after every collision. +//! \param p The particle being tracked +void score_point_tally(Particle& p); +//! Retrieve the position and exclusion sphere (R0) of a detector. +//! +//! This function fetches the spatial position [x, y, z] and the exclusion +//! sphere radius (R0) of the point detector for a specified tally index. It +//! ensures that the tally contains exactly three spatial coordinates and an R0 +//! value. +//! +//! \param det_pos Array to store the detector position [x, y, z] and exclusion +//! sphere radius R0. \param i_tally Index of the tally for which the detector +//! position and R0 are required. +void get_det_pos(double (&det_pos)[4], int i_tally); +//! Score to point tally using data from a source site. +//! +//! This function calculates and scores flux for active point tallies +//! based on the properties of a given source site. It accounts for the source's +//! spatial position, directional distribution, and type to compute the +//! appropriate contribution to the tally. The function handles various source +//! angle distributions, including isotropic and polar-azimuthal types. +//! +//! \param src Pointer to the source site containing particle and source +//! information. \throws RuntimeError If an unexpected angle distribution type +//! is encountered. +void score_point_tally_from_source(const SourceSite* src); +//! Calculate the total mean free path (MFP) for a ghost particle along a +//! specified distance. +//! +//! This function computes the total mean free path (MFP) for a ghost particle +//! traveling from a collision point to a detector located at a given distance. +//! The direction of travel is towards the point detector. The cross-section +//! used in the calculation is computed based on the ghost particle's energy, +//! which remains constant throughout the traversal. The function accounts for +//! shielding by iteratively advancing the particle to boundaries and summing +//! contributions to the MFP along the path. +//! +//! \param ghost_particle The particle whose MFP is being calculated, with its +//! direction towards the point detector. \param total_distance The distance +//! [cm] from the collision point to the detector. \return The total mean free +//! path for the particle over the specified distance. \throws RuntimeError If +//! the particle encounters unexpected behavior during boundary crossing. +double get_MFP(Particle& ghost_particle, double total_distance); +//! Calculate the probability density function (PDF) for elastic and inelastic +//! scattering to a point detector. +//! +//! This function computes the PDF for a particle scattering (elastic or +//! inelastic) to a detector at a specified position. All computations are +//! performed fully relativistically, including energy, momentum, and +//! transformations between the center-of-mass (CM) frame and the lab frame. The +//! Jacobian used in the calculations corresponds to a Lorentz transformation, +//! ensuring accurate mapping between frames of reference. +//! +//! Ghost particles are created to represent possible scattering trajectories. +//! Each ghost particle is initialized with the correct energy and direction as +//! if the particle had scattered towards the point detector. +//! +//! In the case of inelastic scattering, the post-collision energy of the +//! outgoing neutron is set in the CM frame, leading to an effective mass +//! (\(m_4\)) for the residual nucleus. For inelastic collisions, +//! \(m_4 \neq m_2\), as \(m_4\) reflects the excitation state of the residual +//! nucleus. + +//! \param det_pos The position of the point detector in [x, y, z, R0]. +//! \param p The incident particle being analyzed for scattering. +//! \param mu_cm A vector to store the cosine of the scattering angles in the CM +//! frame. \param Js A vector to store the Jacobians for transforming the +//! distributions to the lab frame using Lorentz transformations. \param +//! ghost_particles A vector to store the ghost particles representing +//! scattering trajectories. +//! Each ghost particle is initialized with the correct +//! energy and direction towards the point detector. +//! \param E3k_cm_given The center-of-mass kinetic energy of one outgoing +//! particle, if specified. +//! If set to -1, it is calculated based on kinematics. +void get_pdf_to_point_elastic(double det_pos[4], Particle& p, + std::vector& mu_cm, std::vector& Js, + std::vector& ghost_particles, double E3k_cm_given = -1); +//! Score a ghost particle contribution to a point tally. +//! +//! This function calculates and scores the contribution of a ghost particle +//! to a specified point tally. The scoring process includes computing the flux +//! at the detector location. +//! It accounts for cases where the particle path is fully or partially within +//! the exclusion sphere (\(R_0\)) or outside it. The scoring respects particle +//! type, ensuring only neutrons are considered. +//! +//! \param ghost_p The ghost particle being scored. +//! \param pdf_lab The probability density function value in the lab frame for +//! the particle direction. \param i_tally The index of the tally to which the +//! ghost particle contribution is scored. +//! +//! \note Only continuous energy mode is currently supported; multi-group mode +//! is not implemented. +void score_ghost_particle(Particle& ghost_p, double pdf_lab, int i_tally); +//! Perform a Lorentz boost of a four-vector. +//! +//! This function boosts a four-vector \( B \) (in the lab frame) into the rest +//! frame defined by another four-vector \( A \). The result of the boost is +//! stored in \( X \), representing \( B \) in the rest frame of \( A \). +//! +//! \param A The four-vector defining the rest frame [E, px, py, pz]. +//! \param B The four-vector to be boosted [E, px, py, pz]. +//! \param X The resulting four-vector after the Lorentz boost [E, px, py, pz] +//! (rest frame). +void boostf(double A[4], double B[4], double X[4]); + //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked @@ -76,6 +182,10 @@ void score_analog_tally_ce(Particle& p); //! Score tallies based on a simple count of events (for multigroup). // +void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, + int filter_index, double filter_weight, int i_nuclide, double atom_density, + double flux); + //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index de0767d0af0..d6291c30e5e 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -52,6 +52,8 @@ class ThermalData { //! \param[inout] seed Pseudorandom seed pointer void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out, double* mu, uint64_t* seed); + double get_pdf(const NuclideMicroXS& micro_xs, double E, double& E_out, + double mu, uint64_t* seed); private: struct Reaction { diff --git a/include/openmc/urr.h b/include/openmc/urr.h index 3978c7b86a2..1e603715847 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -27,10 +27,10 @@ class UrrData { double heating; }; - Interpolation interp_; //!< interpolation type - int inelastic_flag_; //!< inelastic competition flag - int absorption_flag_; //!< other absorption flag - bool multiply_smooth_; //!< multiply by smooth cross section? + Interpolation interp_; //!< interpolation type + int inelastic_flag_; //!< inelastic competition flag + int absorption_flag_; //!< other absorption flag + bool multiply_smooth_; //!< multiply by smooth cross section? vector energy_; //!< incident energies auto n_energy() const { return energy_.size(); } diff --git a/openmc/tallies.py b/openmc/tallies.py index 06de0205c9e..2626140e88a 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -7,7 +7,7 @@ import operator from pathlib import Path import lxml.etree as ET - +from re import sub import h5py import numpy as np import pandas as pd @@ -34,7 +34,7 @@ _FILTER_CLASSES = (openmc.Filter, openmc.CrossFilter, openmc.AggregateFilter) # Valid types of estimators -ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog'} +ESTIMATOR_TYPES = ['tracklength', 'collision', 'analog', 'point'] class Tally(IDManagerMixin): @@ -58,7 +58,7 @@ class Tally(IDManagerMixin): multiply_density : bool Whether reaction rates should be multiplied by atom density - .. versionadded:: 0.14.0 + .. versionadded:: 0.13.4 filters : list of openmc.Filter List of specified filters for the tally nuclides : list of str @@ -135,7 +135,7 @@ def __init__(self, tally_id=None, name=''): self._with_batch_statistics = False self._derived = False self._sparse = False - + self._positions = None self._sp_filename = None self._results_read = False @@ -174,14 +174,18 @@ def __repr__(self): parts.append('{: <15}=\t{}'.format('ID', self.id)) parts.append('{: <15}=\t{}'.format('Name', self.name)) if self.derivative is not None: - parts.append('{: <15}=\t{}'.format('Derivative ID', self.derivative.id)) + parts.append('{: <15}=\t{}'.format( + 'Derivative ID', self.derivative.id)) filters = ', '.join(type(f).__name__ for f in self.filters) parts.append('{: <15}=\t{}'.format('Filters', filters)) nuclides = ' '.join(str(nuclide) for nuclide in self.nuclides) parts.append('{: <15}=\t{}'.format('Nuclides', nuclides)) parts.append('{: <15}=\t{}'.format('Scores', self.scores)) parts.append('{: <15}=\t{}'.format('Estimator', self.estimator)) - parts.append('{: <15}=\t{}'.format('Multiply dens.', self.multiply_density)) + parts.append('{: <15}=\t{}'.format( + 'Multiply dens.', self.multiply_density)) + if self.positions is not None: + parts.append('{: <15}=\t{}'.format('positions', self.positions)) return '\n\t'.join(parts) @staticmethod @@ -212,6 +216,14 @@ def name(self, name): cv.check_type('tally name', name, str, none_ok=True) self._name = name + @property + def positions(self): + return self._positions + + @positions.setter + def positions(self, positions): + self._positions = positions + @property def multiply_density(self): return self._multiply_density @@ -239,7 +251,8 @@ def filters(self, filters): raise ValueError(msg) visited_filters.add(f) - self._filters = cv.CheckedList(_FILTER_CLASSES, 'tally filters', filters) + self._filters = cv.CheckedList( + _FILTER_CLASSES, 'tally filters', filters) @property @ensure_results @@ -290,7 +303,7 @@ def scores(self, scores): # Check to see if scores are deprecated before storing for deprecated in ['scatter-', 'nu-scatter-', 'scatter-p', 'nu-scatter-p', 'scatter-y', 'nu-scatter-y', - 'flux-y', 'total-y']: + 'flux-y', 'total-y']: if score.strip().startswith(deprecated): msg = score.strip() + ' is no longer supported.' raise ValueError(msg) @@ -346,7 +359,8 @@ def num_realizations(self): @num_realizations.setter def num_realizations(self, num_realizations): cv.check_type('number of realizations', num_realizations, Integral) - cv.check_greater_than('number of realizations', num_realizations, 0, True) + cv.check_greater_than('number of realizations', + num_realizations, 0, True) self._num_realizations = num_realizations @property @@ -387,8 +401,10 @@ def _read_results(self): # Convert NumPy arrays to SciPy sparse LIL matrices if self.sparse: - self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape) - self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) + self._sum = sps.lil_matrix( + self._sum.flatten(), self._sum.shape) + self._sum_sq = sps.lil_matrix( + self._sum_sq.flatten(), self._sum_sq.shape) # Read simulation time (needed for figure of merit) self._simulation_time = f["runtime"]["simulation"][()] @@ -524,7 +540,8 @@ def sparse(self, sparse): # Convert NumPy arrays to SciPy sparse LIL matrices if sparse and not self.sparse: if self._sum is not None: - self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape) + self._sum = sps.lil_matrix( + self._sum.flatten(), self._sum.shape) if self._sum_sq is not None: self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape) @@ -748,7 +765,8 @@ def can_merge(self, other): can_merge_filters = self._can_merge_filters(other) can_merge_nuclides = self._can_merge_nuclides(other) can_merge_scores = self._can_merge_scores(other) - mergeability = [can_merge_filters, can_merge_nuclides, can_merge_scores] + mergeability = [can_merge_filters, + can_merge_nuclides, can_merge_scores] if not all(mergeability): return False @@ -804,7 +822,8 @@ def merge(self, other): for i, filter1 in enumerate(self.filters): for filter2 in other.filters: if filter1 != filter2 and filter1.can_merge(filter2): - other_copy._swap_filters(other_copy.filters[i], filter2) + other_copy._swap_filters( + other_copy.filters[i], filter2) merged_tally.filters[i] = filter1.merge(filter2) join_right = filter1 < filter2 merge_axis = i @@ -867,7 +886,8 @@ def merge(self, other): merged_sum_sq = np.concatenate((other_sum_sq, self_sum_sq), axis=merge_axis) - merged_tally._sum_sq = np.reshape(merged_sum_sq, merged_tally.shape) + merged_tally._sum_sq = np.reshape( + merged_sum_sq, merged_tally.shape) # Concatenate mean arrays if present in both tallies if self.mean is not None and other.mean is not None: @@ -895,7 +915,8 @@ def merge(self, other): merged_std_dev = np.concatenate((other_std_dev, self_std_dev), axis=merge_axis) - merged_tally._std_dev = np.reshape(merged_std_dev, merged_tally.shape) + merged_tally._std_dev = np.reshape( + merged_std_dev, merged_tally.shape) # Sparsify merged tally if both tallies are sparse merged_tally.sparse = self.sparse and other.sparse @@ -949,6 +970,10 @@ def to_xml_element(self): subelement = ET.SubElement(element, "estimator") subelement.text = self.estimator + if self.positions is not None: + subelement = ET.SubElement(element, "positions") + subelement.text = ' '.join(str(x) for x in self.positions) + # Optional Triggers for trigger in self.triggers: element.append(trigger.to_xml_element()) @@ -1035,6 +1060,10 @@ def from_xml_element(cls, elem, **kwargs): if estimator_elem is not None: tally.estimator = estimator_elem.text + positions_elem = elem.find('positions') + if positions_elem is not None: + tally.positions = positions_elem.text + # Read triggers tally.triggers = [ openmc.Trigger.from_xml_element(trigger_elem) @@ -1215,7 +1244,8 @@ def get_filter_indices(self, filters=[], filter_bins=[]): for j, test_filter in enumerate(filters): if type(self_filter) is test_filter: bins = filter_bins[j] - indices = np.array([self_filter.get_bin_index(b) for b in bins]) + indices = np.array( + [self_filter.get_bin_index(b) for b in bins]) break else: indices = np.arange(self_filter.num_bins) @@ -1720,17 +1750,19 @@ def hybrid_product(self, other, binary_op, filter_product=None, elif binary_op == '*': with np.errstate(divide='ignore', invalid='ignore'): self_rel_err = data['self']['std. dev.'] / data['self']['mean'] - other_rel_err = data['other']['std. dev.'] / data['other']['mean'] + other_rel_err = data['other']['std. dev.'] / \ + data['other']['mean'] new_tally._mean = data['self']['mean'] * data['other']['mean'] new_tally._std_dev = np.abs(new_tally.mean) * \ - np.sqrt(self_rel_err**2 + other_rel_err**2) + np.sqrt(self_rel_err**2 + other_rel_err**2) elif binary_op == '/': with np.errstate(divide='ignore', invalid='ignore'): self_rel_err = data['self']['std. dev.'] / data['self']['mean'] - other_rel_err = data['other']['std. dev.'] / data['other']['mean'] + other_rel_err = data['other']['std. dev.'] / \ + data['other']['mean'] new_tally._mean = data['self']['mean'] / data['other']['mean'] new_tally._std_dev = np.abs(new_tally.mean) * \ - np.sqrt(self_rel_err**2 + other_rel_err**2) + np.sqrt(self_rel_err**2 + other_rel_err**2) elif binary_op == '^': with np.errstate(divide='ignore', invalid='ignore'): mean_ratio = data['other']['mean'] / data['self']['mean'] @@ -1739,7 +1771,7 @@ def hybrid_product(self, other, binary_op, filter_product=None, np.log(data['self']['mean']) * data['other']['std. dev.'] new_tally._mean = data['self']['mean'] ** data['other']['mean'] new_tally._std_dev = np.abs(new_tally.mean) * \ - np.sqrt(first_term**2 + second_term**2) + np.sqrt(first_term**2 + second_term**2) # Convert any infs and nans to zero new_tally._mean[np.isinf(new_tally._mean)] = 0 @@ -1853,14 +1885,16 @@ def _align_tally_data(self, other, filter_product, nuclide_product, for other_filter in other_missing_filters: filter_copy = copy.deepcopy(other_filter) other._mean = np.repeat(other.mean, filter_copy.num_bins, axis=0) - other._std_dev = np.repeat(other.std_dev, filter_copy.num_bins, axis=0) + other._std_dev = np.repeat( + other.std_dev, filter_copy.num_bins, axis=0) other.filters.append(filter_copy) # Add filters present in other but not in self to self for self_filter in self_missing_filters: filter_copy = copy.deepcopy(self_filter) self._mean = np.repeat(self.mean, filter_copy.num_bins, axis=0) - self._std_dev = np.repeat(self.std_dev, filter_copy.num_bins, axis=0) + self._std_dev = np.repeat( + self.std_dev, filter_copy.num_bins, axis=0) self.filters.append(filter_copy) # Align other filters with self filters @@ -1890,7 +1924,8 @@ def _align_tally_data(self, other, filter_product, nuclide_product, # Add nuclides present in self but not in other to other for nuclide in other_missing_nuclides: - other._mean = np.insert(other.mean, other.num_nuclides, 0, axis=1) + other._mean = np.insert( + other.mean, other.num_nuclides, 0, axis=1) other._std_dev = np.insert(other.std_dev, other.num_nuclides, 0, axis=1) other.nuclides.append(nuclide) @@ -1929,14 +1964,17 @@ def _align_tally_data(self, other, filter_product, nuclide_product, # Add scores present in self but not in other to other for score in other_missing_scores: - other._mean = np.insert(other.mean, other.num_scores, 0, axis=2) - other._std_dev = np.insert(other.std_dev, other.num_scores, 0, axis=2) + other._mean = np.insert( + other.mean, other.num_scores, 0, axis=2) + other._std_dev = np.insert( + other.std_dev, other.num_scores, 0, axis=2) other.scores.append(score) # Add scores present in other but not in self to self for score in self_missing_scores: self._mean = np.insert(self.mean, self.num_scores, 0, axis=2) - self._std_dev = np.insert(self.std_dev, self.num_scores, 0, axis=2) + self._std_dev = np.insert( + self.std_dev, self.num_scores, 0, axis=2) self.scores.append(score) # Align other scores with self scores @@ -2880,7 +2918,7 @@ def summation(self, scores=[], filter_type=None, # Add AggregateFilter to the tally sum if not remove_filter: filter_sum = openmc.AggregateFilter(self_filter, - [tuple(filter_bins)], 'sum') + [tuple(filter_bins)], 'sum') tally_sum.filters.append(filter_sum) # Add a copy of each filter not summed across to the tally sum @@ -2893,7 +2931,8 @@ def summation(self, scores=[], filter_type=None, # Sum across any nuclides specified by the user if len(nuclides) != 0: - nuclide_bins = [self.get_nuclide_index(nuclide) for nuclide in nuclides] + nuclide_bins = [self.get_nuclide_index( + nuclide) for nuclide in nuclides] axis_index = self.num_filters mean = np.take(mean, indices=nuclide_bins, axis=axis_index) std_dev = np.take(std_dev, indices=nuclide_bins, axis=axis_index) @@ -3033,7 +3072,7 @@ def average(self, scores=[], filter_type=None, # Add AggregateFilter to the tally avg if not remove_filter: filter_sum = openmc.AggregateFilter(self_filter, - [tuple(filter_bins)], 'avg') + [tuple(filter_bins)], 'avg') tally_avg.filters.append(filter_sum) # Add a copy of each filter not averaged across to the tally avg @@ -3046,7 +3085,8 @@ def average(self, scores=[], filter_type=None, # Sum across any nuclides specified by the user if len(nuclides) != 0: - nuclide_bins = [self.get_nuclide_index(nuclide) for nuclide in nuclides] + nuclide_bins = [self.get_nuclide_index( + nuclide) for nuclide in nuclides] axis_index = self.num_filters mean = np.take(mean, indices=nuclide_bins, axis=axis_index) std_dev = np.take(std_dev, indices=nuclide_bins, axis=axis_index) diff --git a/prints b/prints new file mode 100644 index 00000000000..b72c9356daf --- /dev/null +++ b/prints @@ -0,0 +1,36 @@ + fmt::print("post energy = {0} pre energy = {1}\n",p.E(),p.E_last()); + fmt::print("event mt = {}\n",p.event_mt()); + fmt::print("post weight = {0} pre weight = {1}\n",p.wgt(),p.wgt_last()); + fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z); + fmt::print("u = {0} , {1} , {2}\n",p.u().x,p.u().y,p.u().z); + fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z); + fmt::print("distance = {}\n",distance); + fmt::print("xs = {}\n",p.macro_xs().total); + fmt::print("number of nuclear = {}\n",n); + fmt::print("name of 1# nuclicde = {}\n",nuc->name_); + fmt::print("name of 2# nuclicde = {}\n",nuc2->name_); + fmt::print("reaction MT = {}\n",rx->mt_); + fmt::print("d type = {}\n",typeid(d.get()).name()); + fmt::print("energy distribution at 0 {}\n",d_->angle().get_energy(0)); + fmt::print("energy distribution at 1 {}\n",d_->angle().get_energy(1)); + fmt::print("post energy = {}\n",p.E); + fmt::print(p.E); + + // Neutron velocity in LAB + Direction v_n = vel * p.u_last(); + + // // Sample velocity of target nucleus + // Direction v_t {}; + // if (!p.neutron_xs(i_nuclide).use_ptable) { + // v_t = sample_target_velocity(*nuc, p.E_last(), p.u_last(), v_n, + // p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); + // } + + // Velocity of center-of-mass + Direction v_cm = (v_n ) / (awr + 1.0); //(v_n + awr * v_t) / (awr + 1.0); + + // Transform to CM frame + v_n -= v_cm; + + // Find speed of neutron in CM + vel = v_n.norm(); diff --git a/src/distribution.cpp b/src/distribution.cpp index ca00cbde4eb..1f0eb29c88b 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -12,6 +12,7 @@ #include "openmc/math_functions.h" #include "openmc/random_dist.h" #include "openmc/random_lcg.h" +#include "openmc/search.h" #include "openmc/xml_interface.h" namespace openmc { @@ -134,6 +135,12 @@ double Discrete::sample(uint64_t* seed) const { return x_[di_.sample(seed)]; } +// Not implemented +double Discrete::get_pdf(double x) const +{ + + return -1; +} //============================================================================== // Uniform implementation @@ -155,6 +162,13 @@ double Uniform::sample(uint64_t* seed) const { return a_ + prn(seed) * (b_ - a_); } +double Uniform::get_pdf(double x) const +{ + if (x <= b_ && x >= a_) + return 1 / (b_ - a_); + else + return 0; +} //============================================================================== // PowerLaw implementation @@ -181,6 +195,10 @@ double PowerLaw::sample(uint64_t* seed) const { return std::pow(offset_ + prn(seed) * span_, ninv_); } +double PowerLaw::get_pdf(double x) const +{ + return x / span_ / ninv_; +} //============================================================================== // Maxwell implementation @@ -195,6 +213,11 @@ double Maxwell::sample(uint64_t* seed) const { return maxwell_spectrum(theta_, seed); } +// Not implemented +double Maxwell::get_pdf(double x) const +{ + return -1; +} //============================================================================== // Watt implementation @@ -215,6 +238,11 @@ double Watt::sample(uint64_t* seed) const { return watt_spectrum(a_, b_, seed); } +// Not implemented +double Watt::get_pdf(double x) const +{ + return -1; +} //============================================================================== // Normal implementation @@ -235,6 +263,12 @@ double Normal::sample(uint64_t* seed) const { return normal_variate(mean_value_, std_dev_, seed); } +double Normal::get_pdf(double x) const +{ + double exponent = -0.5 * std::pow((x - mean_value_) / std_dev_, 2); + double coefficient = 1 / (std_dev_ * std::sqrt(2 * PI)); + return coefficient * std::exp(exponent); +} //============================================================================== // Tabular implementation @@ -356,6 +390,41 @@ double Tabular::sample(uint64_t* seed) const } } +double Tabular::get_pdf(double x) const +{ + // get PDF value at x + + int i; + std::size_t n = x_.size(); + if (x < x_[0]) { + return 0; + } else if (x > x_[n - 1]) { + return 0; + } else { + i = lower_bound_index(x_.begin(), x_.end(), x); + } + + // Determine bounding PDF values + double x_i = x_[i]; + double p_i = p_[i]; + + if (interp_ == Interpolation::histogram) { + // Histogram interpolation + return p_i; + } else { + // Linear-linear interpolation + double x_i1 = x_[i + 1]; + double p_i1 = p_[i + 1]; + + double m = (p_i1 - p_i) / (x_i1 - x_i); + if (m == 0.0) { + return p_i; + } else { + return p_i + (x - x_i) * m; + } + } +} + //============================================================================== // Equiprobable implementation //============================================================================== @@ -371,6 +440,10 @@ double Equiprobable::sample(uint64_t* seed) const double xr = x_[i + i]; return xl + ((n - 1) * r - i) * (xr - xl); } +double Equiprobable::get_pdf(double x) const +{ + return -1; +} //============================================================================== // Mixture implementation @@ -420,6 +493,10 @@ double Mixture::sample(uint64_t* seed) const // Sample the chosen distribution return it->second->sample(seed); } +double Mixture::get_pdf(double x) const +{ + return -1; +} //============================================================================== // Helper function diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp index 5e92b794a72..5c35ff9497a 100644 --- a/src/distribution_angle.cpp +++ b/src/distribution_angle.cpp @@ -95,4 +95,37 @@ double AngleDistribution::sample(double E, uint64_t* seed) const return mu; } +double AngleDistribution::get_pdf(double E, double mu, uint64_t* seed) const +{ + // Determine number of incoming energies + auto n = energy_.size(); + + // Find energy bin and calculate interpolation factor -- if the energy is + // outside the range of the tabulated energies, choose the first or last bins + int i; + double r; + if (E < energy_[0]) { + i = 0; + r = 0.0; + } else if (E > energy_[n - 1]) { + i = n - 2; + r = 1.0; + } else { + i = lower_bound_index(energy_.begin(), energy_.end(), E); + r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]); + } + + // Sample between the ith and (i+1)th bin + if (r > prn(seed)) + ++i; + + // Sample i-th distribution + double pdf = distribution_[i]->get_pdf(mu); + + // Make sure pdf > 0 and return + if (std::abs(pdf) < 0) + fmt::print("pdf = {} <1", pdf); + return pdf; +} + } // namespace openmc diff --git a/src/output.cpp b/src/output.cpp index 0a14e8843d8..64a8505c913 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -48,29 +48,30 @@ namespace openmc { void title() { - fmt::print(" %%%%%%%%%%%%%%%\n" - " %%%%%%%%%%%%%%%%%%%%%%%%\n" - " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" - " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" - " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" - " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" - " %%%%%%%%%%%%%%%%%%%%%%%%\n" - " %%%%%%%%%%%%%%%%%%%%%%%%\n" - " ############### %%%%%%%%%%%%%%%%%%%%%%%%\n" - " ################## %%%%%%%%%%%%%%%%%%%%%%%\n" - " ################### %%%%%%%%%%%%%%%%%%%%%%%\n" - " #################### %%%%%%%%%%%%%%%%%%%%%%\n" - " ##################### %%%%%%%%%%%%%%%%%%%%%\n" - " ###################### %%%%%%%%%%%%%%%%%%%%\n" - " ####################### %%%%%%%%%%%%%%%%%%\n" - " ####################### %%%%%%%%%%%%%%%%%\n" - " ###################### %%%%%%%%%%%%%%%%%\n" - " #################### %%%%%%%%%%%%%%%%%\n" - " ################# %%%%%%%%%%%%%%%%%\n" - " ############### %%%%%%%%%%%%%%%%\n" - " ############ %%%%%%%%%%%%%%%\n" - " ######## %%%%%%%%%%%%%%\n" - " %%%%%%%%%%%\n\n"); + fmt::print( + " %%%%%%%%%%%%%%%\n" + " %%%%%%%%%%%%%%%%%%%%%%%%\n" + " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" + " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" + " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" + " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n" + " %%%%%%%%%%%%%%%%%%%%%%%%\n" + " %%%%%%%%%%%%%%%%%%%%%%%%\n" + " ############### %%%%%%%%%%%%%%%%%%%%%%%%\n" + " ################## %%%%%%%%%%%%%%%%%%%%%%%\n" + " ################### %%%%%%%%%%%%%%%%%%%%%%%\n" + " #################### %%%%%%%%%%%%%%%%%%%%%%%\n" + " ##################### %%%%%%%%%%%%%%%%%%%%%%%\n" + " ###################### %%%%%%%%%%%%%%%%%%%%\n" + " ####################### %%%%%%%%%%%%%%%%%%\n" + " ####################### %%%%%%%%%%%%%%%%%\n" + " ###################### %%%%%%%%%%%%%%%%%\n" + " #################### %%%%%%%%%%%%%%%%%\n" + " ################# %%%%%%%%%%%%%%%%%\n" + " ############### %%%%%%%%%%%%%%%%\n" + " ############ %%%%%%%%%%%%%%%\n" + " ######## %%%%%%%%%%%%%%\n" + " %%%%%%%%%%%\n\n"); // Write version information fmt::print( diff --git a/src/particle.cpp b/src/particle.cpp index 118b084b04e..94cded90d07 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -74,6 +74,90 @@ double Particle::speed() const return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma); } } +double Particle::getMass() const +{ + // 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; + } + + return mass; +} + +void Particle::initialize_ghost_particle( + Particle& p, Direction u_new, double E_new) +{ + clear(); + surface() = 0; + cell_born() = C_NONE; + material() = C_NONE; + n_collision() = p.n_collision(); + fission() = false; + zero_flux_derivs(); + + type() = p.type(); + wgt() = p.wgt(); + wgt_last() = p.wgt_last(); + r() = p.r(); + r_last() = p.r(); + u() = u_new; + u_last() = u_new; + // u_last() = p.u_last(); + E() = E_new; // settings::run_CE ? p.E_last() : p.g_last(); + E_last() = E_new; + time() = p.time_last(); + time_last() = p.time_last(); +} + +void Particle::initialize_ghost_particle_from_source( + const SourceSite* src, Direction u_new) +{ + // Reset some attributes + clear(); + surface() = 0; + cell_born() = C_NONE; + material() = C_NONE; + n_collision() = 0; + fission() = false; + zero_flux_derivs(); + + // Copy attributes from source bank site + type() = src->particle; + wgt() = src->wgt; + wgt_last() = src->wgt; + r() = src->r; + u() = u_new; + r_last_current() = src->r; + r_last() = src->r; + u_last() = src->u; + if (settings::run_CE) { + E() = src->E; + g() = 0; + } else { + g() = static_cast(src->E); + g_last() = static_cast(src->E); + E() = data::mg.energy_bin_avg_[g()]; + } + E_last() = E(); + time() = src->time; + time_last() = src->time; + + // fmt::print("==============================particle + // created==============================\n"); fmt::print("u = {0} , {1} , + // {2}\n",u().x,u().y,u().z); fmt::print("u_last = {0} , {1} , + // {2}\n",u_last().x,u_last().y,u_last().z); fmt::print("pos = {0} , {1} , + // {2}\n",r().x,r().y,r().z); +} bool Particle::create_secondary( double wgt, Direction u, double E, ParticleType type) @@ -154,6 +238,14 @@ void Particle::from_source(const SourceSite* src) int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1; surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one; } + if (!model::active_point_tallies.empty()) + score_point_tally_from_source(src); + + // fmt::print("==============================particle + // created==============================\n"); fmt::print("u = {0} , {1} , + // {2}\n",u().x,u().y,u().z); fmt::print("u_last = {0} , {1} , + // {2}\n",u_last().x,u_last().y,u_last().z); fmt::print("pos = {0} , {1} , + // {2}\n",r().x,r().y,r().z); } void Particle::event_calculate_xs() @@ -359,8 +451,11 @@ void Particle::event_collide() // Score collision estimator tallies -- this is done after a collision // has occurred rather than before because we need information on the // outgoing energy for any tallies with an outgoing energy filter + if (!model::active_collision_tallies.empty()) score_collision_tally(*this); + if (!model::active_point_tallies.empty() && exhaustive_find_cell(*this)) + score_point_tally(*this); if (!model::active_analog_tallies.empty()) { if (settings::run_CE) { score_analog_tally_ce(*this); diff --git a/src/particle_data.cpp b/src/particle_data.cpp index 370ca12e469..efa19c4b8f0 100644 --- a/src/particle_data.cpp +++ b/src/particle_data.cpp @@ -47,6 +47,13 @@ void LocalCoord::reset() lattice_index_[2] = 0; rotated_ = false; } +// void ParticleData::reset_cords() +//{ +// Create and clear coordinate levels +// coord_.resize(model::n_coord_levels); +// cell_last_.resize(model::n_coord_levels); +// clear(); +//} GeometryState::GeometryState() { diff --git a/src/physics.cpp b/src/physics.cpp index 3c06e543de1..af894e0febb 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -1,4 +1,5 @@ #include "openmc/physics.h" +#include #include "openmc/bank.h" #include "openmc/bremsstrahlung.h" @@ -25,9 +26,9 @@ #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" - #include #include // for max, min, max_element @@ -202,6 +203,11 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Counter for the number of fission sites successfully stored to the shared // fission bank or the secondary particle bank int n_sites_stored; + // Initialize for point detector + std::vector ghost_particles; + std::vector mu_cm; + std::vector Js; + std::vector pdfs_lab; for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { // Initialize fission site object with particle data @@ -216,7 +222,11 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Sample delayed group and angle/energy for fission reaction sample_fission_neutron(i_nuclide, rx, &site, p); - + p.event_index_mt() = -999; + for (auto i_tally : model::active_point_tallies) { + score_fission_neutron( + i_tally, i_nuclide, rx, &site, p, mu_cm, Js, ghost_particles, pdfs_lab); + } // Store fission site in bank if (use_fission_bank) { int64_t idx = simulation::fission_bank.thread_safe_append(site); @@ -665,7 +675,7 @@ void scatter(Particle& p, int i_nuclide) { // copy incoming direction Direction u_old {p.u()}; - + p.event_index_mt() = 0; // Get pointer to nuclide and grid index/interpolation factor const auto& nuc {data::nuclides[i_nuclide]}; const auto& micro {p.neutron_xs(i_nuclide)}; @@ -702,7 +712,7 @@ void scatter(Particle& p, int i_nuclide) // S(A,B) SCATTERING sab_scatter(i_nuclide, micro.index_sab, p); - + p.event_index_mt() = -1234; // to distinguish from elastic p.event_mt() = ELASTIC; sampled = true; } @@ -724,6 +734,7 @@ void scatter(Particle& p, int i_nuclide) const auto& rx {nuc->reactions_[i]}; inelastic_scatter(*nuc, *rx, p); p.event_mt() = rx->mt_; + p.event_index_mt() = i; } // Set event component @@ -758,8 +769,8 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p) v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n, p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); } - - // Velocity of center-of-mass + p.v_t() = C_LIGHT * std::sqrt(2 / p.getMass()) * v_t; + // Velocity of center-of-mass Direction v_cm = (v_n + awr * v_t) / (awr + 1.0); // Transform to CM frame @@ -797,6 +808,11 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p) // neutron's pre- and post-collision angle p.mu() = p.u().dot(v_n) / vel; + // std::ofstream outfile; + + // outfile.open("scatter", std::ios_base::app); // append instead of overwrite + // outfile << v_cm.dot(v_n)/v_cm.norm()/v_n.norm()<<"\n"; + // Set energy and direction of particle in LAB frame p.u() = v_n / vel; @@ -818,6 +834,22 @@ void sab_scatter(int i_nuclide, int i_sab, Particle& p) data::thermal_scatt[i_sab]->data_[i_temp].sample( micro, p.E(), &E_out, &p.mu(), p.current_seed()); + for (auto i_tally : model::active_point_tallies) { + double E_out_ghost; + double det_pos[4]; + get_det_pos(det_pos, i_tally); + Direction u_lab {det_pos[0] - p.r().x, // towards the detector + det_pos[1] - p.r().y, det_pos[2] - p.r().z}; + Direction u_lab_unit = u_lab / u_lab.norm(); // normalize + double mu_det = + u_lab_unit.dot(p.u_last()); // target velocity is treated as zero + double pdf = data::thermal_scatt[i_sab]->data_[i_temp].get_pdf( + micro, p.E(), E_out_ghost, mu_det, p.current_seed()); + Particle ghost_particle = Particle(); + ghost_particle.initialize_ghost_particle(p, u_lab_unit, E_out_ghost); + score_ghost_particle(ghost_particle, pdf, i_tally); + } + // Set energy to outgoing, change direction of particle p.E() = E_out; p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); @@ -854,12 +886,13 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, // use appropriate target velocity sampling method switch (sampling_method) { case ResScatMethod::cxs: - - // sample target velocity with the constant cross section (cxs) approx. + // fmt::print("ResScatMethod = ResScatMethod::cxs\n"); + // sample target velocity with the constant cross section (cxs) approx. return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); case ResScatMethod::dbrc: case ResScatMethod::rvs: { + fmt::print("ResScatMethod = ResScatMethod::rvs\n"); double E_red = std::sqrt(nuc.awr_ * E / kT); double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc.awr_; double E_up = (E_red + 4.0) * (E_red + 4.0) * kT / nuc.awr_; @@ -1102,8 +1135,96 @@ void sample_fission_neutron( site->u = rotate_angle(p.u(), mu, nullptr, seed); } +void score_fission_neutron(int i_tally, int i_nuclide, const Reaction& rx, + SourceSite* site, Particle& p, std::vector& mu_cm, + std::vector& Js, std::vector& ghost_particles, + std::vector& pdfs_lab) +{ + + // Get attributes of particle + double E_in = p.E(); + uint64_t* seed = p.current_seed(); + + // Determine total nu, delayed nu, and delayed neutron fraction + const auto& nuc {data::nuclides[i_nuclide]}; + double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total); + double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed); + double beta = nu_d / nu_t; + + if (prn(seed) < beta) { + // ==================================================================== + // DELAYED NEUTRON SAMPLED + + // sampled delayed precursor group + double xi = prn(seed) * nu_d; + double prob = 0.0; + int group; + for (group = 1; group < nuc->n_precursor_; ++group) { + // determine delayed neutron precursor yield for group j + double yield = (*rx.products_[group].yield_)(E_in); + + // Check if this group is sampled + prob += yield; + if (xi < prob) + break; + } + + // if the sum of the probabilities is slightly less than one and the + // random number is greater, j will be greater than nuc % + // n_precursor -- check for this condition + group = std::min(group, nuc->n_precursor_); + + // set the delayed group for the particle born from fission + site->delayed_group = group; + + } else { + // ==================================================================== + // PROMPT NEUTRON SAMPLED + + // set the delayed group for the particle born from fission to 0 + site->delayed_group = 0; + } + + // sample from prompt neutron energy distribution + int n_sample = 0; + double mu; + while (true) { + // erase previous elemts bc they were rejected. + mu_cm.clear(); + Js.clear(); + pdfs_lab.clear(); + ghost_particles.clear(); + double E_out; + rx.products_[site->delayed_group].get_pdf( + i_tally, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab); + + // resample if energy is greater than maximum neutron energy + constexpr int neutron = static_cast(ParticleType::neutron); + if (E_out < data::energy_max[neutron]) + break; + + // check for large number of resamples + ++n_sample; + if (n_sample == MAX_SAMPLE) { + // particle_write_restart(p) + fatal_error("Resampled energy distribution maximum number of times " + "for nuclide " + + nuc->name_); + } + } + // starting scoring loop on ghost particles + + for (size_t index = 0; index < ghost_particles.size(); ++index) { + auto& ghost_p = ghost_particles[index]; + double pdf_lab = pdfs_lab[index]; + score_ghost_particle(ghost_p, pdf_lab, i_tally); + + } // for loop on ghost particles +} + void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) { + p.v_t() = {0, 0, 0}; // copy energy of neutron double E_in = p.E(); diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 3ba2c0cfb05..f38753a64c2 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -13,7 +13,9 @@ #include "openmc/secondary_correlated.h" #include "openmc/secondary_kalbach.h" #include "openmc/secondary_nbody.h" +#include "openmc/secondary_thermal.h" #include "openmc/secondary_uncorrelated.h" +#include "openmc/tallies/tally_scoring.h" namespace openmc { @@ -128,5 +130,71 @@ void ReactionProduct::sample( distribution_[0]->sample(E_in, E_out, mu, seed); } } +void ReactionProduct::get_pdf(int i_tally, double E_in, double& E_out, + uint64_t* seed, Particle& p, std::vector& mu_cm, + std::vector& Js, std::vector& ghost_particles, + std::vector& pdfs_lab) const +{ + + double det_pos[4]; + get_det_pos(det_pos, i_tally); + + int distribution_index; + auto n = applicability_.size(); + if (n > 1) { + double prob = 0.0; + double c = prn(seed); + for (int i = 0; i < n; ++i) { + // Determine probability that i-th energy distribution is sampled + prob += applicability_[i](E_in); + + // If i-th distribution is sampled, sample energy from the distribution + if (c <= prob) { + // distribution_[i]->sample(E_in, E_out, mu, seed); + distribution_index = i; + break; + } + } + } else { + // If only one distribution is present, go ahead and sample it + // distribution_[0]->sample(E_in, E_out, mu, seed); + distribution_index = 0; + } + // now extract pdf + + AngleEnergy* angleEnergyPtr = distribution_[distribution_index].get(); + + if (CorrelatedAngleEnergy* correlatedAE = + dynamic_cast(angleEnergyPtr)) { + + (*correlatedAE) + .get_pdf( + det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab); + // Handle CorrelatedAngleEnergy + } else if (KalbachMann* kalbachMann = + dynamic_cast(angleEnergyPtr)) { + + (*kalbachMann) + .get_pdf( + det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab); + + // Handle KalbachMann + } else if (NBodyPhaseSpace* nBodyPS = + dynamic_cast(angleEnergyPtr)) { + + (*nBodyPS).get_pdf( + det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab); + // Handle NBodyPhaseSpace + } else if (UncorrelatedAngleEnergy* uncorrelatedAE = + dynamic_cast(angleEnergyPtr)) { + + (*uncorrelatedAE) + .get_pdf( + det_pos, E_in, E_out, seed, p, mu_cm, Js, ghost_particles, pdfs_lab); + // Handle UncorrelatedAngleEnergy + } else { + std::cout << "Unknown derived type." << std::endl; + } +} } // namespace openmc diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 0e4891dd3e7..84c2fc2eda6 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -10,8 +10,11 @@ #include "openmc/endf.h" #include "openmc/hdf5_interface.h" +#include "openmc/nuclide.h" +#include "openmc/particle.h" #include "openmc/random_lcg.h" #include "openmc/search.h" +#include "openmc/tallies/tally_scoring.h" namespace openmc { @@ -263,4 +266,152 @@ void CorrelatedAngleEnergy::sample( } } +void CorrelatedAngleEnergy::get_pdf(double det_pos[4], double E_in, + double& E_out, uint64_t* seed, Particle& p, std::vector& mu_cm, + std::vector& Js, std::vector& ghost_particles, + std::vector& pdfs_lab) const +{ + // Find energy bin and calculate interpolation factor -- if the energy is + // outside the range of the tabulated energies, choose the first or last bins + auto n_energy_in = energy_.size(); + int i; + double r; + if (E_in < energy_[0]) { + i = 0; + r = 0.0; + } else if (E_in > energy_[n_energy_in - 1]) { + i = n_energy_in - 2; + r = 1.0; + } else { + i = lower_bound_index(energy_.begin(), energy_.end(), E_in); + r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]); + } + + // Sample between the ith and [i+1]th bin + int l = r > prn(seed) ? i + 1 : i; + + // Interpolation for energy E1 and EK + int n_energy_out = distribution_[i].e_out.size(); + int n_discrete = distribution_[i].n_discrete; + double E_i_1 = distribution_[i].e_out[n_discrete]; + double E_i_K = distribution_[i].e_out[n_energy_out - 1]; + + n_energy_out = distribution_[i + 1].e_out.size(); + n_discrete = distribution_[i + 1].n_discrete; + double E_i1_1 = distribution_[i + 1].e_out[n_discrete]; + double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1]; + + double E_1 = E_i_1 + r * (E_i1_1 - E_i_1); + double E_K = E_i_K + r * (E_i1_K - E_i_K); + + // Determine outgoing energy bin + n_energy_out = distribution_[l].e_out.size(); + n_discrete = distribution_[l].n_discrete; + double r1 = prn(seed); + double c_k = distribution_[l].c[0]; + int k = 0; + int end = n_energy_out - 2; + + // Discrete portion + for (int j = 0; j < n_discrete; ++j) { + k = j; + c_k = distribution_[l].c[k]; + if (r1 < c_k) { + end = j; + break; + } + } + + // Continuous portion + double c_k1; + for (int j = n_discrete; j < end; ++j) { + k = j; + c_k1 = distribution_[l].c[k + 1]; + if (r1 < c_k1) + break; + k = j + 1; + c_k = c_k1; + } + + double E_l_k = distribution_[l].e_out[k]; + double p_l_k = distribution_[l].p[k]; + if (distribution_[l].interpolation == Interpolation::histogram) { + // Histogram interpolation + if (p_l_k > 0.0 && k >= n_discrete) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = E_l_k; + } + + } else if (distribution_[l].interpolation == Interpolation::lin_lin) { + // Linear-linear interpolation + double E_l_k1 = distribution_[l].e_out[k + 1]; + double p_l_k1 = distribution_[l].p[k + 1]; + + double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); + if (frac == 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = + E_l_k + + (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - + p_l_k) / + frac; + } + } + + // Now interpolate between incident energy bins i and i + 1 + if (k >= n_discrete) { + if (l == i) { + E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); + } else { + E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); + } + } + + const auto& nuc {data::nuclides[p.event_nuclide()]}; + const auto& rx {nuc->reactions_[p.event_index_mt()]}; + if (rx->scatter_in_cm_) { + get_pdf_to_point_elastic( + det_pos, p, mu_cm, Js, ghost_particles, E_out / 1e6); + for (std::size_t i = 0; i < mu_cm.size(); ++i) { + // Assuming Js.size() is the same as mu_cm.size() + double mu_c = mu_cm[i]; + double derivative = Js[i]; + double pdf_cm; + if (r1 - c_k < c_k1 - r1 || + distribution_[l].interpolation == Interpolation::histogram) { + pdf_cm = distribution_[l].angle[k]->get_pdf(mu_c); + } else { + pdf_cm = distribution_[l].angle[k + 1]->get_pdf(mu_c); + } + pdfs_lab.push_back(pdf_cm / std::abs(derivative)); + } + } + + if (!rx->scatter_in_cm_) { + // fatal_error("Didnt implemt lab"); + Direction u_lab {det_pos[0] - p.r().x, // towards the detector + det_pos[1] - p.r().y, det_pos[2] - p.r().z}; + Direction u_lab_unit = u_lab / u_lab.norm(); // normalize + double E_lab = E_out; + Particle ghost_particle = Particle(); + ghost_particle.initialize_ghost_particle(p, u_lab_unit, E_lab); + ghost_particles.push_back(ghost_particle); + double pdf_mu_lab; + if (r1 - c_k < c_k1 - r1 || + distribution_[l].interpolation == Interpolation::histogram) { + pdf_mu_lab = + distribution_[l].angle[k]->get_pdf(u_lab_unit.dot(p.u_last())); + } else { + pdf_mu_lab = + distribution_[l].angle[k + 1]->get_pdf(u_lab_unit.dot(p.u_last())); + } + + pdfs_lab.push_back(pdf_mu_lab); + + // fatal_error("didn't implement lab"); + } +} + } // namespace openmc diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 4a7878343c6..d2c901bf510 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -9,9 +9,12 @@ #include "xtensor/xview.hpp" #include "openmc/hdf5_interface.h" +#include "openmc/nuclide.h" +#include "openmc/particle.h" #include "openmc/random_dist.h" #include "openmc/random_lcg.h" #include "openmc/search.h" +#include "openmc/tallies/tally_scoring.h" #include "openmc/vector.h" namespace openmc { @@ -237,5 +240,143 @@ void KalbachMann::sample( mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a; } } +void KalbachMann::get_pdf(double det_pos[4], double E_in, double& E_out, + uint64_t* seed, Particle& p, std::vector& mu_cm, + std::vector& Js, std::vector& ghost_particles, + std::vector& pdfs_lab) const +{ + // Find energy bin and calculate interpolation factor -- if the energy is + // outside the range of the tabulated energies, choose the first or last bins + const auto& nuc {data::nuclides[p.event_nuclide()]}; + + // double E_out; + auto n_energy_in = energy_.size(); + int i; + double r; + if (E_in < energy_[0]) { + i = 0; + r = 0.0; + } else if (E_in > energy_[n_energy_in - 1]) { + i = n_energy_in - 2; + r = 1.0; + } else { + i = lower_bound_index(energy_.begin(), energy_.end(), E_in); + r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]); + } + + // Sample between the ith and [i+1]th bin + int l = r > prn(seed) ? i + 1 : i; + + // Interpolation for energy E1 and EK + int n_energy_out = distribution_[i].e_out.size(); + int n_discrete = distribution_[i].n_discrete; + double E_i_1 = distribution_[i].e_out[n_discrete]; + double E_i_K = distribution_[i].e_out[n_energy_out - 1]; + + n_energy_out = distribution_[i + 1].e_out.size(); + n_discrete = distribution_[i + 1].n_discrete; + double E_i1_1 = distribution_[i + 1].e_out[n_discrete]; + double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1]; + + double E_1 = E_i_1 + r * (E_i1_1 - E_i_1); + double E_K = E_i_K + r * (E_i1_K - E_i_K); + + // Determine outgoing energy bin + n_energy_out = distribution_[l].e_out.size(); + n_discrete = distribution_[l].n_discrete; + double r1 = prn(seed); + double c_k = distribution_[l].c[0]; + int k = 0; + int end = n_energy_out - 2; + + // Discrete portion + for (int j = 0; j < n_discrete; ++j) { + k = j; + c_k = distribution_[l].c[k]; + if (r1 < c_k) { + end = j; + break; + } + } + + // Continuous portion + double c_k1; + for (int j = n_discrete; j < end; ++j) { + k = j; + c_k1 = distribution_[l].c[k + 1]; + if (r1 < c_k1) + break; + k = j + 1; + c_k = c_k1; + } + + 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) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = E_l_k; + } + // Determine Kalbach-Mann parameters + km_r = distribution_[l].r[k]; + km_a = distribution_[l].a[k]; + + } else { + // Linear-linear interpolation + double E_l_k1 = distribution_[l].e_out[k + 1]; + double p_l_k1 = distribution_[l].p[k + 1]; + + double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); + if (frac == 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = + E_l_k + + (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - + p_l_k) / + frac; + } + // Linear-linear interpolation + // Determine Kalbach-Mann parameters + km_r = distribution_[l].r[k] + + (E_out - E_l_k) / (E_l_k1 - E_l_k) * + (distribution_[l].r[k + 1] - distribution_[l].r[k]); + km_a = distribution_[l].a[k] + + (E_out - E_l_k) / (E_l_k1 - E_l_k) * + (distribution_[l].a[k + 1] - distribution_[l].a[k]); + } + + // Now interpolate between incident energy bins i and i + 1 + if (k >= n_discrete) { + if (l == i) { + E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); + } else { + E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); + } + } + + get_pdf_to_point_elastic(det_pos, p, mu_cm, Js, ghost_particles, E_out / 1e6); + for (std::size_t i = 0; i < mu_cm.size(); ++i) { + // Assuming Js.size() is the same as mu_cm.size() + double mu_c = mu_cm[i]; + double derivative = Js[i]; + double pdf_cm = km_a / (2 * std::sinh(km_a)) * + (std::cosh(km_a * mu_c) + + km_r * std::sinh(km_a * mu_c)); // center of mass + pdfs_lab.push_back(pdf_cm / std::abs(derivative)); + } + + // https://docs.openmc.org/en/v0.8.0/methods/physics.html#equation-KM-pdf-angle + // double pdf_mu = km_a / (2 * std::sinh(km_a)) * (std::cosh(km_a * mymu) + + // km_r * std::sinh(km_a * mymu)); // center of mass return pdf_mu; + + const auto& rx {nuc->reactions_[p.event_index_mt()]}; + if (!rx->scatter_in_cm_) { + fatal_error("didn't implement lab"); + } +} } // namespace openmc diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp index da0bb81c471..cce2f153bdb 100644 --- a/src/secondary_nbody.cpp +++ b/src/secondary_nbody.cpp @@ -5,8 +5,10 @@ #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/math_functions.h" +#include "openmc/nuclide.h" #include "openmc/random_dist.h" #include "openmc/random_lcg.h" +#include "openmc/tallies/tally_scoring.h" namespace openmc { @@ -67,4 +69,65 @@ void NBodyPhaseSpace::sample( E_out = E_max * v; } +void NBodyPhaseSpace::get_pdf(double det_pos[4], double E_in, double& E_out, + uint64_t* seed, Particle& p, std::vector& mu_cm, + std::vector& Js, std::vector& ghost_particles, + std::vector& pdfs_lab) const +{ + // By definition, the distribution of the angle is isotropic for an N-body + // phase space distribution + + // Determine E_max parameter + double Ap = mass_ratio_; + double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_); + + // x is essentially a Maxwellian distribution + double x = maxwell_spectrum(1.0, seed); + + double y; + double r1, r2, r3, r4, r5, r6; + switch (n_bodies_) { + case 3: + y = maxwell_spectrum(1.0, seed); + break; + case 4: + r1 = prn(seed); + r2 = prn(seed); + r3 = prn(seed); + y = -std::log(r1 * r2 * r3); + break; + case 5: + r1 = prn(seed); + r2 = prn(seed); + r3 = prn(seed); + r4 = prn(seed); + r5 = prn(seed); + r6 = prn(seed); + y = -std::log(r1 * r2 * r3 * r4) - + 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."}; + } + + // Now determine v and E_out + double v = x / (x + y); + E_out = E_max * v; + + get_pdf_to_point_elastic(det_pos, p, mu_cm, Js, ghost_particles, E_out / 1e6); + for (std::size_t i = 0; i < mu_cm.size(); ++i) { + // Assuming Js.size() is the same as mu_cm.size() + double mu_c = mu_cm[i]; + double derivative = Js[i]; + double pdf_cm = 0.5; + pdfs_lab.push_back(pdf_cm / std::abs(derivative)); + } + + const auto& nuc {data::nuclides[p.event_nuclide()]}; + const auto& rx {nuc->reactions_[p.event_index_mt()]}; + if (!rx->scatter_in_cm_) { + fatal_error("didn't implement lab"); + } +} + } // namespace openmc diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 8b9e8737c66..6c790d5a611 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -25,6 +25,299 @@ void get_energy_index( } } +// Structure to hold information about nearest neighbors and their counts +struct Neighbors { + double a_0, a_1, b_0, b_1; + int a_0_index, b_0_index; + int leftCount, rightCount; +}; + +Neighbors findNearestNeighbors( + const std::vector& sortedVector, double mu_0) +{ + Neighbors result; + result.leftCount = 0; + result.rightCount = 0; + result.a_0_index = -1; // Initialize with invalid index + result.b_0_index = -1; // Initialize with invalid index + + // Find the iterator pointing to the first element that is not less than mu_0 + auto it = std::lower_bound(sortedVector.begin(), sortedVector.end(), mu_0); + + // Process neighbors to the left + if (it != sortedVector.begin()) { + --it; + result.a_0 = *it; + result.a_0_index = std::distance(sortedVector.begin(), it); + ++result.leftCount; + + // Process the second left neighbor, if available + if (it != sortedVector.begin()) { + --it; + result.a_1 = *it; + ++result.leftCount; + } else { + result.a_1 = + result + .a_0; // If no second left neighbor, set it to the same as the first + } + } else { + result.a_0 = result.a_1 = mu_0; // If no left neighbors, set both to mu_0 + } + + // Find the iterator pointing to the first element that is greater than mu_0 + it = std::upper_bound(sortedVector.begin(), sortedVector.end(), mu_0); + + // Process neighbors to the right + if (it != sortedVector.end()) { + result.b_0 = *it; + result.b_0_index = std::distance(sortedVector.begin(), it); + ++result.rightCount; + + // Process the second right neighbor, if available + if (++it != sortedVector.end()) { + result.b_1 = *it; + ++result.rightCount; + } else { + result.b_1 = + result + .b_0; // If no second right neighbor, set it to the same as the first + } + } else { + result.b_0 = result.b_1 = mu_0; // If no right neighbors, set both to mu_0 + } + + // Return the result + return result; +} + +int checkNeighborCase(Neighbors neighbors) +{ + + if (neighbors.leftCount == 2 && neighbors.rightCount == 2) + return 1; + if (neighbors.leftCount == 0 && neighbors.rightCount == 2) + return 2; + if (neighbors.leftCount == 1 && neighbors.rightCount == 2) + return 3; + if (neighbors.leftCount == 2 && neighbors.rightCount == 1) + return 4; + if (neighbors.leftCount == 2 && neighbors.rightCount == 0) + return 5; + if (neighbors.leftCount == 1 && neighbors.rightCount == 1) + return 6; + if (neighbors.leftCount == 0 && neighbors.rightCount == 1) + return 7; + if (neighbors.leftCount == 1 && neighbors.rightCount == 0) + return 8; + + return -1; +} + +double get_pdf_discrete(const std::vector& sortedVector, double mu_0, + double nu, const std::vector& cdfVector = std::vector {-1}) +{ + // Use findNearestNeighbors to get information about nearest neighbors + // Make sure mu is in range [-1,1] and return + if (std::abs(mu_0) > 1.0) + mu_0 = std::copysign(1.0, mu_0); + Neighbors neighbors = findNearestNeighbors(sortedVector, mu_0); + int mycase = checkNeighborCase(neighbors); + double pdf; + double a_0 = neighbors.a_0; + double a_1 = neighbors.a_1; + double b_0 = neighbors.b_0; + double b_1 = neighbors.b_1; + int a_0_index = neighbors.a_0_index; + int b_0_index = neighbors.b_0_index; + int mu_index; + double prob; + switch (mycase) { + case 1: { + // Calculate Delta_a and Delta_b + double delta_a = 0.5 * std::min(b_0 - a_0, a_0 - a_1); + double delta_b = 0.5 * std::min(b_1 - b_0, b_0 - a_0); + + // Check conditions and calculate pdf + if (mu_0 >= a_0 && mu_0 < a_0 + delta_a) { + pdf = 1.0 / (2.0 * delta_a); + mu_index = a_0_index; + } else if (mu_0 > b_0 - delta_b && mu_0 <= b_0) { + pdf = 1.0 / (2.0 * delta_b); + mu_index = b_0_index; + } else if (mu_0 > a_0 + delta_a && mu_0 <= b_0 - delta_b) { + pdf = 0.0; + mu_index = 0; + } else { + pdf = -1; // + mu_index = -1; + } + break; + } + case 2: { + // Calculate Delta_b for Case 2 + double delta_b_case2 = 0.5 * std::min(b_1 - b_0, b_0 - (-1)); + // Check conditions and calculate pdf + if (mu_0 > b_0 - delta_b_case2 && mu_0 <= b_0) { + pdf = 1.0 / (2.0 * delta_b_case2); + mu_index = b_0_index; + } else if (mu_0 >= -1 && mu_0 <= b_0 - delta_b_case2) { + pdf = 0.0; + mu_index = 0; + } else { + pdf = -1; // + mu_index = -1; + } + break; + } + case 3: { + + // Calculate Delta_a and Delta_b for Case 3 + double delta_a_case3 = 0.5 * std::min(b_0 - a_0, a_0 - (-1)); + double delta_b_case3 = 0.5 * std::min(b_1 - b_0, b_0 - a_0); + + // Check conditions and calculate pdf + if (mu_0 >= a_0 && mu_0 < a_0 + delta_a_case3) { + pdf = 1.0 / (2.0 * delta_a_case3); + mu_index = a_0_index; + } else if (mu_0 > b_0 - delta_b_case3 && mu_0 <= b_0) { + pdf = 1.0 / (2.0 * delta_b_case3); + mu_index = b_0_index; + } else if (mu_0 > a_0 + delta_a_case3 && mu_0 <= b_0 - delta_b_case3) { + pdf = 0.0; + mu_index = 0; + } else { + pdf = -1; // + mu_index = -1; + } + break; + } + case 4: { + // Calculate Delta_a and Delta_b for Case 4 + double delta_a_case4 = 0.5 * std::min(b_0 - a_0, a_0 - a_1); + double delta_b_case4 = 0.5 * std::min(1 - b_0, b_0 - a_0); + + // Check conditions and calculate pdf + if (mu_0 >= a_0 && mu_0 < a_0 + delta_a_case4) { + pdf = 1.0 / (2.0 * delta_a_case4); + mu_index = a_0_index; + } else if (mu_0 > b_0 - delta_b_case4 && mu_0 <= b_0) { + pdf = 1.0 / (2.0 * delta_b_case4); + mu_index = b_0_index; + } else if (mu_0 > a_0 + delta_a_case4 && mu_0 <= b_0 - delta_b_case4) { + pdf = 0.0; + mu_index = 0; + } else { + pdf = -1; // + mu_index = -1; + } + + break; + } + case 5: { + // Calculate Delta_a for Case 5 + double delta_a_case5 = 0.5 * std::min(1 - a_0, a_0 - a_1); + + // Check conditions and calculate pdf + if (mu_0 >= a_0 && mu_0 < a_0 + delta_a_case5) { + pdf = 1.0 / (2.0 * delta_a_case5); + mu_index = a_0_index; + } else if (mu_0 > a_0 + delta_a_case5 && mu_0 <= 1) { + pdf = 0.0; + mu_index = 0; + } else { + pdf = -1; // + mu_index = -1; + } + + break; + } + case 6: { // mu has exactly one left neighbor and exactly one right neighbor + + // Calculate Delta_a and Delta_b for Case 6 + double delta_a_case6 = 0.5 * std::min(b_0 - a_0, a_0 - (-1)); + double delta_b_case6 = 0.5 * std::min(1 - b_0, b_0 - a_0); + + // Check conditions and calculate pdf + if (mu_0 >= a_0 && mu_0 < a_0 + delta_a_case6) { + pdf = 1.0 / (2.0 * delta_a_case6); + mu_index = a_0_index; + } else if (mu_0 > b_0 - delta_b_case6 && mu_0 <= b_0) { + pdf = 1.0 / (2.0 * delta_b_case6); + mu_index = b_0_index; + } else if (mu_0 > a_0 + delta_a_case6 && mu_0 <= b_0 - delta_b_case6) { + pdf = 0.0; + mu_index = 0; + } else { + pdf = -1; // + mu_index = -1; + } + + break; + } + case 7: { // mu has no left neighbor and exactly one right neighbor + // Calculate Delta_a and Delta_b for Case 7 + double delta_case7 = 0.5 * std::min(1 - b_0, b_0 - (-1)); + + // Check conditions and calculate pdf + if (mu_0 >= b_0 - delta_case7 && mu_0 <= b_0) { + pdf = 1.0 / (2.0 * delta_case7); + mu_index = b_0_index; + } else if (mu_0 >= (-1) && mu_0 <= b_0 - delta_case7) { + pdf = 0.0; + mu_index = 0; + } else { + pdf = -1; // + mu_index = b_0_index; + } + + break; + } + + case 8: { + // Calculate Delta_a and Delta_b for Case 8 + double delta_case8 = 0.5 * std::min(1 - a_0, a_0 - (-1)); + + // Check conditions and calculate pdf + if (mu_0 >= a_0 && mu_0 <= a_0 + delta_case8) { + pdf = 1.0 / (2.0 * delta_case8); + mu_index = a_0_index; + } else if (mu_0 <= 1 && mu_0 >= a_0 + delta_case8) { + pdf = 0.0; + mu_index = 0; + } else { + + pdf = -1; // + mu_index = -1; + } + + break; + } + default: { + + pdf = -1; + } + } + if (cdfVector[0] == -1) + // equiprobable distribution + prob = 1 / nu; + else { + // calculate probability for each mu + std::vector pVector; + if (!cdfVector.empty()) { + double norm = cdfVector.back(); + + for (size_t i = cdfVector.size() - 1; i > 0; --i) { + pVector.push_back((cdfVector[i] - cdfVector[i - 1]) / norm); + } + pVector.push_back(cdfVector[0] / norm); + } + prob = pVector[mu_index]; + } + + return pdf * prob; +} + //============================================================================== // CoherentElasticAE implementation //============================================================================== @@ -55,6 +348,36 @@ void CoherentElasticAE::sample( mu = 1.0 - 2.0 * energies[k] / E_in; } +double CoherentElasticAE::get_pdf( + double E_in, double& E_out, double& mu, 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 i = lower_bound_index(energies.begin(), energies.end(), E_in); + std::vector energies_cut(energies.begin(), energies.begin() + i + 1); + std::vector factors_cut(factors.begin(), factors.begin() + i + 1); + + std::vector mu_vector_rev; + std::transform(energies_cut.begin(), energies_cut.end(), + std::back_inserter(mu_vector_rev), + [E_in](double Ei) { return 1 - 2 * Ei / E_in; }); + std::vector mu_vector(mu_vector_rev.rbegin(), mu_vector_rev.rend()); + // find mu index in mu_vector + const int k_i = lower_bound_index(mu_vector.begin(), mu_vector.end(), mu); + const int k_f = k_i + 1; + + return get_pdf_discrete(mu_vector, mu, 1, factors_cut); +} + //============================================================================== // IncoherentElasticAE implementation //============================================================================== @@ -74,6 +397,19 @@ void IncoherentElasticAE::sample( // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4) E_out = E_in; } +double IncoherentElasticAE::get_pdf( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4 + double c = 2 * E_in * debye_waller_; + E_out = E_in; + + double A = c / (1 - std::exp(-2.0 * c)); // normalization factor + double pdf = A * std::exp(-c * (1 - mu)); + return pdf; + + // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4) +} //============================================================================== // IncoherentElasticAEDiscrete implementation @@ -129,6 +465,27 @@ void IncoherentElasticAEDiscrete::sample( E_out = E_in; } +double IncoherentElasticAEDiscrete::get_pdf( + double E_in, double& E_out, double& mu, 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; + + 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, n_mu); +} + //============================================================================== // IncoherentInelasticAEDiscrete implementation //============================================================================== @@ -201,6 +558,59 @@ void IncoherentInelasticAEDiscrete::sample( // Cosine of angle between incoming and outgoing neutron mu = (1 - f) * mu_ijk + f * mu_i1jk; } +double IncoherentInelasticAEDiscrete::get_pdf( + double E_in, double& E_out, double& mu, uint64_t* seed, int l) const +{ + // Get index and interpolation factor for inelastic grid + int i; + double f; + get_energy_index(energy_, E_in, i, f); + int j; + int n = energy_out_.shape()[1]; + if (!skewed_) { + // All bins equally likely + j = prn(seed) * n; + } else { + // Distribution skewed away from edge points + double r = prn(seed) * (n - 3); + if (r > 1.0) { + // equally likely N-4 middle bins + j = r + 1; + } else if (r > 0.6) { + // second to last bin has relative probability of 0.4 + j = n - 2; + } else if (r > 0.5) { + // last bin has relative probability of 0.1 + j = n - 1; + } else if (r > 0.1) { + // second bin has relative probability of 0.4 + j = 1; + } else { + // first bin has relative probability of 0.1 + j = 0; + } + } + if (l != -1) { + j = l; + } // take j as param + // Determine outgoing energy corresponding to E_in[i] and E_in[i+1] + double E_ij = energy_out_(i, j); + double E_i1j = energy_out_(i + 1, j); + + // Outgoing energy + E_out = (1 - f) * E_ij + f * E_i1j; + int m = mu_out_.shape()[2]; + std::vector mu_vector; + + 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, m); +} //============================================================================== // IncoherentInelasticAE implementation @@ -331,6 +741,76 @@ void IncoherentInelasticAE::sample( mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5); } +double IncoherentInelasticAE::get_pdf( + 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); + + // Pick closer energy based on interpolation factor + int l = f > 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; + 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) + break; + c_j = c_j1; + } + + // check to make sure j is <= n_energy_out - 2 + j = std::min(j, n - 2); + + // Get the data to interpolate between + double E_l_j = distribution_[l].e_out[j]; + double p_l_j = distribution_[l].e_out_pdf[j]; + + // Next part assumes linear-linear interpolation in standard + double E_l_j1 = distribution_[l].e_out[j + 1]; + double p_l_j1 = distribution_[l].e_out_pdf[j + 1]; + + // Find secondary energy (variable E) + double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j); + if (frac == 0.0) { + E_out = E_l_j + (r1 - c_j) / p_l_j; + } else { + E_out = E_l_j + + (std::sqrt(std::max(0.0, p_l_j * p_l_j + 2.0 * frac * (r1 - c_j))) - + p_l_j) / + frac; + } + + // Adjustment of outgoing energy + double E_l = energy_[l]; + if (E_out < 0.5 * E_l) { + E_out *= 2.0 * E_in / E_l - 1.0; + } else { + E_out += E_in - E_l; + } + + // Sample outgoing cosine bin + int n_mu = distribution_[l].mu.shape()[1]; + const auto& mu_l = distribution_[l].mu; + f = (r1 - c_j) / (c_j1 - c_j); + + 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, n_mu); +} + //============================================================================== // MixedElasticAE implementation //============================================================================== diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp index 5cbb76fb9b9..50a26defed7 100644 --- a/src/secondary_uncorrelated.cpp +++ b/src/secondary_uncorrelated.cpp @@ -6,7 +6,10 @@ #include "openmc/error.h" #include "openmc/hdf5_interface.h" +#include "openmc/nuclide.h" +#include "openmc/particle.h" #include "openmc/random_dist.h" +#include "openmc/tallies/tally_scoring.h" namespace openmc { @@ -65,4 +68,60 @@ void UncorrelatedAngleEnergy::sample( E_out = energy_->sample(E_in, seed); } +void UncorrelatedAngleEnergy::get_pdf(double det_pos[4], double E_in, + double& E_out, uint64_t* seed, Particle& p, std::vector& mu_cm, + std::vector& Js, std::vector& ghost_particles, + std::vector& pdfs_lab) const +{ + bool COM = false; + const auto& nuc {data::nuclides[p.event_nuclide()]}; + if (p.event_index_mt() != -999) { + const auto& rx {nuc->reactions_[p.event_index_mt()]}; + COM = rx->scatter_in_cm_; + } + + // Sample outgoing energy + E_out = energy_->sample(E_in, seed); + + if (COM) { + + get_pdf_to_point_elastic( + det_pos, p, mu_cm, Js, ghost_particles, E_out / 1e6); + for (std::size_t i = 0; i < mu_cm.size(); ++i) { + // Assuming Js.size() is the same as mu_cm.size() + double mu_c = mu_cm[i]; + double derivative = Js[i]; + double pdf_cm; + if (!angle_.empty()) { + pdf_cm = angle_.get_pdf(E_in, mu_c, seed); + } else { + // no angle distribution given => assume isotropic for all energies + pdf_cm = 0.5; + } + pdfs_lab.push_back(pdf_cm / std::abs(derivative)); + } + } + + if (!COM) { + // finding mu_cm, E_out is in lab + // fatal_error("Didnt implemt lab"); + Direction u_lab {det_pos[0] - p.r().x, // towards the detector + det_pos[1] - p.r().y, det_pos[2] - p.r().z}; + Direction u_lab_unit = u_lab / u_lab.norm(); // normalize + double E_lab = E_out; + Particle ghost_particle = Particle(); + ghost_particle.initialize_ghost_particle(p, u_lab_unit, E_lab); + ghost_particles.push_back(ghost_particle); + double pdf_mu_lab; + if (!angle_.empty()) { + pdf_mu_lab = angle_.get_pdf(E_in, u_lab_unit.dot(p.u_last()), seed); + } else { + // no angle distribution given => assume isotropic for all energies + pdf_mu_lab = 0.5; + } + + pdfs_lab.push_back(pdf_mu_lab); + } +} + } // namespace openmc diff --git a/src/source.cpp b/src/source.cpp index 12323f7bd7e..22a04753019 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -649,7 +649,8 @@ SourceSite sample_external_source(uint64_t* seed) data::mg.rev_energy_bins_.end(), site.E); site.E = data::mg.num_energy_groups_ - site.E - 1.; } - + site.ext = true; + site.source_index = i; return site; } diff --git a/src/state_point.cpp b/src/state_point.cpp index 8195c486507..93082ee32b9 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -207,6 +207,8 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) write_dataset(tally_group, "estimator", "tracklength"); } else if (tally->estimator_ == TallyEstimator::COLLISION) { write_dataset(tally_group, "estimator", "collision"); + } else if (tally->estimator_ == TallyEstimator::POINT) { + write_dataset(tally_group, "estimator", "point"); } write_dataset(tally_group, "n_realizations", tally->n_realizations_); diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index e737f79d354..dd0523644d7 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -57,6 +57,7 @@ vector active_tallies; vector active_analog_tallies; vector active_tracklength_tallies; vector active_collision_tallies; +vector active_point_tallies; vector active_meshsurf_tallies; vector active_surface_tallies; vector active_pulse_height_tallies; @@ -171,6 +172,10 @@ Tally::Tally(pugi::xml_node node) this->set_nuclides(node); + // READ DATA FOR POSITIONS + + this->set_positions(node); + // ======================================================================= // READ DATA FOR SCORES @@ -342,6 +347,7 @@ Tally::Tally(pugi::xml_node node) // Check if user specified estimator if (check_for_node(node, "estimator")) { std::string est = get_node_value(node, "estimator"); + fmt::print("est = {}\n", est); if (est == "analog") { estimator_ = TallyEstimator::ANALOG; } else if (est == "tracklength" || est == "track-length" || @@ -370,6 +376,18 @@ Tally::Tally(pugi::xml_node node) // Set estimator to collision estimator estimator_ = TallyEstimator::COLLISION; + } else if (est == "point") { + // If the estimator was set to an analog estimator, this means the + // tally needs post-collision information + if (estimator_ == TallyEstimator::ANALOG) { + throw std::runtime_error {fmt::format("Cannot use collision estimator " + "for tally ", + id_)}; + } + + // Set estimator to collision estimator + estimator_ = TallyEstimator::POINT; + } else { throw std::runtime_error { fmt::format("Invalid estimator '{}' on tally {}", est, id_)}; @@ -720,6 +738,31 @@ void Tally::set_nuclides(const vector& nuclides) } } +void Tally::set_positions(pugi::xml_node node) +{ + positions_.clear(); + + // By default, we tally just the total material rates. + if (!check_for_node(node, "positions")) { + positions_.push_back(""); + return; + } + + // The user provided specifics nuclides. Parse it as an array with either + // "total" or a nuclide name like "U235" in each position. + auto words = get_node_array(node, "positions"); + this->set_positions(words); +} + +void Tally::set_positions(const vector& positions) +{ + positions_.clear(); + + for (const auto& pos : positions) { + positions_.push_back(pos); + } +} + void Tally::init_triggers(pugi::xml_node node) { for (auto trigger_node : node.children("trigger")) { @@ -1074,6 +1117,7 @@ void setup_active_tallies() model::active_analog_tallies.clear(); model::active_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(); @@ -1093,6 +1137,9 @@ void setup_active_tallies() case TallyEstimator::TRACKLENGTH: model::active_tracklength_tallies.push_back(i); break; + case TallyEstimator::POINT: + model::active_point_tallies.push_back(i); + break; case TallyEstimator::COLLISION: model::active_collision_tallies.push_back(i); } @@ -1128,6 +1175,7 @@ void free_memory_tally() model::active_analog_tallies.clear(); model::active_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(); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e73fb90f311..55819244c80 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1,9 +1,10 @@ #include "openmc/tallies/tally_scoring.h" - #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/constants.h" +#include "openmc/distribution_multi.h" #include "openmc/error.h" +#include "openmc/geometry.h" #include "openmc/ifp.h" #include "openmc/material.h" #include "openmc/mgxs_interface.h" @@ -11,17 +12,24 @@ #include "openmc/photon.h" #include "openmc/reaction_product.h" #include "openmc/search.h" +#include "openmc/secondary_correlated.h" +#include "openmc/secondary_uncorrelated.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" #include "openmc/tallies/filter_cell.h" #include "openmc/tallies/filter_delayedgroup.h" #include "openmc/tallies/filter_energy.h" - +#include +#include +#include #include - +#include +int ghost_counter = 0; +int col_counter = 0; namespace openmc { //============================================================================== @@ -1614,6 +1622,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, double filter_weight, int i_nuclide, double atom_density, double flux) { + auto& tally {*model::tallies[i_tally]}; // Set the direction and group to use with get_xs @@ -2559,6 +2568,147 @@ void score_collision_tally(Particle& p) match.bins_present_ = false; } +void score_point_tally(Particle& p) +{ + if ((p.event_mt() == 101) || (p.event_mt() == 18) || + (p.event_index_mt() == -1234)) + return; // absorption or fission or s(a,b) + + for (auto i_tally : model::active_point_tallies) { + + double yield = 1; + col_counter++; + const auto& nuc {data::nuclides[p.event_nuclide()]}; + const auto& rx {nuc->reactions_[p.event_index_mt()]}; + auto& d = rx->products_[0].distribution_[0]; + auto d_ = dynamic_cast(d.get()); + // Initialize + std::vector mu_cm; + std::vector Js; + std::vector ghost_particles; + std::vector pdfs_lab; + // std::vector fluxes; + if (std::isnan(p.r().x)) { + return; // sometimes the get_pdf function turns p.r() into nan + } + + // Get position (x,y,z) of detector + double det_pos[4]; + get_det_pos(det_pos, i_tally); + + if (p.event_mt() == 2 && p.event_index_mt() != -1234) { + get_pdf_to_point_elastic(det_pos, p, mu_cm, Js, ghost_particles); + for (std::size_t i = 0; i < mu_cm.size(); ++i) { + // Assuming Js.size() is the same as mu_cm.size() + double mu_c = mu_cm[i]; + double derivative = Js[i]; + double pdf_cm = d_->angle().get_pdf(p.E_last(), mu_c, p.current_seed()); + pdfs_lab.push_back(pdf_cm / std::abs(derivative)); + } + } + if (p.event_mt() != 2) { // Inelastic + // make sure v_t is 0 + // copy energy of neutron + double E_in = p.E_last(); + + double E_out; + rx->products_[0].get_pdf(i_tally, E_in, E_out, p.current_seed(), p, mu_cm, + Js, ghost_particles, pdfs_lab); + + yield = (*rx->products_[0].yield_)(p.E_last()); + if (std::floor(yield) != yield && yield > 0) { + yield = 1; + } + for (auto& p_ghost : ghost_particles) { + p_ghost.wgt() = p_ghost.wgt() * yield; + ; + } + } + + for (size_t index = 0; index < ghost_particles.size(); ++index) { + auto& ghost_p = ghost_particles[index]; + double pdf_lab = pdfs_lab[index]; + score_ghost_particle(ghost_p, pdf_lab, i_tally); + // calculate shielding + } // for loop on ghost particles + } +} + +void get_det_pos(double (&det_pos)[4], int i_tally) +{ + const Tally& tally {*model::tallies[i_tally]}; + if (tally.positions_.size() == 4) { + for (auto i = 0; i < 3; ++i) { + auto pos_coord = tally.positions_[i]; + det_pos[i] = std::stod(pos_coord); + } + auto R0 = tally.positions_[3]; + det_pos[4] = std::stod(R0); + } else { + fatal_error("user must use 3 positions and R0"); + } +} + +void score_point_tally_from_source(const SourceSite* src) +{ + if (!src->ext) { + return; + } + double flux = 0; + for (auto i_tally : model::active_point_tallies) { + double det_pos[4]; // Get position (x,y,z) of detector + get_det_pos(det_pos, i_tally); + Direction u_lab {det_pos[0] - src->r.x, // towards the detector + det_pos[1] - src->r.y, det_pos[2] - src->r.z}; + Direction u_lab_unit = u_lab / u_lab.norm(); // normalize + // Get the angle distribution + double pdf_mu_det = 0; + Source& mysource = *(model::external_sources[src->source_index]); + if (auto* independentSource = dynamic_cast(&mysource)) + // Check if the actual type is IndependentSource + { + UnitSphereDistribution* angleDistribution = independentSource->angle(); + if (angleDistribution) { + Direction my_u_ref = angleDistribution->u_ref_; + Direction my_u_ref_unit = my_u_ref / my_u_ref.norm(); + if (typeid(*angleDistribution) == typeid(PolarAzimuthal)) + // polar and assuming phi is isotropic + { + double my_det_mu = my_u_ref_unit.dot(u_lab_unit); + if (std::abs(my_det_mu) > 1.0) { + my_det_mu = std::copysign(1.0, my_det_mu); + } + auto* polarAzimuthalDistribution = + dynamic_cast(angleDistribution); + Distribution* muDistribution = polarAzimuthalDistribution->mu(); + Distribution* phiDistribution = polarAzimuthalDistribution->phi(); + pdf_mu_det = muDistribution->get_pdf(my_det_mu); + } else if (typeid(*angleDistribution) == typeid(Isotropic)) + // Isotropic + { + pdf_mu_det = 0.5; + } else if (typeid(*angleDistribution) == typeid(Monodirectional)) + // pencil + { + if (my_u_ref_unit.dot(u_lab_unit) == 1) { + double total_distance = u_lab.norm(); + pdf_mu_det = + (2 * PI * total_distance * total_distance); // to cancel the Rs + } // pdf is not defined for discrete case + } else // unexpected type + { + fatal_error("unexpected type"); + } + } + + Particle ghost_particle = Particle(); + ghost_particle.initialize_ghost_particle_from_source(src, u_lab_unit); + + score_ghost_particle(ghost_particle, pdf_mu_det, i_tally); + } + } +} + void score_surface_tally(Particle& p, const vector& tallies) { double current = p.wgt_last(); @@ -2603,6 +2753,326 @@ void score_surface_tally(Particle& p, const vector& tallies) match.bins_present_ = false; } +void boostf(double A[4], double B[4], double X[4]) +{ + // + // boosts B(labfram) to A rest frame and gives output in X + // + double W; + int j; + + if ((A[0] * A[0] - A[1] * A[1] - A[2] * A[2] - A[3] * A[3]) <= 0) { + } + + W = sqrt(A[0] * A[0] - A[1] * A[1] - A[2] * A[2] - A[3] * A[3]); + + if (W == 0 || W == (-A[0])) + + X[0] = (B[0] * A[0] - B[3] * A[3] - B[2] * A[2] - B[1] * A[1]) / W; + for (j = 1; j <= 3; j++) { + X[j] = B[j] - A[j] * (B[0] + X[0]) / (A[0] + W); + } + + return; +} + +double get_MFP(Particle& ghost_particle, double total_distance) +{ + // calculate shilding + double remaining_distance = total_distance; + double total_MFP = 0; + ghost_particle.event_calculate_xs(); + ghost_particle.boundary() = distance_to_boundary(ghost_particle); + double advance_distance = ghost_particle.boundary().distance(); + + while (advance_distance < remaining_distance) // Advance to next boundary + { + total_MFP += advance_distance * ghost_particle.macro_xs().total; + // Advance particle in space and time + for (int j = 0; j < ghost_particle.n_coord(); ++j) { + ghost_particle.coord(j).r() += + advance_distance * ghost_particle.coord(j).u(); + } + remaining_distance -= advance_distance; + ghost_particle.time() += advance_distance / ghost_particle.speed(); + ghost_particle.event_cross_surface(); + ghost_particle.event_calculate_xs(); + ghost_particle.boundary() = distance_to_boundary(ghost_particle); + advance_distance = ghost_particle.boundary().distance(); + } + total_MFP += remaining_distance * + ghost_particle.macro_xs().total; // advance to next boundary + for (int j = 0; j < ghost_particle.n_coord(); ++j) { + ghost_particle.coord(j).r() += + remaining_distance * ghost_particle.coord(j).u(); + } + ghost_particle.time() += remaining_distance / ghost_particle.speed(); + return total_MFP; +} + +void get_pdf_to_point_elastic(double det_pos[4], Particle& p, + std::vector& mu_cm, std::vector& Js, + std::vector& ghost_particles, double E3k_cm_given) +{ + Direction u_lab {det_pos[0] - p.r().x, // towards the detector + det_pos[1] - p.r().y, det_pos[2] - p.r().z}; + Direction u_lab_unit = u_lab / u_lab.norm(); // normalize + + double m1 = p.getMass() / 1e6; // mass of incoming particle in MeV + const auto& nuc {data::nuclides[p.event_nuclide()]}; + double awr = nuc->awr_; + double m2 = m1 * awr; // mass of target + double m3 = m1; // mass of outgoing particle to detector + double m4 = m2; // mass of recoil target system + + double E1_tot = + p.E_last() / 1e6 + m1; // total Energy of incoming particle in MeV + double p1_tot = std::sqrt( + E1_tot * E1_tot - m1 * m1); // total momenta of incoming particle in MeV + // without this the get_pdf function turns p.r() into nan + Direction p1 = p1_tot * p.u_last(); // 3 momentum of incoming particle + Direction p2 = p.v_t() * m2 / C_LIGHT; // 3 momentum of target in lab + double E2_tot = std::sqrt(p2.norm() * p2.norm() + m2 * m2); + double E_cm = E1_tot + E2_tot; + Direction p_cm = p1 + p2; + double p_tot_cm = p_cm.norm(); + + double cos_lab = u_lab_unit.dot(p_cm) / (p_tot_cm); // between cm and p3 + if (std::abs(cos_lab) > 1.0) { + cos_lab = std::copysign(1.0, cos_lab); + } + + double theta = std::acos(cos_lab); + double sin_lab_sq = 1 - cos_lab * cos_lab; + + double M_cm = std::sqrt( + E_cm * E_cm - + p_tot_cm * p_tot_cm); // mass of the center of mass (incoming and target) + double gamma = E_cm / M_cm; + double p1_cm[4]; + double A[4] = {E_cm, p_cm.x, p_cm.y, p_cm.z}; + // double invA[4] = {E_cm, -p_cm.x , -p_cm.y , -p_cm.z}; + // boostf( invA ,p1_cm, maybe_p1_lab); boost back to lab + double B[4] = {E1_tot, p1.x, p1.y, p1.z}; + boostf(A, B, p1_cm); + double p1_tot_cm = + std::sqrt(p1_cm[1] * p1_cm[1] + p1_cm[2] * p1_cm[2] + p1_cm[3] * p1_cm[3]); + double E3_cm = (M_cm * M_cm + m3 * m3 - m4 * m4) / (2 * M_cm); + if (E3k_cm_given >= 0.0) { + E3_cm = E3k_cm_given + m3; + m4 = std::sqrt(M_cm * M_cm + m3 * m3 - 2 * M_cm * E3_cm); + } + double p3_tot_cm = std::sqrt(E3_cm * E3_cm - m3 * m3); + double cond = (M_cm / p_tot_cm) * (p3_tot_cm / m3); + double insq = (M_cm * M_cm * p3_tot_cm * p3_tot_cm - + m3 * m3 * p_tot_cm * p_tot_cm * sin_lab_sq); + double p3_tot_1 = 0; + double p3_tot_2 = 0; + double E3k_1 = 0; + double E3k_2 = 0; + Direction p3_1 = {0, 0, 0}; + Direction p3_2 = {0, 0, 0}; + double Fp3cm_1[4]; + double Fp3cm_2[4]; + const auto& rx {nuc->reactions_[0]}; + // auto& d = rx->products_[0].distribution_[0]; + // auto d_ = dynamic_cast(d.get()); + + double q = (p_tot_cm / E_cm) * (E3_cm / p3_tot_cm); + double approx_tol = 0.0001; + + if (insq >= 0) //( (cond > 1) || ( (cond < 1) && (theta < std::asin(cond)) ) ) + { + // first solution + + p3_tot_1 = ((M_cm * M_cm + m3 * m3 - m4 * m4) * p_tot_cm * cos_lab + + 2 * E_cm * std::sqrt(insq)) / + 2 / (M_cm * M_cm + p_tot_cm * p_tot_cm * sin_lab_sq); + if (p3_tot_1 <= 0) + return; + p3_1 = u_lab_unit * p3_tot_1; + double E3_tot_1 = std::sqrt(p3_tot_1 * p3_tot_1 + m3 * m3); + E3k_1 = (E3_tot_1 - m3) * 1e6; // back to eV + double B1[4] = {E3_tot_1, p3_1.x, p3_1.y, p3_1.z}; + boostf(A, B1, Fp3cm_1); + + double p3cm_tot_1 = + std::sqrt(Fp3cm_1[1] * Fp3cm_1[1] + Fp3cm_1[2] * Fp3cm_1[2] + + Fp3cm_1[3] * Fp3cm_1[3]); + double mucm_1 = + (Fp3cm_1[1] * p1_cm[1] + Fp3cm_1[2] * p1_cm[2] + Fp3cm_1[3] * p1_cm[3]) / + (p1_tot_cm * p3cm_tot_1); // good until here + if (std::abs(mucm_1) > 1.0) { + mucm_1 = std::copysign(1.0, mucm_1); + } + // double pdf1cm = d_->angle().get_pdf(p.E_last(),mucm_1,p.current_seed()); + // pdfs_cm.push_back(pdf1cm); + mu_cm.push_back(mucm_1); + + double mucm03_1 = + (Fp3cm_1[1] * p_cm.x + Fp3cm_1[2] * p_cm.y + Fp3cm_1[3] * p_cm.z) / + (p_tot_cm * p3cm_tot_1); + + if (std::abs(mucm03_1) > 1.0) { + mucm03_1 = std::copysign(1.0, mucm03_1); + } + double sincm1 = std::sqrt( + 1 - mucm03_1 * mucm03_1); // if this is zero derivative is inf so pdf is 0 + double sin_ratio1 = std::sqrt(sin_lab_sq) / sincm1; + double derivative1 = + gamma * (1 + q * mucm03_1) * (sin_ratio1 * sin_ratio1 * sin_ratio1); + if (sin_lab_sq < approx_tol) { + derivative1 = ((cos_lab) / (gamma * mucm03_1 * (1 + q * mucm03_1))) * + ((cos_lab) / (gamma * mucm03_1 * (1 + q * mucm03_1))); + } + Js.push_back(derivative1); + + Particle ghost_particle = Particle(); + ghost_particle.initialize_ghost_particle(p, u_lab_unit, E3k_1); + ghost_particles.push_back(ghost_particle); + + if (true) //((cond < 1) && (theta < std::asin(cond))) + { + // second solution + + p3_tot_2 = ((M_cm * M_cm + m3 * m3 - m4 * m4) * p_tot_cm * cos_lab - + 2 * E_cm * std::sqrt(insq)) / + 2 / (M_cm * M_cm + p_tot_cm * p_tot_cm * sin_lab_sq); + if (p3_tot_2 < 0) + return; + p3_2 = u_lab_unit * p3_tot_2; + double E3_tot_2 = std::sqrt(p3_tot_2 * p3_tot_2 + m3 * m3); + E3k_2 = (E3_tot_2 - m3) * 1e6; + if (p3_tot_2 < 0 || E3k_2 < 0) + return; + double B2[4] = {E3_tot_2, p3_2.x, p3_2.y, p3_2.z}; + boostf(A, B2, Fp3cm_2); + double p3cm_tot_2 = + std::sqrt(Fp3cm_2[1] * Fp3cm_2[1] + Fp3cm_2[2] * Fp3cm_2[2] + + Fp3cm_2[3] * Fp3cm_2[3]); + double mucm_2 = (Fp3cm_2[1] * p1_cm[1] + Fp3cm_2[2] * p1_cm[2] + + Fp3cm_2[3] * p1_cm[3]) / + (p1_tot_cm * p3cm_tot_2); + if (std::abs(mucm_2) > 1) { + mucm_2 = std::copysign(1.0, mucm_2); + } + // double pdf2cm = + // d_->angle().get_pdf(p.E_last(),mucm_2,p.current_seed()); + // pdfs_cm.push_back(pdf2cm); + mu_cm.push_back(mucm_2); + + double mucm03_2 = + (Fp3cm_2[1] * p_cm.x + Fp3cm_2[2] * p_cm.y + Fp3cm_2[3] * p_cm.z) / + (p_tot_cm * p3cm_tot_1); + if (std::abs(mucm03_2) > 1.0) { + mucm03_2 = std::copysign(1.0, mucm03_2); + } + double sincm2 = std::sqrt(1 - mucm03_2 * mucm03_2); + double sin_ratio2 = std::sqrt(sin_lab_sq) / sincm2; + double derivative2 = + gamma * (1 + q * mucm03_2) * (sin_ratio2 * sin_ratio2 * sin_ratio2); + if (sin_lab_sq < approx_tol) { + derivative2 = ((cos_lab) / (gamma * mucm03_2 * (1 + q * mucm03_2))) * + ((cos_lab) / (gamma * mucm03_2 * (1 + q * mucm03_2))); + } + Js.push_back(derivative2); + + Particle ghost_particle = Particle(); + ghost_particle.initialize_ghost_particle(p, u_lab_unit, E3k_2); + ghost_particles.push_back(ghost_particle); + } + } +} +void score_ghost_particle(Particle& ghost_p, double pdf_lab, int i_tally) +{ + + double myflux; + double det_pos[4]; + get_det_pos(det_pos, i_tally); + Direction u_lab {det_pos[0] - ghost_p.r().x, // towards the detector + det_pos[1] - ghost_p.r().y, det_pos[2] - ghost_p.r().z}; + Direction u_lab_unit = u_lab / u_lab.norm(); // normalize + double total_distance = u_lab.norm(); + double R0 = det_pos[3]; + double total_MFP1 = get_MFP(ghost_p, total_distance); + + if (total_distance < R0) { + if (ghost_p.macro_xs().total == 0) { + myflux = (ghost_p.wgt() * pdf_lab) / (2 / 3 * PI * R0 * R0); + + } else { + // mutliplying in 1=10000/10000 to avoid float point error + myflux = (10000 * ghost_p.wgt() * pdf_lab * + (1 - exp(-ghost_p.macro_xs().total * R0))) / + (10000 * 2 / 3 * PI * R0 * R0 * R0 * ghost_p.macro_xs().total); + } + } else { + myflux = (ghost_p.wgt()) * exp(-total_MFP1) / + (2 * PI * total_distance * total_distance) * pdf_lab; + } + + if (std::isnan(myflux) || std::abs(myflux) > 1e10) { + + // Cause a segmentation fault to crash the program + int* ptr = nullptr; + *ptr = 42; // This will trigger a segmentation fault + } + + if (myflux < 0) { + fatal_error("negetive flux"); + } + if (ghost_p.type() != ParticleType::neutron) { + myflux = 0; + } + 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, ghost_p); + auto end = FilterBinIter(tally, true, &ghost_p.filter_matches()); + if (filter_iter == end) + return; + + // 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]; + + double atom_density = 0.; + if (i_nuclide >= 0) { + auto j = + model::materials[ghost_p.material()]->mat_nuclide_index_[i_nuclide]; + if (j == C_NONE) + continue; + atom_density = model::materials[ghost_p.material()]->atom_density_(j); + } + // TODO: consider replacing this "if" with pointers or templates + if (settings::run_CE) { + score_general_ce_nonanalog(ghost_p, i_tally, i * tally.scores_.size(), + filter_index, filter_weight, i_nuclide, atom_density, myflux); + } else { + fatal_error("multi group not implemnted for point tally"); + } + } + } + + // 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) + return; + + // Reset all the filter matches for the next tally event. + for (auto& match : ghost_p.filter_matches()) + match.bins_present_ = false; +} + void score_pulse_height_tally(Particle& p, const vector& tallies) { // The pulse height tally in OpenMC hijacks the logic of CellFilter and @@ -2678,4 +3148,5 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) p.E_last() = orig_E_last; } } +// } // namespace openmc diff --git a/src/thermal.cpp b/src/thermal.cpp index cbe0983ed65..a4206f31e41 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -1,15 +1,15 @@ #include "openmc/thermal.h" - -#include // for sort, move, min, max, find -#include // for round, sqrt, abs - #include "xtensor/xarray.hpp" #include "xtensor/xbuilder.hpp" #include "xtensor/xmath.hpp" #include "xtensor/xsort.hpp" #include "xtensor/xtensor.hpp" #include "xtensor/xview.hpp" +#include // for sort, move, min, max, find +#include // for round, sqrt, abs #include +#include +#include #include "openmc/constants.h" #include "openmc/endf.h" @@ -313,6 +313,52 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E, *mu = std::copysign(1.0, *mu); } +double ThermalData::get_pdf(const NuclideMicroXS& micro_xs, double E, + double& E_out, double mu, uint64_t* seed) +{ + double pdf = -1; + AngleEnergy* angleEnergyPtr; + // Determine whether inelastic or elastic scattering occured + if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) { + // elastic_.distribution->get_pdf(E, *E_out, *mu, seed); + angleEnergyPtr = elastic_.distribution.get(); + } else { + // inelastic_.distribution->get_pdf(E, *E_out, *mu, seed); + angleEnergyPtr = inelastic_.distribution.get(); + } + + if (CoherentElasticAE* coherentElasticAE = + dynamic_cast(angleEnergyPtr)) { + pdf = (*coherentElasticAE).get_pdf(E, E_out, mu, seed); + + bool creat_pdf_file = false; + + // Handle CoherentElasticAE + } else if (IncoherentElasticAE* incoherentElasticAE = + dynamic_cast(angleEnergyPtr)) { + // Handle IncoherentElasticAE + } else if (IncoherentElasticAEDiscrete* incoherentElasticAEDiscrete = + dynamic_cast(angleEnergyPtr)) { + // Handle IncoherentElasticAEDiscrete + } else if (IncoherentInelasticAEDiscrete* incoherentInelasticAEDiscrete = + dynamic_cast(angleEnergyPtr)) { + pdf = (*incoherentInelasticAEDiscrete).get_pdf(E, E_out, mu, seed, -1); + bool creat_pdf_file = false; + // Handle IncoherentInelasticAEDiscrete + } else if (IncoherentInelasticAE* incoherentInelasticAE = + dynamic_cast(angleEnergyPtr)) { + + // Handle IncoherentInelasticAE + } else if (MixedElasticAE* mixedElasticAE = + dynamic_cast(angleEnergyPtr)) { + // Handle MixedElasticAE + } else { + std::cout << "Unknown derived type." << std::endl; + } + + return pdf; +} + void free_memory_thermal() { data::thermal_scatt.clear(); diff --git a/vendor/gsl-lite b/vendor/gsl-lite new file mode 160000 index 00000000000..913e86d49c6 --- /dev/null +++ b/vendor/gsl-lite @@ -0,0 +1 @@ +Subproject commit 913e86d49c6a1acca980f4e325378f9dc393493a