diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 7372f9e84a..d05dcdf71b 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1392,13 +1392,15 @@ 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`` - 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. @@ -1484,6 +1486,14 @@ 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..375c517e76 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -60,20 +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 - 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 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 + 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 @@ -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 627889db5c..58e4013355 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3128,6 +3128,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 bfa5108a00..f2c261ec75 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -654,6 +654,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 4d4d6bf938..0586860191 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 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 + ! 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(TWO * p % E / (MASS_NEUTRON_MEV)) * C_LIGHT) + + else + ! For inverse velocity, we need no cross section + score = flux / (sqrt(TWO * p % E / (MASS_NEUTRON_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()