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, &