From b13061b5c079007e9755ca028c31c70ba14de255 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 22 Oct 2015 18:39:49 -0400 Subject: [PATCH 1/4] added inverse-velocity score and accompanying tests --- docs/source/usersguide/input.rst | 9 ++++-- src/constants.F90 | 7 ++-- src/input_xml.F90 | 2 ++ src/output.F90 | 1 + src/state_point.F90 | 2 ++ src/summary.F90 | 2 ++ src/tally.F90 | 21 ++++++++++++ .../inputs_true.dat | 1 + .../results_true.dat | 29 +++++++++++++++++ .../test_score_inversevelocity.py | 32 +++++++++++++++++++ 10 files changed, 102 insertions(+), 4 deletions(-) create mode 100644 tests/test_score_inverse_velocity/inputs_true.dat create mode 100644 tests/test_score_inverse_velocity/results_true.dat create mode 100644 tests/test_score_inverse_velocity/test_score_inversevelocity.py diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 0c10756ba1..8ddf7082c5 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1386,8 +1386,8 @@ The ```` element accepts the following sub-elements: options are "flux", "total", "scatter", "absorption", "fission", "nu-fission", "delayed-nu-fission", "kappa-fission", "nu-scatter", "scatter-N", "scatter-PN", "scatter-YN", "nu-scatter-N", "nu-scatter-PN", - "nu-scatter-YN", "flux-YN", "total-YN", "current", and "events". These - correspond to the following physical quantities: + "nu-scatter-YN", "flux-YN", "total-YN", "current", "inverse-velocity" and + "events". These correspond to the following physical quantities: :flux: Total flux in particle-cm per source particle. Note: The ``analog`` @@ -1478,6 +1478,11 @@ The ```` element accepts the following sub-elements: specified. Furthermore, it may not be used in conjunction with any other score. + :inverse-velocity: + The flux-weighted inverse velocity where the velocity is in units of + meters per second. Note: The ``analog`` estimator is actually identical + to the ``collision`` estimator for the inverse-velocity score. + :events: Number of scoring events. Units are events per source particle. diff --git a/src/constants.F90 b/src/constants.F90 index 9e973f24af..7c8bc8fc5d 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -65,6 +65,8 @@ module constants MASS_NEUTRON = 1.008664916_8, & ! mass of a neutron in amu MASS_PROTON = 1.007276466812_8, & ! mass of a proton in amu AMU = 1.660538921e-27_8, & ! 1 amu in kg + AMU_MEV = 931.494061_8, & ! 1 amu in MeV/c^2 + C_LIGHT = 2.99792458e8_8, & ! speed of light in a vacuum N_AVOGADRO = 0.602214129_8, & ! Avogadro's number in 10^24/mol K_BOLTZMANN = 8.6173324e-11_8, & ! Boltzmann constant in MeV/K INFINITY = huge(0.0_8), & ! positive infinity @@ -257,7 +259,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 21 + integer, parameter :: N_SCORE_TYPES = 22 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -279,7 +281,8 @@ module constants SCORE_SCATTER_YN = -18, & ! angular flux-weighted scattering moment (0:N) SCORE_NU_SCATTER_YN = -19, & ! angular flux-weighted nu-scattering moment (0:N) SCORE_EVENTS = -20, & ! number of events - SCORE_DELAYED_NU_FISSION = -21 ! delayed neutron production rate + SCORE_DELAYED_NU_FISSION = -21, & ! delayed neutron production rate + SCORE_INVERSE_VELOCITY = -22 ! flux weighted inverse velocity ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 528ea40f51..3fe344ed92 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3110,6 +3110,8 @@ contains end if case ('kappa-fission') t % score_bins(j) = SCORE_KAPPA_FISSION + case ('inverse-velocity') + t % score_bins(j) = SCORE_INVERSE_VELOCITY case ('current') t % score_bins(j) = SCORE_CURRENT t % type = TALLY_SURFACE_CURRENT diff --git a/src/output.F90 b/src/output.F90 index 3463244e31..55cd5b2a5b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -985,6 +985,7 @@ contains score_names(abs(SCORE_NU_SCATTER_PN)) = "Scattering Prod. Rate Moment" score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment" score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate" + score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/state_point.F90 b/src/state_point.F90 index 6d1c6f3129..046e7e6669 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -353,6 +353,8 @@ contains str_array(j) = "nu-scatter-yn" case (SCORE_EVENTS) str_array(j) = "events" + case (SCORE_INVERSE_VELOCITY) + str_array(j) = "inverse-velocity" case default str_array(j) = reaction_name(tally%score_bins(j)) end select diff --git a/src/summary.F90 b/src/summary.F90 index f088cb7bbf..bf7621880b 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -649,6 +649,8 @@ contains str_array(j) = "nu-scatter-yn" case (SCORE_EVENTS) str_array(j) = "events" + case (SCORE_INVERSE_VELOCITY) + str_array(j) = "inverse-velocity" case default str_array(j) = reaction_name(t%score_bins(j)) end select diff --git a/src/tally.F90 b/src/tally.F90 index 7ad174c9bc..6e2e279120 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -127,6 +127,27 @@ contains end if + case (SCORE_INVERSE_VELOCITY) + if (t % estimator == ESTIMATOR_ANALOG) then + ! All events score to a inverse velocity bin. We actually use a + ! collision estimator in place of an analog one since there is no way + ! to count 'events' exactly for the inverse velocity + if (survival_biasing) then + ! We need to account for the fact that some weight was already + ! absorbed + score = p % last_wgt + p % absorb_wgt + else + score = p % last_wgt + end if + score = score / material_xs % total & + / (sqrt(2 * p % E / (MASS_NEUTRON * AMU_MEV)) * C_LIGHT) + + else + ! For inverse velocity, we need no cross section + score = flux / (sqrt(2 * p % E / (MASS_NEUTRON * AMU_MEV)) * C_LIGHT) + end if + + case (SCORE_SCATTER, SCORE_SCATTER_N) if (t % estimator == ESTIMATOR_ANALOG) then ! Skip any event where the particle didn't scatter diff --git a/tests/test_score_inverse_velocity/inputs_true.dat b/tests/test_score_inverse_velocity/inputs_true.dat new file mode 100644 index 0000000000..0546cf828f --- /dev/null +++ b/tests/test_score_inverse_velocity/inputs_true.dat @@ -0,0 +1 @@ +ba1010f940c50314d61aae9f729b7bb476a6b35c5556fbc689ba3fde70ff50b0ba5dad0db3b38349a1562d67817047090ac450a64d895c8f3393163fe7914763 \ No newline at end of file diff --git a/tests/test_score_inverse_velocity/results_true.dat b/tests/test_score_inverse_velocity/results_true.dat new file mode 100644 index 0000000000..f86a8beb5b --- /dev/null +++ b/tests/test_score_inverse_velocity/results_true.dat @@ -0,0 +1,29 @@ +k-combined: +9.903196E-01 4.279617E-02 +tally 1: +1.049628E-03 +2.261930E-07 +4.056389E-04 +3.411247E-08 +2.243766E-03 +1.069671E-06 +6.354432E-04 +8.370608E-08 +tally 2: +1.048031E-03 +2.263789E-07 +4.276093E-04 +4.136358E-08 +2.667795E-03 +1.528407E-06 +6.199140E-04 +7.840402E-08 +tally 3: +1.048031E-03 +2.263789E-07 +4.276093E-04 +4.136358E-08 +2.667795E-03 +1.528407E-06 +6.199140E-04 +7.840402E-08 diff --git a/tests/test_score_inverse_velocity/test_score_inversevelocity.py b/tests/test_score_inverse_velocity/test_score_inversevelocity.py new file mode 100644 index 0000000000..787e9df091 --- /dev/null +++ b/tests/test_score_inverse_velocity/test_score_inversevelocity.py @@ -0,0 +1,32 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness, PyAPITestHarness +import openmc + + +class ScoreInverseVelocityTestHarness(PyAPITestHarness): + def _build_inputs(self): + filt = openmc.Filter(type='cell', bins=(21, 22, 23, 27)) + tallies = [openmc.Tally(tally_id=i) for i in range(1, 4)] + [t.add_filter(filt) for t in tallies] + [t.add_score('inverse-velocity') for t in tallies] + tallies[0].estimator = 'tracklength' + tallies[1].estimator = 'analog' + tallies[2].estimator = 'collision' + self._input_set.tallies = openmc.TalliesFile() + [self._input_set.tallies.add_tally(t) for t in tallies] + + super(ScoreInverseVelocityTestHarness, self)._build_inputs() + + def _cleanup(self): + super(ScoreInverseVelocityTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + + +if __name__ == '__main__': + harness = ScoreInverseVelocityTestHarness('statepoint.10.*', True) + harness.main() From 9828b34ec182ca1e3c210cb1edbd76eb5668a614 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 22 Oct 2015 21:28:28 -0400 Subject: [PATCH 2/4] added dash in comment for score inverse velocity --- src/constants.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constants.F90 b/src/constants.F90 index 7c8bc8fc5d..2aceea15cf 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -282,7 +282,7 @@ module constants SCORE_NU_SCATTER_YN = -19, & ! angular flux-weighted nu-scattering moment (0:N) SCORE_EVENTS = -20, & ! number of events SCORE_DELAYED_NU_FISSION = -21, & ! delayed neutron production rate - SCORE_INVERSE_VELOCITY = -22 ! flux weighted inverse velocity + SCORE_INVERSE_VELOCITY = -22 ! flux-weighted inverse velocity ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 From 6038733a2234d8889503846fdda2e4a5aea1eceb Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Fri, 23 Oct 2015 18:54:36 -0400 Subject: [PATCH 3/4] changed constant for mass of neutron in MeV/c^2 and modified inverse-velocity documentation --- docs/source/usersguide/input.rst | 15 ++++++++++----- src/constants.F90 | 32 ++++++++++++++++---------------- src/tally.F90 | 4 ++-- 3 files changed, 28 insertions(+), 23 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 8e52813bc2..726624cc85 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1390,9 +1390,11 @@ The ```` element accepts the following sub-elements: "events". These correspond to the following physical quantities: :flux: - Total flux in particle-cm per source particle. Note: The ``analog`` - estimator is actually identical to the ``collision`` estimator for the - flux score. + Total flux in particle-cm per source particle. + + .. note:: + The ``analog`` estimator is actually identical to the ``collision`` + estimator for the flux score. :total: Total reaction rate in reactions per source particle. @@ -1480,8 +1482,11 @@ The ```` element accepts the following sub-elements: :inverse-velocity: The flux-weighted inverse velocity where the velocity is in units of - meters per second. Note: The ``analog`` estimator is actually identical - to the ``collision`` estimator for the inverse-velocity score. + meters per second. + + .. note:: + The ``analog`` estimator is actually identical to the ``collision`` + estimator for the inverse-velocity score. :events: Number of scoring events. Units are events per source particle. diff --git a/src/constants.F90 b/src/constants.F90 index 2aceea15cf..19fccf2677 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -60,22 +60,22 @@ module constants ! Values here are from the Committee on Data for Science and Technology ! (CODATA) 2010 recommendation (doi:10.1103/RevModPhys.84.1527). - real(8), parameter :: & - PI = 3.1415926535898_8, & ! pi - MASS_NEUTRON = 1.008664916_8, & ! mass of a neutron in amu - MASS_PROTON = 1.007276466812_8, & ! mass of a proton in amu - AMU = 1.660538921e-27_8, & ! 1 amu in kg - AMU_MEV = 931.494061_8, & ! 1 amu in MeV/c^2 - C_LIGHT = 2.99792458e8_8, & ! speed of light in a vacuum - N_AVOGADRO = 0.602214129_8, & ! Avogadro's number in 10^24/mol - K_BOLTZMANN = 8.6173324e-11_8, & ! Boltzmann constant in MeV/K - INFINITY = huge(0.0_8), & ! positive infinity - ZERO = 0.0_8, & - HALF = 0.5_8, & - ONE = 1.0_8, & - TWO = 2.0_8, & - THREE = 3.0_8, & - FOUR = 4.0_8 + real(8), parameter :: & + PI = 3.1415926535898_8, & ! pi + MASS_NEUTRON = 1.008664916_8, & ! mass of a neutron in amu + MASS_NEUTRON_MEV = 939.565379_8, & ! mass of a neutron in MeV/c^2 + MASS_PROTON = 1.007276466812_8, & ! mass of a proton in amu + AMU = 1.660538921e-27_8, & ! 1 amu in kg + C_LIGHT = 2.99792458e8_8, & ! speed of light in a vacuum + N_AVOGADRO = 0.602214129_8, & ! Avogadro's number in 10^24/mol + K_BOLTZMANN = 8.6173324e-11_8, & ! Boltzmann constant in MeV/K + INFINITY = huge(0.0_8), & ! positive infinity + ZERO = 0.0_8, & + HALF = 0.5_8, & + ONE = 1.0_8, & + TWO = 2.0_8, & + THREE = 3.0_8, & + FOUR = 4.0_8 ! ============================================================================ ! GEOMETRY-RELATED CONSTANTS diff --git a/src/tally.F90 b/src/tally.F90 index 475cc09b10..ee07db56c1 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -140,11 +140,11 @@ contains score = p % last_wgt end if score = score / material_xs % total & - / (sqrt(2 * p % E / (MASS_NEUTRON * AMU_MEV)) * C_LIGHT) + / (sqrt(TWO * p % E / (MASS_NEUTRON_MEV)) * C_LIGHT) else ! For inverse velocity, we need no cross section - score = flux / (sqrt(2 * p % E / (MASS_NEUTRON * AMU_MEV)) * C_LIGHT) + score = flux / (sqrt(TWO * p % E / (MASS_NEUTRON_MEV)) * C_LIGHT) end if From 5799c05be4deba453ce3c24c22586c4f0fdadc7b Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Sun, 25 Oct 2015 17:54:15 -0400 Subject: [PATCH 4/4] added units to speed of light constant --- src/constants.F90 | 2 +- src/tally.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 19fccf2677..375c517e76 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -66,7 +66,7 @@ module constants MASS_NEUTRON_MEV = 939.565379_8, & ! mass of a neutron in MeV/c^2 MASS_PROTON = 1.007276466812_8, & ! mass of a proton in amu AMU = 1.660538921e-27_8, & ! 1 amu in kg - C_LIGHT = 2.99792458e8_8, & ! speed of light in a vacuum + C_LIGHT = 2.99792458e8_8, & ! speed of light in m/s N_AVOGADRO = 0.602214129_8, & ! Avogadro's number in 10^24/mol K_BOLTZMANN = 8.6173324e-11_8, & ! Boltzmann constant in MeV/K INFINITY = huge(0.0_8), & ! positive infinity diff --git a/src/tally.F90 b/src/tally.F90 index fa1c13965a..0586860191 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -129,7 +129,7 @@ contains case (SCORE_INVERSE_VELOCITY) if (t % estimator == ESTIMATOR_ANALOG) then - ! All events score to a inverse velocity bin. We actually use a + ! All events score to an inverse velocity bin. We actually use a ! collision estimator in place of an analog one since there is no way ! to count 'events' exactly for the inverse velocity if (survival_biasing) then