crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Added the capability to handle an arbitrary number of legendre moments. This required changes to the input, the output, and the tally handling. The input and math routines are not yet complete, as it still will not allow for any value for the scattering order. I do not expect to extend this past n=5 right now, but I want to have the framework in place so that an extension only requires changing an if-statement in input.xml and adding the additional legendre polynomials to math.F90

Browse source
This commit is contained in:
Adam Nelson 2012-12-03 20:48:14 -05:00
parent 4d8b8b9fb8
commit fa5934cb03
6 changed files with 578 additions and 508 deletions
Download patch file Download diff file Expand all files Collapse all files

20
src/cmfd_input.F90
View file

@ -369,12 +369,16 @@ contains
! allocate scoring bins
allocate(t % score_bins(4))
t % n_score_bins = 4
! allocate scattering order data
allocate(t % scatt_order(4))
t % scatt_order = 0
! set macro_bins
t % score_bins(1) = SCORE_FLUX
t % score_bins(2) = SCORE_TOTAL
t % score_bins(3) = SCORE_SCATTER_1
t % score_bins(4) = SCORE_DIFFUSION
t % score_bins(1) = SCORE_FLUX
t % score_bins(2) = SCORE_TOTAL
t % score_bins(3) = SCORE_SCATTER_N
t % scatt_order(3) = 1
t % score_bins(4) = SCORE_DIFFUSION
! Increment the appropriate index and set pointer
analog_tallies(n_user_analog_tallies + 1) = i
@ -413,6 +417,9 @@ contains
! allocate macro reactions
allocate(t % score_bins(2))
t % n_score_bins = 2
! allocate scattering order data
allocate(t % scatt_order(2))
t % scatt_order = 0
! set macro_bins
t % score_bins(1) = SCORE_NU_SCATTER
@ -454,6 +461,9 @@ contains
! allocate macro reactions
allocate(t % score_bins(1))
t % n_score_bins = 1
! allocate scattering order data
allocate(t % scatt_order(1))
t % scatt_order = 0
! set macro bins
t % score_bins(1) = SCORE_CURRENT

31
src/constants.F90
View file

@ -251,28 +251,25 @@ module constants
EVENT_FISSION = 3
! Tally score type
integer, parameter :: N_SCORE_TYPES = 20
integer, parameter :: N_SCORE_TYPES = 17
integer, parameter :: &
SCORE_FLUX = -1, & ! flux
SCORE_TOTAL = -2, & ! total reaction rate
SCORE_SCATTER = -3, & ! scattering rate
SCORE_NU_SCATTER = -4, & ! scattering production rate
SCORE_SCATTER_1 = -5, & ! first scattering moment
SCORE_SCATTER_2 = -6, & ! second scattering moment
SCORE_SCATTER_3 = -7, & ! third scattering moment
SCORE_SCATTER_4 = -8, & ! fourth scattering moment
SCORE_SCATTER_5 = -9, & ! fifth scattering moment
SCORE_TRANSPORT = -10, & ! transport reaction rate
SCORE_DIFFUSION = -11, & ! diffusion coefficient
SCORE_N_1N = -12, & ! (n,1n) rate
SCORE_N_2N = -13, & ! (n,2n) rate
SCORE_N_3N = -14, & ! (n,3n) rate
SCORE_N_4N = -15, & ! (n,4n) rate
SCORE_ABSORPTION = -16, & ! absorption rate
SCORE_FISSION = -17, & ! fission rate
SCORE_NU_FISSION = -18, & ! neutron production rate
SCORE_CURRENT = -19, & ! partial current
SCORE_EVENTS = -20 ! number of events
SCORE_SCATTER_N = -5, & ! arbitrary scattering moment
SCORE_SCATTER_PN = -6, & ! system for scoring 0th through nth moment
SCORE_TRANSPORT = -7, & ! transport reaction rate
SCORE_DIFFUSION = -8, & ! diffusion coefficient
SCORE_N_1N = -9, & ! (n,1n) rate
SCORE_N_2N = -10, & ! (n,2n) rate
SCORE_N_3N = -11, & ! (n,3n) rate
SCORE_N_4N = -12, & ! (n,4n) rate
SCORE_ABSORPTION = -13, & ! absorption rate
SCORE_FISSION = -14, & ! fission rate
SCORE_NU_FISSION = -15, & ! neutron production rate
SCORE_CURRENT = -16, & ! partial current
SCORE_EVENTS = -17 ! number of events
! Tally map bin finding
integer, parameter :: NO_BIN_FOUND = -1

490
src/input_xml.F90
View file

File diff suppressed because it is too large Load diff

352
src/output.F90
View file

@ -32,35 +32,35 @@ contains
subroutine title()
write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') &
' .d88888b. 888b d888 .d8888b.', &
' d88P" "Y88b 8888b d8888 d88P Y88b', &
' 888 888 88888b.d88888 888 888', &
' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', &
' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', &
' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', &
' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', &
' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', &
'__________________888______________________________________________________', &
' 888', &
' 888'
' .d88888b. 888b d888 .d8888b.', &
' d88P" "Y88b 8888b d8888 d88P Y88b', &
' 888 888 88888b.d88888 888 888', &
' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', &
' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', &
' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', &
' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', &
' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', &
'__________________888______________________________________________________', &
' 888', &
' 888'
! Write version information
write(UNIT=OUTPUT_UNIT, FMT=*) &
' Developed At: Massachusetts Institute of Technology'
' Developed At: Massachusetts Institute of Technology'
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') &
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1
#endif
! Write the date and time
write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') &
time_stamp()
time_stamp()
#ifdef MPI
! Write number of processors
write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') &
trim(to_str(n_procs))
trim(to_str(n_procs))
#endif
end subroutine title
@ -77,7 +77,7 @@ contains
call date_and_time(DATE=date_, TIME=time_)
current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // &
" " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
" " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6)
end function time_stamp
@ -125,13 +125,13 @@ contains
select case (header_level)
case (1)
write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), &
repeat('=', n) // '> ' // trim(line) // ' <' // &
repeat('=', m), repeat('=', 75)
repeat('=', n) // '> ' // trim(line) // ' <' // &
repeat('=', m), repeat('=', 75)
case (2)
write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75)
case (3)
write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // &
trim(line) // ' <' // repeat('=', m)
trim(line) // ' <' // repeat('=', m)
end select
end subroutine header
@ -145,11 +145,11 @@ contains
if (master) then
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') &
"OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
"OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2012 &
&Massachusetts Institute of Technology"
&Massachusetts Institute of Technology"
write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at &
&<http://mit-crpg.github.com/openmc/license.html>"
&<http://mit-crpg.github.com/openmc/license.html>"
end if
end subroutine print_version
@ -207,7 +207,7 @@ contains
else
! Determine last space in current line
i_end = i_start + index(message(i_start+1:i_start+line_wrap), &
' ', BACK=.true.)
' ', BACK=.true.)
! Write up to last space
write(ou, fmt='(1X,A)') message(i_start+1:i_end-1)
@ -270,8 +270,8 @@ contains
l => lattices(coord % lattice)
write(ou,*) ' Lattice = ' // trim(to_str(l % id))
write(ou,*) ' Lattice position = (' // trim(to_str(&
p % coord % lattice_x)) // ',' // trim(to_str(&
p % coord % lattice_y)) // ')'
p % coord % lattice_x)) // ',' // trim(to_str(&
p % coord % lattice_y)) // ')'
end if
! Print local coordinates
@ -390,7 +390,7 @@ contains
case default
index_surf = abs(c % surfaces(i))
string = trim(string) // ' ' // to_str(sign(&
surfaces(index_surf) % id, c % surfaces(i)))
surfaces(index_surf) % id, c % surfaces(i)))
end select
end do
write(unit_,*) ' Surface Specification:' // trim(string)
@ -595,7 +595,7 @@ contains
! Write total atom density in atom/b-cm
write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) &
// ' atom/b-cm'
// ' atom/b-cm'
! Write atom density for each nuclide in material
write(unit_,*) ' Nuclides:'
@ -603,7 +603,7 @@ contains
nuc => nuclides(mat % nuclide(i))
density = mat % atom_density(i)
string = ' ' // trim(nuc % name) // ' = ' // &
trim(to_str(density)) // ' atom/b-cm'
trim(to_str(density)) // ' atom/b-cm'
write(unit_,*) trim(string)
end do
@ -628,7 +628,9 @@ contains
integer :: j ! index in filters array
integer :: id ! user-specified id
integer :: unit_ ! unit to write to
integer :: n ! scattering order to include in name
character(MAX_LINE_LEN) :: string
character(MAX_WORD_LEN) :: pn_string
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Universe), pointer :: u => null()
@ -740,7 +742,7 @@ contains
string = ""
do i = 1, t % filters(j) % n_bins + 1
string = trim(string) // ' ' // trim(to_str(&
t % filters(j) % real_bins(i)))
t % filters(j) % real_bins(i)))
end do
write(unit_,*) ' Incoming Energy Bins:' // trim(string)
end if
@ -751,7 +753,7 @@ contains
string = ""
do i = 1, t % filters(j) % n_bins + 1
string = trim(string) // ' ' // trim(to_str(&
t % filters(j) % real_bins(i)))
t % filters(j) % real_bins(i)))
end do
write(unit_,*) ' Outgoing Energy Bins:' // trim(string)
end if
@ -763,18 +765,20 @@ contains
write(unit_,fmt='(A)',advance='no') ' total'
else
write(unit_,fmt='(A)',advance='no') ' ' // trim(adjustl(&
nuclides(t % nuclide_bins(i)) % name))
nuclides(t % nuclide_bins(i)) % name))
end if
if (mod(i,4) == 0 .and. i /= t % n_nuclide_bins) &
write(unit_,'(/18X)',advance='no')
write(unit_,'(/18X)',advance='no')
end do
write(unit_,*)
! Write score bins
string = ""
do i = 1, t % n_score_bins
select case (t % score_bins(i))
string = ""
j = 0
do i = 1, t % n_nonPN_score_bins
j = j + 1
select case (t % score_bins(j))
case (SCORE_FLUX)
string = trim(string) // ' flux'
case (SCORE_TOTAL)
@ -783,16 +787,17 @@ contains
string = trim(string) // ' scatter'
case (SCORE_NU_SCATTER)
string = trim(string) // ' nu-scatter'
case (SCORE_SCATTER_1)
string = trim(string) // ' scatter-1'
case (SCORE_SCATTER_2)
string = trim(string) // ' scatter-2'
case (SCORE_SCATTER_3)
string = trim(string) // ' scatter-3'
case (SCORE_SCATTER_4)
string = trim(string) // ' scatter-4'
case (SCORE_SCATTER_5)
string = trim(string) // ' scatter-5'
case (SCORE_SCATTER_N)
pn_string = ' scatter-' // trim(to_str(t % scatt_order(j)))
string = trim(string) // pn_string
case (SCORE_SCATTER_PN)
pn_string = ' scatter'
string = trim(string) // pn_string
do n = 1, t % scatt_order(j)
pn_string = ' scatter-' // trim(to_str(n))
string = trim(string) // pn_string
end do
j = j + n - 1
case (SCORE_TRANSPORT)
string = trim(string) // ' transport'
case (SCORE_DIFFUSION)
@ -933,8 +938,8 @@ contains
end if
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
! Accumulate data size
size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8
@ -971,11 +976,11 @@ contains
write(unit_,*) ' Memory Requirements'
write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes'
write(unit_,*) ' Secondary angle distributions = ' // &
trim(to_str(size_angle_total)) // ' bytes'
trim(to_str(size_angle_total)) // ' bytes'
write(unit_,*) ' Secondary energy distributions = ' // &
trim(to_str(size_energy_total)) // ' bytes'
trim(to_str(size_energy_total)) // ' bytes'
write(unit_,*) ' Probability Tables = ' // &
trim(to_str(size_urr)) // ' bytes'
trim(to_str(size_urr)) // ' bytes'
write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes'
! Blank line at end of nuclide
@ -1011,28 +1016,28 @@ contains
! Inelastic data
write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_in))
trim(to_str(sab % n_inelastic_e_in))
write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_e_out))
trim(to_str(sab % n_inelastic_e_out))
write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // &
trim(to_str(sab % n_inelastic_mu))
trim(to_str(sab % n_inelastic_mu))
write(unit_,*) ' Threshold for Inelastic = ' // &
trim(to_str(sab % threshold_inelastic))
trim(to_str(sab % threshold_inelastic))
! Elastic data
if (sab % n_elastic_e_in > 0) then
write(unit_,*) ' # of Incoming Energies (Elastic) = ' // &
trim(to_str(sab % n_elastic_e_in))
trim(to_str(sab % n_elastic_e_in))
write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // &
trim(to_str(sab % n_elastic_mu))
trim(to_str(sab % n_elastic_mu))
write(unit_,*) ' Threshold for Elastic = ' // &
trim(to_str(sab % threshold_elastic))
trim(to_str(sab % threshold_elastic))
end if
! Determine memory used by S(a,b) table and write out
size_sab = 8 * (sab % n_inelastic_e_in * (2 + sab % n_inelastic_e_out * &
(1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * &
(2 + sab % n_elastic_mu))
(1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * &
(2 + sab % n_elastic_mu))
write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes'
! Blank line at end
@ -1060,19 +1065,19 @@ contains
call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1)
write(UNIT=UNIT_SUMMARY, FMT=*) &
"Copyright: 2011-2012 Massachusetts Institute of Technology"
"Copyright: 2011-2012 Massachusetts Institute of Technology"
write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') &
"Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
"Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
#ifdef GIT_SHA1
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1
#endif
write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A)') &
time_stamp()
time_stamp()
! Write information on number of processors
#ifdef MPI
write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') &
trim(to_str(n_procs))
trim(to_str(n_procs))
#endif
! Display problem summary
@ -1117,7 +1122,7 @@ contains
call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY)
write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid))
write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(&
n_grid*n_nuclides_total*4)) // " bytes"
n_grid*n_nuclides_total*4)) // " bytes"
! print summary of variance reduction
call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY)
@ -1227,16 +1232,16 @@ contains
! write out information batch and option independent output
write(UNIT=OUTPUT_UNIT, FMT='(2X,I5)', ADVANCE='NO') current_batch
write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') &
k_batch(current_batch)
k_batch(current_batch)
! write out entropy info
if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
entropy(current_batch)
entropy(current_batch)
! write out accumulated k-effective if after first active batch
if (current_batch > n_inactive + 1) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
keff, keff_std
keff, keff_std
else
write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO')
end if
@ -1244,11 +1249,11 @@ contains
! write out cmfd keff if it is active
if (cmfd_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % keff
cmfd % keff
! write out cmfd entopy
if (cmfd_on .and. entropy_on) write(UNIT=OUTPUT_UNIT, &
FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy
FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy
! next line
write(UNIT=OUTPUT_UNIT, FMT=*)
@ -1275,18 +1280,18 @@ contains
! Write plotting origin
write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
" " // trim(to_str(pl % origin(2))) // " " // &
trim(to_str(pl % origin(3)))
! Write plotting width
if (pl % type == PLOT_TYPE_SLICE) then
write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
" " // trim(to_str(pl % width(2)))
" " // trim(to_str(pl % width(2)))
write(ou,100) "Coloring:", trim(to_str(pl % color_by))
write(ou,100) "Basis:", trim(to_str(pl % basis))
write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
trim(to_str(pl % pixels(2)))
trim(to_str(pl % pixels(2)))
end if
write(ou,*)
@ -1319,7 +1324,7 @@ contains
write(ou,100) " Reading cross sections", time_read_xs % elapsed
write(ou,100) " Unionizing energy grid", time_unionize % elapsed
write(ou,100) "Total time in simulation", time_inactive % elapsed + &
time_active % elapsed
time_active % elapsed
write(ou,100) " Time in transport only", time_transport % elapsed
if(cmfd_run) write(ou,100) "Total CMFD time", time_cmfd % elapsed
write(ou,100) " Time in inactive batches", time_inactive % elapsed
@ -1333,14 +1338,14 @@ contains
if (restart_run) then
total_particles = n_particles * (n_batches - &
restart_batch) * gen_per_batch
restart_batch) * gen_per_batch
else
total_particles = n_particles * n_batches * gen_per_batch
end if
! display calculate rate
speed = real(total_particles) / (time_inactive % elapsed + &
time_active % elapsed)
time_active % elapsed)
string = to_str(speed)
write(ou,101) "Calculation Rate", trim(string)
@ -1358,13 +1363,13 @@ contains
! write global tallies
write(ou,102) "k-effective (Analog)", global_tallies(K_ANALOG) % sum, &
global_tallies(K_ANALOG) % sum_sq
global_tallies(K_ANALOG) % sum_sq
write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) % sum, &
global_tallies(K_COLLISION) % sum_sq
global_tallies(K_COLLISION) % sum_sq
write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum, &
global_tallies(K_TRACKLENGTH) % sum_sq
global_tallies(K_TRACKLENGTH) % sum_sq
write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, &
global_tallies(LEAKAGE) % sum_sq
global_tallies(LEAKAGE) % sum_sq
write(ou,*)
! format for write statements
@ -1385,17 +1390,21 @@ contains
integer :: j ! level in tally hierarchy
integer :: k ! loop index for scoring bins
integer :: n ! loop index for nuclides
integer :: l ! loop index for nonPN scores
integer :: type ! type of tally filter
integer :: indent ! number of spaces to preceed output
integer :: filter_index ! index in scores array for filters
integer :: score_index ! scoring bin index
integer :: i_nuclide ! index in nuclides array
integer :: i_listing ! index in xs_listings array
integer :: nOrder ! loop index for scattering orders
real(8) :: t_value ! t-values for confidence intervals
real(8) :: alpha ! significance level for CI
character(MAX_FILE_LEN) :: filename ! name of output file
character(15) :: filter_name(N_FILTER_TYPES) ! names of tally filters
character(27) :: score_name(N_SCORE_TYPES) ! names of scoring function
character(27) :: score_names(N_SCORE_TYPES) ! names of scoring function
character(27) :: score_name ! names of scoring function
! to be applied at write-time
type(TallyObject), pointer :: t
! Skip if there are no tallies
@ -1412,25 +1421,22 @@ contains
filter_name(FILTER_ENERGYOUT) = "Outgoing Energy"
! Initialize names for scores
score_name(abs(SCORE_FLUX)) = "Flux"
score_name(abs(SCORE_TOTAL)) = "Total Reaction Rate"
score_name(abs(SCORE_SCATTER)) = "Scattering Rate"
score_name(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate"
score_name(abs(SCORE_SCATTER_1)) = "First Scattering Moment"
score_name(abs(SCORE_SCATTER_2)) = "Second Scattering Moment"
score_name(abs(SCORE_SCATTER_3)) = "Third Scattering Moment"
score_name(abs(SCORE_SCATTER_4)) = "Fourth Scattering Moment"
score_name(abs(SCORE_SCATTER_5)) = "Fifth Scattering Moment"
score_name(abs(SCORE_TRANSPORT)) = "Transport Rate"
score_name(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient"
score_name(abs(SCORE_N_1N)) = "(n,1n) Rate"
score_name(abs(SCORE_N_2N)) = "(n,2n) Rate"
score_name(abs(SCORE_N_3N)) = "(n,3n) Rate"
score_name(abs(SCORE_N_4N)) = "(n,4n) Rate"
score_name(abs(SCORE_ABSORPTION)) = "Absorption Rate"
score_name(abs(SCORE_FISSION)) = "Fission Rate"
score_name(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate"
score_name(abs(SCORE_EVENTS)) = "Events"
score_names(abs(SCORE_FLUX)) = "Flux"
score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate"
score_names(abs(SCORE_SCATTER)) = "Scattering Rate"
score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate"
score_names(abs(SCORE_SCATTER_N)) = ""
score_names(abs(SCORE_SCATTER_PN)) = ""
score_names(abs(SCORE_TRANSPORT)) = "Transport Rate"
score_names(abs(SCORE_DIFFUSION)) = "Diffusion Coefficient"
score_names(abs(SCORE_N_1N)) = "(n,1n) Rate"
score_names(abs(SCORE_N_2N)) = "(n,2n) Rate"
score_names(abs(SCORE_N_3N)) = "(n,3n) Rate"
score_names(abs(SCORE_N_4N)) = "(n,4n) Rate"
score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate"
score_names(abs(SCORE_FISSION)) = "Fission Rate"
score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate"
score_names(abs(SCORE_EVENTS)) = "Events"
! Create filename for tally output
if (run_mode == MODE_TALLIES) then
@ -1468,10 +1474,10 @@ contains
! Write header block
if (t % label == "") then
call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, &
level=3)
level=3)
else
call header("TALLY " // trim(to_str(t % id)) // ": " &
// trim(t % label), unit=UNIT_TALLY, level=3)
// trim(t % label), unit=UNIT_TALLY, level=3)
endif
! Handle surface current tallies separately
@ -1520,7 +1526,7 @@ contains
! Print current filter information
type = t % filters(j) % type
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
trim(filter_name(type)), trim(get_label(t, j))
indent = indent + 2
j = j + 1
end if
@ -1531,7 +1537,7 @@ contains
if (t % n_filters > 0) then
type = t % filters(j) % type
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(filter_name(type)), trim(get_label(t, j))
trim(filter_name(type)), trim(get_label(t, j))
end if
! Determine scoring index for this bin combination -- note that unlike
@ -1543,7 +1549,7 @@ contains
else
filter_index = 1
end if
! Write scores for this filter bin combination
score_index = 0
if (t % n_filters > 0) indent = indent + 2
@ -1552,20 +1558,52 @@ contains
i_nuclide = t % nuclide_bins(n)
if (i_nuclide == -1) then
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
"Total Material"
"Total Material"
else
i_listing = nuclides(i_nuclide) % listing
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(xs_listings(i_listing) % alias)
trim(xs_listings(i_listing) % alias)
end if
indent = indent + 2
do k = 1, t % n_score_bins
k = 0
do l = 1, t % n_nonPN_score_bins
k = k + 1
score_index = score_index + 1
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name(abs(t % score_bins(k))), &
to_str(t % scores(score_index,filter_index) % sum), &
trim(to_str(t % scores(score_index,filter_index) % sum_sq))
if (t % score_bins(k) == SCORE_SCATTER_N) then
if (t % scatt_order(k) == 0) then
score_name = "Scattering Rate"
else
score_name = 'Order-' // trim(to_str(t % scatt_order(k))) // &
' Scattering Moment'
end if
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % scores(score_index,filter_index) % sum), &
trim(to_str(t % scores(score_index,filter_index) % sum_sq))
else if (t % score_bins(k) == SCORE_SCATTER_PN) then
score_name = "Scattering Rate"
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % scores(score_index,filter_index) % sum), &
trim(to_str(t % scores(score_index,filter_index) % sum_sq))
do nOrder = 1, t % scatt_order(k)
score_index = score_index + 1
score_name = 'Order-' // trim(to_str(nOrder)) // &
' Scattering Moment'
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % scores(score_index,filter_index) % sum), &
trim(to_str(t % scores(score_index,filter_index) % sum_sq))
end do
k = k + nOrder - 1
else
score_name = score_names(abs(t % score_bins(k)))
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % scores(score_index,filter_index) % sum), &
trim(to_str(t % scores(score_index,filter_index) % sum_sq))
end if
end do
indent = indent - 2
@ -1642,110 +1680,110 @@ contains
! Write incoming energy bin
write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') &
"Incoming Energy", trim(get_label(t, i_filter_ein))
"Incoming Energy", trim(get_label(t, i_filter_ein))
end if
! Left Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Left", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Left", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Left", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Left", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
! Right Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Right", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Right", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Right", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Right", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
! Back Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Back", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Back", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Back", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Back", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
! Front Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Front", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Front", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_FRONT
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Front", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Front", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
! Bottom Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Bottom", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Bottom", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Bottom", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Bottom", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
! Top Surface
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
t % matching_bins(i_filter_surf) = IN_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Incoming Current from Top", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Incoming Current from Top", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
t % matching_bins(i_filter_surf) = OUT_TOP
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') &
"Outgoing Current to Top", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
"Outgoing Current to Top", &
to_str(t % scores(1,filter_index) % sum), &
trim(to_str(t % scores(1,filter_index) % sum_sq))
end do
end do
@ -1793,10 +1831,10 @@ contains
call bin_to_mesh_indices(m, bin, ijk)
if (m % n_dimension == 2) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ")"
trim(to_str(ijk(2))) // ")"
elseif (m % n_dimension == 3) then
label = "Index (" // trim(to_str(ijk(1))) // ", " // &
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
end if
case (FILTER_ENERGYIN, FILTER_ENERGYOUT)
E0 = t % filters(i_filter) % real_bins(bin)

188
src/tally.F90
View file

@ -6,7 +6,7 @@ module tally
use global
use math, only: t_percentile, calc_pn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, &
mesh_indices_to_bin, mesh_intersects
mesh_indices_to_bin, mesh_intersects
use mesh_header, only: StructuredMesh
use output, only: header
use search, only: binary_search
@ -34,6 +34,9 @@ contains
integer :: j ! loop index for scoring bins
integer :: k ! loop index for nuclide bins
integer :: n ! loop index for scattering order
integer :: l ! scoring bin loop index, allowing for changing
! position during the loop
integer :: filter_index ! single index for single bin
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
integer :: i_nuclide ! index in nuclides array
@ -119,7 +122,9 @@ contains
end if
! Determine score for each bin
SCORE_LOOP: do j = 1, t % n_score_bins
j = 0
SCORE_LOOP: do l = 1, t % n_nonPN_score_bins
j = j + 1
! determine what type of score bin
score_bin = t % score_bins(j)
@ -160,52 +165,36 @@ contains
! reaction with neutrons in the exit channel
score = wgt
case (SCORE_SCATTER_1)
case (SCORE_SCATTER_N)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
! The first scattering moment can be determined by using the
! rate of scattering reactions multiplied by the cosine of the
! change in neutron's angle due to the collision
! Find the scattering order for a singly requested moment, and
! store its moment contribution.
score = last_wgt * mu
if (t % scatt_order(j) == 1) then
score = last_wgt * mu ! avoid function call overhead
else
score = last_wgt * calc_pn(t % scatt_order(j), mu)
endif
case (SCORE_SCATTER_2)
case (SCORE_SCATTER_PN)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
score_index = score_index - 1
! Find the scattering order for a collection of requested moments
! and store the moment contribution of each
do n = 0, t % scatt_order(j)
! determine scoring bin index
score_index = score_index + 1
! get the score and tally it
score = last_wgt * calc_pn(n, mu)
if (n < t % scatt_order(j)) &
call add_to_score(t % scores(score_index, filter_index), score)
end do
j = j + t % scatt_order(j)
! The second scattering moment can be determined in a similar
! manner to the first scattering moment
score = last_wgt * calc_pn(2, mu)
case (SCORE_SCATTER_3)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
! The third scattering moment can be determined in a similar
! manner to the first scattering moment
score = last_wgt * calc_pn(3, mu)
case (SCORE_SCATTER_4)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
! The fourth scattering moment can be determined in a similar
! manner to the first scattering moment
score = last_wgt * calc_pn(4, mu)
case (SCORE_SCATTER_5)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
! The fifth scattering moment can be determined in a similar
! manner to the first scattering moment
score = last_wgt * calc_pn(5, mu)
case (SCORE_TRANSPORT)
! Skip any event where the particle didn't scatter
if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP
@ -233,7 +222,7 @@ contains
! mu*Sigma_s)).
score = last_wgt / (3.0_8 * score * (material_xs % total - &
mu * score))
mu * score))
case (SCORE_N_1N)
! Skip any event where the particle didn't scatter
@ -326,6 +315,7 @@ contains
case default
message = "Invalid score type on tally " // to_str(t % id) // "."
write(*,*) score_index, score_bin, l,j, t % n_score_bins
call fatal_error()
end select
@ -393,7 +383,7 @@ contains
! change outgoing energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
size(t % filters(i) % real_bins), E_out)
size(t % filters(i) % real_bins), E_out)
! determine scoring index
i_filter = sum((t % matching_bins - 1) * t % stride) + 1
@ -474,7 +464,7 @@ contains
if (t % all_nuclides) then
call score_all_nuclides(tracklength_tallies(curr_ptr % data), flux, &
filter_index)
filter_index)
else
NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins
@ -511,20 +501,20 @@ contains
select case(score_bin)
case (SCORE_TOTAL)
score = micro_xs(i_nuclide) % total * &
atom_density * flux
atom_density * flux
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
case (SCORE_ABSORPTION)
score = micro_xs(i_nuclide) % absorption * &
atom_density * flux
atom_density * flux
case (SCORE_FISSION)
score = micro_xs(i_nuclide) % fission * &
atom_density * flux
atom_density * flux
case (SCORE_NU_FISSION)
score = micro_xs(i_nuclide) % nu_fission * &
atom_density * flux
atom_density * flux
case (SCORE_EVENTS)
score = ONE
case default
@ -634,7 +624,7 @@ contains
score = micro_xs(i_nuclide) % total * atom_density * flux
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * atom_density * flux
micro_xs(i_nuclide) % absorption) * atom_density * flux
case (SCORE_ABSORPTION)
score = micro_xs(i_nuclide) % absorption * atom_density * flux
case (SCORE_FISSION)
@ -777,27 +767,27 @@ contains
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
p % material, i_tally)
case (FILTER_CELL)
! determine next cell bin
! TODO: Account for cells in multiple levels when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
p % coord % cell, i_tally)
p % coord % cell, i_tally)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
p % surface, i_tally)
case (FILTER_ENERGYIN)
! determine incoming energy bin
@ -805,12 +795,12 @@ contains
! check if energy of the particle is within energy bins
if (p % E < t % filters(i) % real_bins(1) .or. &
p % E > t % filters(i) % real_bins(k + 1)) then
p % E > t % filters(i) % real_bins(k + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
k + 1, p % E)
k + 1, p % E)
end if
end select
@ -931,25 +921,25 @@ contains
select case(score_bin)
case (SCORE_TOTAL)
score = micro_xs(i_nuclide) % total * &
atom_density * flux
atom_density * flux
case (SCORE_SCATTER)
score = (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
micro_xs(i_nuclide) % absorption) * &
atom_density * flux
case (SCORE_ABSORPTION)
score = micro_xs(i_nuclide) % absorption * &
atom_density * flux
atom_density * flux
case (SCORE_FISSION)
score = micro_xs(i_nuclide) % fission * &
atom_density * flux
atom_density * flux
case (SCORE_NU_FISSION)
score = micro_xs(i_nuclide) % nu_fission * &
atom_density * flux
atom_density * flux
case (SCORE_EVENTS)
score = ONE
case default
message = "Invalid score type on tally " // &
to_str(t % id) // "."
to_str(t % id) // "."
call fatal_error()
end select
@ -972,7 +962,7 @@ contains
score = ONE
case default
message = "Invalid score type on tally " // &
to_str(t % id) // "."
to_str(t % id) // "."
call fatal_error()
end select
end if
@ -1030,27 +1020,27 @@ contains
! determine next universe bin
! TODO: Account for multiple universes when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, &
p % coord % universe, i_tally)
p % coord % universe, i_tally)
case (FILTER_MATERIAL)
t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, &
p % material, i_tally)
p % material, i_tally)
case (FILTER_CELL)
! determine next cell bin
! TODO: Account for cells in multiple levels when performing this filter
t % matching_bins(i) = get_next_bin(FILTER_CELL, &
p % coord % cell, i_tally)
p % coord % cell, i_tally)
case (FILTER_CELLBORN)
! determine next cellborn bin
t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, &
p % cell_born, i_tally)
p % cell_born, i_tally)
case (FILTER_SURFACE)
! determine next surface bin
t % matching_bins(i) = get_next_bin(FILTER_SURFACE, &
p % surface, i_tally)
p % surface, i_tally)
case (FILTER_ENERGYIN)
! determine incoming energy bin
@ -1065,12 +1055,12 @@ contains
! check if energy of the particle is within energy bins
if (E < t % filters(i) % real_bins(1) .or. &
E > t % filters(i) % real_bins(n + 1)) then
E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, E)
n + 1, E)
end if
case (FILTER_ENERGYOUT)
@ -1079,12 +1069,12 @@ contains
! check if energy of the particle is within energy bins
if (E < t % filters(i) % real_bins(1) .or. &
E > t % filters(i) % real_bins(n + 1)) then
E > t % filters(i) % real_bins(n + 1)) then
t % matching_bins(i) = NO_BIN_FOUND
else
! search to find incoming energy bin
t % matching_bins(i) = binary_search(t % filters(i) % real_bins, &
n + 1, p % E)
n + 1, p % E)
end if
end select
@ -1175,14 +1165,14 @@ contains
n = t % filters(j) % n_bins
! check if energy of the particle is within energy bins
if (p % E < t % filters(j) % real_bins(1) .or. &
p % E > t % filters(j) % real_bins(n + 1)) then
p % E > t % filters(j) % real_bins(n + 1)) then
curr_ptr => curr_ptr % next ! select next tally
cycle
end if
! search to find incoming energy bin
t % matching_bins(j) = binary_search(t % filters(j) % real_bins, &
n + 1, p % E)
n + 1, p % E)
end if
! =======================================================================
@ -1200,7 +1190,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
end if
@ -1211,7 +1201,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
end if
@ -1227,7 +1217,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
end if
@ -1238,7 +1228,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
end if
@ -1254,7 +1244,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
end if
@ -1265,7 +1255,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
filter_index = sum((t % matching_bins - 1) * t % stride) + 1
call add_to_score(t % scores(1, filter_index), p % wgt)
end if
@ -1318,7 +1308,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(1) = ijk0(1) + 1
xyz_cross(1) = xyz_cross(1) + m % width(1)
@ -1330,7 +1320,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_RIGHT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
elseif (distance == d(2)) then
@ -1340,7 +1330,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(2) = ijk0(2) + 1
xyz_cross(2) = xyz_cross(2) + m % width(2)
@ -1352,7 +1342,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_FRONT
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
else if (distance == d(3)) then
@ -1362,7 +1352,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = OUT_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
ijk0(3) = ijk0(3) + 1
xyz_cross(3) = xyz_cross(3) + m % width(3)
@ -1374,7 +1364,7 @@ contains
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
t % matching_bins(i_filter_surf) = IN_TOP
t % matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, ijk0 + 1, .true.)
mesh_indices_to_bin(m, ijk0 + 1, .true.)
end if
end if
end if
@ -1385,7 +1375,7 @@ contains
! Check for errors
if (filter_index <= 0 .or. filter_index > &
t % total_filter_bins) then
t % total_filter_bins) then
message = "Score index outside range."
call fatal_error()
end if
@ -1442,7 +1432,7 @@ contains
end if
i_tally_check = tally_maps(filter_type) % items(filter_value) % &
elements(position(filter_type)) % index_tally
elements(position(filter_type)) % index_tally
if (i_tally_check > i_tally) then
! Since the index being checked against is greater than the index we
@ -1454,7 +1444,7 @@ contains
elseif (i_tally_check == i_tally) then
! Found a match
bin = tally_maps(filter_type) % items(filter_value) % &
elements(position(filter_type)) % index_bin
elements(position(filter_type)) % index_bin
return
end if
@ -1539,14 +1529,14 @@ contains
! The MPI_IN_PLACE specifier allows the master to copy values into
! a receive buffer without having a temporary variable
call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Transfer values to value on master
t % scores(:,:) % value = tally_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, &
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Reset value on other processors
t % scores(:,:) % value = 0
@ -1562,14 +1552,14 @@ contains
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Transfer values back to global_tallies on master
global_tallies(:) % value = global_temp
else
! Receive buffer not significant at other processors
call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, &
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err)
! Reset value on other processors
global_tallies(:) % value = ZERO
@ -1580,11 +1570,11 @@ contains
! last realization
if (master) then
call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, MPI_COMM_WORLD, mpi_err)
else
! Receive buffer not significant at other processors
call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, &
0, MPI_COMM_WORLD, mpi_err)
0, MPI_COMM_WORLD, mpi_err)
end if
if (associated(curr_ptr)) nullify(curr_ptr)
@ -1690,7 +1680,7 @@ contains
score % sum = score % sum/n
score % sum_sq = sqrt((score % sum_sq/n - score % sum * &
score % sum) / (n - 1))
score % sum) / (n - 1))
end subroutine statistics_score
@ -1911,7 +1901,7 @@ contains
end if
do i = n_cmfd_current_tallies + n_user_current_tallies, &
n_user_current_tallies + 1, -1
n_user_current_tallies + 1, -1
! allocate node
allocate(curr_ptr)

5
src/tally_header.F90
View file

@ -97,8 +97,13 @@ module tally_header
logical :: all_nuclides = .false.
! Values to score, e.g. flux, absorption, etc.
! scat_order is the scattering order for each score.
! It is to be 0 if the scattering order is 0, or if the score is not a
! scattering response.
integer :: n_score_bins = 0
integer, allocatable :: score_bins(:)
integer, allocatable :: scatt_order(:)
integer :: n_nonPN_score_bins = 0
! Scores for each bin -- the first dimenion of the array is for scores
! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the