Skip to content

Commit 8a07fee

Browse files
paulromanopaulromano
committed
Sample isotropic direction consistently
1 parent f8241f5 commit 8a07fee

4 files changed

Lines changed: 17 additions & 31 deletions

File tree

‎include/openmc/distribution_multi.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -57,6 +57,8 @@ class PolarAzimuthal : public UnitSphereDistribution {
5757
//! Uniform distribution on the unit sphere
5858
//==============================================================================
5959

60+
Direction isotropic_direction(uint64_t* seed);
61+
6062
class Isotropic : public UnitSphereDistribution {
6163
public:
6264
Isotropic() { };

‎src/distribution_multi.cpp‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -69,14 +69,19 @@ Direction PolarAzimuthal::sample(uint64_t* seed) const
6969
// Isotropic implementation
7070
//==============================================================================
7171

72-
Direction Isotropic::sample(uint64_t* seed) const
72+
Direction isotropic_direction(uint64_t* seed)
7373
{
7474
double phi = 2.0*PI*prn(seed);
7575
double mu = 2.0*prn(seed) - 1.0;
7676
return {mu, std::sqrt(1.0 - mu*mu) * std::cos(phi),
7777
std::sqrt(1.0 - mu*mu) * std::sin(phi)};
7878
}
7979

80+
Direction Isotropic::sample(uint64_t* seed) const
81+
{
82+
return isotropic_direction(seed);
83+
}
84+
8085
//==============================================================================
8186
// Monodirectional implementation
8287
//==============================================================================

‎src/photon.cpp‎

Lines changed: 3 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
#include "openmc/bremsstrahlung.h"
44
#include "openmc/constants.h"
5+
#include "openmc/distribution_multi.h"
56
#include "openmc/hdf5_interface.h"
67
#include "openmc/nuclide.h"
78
#include "openmc/particle.h"
@@ -662,12 +663,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
662663
{
663664
// If no transitions, assume fluorescent photon from captured free electron
664665
if (shell.n_transitions == 0) {
665-
double mu = 2.0*prn(p.current_seed()) - 1.0;
666-
double phi = 2.0*PI*prn(p.current_seed());
667-
Direction u;
668-
u.x = mu;
669-
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
670-
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
666+
Direction u = isotropic_direction(p.current_seed());
671667
double E = shell.binding_energy;
672668
p.create_secondary(p.wgt_, u, E, Particle::Type::photon);
673669
return;
@@ -687,12 +683,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
687683
int secondary = shell.transition_subshells(i_transition, 1);
688684

689685
// Sample angle isotropically
690-
double mu = 2.0*prn(p.current_seed()) - 1.0;
691-
double phi = 2.0*PI*prn(p.current_seed());
692-
Direction u;
693-
u.x = mu;
694-
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
695-
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
686+
Direction u = isotropic_direction(p.current_seed());
696687

697688
// Get the transition energy
698689
double E = shell.transition_energy(i_transition);

‎src/physics.cpp‎

Lines changed: 6 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "openmc/bank.h"
44
#include "openmc/bremsstrahlung.h"
55
#include "openmc/constants.h"
6+
#include "openmc/distribution_multi.h"
67
#include "openmc/eigenvalue.h"
78
#include "openmc/endf.h"
89
#include "openmc/error.h"
@@ -425,12 +426,7 @@ void sample_positron_reaction(Particle& p)
425426
}
426427

427428
// Sample angle isotropically
428-
double mu = 2.0*prn(p.current_seed()) - 1.0;
429-
double phi = 2.0*PI*prn(p.current_seed());
430-
Direction u;
431-
u.x = mu;
432-
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
433-
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
429+
Direction u = isotropic_direction(p.current_seed());
434430

435431
// Create annihilation photon pair traveling in opposite directions
436432
p.create_secondary(p.wgt_, u, MASS_ELECTRON_EV, Particle::Type::photon);
@@ -712,13 +708,7 @@ void scatter(Particle& p, int i_nuclide)
712708
int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
713709
if (mat->p0_[i_nuc_mat]) {
714710
// Sample isotropic-in-lab outgoing direction
715-
double mu = 2.0*prn(p.current_seed()) - 1.0;
716-
double phi = 2.0*PI*prn(p.current_seed());
717-
718-
// Change direction of particle
719-
p.u().x = mu;
720-
p.u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
721-
p.u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
711+
p.u() = isotropic_direction(p.current_seed());
722712
p.mu_ = u_old.dot(p.u());
723713
}
724714
}
@@ -1068,11 +1058,9 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part
10681058
}
10691059

10701060
// Sample azimuthal angle uniformly in [0, 2*pi) and assign angle
1071-
// TODO: use rotate_angle instead
1072-
double phi = 2.0*PI*prn(seed);
1073-
site->u.x = mu;
1074-
site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi);
1075-
site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi);
1061+
// TODO: account for dependence on incident neutron?
1062+
Direction ref(1., 0., 0.);
1063+
site->u = rotate_angle(ref, mu, nullptr, seed);
10761064
}
10771065

10781066
void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)

0 commit comments

Comments
 (0)