mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
First shot at tallies with outgoing energy filter (Answers not correct). Improved tally output and changed header subroutine to take an optional argument specifying unit.
This commit is contained in:
Paul Romano
parent
cb90da9408
commit
c9b2967c64
4 changed files with 95 additions and 37 deletions
|
|
@ -92,13 +92,15 @@ contains
|
|||
! specified, it is assumed to be a minor header block (H3).
|
||||
!===============================================================================
|
||||
|
||||
subroutine header(msg, level)
|
||||
subroutine header(msg, level, unit_file)
|
||||
|
||||
character(*), intent(in) :: msg
|
||||
integer, optional :: level
|
||||
integer, optional :: unit_file
|
||||
|
||||
integer :: header_level
|
||||
integer :: n, m
|
||||
integer :: unit
|
||||
character(75) :: line
|
||||
|
||||
! set default header level
|
||||
|
|
@ -108,8 +110,15 @@ contains
|
|||
header_level = level
|
||||
end if
|
||||
|
||||
! set default unit
|
||||
if (present(unit_file)) then
|
||||
unit = unit_file
|
||||
else
|
||||
unit = ou
|
||||
end if
|
||||
|
||||
! Print first blank line
|
||||
write(ou,*)
|
||||
write(unit,*)
|
||||
|
||||
! determine how many times to repeat '=' character
|
||||
n = (63 - len_trim(msg))/2
|
||||
|
|
@ -124,21 +133,21 @@ contains
|
|||
select case (header_level)
|
||||
case (1)
|
||||
! determine number of spaces to put in from of header
|
||||
write(ou,*) repeat('=', 75)
|
||||
write(ou,*) repeat('=', n) // '> ' // trim(line) // ' <' // &
|
||||
write(unit,*) repeat('=', 75)
|
||||
write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // &
|
||||
& repeat('=', m)
|
||||
write(ou,*) repeat('=', 75)
|
||||
write(unit,*) repeat('=', 75)
|
||||
case (2)
|
||||
write(ou,*) trim(line)
|
||||
write(ou,*) repeat('-', 75)
|
||||
write(unit,*) trim(line)
|
||||
write(unit,*) repeat('-', 75)
|
||||
case (3)
|
||||
n = (63 - len_trim(line))/2
|
||||
write(ou,*) repeat('=', n) // '> ' // trim(line) // ' <' // &
|
||||
write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // &
|
||||
& repeat('=', m)
|
||||
end select
|
||||
|
||||
! Print trailing blank line
|
||||
write(ou, *)
|
||||
write(unit, *)
|
||||
|
||||
end subroutine header
|
||||
|
||||
|
|
|
|||
|
|
@ -20,6 +20,10 @@ module particle_header
|
|||
real(8) :: E ! energy
|
||||
logical :: alive ! is particle alive?
|
||||
|
||||
! Pre-collision physical data
|
||||
real(8) :: last_wgt ! last particle weight
|
||||
real(8) :: last_E ! last energy
|
||||
|
||||
! Energy grid data
|
||||
integer :: IE ! index on energy grid
|
||||
real(8) :: interp ! interpolation factor for energy grid
|
||||
|
|
|
|||
|
|
@ -283,23 +283,25 @@ contains
|
|||
real(8) :: prob ! cumulative probability
|
||||
real(8) :: cutoff ! random number
|
||||
real(8) :: atom_density ! atom density of nuclide in atom/b-cm
|
||||
logical :: scatter ! was this a scattering reaction?
|
||||
character(MAX_LINE_LEN) :: msg ! output/error message
|
||||
type(Material), pointer :: mat
|
||||
type(Nuclide), pointer :: nuc
|
||||
type(Reaction), pointer :: rxn
|
||||
|
||||
! Set scatter to false by default
|
||||
scatter = .false.
|
||||
|
||||
! Store pre-collision particle properties
|
||||
p % last_wgt = p % wgt
|
||||
p % last_E = p % E
|
||||
|
||||
! Add to collision counter for particle
|
||||
p % n_collision = p % n_collision + 1
|
||||
|
||||
! Get pointer to current material
|
||||
mat => materials(p % material)
|
||||
|
||||
! Score collision estimator tallies for any macro tallies -- we can do this
|
||||
! before sampling the nuclide since
|
||||
if (tallies_on) then
|
||||
call score_tally(p)
|
||||
end if
|
||||
|
||||
! sample nuclide
|
||||
i = 0
|
||||
total = material_xs % total
|
||||
|
|
@ -353,10 +355,12 @@ contains
|
|||
select case (rxn % MT)
|
||||
case (ELASTIC)
|
||||
call elastic_scatter(p, nuc, rxn)
|
||||
scatter = .true.
|
||||
case (N_NA, N_N3A, N_NP, N_N2A, N_ND, N_NT, N_N3HE, N_NT2A, N_N2P, &
|
||||
N_NPA, N_N1 : N_NC, N_2ND, N_2N, N_3N, N_2NA, N_3NA, N_2N2A, &
|
||||
N_4N, N_2NP, N_3NP)
|
||||
call inelastic_scatter(p, nuc, rxn)
|
||||
scatter = .true.
|
||||
case (N_FISSION, N_F, N_NF, N_2NF, N_3NF)
|
||||
call n_fission_event(p, nuc, rxn)
|
||||
case (N_GAMMA : N_DA)
|
||||
|
|
@ -373,6 +377,12 @@ contains
|
|||
! call warning(msg)
|
||||
end if
|
||||
|
||||
! Score collision estimator tallies for any macro tallies -- we can do this
|
||||
! before sampling the nuclide since
|
||||
if (tallies_on) then
|
||||
call score_tally(p, scatter)
|
||||
end if
|
||||
|
||||
! find energy index, interpolation factor
|
||||
call find_energy_index(p)
|
||||
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@ module tally
|
|||
use error, only: fatal_error
|
||||
use global
|
||||
use mesh, only: get_mesh_bin, bin_to_mesh_indices
|
||||
use output, only: message
|
||||
use output, only: message, header
|
||||
use search, only: binary_search
|
||||
use string, only: int_to_str, real_to_str
|
||||
use tally_header, only: TallyScore, TallyMapItem, TallyMapElement
|
||||
|
|
@ -258,9 +258,10 @@ contains
|
|||
! SCORE_TALLY contains the main logic for scoring user-specified tallies
|
||||
!===============================================================================
|
||||
|
||||
subroutine score_tally(p)
|
||||
subroutine score_tally(p, scatter)
|
||||
|
||||
type(Particle), pointer :: p ! particle
|
||||
type(Particle), pointer :: p
|
||||
logical, intent(in) :: scatter
|
||||
|
||||
integer :: i
|
||||
integer :: j
|
||||
|
|
@ -268,10 +269,17 @@ contains
|
|||
integer :: bins(TALLY_TYPES)
|
||||
integer :: score_index
|
||||
real(8) :: score
|
||||
real(8) :: last_wgt
|
||||
real(8) :: wgt
|
||||
logical :: in_mesh
|
||||
logical :: has_outgoing_energy
|
||||
type(TallyObject), pointer :: t
|
||||
type(StructuredMesh), pointer :: m
|
||||
|
||||
! Copy particle's pre- and post-collision weight
|
||||
last_wgt = p % last_wgt
|
||||
wgt = p % wgt
|
||||
|
||||
! A loop over all tallies is necessary because we need to simultaneously
|
||||
! determine different filter bins for the same tally in order to score to it
|
||||
|
||||
|
|
@ -334,10 +342,11 @@ contains
|
|||
n = t % n_bins(T_ENERGYIN)
|
||||
if (n > 0) then
|
||||
! check if energy of the particle is within energy bins
|
||||
if (p % E < t % energy_in(1) .or. p % E > t % energy_in(n + 1)) cycle
|
||||
if (p % last_E < t % energy_in(1) .or. &
|
||||
p % last_E > t % energy_in(n + 1)) cycle
|
||||
|
||||
! search to find incoming energy bin
|
||||
bins(T_ENERGYIN) = binary_search(t % energy_in, n + 1, p % E)
|
||||
bins(T_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E)
|
||||
else
|
||||
bins(T_ENERGYIN) = 1
|
||||
end if
|
||||
|
|
@ -350,8 +359,10 @@ contains
|
|||
|
||||
! search to find incoming energy bin
|
||||
bins(T_ENERGYOUT) = binary_search(t % energy_out, n + 1, p % E)
|
||||
has_outgoing_energy = .true.
|
||||
else
|
||||
bins(T_ENERGYOUT) = 1
|
||||
has_outgoing_energy = .false.
|
||||
end if
|
||||
|
||||
! =======================================================================
|
||||
|
|
@ -364,24 +375,44 @@ contains
|
|||
! Determine scoring index for this filter combination
|
||||
score_index = sum((bins - 1) * t % stride) + 1
|
||||
|
||||
! Determine score for each bin
|
||||
do j = 1, t % n_macro_bins
|
||||
! Determine score
|
||||
select case(t % macro_bins(j) % scalar)
|
||||
case (MACRO_FLUX)
|
||||
score = p % wgt / material_xs % total
|
||||
case (MACRO_TOTAL)
|
||||
score = p % wgt
|
||||
case (MACRO_SCATTER)
|
||||
score = p % wgt * (material_xs % total - material_xs % absorption) &
|
||||
/ material_xs % total
|
||||
case (MACRO_ABSORPTION)
|
||||
score = p % wgt * material_xs % absorption / material_xs % total
|
||||
case (MACRO_FISSION)
|
||||
score = p % wgt * material_xs % fission / material_xs % total
|
||||
case (MACRO_NU_FISSION)
|
||||
score = p % wgt * material_xs % nu_fission / material_xs % total
|
||||
end select
|
||||
if (has_outgoing_energy) then
|
||||
! If this tally has an outgoing energy filter, the only supported
|
||||
! reaction is scattering. For all other reactions, about
|
||||
|
||||
if (.not. scatter) return
|
||||
|
||||
! Make sure bin is scattering -- since scattering has already
|
||||
! occured, we do not need to multiply by the scattering cross
|
||||
! section
|
||||
|
||||
if (t % macro_bins(j) % scalar == MACRO_SCATTER) then
|
||||
score = last_wgt
|
||||
else
|
||||
! call fatal_error
|
||||
end if
|
||||
else
|
||||
! For tallies with no outgoing energy filter, the score is
|
||||
! calculated normally depending on the quantity specified
|
||||
|
||||
select case(t % macro_bins(j) % scalar)
|
||||
case (MACRO_FLUX)
|
||||
score = last_wgt / material_xs % total
|
||||
case (MACRO_TOTAL)
|
||||
score = last_wgt
|
||||
case (MACRO_SCATTER)
|
||||
score = last_wgt * (material_xs % total - material_xs % absorption) &
|
||||
/ material_xs % total
|
||||
case (MACRO_ABSORPTION)
|
||||
score = last_wgt * material_xs % absorption / material_xs % total
|
||||
case (MACRO_FISSION)
|
||||
score = last_wgt * material_xs % fission / material_xs % total
|
||||
case (MACRO_NU_FISSION)
|
||||
score = last_wgt * material_xs % nu_fission / material_xs % total
|
||||
end select
|
||||
end if
|
||||
|
||||
! Add score to tally
|
||||
call add_to_score(t % scores(score_index, j), score)
|
||||
|
||||
|
|
@ -561,6 +592,9 @@ contains
|
|||
do i = 1, n_tallies
|
||||
t => tallies(i)
|
||||
|
||||
! Write header block
|
||||
call header("TALLY " // trim(int_to_str(t % uid)), 3, UNIT_TALLY)
|
||||
|
||||
! First determine which filters this tally has
|
||||
do j = 1, TALLY_TYPES
|
||||
if (t % n_bins(j) > 0) then
|
||||
|
|
@ -577,7 +611,8 @@ contains
|
|||
! loop for each filter variable (given that only a few filters are likely
|
||||
! to be used for a given tally.
|
||||
|
||||
! Initialize filter level and indentation
|
||||
! Initialize bins, filter level, and indentation
|
||||
bins = 0
|
||||
j = 1
|
||||
indent = 0
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue