From b803032267c71bec6ea04050528622329c94ba6d Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 6 May 2026 22:43:47 +0300 Subject: [PATCH] implement jacobian calculation and energy and angle transformation --- include/openmc/angle_energy.h | 48 +++++++++++++++-- include/openmc/chain.h | 8 +-- include/openmc/reaction_product.h | 8 +-- include/openmc/secondary_correlated.h | 8 +-- include/openmc/secondary_kalbach.h | 8 +-- include/openmc/secondary_nbody.h | 8 +-- include/openmc/secondary_thermal.h | 48 ++++++++++------- include/openmc/secondary_uncorrelated.h | 8 +-- include/openmc/thermal.h | 6 ++- src/chain.cpp | 8 ++- src/reaction_product.cpp | 7 +-- src/secondary_correlated.cpp | 11 ++-- src/secondary_kalbach.cpp | 10 ++-- src/secondary_nbody.cpp | 9 ++-- src/secondary_thermal.cpp | 68 +++++++++++++++++-------- src/secondary_uncorrelated.cpp | 12 +++-- src/thermal.cpp | 5 +- 17 files changed, 197 insertions(+), 83 deletions(-) diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h index 55deb5d415e..1858d7a3b1c 100644 --- a/include/openmc/angle_energy.h +++ b/include/openmc/angle_energy.h @@ -1,8 +1,11 @@ #ifndef OPENMC_ANGLE_ENERGY_H #define OPENMC_ANGLE_ENERGY_H +#include // for sqrt #include +#include "openmc/random_lcg.h" + namespace openmc { //============================================================================== @@ -27,12 +30,51 @@ class AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - virtual double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const = 0; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + virtual double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const = 0; virtual ~AngleEnergy() = default; }; +inline double get_jac_and_transform( + double E_in, double& mu, double& E_out, uint64_t* seed, double awr) +{ + double E_cm = E_out; + double mu_lab = mu; + double D = mu_lab * mu_lab - 1.0 + (awr + 1.0) * (awr + 1.0) * E_cm / E_in; + if (D < 0.0) + return 0.0; + D = std::sqrt(D); + double E_out1 = + E_in * ((mu_lab + D) / (awr + 1.0)) * ((mu_lab + D) / (awr + 1.0)); + double E_out2 = + E_in * ((mu_lab - D) / (awr + 1.0)) * ((mu_lab - D) / (awr + 1.0)); + double mu_cm1 = + mu_lab * std::sqrt(E_out1 / E_cm) - std::sqrt(E_in / E_cm) / (awr + 1.0); + double mu_cm2 = + mu_lab * std::sqrt(E_out1 / E_cm) - std::sqrt(E_in / E_cm) / (awr + 1.0); + double mult = 1.0; + if (std::abs(mu_cm1) > 1.0) { + mu = mu_cm2; + E_out = E_out2; + } else if (std::abs(mu_cm2) > 1.0) { + mu = mu_cm1; + E_out = E_out1; + } else { + mult = 2.0; + if (prn(seed) < 0.5) { + mu = mu_cm2; + E_out = E_out2; + } else { + mu = mu_cm2; + E_out = E_out2; + } + } + return mult * E_out * (awr + 1.0) / (D * std::sqrt(E_cm * E_in)); +} + } // namespace openmc #endif // OPENMC_ANGLE_ENERGY_H diff --git a/include/openmc/chain.h b/include/openmc/chain.h index 6f58303585a..052dad486ce 100644 --- a/include/openmc/chain.h +++ b/include/openmc/chain.h @@ -76,9 +76,11 @@ class DecayPhotonAngleEnergy : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: const Distribution* photon_energy_; diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index 79a6160e25b..ba464633650 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -60,9 +60,11 @@ class ReactionProduct { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const; ParticleType particle_; //!< Particle type EmissionMode emission_mode_; //!< Emission mode diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h index 69b22981aff..df5bd9e9c35 100644 --- a/include/openmc/secondary_correlated.h +++ b/include/openmc/secondary_correlated.h @@ -53,9 +53,11 @@ class CorrelatedAngleEnergy : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; // energy property vector& energy() { return energy_; } diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index b25352be9f0..22e55285edb 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -46,9 +46,11 @@ class KalbachMann : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: //! Outgoing energy/angle at a single incoming energy diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h index 9d033a6b869..eb1b53b9dba 100644 --- a/include/openmc/secondary_nbody.h +++ b/include/openmc/secondary_nbody.h @@ -39,9 +39,11 @@ class NBodyPhaseSpace : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: int n_bodies_; //!< Number of particles distributed diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 45d2c5260f1..54d791d8498 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -39,9 +39,11 @@ class CoherentElasticAE : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section @@ -74,9 +76,11 @@ class IncoherentElasticAE : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: double debye_waller_; @@ -108,9 +112,11 @@ class IncoherentElasticAEDiscrete : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: const vector& energy_; //!< Energies at which cosines are tabulated @@ -149,9 +155,11 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: const vector& energy_; //!< Incident energies @@ -196,9 +204,11 @@ class IncoherentInelasticAE : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: //! Secondary energy/angle distribution @@ -246,9 +256,11 @@ class MixedElasticAE : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; private: CoherentElasticAE coherent_dist_; //!< Coherent distribution diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h index f895ae77f63..94301e7e16d 100644 --- a/include/openmc/secondary_uncorrelated.h +++ b/include/openmc/secondary_uncorrelated.h @@ -37,9 +37,11 @@ class UncorrelatedAngleEnergy : public AngleEnergy { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine - double sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const override; + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine + double sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const override; // Accessors AngleDistribution& angle() { return angle_; } diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index c06a2ee0dde..a50af6cfd0b 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -67,9 +67,11 @@ class ThermalData { //! \param[in] mu Scattering cosine with respect to current direction //! \param[out] E_out Outgoing energy in [eV] //! \param[inout] seed Pseudorandom seed pointer - //! \return Probability density for the scattering cosine + //! \param[in] is_com Is scattering cosine is given in center of mass + //! coordinates \param[in] awr Weight of nucleus in neutron masses \return + //! Probability density for the scattering cosine double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in, - double mu, double& E_out, uint64_t* seed) const; + double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const; private: struct Reaction { diff --git a/src/chain.cpp b/src/chain.cpp index e4d0324d3b5..ab7be11235c 100644 --- a/src/chain.cpp +++ b/src/chain.cpp @@ -74,9 +74,13 @@ void DecayPhotonAngleEnergy::sample( mu = Uniform(-1., 1.).sample(seed).first; } -double DecayPhotonAngleEnergy::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double DecayPhotonAngleEnergy::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { + if (is_com) + fatal_error( + "DecayPhotonAngleEnergy distribution must be in lab coordinates"); + E_out = photon_energy_->sample(seed).first; return 0.5; } diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index a1c93786166..69fc959b944 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -134,10 +134,11 @@ void ReactionProduct::sample( sample_dist(E_in, seed).sample(E_in, E_out, mu, seed); } -double ReactionProduct::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double ReactionProduct::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { - return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed); + return sample_dist(E_in, seed) + .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr); } } // namespace openmc diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 32701791ade..99b49cb1970 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -260,10 +260,15 @@ void CorrelatedAngleEnergy::sample( mu = sample_dist(E_in, E_out, seed).sample(seed).first; } -double CorrelatedAngleEnergy::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double CorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { - return sample_dist(E_in, E_out, seed).evaluate(mu); + auto& dist = sample_dist(E_in, E_out, seed); + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + + return jac * dist.evaluate(mu); } } // namespace openmc diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 018ce1c8a9c..3472d9e7391 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -232,14 +232,18 @@ void KalbachMann::sample( mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a; } } -double KalbachMann::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double KalbachMann::sample_energy_and_pdf(double E_in, double mu, double& E_out, + uint64_t* seed, bool is_com, double awr) const { double km_r, km_a; sample_params(E_in, E_out, km_a, km_r, seed); + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + // https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle - return km_a / (2 * std::sinh(km_a)) * + return jac * km_a / (2 * std::sinh(km_a)) * (std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu)); } diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp index 72f0b0b92d0..9ae7a5152f4 100644 --- a/src/secondary_nbody.cpp +++ b/src/secondary_nbody.cpp @@ -72,11 +72,14 @@ void NBodyPhaseSpace::sample( E_out = sample_energy(E_in, seed); } -double NBodyPhaseSpace::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double NBodyPhaseSpace::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { E_out = sample_energy(E_in, seed); - return 0.5; + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + return jac * 0.5; } } // namespace openmc diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index b0f601809df..7634c15ed16 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -53,11 +53,16 @@ void CoherentElasticAE::sample( mu = 1.0 - 2.0 * energies[k] / E_in; } -double CoherentElasticAE::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double CoherentElasticAE::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1) E_out = E_in; + + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + const auto& factors = xs_.factors(); if (E_in < bragg_edges_.front()) @@ -68,8 +73,9 @@ double CoherentElasticAE::sample_energy_and_pdf( double E = 0.5 * (1 - mu) * E_in; double C = 0.5 * E_in / factors[i]; - return C * get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)), - factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in); + return jac * C * + get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)), + factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in); } //============================================================================== @@ -90,15 +96,20 @@ void IncoherentElasticAE::sample( double c = 2 * E_in * debye_waller_; mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0; } -double IncoherentElasticAE::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const + +double IncoherentElasticAE::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { E_out = E_in; + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4 double c = 2 * E_in * debye_waller_; double A = c / (1 - std::exp(-2.0 * c)); // normalization factor - return A * std::exp(-c * (1 - mu)); + return jac * A * std::exp(-c * (1 - mu)); } //============================================================================== @@ -155,8 +166,8 @@ void IncoherentElasticAEDiscrete::sample( E_out = E_in; } -double IncoherentElasticAEDiscrete::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double IncoherentElasticAEDiscrete::sample_energy_and_pdf(double E_in, + double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const { // Get index and interpolation factor for elastic grid int i; @@ -165,8 +176,12 @@ double IncoherentElasticAEDiscrete::sample_energy_and_pdf( // Energy doesn't change in elastic scattering E_out = E_in; - return get_pdf_discrete_interpolated( - mu_out_.slice(i, tensor::all), mu_out_.slice(i + 1, tensor::all), f, mu); + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + + return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, tensor::all), + mu_out_.slice(i + 1, tensor::all), f, mu); } //============================================================================== @@ -253,8 +268,8 @@ void IncoherentInelasticAEDiscrete::sample( mu = (1 - f) * mu_ijk + f * mu_i1jk; } -double IncoherentInelasticAEDiscrete::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(double E_in, + double mu, double& E_out, uint64_t* seed, bool is_com, double awr) const { // Get index and interpolation factor for inelastic grid int i; @@ -263,8 +278,12 @@ double IncoherentInelasticAEDiscrete::sample_energy_and_pdf( int j; sample_params(E_in, E_out, j, seed); - return get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all), - mu_out_.slice(i + 1, j, tensor::all), f, mu); + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + + return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all), + mu_out_.slice(i + 1, j, tensor::all), f, mu); } //============================================================================== @@ -403,8 +422,8 @@ void IncoherentInelasticAE::sample( mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5); } -double IncoherentInelasticAE::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double IncoherentInelasticAE::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { double f; int l, j; @@ -412,8 +431,12 @@ double IncoherentInelasticAE::sample_energy_and_pdf( const auto& mu_l = distribution_[l].mu; - return get_pdf_discrete_interpolated( - mu_l.slice(j, tensor::all), mu_l.slice(j + 1, tensor::all), f, mu); + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + + return jac * get_pdf_discrete_interpolated(mu_l.slice(j, tensor::all), + mu_l.slice(j + 1, tensor::all), f, mu); } //============================================================================== @@ -458,10 +481,11 @@ void MixedElasticAE::sample( sample_dist(E_in, seed).sample(E_in, E_out, mu, seed); } -double MixedElasticAE::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double MixedElasticAE::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { - return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed); + return sample_dist(E_in, seed) + .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr); } } // namespace openmc diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp index ec2af710258..59818832247 100644 --- a/src/secondary_uncorrelated.cpp +++ b/src/secondary_uncorrelated.cpp @@ -65,8 +65,8 @@ void UncorrelatedAngleEnergy::sample( E_out = energy_->sample(E_in, seed); } -double UncorrelatedAngleEnergy::sample_energy_and_pdf( - double E_in, double mu, double& E_out, uint64_t* seed) const +double UncorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu, + double& E_out, uint64_t* seed, bool is_com, double awr) const { // Sample outgoing energy if (energy_ != nullptr) { @@ -75,11 +75,15 @@ double UncorrelatedAngleEnergy::sample_energy_and_pdf( E_out = E_in; } + double jac = 1.0; + if (is_com) + jac = get_jac_and_transform(E_in, mu, E_out, seed, awr); + if (!angle_.empty()) { - return angle_.evaluate(E_in, mu); + return jac * angle_.evaluate(E_in, mu); } else { // no angle distribution given => assume isotropic for all energies - return 0.5; + return jac * 0.5; } } diff --git a/src/thermal.cpp b/src/thermal.cpp index edfbddf23e8..fd753165a39 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -314,10 +314,11 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E, } double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs, - double E_in, double mu, double& E_out, uint64_t* seed) const + double E_in, double mu, double& E_out, uint64_t* seed, bool is_com, + double awr) const { return sample_dist(micro_xs, E_in, seed) - .sample_energy_and_pdf(E_in, mu, E_out, seed); + .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr); } void free_memory_thermal()