diff --git a/src/output.f90 b/src/output.f90 index 7b23d636ee..554c369943 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -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 diff --git a/src/particle_header.f90 b/src/particle_header.f90 index d8af8e0dff..ba1da0d403 100644 --- a/src/particle_header.f90 +++ b/src/particle_header.f90 @@ -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 diff --git a/src/physics.f90 b/src/physics.f90 index fae0bfbc01..0b08b2cd0a 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -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) diff --git a/src/tally.f90 b/src/tally.f90 index 2310550469..db1f381b40 100644 --- a/src/tally.f90 +++ b/src/tally.f90 @@ -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