mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Added create_fission_sites and absorption to physics_mg.cpp
This commit is contained in:
Adam G Nelson
parent
960f71e276
commit
ef9908e774
5 changed files with 183 additions and 164 deletions
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,12 +1,23 @@
|
|||
#include "openmc/physics_mg.h"
|
||||
|
||||
#include <stdexcept>
|
||||
#include <sstream>
|
||||
|
||||
#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<int>(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_t>(size_bank);
|
||||
i < static_cast<size_t>(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<int>(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<double>(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
|
||||
|
|
|
|||
|
|
@ -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, &
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue