crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Convert sample_photon_product to C++

Browse source
This commit is contained in:
Paul Romano 2018-11-20 07:35:30 -06:00
parent 85b60badc2
commit 9105cd195b
3 changed files with 31 additions and 97 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/physics.h
View file

@ -47,7 +47,7 @@ void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
Reaction* sample_fission(int i_nuclide, double E);
extern "C" void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
void absorption(Particle* p, int i_nuclide);

59
src/physics.F90
View file

@ -254,65 +254,6 @@ contains
end function sample_element
!===============================================================================
! SAMPLE_PHOTON_PRODUCT
!===============================================================================
subroutine sample_photon_product(i_nuclide, E, i_reaction, i_product) bind(C)
integer(C_INT), value :: i_nuclide ! index in nuclides array
real(C_DOUBLE), value :: E ! energy of neutron
integer(C_INT), intent(out) :: i_reaction ! index in nuc % reactions array
integer(C_INT), intent(out) :: i_product ! index in reaction % products array
integer :: i_grid
integer :: i_temp
integer :: threshold
integer :: last_valid_reaction
integer :: last_valid_product
real(8) :: f
real(8) :: prob
real(8) :: cutoff
real(8) :: yield
! Get pointer to nuclide
associate (nuc => nuclides(i_nuclide))
! Get grid index and interpolation factor and sample photon production cdf
i_temp = micro_xs(i_nuclide) % index_temp
i_grid = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
cutoff = prn() * micro_xs(i_nuclide) % photon_prod
prob = ZERO
! Loop through each reaction type
REACTION_LOOP: do i_reaction = 1, size(nuc % reactions)
associate (rx => nuc % reactions(i_reaction))
threshold = rx % xs_threshold(i_temp)
! if energy is below threshold for this reaction, skip it
if (i_grid < threshold) cycle
do i_product = 1, rx % products_size()
if (rx % product_particle(i_product) == PHOTON) then
! add to cumulative probability
yield = rx % product_yield(i_product, E)
prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) &
+ f*(rx % xs(i_temp, i_grid - threshold + 2))) * yield
if (prob > cutoff) return
last_valid_reaction = i_reaction
last_valid_product = i_product
end if
end do
end associate
end do REACTION_LOOP
end associate
i_reaction = last_valid_reaction
i_product = last_valid_product
end subroutine sample_photon_product
!===============================================================================
! SCATTER
!===============================================================================

67
src/physics.cpp
View file

@ -530,45 +530,38 @@ Reaction* sample_fission(int i_nuclide, double E)
}
}
// void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
// {
// // Get pointer to nuclide
// associate (nuc => nuclides(i_nuclide))
void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product)
{
// Get grid index and interpolation factor and sample photon production cdf
int i_temp = simulation::micro_xs[i_nuclide-1].index_temp;
int i_grid = simulation::micro_xs[i_nuclide-1].index_grid;
double f = simulation::micro_xs[i_nuclide-1].interp_factor;
double cutoff = prn() * simulation::micro_xs[i_nuclide-1].photon_prod;
double prob = 0.0;
// // Get grid index and interpolation factor and sample photon production cdf
// i_temp = simulation::micro_xs[i_nuclide-1].index_temp
// i_grid = simulation::micro_xs[i_nuclide-1].index_grid
// f = simulation::micro_xs[i_nuclide-1].interp_factor
// cutoff = prn() * simulation::micro_xs[i_nuclide-1].photon_prod
// prob = 0.0
// Loop through each reaction type
const auto& nuc {data::nuclides[i_nuclide-1]};
for (int i = 0; i < nuc->reactions_.size(); ++i) {
const auto& rx = nuc->reactions_[i];
int threshold = rx->xs_[i_temp-1].threshold;
// // Loop through each reaction type
// REACTION_LOOP: do i_reaction = 1, size(nuc % reactions)
// associate (rx => nuc % reactions(i_reaction))
// threshold = rx % xs_threshold(i_temp)
// if energy is below threshold for this reaction, skip it
if (i_grid < threshold) continue;
// // if energy is below threshold for this reaction, skip it
// if (i_grid < threshold) cycle
for (int j = 0; j < rx->products_.size(); ++j) {
if (rx->products_[j].particle_ == ParticleType::photon) {
// add to cumulative probability
double yield = (*rx->products_[j].yield_)(E);
prob += ((1.0 - f) * rx->xs_[i_temp-1].value[i_grid - threshold]
+ f*(rx->xs_[i_temp-1].value[i_grid - threshold + 1])) * yield;
// do i_product = 1, rx % products_size()
// if (rx % product_particle(i_product) == PHOTON) {
// // add to cumulative probability
// yield = rx % product_yield(i_product, E)
// prob = prob + ((1.0 - f) * rx % xs(i_temp, i_grid - threshold + 1) &
// + f*(rx % xs(i_temp, i_grid - threshold + 2))) * yield
// if (prob > cutoff) return
// last_valid_reaction = i_reaction
// last_valid_product = i_product
// }
// end do
// end associate
// end do REACTION_LOOP
// end associate
// i_reaction = last_valid_reaction
// i_product = last_valid_product
// }
*i_rx = i;
*i_product = j;
if (prob > cutoff) return;
}
}
}
}
void absorption(Particle* p, int i_nuclide)
{
@ -1144,10 +1137,10 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
sample_photon_product(i_nuclide, p->E, &i_rx, &i_product);
// Sample the outgoing energy and angle
auto& rx = data::nuclides[i_nuclide-1]->reactions_[i_rx-1];
auto& rx = data::nuclides[i_nuclide-1]->reactions_[i_rx];
double E;
double mu;
rx->products_[i_product-1].sample(p->E, E, mu);
rx->products_[i_product].sample(p->E, E, mu);
// Sample the new direction
double uvw[3];