From cd493057bf679a01ac601ad19374ad39aba96f99 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Sun, 7 Oct 2018 07:18:29 -0400 Subject: [PATCH 1/7] Converted physics_common to C++ --- CMakeLists.txt | 1 + include/openmc/physics_common.h | 16 ++++++++++++++++ src/physics_common.F90 | 31 +++++++++---------------------- src/physics_common.cpp | 26 ++++++++++++++++++++++++++ 4 files changed, 52 insertions(+), 22 deletions(-) create mode 100644 include/openmc/physics_common.h create mode 100644 src/physics_common.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 9facdf7ba5..e62c413476 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -423,6 +423,7 @@ add_library(libopenmc SHARED src/nuclide.cpp src/output.cpp src/particle.cpp + src/physics_common.cpp src/plot.cpp src/position.cpp src/pugixml/pugixml_c.cpp diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h new file mode 100644 index 0000000000..33fb8ed614 --- /dev/null +++ b/include/openmc/physics_common.h @@ -0,0 +1,16 @@ +//! \file physics_common.h +//! A collection of physics methods common to MG, CE, photon, etc. + +#ifndef OPENMC_PHYSICS_COMMON_H +#define OPENMC_PHYSICS_COMMON_H + +#include "openmc/particle.h" + +namespace openmc { + +//! \brief Performs the russian roulette operation for a particle +extern "C" void +russian_roulette(Particle* p); + +} // namespace openmc +#endif // OPENMC_PHYSICS_COMMON_H diff --git a/src/physics_common.F90 b/src/physics_common.F90 index 63bf62ce62..0d5b9fb324 100644 --- a/src/physics_common.F90 +++ b/src/physics_common.F90 @@ -1,33 +1,20 @@ module physics_common - use constants + use, intrinsic :: ISO_C_BINDING + use particle_header, only: Particle - use random_lcg, only: prn - use settings, only: weight_cutoff, weight_survive implicit none -contains - !=============================================================================== -! RUSSIAN_ROULETTE +! RUSSIAN_ROULETTE FROM C !=============================================================================== - subroutine russian_roulette(p) - - type(Particle), intent(inout) :: p - - if (p % wgt < weight_cutoff) then - if (prn() < p % wgt / weight_survive) then - p % wgt = weight_survive - p % last_wgt = p % wgt - else - p % wgt = ZERO - p % last_wgt = ZERO - p % alive = .false. - end if - end if - - end subroutine russian_roulette + interface + subroutine russian_roulette(p) bind(C) + import Particle + type(Particle), intent(inout) :: p + end subroutine russian_roulette + end interface end module physics_common diff --git a/src/physics_common.cpp b/src/physics_common.cpp new file mode 100644 index 0000000000..02c77ae7d9 --- /dev/null +++ b/src/physics_common.cpp @@ -0,0 +1,26 @@ +#include "openmc/physics_common.h" + +#include "openmc/settings.h" +#include "openmc/random_lcg.h" + +namespace openmc { + +//============================================================================== +// RUSSIAN_ROULETTE +//============================================================================== + +void russian_roulette(Particle* p) +{ + if (p->wgt < settings::weight_cutoff) { + if (prn() < p->wgt / settings::weight_survive) { + p->wgt = settings::weight_survive; + p->last_wgt = p->wgt; + } else { + p->wgt = 0.; + p->last_wgt = 0.; + p->alive = false; + } + } +} + +} //namespace openmc From f3fb3d197991a1e6d5c0f754a73d1956a6cf7124 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Mon, 8 Oct 2018 14:01:14 -0400 Subject: [PATCH 2/7] replaced physics_mg.F90:scatter with the c++ version --- CMakeLists.txt | 1 + src/mgxs_interface.F90 | 2 +- src/physics_mg.F90 | 29 +++++++++-------------------- 3 files changed, 11 insertions(+), 21 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e62c413476..135a724311 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -424,6 +424,7 @@ add_library(libopenmc SHARED src/output.cpp src/particle.cpp src/physics_common.cpp + src/physics_mg.cpp src/plot.cpp src/position.cpp src/pugixml/pugixml_c.cpp diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90 index 9dccc92301..55a18afaba 100644 --- a/src/mgxs_interface.F90 +++ b/src/mgxs_interface.F90 @@ -156,7 +156,7 @@ module mgxs_interface real(8), allocatable :: energy_bins(:) ! Midpoint of the energy group structure - real(8), allocatable :: energy_bin_avg(:) + real(C_DOUBLE), allocatable :: energy_bin_avg(:) ! Energy group structure with increasing energy real(C_DOUBLE), allocatable, target :: rev_energy_bins(:) diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 1d60b371ee..8d0c976180 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -20,6 +20,14 @@ module physics_mg implicit none + interface + subroutine scatter(p, energy_bin_avg) bind(C) + import Particle, C_DOUBLE + type(Particle), intent(inout) :: p + real(C_DOUBLE), intent(in) :: energy_bin_avg(*) + end subroutine scatter + end interface + contains !=============================================================================== @@ -85,7 +93,7 @@ contains ! Sample a scattering reaction and determine the secondary energy of the ! exiting neutron - call scatter(p) + call scatter(p, energy_bin_avg) ! Play russian roulette if survival biasing is turned on if (survival_biasing) then @@ -133,25 +141,6 @@ contains end subroutine absorption -!=============================================================================== -! SCATTER -!=============================================================================== - - subroutine scatter(p) - - type(Particle), intent(inout) :: p - - call sample_scatter_c(p % material, p % last_g, p % g, p % mu, & - p % wgt, p % coord(1) % uvw) - - ! Update energy value for downstream compatability (in tallying) - p % E = energy_bin_avg(p % g) - - ! Set event component - p % event = EVENT_SCATTER - - end subroutine scatter - !=============================================================================== ! CREATE_FISSION_SITES determines the average total, prompt, and delayed ! neutrons produced from fission and creates appropriate bank sites. From 960f71e276c753ec056978b182cff751728f81cc Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Mon, 8 Oct 2018 14:01:52 -0400 Subject: [PATCH 3/7] whoops, including files --- include/openmc/physics_mg.h | 20 ++++++++++++++++++++ src/physics_mg.cpp | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+) create mode 100644 include/openmc/physics_mg.h create mode 100644 src/physics_mg.cpp diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h new file mode 100644 index 0000000000..e570de917a --- /dev/null +++ b/include/openmc/physics_mg.h @@ -0,0 +1,20 @@ +//! \file physics_mg.h +//! Methods needed to perform the collision physics for multi-group mode + +#ifndef OPENMC_PHYSICS_MG_H +#define OPENMC_PHYSICS_MG_H + +#include "openmc/particle.h" + +namespace openmc { + +//============================================================================== +// SCATTER +//============================================================================== + +//! \brief Samples the scattering event +extern "C" void +scatter(Particle* p, const double* energy_bin_avg); + +} // namespace openmc +#endif // OPENMC_PHYSICS_MG_H diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp new file mode 100644 index 0000000000..27379a2d6c --- /dev/null +++ b/src/physics_mg.cpp @@ -0,0 +1,32 @@ +#include "openmc/physics_mg.h" + +#include "openmc/constants.h" +#include "openmc/math_functions.h" +#include "openmc/mgxs_interface.h" + +namespace openmc { + +void scatter(Particle* p, const double* energy_bin_avg) +{ + // Adjust indices for Fortran to C++ indexing + // TODO: Remove when no longer needed + int gin = p->last_g - 1; + int gout = p->g - 1; + int i_mat = p->material - 1; + macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt); + + // Adjust return value for fortran indexing + // TODO: Remove when no longer needed + p->g = gout + 1; + + // Rotate the angle + rotate_angle_c(p->coord[0].uvw, p->mu, nullptr); + + // Update energy value for downstream compatability (in tallying) + p->E = energy_bin_avg[gout]; + + // Set event component + p->event = EVENT_SCATTER; +} + +} //namespace openmc From ef9908e774a9c1ef6105205548ed810a3c3ee5cb Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Tue, 9 Oct 2018 20:35:25 -0400 Subject: [PATCH 4/7] Added create_fission_sites and absorption to physics_mg.cpp --- include/openmc/capi.h | 2 + include/openmc/physics_mg.h | 12 +++ src/physics_mg.F90 | 186 +++++------------------------------ src/physics_mg.cpp | 145 ++++++++++++++++++++++++++- src/tallies/tally_header.F90 | 2 +- 5 files changed, 183 insertions(+), 164 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index c48b212bd7..57d817d23d 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -149,6 +149,8 @@ extern "C" { extern int32_t n_tallies; extern int32_t n_universes; extern bool openmc_simulation_initialized; + extern double global_tally_absorption; +#pragma omp threadprivate(global_tally_absorption) // Variables that are shared by necessity (can be removed from public header // later) diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index e570de917a..7fed561f34 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -4,7 +4,9 @@ #ifndef OPENMC_PHYSICS_MG_H #define OPENMC_PHYSICS_MG_H +#include "openmc/capi.h" #include "openmc/particle.h" +#include "openmc/nuclide.h" namespace openmc { @@ -16,5 +18,15 @@ namespace openmc { extern "C" void scatter(Particle* p, const double* energy_bin_avg); +//! \brief Determines the average total, prompt and delayed neutrons produced +//! from fission and creates the appropriate bank sites. +extern "C" void +create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, + int64_t& bank_array_size, MaterialMacroXS& material_xs); + +//! \brief Handles an absorption event +extern "C" void +absorption(Particle* p, MaterialMacroXS& material_xs); + } // namespace openmc #endif // OPENMC_PHYSICS_MG_H diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 8d0c976180..d0a69c8631 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -9,7 +9,7 @@ module physics_mg use math, only: rotate_angle use mgxs_interface use message_passing - use nuclide_header, only: material_xs + use nuclide_header, only: MaterialMacroXS, material_xs use particle_header use physics_common use random_lcg, only: prn @@ -26,6 +26,22 @@ module physics_mg type(Particle), intent(inout) :: p real(C_DOUBLE), intent(in) :: energy_bin_avg(*) end subroutine scatter + + subroutine create_fission_sites(p, bank_array, size_bank, bank_array_size, & + material_xs) bind(C) + import Particle, Bank, C_INT64_T, MaterialMacroXS + type(Particle), intent(inout) :: p + type(Bank), intent(inout) :: bank_array(*) + integer(C_INT64_T), intent(inout) :: size_bank + integer(C_INT64_T), intent(in) :: bank_array_size + type(MaterialMacroXS), intent(in) :: material_xs + end subroutine create_fission_sites + + subroutine absorption(p, material_xs) bind(C) + import Particle, MaterialMacroXS + type(Particle), intent(inout) :: p + type(MaterialMacroXS), intent(in) :: material_xs + end subroutine absorption end interface contains @@ -75,9 +91,13 @@ contains if (mat % fissionable) then if (run_mode == MODE_EIGENVALUE) then - call create_fission_sites(p, fission_bank, n_bank) + call create_fission_sites(p, fission_bank, n_bank, & + size(fission_bank, KIND=C_INT64_T), & + material_xs) elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then - call create_fission_sites(p, p % secondary_bank, p % n_secondary) + call create_fission_sites(p, p % secondary_bank, p % n_secondary, & + size(p % secondary_bank, KIND=C_INT64_T), & + material_xs) end if end if @@ -85,7 +105,7 @@ contains ! weight of the particle. Otherwise, it checks to see if absorption occurs if (material_xs % absorption > ZERO) then - call absorption(p) + call absorption(p, material_xs) else p % absorb_wgt = ZERO end if @@ -103,162 +123,4 @@ contains end subroutine sample_reaction -!=============================================================================== -! ABSORPTION -!=============================================================================== - - subroutine absorption(p) - - type(Particle), intent(inout) :: p - - if (survival_biasing) then - ! Determine weight absorbed in survival biasing - p % absorb_wgt = (p % wgt * & - material_xs % absorption / material_xs % total) - - ! Adjust weight of particle by probability of absorption - p % wgt = p % wgt - p % absorb_wgt - p % last_wgt = p % wgt - - ! Score implicit absorption estimate of keff -!$omp atomic - global_tallies(RESULT_VALUE, K_ABSORPTION) = & - global_tallies(RESULT_VALUE, K_ABSORPTION) + p % absorb_wgt * & - material_xs % nu_fission / material_xs % absorption - else - ! See if disappearance reaction happens - if (material_xs % absorption > prn() * material_xs % total) then - ! Score absorption estimate of keff -!$omp atomic - global_tallies(RESULT_VALUE, K_ABSORPTION) = & - global_tallies(RESULT_VALUE, K_ABSORPTION) + p % wgt * & - material_xs % nu_fission / material_xs % absorption - - p % alive = .false. - p % event = EVENT_ABSORB - end if - end if - - end subroutine absorption - -!=============================================================================== -! CREATE_FISSION_SITES determines the average total, prompt, and delayed -! neutrons produced from fission and creates appropriate bank sites. -!=============================================================================== - - subroutine create_fission_sites(p, bank_array, size_bank) - type(Particle), intent(inout) :: p - type(Bank), intent(inout) :: bank_array(:) - integer(8), intent(inout) :: size_bank - - integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born - integer :: i ! loop index - integer :: dg ! delayed group - integer :: gout ! group out - integer :: nu ! actual number of neutrons produced - real(8) :: nu_t ! total nu - real(8) :: mu ! fission neutron angular cosine - real(8) :: phi ! fission neutron azimuthal angle - real(8) :: weight ! weight adjustment for ufs method - - interface - function ufs_get_weight(p) result(weight) bind(C) - import Particle, C_DOUBLE - type(Particle), intent(in) :: p - real(C_DOUBLE) :: WEIGHT - end function - end interface - - ! TODO: Heat generation from fission - - ! If uniform fission source weighting is turned on, we increase of decrease - ! the expected number of fission sites produced - - if (ufs) then - weight = ufs_get_weight(p) - else - weight = ONE - end if - - ! Determine expected number of neutrons produced - nu_t = p % wgt / keff * weight * & - material_xs % nu_fission / material_xs % total - - ! Sample number of neutrons produced - if (prn() > nu_t - int(nu_t)) then - nu = int(nu_t) - else - nu = int(nu_t) + 1 - end if - - ! Check for bank size getting hit. For fixed source calculations, this is a - ! fatal error. For eigenvalue calculations, it just means that k-effective - ! was too high for a single batch. - if (size_bank + nu > size(bank_array)) then - if (run_mode == MODE_FIXEDSOURCE) then - call fatal_error("Secondary particle bank size limit reached. If you & - &are running a subcritical multiplication problem, k-effective & - &may be too close to one.") - else - if (master) call warning("Maximum number of sites in fission bank & - &reached. This can result in irreproducible results using different & - &numbers of processes/threads.") - end if - end if - - ! Bank source neutrons - if (nu == 0 .or. size_bank == size(bank_array)) return - - ! Initialize counter of delayed neutrons encountered for each delayed group - ! to zero. - nu_d(:) = 0 - - p % fission = .true. ! Fission neutrons will be banked - do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4) - ! Bank source neutrons by copying particle data - bank_array(i) % xyz = p % coord(1) % xyz - - ! Set particle as neutron - bank_array(i) % particle = NEUTRON - - ! Set weight of fission bank site - bank_array(i) % wgt = ONE/weight - - ! Sample cosine of angle -- fission neutrons are treated as being emitted - ! isotropically. - mu = TWO * prn() - ONE - - ! Sample azimuthal angle uniformly in [0,2*pi) - phi = TWO * PI * prn() - bank_array(i) % uvw(1) = mu - bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) - - ! Sample secondary energy distribution for fission reaction and set energy - ! in fission bank - call sample_fission_energy_c(p % material, p % g, dg, gout) - - bank_array(i) % E = real(gout, 8) - bank_array(i) % delayed_group = dg - - ! Set delayed group on particle too - p % delayed_group = dg - - ! Increment the number of neutrons born delayed - if (p % delayed_group > 0) then - nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1 - end if - - end do - - ! increment number of bank sites - size_bank = min(size_bank + nu, int(size(bank_array),8)) - - ! Store total weight banked for analog fission tallies - p % n_bank = nu - p % wgt_bank = nu/weight - p % n_delayed_bank(:) = nu_d(:) - - end subroutine create_fission_sites - end module physics_mg diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 27379a2d6c..bd7bcec11b 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -1,12 +1,23 @@ #include "openmc/physics_mg.h" +#include +#include + +#include "xtensor/xarray.hpp" + #include "openmc/constants.h" +#include "openmc/eigenvalue.h" +#include "openmc/error.h" #include "openmc/math_functions.h" +#include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" +#include "openmc/random_lcg.h" +#include "openmc/settings.h" namespace openmc { -void scatter(Particle* p, const double* energy_bin_avg) +void +scatter(Particle* p, const double* energy_bin_avg) { // Adjust indices for Fortran to C++ indexing // TODO: Remove when no longer needed @@ -29,4 +40,136 @@ void scatter(Particle* p, const double* energy_bin_avg) p->event = EVENT_SCATTER; } +void +create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, + int64_t& bank_array_size, MaterialMacroXS& material_xs) +{ + // TODO: Heat generation from fission + + // If uniform fission source weighting is turned on, we increase or decrease + // the expected number of fission sites produced + double weight; + if (settings::ufs_on) { + weight = ufs_get_weight(p); + } else { + weight = 1.; + } + + // Determine the expected number of neutrons produced + double nu_t = p->wgt / openmc_keff * weight * material_xs.nu_fission / + material_xs.total; + + // Sample the number of neutrons produced + int nu = static_cast(nu_t); + if (prn() <= (nu_t - int(nu_t))) { + nu++; + } + + // Check for the bank size getting hit. For fixed source calculations, this + // is a fatal error; for eigenvalue calculations, it just means that k-eff + // was too high for a single batch. + if (size_bank + nu > bank_array_size) { + if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { + throw std::runtime_error{"Secondary particle bank size limit reached." + " If you are running a subcritical multiplication problem," + " k-effective may be too close to one."}; + } else { + if (mpi::master) { + std::stringstream msg; + msg << "Maximum number of sites in fission bank reached. This can" + " result in irreproducible results using different numbers of" + " processes/threads."; + warning(msg); + } + } + } + + // Begin banking the source neutrons + // First, if our bank is full then don't continue + if ((nu == 0) || (size_bank == bank_array_size)) return; + + // Initialize the counter of delayed neutrons encountered for each delayed + // group. + double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + p->fission = true; + // TODO: +1 for start? + for (size_t i = static_cast(size_bank); + i < static_cast(std::min(size_bank + nu, bank_array_size)); i++) { + // Bank source neutrons by copying the particle data + bank_array[i].xyz[0] = p->coord[0].xyz[0]; + bank_array[i].xyz[1] = p->coord[0].xyz[1]; + bank_array[i].xyz[2] = p->coord[0].xyz[2]; + + // Set that the bank particle is a neutron + bank_array[i].particle = static_cast(ParticleType::neutron); + + // Set the weight of the fission bank site + bank_array[i].wgt = 1. / weight; + + // Sample the cosine of the angle, assuming fission neutrons are emitted + // isotropically + double mu = 2. * prn() - 1.; + + // Sample the azimuthal angle uniformly in [0, 2.pi) + double phi = 2. * PI * prn(); + bank_array[i].uvw[0] = mu; + bank_array[i].uvw[1] = std::sqrt(1. - mu * mu) * std::cos(phi); + bank_array[i].uvw[2] = std::sqrt(1. - mu * mu) * std::sin(phi); + + // Sample secondary energy distribution for the fission reaction and set + // the energy in the fission bank + int dg; + int gout; + macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout); + bank_array[i].E = static_cast(gout + 1); + bank_array[i].delayed_group = dg + 1; + + // Set the delayed group on the particle as well + p->delayed_group = dg + 1; + + // Increment the number of neutrons born delayed + if (p->delayed_group > 0) { + nu_d[dg]++; + } + } + + // Increment number of bank sites + size_bank = std::min(size_bank + nu, bank_array_size); + + // Store the total weight banked for analog fission tallies + p->n_bank = nu; + p->wgt_bank = nu / weight; + for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { + p->n_delayed_bank[d] = nu_d[d]; + } +} + +void +absorption(Particle* p, MaterialMacroXS& material_xs) +{ + if (settings::survival_biasing) { + // Determine weight absorbed in survival biasing + p->absorb_wgt = p->wgt * material_xs.absorption / material_xs.total; + + // Adjust weight of particle by the probability of absorption + p->wgt -= p->absorb_wgt; + p->last_wgt = p->wgt; + + // Score implicit absorpion estimate of keff +#pragma omp atomic + global_tally_absorption += p->absorb_wgt * material_xs.nu_fission / + material_xs.absorption; + } else { + if (material_xs.absorption > prn() * material_xs.total) { +#pragma omp atomic + global_tally_absorption += p->wgt * material_xs.nu_fission / + material_xs.absorption; + p->alive = false; + p->event = EVENT_ABSORB; + } + + } +} + } //namespace openmc diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 367214a930..e9b8e53de9 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -137,7 +137,7 @@ module tally_header ! invalidation. Thus, we use threadprivate variables to accumulate global ! tallies and then reduce at the end of a generation. real(C_DOUBLE), public :: global_tally_collision = ZERO - real(C_DOUBLE), public :: global_tally_absorption = ZERO + real(C_DOUBLE), public, bind(C) :: global_tally_absorption = ZERO real(C_DOUBLE), public :: global_tally_tracklength = ZERO real(C_DOUBLE), public :: global_tally_leakage = ZERO !$omp threadprivate(global_tally_collision, global_tally_absorption, & From e3e388bb6d16ef6c88f1ee50290a72c720565889 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Fri, 12 Oct 2018 10:40:35 -0400 Subject: [PATCH 5/7] Completed replacement of physics_mg.F90 with C++ --- CMakeLists.txt | 1 - include/openmc/material.h | 1 + include/openmc/physics_mg.h | 24 +++++-- src/input_xml.F90 | 2 +- src/material.cpp | 7 ++ src/material_header.F90 | 31 ++++++++- src/mgxs_data.F90 | 3 +- src/physics_mg.F90 | 126 ------------------------------------ src/physics_mg.cpp | 67 ++++++++++++++++++- src/tracking.F90 | 18 +++++- 10 files changed, 140 insertions(+), 140 deletions(-) delete mode 100644 src/physics_mg.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 135a724311..b72bb57eb9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -345,7 +345,6 @@ add_library(libopenmc SHARED src/photon_physics.F90 src/physics_common.F90 src/physics.F90 - src/physics_mg.F90 src/plot.F90 src/plot_header.F90 src/progress_header.F90 diff --git a/include/openmc/material.h b/include/openmc/material.h index 1af6f85976..909215ee37 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -26,6 +26,7 @@ class Material public: int32_t id; //!< Unique ID double volume_ {-1.0}; //!< Volume in [cm^3] + bool fissionable {false}; //!< Does this material contain fissionable nuclides //! \brief Default temperature for cells containing this material. //! diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 7fed561f34..fdfebf71e4 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -14,19 +14,33 @@ namespace openmc { // SCATTER //============================================================================== -//! \brief Samples the scattering event +//! \brief samples particle behavior after a collision event. extern "C" void +collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size, + const double* energy_bin_avg, const MaterialMacroXS& material_xs); + +//! \brief samples a reaction type. +//! +//! Note that there is special logic when suvival biasing is turned on since +//! fission and disappearance are treated implicitly. +void +sample_reaction(Particle* p, Bank* fission_bank, + const int64_t fission_bank_size, const double* energy_bin_avg, + const MaterialMacroXS& material_xs); + +//! \brief Samples the scattering event +void scatter(Particle* p, const double* energy_bin_avg); //! \brief Determines the average total, prompt and delayed neutrons produced //! from fission and creates the appropriate bank sites. -extern "C" void +void create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, - int64_t& bank_array_size, MaterialMacroXS& material_xs); + const int64_t& bank_array_size, const MaterialMacroXS& material_xs); //! \brief Handles an absorption event -extern "C" void -absorption(Particle* p, MaterialMacroXS& material_xs); +void +absorption(Particle* p, const MaterialMacroXS& material_xs); } // namespace openmc #endif // OPENMC_PHYSICS_MG_H diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 420ae5bbb8..12b97d9407 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2986,7 +2986,7 @@ contains ! Check if material is fissionable if (nuclides(materials(i) % nuclide(j)) % fissionable) then - materials(i) % fissionable = .true. + call materials(i) % set_fissionable(logical(.true., C_BOOL)) end if end do diff --git a/src/material.cpp b/src/material.cpp index 1a0c2d2078..2bd9228753 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -122,6 +122,13 @@ extern "C" { material_map[id] = index - 1; } + bool material_fissionable(Material* mat) {return mat->fissionable;} + + void material_set_fissionable(Material* mat, bool fissionable, int32_t index) + { + mat->fissionable = fissionable; + } + void extend_materials_c(int32_t n) { materials.reserve(materials.size() + n); diff --git a/src/material_header.F90 b/src/material_header.F90 index 71e6ec5f1a..a067269ce2 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -51,6 +51,20 @@ module material_header integer(C_INT32_T), intent(in), value :: index end subroutine material_set_id_c + function material_fissionable_c(mat_ptr) & + bind(C, name='material_fissionable') result(fissionable) + import C_PTR, C_BOOL + type(C_PTR), intent(in), value :: mat_ptr + logical(C_BOOL) :: fissionable + end function material_fissionable_c + + subroutine material_set_fissionable_c(mat_ptr, fissionable) & + bind(C, name='material_set_fissionable') + import C_PTR, C_BOOL + type(C_PTR), intent(in), value :: mat_ptr + logical(C_BOOL), intent(in), value :: fissionable + end subroutine material_set_fissionable_c + subroutine extend_materials_c(n) bind(C) import C_INT32_T integer(C_INT32_T), intent(in), value :: n @@ -95,7 +109,6 @@ module material_header character(20), allocatable :: sab_names(:) ! name of S(a,b) table ! Does this material contain fissionable nuclides? Is it depletable? - logical :: fissionable = .false. logical :: depletable = .false. ! enforce isotropic scattering in lab for specific nuclides @@ -105,6 +118,8 @@ module material_header contains procedure :: id => material_id procedure :: set_id => material_set_id + procedure :: fissionable => material_fissionable + procedure :: set_fissionable => material_set_fissionable procedure :: set_density => material_set_density procedure :: init_nuclide_index => material_init_nuclide_index procedure :: assign_sab_tables => material_assign_sab_tables @@ -139,6 +154,18 @@ contains call material_set_id_c(this % ptr, id, index) end subroutine material_set_id + function material_fissionable(this) result(fissionable) + class(Material), intent(in) :: this + logical(C_BOOL) :: fissionable + fissionable = material_fissionable_c(this % ptr) + end function material_fissionable + + subroutine material_set_fissionable(this, fissionable) + class(Material), intent(in) :: this + logical(C_BOOL), intent(in) :: fissionable + call material_set_fissionable_c(this % ptr, fissionable) + end subroutine material_set_fissionable + function material_set_density(this, density) result(err) class(Material), intent(inout) :: this real(8), intent(in) :: density @@ -685,7 +712,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= size(materials)) then - fissionable = materials(index) % fissionable + fissionable = materials(index) % fissionable() err = 0 else err = E_OUT_OF_BOUNDS diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 03a0e9402f..e6eb5a3da5 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -90,7 +90,8 @@ contains end if end do NUCLIDE_LOOP - mat % fissionable = query_fissionable_c(mat % n_nuclides, mat % nuclide) + call mat % set_fissionable(query_fissionable_c(mat % n_nuclides, & + mat % nuclide)) end do MATERIAL_LOOP diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 deleted file mode 100644 index d0a69c8631..0000000000 --- a/src/physics_mg.F90 +++ /dev/null @@ -1,126 +0,0 @@ -module physics_mg - ! This module contains the multi-group specific physics routines so as to not - ! hinder performance of the CE versions with multiple if-thens. - - use bank_header - use constants - use error, only: fatal_error, warning, write_message - use material_header, only: Material, materials - use math, only: rotate_angle - use mgxs_interface - use message_passing - use nuclide_header, only: MaterialMacroXS, material_xs - use particle_header - use physics_common - use random_lcg, only: prn - use settings - use simulation_header - use string, only: to_str - use tally_header - - implicit none - - interface - subroutine scatter(p, energy_bin_avg) bind(C) - import Particle, C_DOUBLE - type(Particle), intent(inout) :: p - real(C_DOUBLE), intent(in) :: energy_bin_avg(*) - end subroutine scatter - - subroutine create_fission_sites(p, bank_array, size_bank, bank_array_size, & - material_xs) bind(C) - import Particle, Bank, C_INT64_T, MaterialMacroXS - type(Particle), intent(inout) :: p - type(Bank), intent(inout) :: bank_array(*) - integer(C_INT64_T), intent(inout) :: size_bank - integer(C_INT64_T), intent(in) :: bank_array_size - type(MaterialMacroXS), intent(in) :: material_xs - end subroutine create_fission_sites - - subroutine absorption(p, material_xs) bind(C) - import Particle, MaterialMacroXS - type(Particle), intent(inout) :: p - type(MaterialMacroXS), intent(in) :: material_xs - end subroutine absorption - end interface - -contains - -!=============================================================================== -! COLLISION_MG samples a nuclide and reaction and then calls the appropriate -! routine for that reaction -!=============================================================================== - - subroutine collision_mg(p) - - type(Particle), intent(inout) :: p - - ! Add to collision counter for particle - p % n_collision = p % n_collision + 1 - - ! Sample nuclide/reaction for the material the particle is in - call sample_reaction(p) - - ! Display information about collision - if (verbosity >= 10 .or. trace) then - call write_message(" " // "Energy Group = " // trim(to_str(p % g))) - end if - - end subroutine collision_mg - -!=============================================================================== -! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for -! each nuclide within a material and then samples a reaction for that nuclide -! and calls the appropriate routine to process the physics. Note that there is -! special logic when suvival biasing is turned on since fission and -! disappearance are treated implicitly. -!=============================================================================== - - subroutine sample_reaction(p) - - type(Particle), intent(inout) :: p - - type(Material), pointer :: mat - - mat => materials(p % material) - - ! Create fission bank sites. Note that while a fission reaction is sampled, - ! it never actually "happens", i.e. the weight of the particle does not - ! change when sampling fission sites. The following block handles all - ! absorption (including fission) - - if (mat % fissionable) then - if (run_mode == MODE_EIGENVALUE) then - call create_fission_sites(p, fission_bank, n_bank, & - size(fission_bank, KIND=C_INT64_T), & - material_xs) - elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then - call create_fission_sites(p, p % secondary_bank, p % n_secondary, & - size(p % secondary_bank, KIND=C_INT64_T), & - material_xs) - end if - end if - - ! If survival biasing is being used, the following subroutine adjusts the - ! weight of the particle. Otherwise, it checks to see if absorption occurs - - if (material_xs % absorption > ZERO) then - call absorption(p, material_xs) - else - p % absorb_wgt = ZERO - end if - if (.not. p % alive) return - - ! Sample a scattering reaction and determine the secondary energy of the - ! exiting neutron - call scatter(p, energy_bin_avg) - - ! Play russian roulette if survival biasing is turned on - if (survival_biasing) then - call russian_roulette(p) - if (.not. p % alive) return - end if - - end subroutine sample_reaction - -end module physics_mg diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index bd7bcec11b..8651e29f18 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -8,14 +8,76 @@ #include "openmc/constants.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/material.h" #include "openmc/math_functions.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" +#include "openmc/physics_common.h" #include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/simulation.h" namespace openmc { +void +collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size, + const double* energy_bin_avg, const MaterialMacroXS& material_xs) +{ + // Add to the collision counter for the particle + p->n_collision++; + + // Sample the reaction type + sample_reaction(p, fission_bank, fission_bank_size, energy_bin_avg, + material_xs); + + // Display information about collision + if ((settings::verbosity >= 10) || (openmc_trace)) { + std::stringstream msg; + msg << " Energy Group = " << p->g; + write_message(msg, 1); + } +} + +void +sample_reaction(Particle* p, Bank* fission_bank, + const int64_t fission_bank_size, const double* energy_bin_avg, + const MaterialMacroXS& material_xs) +{ + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) + + if (materials[p->material - 1]->fissionable) { + if (settings::run_mode == RUN_MODE_EIGENVALUE) { + create_fission_sites(p, fission_bank, n_bank, fission_bank_size, + material_xs); + } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && + (settings::create_fission_neutrons)) { + create_fission_sites(p, p->secondary_bank, p->n_secondary, + MAX_SECONDARY, material_xs); + } + } + + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs. + if (material_xs.absorption > 0.) { + absorption(p, material_xs); + } else { + p->absorb_wgt = 0.; + } + if (!p->alive) return; + + // Sample a scattering event to determine the energy of the exiting neutron + scatter(p, energy_bin_avg); + + // Play Russian roulette if survival biasing is turned on + if (settings::survival_biasing) { + russian_roulette(p); + if (!p->alive) return; + } +} + void scatter(Particle* p, const double* energy_bin_avg) { @@ -42,7 +104,8 @@ scatter(Particle* p, const double* energy_bin_avg) void create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, - int64_t& bank_array_size, MaterialMacroXS& material_xs) + const int64_t& bank_array_size, + const MaterialMacroXS& material_xs) { // TODO: Heat generation from fission @@ -146,7 +209,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, } void -absorption(Particle* p, MaterialMacroXS& material_xs) +absorption(Particle* p, const MaterialMacroXS& material_xs) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing diff --git a/src/tracking.F90 b/src/tracking.F90 index 8cf2fc05c7..02f10f7cc8 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -17,7 +17,7 @@ module tracking use nuclide_header use particle_header use physics, only: collision - use physics_mg, only: collision_mg + ! use physics_mg, only: collision_mg use random_lcg, only: prn, prn_set_stream use settings use simulation_header @@ -33,6 +33,19 @@ module tracking implicit none + interface + subroutine collision_mg(p, fission_bank, fission_bank_size, & + energy_bin_avg, material_xs) bind(C) + import Particle, Bank, C_INT64_T, C_DOUBLE, MaterialMacroXS + type(Particle), intent(inout) :: p + type(Bank), intent(inout) :: fission_bank(*) + integer(C_INT64_T), value, intent(in) :: fission_bank_size + real(C_DOUBLE), intent(in) :: energy_bin_avg(*) + type(MaterialMacroXS), intent(in) :: material_xs + end subroutine collision_mg + + end interface + contains !=============================================================================== @@ -226,7 +239,8 @@ contains if (run_CE) then call collision(p) else - call collision_mg(p) + call collision_mg(p, fission_bank, & + size(fission_bank, kind=C_INT64_T), energy_bin_avg, material_xs) end if ! Score collision estimator tallies -- this is done after a collision From a9641b8e2adae6e69dcbcd95c6cdad37a9347f9d Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Fri, 12 Oct 2018 11:51:23 -0400 Subject: [PATCH 6/7] Cleaned up the C and Fortran interfaces; includes removing interfaces which are no longer needed and removing the _c from C functions which no longer have Fortran links --- include/openmc/math_functions.h | 8 +++--- include/openmc/mgxs_interface.h | 10 +++----- include/openmc/physics_mg.h | 13 +++++----- openmc/capi/math.py | 23 +++++++++-------- src/distribution.cpp | 4 +-- src/distribution_energy.cpp | 4 +-- src/math.F90 | 44 --------------------------------- src/math_functions.cpp | 12 ++++----- src/mgxs_interface.F90 | 20 --------------- src/mgxs_interface.cpp | 30 ---------------------- src/physics_mg.cpp | 16 ++++++------ src/scattdata.cpp | 12 ++++----- src/secondary_nbody.cpp | 4 +-- src/tracking.F90 | 10 +++----- 14 files changed, 54 insertions(+), 156 deletions(-) diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index a3c9c07cf0..8bfa118995 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -22,7 +22,7 @@ namespace openmc { //! @return The requested percentile //============================================================================== -extern "C" double normal_percentile_c(double p); +extern "C" double normal_percentile(double p); //============================================================================== //! Calculate the percentile of the Student's t distribution with a specified @@ -58,7 +58,7 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]); //! evaluated at x //============================================================================== -extern "C" double evaluate_legendre_c(int n, const double data[], double x); +extern "C" double evaluate_legendre(int n, const double data[], double x); //============================================================================== //! Calculate the n-th order real spherical harmonics for a given angle (in @@ -144,7 +144,7 @@ Direction rotate_angle(Direction u, double mu, double* phi); //! @result The sampled outgoing energy //============================================================================== -extern "C" double maxwell_spectrum_c(double T); +extern "C" double maxwell_spectrum(double T); //============================================================================== //! Samples an energy from a Watt energy-dependent fission distribution. @@ -159,7 +159,7 @@ extern "C" double maxwell_spectrum_c(double T); //! @result The sampled outgoing energy //============================================================================== -extern "C" double watt_spectrum_c(double a, double b); +extern "C" double watt_spectrum(double a, double b); //============================================================================== //! Doppler broadens the windowed multipole curvefit. diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index 02a43fece4..e0addf5bc0 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -18,6 +18,9 @@ extern std::vector nuclides_MG; extern std::vector macro_xs; extern "C" int num_energy_groups; +//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also +// bring over energy_bin_avg, energy_bins, etc, as vectors) + //============================================================================== // Mgxs data loading interface methods //============================================================================== @@ -44,13 +47,6 @@ extern "C" void calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs); -extern "C" void -sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt, - double uvw[3]); - -extern "C" void -sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout); - extern "C" double get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg); diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index fdfebf71e4..09a8cc617d 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -10,22 +10,21 @@ namespace openmc { -//============================================================================== -// SCATTER -//============================================================================== +//TODO: Remove energy_bin_avg and material_xs parameters when they reside on +// the C-side this should happen after materials, physics, input, and tallies +// are brought over //! \brief samples particle behavior after a collision event. extern "C" void -collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size, - const double* energy_bin_avg, const MaterialMacroXS& material_xs); +collision_mg(Particle* p, const double* energy_bin_avg, + const MaterialMacroXS& material_xs); //! \brief samples a reaction type. //! //! Note that there is special logic when suvival biasing is turned on since //! fission and disappearance are treated implicitly. void -sample_reaction(Particle* p, Bank* fission_bank, - const int64_t fission_bank_size, const double* energy_bin_avg, +sample_reaction(Particle* p, const double* energy_bin_avg, const MaterialMacroXS& material_xs); //! \brief Samples the scattering event diff --git a/openmc/capi/math.py b/openmc/capi/math.py index 390ee466ff..a8ea849b3e 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -12,8 +12,8 @@ _dll.t_percentile_c.argtypes = [c_double, c_int] _dll.calc_pn_c.restype = None _dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)] -_dll.evaluate_legendre_c.restype = c_double -_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double] +_dll.evaluate_legendre.restype = c_double +_dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double] _dll.calc_rn_c.restype = None _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)] @@ -27,11 +27,11 @@ _dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)] _dll.rotate_angle_c.restype = None _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, POINTER(c_double)] -_dll.maxwell_spectrum_c.restype = c_double -_dll.maxwell_spectrum_c.argtypes = [c_double] +_dll.maxwell_spectrum.restype = c_double +_dll.maxwell_spectrum.argtypes = [c_double] -_dll.watt_spectrum_c.restype = c_double -_dll.watt_spectrum_c.argtypes = [c_double, c_double] +_dll.watt_spectrum.restype = c_double +_dll.watt_spectrum.argtypes = [c_double, c_double] _dll.broaden_wmp_polynomials_c.restype = None _dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int, @@ -100,9 +100,8 @@ def evaluate_legendre(data, x): """ data_arr = np.array(data, dtype=np.float64) - return _dll.evaluate_legendre_c(len(data), - data_arr.ctypes.data_as(POINTER(c_double)), - x) + return _dll.evaluate_legendre(len(data), + data_arr.ctypes.data_as(POINTER(c_double)), x) def calc_rn(n, uvw): @@ -182,7 +181,7 @@ def calc_zn_rad(n, rho): zn_rad = np.zeros(num_bins, dtype=np.float64) _dll.calc_zn_rad_c(n, rho, zn_rad) return zn_rad - + def rotate_angle(uvw0, mu, phi=None): """ Rotates direction cosines through a polar angle whose cosine is @@ -231,7 +230,7 @@ def maxwell_spectrum(T): """ - return _dll.maxwell_spectrum_c(T) + return _dll.maxwell_spectrum(T) def watt_spectrum(a, b): @@ -251,7 +250,7 @@ def watt_spectrum(a, b): """ - return _dll.watt_spectrum_c(a, b) + return _dll.watt_spectrum(a, b) def broaden_wmp_polynomials(E, dopp, n): diff --git a/src/distribution.cpp b/src/distribution.cpp index 01b09e65db..e546b11ddd 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -90,7 +90,7 @@ Maxwell::Maxwell(pugi::xml_node node) double Maxwell::sample() const { - return maxwell_spectrum_c(theta_); + return maxwell_spectrum(theta_); } //============================================================================== @@ -110,7 +110,7 @@ Watt::Watt(pugi::xml_node node) double Watt::sample() const { - return watt_spectrum_c(a_, b_); + return watt_spectrum(a_, b_); } //============================================================================== diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index 863ffd8b08..e311e322d2 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -279,7 +279,7 @@ double MaxwellEnergy::sample(double E) const while (true) { // Sample maxwell fission spectrum - double E_out = maxwell_spectrum_c(theta); + double E_out = maxwell_spectrum(theta); // Accept energy based on restriction energy if (E_out <= E - u_) return E_out; @@ -343,7 +343,7 @@ double WattEnergy::sample(double E) const while (true) { // Sample energy-dependent Watt fission spectrum - double E_out = watt_spectrum_c(a, b); + double E_out = watt_spectrum(a, b); // Accept energy based on restriction energy if (E_out <= E - u_) return E_out; diff --git a/src/math.F90 b/src/math.F90 index 9964922c17..143e21e030 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -12,10 +12,7 @@ module math public :: calc_rn public :: calc_zn public :: calc_zn_rad - public :: evaluate_legendre public :: rotate_angle - public :: maxwell_spectrum - public :: watt_spectrum public :: faddeeva public :: w_derivative public :: broaden_wmp_polynomials @@ -42,16 +39,6 @@ module math real(C_DOUBLE), intent(out) :: pnx(n + 1) end subroutine calc_pn - pure function evaluate_legendre_c_intfc(n, data, x) & - bind(C, name='evaluate_legendre_c') result(val) - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), intent(in) :: data(n) - real(C_DOUBLE), value, intent(in) :: x - real(C_DOUBLE) :: val - end function evaluate_legendre_c_intfc - pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c') use ISO_C_BINDING implicit none @@ -85,23 +72,6 @@ module math real(C_DOUBLE), optional, intent(in) :: phi end subroutine rotate_angle_c_intfc - function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') & - result(E_out) - use ISO_C_BINDING - implicit none - real(C_DOUBLE), value, intent(in) :: T - real(C_DOUBLE) :: E_out - end function maxwell_spectrum - - function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') & - result(E_out) - use ISO_C_BINDING - implicit none - real(C_DOUBLE), value, intent(in) :: a - real(C_DOUBLE), value, intent(in) :: b - real(C_DOUBLE) :: E_out - end function watt_spectrum - subroutine broaden_wmp_polynomials(E, dopp, n, factors) & bind(C, name='broaden_wmp_polynomials_c') use ISO_C_BINDING @@ -157,20 +127,6 @@ module math contains -!=============================================================================== -! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients -! and the value of x -!=============================================================================== - - pure function evaluate_legendre(data, x) result(val) bind(C) - real(C_DOUBLE), intent(in) :: data(:) - real(C_DOUBLE), intent(in) :: x - real(C_DOUBLE) :: val - - val = evaluate_legendre_c_intfc(size(data) - 1, data, x) - - end function evaluate_legendre - !=============================================================================== ! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is ! mu and through an azimuthal angle sampled uniformly. Note that this is done diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 5c1a2f4dce..8d30ea84d7 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -6,7 +6,7 @@ namespace openmc { // Mathematical methods //============================================================================== -double normal_percentile_c(double p) { +double normal_percentile(double p) { constexpr double p_low = 0.02425; constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2, -2.759285104469687e2, 1.383577518672690e2, @@ -79,7 +79,7 @@ double t_percentile_c(double p, int df){ // 16 (4), pp. 1123-1132 (1987). double n = df; double k = 1. / (n - 2.); - double z = normal_percentile_c(p); + double z = normal_percentile(p); double z2 = z * z; t = std::sqrt(n * k) * (z + (z2 - 3.) * z * k / 4. + ((5. * z2 - 56.) * z2 + 75.) * z * k * k / 96. + (((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) * @@ -103,7 +103,7 @@ void calc_pn_c(int n, double x, double pnx[]) { } -double evaluate_legendre_c(int n, const double data[], double x) { +double evaluate_legendre(int n, const double data[], double x) { double pnx[n + 1]; double val = 0.0; calc_pn_c(n, x, pnx); @@ -658,7 +658,7 @@ Direction rotate_angle(Direction u, double mu, double* phi) } -double maxwell_spectrum_c(double T) { +double maxwell_spectrum(double T) { // Set the random numbers double r1 = prn(); double r2 = prn(); @@ -674,8 +674,8 @@ double maxwell_spectrum_c(double T) { } -double watt_spectrum_c(double a, double b) { - double w = maxwell_spectrum_c(a); +double watt_spectrum(double a, double b) { + double w = maxwell_spectrum(a); double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w); return E_out; diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90 index 55a18afaba..77f629b80e 100644 --- a/src/mgxs_interface.F90 +++ b/src/mgxs_interface.F90 @@ -62,26 +62,6 @@ module mgxs_interface real(C_DOUBLE), intent(inout) :: nu_fiss_xs end subroutine calculate_xs_c - subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C) - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: i_mat - integer(C_INT), value, intent(in) :: gin - integer(C_INT), intent(inout) :: gout - real(C_DOUBLE), intent(inout) :: mu - real(C_DOUBLE), intent(inout) :: wgt - real(C_DOUBLE), intent(inout) :: uvw(1:3) - end subroutine sample_scatter_c - - subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C) - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: i_mat - integer(C_INT), value, intent(in) :: gin - integer(C_INT), intent(inout) :: dg - integer(C_INT), intent(inout) :: gout - end subroutine sample_fission_energy_c - subroutine get_name_c(index, name_len, name) bind(C) use ISO_C_BINDING implicit none diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 56f3392ff2..768ffe9790 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -97,36 +97,6 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3], //============================================================================== -void -sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt, - double uvw[3]) -{ - int gout_c = gout - 1; - macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt); - - // adjust return value for fortran indexing - gout = gout_c + 1; - - // Rotate the angle - rotate_angle_c(uvw, mu, nullptr); -} - -//============================================================================== - -void -sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout) -{ - int dg_c = 0; - int gout_c = 0; - macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c); - - // adjust return values for fortran indexing - dg = dg_c + 1; - gout = gout_c + 1; -} - -//============================================================================== - double get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg) { diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 8651e29f18..10e6878b2c 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -20,15 +20,14 @@ namespace openmc { void -collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size, - const double* energy_bin_avg, const MaterialMacroXS& material_xs) +collision_mg(Particle* p, const double* energy_bin_avg, + const MaterialMacroXS& material_xs) { // Add to the collision counter for the particle p->n_collision++; // Sample the reaction type - sample_reaction(p, fission_bank, fission_bank_size, energy_bin_avg, - material_xs); + sample_reaction(p, energy_bin_avg, material_xs); // Display information about collision if ((settings::verbosity >= 10) || (openmc_trace)) { @@ -39,8 +38,7 @@ collision_mg(Particle* p, Bank* fission_bank, const int64_t fission_bank_size, } void -sample_reaction(Particle* p, Bank* fission_bank, - const int64_t fission_bank_size, const double* energy_bin_avg, +sample_reaction(Particle* p, const double* energy_bin_avg, const MaterialMacroXS& material_xs) { // Create fission bank sites. Note that while a fission reaction is sampled, @@ -50,7 +48,11 @@ sample_reaction(Particle* p, Bank* fission_bank, if (materials[p->material - 1]->fissionable) { if (settings::run_mode == RUN_MODE_EIGENVALUE) { - create_fission_sites(p, fission_bank, n_bank, fission_bank_size, + Bank* result_bank; + int64_t result_bank_size; + // Get pointer to fission bank from Fortran side + openmc_fission_bank(&result_bank, &result_bank_size); + create_fission_sites(p, result_bank, n_bank, result_bank_size, material_xs); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 8d71386243..21e7f9df7a 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -316,8 +316,8 @@ ScattDataLegendre::update_max_val() } // Calculate probability - double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1, - dist[gin][i_gout].data(), mu); + double f = evaluate_legendre(dist[gin][i_gout].size() - 1, + dist[gin][i_gout].data(), mu); // if this is a new maximum, store it if (f > max_val[gin][i_gout]) max_val[gin][i_gout] = f; @@ -339,8 +339,8 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu) f = 0.; } else { int i_gout = gout - gmin[gin]; - f = evaluate_legendre_c(dist[gin][i_gout].size() - 1, - dist[gin][i_gout].data(), mu); + f = evaluate_legendre(dist[gin][i_gout].size() - 1, + dist[gin][i_gout].data(), mu); } return f; } @@ -881,8 +881,8 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, tab.fmu[gin][i_gout].resize(n_mu); for (int imu = 0; imu < n_mu; imu++) { tab.fmu[gin][i_gout][imu] = - evaluate_legendre_c(leg.dist[gin][i_gout].size() - 1, - leg.dist[gin][i_gout].data(), tab.mu[imu]); + evaluate_legendre(leg.dist[gin][i_gout].size() - 1, + leg.dist[gin][i_gout].data(), tab.mu[imu]); } // Ensure positivity diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp index 45f6d3cb86..dabd9de3a8 100644 --- a/src/secondary_nbody.cpp +++ b/src/secondary_nbody.cpp @@ -32,13 +32,13 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_); // x is essentially a Maxwellian distribution - double x = maxwell_spectrum_c(1.0); + double x = maxwell_spectrum(1.0); double y; double r1, r2, r3, r4, r5, r6; switch (n_bodies_) { case 3: - y = maxwell_spectrum_c(1.0); + y = maxwell_spectrum(1.0); break; case 4: r1 = prn(); diff --git a/src/tracking.F90 b/src/tracking.F90 index 02f10f7cc8..982daeccce 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -34,12 +34,9 @@ module tracking implicit none interface - subroutine collision_mg(p, fission_bank, fission_bank_size, & - energy_bin_avg, material_xs) bind(C) - import Particle, Bank, C_INT64_T, C_DOUBLE, MaterialMacroXS + subroutine collision_mg(p, energy_bin_avg, material_xs) bind(C) + import Particle, C_DOUBLE, MaterialMacroXS type(Particle), intent(inout) :: p - type(Bank), intent(inout) :: fission_bank(*) - integer(C_INT64_T), value, intent(in) :: fission_bank_size real(C_DOUBLE), intent(in) :: energy_bin_avg(*) type(MaterialMacroXS), intent(in) :: material_xs end subroutine collision_mg @@ -239,8 +236,7 @@ contains if (run_CE) then call collision(p) else - call collision_mg(p, fission_bank, & - size(fission_bank, kind=C_INT64_T), energy_bin_avg, material_xs) + call collision_mg(p, energy_bin_avg, material_xs) end if ! Score collision estimator tallies -- this is done after a collision From 5af41ce3c350a8219fbf242b9387ff5394bdc310 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Sat, 13 Oct 2018 13:33:06 -0400 Subject: [PATCH 7/7] Resolving review comments: adding doxygen param lists, cleaning up references/pointers --- include/openmc/physics_mg.h | 25 ++++++++++++++---- src/input_xml.F90 | 2 +- src/material.cpp | 2 +- src/material_header.F90 | 6 ++--- src/mgxs_data.F90 | 4 +-- src/physics_mg.cpp | 51 ++++++++++++++++--------------------- src/tracking.F90 | 1 - 7 files changed, 49 insertions(+), 42 deletions(-) diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 09a8cc617d..14b7e6cbd5 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -15,31 +15,46 @@ namespace openmc { // are brought over //! \brief samples particle behavior after a collision event. +//! \param p Particle to operate on +//! \param energy_bin_avg Average energy within each energy bin +//! \param material_xs The cross section cache for the current material extern "C" void collision_mg(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS& material_xs); + const MaterialMacroXS* material_xs); //! \brief samples a reaction type. //! //! Note that there is special logic when suvival biasing is turned on since //! fission and disappearance are treated implicitly. +//! \param p Particle to operate on +//! \param energy_bin_avg Average energy within each energy bin +//! \param material_xs The cross section cache for the current material void sample_reaction(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS& material_xs); + const MaterialMacroXS* material_xs); //! \brief Samples the scattering event +//! \param p Particle to operate on +//! \param energy_bin_avg Average energy within each energy bin void scatter(Particle* p, const double* energy_bin_avg); //! \brief Determines the average total, prompt and delayed neutrons produced //! from fission and creates the appropriate bank sites. +//! \param p Particle to operate on +//! \param bank_array The particle bank to populate +//! \param size_bank Number of particles currently in the bank +//! \param bank_array_size Allocated size of the bank +//! \param material_xs The cross section cache for the current material void -create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, - const int64_t& bank_array_size, const MaterialMacroXS& material_xs); +create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, + int64_t bank_array_size, const MaterialMacroXS* material_xs); //! \brief Handles an absorption event +//! \param p Particle to operate on +//! \param material_xs The cross section cache for the current material void -absorption(Particle* p, const MaterialMacroXS& material_xs); +absorption(Particle* p, const MaterialMacroXS* material_xs); } // namespace openmc #endif // OPENMC_PHYSICS_MG_H diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 12b97d9407..989f7b357f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2986,7 +2986,7 @@ contains ! Check if material is fissionable if (nuclides(materials(i) % nuclide(j)) % fissionable) then - call materials(i) % set_fissionable(logical(.true., C_BOOL)) + call materials(i) % set_fissionable(.true.) end if end do diff --git a/src/material.cpp b/src/material.cpp index 2bd9228753..ecfc7f55e2 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -124,7 +124,7 @@ extern "C" { bool material_fissionable(Material* mat) {return mat->fissionable;} - void material_set_fissionable(Material* mat, bool fissionable, int32_t index) + void material_set_fissionable(Material* mat, bool fissionable) { mat->fissionable = fissionable; } diff --git a/src/material_header.F90 b/src/material_header.F90 index a067269ce2..b4d0617255 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -161,9 +161,9 @@ contains end function material_fissionable subroutine material_set_fissionable(this, fissionable) - class(Material), intent(in) :: this - logical(C_BOOL), intent(in) :: fissionable - call material_set_fissionable_c(this % ptr, fissionable) + class(Material),intent(in) :: this + logical, intent(in) :: fissionable + call material_set_fissionable_c(this % ptr, logical(fissionable, C_BOOL)) end subroutine material_set_fissionable function material_set_density(this, density) result(err) diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index e6eb5a3da5..90178a015d 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -90,8 +90,8 @@ contains end if end do NUCLIDE_LOOP - call mat % set_fissionable(query_fissionable_c(mat % n_nuclides, & - mat % nuclide)) + call mat % set_fissionable( & + logical(query_fissionable_c(mat % n_nuclides, mat % nuclide))) end do MATERIAL_LOOP diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 10e6878b2c..38eb9a5b95 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -21,7 +21,7 @@ namespace openmc { void collision_mg(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS& material_xs) + const MaterialMacroXS* material_xs) { // Add to the collision counter for the particle p->n_collision++; @@ -39,7 +39,7 @@ collision_mg(Particle* p, const double* energy_bin_avg, void sample_reaction(Particle* p, const double* energy_bin_avg, - const MaterialMacroXS& material_xs) + const MaterialMacroXS* material_xs) { // Create fission bank sites. Note that while a fission reaction is sampled, // it never actually "happens", i.e. the weight of the particle does not @@ -52,18 +52,18 @@ sample_reaction(Particle* p, const double* energy_bin_avg, int64_t result_bank_size; // Get pointer to fission bank from Fortran side openmc_fission_bank(&result_bank, &result_bank_size); - create_fission_sites(p, result_bank, n_bank, result_bank_size, + create_fission_sites(p, result_bank, &n_bank, result_bank_size, material_xs); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { - create_fission_sites(p, p->secondary_bank, p->n_secondary, + create_fission_sites(p, p->secondary_bank, &(p->n_secondary), MAX_SECONDARY, material_xs); } } // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (material_xs.absorption > 0.) { + if (material_xs->absorption > 0.) { absorption(p, material_xs); } else { p->absorb_wgt = 0.; @@ -105,24 +105,18 @@ scatter(Particle* p, const double* energy_bin_avg) } void -create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, - const int64_t& bank_array_size, - const MaterialMacroXS& material_xs) +create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, + int64_t bank_array_size, const MaterialMacroXS* material_xs) { // TODO: Heat generation from fission // If uniform fission source weighting is turned on, we increase or decrease // the expected number of fission sites produced - double weight; - if (settings::ufs_on) { - weight = ufs_get_weight(p); - } else { - weight = 1.; - } + double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p->wgt / openmc_keff * weight * material_xs.nu_fission / - material_xs.total; + double nu_t = p->wgt / openmc_keff * weight * material_xs->nu_fission / + material_xs->total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -133,7 +127,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, // Check for the bank size getting hit. For fixed source calculations, this // is a fatal error; for eigenvalue calculations, it just means that k-eff // was too high for a single batch. - if (size_bank + nu > bank_array_size) { + if (*size_bank + nu > bank_array_size) { if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { throw std::runtime_error{"Secondary particle bank size limit reached." " If you are running a subcritical multiplication problem," @@ -151,16 +145,15 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, // Begin banking the source neutrons // First, if our bank is full then don't continue - if ((nu == 0) || (size_bank == bank_array_size)) return; + if ((nu == 0) || (*size_bank == bank_array_size)) return; // Initialize the counter of delayed neutrons encountered for each delayed // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; p->fission = true; - // TODO: +1 for start? - for (size_t i = static_cast(size_bank); - i < static_cast(std::min(size_bank + nu, bank_array_size)); i++) { + for (size_t i = static_cast(*size_bank); + i < static_cast(std::min(*size_bank + nu, bank_array_size)); i++) { // Bank source neutrons by copying the particle data bank_array[i].xyz[0] = p->coord[0].xyz[0]; bank_array[i].xyz[1] = p->coord[0].xyz[1]; @@ -200,7 +193,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, } // Increment number of bank sites - size_bank = std::min(size_bank + nu, bank_array_size); + *size_bank = std::min(*size_bank + nu, bank_array_size); // Store the total weight banked for analog fission tallies p->n_bank = nu; @@ -211,11 +204,11 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t& size_bank, } void -absorption(Particle* p, const MaterialMacroXS& material_xs) +absorption(Particle* p, const MaterialMacroXS* material_xs) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->absorb_wgt = p->wgt * material_xs.absorption / material_xs.total; + p->absorb_wgt = p->wgt * material_xs->absorption / material_xs->total; // Adjust weight of particle by the probability of absorption p->wgt -= p->absorb_wgt; @@ -223,13 +216,13 @@ absorption(Particle* p, const MaterialMacroXS& material_xs) // Score implicit absorpion estimate of keff #pragma omp atomic - global_tally_absorption += p->absorb_wgt * material_xs.nu_fission / - material_xs.absorption; + global_tally_absorption += p->absorb_wgt * material_xs->nu_fission / + material_xs->absorption; } else { - if (material_xs.absorption > prn() * material_xs.total) { + if (material_xs->absorption > prn() * material_xs->total) { #pragma omp atomic - global_tally_absorption += p->wgt * material_xs.nu_fission / - material_xs.absorption; + global_tally_absorption += p->wgt * material_xs->nu_fission / + material_xs->absorption; p->alive = false; p->event = EVENT_ABSORB; } diff --git a/src/tracking.F90 b/src/tracking.F90 index 982daeccce..94c009be00 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -17,7 +17,6 @@ module tracking use nuclide_header use particle_header use physics, only: collision - ! use physics_mg, only: collision_mg use random_lcg, only: prn, prn_set_stream use settings use simulation_header