crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Sample isotropic direction consistently

Browse source
This commit is contained in:
Paul Romano 2021-03-24 08:11:40 -05:00
parent f8241f567f
commit 8a07feebcb
4 changed files with 17 additions and 31 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/distribution_multi.h
View file

@ -57,6 +57,8 @@ private:
//! Uniform distribution on the unit sphere
//==============================================================================
Direction isotropic_direction(uint64_t* seed);
class Isotropic : public UnitSphereDistribution {
public:
Isotropic() { };

7
src/distribution_multi.cpp
View file

@ -69,7 +69,7 @@ Direction PolarAzimuthal::sample(uint64_t* seed) const
// Isotropic implementation
//==============================================================================
Direction Isotropic::sample(uint64_t* seed) const
Direction isotropic_direction(uint64_t* seed)
{
double phi = 2.0*PI*prn(seed);
double mu = 2.0*prn(seed) - 1.0;
@ -77,6 +77,11 @@ Direction Isotropic::sample(uint64_t* seed) const
std::sqrt(1.0 - mu*mu) * std::sin(phi)};
}
Direction Isotropic::sample(uint64_t* seed) const
{
return isotropic_direction(seed);
}
//==============================================================================
// Monodirectional implementation
//==============================================================================

15
src/photon.cpp
View file

@ -2,6 +2,7 @@
#include "openmc/bremsstrahlung.h"
#include "openmc/constants.h"
#include "openmc/distribution_multi.h"
#include "openmc/hdf5_interface.h"
#include "openmc/nuclide.h"
#include "openmc/particle.h"
@ -662,12 +663,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
{
// If no transitions, assume fluorescent photon from captured free electron
if (shell.n_transitions == 0) {
double mu = 2.0*prn(p.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
Direction u = isotropic_direction(p.current_seed());
double E = shell.binding_energy;
p.create_secondary(p.wgt_, u, E, Particle::Type::photon);
return;
@ -687,12 +683,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
int secondary = shell.transition_subshells(i_transition, 1);
// Sample angle isotropically
double mu = 2.0*prn(p.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
Direction u = isotropic_direction(p.current_seed());
// Get the transition energy
double E = shell.transition_energy(i_transition);

24
src/physics.cpp
View file

@ -3,6 +3,7 @@
#include "openmc/bank.h"
#include "openmc/bremsstrahlung.h"
#include "openmc/constants.h"
#include "openmc/distribution_multi.h"
#include "openmc/eigenvalue.h"
#include "openmc/endf.h"
#include "openmc/error.h"
@ -425,12 +426,7 @@ void sample_positron_reaction(Particle& p)
}
// Sample angle isotropically
double mu = 2.0*prn(p.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
Direction u = isotropic_direction(p.current_seed());
// Create annihilation photon pair traveling in opposite directions
p.create_secondary(p.wgt_, u, MASS_ELECTRON_EV, Particle::Type::photon);
@ -712,13 +708,7 @@ void scatter(Particle& p, int i_nuclide)
int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
if (mat->p0_[i_nuc_mat]) {
// Sample isotropic-in-lab outgoing direction
double mu = 2.0*prn(p.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
// Change direction of particle
p.u().x = mu;
p.u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
p.u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
p.u() = isotropic_direction(p.current_seed());
p.mu_ = u_old.dot(p.u());
}
}
@ -1068,11 +1058,9 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part
}
// Sample azimuthal angle uniformly in [0, 2*pi) and assign angle
// TODO: use rotate_angle instead
double phi = 2.0*PI*prn(seed);
site->u.x = mu;
site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi);
site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi);
// TODO: account for dependence on incident neutron?
Direction ref(1., 0., 0.);
site->u = rotate_angle(ref, mu, nullptr, seed);
}
void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)