From f2674c28abf40b235fcb8a27e2dde9ee5d981cb8 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 2 Oct 2016 11:58:39 -0400 Subject: [PATCH] Move threadprivate flux derivatives to Particle --- src/global.F90 | 1 - src/input_xml.F90 | 12 +--- src/particle_header.F90 | 3 + src/simulation.F90 | 5 +- src/tally.F90 | 154 +++++++++++++++++++--------------------- src/tally_header.F90 | 1 - src/tracking.F90 | 4 +- 7 files changed, 84 insertions(+), 96 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index 46508ab5b7..7ff8a26b62 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -181,7 +181,6 @@ module global ! Tally derivatives type(TallyDerivative), allocatable :: tally_derivs(:) -!$omp threadprivate(tally_derivs) ! Normalization for statistics integer :: n_realizations = 0 ! # of independent realizations diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cd3e66cefa..ade8834111 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2630,10 +2630,8 @@ contains filename = trim(path_input) // "tallies.xml" inquire(FILE=filename, EXIST=file_exists) if (.not. file_exists) then - ! We need each thread to allocate tally_derivs to avoid segfaults -!$omp parallel + ! We need to allocate tally_derivs to avoid segfaults allocate(tally_derivs(0)) -!$omp end parallel ! Since a tallies.xml file is optional, no error is issued here return @@ -2816,15 +2814,9 @@ contains ! ========================================================================== ! READ DATA FOR DERIVATIVES - ! Get pointer list to XML . + ! Get pointer list to XML nodes and allocate global array. call get_node_list(doc, "derivative", node_deriv_list) - - ! Allocate TallyDerivative array on each thread. The attributes of the - ! TallyDerivatives will be set on the master thread and then 'copyin'ed - ! in simulate.F90 -!$omp parallel allocate(tally_derivs(get_list_size(node_deriv_list))) -!$omp end parallel ! Read derivative attributes. do i = 1, get_list_size(node_deriv_list) diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 53a503ce2e..ffb0617e49 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -100,6 +100,9 @@ module particle_header integer(8) :: n_secondary = 0 type(Bank) :: secondary_bank(MAX_SECONDARY) + ! Flux (weight) derivatives for differential tallies + real(8), allocatable :: flux_derivs(:) + contains procedure :: initialize => initialize_particle procedure :: clear => clear_particle diff --git a/src/simulation.F90 b/src/simulation.F90 index 93950b1f9d..312b56d824 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -48,6 +48,9 @@ contains if (.not. restart_run) call initialize_source() + ! Allocate flux derivative array for differential tallies + allocate(p % flux_derivs(size(tally_derivs))) + ! Display header if (master) then if (run_mode == MODE_FIXEDSOURCE) then @@ -84,7 +87,7 @@ contains ! ==================================================================== ! LOOP OVER PARTICLES -!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs) +!$omp parallel do schedule(static) firstprivate(p) PARTICLE_LOOP: do i_work = 1, work current_work = i_work diff --git a/src/tally.F90 b/src/tally.F90 index 10ab35a797..889adb406d 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2743,24 +2743,27 @@ contains integer, intent(in) :: score_bin real(8), intent(inout) :: score - integer :: l ! loop index for nuclides in material + integer :: l logical :: scoring_diff_nuclide + real(8) :: flux_deriv real(8) :: dsigT, dsigA, dsigF, cum_dsig if (score == ZERO) return + flux_deriv = p % flux_derivs(t % deriv) + associate(deriv => tally_derivs(t % deriv)) - select case (deriv % variable) + select case (tally_derivs(t % deriv) % variable) case (DIFF_DENSITY) select case (t % estimator) case (ESTIMATOR_ANALOG) if (materials(p % material) % id == deriv % diff_material) then - score = score * (deriv % flux_deriv + ONE & + score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (ESTIMATOR_COLLISION) @@ -2768,52 +2771,52 @@ contains select case (score_bin) case (SCORE_FLUX) - score = score * deriv % flux_deriv + score = score * flux_deriv case (SCORE_TOTAL) if (materials(p % material) % id == deriv % diff_material & .and. material_xs % total /= ZERO) then - score = score * (deriv % flux_deriv + ONE & + score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_SCATTER) if (materials(p % material) % id == deriv % diff_material & .and. material_xs % total - material_xs % absorption /= ZERO) & then - score = score * (deriv % flux_deriv + ONE & + score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_ABSORPTION) if (materials(p % material) % id == deriv % diff_material & .and. material_xs % absorption /= ZERO) then - score = score * (deriv % flux_deriv + ONE & + score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_FISSION) if (materials(p % material) % id == deriv % diff_material & .and. material_xs % fission /= ZERO) then - score = score * (deriv % flux_deriv + ONE & + score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_NU_FISSION) if (materials(p % material) % id == deriv % diff_material & .and. material_xs % nu_fission /= ZERO) then - score = score * (deriv % flux_deriv + ONE & + score = score * (flux_deriv + ONE & / materials(p % material) % density_gpcc) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case default @@ -2836,11 +2839,11 @@ contains do l = 1, mat % n_nuclides if (mat % nuclide(l) == deriv % diff_nuclide) exit end do - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + ONE / mat % atom_density(l)) end associate else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (ESTIMATOR_COLLISION) @@ -2851,78 +2854,78 @@ contains select case (score_bin) case (SCORE_FLUX) - score = score * deriv % flux_deriv + score = score * flux_deriv case (SCORE_TOTAL) if (i_nuclide == -1 .and. & materials(p % material) % id == deriv % diff_material .and. & material_xs % total /= ZERO) then - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % total & / material_xs % total) else if (scoring_diff_nuclide .and. & micro_xs(deriv % diff_nuclide) % total /= ZERO) then - score = score * (deriv % flux_deriv + ONE / atom_density) + score = score * (flux_deriv + ONE / atom_density) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_SCATTER) if (i_nuclide == -1 .and. & materials(p % material) % id == deriv % diff_material .and. & material_xs % total - material_xs % absorption /= ZERO) then - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + (micro_xs(deriv % diff_nuclide) % total & - micro_xs(deriv % diff_nuclide) % absorption) & / (material_xs % total - material_xs % absorption)) else if (scoring_diff_nuclide .and. & (micro_xs(deriv % diff_nuclide) % total & - micro_xs(deriv % diff_nuclide) % absorption) /= ZERO) then - score = score * (deriv % flux_deriv + ONE / atom_density) + score = score * (flux_deriv + ONE / atom_density) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_ABSORPTION) if (i_nuclide == -1 .and. & materials(p % material) % id == deriv % diff_material .and. & material_xs % absorption /= ZERO) then - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % absorption & / material_xs % absorption ) else if (scoring_diff_nuclide .and. & micro_xs(deriv % diff_nuclide) % absorption /= ZERO) then - score = score * (deriv % flux_deriv + ONE / atom_density) + score = score * (flux_deriv + ONE / atom_density) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_FISSION) if (i_nuclide == -1 .and. & materials(p % material) % id == deriv % diff_material .and. & material_xs % fission /= ZERO) then - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % fission & / material_xs % fission) else if (scoring_diff_nuclide .and. & micro_xs(deriv % diff_nuclide) % fission /= ZERO) then - score = score * (deriv % flux_deriv + ONE / atom_density) + score = score * (flux_deriv + ONE / atom_density) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_NU_FISSION) if (i_nuclide == -1 .and. & materials(p % material) % id == deriv % diff_material .and. & material_xs % nu_fission /= ZERO) then - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + micro_xs(deriv % diff_nuclide) % nu_fission & / material_xs % nu_fission) else if (scoring_diff_nuclide .and. & micro_xs(deriv % diff_nuclide) % nu_fission /= ZERO) then - score = score * (deriv % flux_deriv + ONE / atom_density) + score = score * (flux_deriv + ONE / atom_density) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case default @@ -2943,7 +2946,7 @@ contains select case (score_bin) case (SCORE_FLUX) - score = score * deriv % flux_deriv + score = score * flux_deriv case (SCORE_TOTAL) if (materials(p % material) % id == deriv % diff_material .and. & @@ -2961,11 +2964,11 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigT * mat % atom_density(l) / material_xs % total) end associate else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_SCATTER) @@ -2986,12 +2989,12 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv + (dsigT - dsigA) & + score = score * (flux_deriv + (dsigT - dsigA) & * mat % atom_density(l) / & (material_xs % total - material_xs % absorption)) end associate else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_ABSORPTION) @@ -3010,11 +3013,11 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigA * mat % atom_density(l) / material_xs % absorption) end associate else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_FISSION) @@ -3033,11 +3036,11 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigF * mat % atom_density(l) / material_xs % fission) end associate else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_NU_FISSION) @@ -3056,13 +3059,13 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigF * mat % atom_density(l) / material_xs % nu_fission& * micro_xs(p % event_nuclide) % nu_fission & / micro_xs(p % event_nuclide) % fission) end associate else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case default @@ -3075,7 +3078,7 @@ contains select case (score_bin) case (SCORE_FLUX) - score = score * deriv % flux_deriv + score = score * flux_deriv case (SCORE_TOTAL) if (i_nuclide == -1 .and. & @@ -3096,7 +3099,7 @@ contains end associate end do end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + cum_dsig / material_xs % total) else if (materials(p % material) % id == deriv % diff_material & .and. material_xs % total > ZERO) then @@ -3109,10 +3112,10 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigT / micro_xs(i_nuclide) % total) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_SCATTER) @@ -3136,7 +3139,7 @@ contains end associate end do end associate - score = score * (deriv % flux_deriv + cum_dsig & + score = score * (flux_deriv + cum_dsig & / (material_xs % total - material_xs % absorption)) else if ( materials(p % material) % id == deriv % diff_material & .and. (material_xs % total - material_xs % absorption) > ZERO)& @@ -3151,11 +3154,11 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv + (dsigT - dsigA) & + score = score * (flux_deriv + (dsigT - dsigA) & / (micro_xs(i_nuclide) % total & - micro_xs(i_nuclide) % absorption)) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_ABSORPTION) @@ -3177,7 +3180,7 @@ contains end associate end do end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + cum_dsig / material_xs % absorption) else if (materials(p % material) % id == deriv % diff_material & .and. material_xs % absorption > ZERO) then @@ -3190,10 +3193,10 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigA / micro_xs(i_nuclide) % absorption) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_FISSION) @@ -3215,7 +3218,7 @@ contains end associate end do end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + cum_dsig / material_xs % fission) else if (materials(p % material) % id == deriv % diff_material & .and. material_xs % fission > ZERO) then @@ -3228,10 +3231,10 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigF / micro_xs(i_nuclide) % fission) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case (SCORE_NU_FISSION) @@ -3255,7 +3258,7 @@ contains end associate end do end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + cum_dsig / material_xs % nu_fission) else if (materials(p % material) % id == deriv % diff_material & .and. material_xs % nu_fission > ZERO) then @@ -3268,10 +3271,10 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) end if end associate - score = score * (deriv % flux_deriv & + score = score * (flux_deriv & + dsigF / micro_xs(i_nuclide) % fission) else - score = score * deriv % flux_deriv + score = score * flux_deriv end if case default @@ -3293,8 +3296,8 @@ contains !=============================================================================== subroutine score_track_derivative(p, distance) - type(particle), intent(in) :: p - real(8), intent(in) :: distance ! Neutron flight distance + type(particle), intent(inout) :: p + real(8), intent(in) :: distance ! Neutron flight distance integer :: i, l real(8) :: dsigT, dsigA, dsigF @@ -3312,7 +3315,7 @@ contains ! phi = e^(-Sigma_tot * dist) ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist - deriv % flux_deriv = deriv % flux_deriv & + p % flux_derivs(i) = p % flux_derivs(i) & - distance * material_xs % total / mat % density_gpcc end if end associate @@ -3323,7 +3326,7 @@ contains ! phi = e^(-Sigma_tot * dist) ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist - deriv % flux_deriv = deriv % flux_deriv & + p % flux_derivs(i) = p % flux_derivs(i) & - distance * micro_xs(deriv % diff_nuclide) % total end if end associate @@ -3338,7 +3341,7 @@ contains p % E <= nuc % multipole % end_E/1.0e6_8) then call multipole_deriv_eval(nuc % multipole, p % E, & p % sqrtkT, dsigT, dsigA, dsigF) - deriv % flux_deriv = deriv % flux_deriv & + p % flux_derivs(i) = p % flux_derivs(i) & - distance * dsigT * mat % atom_density(l) end if end associate @@ -3356,7 +3359,7 @@ contains !=============================================================================== subroutine score_collision_derivative(p) - type(particle), intent(in) :: p + type(particle), intent(inout) :: p integer :: i, j, l real(8) :: dsigT, dsigA, dsigF @@ -3374,7 +3377,7 @@ contains ! phi = Sigma_MT ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_MT / d_rho) / Sigma_MT ! (1 / phi) * (d_phi / d_rho) = 1 / rho - deriv % flux_deriv = deriv % flux_deriv & + p % flux_derivs(i) = p % flux_derivs(i) & + ONE / mat % density_gpcc end if end associate @@ -3395,7 +3398,7 @@ contains ! (1 / phi) * (d_phi / d_N) = (d_Sigma_MT / d_N) / Sigma_MT ! (1 / phi) * (d_phi / d_N) = sigma_MT / Sigma_MT ! (1 / phi) * (d_phi / d_N) = 1 / N - deriv % flux_deriv = deriv % flux_deriv & + p % flux_derivs(i) = p % flux_derivs(i) & + ONE / mat % atom_density(j) end if end associate @@ -3416,11 +3419,11 @@ contains p % sqrtkT, dsigT, dsigA, dsigF) select case(p % event) case (EVENT_SCATTER) - deriv % flux_deriv = deriv % flux_deriv + (dsigT - dsigA)& + p % flux_derivs(i) = p % flux_derivs(i) + (dsigT - dsigA)& / (micro_xs(mat % nuclide(l)) % total & - micro_xs(mat % nuclide(l)) % absorption) case (EVENT_ABSORB) - deriv % flux_deriv = deriv % flux_deriv & + p % flux_derivs(i) = p % flux_derivs(i) & + dsigA / micro_xs(mat % nuclide(l)) % absorption end select end if @@ -3433,17 +3436,6 @@ contains end do end subroutine score_collision_derivative -!=============================================================================== -! ZERO_FLUX_DERIVS Set the flux derivatives on differential tallies to zero. -!=============================================================================== - - subroutine zero_flux_derivs() - integer :: i - do i = 1, size(tally_derivs) - tally_derivs(i) % flux_deriv = ZERO - end do - end subroutine zero_flux_derivs - !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable diff --git a/src/tally_header.F90 b/src/tally_header.F90 index f30a40debe..5a80f249fc 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -25,7 +25,6 @@ module tally_header type TallyDerivative integer :: id - real(8) :: flux_deriv integer :: variable integer :: diff_material integer :: diff_nuclide diff --git a/src/tracking.F90 b/src/tracking.F90 index b0ba615ee5..7bdd32c227 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -16,7 +16,7 @@ module tracking use tally, only: score_analog_tally, score_tracklength_tally, & score_collision_tally, score_surface_current, & score_track_derivative, & - score_collision_derivative, zero_flux_derivs + score_collision_derivative use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -66,7 +66,7 @@ contains endif ! Every particle starts with no accumulated flux derivative. - if (active_tallies % size() > 0) call zero_flux_derivs() + p % flux_derivs(:) = ZERO EVENT_LOOP: do ! If the cell hasn't been determined based on the particle's location,