From b2d679c996800a4d272da26e89c72a61f52d5115 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 24 Feb 2014 11:44:48 -0500 Subject: [PATCH 1/8] Avoided NaNs in keff and fission analog tallies when micro_xs(i)%absorption was zero as it is for He-4. --- src/physics.F90 | 6 +++++- src/tally.F90 | 54 ++++++++++++++++++++++++++++++------------------- 2 files changed, 38 insertions(+), 22 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 2d3cfeb39f..1bc975f968 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -98,7 +98,11 @@ contains ! If survival biasing is being used, the following subroutine adjusts the ! weight of the particle. Otherwise, it checks to see if absorption occurs - call absorption(p, i_nuclide) + if (micro_xs(i_nuclide) % absorption > ZERO) then + call absorption(p, i_nuclide) + else + p % absorb_wgt = ZERO + end if if (.not. p % alive) return ! Play russian roulette if survival biasing is turned on diff --git a/src/tally.F90 b/src/tally.F90 index b731c1ae5e..9e00ea2bd1 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -42,7 +42,7 @@ 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 + 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 @@ -164,7 +164,7 @@ contains score = last_wgt - case (SCORE_NU_SCATTER) + case (SCORE_NU_SCATTER) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -173,7 +173,7 @@ contains ! reaction with neutrons in the exit channel score = wgt - + case (SCORE_SCATTER_N) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -201,7 +201,7 @@ contains score_index = score_index + 1 ! get the score and tally it score = last_wgt * calc_pn(n, mu) - + !$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score @@ -257,8 +257,12 @@ contains ! calculate fraction of absorptions that would have resulted in ! fission - score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission / & + micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if else ! Skip any non-fission events @@ -277,8 +281,12 @@ contains ! calculate fraction of absorptions that would have resulted in ! nu-fission - score = p % absorb_wgt * micro_xs(p % event_nuclide) % & - nu_fission / micro_xs(p % event_nuclide) % absorption + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * micro_xs(p % event_nuclide) % & + nu_fission / micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if else ! Skip any non-fission events @@ -305,17 +313,21 @@ contains end if end if - + case (SCORE_KAPPA_FISSION) if (survival_biasing) then ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in ! fission and multiply by Q - score = p % absorb_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption - + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * & + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if + else ! Skip any non-fission events if (.not. p % fission) cycle SCORE_LOOP @@ -323,7 +335,7 @@ contains ! All fission events will contribute, so again we can use ! particle's weight entering the collision as the estimate for ! the fission energy production rate - + n = nuclides(p % event_nuclide) % index_fission(1) score = last_wgt * & nuclides(p % event_nuclide) % reactions(n) % Q_value @@ -665,7 +677,7 @@ contains do l = 1, mat % n_nuclides ! Get atom density atom_density = mat % atom_density(l) - + ! Get index in nuclides array i_nuc = mat % nuclide(l) @@ -780,7 +792,7 @@ contains ! determine what type of score bin score_bin = t % score_bins(j) - ! Determine macroscopic nuclide cross section + ! Determine macroscopic nuclide cross section select case(score_bin) case (SCORE_FLUX) score = flux @@ -812,7 +824,7 @@ contains ! sections that are used often (e.g. n2n, ngamma, etc. for depletion), ! it might make sense to optimize this section or pre-calculate cross ! sections - + if (score_bin > 1) then ! Set default score score = ZERO @@ -869,7 +881,7 @@ contains ! determine what type of score bin score_bin = t % score_bins(j) - ! Determine macroscopic material cross section + ! Determine macroscopic material cross section select case(score_bin) case (SCORE_FLUX) score = flux @@ -899,7 +911,7 @@ contains ! Any other cross section has to be calculated on-the-fly. This is ! somewhat costly since it requires a loop over each nuclide in a ! material and each reaction in the nuclide - + if (score_bin > 1) then ! Set default score score = ZERO @@ -1202,7 +1214,7 @@ contains score_bin = t % score_bins(j) if (i_nuclide > 0) then - ! Determine macroscopic nuclide cross section + ! Determine macroscopic nuclide cross section select case(score_bin) case (SCORE_FLUX) score = flux @@ -1233,7 +1245,7 @@ contains end select else - ! Determine macroscopic material cross section + ! Determine macroscopic material cross section select case(score_bin) case (SCORE_FLUX) score = flux From 594bd66c86dc962007bf9c29e59e1736dffd31ef Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 24 Feb 2014 14:52:40 -0500 Subject: [PATCH 2/8] Fixed nu-scatter (w/ energyout filter) survival biasing issue. --- src/tally.F90 | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 9e00ea2bd1..cd65d6dadb 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -220,7 +220,7 @@ contains ! Score total rate - p1 scatter rate Note estimator needs to be ! adjusted since tallying is only occuring when a scatter has - ! happend. Effectively this means multiplying the estimator by + ! happened. Effectively this means multiplying the estimator by ! total/scatter macro score = (macro_total - mu*macro_scatt)*(ONE/macro_scatt) @@ -281,11 +281,24 @@ contains ! calculate fraction of absorptions that would have resulted in ! nu-fission - if (micro_xs(p % event_nuclide) % absorption > ZERO) then - score = p % absorb_wgt * micro_xs(p % event_nuclide) % & - nu_fission / micro_xs(p % event_nuclide) % absorption + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + ! Normally, we only need to make contributions to one scoring + ! bin. However, in the case of fission, since multiple fission + ! neutrons were emitted with different energies, multiple + ! outgoing energy bins may have been scored to. The following + ! logic treats this special case and results to multiple bins + + call score_fission_eout(p, t, score_index) + cycle SCORE_LOOP + else - score = ZERO + + if (micro_xs(p % event_nuclide) % absorption > ZERO) then + score = p % absorb_wgt * micro_xs(p % event_nuclide) % & + nu_fission / micro_xs(p % event_nuclide) % absorption + else + score = ZERO + end if end if else @@ -395,7 +408,7 @@ contains integer :: k ! loop index for bank sites integer :: bin_energyout ! original outgoing energy bin integer :: i_filter ! index for matching filter bin combination - real(8) :: score ! actualy score + real(8) :: score ! actual score real(8) :: E_out ! energy of fission bank site ! save original outgoing energy bin and score index From 3ef6249f5c03a5a7ee3734160196cc7310e8a66d Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 10 Apr 2014 13:49:58 -0400 Subject: [PATCH 3/8] Fixed erroneous tallying of fission and kappa-fission in cases without survival biasing. --- src/tally.F90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index cd65d6dadb..ef5ae811ee 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -259,20 +259,21 @@ contains if (micro_xs(p % event_nuclide) % absorption > ZERO) then score = p % absorb_wgt * micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption + micro_xs(p % event_nuclide) % absorption else score = ZERO end if else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP + ! Skip any non-absorption events + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP ! All fission events will contribute, so again we can use ! particle's weight entering the collision as the estimate for the ! fission reaction rate - score = last_wgt + score = last_wgt * micro_xs(p % event_nuclide) % fission / & + micro_xs(p % event_nuclide) % absorption end if case (SCORE_NU_FISSION) @@ -335,23 +336,22 @@ contains if (micro_xs(p % event_nuclide) % absorption > ZERO) then score = p % absorb_wgt * & - micro_xs(p % event_nuclide) % kappa_fission / & - micro_xs(p % event_nuclide) % absorption + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption else score = ZERO end if else - ! Skip any non-fission events - if (.not. p % fission) cycle SCORE_LOOP + ! Skip any non-absorption events + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP ! All fission events will contribute, so again we can use ! particle's weight entering the collision as the estimate for ! the fission energy production rate - - n = nuclides(p % event_nuclide) % index_fission(1) score = last_wgt * & - nuclides(p % event_nuclide) % reactions(n) % Q_value + micro_xs(p % event_nuclide) % kappa_fission / & + micro_xs(p % event_nuclide) % absorption end if case (SCORE_EVENTS) ! Simply count number of scoring events From 047e7fe3ef371e760848b4de804467600681dc45 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 10 Apr 2014 14:13:26 -0400 Subject: [PATCH 4/8] Fixed flux tallying for case of survival biasing. --- src/tally.F90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/tally.F90 b/src/tally.F90 index ef5ae811ee..1e7e2e5998 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -139,7 +139,15 @@ contains ! estimator since there is no way to count 'events' exactly for ! the flux - score = last_wgt / material_xs % total + if (survival_biasing) then + ! We need to account for the fact that some weight was already + ! absorbed + score = last_wgt + p % absorb_wgt + else + score = last_wgt + end if + + score = score / material_xs % total case (SCORE_TOTAL) ! All events will score to the total reaction rate. We can just From 56d0b3b89f2018554e91d43a4630e41504dbfe99 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 26 Jun 2014 13:06:46 -0400 Subject: [PATCH 5/8] Moved scattering to occur before russian roulette, and set last_wgt to 0 as well when a particle is killed via russian roulette. --- src/physics.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 940239fc62..71a5cc6cc3 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -105,6 +105,10 @@ contains end if if (.not. p % alive) return + ! Sample a scattering reaction and determine the secondary energy of the + ! exiting neutron + call scatter(p, i_nuclide) + ! Play russian roulette if survival biasing is turned on if (survival_biasing) then @@ -112,11 +116,6 @@ contains if (.not. p % alive) return end if - ! Sample a scattering reaction and determine the secondary energy of the - ! exiting neutron - - call scatter(p, i_nuclide) - end subroutine sample_reaction !=============================================================================== @@ -297,6 +296,7 @@ contains p % last_wgt = p % wgt else p % wgt = ZERO + p % last_wgt = ZERO p % alive = .false. end if end if From c9561a8884890555ce14aa867d5448469275db9c Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 27 Jun 2014 06:51:37 -0400 Subject: [PATCH 6/8] Fixed a merge conflict issue. --- src/tally.F90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 0b4c23daed..f86d67d168 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -4,7 +4,7 @@ module tally use constants use error, only: fatal_error use global - use math, only: t_percentile, calc_pn + use math, only: t_percentile, calc_pn, calc_rn use mesh, only: get_mesh_bin, bin_to_mesh_indices, & get_mesh_indices, mesh_indices_to_bin, & mesh_intersects_2d, mesh_intersects_3d @@ -244,7 +244,7 @@ contains score = last_wgt - case (SCORE_NU_SCATTER) + case (SCORE_NU_SCATTER) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -253,7 +253,7 @@ contains ! reaction with neutrons in the exit channel score = wgt - + case (SCORE_SCATTER_N) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -281,7 +281,7 @@ contains score_index = score_index + 1 ! get the score and tally it score = last_wgt * calc_pn(n, mu) - + !$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score @@ -503,7 +503,7 @@ contains end if end if - + case (SCORE_KAPPA_FISSION) if (survival_biasing) then ! No fission events occur if survival biasing is on -- need to @@ -977,7 +977,7 @@ contains do l = 1, mat % n_nuclides ! Get atom density atom_density = mat % atom_density(l) - + ! Get index in nuclides array i_nuc = mat % nuclide(l) @@ -1179,7 +1179,7 @@ contains ! sections that are used often (e.g. n2n, ngamma, etc. for depletion), ! it might make sense to optimize this section or pre-calculate cross ! sections - + if (score_bin > 1) then ! Set default score score = ZERO @@ -1318,7 +1318,7 @@ contains ! Any other cross section has to be calculated on-the-fly. This is ! somewhat costly since it requires a loop over each nuclide in a ! material and each reaction in the nuclide - + if (score_bin > 1) then ! Set default score score = ZERO @@ -1630,7 +1630,7 @@ contains score_index = (b - 1)*t % n_score_bins + j if (i_nuclide > 0) then - ! Determine macroscopic nuclide cross section + ! Determine macroscopic nuclide cross section select case(score_bin) case (SCORE_FLUX) score = flux @@ -1713,7 +1713,7 @@ contains end select else - ! Determine macroscopic material cross section + ! Determine macroscopic material cross section select case(score_bin) case (SCORE_FLUX) score = flux From 72b9cb7a55fd1effc207fcf7fa16ac6fd6b4efd7 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 27 Jun 2014 07:01:25 -0400 Subject: [PATCH 7/8] Updated survival biasing test results. --- tests/test_survival_biasing/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_survival_biasing/results_true.dat b/tests/test_survival_biasing/results_true.dat index 004b3ac613..78f87d0227 100644 --- a/tests/test_survival_biasing/results_true.dat +++ b/tests/test_survival_biasing/results_true.dat @@ -1,2 +1,2 @@ k-combined: -9.681467E-01 8.159286E-03 +1.043266E+00 1.481033E-02 From fb3c78233917ff151564685bd45d609871bbdc7b Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 27 Jun 2014 07:50:07 -0400 Subject: [PATCH 8/8] Updating test score fission results based on survival biasing updates. --- tests/test_score_fission/results_true.dat | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_score_fission/results_true.dat b/tests/test_score_fission/results_true.dat index dfcab1a2a5..e9440264c9 100644 --- a/tests/test_score_fission/results_true.dat +++ b/tests/test_score_fission/results_true.dat @@ -10,11 +10,11 @@ tally 1: 1.114257E+00 2.508773E-01 tally 2: -1.180000E+00 -2.886000E-01 +7.809260E-01 +1.258561E-01 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 -1.680000E+00 -5.682000E-01 +1.056654E+00 +2.244276E-01