From 934e26ca533a2d77e431a3db6b8b83ef76297af4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 17 Nov 2018 16:08:36 -0600 Subject: [PATCH] Move sample_nuclide to C++ --- include/openmc/physics.h | 2 +- src/physics.F90 | 62 --------------------------- src/physics.cpp | 92 +++++++++++++++++++++++----------------- 3 files changed, 53 insertions(+), 103 deletions(-) diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 96d6d10e9e..6ea043c36e 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -36,7 +36,7 @@ void sample_electron_reaction(Particle* p); //! MeV) are created and travel in opposite directions. void sample_positron_reaction(Particle* p); -extern "C" void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat); +void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat); //! Determine the average total, prompt, and delayed neutrons produced from //! fission and creates appropriate bank sites. diff --git a/src/physics.F90 b/src/physics.F90 index fdef8ecd7f..7f956be907 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -211,68 +211,6 @@ contains end subroutine sample_photon_reaction -!=============================================================================== -! SAMPLE_NUCLIDE -!=============================================================================== - - subroutine sample_nuclide(p, base, i_nuclide, i_nuc_mat) bind(C) - - type(Particle), intent(in) :: p - integer(C_INT), value :: base ! which reaction to sample based on - integer(C_INT), intent(out) :: i_nuclide - integer(C_INT), intent(out) :: i_nuc_mat - - real(8) :: prob - real(8) :: cutoff - real(8) :: atom_density ! atom density of nuclide in atom/b-cm - real(8) :: sigma ! microscopic total xs for nuclide - type(Material), pointer :: mat - - ! Get pointer to current material - mat => materials(p % material) - - ! Sample cumulative distribution function - select case (base) - case (SCORE_TOTAL) - cutoff = prn() * material_xs % total - case (SCORE_SCATTER) - cutoff = prn() * (material_xs % total - material_xs % absorption) - case (SCORE_FISSION) - cutoff = prn() * material_xs % fission - end select - - i_nuc_mat = 0 - prob = ZERO - do while (prob < cutoff) - i_nuc_mat = i_nuc_mat + 1 - - ! Check to make sure that a nuclide was sampled - if (i_nuc_mat > mat % n_nuclides) then - call particle_write_restart(p) - call fatal_error("Did not sample any nuclide during collision.") - end if - - ! Find atom density - i_nuclide = mat % nuclide(i_nuc_mat) - atom_density = mat % atom_density(i_nuc_mat) - - ! Determine microscopic cross section - select case (base) - case (SCORE_TOTAL) - sigma = atom_density * micro_xs(i_nuclide) % total - case (SCORE_SCATTER) - sigma = atom_density * (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) - case (SCORE_FISSION) - sigma = atom_density * micro_xs(i_nuclide) % fission - end select - - ! Increment probability to compare to cutoff - prob = prob + sigma - end do - - end subroutine sample_nuclide - !=============================================================================== ! SAMPLE_ELEMENT !=============================================================================== diff --git a/src/physics.cpp b/src/physics.cpp index db160c3a25..44b141b0dc 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -397,57 +397,69 @@ void sample_positron_reaction(Particle* p) p->alive = false; } -// void sample_nuclide(Particle* p, int mt, int i_nuclide, int i_nuc_mat) -// { -// // Get pointer to current material -// mat => materials(p->material) +void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat) +{ + // Sample cumulative distribution function + double cutoff; + switch (mt) { + case SCORE_TOTAL: + cutoff = prn() * simulation::material_xs.total; + break; + case SCORE_SCATTER: + cutoff = prn() * (simulation::material_xs.total - + simulation::material_xs.absorption); + break; + case SCORE_FISSION: + cutoff = prn() * simulation::material_xs.fission; + break; + } -// // Sample cumulative distribution function -// select case (base) -// case ('total') -// cutoff = prn() * material_xs % total -// case ('scatter') -// cutoff = prn() * (material_xs % total - material_xs % absorption) -// case ('fission') -// cutoff = prn() * material_xs % fission -// end select + // Get pointers to nuclide/density arrays + int* nuclides; + double* densities; + int n; + openmc_material_get_densities(p->material, &nuclides, &densities, &n); -// i_nuc_mat = 0 -// prob = 0.0 -// do while (prob < cutoff) -// i_nuc_mat = i_nuc_mat + 1 + *i_nuc_mat = 0; + double prob = 0.0; + while (prob < cutoff) { + // Check to make sure that a nuclide was sampled + if (*i_nuc_mat > n) { + p->write_restart(); + fatal_error("Did not sample any nuclide during collision."); + } -// // Check to make sure that a nuclide was sampled -// if (i_nuc_mat > mat % n_nuclides) { -// particle_write_restart(p) -// fatal_error("Did not sample any nuclide during collision.") -// } + // Find atom density + *i_nuclide = nuclides[*i_nuc_mat]; + double atom_density = densities[*i_nuc_mat]; -// // Find atom density -// i_nuclide = mat % nuclide(i_nuc_mat) -// atom_density = mat % atom_density(i_nuc_mat) + // Determine microscopic cross section + double sigma; + switch (mt) { + case SCORE_TOTAL: + sigma = atom_density * simulation::micro_xs[*i_nuclide-1].total; + break; + case SCORE_SCATTER: + sigma = atom_density * (simulation::micro_xs[*i_nuclide-1].total - + simulation::micro_xs[*i_nuclide-1].absorption); + break; + case SCORE_FISSION: + sigma = atom_density * simulation::micro_xs[*i_nuclide-1].fission; + break; + } -// // Determine microscopic cross section -// select case (base) -// case ('total') -// sigma = atom_density * micro_xs(i_nuclide) % total -// case ('scatter') -// sigma = atom_density * (micro_xs(i_nuclide) % total - & -// micro_xs(i_nuclide) % absorption) -// case ('fission') -// sigma = atom_density * micro_xs(i_nuclide) % fission -// end select + // Increment probability to compare to cutoff + prob += sigma; -// // Increment probability to compare to cutoff -// prob = prob + sigma -// end do -// } + ++(*i_nuc_mat); + } +} // void sample_element(Particle* p) // { // associate (mat => materials(p->material)) // // Sample cumulative distribution function -// cutoff = prn() * material_xs % total +// cutoff = prn() * simulation::material_xs.total // i = 0 // prob = 0.0