diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 6e29e7ee5f..c9aa8cc866 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -41,7 +41,7 @@ for the geometry, materials, and settings. Additionally, there are three optiona input files. The first is a tallies XML file that specifies physical quantities to be tallied. The second is a plots XML file that specifies regions of geometry which should be plotted. The third is a CMFD XML file that specifies coarse mesh -acceleration geometry and execution parameters. OpenMC expects that these +acceleration geometry and execution parameters. OpenMC expects that these files are called: * ``geometry.xml`` @@ -105,7 +105,7 @@ default. This element has the following attributes/sub-elements: The ```` element indicates that a :math:`k`-eigenvalue calculation should be performed. It has the following attributes/sub-elements: - :batches: + :batches: The total number of batches, where each batch corresponds to multiple fission source iterations. Batching is done to eliminate correlation between realizations of random variables. @@ -171,7 +171,7 @@ problem. It has the following attributes/sub-elements: The ```` element indicates that a fixed source calculation should be performed. It has the following attributes/sub-elements: - :batches: + :batches: The total number of batches. For fixed source calculations, each batch represents a realization of random variables for tallies. @@ -229,7 +229,7 @@ out the file and "false" will not. *Default*: false - :summary: + :summary: Writes out an ASCII summary file describing all of the user input files that were read in. @@ -297,7 +297,7 @@ attributes/sub-elements: An element specifying the spatial distribution of source sites. This element has the following attributes: - :type: + :type: The type of spatial distribution. Valid options are "box" and "point". A "box" spatial distribution has coordinates sampled uniformly in a parallelepiped. A "point" spatial distribution has coordinates specified @@ -321,7 +321,7 @@ attributes/sub-elements: An element specifying the angular distribution of source sites. This element has the following attributes: - :type: + :type: The type of angular distribution. Valid options are "isotropic" and "monodirectional". The angle of the particle emitted from a source site is isotropic if the "isotropic" option is given. The angle of the particle @@ -344,7 +344,7 @@ attributes/sub-elements: An element specifying the energy distribution of source sites. This element has the following attributes: - :type: + :type: The type of energy distribution. Valid options are "monoenergetic", "watt", and "maxwell". The "monoenergetic" option produces source sites at @@ -373,8 +373,8 @@ attributes/sub-elements: The ```` element indicates at what batches a state point file should be written. A state point file can be used to restart a run or to get -tally results at any batch. The default behavior when using this tag is to -write out the source bank in the state_point file. This behavior can be +tally results at any batch. The default behavior when using this tag is to +write out the source bank in the state_point file. This behavior can be customized by using the ```` element. This element has the following attributes/sub-elements: @@ -394,8 +394,8 @@ following attributes/sub-elements: ```` Element -------------------------- -The ```` element indicates at what batches the source bank -should be written. The source bank can be either written out within a state +The ```` element indicates at what batches the source bank +should be written. The source bank can be either written out within a state point file or separately in a source point file. This element has the following attributes/sub-elements: @@ -425,17 +425,17 @@ attributes/sub-elements: :source_write: If this element is set to "false", source sites are not written - to the state point or source point file. This can substantially reduce the + to the state point or source point file. This can substantially reduce the size of state points if large numbers of particles per batch are used. *Default*: true :overwrite_latest: If this element is set to "true", a source point file containing - the source bank will be written out to a separate file named - ``source.binary`` or ``source.h5`` depending on if HDF5 is enabled. + the source bank will be written out to a separate file named + ``source.binary`` or ``source.h5`` depending on if HDF5 is enabled. This file will be overwritten at every single batch so that the latest - source bank will be available. It should be noted that a user can set both + source bank will be available. It should be noted that a user can set both this element to "true" and specify batches to write a permanent source bank. *Default*: false @@ -879,7 +879,7 @@ The ```` element accepts the following sub-elements: :label: This is an optional sub-element specifying the name of this tally to be used - for output purposes. This string is limited to 52 characters for formatting + for output purposes. This string is limited to 52 characters for formatting purposes. :filter: @@ -949,10 +949,11 @@ The ```` element accepts the following sub-elements: :scores: A space-separated list of the desired responses to be accumulated. Accepted - options are "flux", "total", "scatter", "nu-scatter", "scatter-N", - "scatter-PN", "absorption", "fission", "nu-fission", "kappa-fission", - "current", and "events". These corresponding to the following physical - quantities. + options are "flux", "total", "scatter", "absorption", "fission", + "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 corresponding to the following + physical quantities: :flux: Total flux @@ -964,24 +965,6 @@ The ```` element accepts the following sub-elements: Total scattering rate. Can also be identified with the ``scatter-0`` response type. - :nu-scatter: - Total production of neutrons due to scattering. This accounts for - multiplicity from (n,2n), (n,3n), and (n,4n) reactions and should be - slightly higher than the scattering rate. - - :scatter-N: - Tally the N\ :sup:`th` \ scattering moment, where N is the Legendre - expansion order. N must be between 0 and 10. As an example, tallying the - 2\ :sup:`nd` \ scattering moment would be specified as `` - scatter-2 ``. - - :scatter-PN: - Tally all of the scattering moments from order 0 to N, where N is the - Legendre expansion order. That is, ``scatter-P1`` is equivalent to - requesting tallies of ``scatter-0`` and ``scatter-1``. N must be between - 0 and 10. As an example, tallying up to the 2\ :sup:`nd` \ scattering - moment would be specified as `` scatter-P2 ``. - :absorption: Total absorption rate. This accounts for all reactions which do not produce secondary neutrons. @@ -991,7 +974,7 @@ The ```` element accepts the following sub-elements: :nu-fission: Total production of neutrons due to fission - + :kappa-fission: The recoverable energy production rate due to fission. The recoverable energy is defined as the fission product kinetic energy, prompt and @@ -1001,6 +984,47 @@ The ```` element accepts the following sub-elements: prompt and delayed :math:`\gamma`-rays are assumed to deposit their energy locally. + :scatter-N: + Tally the N\ :sup:`th` \ scattering moment, where N is the Legendre + expansion order of the change in particle angle :math:`\left(\mu\right)`. + N must be between 0 and 10. As an example, tallying the + 2\ :sup:`nd` \ scattering moment would be specified as + `` scatter-2 ``. + + :scatter-PN: + Tally all of the scattering moments from order 0 to N, where N is the + Legendre expansion order of the change in particle angle :math:`\left(\mu\right)`. + That is, ``scatter-P1`` is equivalent to requesting tallies of + ``scatter-0`` and ``scatter-1``. Like for ``scatter-N``, + N must be between 0 and 10. As an example, tallying up to the + 2\ :sup:`nd` \ scattering moment would be specified as + `` scatter-P2 ``. + + :scatter-YN: + ``scatter-YN`` is similar to ``scatter-PN`` except an additional + expansion is performed for the incoming particle direction + :math:`\left(\Omega\right)` using the real spherical harmonics. This is useful + for performing angular flux moment weighting of the scattering moments. + Like ``scatter-PN``, ``scatter-YN`` will tally all of the moments from + order 0 to N; N again must be between 0 and 10. + + :nu-scatter, nu-scatter-N, nu-scatter-PN, nu-scatter-YN: + These scores are similar in functionality to their ``scatter*`` + equivalents except the total production of neutrons due to + scattering is scored vice simply the scattering rate. This accounts for + multiplicity from (n,2n), (n,3n), and (n,4n) reactions. + + :flux-YN: + Spherical harmonic expansion of the direction of motion + :math:`\left(\Omega\right)` of the total flux. This score will tally + all of the harmonic moments of order 0 to N. N must be between 0 and 10. + + :total-YN: + Spherical harmonic expansion of the incoming particle's direction of + motion :math:`\left(\Omega\right)` of the total flux. This score will + tally all of the harmonic moments of order 0 to N. N must be between 0 + and 10. + :current: Partial currents on the boundaries of each cell in a mesh. @@ -1171,7 +1195,7 @@ attributes or sub-elements. These are not used in "voxel" plots: Any number of this optional tag may be included in each ```` element, which can override the default random colors for cells or materials. Each ``col_spec`` element must contain ``id`` and ``rgb`` sub-elements. - + :id: Specifies the cell or material unique id for the color specification. @@ -1208,20 +1232,20 @@ attributes or sub-elements. These are not used in "voxel" plots: ------------------------------ CMFD Specification -- cmfd.xml ------------------------------ - + Coarse mesh finite difference acceleration method has been implemented in OpenMC. -Currently, it allows users to accelerate fission source convergence during -inactive neutron batches. To run CMFD, the ```` element in +Currently, it allows users to accelerate fission source convergence during +inactive neutron batches. To run CMFD, the ```` element in ``settings.xml`` should be set to "true". ```` Element -------------------------- The ```` element controls the batch where CMFD tallies should be -reset. CMFD tallies should be reset before active batches so they are accumulated +reset. CMFD tallies should be reset before active batches so they are accumulated without bias. - *Default*: 0 + *Default*: 0 ```` Element ------------------- @@ -1258,9 +1282,9 @@ It can be turned on with "true" and off with "false". ```` Element ---------------------- -The ```` element controls if cmfd tallies should be accumulated -during inactive batches. For some applications, CMFD tallies may not be -needed until the start of active batches. This option can be turned on +The ```` element controls if cmfd tallies should be accumulated +during inactive batches. For some applications, CMFD tallies may not be +needed until the start of active batches. This option can be turned on with "true" and off with "false" *Default*: true @@ -1277,12 +1301,12 @@ occurs. The amout of resets is controlled with the ```` element. ```` Element ------------------------- -The ```` element is used to view the convergence of linear GMRES -iterations in PETSc. This option can be turned on with "true" and turned off +The ```` element is used to view the convergence of linear GMRES +iterations in PETSc. This option can be turned on with "true" and turned off with "false". - *Default*: false + *Default*: false ```` Element ------------------ @@ -1295,7 +1319,7 @@ attributes/sub-elements: given, it is assumed that the mesh is an x-y mesh. :upper_right: - The upper-right corner of the structrued mesh. If only two coordinate are + The upper-right corner of the structrued mesh. If only two coordinate are given, it is assumed that the mesh is an x-y mesh. :dimension: @@ -1306,8 +1330,8 @@ attributes/sub-elements: :energy: Energy bins [in MeV], listed in ascending order (e.g. 0.0 0.625e-7 20.0) - for CMFD tallies and acceleration. If no energy bins are listed, OpenMC - automatically assumes a one energy group calculation over the entire + for CMFD tallies and acceleration. If no energy bins are listed, OpenMC + automatically assumes a one energy group calculation over the entire energy range. :albedo: @@ -1317,10 +1341,10 @@ attributes/sub-elements: *Default*: 1.0 1.0 1.0 1.0 1.0 1.0 :map: - An optional acceleration map can be specified to overlay on the coarse - mesh spatial grid. If this option is used a ``1`` is used for a + An optional acceleration map can be specified to overlay on the coarse + mesh spatial grid. If this option is used a ``1`` is used for a non-accelerated region and a ``2`` is used for an accelerated region. - For a simple 4x4 coarse mesh with a 2x2 fuel lattice surrounded by + For a simple 4x4 coarse mesh with a 2x2 fuel lattice surrounded by reflector, the map is: ``1 1 1 1`` @@ -1331,24 +1355,24 @@ attributes/sub-elements: ``1 1 1 1`` - Therefore a 2x2 system of equations is solved rather than a 4x4. This - is extremely important to use in reflectors as neutrons will not + Therefore a 2x2 system of equations is solved rather than a 4x4. This + is extremely important to use in reflectors as neutrons will not contribute to any tallies far away from fission source neutron regions. A ``2`` must be used to identify any fission source region. - .. note:: Only two of the following three sub-elements are needed: - ``lower_left``, ``upper_right`` and ``width``. Any combination + .. note:: Only two of the following three sub-elements are needed: + ``lower_left``, ``upper_right`` and ``width``. Any combination of two of these will yield the third. ```` Element ------------------ -The ```` element is used to normalize the CMFD fission source distribution -to a particular value. For example, if a fission source is calculated for a -17 x 17 lattice of pins, the fission source may be normalized to the number of -fission source regions, in this case 289. This is useful when visualizing this -distribution as the average peaking factor will be unity. This parameter will -not impact the calculation. +The ```` element is used to normalize the CMFD fission source distribution +to a particular value. For example, if a fission source is calculated for a +17 x 17 lattice of pins, the fission source may be normalized to the number of +fission source regions, in this case 289. This is useful when visualizing this +distribution as the average peaking factor will be unity. This parameter will +not impact the calculation. *Default*: 1.0 @@ -1363,7 +1387,7 @@ occur during inactive CMFD batches. ```` Element --------------------------- -The ```` element is used to view the convergence of power iteration. +The ```` element is used to view the convergence of power iteration. This option can be turned on with "true" and turned off with "false". *Default*: false @@ -1389,7 +1413,7 @@ function in PETSc. This option can be turned on with "true" and turned off with -------------------- The ```` element controls whether the CMFD eigenproblem is solved with -standard power iteration or nonlinear Jacobian-free Newton Krylov (JFNK). +standard power iteration or nonlinear Jacobian-free Newton Krylov (JFNK). By setting "power", power iteration is used and by setting "jfnk", JFNK is used. *Default*: power diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index ad02b542a3..7e6ca18e53 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -8,7 +8,7 @@ module cmfd_input implicit none private - public :: configure_cmfd + public :: configure_cmfd contains @@ -119,7 +119,7 @@ contains cmfd % egrid = (/0.0_8,20.0_8/) cmfd % indices(4) = 1 ! one energy group end if - + ! Set global albedo if (check_for_node(node_mesh, "albedo")) then call get_node_array(node_mesh, "albedo", cmfd % albedo) @@ -134,7 +134,7 @@ contains if (get_arraysize_integer(node_mesh, "map") /= & product(cmfd % indices(1:3))) then message = 'FATAL==>CMFD coremap not to correct dimensions' - call fatal_error() + call fatal_error() end if allocate(iarray(get_arraysize_integer(node_mesh, "map"))) call get_node_array(node_mesh, "map", iarray) @@ -206,7 +206,7 @@ contains ! Batch to begin cmfd if (check_for_node(doc, "begin")) & - call get_node_value(doc, "begin", cmfd_begin) + call get_node_value(doc, "begin", cmfd_begin) ! Tally during inactive batches if (check_for_node(doc, "inactive")) then @@ -288,7 +288,7 @@ contains m => meshes(n_user_meshes+1) ! Set mesh id - m % id = n_user_meshes + 1 + m % id = n_user_meshes + 1 ! Set mesh type to rectangular m % type = LATTICE_RECT @@ -422,7 +422,7 @@ contains if (check_for_node(node_mesh, "energy")) then n_filters = n_filters + 1 filters(n_filters) % type = FILTER_ENERGYIN - ng = get_arraysize_double(node_mesh, "energy") + ng = get_arraysize_double(node_mesh, "energy") filters(n_filters) % n_bins = ng - 1 allocate(filters(n_filters) % real_bins(ng)) call get_node_array(node_mesh, "energy", & @@ -454,20 +454,20 @@ contains allocate(t % filters(n_filters)) t % filters = filters(1:n_filters) - ! Allocate scoring bins + ! Allocate scoring bins allocate(t % score_bins(3)) t % n_score_bins = 3 t % n_user_score_bins = 3 ! Allocate scattering order data - allocate(t % scatt_order(3)) - t % scatt_order = 0 - + allocate(t % moment_order(3)) + t % moment_order = 0 + ! Set macro_bins t % score_bins(1) = SCORE_FLUX t % score_bins(2) = SCORE_TOTAL t % score_bins(3) = SCORE_SCATTER_N - t % scatt_order(3) = 1 + t % moment_order(3) = 1 else if (i == 2) then @@ -507,8 +507,8 @@ contains t % n_user_score_bins = 2 ! Allocate scattering order data - allocate(t % scatt_order(2)) - t % scatt_order = 0 + allocate(t % moment_order(2)) + t % moment_order = 0 ! Set macro_bins t % score_bins(1) = SCORE_NU_SCATTER @@ -550,8 +550,8 @@ contains t % n_user_score_bins = 1 ! Allocate scattering order data - allocate(t % scatt_order(1)) - t % scatt_order = 0 + allocate(t % moment_order(1)) + t % moment_order = 0 ! Set macro bins t % score_bins(1) = SCORE_CURRENT diff --git a/src/constants.F90 b/src/constants.F90 index db22f74cda..69954618d0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -11,7 +11,7 @@ module constants integer, parameter :: VERSION_RELEASE = 4 ! Revision numbers for binary files - integer, parameter :: REVISION_STATEPOINT = 11 + integer, parameter :: REVISION_STATEPOINT = 12 integer, parameter :: REVISION_PARTICLE_RESTART = 1 ! Binary file types @@ -276,7 +276,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 14 + integer, parameter :: N_SCORE_TYPES = 20 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -284,18 +284,37 @@ module constants SCORE_NU_SCATTER = -4, & ! scattering production rate SCORE_SCATTER_N = -5, & ! arbitrary scattering moment SCORE_SCATTER_PN = -6, & ! system for scoring 0th through nth moment - SCORE_TRANSPORT = -7, & ! transport reaction rate - SCORE_N_1N = -8, & ! (n,1n) rate - SCORE_ABSORPTION = -9, & ! absorption rate - SCORE_FISSION = -10, & ! fission rate - SCORE_NU_FISSION = -11, & ! neutron production rate - SCORE_KAPPA_FISSION = -12, & ! fission energy production rate - SCORE_CURRENT = -13, & ! partial current - SCORE_EVENTS = -14 ! number of events + SCORE_NU_SCATTER_N = -7, & ! arbitrary nu-scattering moment + SCORE_NU_SCATTER_PN = -8, & ! system for scoring 0th through nth nu-scatter moment + SCORE_TRANSPORT = -9, & ! transport reaction rate + SCORE_N_1N = -10, & ! (n,1n) rate + SCORE_ABSORPTION = -11, & ! absorption rate + SCORE_FISSION = -12, & ! fission rate + SCORE_NU_FISSION = -13, & ! neutron production rate + SCORE_KAPPA_FISSION = -14, & ! fission energy production rate + SCORE_CURRENT = -15, & ! partial current + SCORE_FLUX_YN = -16, & ! angular moment of flux + SCORE_TOTAL_YN = -17, & ! angular moment of total reaction rate + 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 ! Maximum scattering order supported - integer, parameter :: SCATT_ORDER_MAX = 10 - character(len=*), parameter :: SCATT_ORDER_MAX_PNSTR = "scatter-p10" + integer, parameter :: MAX_ANG_ORDER = 10 + + ! Names of *-PN & *-YN scores (MOMENT_STRS) and *-N moment scores + character(*), parameter :: & + MOMENT_STRS(6) = (/ "scatter-p ", & + "nu-scatter-p", & + "flux-y ", & + "total-y ", & + "scatter-y ", & + "nu-scatter-y"/), & + MOMENT_N_STRS(2) = (/ "scatter- ", & + "nu-scatter- "/) + + ! Location in MOMENT_STRS where the YN data begins + integer, parameter :: YN_LOC = 3 ! Tally map bin finding integer, parameter :: NO_BIN_FOUND = -1 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2d88871db6..46461b78da 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -185,7 +185,7 @@ contains ! If the number of particles was specified as a command-line argument, we ! don't set it here - if (n_particles == 0) n_particles = temp_long + if (n_particles == 0) n_particles = temp_long ! Copy batch information call get_node_value(node_mode, "batches", n_batches) @@ -278,10 +278,10 @@ contains else ! Spatial distribution for external source - if (check_for_node(node_source, "space")) then + if (check_for_node(node_source, "space")) then ! Get pointer to spatial distribution - call get_node_ptr(node_source, "space", node_dist) + call get_node_ptr(node_source, "space", node_dist) ! Check for type of spatial distribution type = '' @@ -624,7 +624,7 @@ contains elseif (check_for_node(node_sp, "interval")) then ! User gave an interval for writing state points call get_node_value(node_sp, "interval", temp_int) - n_state_points = n_batches / temp_int + n_state_points = n_batches / temp_int do i = 1, n_state_points call statepoint_batch % add(temp_int * i) end do @@ -648,7 +648,7 @@ contains ! Determine number of batches at which to store source points if (check_for_node(node_sp, "batches")) then - n_source_points = get_arraysize_integer(node_sp, "batches") + n_source_points = get_arraysize_integer(node_sp, "batches") else n_source_points = 0 end if @@ -700,7 +700,7 @@ contains end if else ! If no tag was present, by default we keep source bank in - ! statepoint file and write it out at statepoints intervals + ! statepoint file and write it out at statepoints intervals source_separate = .false. n_source_points = n_state_points do i = 1, n_state_points @@ -944,7 +944,7 @@ contains ! Check to make sure that either material or fill was specified if (c % material == NONE .and. c % fill == NONE) then - message = "Neither material nor fill was specified for cell " // & + message = "Neither material nor fill was specified for cell " // & trim(to_str(c % id)) call fatal_error() end if @@ -1139,7 +1139,7 @@ contains n = get_arraysize_double(node_surf, "coeffs") if (n < coeffs_reqd) then - message = "Not enough coefficients specified for surface: " // & + message = "Not enough coefficients specified for surface: " // & trim(to_str(s % id)) call fatal_error() elseif (n > coeffs_reqd) then @@ -1615,7 +1615,7 @@ contains ! Check to make sure cross-section is continuous energy neutron table n = len_trim(name) if (name(n:n) /= 'c') then - message = "Cross-section table " // trim(name) // & + message = "Cross-section table " // trim(name) // & " is not a continuous-energy neutron table." call fatal_error() end if @@ -1644,7 +1644,7 @@ contains ! Check to make sure either all atom percents or all weight percents are ! given - if (.not. (all(mat % atom_density > ZERO) .or. & + if (.not. (all(mat % atom_density > ZERO) .or. & all(mat % atom_density < ZERO))) then message = "Cannot mix atom and weight percents in material " // & to_str(mat % id) @@ -1745,12 +1745,14 @@ contains integer :: n ! size of arrays in mesh specification integer :: n_words ! number of words read integer :: n_filters ! number of filters - integer :: n_new ! number of new scores to add based on Pn tally - integer :: n_scores ! number of tot scores after adjusting for Pn tally - integer :: n_order ! Scattering order requested + integer :: n_new ! number of new scores to add based on Yn/Pn tally + integer :: n_scores ! number of tot scores after adjusting for Yn/Pn tally + integer :: n_bins ! Total new bins for this score + integer :: n_order ! Moment order requested integer :: n_order_pos ! Position of Scattering order in score name string integer :: MT ! user-specified MT for score integer :: iarray3(3) ! temporary integer array + integer :: imomstr ! Index of MOMENT_STRS & MOMENT_N_STRS logical :: file_exists ! does tallies.xml file exist? real(8) :: rarray3(3) ! temporary double prec. array character(MAX_LINE_LEN) :: filename @@ -2168,7 +2170,7 @@ contains case default ! Specified tally filter is invalid, raise error - message = "Unknown filter type '" // & + message = "Unknown filter type '" // & trim(temp_str) // "' on tally " // & trim(to_str(t % id)) // "." call fatal_error() @@ -2220,7 +2222,7 @@ contains t % all_nuclides = .true. else ! Any other case, e.g. U-235 Pu-239 - n_words = get_arraysize_string(node_tal, "nuclides") + n_words = get_arraysize_string(node_tal, "nuclides") allocate(t % nuclide_bins(n_words)) do j = 1, n_words ! Check if total material was specified @@ -2242,7 +2244,7 @@ contains word = pair_list % key(1:150) exit end if - + ! Advance to next pair_list => pair_list % next end do @@ -2291,79 +2293,110 @@ contains ! READ DATA FOR SCORES if (check_for_node(node_tal, "scores")) then - ! Loop through scores and determine if a scatter-p# input was used - ! to allow for proper pre-allocating of t % score_bins - ! This scheme allows multiple scatter-p# to be requested by the user - ! if so desired n_words = get_arraysize_string(node_tal, "scores") allocate(sarray(n_words)) call get_node_array(node_tal, "scores", sarray) + + ! Before we can allocate storage for scores, we must determine the + ! number of additional scores required due to the moment scores + ! (i.e., scatter-p#, flux-y#) n_new = 0 do j = 1, n_words call lower_case(sarray(j)) - ! Find if scores(j) is of the form 'scatter-p' - ! If so, get the number and do a select case on that. + ! Find if scores(j) is of the form 'moment-p' or 'moment-y' present in + ! MOMENT_STRS(:) + ! If so, check the order, store if OK, then reset the number to 'n' score_name = trim(sarray(j)) - if (starts_with(score_name,'scatter-p')) then - n_order_pos = scan(score_name,'0123456789') - n_order = int(str_to_int( & - score_name(n_order_pos:(len_trim(score_name)))),4) - if (n_order > SCATT_ORDER_MAX) then - ! Throw a warning. Set to the maximum number. - ! The above scheme will essentially take the absolute value - message = "Invalid scattering order of " // trim(to_str(n_order)) // & - " requested. Setting to the maximum permissible value, " // & - trim(to_str(SCATT_ORDER_MAX)) - call warning() - n_order = SCATT_ORDER_MAX - sarray(j) = SCATT_ORDER_MAX_PNSTR + do imomstr = 1, size(MOMENT_STRS) + if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then + n_order_pos = scan(score_name,'0123456789') + n_order = int(str_to_int( & + score_name(n_order_pos:(len_trim(score_name)))),4) + if (n_order > MAX_ANG_ORDER) then + ! User requested too many orders; throw a warning and set to the + ! maximum order. + ! The above scheme will essentially take the absolute value + message = "Invalid scattering order of " // trim(to_str(n_order)) // & + " requested. Setting to the maximum permissible value, " // & + trim(to_str(MAX_ANG_ORDER)) + call warning() + n_order = MAX_ANG_ORDER + sarray(j) = trim(MOMENT_STRS(imomstr)) // & + trim(to_str(MAX_ANG_ORDER)) + end if + ! Find total number of bins for this case + if (imomstr >= YN_LOC) then + n_bins = (n_order + 1)**2 + else + n_bins = n_order + 1 + end if + ! We subtract one since n_words already included + n_new = n_new + n_bins - 1 + exit end if - n_new = n_new + n_order - end if + end do end do n_scores = n_words + n_new - - ! Allocate accordingly + + ! Allocate score storage accordingly allocate(t % score_bins(n_scores)) - allocate(t % scatt_order(n_scores)) - t % scatt_order = 0 + allocate(t % moment_order(n_scores)) + t % moment_order = 0 j = 0 do l = 1, n_words j = j + 1 - ! Get the input string in scores(l) but if scatter-n or scatter-pn - ! then strip off the n, and store it as an integer to be used later - ! Peform the select case on this modified (number removed) string + ! Get the input string in scores(l) but if score is one of the moment + ! scores then strip off the n and store it as an integer to be used + ! later. Then perform the select case on this modified (number removed) string score_name = sarray(l) - if (starts_with(score_name,'scatter-p')) then - n_order_pos = scan(score_name,'0123456789') - n_order = int(str_to_int( & - score_name(n_order_pos:(len_trim(score_name)))),4) - if (n_order > SCATT_ORDER_MAX) then - ! Throw a warning. Set to the maximum number. - ! The above scheme will essentially take the absolute value - message = "Invalid scattering order of " // trim(to_str(n_order)) // & - " requested. Setting to the maximum permissible value, " // & - trim(to_str(SCATT_ORDER_MAX)) - call warning() - n_order = SCATT_ORDER_MAX + do imomstr = 1, size(MOMENT_STRS) + if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then + n_order_pos = scan(score_name,'0123456789') + n_order = int(str_to_int( & + score_name(n_order_pos:(len_trim(score_name)))),4) + if (n_order > MAX_ANG_ORDER) then + ! User requested too many orders; throw a warning and set to the + ! maximum order. + ! The above scheme will essentially take the absolute value + message = "Invalid scattering order of " // trim(to_str(n_order)) // & + " requested. Setting to the maximum permissible value, " // & + trim(to_str(MAX_ANG_ORDER)) + call warning() + n_order = MAX_ANG_ORDER + end if + score_name = trim(MOMENT_STRS(imomstr)) // "n" + ! Find total number of bins for this case + if (imomstr >= YN_LOC) then + n_bins = (n_order + 1)**2 + else + n_bins = n_order + 1 + end if + exit end if - score_name = "scatter-pn" - else if (starts_with(score_name,'scatter-')) then - n_order_pos = scan(score_name,'0123456789') - n_order = int(str_to_int( & - score_name(n_order_pos:(len_trim(score_name)))),4) - if (n_order > SCATT_ORDER_MAX) then - ! Throw a warning. Set to the maximum number. - ! The above scheme will essentially take the absolute value - message = "Invalid scattering order of " // trim(to_str(n_order)) // & - " requested. Setting to the maximum permissible value, " // & - trim(to_str(SCATT_ORDER_MAX)) - call warning() - n_order = SCATT_ORDER_MAX - end if - score_name = "scatter-n" + end do + ! Now check the Moment_N_Strs, but only if we werent successful above + if (imomstr > size(MOMENT_STRS)) then + do imomstr = 1, size(MOMENT_N_STRS) + if (starts_with(score_name,trim(MOMENT_N_STRS(imomstr)))) then + n_order_pos = scan(score_name,'0123456789') + n_order = int(str_to_int( & + score_name(n_order_pos:(len_trim(score_name)))),4) + if (n_order > MAX_ANG_ORDER) then + ! User requested too many orders; throw a warning and set to the + ! maximum order. + ! The above scheme will essentially take the absolute value + message = "Invalid scattering order of " // trim(to_str(n_order)) // & + " requested. Setting to the maximum permissible value, " // & + trim(to_str(MAX_ANG_ORDER)) + call warning() + n_order = MAX_ANG_ORDER + end if + score_name = trim(MOMENT_N_STRS(imomstr)) // "n" + exit + end if + end do end if - + select case (trim(score_name)) case ('flux') ! Prohibit user from tallying flux for an individual nuclide @@ -2378,6 +2411,23 @@ contains message = "Cannot tally flux with an outgoing energy filter." call fatal_error() end if + case ('flux-yn') + ! Prohibit user from tallying flux for an individual nuclide + if (.not. (t % n_nuclide_bins == 1 .and. & + t % nuclide_bins(1) == -1)) then + message = "Cannot tally flux for an individual nuclide." + call fatal_error() + end if + + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + message = "Cannot tally flux with an outgoing energy filter." + call fatal_error() + end if + + t % score_bins(j : j + n_bins - 1) = SCORE_FLUX_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 + case ('total') t % score_bins(j) = SCORE_TOTAL if (t % find_filter(FILTER_ENERGYOUT) > 0) then @@ -2385,9 +2435,21 @@ contains &outgoing energy filter." call fatal_error() end if + + case ('total-yn') + if (t % find_filter(FILTER_ENERGYOUT) > 0) then + message = "Cannot tally total reaction rate with an & + &outgoing energy filter." + call fatal_error() + end if + + t % score_bins(j : j + n_bins - 1) = SCORE_TOTAL_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 + case ('scatter') t % score_bins(j) = SCORE_SCATTER - + case ('nu-scatter') t % score_bins(j) = SCORE_NU_SCATTER @@ -2401,22 +2463,53 @@ contains ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG end if - t % scatt_order(j) = n_order - + t % moment_order(j) = n_order + + case ('nu-scatter-n') + ! Set tally estimator to analog + t % estimator = ESTIMATOR_ANALOG + if (n_order == 0) then + t % score_bins(j) = SCORE_NU_SCATTER + else + t % score_bins(j) = SCORE_NU_SCATTER_N + end if + t % moment_order(j) = n_order + case ('scatter-pn') t % estimator = ESTIMATOR_ANALOG ! Setup P0:Pn - t % score_bins(j : j + n_order) = SCORE_SCATTER_PN - t % scatt_order(j : j + n_order) = n_order - j = j + n_order - + t % score_bins(j : j + n_bins - 1) = SCORE_SCATTER_PN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 + + case ('nu-scatter-pn') + t % estimator = ESTIMATOR_ANALOG + ! Setup P0:Pn + t % score_bins(j : j + n_bins - 1) = SCORE_NU_SCATTER_PN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 + + case ('scatter-yn') + t % estimator = ESTIMATOR_ANALOG + ! Setup P0:Pn + t % score_bins(j : j + n_bins - 1) = SCORE_SCATTER_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 + + case ('nu-scatter-yn') + t % estimator = ESTIMATOR_ANALOG + ! Setup P0:Pn + t % score_bins(j : j + n_bins - 1) = SCORE_NU_SCATTER_YN + t % moment_order(j : j + n_bins - 1) = n_order + j = j + n_bins - 1 + case('transport') t % score_bins(j) = SCORE_TRANSPORT ! Set tally estimator to analog t % estimator = ESTIMATOR_ANALOG case ('diffusion') - message = "Diffusion score no longer supported for tallies, & + message = "Diffusion score no longer supported for tallies, & &please remove" call fatal_error() case ('n1n') @@ -2519,14 +2612,14 @@ contains t % score_bins(j) = MT else message = "Invalid MT on : " // & - trim(sarray(j)) + trim(sarray(l)) call fatal_error() end if else ! Specified score was not an integer message = "Unknown scoring function: " // & - trim(sarray(j)) + trim(sarray(l)) call fatal_error() end if @@ -2675,7 +2768,7 @@ contains pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & "_" // trim(filename) // ".voxel" end select - + ! Copy plot pixel size if (pl % type == PLOT_TYPE_SLICE) then if (get_arraysize_integer(node_plot, "pixels") == 2) then @@ -2698,7 +2791,7 @@ contains ! Copy plot background color if (check_for_node(node_plot, "background")) then if (pl % type == PLOT_TYPE_VOXEL) then - message = "Background color ignored in voxel plot " // & + message = "Background color ignored in voxel plot " // & trim(to_str(pl % id)) call warning() end if @@ -2712,7 +2805,7 @@ contains else pl % not_found % rgb = (/ 255, 255, 255 /) end if - + ! Copy plot basis if (pl % type == PLOT_TYPE_SLICE) then temp_str = 'xy' @@ -2727,12 +2820,12 @@ contains case ("yz") pl % basis = PLOT_BASIS_YZ case default - message = "Unsupported plot basis '" // trim(temp_str) & + message = "Unsupported plot basis '" // trim(temp_str) & // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select end if - + ! Copy plotting origin if (get_arraysize_double(node_plot, "origin") == 3) then call get_node_array(node_plot, "origin", pl % origin) @@ -2799,13 +2892,13 @@ contains ! Copy user specified colors if (n_cols /= 0) then - + if (pl % type == PLOT_TYPE_VOXEL) then - message = "Color specifications ignored in voxel plot " // & + message = "Color specifications ignored in voxel plot " // & trim(to_str(pl % id)) call warning() end if - + do j = 1, n_cols ! Get pointer to color spec XML node @@ -2815,7 +2908,7 @@ contains if (get_arraysize_double(node_col, "rgb") /= 3) then message = "Bad RGB " & // "in plot " // trim(to_str(pl % id)) - call fatal_error() + call fatal_error() end if ! Ensure that there is an id for this color specification @@ -2858,13 +2951,13 @@ contains call get_node_list(node_plot, "mask", node_mask_list) n_masks = get_list_size(node_mask_list) if (n_masks /= 0) then - + if (pl % type == PLOT_TYPE_VOXEL) then - message = "Mask ignored in voxel plot " // & + message = "Mask ignored in voxel plot " // & trim(to_str(pl % id)) call warning() end if - + select case(n_masks) case default message = "Mutliple masks" // & @@ -2885,14 +2978,14 @@ contains end if allocate(iarray(n_comp)) call get_node_array(node_mask, "components", iarray) - + ! First we need to change the user-specified identifiers to indices ! in the cell and material arrays do j=1, n_comp col_id = iarray(j) - + if (pl % color_by == PLOT_COLOR_CELLS) then - + if (cell_dict % has_key(col_id)) then iarray(j) = cell_dict % get_key(col_id) else @@ -2900,9 +2993,9 @@ contains " specified in the mask in plot " // trim(to_str(pl % id)) call fatal_error() end if - + else if (pl % color_by == PLOT_COLOR_MATS) then - + if (material_dict % has_key(col_id)) then iarray(j) = material_dict % get_key(col_id) else @@ -2910,10 +3003,10 @@ contains " specified in the mask in plot " // trim(to_str(pl % id)) call fatal_error() end if - - end if + + end if end do - + ! Alter colors based on mask information do j=1,size(pl % colors) if (.not. any(j .eq. iarray)) then @@ -2928,9 +3021,9 @@ contains end do deallocate(iarray) - + end select - + end if ! Add plot to dictionary @@ -2970,7 +3063,7 @@ contains "' does not exist!" call fatal_error() end if - + message = "Reading cross sections XML file..." call write_message(5) diff --git a/src/math.F90 b/src/math.F90 index 61bb8e5f35..b3046ccc4a 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -1,6 +1,6 @@ module math - use constants, only: PI, ONE, TWO + use constants, only: PI, ONE, TWO, ZERO use random_lcg, only: prn implicit none @@ -116,20 +116,20 @@ contains !=============================================================================== ! CALC_PN calculates the n-th order Legendre polynomial at the value of x. ! Since this function is called repeatedly during the neutron transport process, -! neither n or x is checked to see if they are in the applicable range. +! neither n or x is checked to see if they are in the applicable range. ! This is left to the client developer to use where applicable. x is to be in ! the domain of [-1,1], and 0<=n<=5. If x is outside of the range, the return -! value will be outside the expected range; if n is outside the stated range, +! value will be outside the expected range; if n is outside the stated range, ! the return value will be 1.0. !=============================================================================== - + pure function calc_pn(n,x) result(pnx) integer, intent(in) :: n ! Legendre order requested - real(8), intent(in) :: x ! Independent variable the Legendre is to be + real(8), intent(in) :: x ! Independent variable the Legendre is to be ! evaluated at; x must be in the domain [-1,1] real(8) :: pnx ! The Legendre poly of order n evaluated at x - + select case(n) case(1) pnx = x @@ -144,7 +144,7 @@ contains case(6) pnx = 14.4375_8 * (x ** 6) - 19.6875_8 * (x ** 4) + & 6.5625_8 * x * x - 0.3125_8 - case(7) + case(7) pnx = 26.8125_8 * (x ** 7) - 43.3125_8 * (x ** 5) + & 19.6875_8 * x * x * x - 2.1875_8 * x case(8) @@ -160,9 +160,393 @@ contains case default pnx = ONE ! correct for case(0), incorrect for the rest end select - + end function calc_pn +!=============================================================================== +! CALC_RN calculates the n-th order spherical harmonics for a given angle +! (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n) +!=============================================================================== + + pure function calc_rn(n,uvw) result(rn) + + integer, intent(in) :: n ! Order requested + real(8), intent(in) :: uvw(3) ! Direction of travel, assumed to be on unit sphere + real(8) :: rn(2*n + 1) ! The resultant R_n(uvw) + + real(8) :: phi, w ! Azimuthal and Cosine of Polar angles (from uvw) + real(8) :: w2m1 ! (w^2 - 1), frequently used in these + + w = uvw(3) ! z = cos(polar) + if (uvw(1) == ZERO) then + phi = ZERO + else +! phi = atan(uvw(2) / uvw(1)) + phi = atan2(uvw(2), uvw(1)) + end if + + w2m1 = (ONE - w**2) + select case(n) + case (0) + ! l = 0, m = 0 + rn(1) = ONE + case (1) + ! l = 1, m = -1 + rn(1) = ONE*sqrt(w2m1) * sin(phi) + ! l = 1, m = 0 + rn(2) = ONE * w + ! l = 1, m = 1 + rn(3) = ONE*sqrt(w2m1) * cos(phi) + case (2) + ! l = 2, m = -2 + rn(1) = 0.288675134594813_8 * (-3.0_8 * w**2 + 3.0_8) * sin(TWO*phi) + ! l = 2, m = -1 + rn(2) = 1.73205080756888_8 * w*sqrt(w2m1) * sin(phi) + ! l = 2, m = 0 + rn(3) = 1.5_8 * w**2 - 0.5_8 + ! l = 2, m = 1 + rn(4) = 1.73205080756888_8 * w*sqrt(w2m1) * cos(phi) + ! l = 2, m = 2 + rn(5) = 0.288675134594813_8 * (-3.0_8 * w**2 + 3.0_8) * cos(TWO*phi) + case (3) + ! l = 3, m = -3 + rn(1) = 0.790569415042095_8 * (w2m1)**(3.0_8/TWO) * sin(3.0_8 * phi) + ! l = 3, m = -2 + rn(2) = 1.93649167310371_8 * w*(w2m1) * sin(TWO*phi) + ! l = 3, m = -1 + rn(3) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - 3.0_8/TWO) * & + sin(phi) + ! l = 3, m = 0 + rn(4) = 2.5_8 * w**3 - 1.5_8 * w + ! l = 3, m = 1 + rn(5) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - 3.0_8/TWO) * & + cos(phi) + ! l = 3, m = 2 + rn(6) = 1.93649167310371_8 * w*(w2m1) * cos(TWO*phi) + ! l = 3, m = 3 + rn(7) = 0.790569415042095_8 * (w2m1)**(3.0_8/TWO) * cos(3.0_8* phi) + case (4) + ! l = 4, m = -4 + rn(1) = 0.739509972887452_8 * (w2m1)**2 * sin(4.0_8*phi) + ! l = 4, m = -3 + rn(2) = 2.09165006633519_8 * w*(w2m1)**(3.0_8/TWO) * sin(3.0_8* phi) + ! l = 4, m = -2 + rn(3) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * & + sin(TWO*phi) + ! l = 4, m = -1 + rn(4) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * & + sin(phi) + ! l = 4, m = 0 + rn(5) = 4.375_8 * w**4 - 3.75_8 * w**2 + 0.375_8 + ! l = 4, m = 1 + rn(6) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * & + cos(phi) + ! l = 4, m = 2 + rn(7) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * & + cos(TWO*phi) + ! l = 4, m = 3 + rn(8) = 2.09165006633519_8 * w*(w2m1)**(3.0_8/TWO) * cos(3.0_8* phi) + ! l = 4, m = 4 + rn(9) = 0.739509972887452_8 * (w2m1)**2 * cos(4.0_8*phi) + case (5) + ! l = 5, m = -5 + rn(1) = 0.701560760020114_8 * (w2m1)**(5.0_8/TWO) * sin(5.0_8*phi) + ! l = 5, m = -4 + rn(2) = 2.21852991866236_8 * w*(w2m1)**2 * sin(4.0_8*phi) + ! l = 5, m = -3 + rn(3) = 0.00996023841111995_8 * (w2m1)**(3.0_8/TWO)* & + ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(3.0_8*phi) + ! l = 5, m = -2 + rn(4) = 0.0487950036474267_8 * (w2m1)*((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * & + sin(TWO*phi) + ! l = 5, m = -1 + rn(5) = 0.258198889747161_8*sqrt(w2m1)* & + ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * sin(phi) + ! l = 5, m = 0 + rn(6) = 7.875_8 * w**5 - 8.75_8 * w**3 + 1.875_8 * w + ! l = 5, m = 1 + rn(7) = 0.258198889747161_8*sqrt(w2m1)* & + ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * cos(phi) + ! l = 5, m = 2 + rn(8) = 0.0487950036474267_8 * (w2m1)* & + ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi) + ! l = 5, m = 3 + rn(9) = 0.00996023841111995_8 * (w2m1)**(3.0_8/TWO)* & + ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(3.0_8*phi) + ! l = 5, m = 4 + rn(10) = 2.21852991866236_8 * w*(w2m1)**2 * cos(4.0_8*phi) + ! l = 5, m = 5 + rn(11) = 0.701560760020114_8 * (w2m1)**(5.0_8/TWO) * cos(5.0_8* phi) + case (6) + ! l = 6, m = -6 + rn(1) = 0.671693289381396_8 * (w2m1)**3 * sin(6.0_8*phi) + ! l = 6, m = -5 + rn(2) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * sin(5.0_8*phi) + ! l = 6, m = -4 + rn(3) = 0.00104990131391452_8 * (w2m1)**2 * & + ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi) + ! l = 6, m = -3 + rn(4) = 0.00575054632785295_8 * (w2m1)**(3.0_8/TWO) * & + ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(3.0_8*phi) + ! l = 6, m = -2 + rn(5) = 0.0345032779671177_8 * (w2m1) * & + ((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * sin(TWO*phi) + ! l = 6, m = -1 + rn(6) = 0.218217890235992_8*sqrt(w2m1) * & + ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * sin(phi) + ! l = 6, m = 0 + rn(7) = 14.4375_8 * w**6 - 19.6875_8 * w**4 + 6.5625_8 * w**2 - 0.3125_8 + ! l = 6, m = 1 + rn(8) = 0.218217890235992_8*sqrt(w2m1) * & + ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * cos(phi) + ! l = 6, m = 2 + rn(9) = 0.0345032779671177_8 * (w2m1) * & + ((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * cos(TWO*phi) + ! l = 6, m = 3 + rn(10) = 0.00575054632785295_8 * (w2m1)**(3.0_8/TWO) * & + ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(3.0_8*phi) + ! l = 6, m = 4 + rn(11) = 0.00104990131391452_8 * (w2m1)**2 * & + ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi) + ! l = 6, m = 5 + rn(12) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * cos(5.0_8*phi) + ! l = 6, m = 6 + rn(13) = 0.671693289381396_8 * (w2m1)**3 * cos(6.0_8*phi) + case (7) + ! l = 7, m = -7 + rn(1) = 0.647259849287749_8 * (w2m1)**(7.0_8/TWO) * sin(7.0_8*phi) + ! l = 7, m = -6 + rn(2) = 2.42182459624969_8 * w*(w2m1)**3 * sin(6.0_8*phi) + ! l = 7, m = -5 + rn(3) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* & + ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi) + ! l = 7, m = -4 + rn(4) = 0.000548293079133141_8 * (w2m1)**2* & + ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi) + ! l = 7, m = -3 + rn(5) = 0.00363696483726654_8 * (w2m1)**(3.0_8/TWO)* & + ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & + sin(3.0_8*phi) + ! l = 7, m = -2 + rn(6) = 0.025717224993682_8 * (w2m1)* & + ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & + sin(TWO*phi) + ! l = 7, m = -1 + rn(7) = 0.188982236504614_8*sqrt(w2m1)* & + ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & + (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi) + ! l = 7, m = 0 + rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 * w + ! l = 7, m = 1 + rn(9) = 0.188982236504614_8*sqrt(w2m1)* & + ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & + (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi) + ! l = 7, m = 2 + rn(10) = 0.025717224993682_8 * (w2m1)* & + ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & + cos(TWO*phi) + ! l = 7, m = 3 + rn(11) = 0.00363696483726654_8 * (w2m1)**(3.0_8/TWO)* & + ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & + cos(3.0_8*phi) + ! l = 7, m = 4 + rn(12) = 0.000548293079133141_8 * (w2m1)**2 * & + ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi) + ! l = 7, m = 5 + rn(13) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* & + ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi) + ! l = 7, m = 6 + rn(14) = 2.42182459624969_8 * w*(w2m1)**3 * cos(6.0_8*phi) + ! l = 7, m = 7 + rn(15) = 0.647259849287749_8 * (w2m1)**(7.0_8/TWO) * cos(7.0_8*phi) + case (8) + ! l = 8, m = -8 + rn(1) = 0.626706654240044_8 * (w2m1)**4 * sin(8.0_8*phi) + ! l = 8, m = -7 + rn(2) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * sin(7.0_8*phi) + ! l = 8, m = -6 + rn(3) = 6.77369783729086d-6*(w2m1)**3* & + ((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi) + ! l = 8, m = -5 + rn(4) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)* & + ((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi) + ! l = 8, m = -4 + rn(5) = 0.000316557156832328_8 * (w2m1)**2* & + ((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * sin(4.0_8*phi) + ! l = 8, m = -3 + rn(6) = 0.00245204119306875_8 * (w2m1)**(3.0_8/TWO)* & + ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + (10395.0_8/8.0_8)*w) * sin(3.0_8*phi) + ! l = 8, m = -2 + rn(7) = 0.0199204768222399_8 * (w2m1)* & + ((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + & + (10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi) + ! l = 8, m = -1 + rn(8) = 0.166666666666667_8*sqrt(w2m1)* & + ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & + (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi) + ! l = 8, m = 0 + rn(9) = 50.2734375_8 * w**8 - 93.84375_8 * w**6 + 54.140625_8 * w**4 -& + 9.84375_8 * w**2 + 0.2734375_8 + ! l = 8, m = 1 + rn(10) = 0.166666666666667_8*sqrt(w2m1)* & + ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & + (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi) + ! l = 8, m = 2 + rn(11) = 0.0199204768222399_8 * (w2m1)*((45045.0_8/16.0_8)*w**6- & + 45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - & + 315.0_8/16.0_8) * cos(TWO*phi) + ! l = 8, m = 3 + rn(12) = 0.00245204119306875_8 * (w2m1)**(3.0_8/TWO)* & + ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + & + (10395.0_8/8.0_8)*w) * cos(3.0_8*phi) + ! l = 8, m = 4 + rn(13) = 0.000316557156832328_8 * (w2m1)**2*((675675.0_8/8.0_8)*w**4 - & + 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi) + ! l = 8, m = 5 + rn(14) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 - & + 135135.0_8/TWO*w) * cos(5.0_8*phi) + ! l = 8, m = 6 + rn(15) = 6.77369783729086d-6*(w2m1)**3*((2027025.0_8/TWO)*w**2 - & + 135135.0_8/TWO) * cos(6.0_8*phi) + ! l = 8, m = 7 + rn(16) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * cos(7.0_8*phi) + ! l = 8, m = 8 + rn(17) = 0.626706654240044_8 * (w2m1)**4 * cos(8.0_8*phi) + case (9) + ! l = 9, m = -9 + rn(1) = 0.609049392175524_8 * (w2m1)**(9.0_8/TWO) * sin(9.0_8*phi) + ! l = 9, m = -8 + rn(2) = 2.58397773170915_8 * w*(w2m1)**4 * sin(8.0_8*phi) + ! l = 9, m = -7 + rn(3) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* & + ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi) + ! l = 9, m = -6 + rn(4) = 3.02928976464514d-6*(w2m1)**3* & + ((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi) + ! l = 9, m = -5 + rn(5) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)* & + ((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + & + 135135.0_8/8.0_8) * sin(5.0_8*phi) + ! l = 9, m = -4 + rn(6) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - & + 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi) + ! l = 9, m = -3 + rn(7) = 0.00173385495536766_8 * (w2m1)**(3.0_8/TWO)* & + ((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + & + (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(3.0_8*phi) + ! l = 9, m = -2 + rn(8) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7- & + 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/16.0_8 * w)* & + sin(TWO*phi) + ! l = 9, m = -1 + rn(9) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - & + 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * sin(phi) + ! l = 9, m = 0 + rn(10) = 94.9609375_8 * w**9 - 201.09375_8 * w**7 + 140.765625_8 * w**5- & + 36.09375_8 * w**3 + 2.4609375_8 * w + ! l = 9, m = 1 + rn(11) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - & + 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * cos(phi) + ! l = 9, m = 2 + rn(12) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7 - & + 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/ 16.0_8 * w) * & + cos(TWO*phi) + ! l = 9, m = 3 + rn(13) = 0.00173385495536766_8 * (w2m1)**(3.0_8/TWO)*((765765.0_8/16.0_8)*w**6 - & + 675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8)* & + cos(3.0_8*phi) + ! l = 9, m = 4 + rn(14) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - & + 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi) + ! l = 9, m = 5 + rn(15) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)*((11486475.0_8/8.0_8)* & + w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi) + ! l = 9, m = 6 + rn(16) = 3.02928976464514d-6*(w2m1)**3*((11486475.0_8/TWO)*w**3 - & + 2027025.0_8/TWO*w) * cos(6.0_8*phi) + ! l = 9, m = 7 + rn(17) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* & + ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi) + ! l = 9, m = 8 + rn(18) = 2.58397773170915_8 * w*(w2m1)**4 * cos(8.0_8*phi) + ! l = 9, m = 9 + rn(19) = 0.609049392175524_8 * (w2m1)**(9.0_8/TWO) * cos(9.0_8*phi) + case (10) + ! l = 10, m = -10 + rn(1) = 0.593627917136573_8 * (w2m1)**5 * sin(10.0_8*phi) + ! l = 10, m = -9 + rn(2) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * sin(9.0_8*phi) + ! l = 10, m = -8 + rn(3) = 2.49953651452314d-8*(w2m1)**4*((654729075.0_8/TWO)*w**2 - & + 34459425.0_8/TWO) * sin(8.0_8*phi) + ! l = 10, m = -7 + rn(4) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* & + ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi) + ! l = 10, m = -6 + rn(5) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - & + 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi) + ! l = 10, m = -5 + rn(6) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* & + ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & + (2027025.0_8/8.0_8)*w) * sin(5.0_8*phi) + ! l = 10, m = -4 + rn(7) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - & + 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & + 45045.0_8/16.0_8) * sin(4.0_8*phi) + ! l = 10, m = -3 + rn(8) = 0.00127230170115096_8 * (w2m1)**(3.0_8/TWO)* & + ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & + (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(3.0_8*phi) + ! l = 10, m = -2 + rn(9) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - & + 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - & + 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi) + ! l = 10, m = -1 + rn(10) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - & + 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - & + 15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi) + ! l = 10, m = 0 + rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 * w**6 - & + 117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8 + ! l = 10, m = 1 + rn(12) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - & + 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ & + 32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi) + ! l = 10, m = 2 + rn(13) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - & + 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -& + 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi) + ! l = 10, m = 3 + rn(14) = 0.00127230170115096_8 * (w2m1)**(3.0_8/TWO)* & + ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & + (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(3.0_8*phi) + ! l = 10, m = 4 + rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - & + 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & + 45045.0_8/16.0_8) * cos(4.0_8*phi) + ! l = 10, m = 5 + rn(16) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* & + ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & + (2027025.0_8/8.0_8)*w) * cos(5.0_8*phi) + ! l = 10, m = 6 + rn(17) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - & + 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi) + ! l = 10, m = 7 + rn(18) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* & + ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi) + ! l = 10, m = 8 + rn(19) = 2.49953651452314d-8*(w2m1)**4* & + ((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi) + ! l = 10, m = 9 + rn(20) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * cos(9.0_8*phi) + ! l = 10, m = 10 + rn(21) = 0.593627917136573_8 * (w2m1)**5 * cos(10.0_8*phi) + case default + rn = ONE + end select + + end function calc_rn + !=============================================================================== ! MAXWELL_SPECTRUM samples an energy from the Maxwell fission distribution based ! on a direct sampling scheme. The probability distribution function for a diff --git a/src/output.F90 b/src/output.F90 index fab3812631..834c28e188 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -672,7 +672,7 @@ contains integer :: j ! index in filters array integer :: id ! user-specified id integer :: unit_ ! unit to write to - integer :: n ! scattering order to include in name + integer :: n ! moment order to include in name character(MAX_LINE_LEN) :: string character(MAX_WORD_LEN) :: pn_string type(Cell), pointer :: c => null() @@ -825,20 +825,63 @@ contains select case (t % score_bins(j)) case (SCORE_FLUX) string = trim(string) // ' flux' + case (SCORE_FLUX_YN) + pn_string = ' flux' + string = trim(string) // pn_string + do n = 1, t % moment_order(j) + pn_string = ' flux-y' // trim(to_str(n)) + string = trim(string) // pn_string + end do + j = j + n - 1 case (SCORE_TOTAL) string = trim(string) // ' total' + case (SCORE_TOTAL_YN) + pn_string = ' total' + string = trim(string) // pn_string + do n = 1, t % moment_order(j) + pn_string = ' total-y' // trim(to_str(n)) + string = trim(string) // pn_string + end do + j = j + n - 1 case (SCORE_SCATTER) string = trim(string) // ' scatter' case (SCORE_NU_SCATTER) string = trim(string) // ' nu-scatter' case (SCORE_SCATTER_N) - pn_string = ' scatter-' // trim(to_str(t % scatt_order(j))) + pn_string = ' scatter-' // trim(to_str(t % moment_order(j))) string = trim(string) // pn_string case (SCORE_SCATTER_PN) pn_string = ' scatter' string = trim(string) // pn_string - do n = 1, t % scatt_order(j) - pn_string = ' scatter-' // trim(to_str(n)) + do n = 1, t % moment_order(j) + pn_string = ' scatter-p' // trim(to_str(n)) + string = trim(string) // pn_string + end do + j = j + n - 1 + case (SCORE_NU_SCATTER_N) + pn_string = ' nu-scatter-' // trim(to_str(t % moment_order(j))) + string = trim(string) // pn_string + case (SCORE_NU_SCATTER_PN) + pn_string = ' nu-scatter' + string = trim(string) // pn_string + do n = 1, t % moment_order(j) + pn_string = ' nu-scatter-p' // trim(to_str(n)) + string = trim(string) // pn_string + end do + j = j + n - 1 + case (SCORE_SCATTER_YN) + pn_string = ' scatter' + string = trim(string) // pn_string + do n = 1, t % moment_order(j) + pn_string = ' scatter-y' // trim(to_str(n)) + string = trim(string) // pn_string + end do + j = j + n - 1 + case (SCORE_NU_SCATTER_YN) + pn_string = ' nu-scatter' + string = trim(string) // pn_string + do n = 1, t % moment_order(j) + pn_string = ' nu-scatter-y' // trim(to_str(n)) string = trim(string) // pn_string end do j = j + n - 1 @@ -1614,13 +1657,14 @@ contains integer :: score_index ! scoring bin index integer :: i_nuclide ! index in nuclides array integer :: i_listing ! index in xs_listings array - integer :: n_order ! loop index for scattering orders + integer :: n_order ! loop index for moment orders + integer :: nm_order ! loop index for Ynm moment orders real(8) :: t_value ! t-values for confidence intervals real(8) :: alpha ! significance level for CI character(MAX_FILE_LEN) :: filename ! name of output file character(15) :: filter_name(N_FILTER_TYPES) ! names of tally filters - character(27) :: score_names(N_SCORE_TYPES) ! names of scoring function - character(27) :: score_name ! names of scoring function + character(36) :: score_names(N_SCORE_TYPES) ! names of scoring function + character(36) :: score_name ! names of scoring function ! to be applied at write-time type(TallyObject), pointer :: t @@ -1642,8 +1686,6 @@ contains score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate" score_names(abs(SCORE_SCATTER)) = "Scattering Rate" score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" - score_names(abs(SCORE_SCATTER_N)) = "" - score_names(abs(SCORE_SCATTER_PN)) = "" score_names(abs(SCORE_TRANSPORT)) = "Transport Rate" score_names(abs(SCORE_N_1N)) = "(n,1n) Rate" score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate" @@ -1651,6 +1693,14 @@ contains score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" score_names(abs(SCORE_KAPPA_FISSION)) = "Kappa-Fission Rate" score_names(abs(SCORE_EVENTS)) = "Events" + score_names(abs(SCORE_FLUX_YN)) = "Flux Moment" + score_names(abs(SCORE_TOTAL_YN)) = "Total Reaction Rate Moment" + score_names(abs(SCORE_SCATTER_N)) = "Scattering Rate Moment" + score_names(abs(SCORE_SCATTER_PN)) = "Scattering Rate Moment" + score_names(abs(SCORE_SCATTER_YN)) = "Scattering Rate Moment" + score_names(abs(SCORE_NU_SCATTER_N)) = "Scattering Prod. Rate Moment" + score_names(abs(SCORE_NU_SCATTER_PN)) = "Scattering Prod. Rate Moment" + score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment" ! Create filename for tally output filename = trim(path_output) // "tallies.out" @@ -1777,34 +1827,42 @@ contains do l = 1, t % n_user_score_bins k = k + 1 score_index = score_index + 1 - if (t % score_bins(k) == SCORE_SCATTER_N) then - if (t % scatt_order(k) == 0) then - score_name = "Scattering Rate" - else - score_name = 'P' // trim(to_str(t % scatt_order(k))) // & - ' Scattering Moment' - end if + select case(t % score_bins(k)) + case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) + score_name = 'P' // trim(to_str(t % moment_order(k))) // " " // & + score_names(abs(t % score_bins(k))) write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) - else if (t % score_bins(k) == SCORE_SCATTER_PN) then - score_name = "Scattering Rate" - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) - do n_order = 1, t % scatt_order(k) + case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) + score_index = score_index - 1 + do n_order = 0, t % moment_order(k) score_index = score_index + 1 - score_name = 'P' // trim(to_str(n_order)) // & - ' Scattering Moment' + score_name = 'P' // trim(to_str(n_order)) // " " //& + score_names(abs(t % score_bins(k))) write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) end do - k = k + n_order - 1 - else + k = k + t % moment_order(k) + case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & + SCORE_TOTAL_YN) + score_index = score_index - 1 + do n_order = 0, t % moment_order(k) + do nm_order = -n_order, n_order + score_index = score_index + 1 + score_name = 'Y' // trim(to_str(n_order)) // ',' // & + trim(to_str(nm_order)) // " " // score_names(abs(t % score_bins(k))) + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) % sum_sq)) + end do + end do + k = k + (t % moment_order(k) + 1)**2 - 1 + case default if (t % score_bins(k) > 0) then score_name = reaction_name(t % score_bins(k)) else @@ -1814,7 +1872,7 @@ contains repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) - end if + end select end do indent = indent - 2 diff --git a/src/particle_header.F90 b/src/particle_header.F90 index e1c696ff16..3e9dca7d56 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -30,7 +30,7 @@ module particle_header ! Pointer to next (more local) set of coordinates type(LocalCoord), pointer :: next => null() end type LocalCoord - + !=============================================================================== ! PARTICLE describes the state of a particle being transported through the ! geometry @@ -53,6 +53,7 @@ module particle_header ! Pre-collision physical data real(8) :: last_xyz(3) ! previous coordinates + real(8) :: last_uvw(3) ! previous direction coordinates real(8) :: last_wgt ! pre-collision particle weight real(8) :: last_E ! pre-collision energy real(8) :: absorb_wgt ! weight absorbed for survival biasing @@ -96,7 +97,7 @@ contains type(LocalCoord), pointer :: coord - if (associated(coord)) then + if (associated(coord)) then ! recursively deallocate lower coordinates if (associated(coord % next)) call deallocate_coord(coord%next) diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index c524170a69..ab7b1a1a2f 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -1,4 +1,4 @@ -module particle_restart +module particle_restart use, intrinsic :: ISO_FORTRAN_ENV @@ -88,6 +88,7 @@ contains ! Set particle last attributes p % last_wgt = p % wgt p % last_xyz = p % coord % xyz + p % last_uvw = p % coord % uvw p % last_E = p % E ! Close hdf5 file diff --git a/src/physics.F90 b/src/physics.F90 index 71a5cc6cc3..47da60f2b9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -36,6 +36,7 @@ contains ! Store pre-collision particle properties p % last_wgt = p % wgt p % last_E = p % E + p % last_uvw = p % coord0 % uvw ! Add to collision counter for particle p % n_collision = p % n_collision + 1 diff --git a/src/source.F90 b/src/source.F90 index ced7fe1dc5..79bc49fd57 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -256,6 +256,7 @@ contains p % coord % xyz = src % xyz p % coord % uvw = src % uvw p % last_xyz = src % xyz + p % last_uvw = src % uvw p % E = src % E p % last_E = src % E diff --git a/src/state_point.F90 b/src/state_point.F90 index 6884620e2f..83d768a1ba 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -58,7 +58,7 @@ contains call write_message(1) if (master) then - ! Create statepoint file + ! Create statepoint file call sp % file_create(filename) ! Write file type @@ -216,12 +216,12 @@ contains group="tallies/tally" // to_str(i), length=t % n_nuclide_bins) deallocate(temp_array) - ! Write number of score bins, score bins, and scatt order + ! Write number of score bins, score bins, and moment order call sp % write_data(t % n_score_bins, "n_score_bins", & group="tallies/tally" // to_str(i)) call sp % write_data(t % score_bins, "score_bins", & group="tallies/tally" // to_str(i), length=t % n_score_bins) - call sp % write_data(t % scatt_order, "scatt_order", & + call sp % write_data(t % moment_order, "moment_order", & group="tallies/tally" // to_str(i), length=t % n_score_bins) ! Write number of user score bins @@ -409,7 +409,7 @@ contains ! Copy global tallies into temporary array for reducing n_bins = 2 * N_GLOBAL_TALLIES global_temp(1,:) = global_tallies(:) % sum - global_temp(2,:) = global_tallies(:) % sum_sq + global_temp(2,:) = global_tallies(:) % sum_sq if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into a @@ -430,7 +430,7 @@ contains tallyresult_temp(:,1) % sum = global_temp(1,:) tallyresult_temp(:,1) % sum_sq = global_temp(2,:) - + ! Write out global tallies sum and sum_sq call sp % write_tally_result(tallyresult_temp, "global_tallies", & n1=N_GLOBAL_TALLIES, n2=1) @@ -484,7 +484,7 @@ contains allocate(tallyresult_temp(m,n)) tallyresult_temp(:,:) % sum = tally_temp(1,:,:) tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:) - + ! Write reduced tally results to file call sp % write_tally_result(t % results, "results", & group="tallies/tally" // to_str(i), n1=m, n2=n) @@ -525,7 +525,7 @@ contains integer :: int_array(3) integer, allocatable :: temp_array(:) logical :: source_present - real(8) :: real_array(3) + real(8) :: real_array(3) type(TallyObject), pointer :: t => null() ! Write message @@ -720,7 +720,7 @@ contains group="tallies/tally" // to_str(i)) call sp % read_data(t % score_bins, "score_bins", & group="tallies/tally" // to_str(i), length=t % n_score_bins) - call sp % read_data(t % scatt_order, "scatt_order", & + call sp % read_data(t % moment_order, "moment_order", & group="tallies/tally" // to_str(i), length=t % n_score_bins) ! Write number of user score bins @@ -741,7 +741,7 @@ contains if (path_source_point == path_state_point .and. .not. source_present) then message = "Source bank must be contained in statepoint restart file" call fatal_error() - end if + end if ! Read tallies to master if (master) then @@ -779,20 +779,20 @@ contains end if end if - ! Read source if in eigenvalue mode + ! Read source if in eigenvalue mode if (run_mode == MODE_EIGENVALUE) then ! Check if source was written out separately if (.not. source_present) then - ! Close statepoint file + ! Close statepoint file call sp % file_close() ! Write message message = "Loading source file " // trim(path_source_point) // "..." call write_message(1) - ! Open source file + ! Open source file call sp % file_open(path_source_point, 'r', serial = .false.) ! Read file type diff --git a/src/tally.F90 b/src/tally.F90 index 5fceeb1f75..0b4c23daed 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -41,7 +41,8 @@ contains integer :: i_tally integer :: j ! loop index for scoring bins integer :: k ! loop index for nuclide bins - integer :: n ! loop index for scattering order + integer :: n ! loop index for legendre order + integer :: num_nm ! Number of N,M orders in harmonic integer :: l ! scoring bin loop index, allowing for changing ! position during the loop integer :: filter_index ! single index for single bin @@ -148,6 +149,42 @@ contains end if score = score / material_xs % total + case (SCORE_FLUX_YN) + ! All events score to a flux bin. We actually use a collision + ! estimator since there is no way to count 'events' exactly for + ! the flux + + score_index = score_index - 1 + + ! get the score + 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 + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP case (SCORE_TOTAL) ! All events will score to the total reaction rate. We can just @@ -162,6 +199,41 @@ contains score = last_wgt end if + case (SCORE_TOTAL_YN) + ! All events will score to the total reaction rate. We can just + ! use the weight of the particle entering the collision as the + ! score + + score_index = score_index - 1 + + ! get the score + 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 + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_SCATTER) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP @@ -172,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 @@ -181,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 @@ -189,33 +261,129 @@ contains ! Find the scattering order for a singly requested moment, and ! store its moment contribution. - if (t % scatt_order(j) == 1) then + if (t % moment_order(j) == 1) then score = last_wgt * mu ! avoid function call overhead else - score = last_wgt * calc_pn(t % scatt_order(j), mu) + score = last_wgt * calc_pn(t % moment_order(j), mu) endif case (SCORE_SCATTER_PN) ! Skip any event where the particle didn't scatter if (p % event /= EVENT_SCATTER) then - j = j + t % scatt_order(j) + j = j + t % moment_order(j) cycle SCORE_LOOP end if score_index = score_index - 1 ! Find the scattering order for a collection of requested moments ! and store the moment contribution of each - do n = 0, t % scatt_order(j) + do n = 0, t % moment_order(j) ! determine scoring bin index 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 +!$omp end critical + end do + j = j + t % moment_order(j) + cycle SCORE_LOOP + + case (SCORE_SCATTER_YN) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + j = j + t % moment_order(j) + cycle SCORE_LOOP + end if + score_index = score_index - 1 + + ! Calculate the number of moments from t % moment_order + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + ! get the score of the scattering moment + score = last_wgt * calc_pn(n, mu) + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + + case (SCORE_NU_SCATTER_N) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) cycle SCORE_LOOP + + ! Find the scattering order for a singly requested moment, and + ! store its moment contribution. + + if (t % moment_order(j) == 1) then + score = wgt * mu ! avoid function call overhead + else + score = wgt * calc_pn(t % moment_order(j), mu) + endif + + case (SCORE_NU_SCATTER_PN) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + j = j + t % moment_order(j) + cycle SCORE_LOOP + end if + score_index = score_index - 1 + ! Find the scattering order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + 1 + ! get the score and tally it + score = wgt * calc_pn(n, mu) !$omp critical t % results(score_index, filter_index) % value = & t % results(score_index, filter_index) % value + score !$omp end critical end do - j = j + t % scatt_order(j) + j = j + t % moment_order(j) + cycle SCORE_LOOP + + case (SCORE_NU_SCATTER_YN) + ! Skip any event where the particle didn't scatter + if (p % event /= EVENT_SCATTER) then + j = j + t % moment_order(j) + cycle SCORE_LOOP + end if + score_index = score_index - 1 + + ! Calculate the number of moments from t % moment_order + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + ! get the score of the scattering moment + score = wgt * calc_pn(n, mu) + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % last_uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 cycle SCORE_LOOP case (SCORE_TRANSPORT) @@ -335,12 +503,12 @@ 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 + ! fission scale by kappa-fission if (micro_xs(p % event_nuclide) % absorption > ZERO) then score = p % absorb_wgt * & @@ -479,6 +647,9 @@ contains integer :: k ! loop index for nuclide bins integer :: l ! loop index for nuclides in material integer :: m ! loop index for reactions + integer :: n ! loop index for legendre order + integer :: num_nm ! Number of N,M orders in harmonic + integer :: q ! loop index for scoring bins integer :: filter_index ! single index for single bin integer :: i_nuclide ! index in nuclides array (from bins) integer :: i_nuc ! index in nuclides array (from material) @@ -565,10 +736,15 @@ contains end if ! Determine score for each bin - SCORE_LOOP: do j = 1, t % n_score_bins + j = 0 + SCORE_LOOP: do q = 1, t % n_user_score_bins + j = j + 1 ! determine what type of score bin score_bin = t % score_bins(j) + ! determine scoring bin index + score_index = (k - 1)*t % n_score_bins + j + if (i_nuclide > 0) then ! ================================================================ ! DETERMINE NUCLIDE CROSS SECTION @@ -578,11 +754,58 @@ contains ! For flux, we need no cross section score = flux + case (SCORE_FLUX_YN) + score_index = score_index - 1 + + ! For flux, we need no cross section + score = flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_TOTAL) ! Total cross section is pre-calculated score = micro_xs(i_nuclide) % total * & atom_density * flux + case (SCORE_TOTAL_YN) + score_index = score_index - 1 + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_SCATTER) ! Scattering cross section is pre-calculated score = (micro_xs(i_nuclide) % total - & @@ -662,10 +885,60 @@ contains ! For flux, we need no cross section score = flux + case (SCORE_FLUX_YN) + score_index = score_index - 1 + + ! For flux, we need no cross section + score = flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_TOTAL) ! Total cross section is pre-calculated score = material_xs % total * flux + case (SCORE_TOTAL_YN) + score_index = score_index - 1 + + ! Total cross section is pre-calculated + score = material_xs % total * flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_SCATTER) ! Scattering cross section is pre-calculated score = (material_xs % total - material_xs % absorption) * flux @@ -704,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) @@ -742,9 +1015,6 @@ contains end select end if - ! Determine scoring bin index - score_index = (k - 1)*t % n_score_bins + j - ! Add score to tally !$omp critical t % results(score_index, filter_index) % value = & @@ -785,6 +1055,9 @@ contains integer :: i ! loop index for nuclides in material integer :: j ! loop index for scoring bin types integer :: m ! loop index for reactions in nuclide + integer :: n ! loop index for legendre order + integer :: num_nm ! Number of N,M orders in harmonic + integer :: q ! loop index for scoring bins integer :: i_nuclide ! index in nuclides array integer :: score_bin ! type of score, e.g. SCORE_FLUX integer :: score_index ! scoring bin index @@ -815,18 +1088,73 @@ contains atom_density = mat % atom_density(i) ! Loop over score types for each bin - SCORE_LOOP: do j = 1, t % n_score_bins + j = 0 + SCORE_LOOP: do q = 1, t % n_user_score_bins + j = j + 1 ! determine what type of score bin score_bin = t % score_bins(j) + ! Determine scoring bin index based on what the index of the nuclide + ! is in the nuclides array + score_index = (i_nuclide - 1)*t % n_score_bins + j + ! Determine macroscopic nuclide cross section select case(score_bin) case (SCORE_FLUX) score = flux + case (SCORE_FLUX_YN) + score_index = score_index - 1 + + ! For flux, we need no cross section + score = flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * atom_density * flux + case (SCORE_TOTAL_YN) + score_index = score_index - 1 + + score = micro_xs(i_nuclide) % total * atom_density * flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & micro_xs(i_nuclide) % absorption) * atom_density * flux @@ -851,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 @@ -886,10 +1214,6 @@ contains end if end select - ! Determine scoring bin index based on what the index of the nuclide - ! is in the nuclides array - score_index = (i_nuclide - 1)*t % n_score_bins + j - ! Add score to tally !$omp critical t % results(score_index, filter_index) % value = & @@ -904,18 +1228,74 @@ contains ! SCORE TOTAL MATERIAL REACTION RATES ! Loop over score types for each bin - MATERIAL_SCORE_LOOP: do j = 1, t % n_score_bins + j = 0 + MATERIAL_SCORE_LOOP: do q = 1, t % n_user_score_bins + j = j + 1 ! determine what type of score bin score_bin = t % score_bins(j) + ! Determine scoring bin index based on what the index of the nuclide + ! is in the nuclides array + score_index = n_nuclides_total*t % n_score_bins + j + ! Determine macroscopic material cross section select case(score_bin) case (SCORE_FLUX) score = flux + case (SCORE_FLUX_YN) + score_index = score_index - 1 + + ! For flux, we need no cross section + score = flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle MATERIAL_SCORE_LOOP + case (SCORE_TOTAL) score = material_xs % total * flux + case (SCORE_TOTAL_YN) + score_index = score_index - 1 + + ! Total cross section is pre-calculated + score = material_xs % total * flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle MATERIAL_SCORE_LOOP + case (SCORE_SCATTER) score = (material_xs % total - material_xs % absorption) * flux @@ -938,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 @@ -986,10 +1366,6 @@ contains end if end select - ! Determine scoring bin index based on what the index of the nuclide - ! is in the nuclides array - score_index = n_nuclides_total*t % n_score_bins + j - ! Add score to tally !$omp critical t % results(score_index, filter_index) % value = & @@ -1016,6 +1392,9 @@ contains integer :: j ! loop index for direction integer :: k ! loop index for mesh cell crossings integer :: b ! loop index for nuclide bins + integer :: n ! loop index for legendre order + integer :: num_nm ! Number of N,M orders in harmonic + integer :: q ! loop index for scoring bins integer :: ijk0(3) ! indices of starting coordinates integer :: ijk1(3) ! indices of ending coordinates integer :: ijk_cross(3) ! indices of mesh cell crossed @@ -1241,18 +1620,75 @@ contains end if ! Determine score for each bin - SCORE_LOOP: do j = 1, t % n_score_bins + j = 0 + SCORE_LOOP: do q = 1, t % n_user_score_bins + j = j + 1 ! determine what type of score bin score_bin = t % score_bins(j) + ! Determine scoring bin index + 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 + + case (SCORE_FLUX_YN) + score_index = score_index - 1 + + score = flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * & atom_density * flux + + case (SCORE_TOTAL_YN) + score_index = score_index - 1 + + ! Total cross section is pre-calculated + score = micro_xs(i_nuclide) % total * & + atom_density * flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & micro_xs(i_nuclide) % absorption) * & @@ -1277,12 +1713,63 @@ contains end select else - ! Determine macroscopic material cross section + ! Determine macroscopic material cross section select case(score_bin) case (SCORE_FLUX) score = flux + + case (SCORE_FLUX_YN) + score_index = score_index - 1 + + score = flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_TOTAL) score = material_xs % total * flux + + case (SCORE_TOTAL_YN) + score_index = score_index - 1 + + ! Total cross section is pre-calculated + score = material_xs % total * flux + + num_nm = 1 + ! Find the order for a collection of requested moments + ! and store the moment contribution of each + do n = 0, t % moment_order(j) + ! determine scoring bin index + score_index = score_index + num_nm + ! Update number of total n,m bins for this n (m = [-n: n]) + num_nm = 2 * n + 1 + + ! multiply score by the angular flux moments and store +!$omp critical + t % results(score_index: score_index + num_nm - 1, filter_index) % value = & + t % results(score_index: score_index + num_nm - 1, filter_index) % value + & + score * calc_rn(n, p % coord0 % uvw) +!$omp end critical + end do + j = j + (t % moment_order(j) + 1)**2 - 1 + cycle SCORE_LOOP + case (SCORE_SCATTER) score = (material_xs % total - material_xs % absorption) * flux case (SCORE_ABSORPTION) @@ -1302,9 +1789,6 @@ contains end select end if - ! Determine scoring bin index - score_index = (b - 1)*t % n_score_bins + j - ! Add score to tally !$omp critical t % results(score_index, filter_index) % value = & diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 17cb1cb0de..e440f472e8 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -47,7 +47,7 @@ module tally_header !=============================================================================== ! TALLYFILTER describes a filter that limits what events score to a tally. For ! example, a cell filter indicates that only particles in a specified cell -! should score to the tally. +! should score to the tally. !=============================================================================== type TallyFilter @@ -55,7 +55,7 @@ module tally_header integer :: n_bins = 0 integer, allocatable :: int_bins(:) real(8), allocatable :: real_bins(:) ! Only used for energy filters - + ! Type-Bound procedures contains procedure :: clear => tallyfilter_clear ! Deallocates TallyFilter @@ -105,7 +105,7 @@ module tally_header ! scattering response. integer :: n_score_bins = 0 integer, allocatable :: score_bins(:) - integer, allocatable :: scatt_order(:) + integer, allocatable :: moment_order(:) integer :: n_user_score_bins = 0 ! Results for each bin -- the first dimension of the array is for scores @@ -122,14 +122,14 @@ module tally_header ! Number of realizations of tally random variables integer :: n_realizations = 0 - + ! Type-Bound procedures contains procedure :: clear => tallyobject_clear ! Deallocates TallyObject end type TallyObject - + contains - + !=============================================================================== ! TALLYFILTER_CLEAR deallocates a TallyFilter element and sets it to its as ! initialized state. @@ -137,16 +137,16 @@ module tally_header subroutine tallyfilter_clear(this) class(TallyFilter), intent(inout) :: this ! The TallyFilter to be cleared - + this % type = NONE this % n_bins = 0 if (allocated(this % int_bins)) & deallocate(this % int_bins) if (allocated(this % real_bins)) & deallocate(this % real_bins) - + end subroutine tallyfilter_clear - + !=============================================================================== ! TALLYOBJECT_CLEAR deallocates a TallyObject element and sets it to its as ! initialized state. @@ -154,44 +154,44 @@ module tally_header subroutine tallyobject_clear(this) class(TallyObject), intent(inout) :: this ! The TallyObject to be cleared - + integer :: i ! Loop Index - + ! This routine will go through each item in TallyObject and set the value ! to its default, as-initialized values, including deallocations. this % label = "" - + if (allocated(this % filters)) then do i = 1, size(this % filters) call this % filters(i) % clear() end do deallocate(this % filters) end if - + if (allocated(this % stride)) & deallocate(this % stride) - + this % find_filter = 0 - + this % n_nuclide_bins = 0 if (allocated(this % nuclide_bins)) & deallocate(this % nuclide_bins) this % all_nuclides = .false. - + this % n_score_bins = 0 if (allocated(this % score_bins)) & deallocate(this % score_bins) - if (allocated(this % scatt_order)) & - deallocate(this % scatt_order) + if (allocated(this % moment_order)) & + deallocate(this % moment_order) this % n_user_score_bins = 0 - + if (allocated(this % results)) & deallocate(this % results) - + this % reset = .false. - + this % n_realizations = 0 - + end subroutine tallyobject_clear end module tally_header diff --git a/src/tracking.F90 b/src/tracking.F90 index c7ee31f0b2..68a0e3d09a 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -66,7 +66,7 @@ contains total_weight = total_weight + p % wgt !$omp end critical - ! Force calculation of cross-sections by setting last energy to zero + ! Force calculation of cross-sections by setting last energy to zero micro_xs % last_E = ZERO ! Prepare to write out particle track. diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 9613421e07..3799aa7253 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -6,81 +6,89 @@ from collections import OrderedDict import numpy as np import scipy.stats +REVISION_STATEPOINT = 12 + filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', 5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'} score_types = {-1: 'flux', - -2: 'total', - -3: 'scatter', - -4: 'nu-scatter', - -5: 'scatter-n', - -6: 'scatter-pn', - -7: 'transport', - -8: 'n1n', - -9: 'absorption', - -10: 'fission', - -11: 'nu-fission', - -12: 'kappa-fission', - -13: 'current', - -14: 'events', - 1: '(n,total)', - 2: '(n,elastic)', - 4: '(n,level)', - 11: '(n,2nd)', - 16: '(n,2n)', - 17: '(n,3n)', - 18: '(n,fission)', - 19: '(n,f)', - 20: '(n,nf)', - 21: '(n,2nf)', - 22: '(n,na)', - 23: '(n,n3a)', - 24: '(n,2na)', - 25: '(n,3na)', - 28: '(n,np)', - 29: '(n,n2a)', - 30: '(n,2n2a)', - 32: '(n,nd)', - 33: '(n,nt)', - 34: '(n,nHe-3)', - 35: '(n,nd2a)', - 36: '(n,nt2a)', - 37: '(n,4n)', - 38: '(n,3nf)', - 41: '(n,2np)', - 42: '(n,3np)', - 44: '(n,n2p)', - 45: '(n,npa)', - 91: '(n,nc)', - 101: '(n,disappear)', - 102: '(n,gamma)', - 103: '(n,p)', - 104: '(n,d)', - 105: '(n,t)', - 106: '(n,3He)', - 107: '(n,a)', - 108: '(n,2a)', - 109: '(n,3a)', - 111: '(n,2p)', - 112: '(n,pa)', - 113: '(n,t2a)', - 114: '(n,d2a)', - 115: '(n,pd)', - 116: '(n,pt)', - 117: '(n,da)', - 201: '(n,Xn)', - 202: '(n,Xgamma)', - 203: '(n,Xp)', - 204: '(n,Xd)', - 205: '(n,Xt)', - 206: '(n,X3He)', - 207: '(n,Xa)', - 444: '(damage)', - 649: '(n,pc)', - 699: '(n,dc)', - 749: '(n,tc)', - 799: '(n,3Hec)', - 849: '(n,tc)'} + -2: 'total', + -3: 'scatter', + -4: 'nu-scatter', + -5: 'scatter-n', + -6: 'scatter-pn', + -7: 'nu-scatter-n', + -8: 'nu-scatter-pn', + -9: 'transport', + -10: 'n1n', + -11: 'absorption', + -12: 'fission', + -13: 'nu-fission', + -14: 'kappa-fission', + -15: 'current', + -16: 'flux-yn', + -17: 'total-yn', + -18: 'scatter-yn', + -19: 'nu-scatter-yn', + -20: 'events', + 1: '(n,total)', + 2: '(n,elastic)', + 4: '(n,level)', + 11: '(n,2nd)', + 16: '(n,2n)', + 17: '(n,3n)', + 18: '(n,fission)', + 19: '(n,f)', + 20: '(n,nf)', + 21: '(n,2nf)', + 22: '(n,na)', + 23: '(n,n3a)', + 24: '(n,2na)', + 25: '(n,3na)', + 28: '(n,np)', + 29: '(n,n2a)', + 30: '(n,2n2a)', + 32: '(n,nd)', + 33: '(n,nt)', + 34: '(n,nHe-3)', + 35: '(n,nd2a)', + 36: '(n,nt2a)', + 37: '(n,4n)', + 38: '(n,3nf)', + 41: '(n,2np)', + 42: '(n,3np)', + 44: '(n,n2p)', + 45: '(n,npa)', + 91: '(n,nc)', + 101: '(n,disappear)', + 102: '(n,gamma)', + 103: '(n,p)', + 104: '(n,d)', + 105: '(n,t)', + 106: '(n,3He)', + 107: '(n,a)', + 108: '(n,2a)', + 109: '(n,3a)', + 111: '(n,2p)', + 112: '(n,pa)', + 113: '(n,t2a)', + 114: '(n,d2a)', + 115: '(n,pd)', + 116: '(n,pt)', + 117: '(n,da)', + 201: '(n,Xn)', + 202: '(n,Xgamma)', + 203: '(n,Xp)', + 204: '(n,Xd)', + 205: '(n,Xt)', + 206: '(n,X3He)', + 207: '(n,Xa)', + 444: '(damage)', + 649: '(n,pc)', + 699: '(n,dc)', + 749: '(n,tc)', + 799: '(n,3Hec)', + 849: '(n,tc)'} score_types.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)}) score_types.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)}) score_types.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)}) @@ -151,7 +159,7 @@ class StatePoint(object): # Read statepoint revision self.revision = self._get_int(path='revision')[0] - if self.revision != 11: + if self.revision != REVISION_STATEPOINT: raise Exception('Statepoint Revision is not consistent.') # Read OpenMC version @@ -287,7 +295,7 @@ class StatePoint(object): t.n_scores = self._get_int(path=base+'n_score_bins')[0] t.scores = [score_types[j] for j in self._get_int( t.n_scores, path=base+'score_bins')] - t.scatt_order = self._get_int(t.n_scores, path=base+'scatt_order') + t.moment_order = self._get_int(t.n_scores, path=base+'moment_order') # Read number of user score bins t.n_user_scores = self._get_int(path=base+'n_user_score_bins')[0] @@ -301,7 +309,7 @@ class StatePoint(object): # Source bank present source_present = self._get_int(path='source_present')[0] if source_present == 1: - self.source_present = True + self.source_present = True else: self.source_present = False diff --git a/tests/test_score_flux_yn/geometry.xml b/tests/test_score_flux_yn/geometry.xml new file mode 100644 index 0000000000..2407bf74b3 --- /dev/null +++ b/tests/test_score_flux_yn/geometry.xml @@ -0,0 +1,185 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + + + + \ No newline at end of file diff --git a/tests/test_score_flux_yn/materials.xml b/tests/test_score_flux_yn/materials.xml new file mode 100644 index 0000000000..9c0b74f3f1 --- /dev/null +++ b/tests/test_score_flux_yn/materials.xml @@ -0,0 +1,272 @@ + + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_score_flux_yn/results.py b/tests/test_score_flux_yn/results.py new file mode 100644 index 0000000000..f4c32a893f --- /dev/null +++ b/tests/test_score_flux_yn/results.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +import sys +import numpy as np + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# extract tally results and convert to vector +results1 = sp.tallies[0].results +shape1 = results1.shape +size1 = (np.product(shape1)) +results1 = np.reshape(results1, size1) +results2 = sp.tallies[1].results +shape2 = results2.shape +size2 = (np.product(shape2)) +results2 = np.reshape(results2, size2) + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write out tally results +outstr += 'tally 1:\n' +for item in results1: + outstr += "{0:12.6E}\n".format(item) +outstr += 'tally 2:\n' +for item in results2: + outstr += "{0:12.6E}\n".format(item) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_score_flux_yn/results_true.dat b/tests/test_score_flux_yn/results_true.dat new file mode 100644 index 0000000000..e586c894ac --- /dev/null +++ b/tests/test_score_flux_yn/results_true.dat @@ -0,0 +1,448 @@ +k-combined: +9.785262E-01 2.895177E-02 +tally 1: +2.435729E+01 +1.203987E+02 +8.194681E+00 +1.380086E+01 +3.940735E+01 +3.173417E+02 +3.867197E+01 +3.002989E+02 +1.337306E+01 +3.597709E+01 +6.603960E+01 +8.773418E+02 +tally 2: +2.435729E+01 +1.203987E+02 +-1.312433E-01 +2.122879E-01 +2.116196E-01 +6.066255E-02 +-5.313614E-01 +2.706195E-01 +-1.075684E-01 +7.055293E-02 +1.243217E-01 +4.159849E-02 +2.468823E-01 +1.768096E-01 +-1.123588E+00 +3.763177E-01 +-1.648770E-02 +1.832803E-01 +2.802677E-01 +5.082361E-02 +-2.061149E-01 +3.165317E-02 +-1.807522E-01 +5.391268E-02 +4.530400E-01 +1.087546E-01 +-6.668795E-01 +1.641587E-01 +-2.043175E-02 +5.328027E-02 +-4.884041E-01 +5.671422E-02 +1.284666E-01 +1.392098E-02 +-1.432327E-01 +4.612313E-02 +2.306910E-01 +5.105478E-02 +1.541349E-01 +3.338018E-02 +1.058470E-01 +5.362792E-02 +-1.047112E-01 +4.224258E-02 +2.943176E-01 +4.868725E-02 +2.587610E-02 +4.315807E-02 +-3.523267E-02 +3.992323E-02 +2.858380E-02 +2.103937E-02 +-8.696462E-02 +2.368855E-02 +-1.389638E-01 +2.956421E-02 +-3.501207E-01 +3.399027E-02 +1.902988E-01 +3.608903E-02 +5.327172E-02 +3.572644E-02 +-2.426249E-01 +1.431615E-02 +-3.050565E-01 +5.036758E-02 +-2.450581E-01 +3.287240E-02 +3.917853E-01 +5.203000E-02 +9.280421E-02 +3.986423E-02 +8.194681E+00 +1.380086E+01 +-1.198609E-01 +2.695675E-02 +1.383659E-01 +2.636170E-02 +-1.749412E-01 +2.760502E-02 +-9.940427E-02 +1.568240E-02 +6.891264E-02 +1.026583E-02 +1.590465E-01 +1.034971E-02 +-3.138227E-01 +2.840189E-02 +-2.863762E-02 +1.621002E-02 +1.352890E-01 +6.206260E-03 +-4.004307E-02 +3.615784E-03 +-7.902125E-02 +1.613512E-02 +1.951944E-01 +2.062619E-02 +-2.509090E-01 +2.990177E-02 +-7.269079E-03 +3.049242E-03 +-5.414351E-02 +3.014519E-03 +6.260478E-02 +1.067701E-03 +-6.869296E-02 +3.560172E-03 +7.374905E-02 +8.314559E-03 +8.279954E-02 +3.348473E-03 +5.763034E-02 +1.854541E-02 +-4.318960E-02 +3.464200E-03 +7.982046E-02 +1.229306E-02 +1.040014E-02 +1.053815E-02 +2.346513E-02 +3.165815E-03 +1.465521E-02 +1.948357E-03 +-1.095875E-02 +2.265977E-03 +-3.900368E-02 +3.842750E-03 +-2.042671E-06 +1.832329E-03 +5.396847E-02 +3.217326E-03 +2.982176E-02 +1.176203E-02 +-2.814986E-02 +1.536677E-03 +-1.081965E-01 +5.407495E-03 +-3.828136E-02 +2.878498E-03 +7.710669E-02 +3.209292E-03 +6.098258E-02 +2.186137E-03 +3.940735E+01 +3.173417E+02 +-4.295696E-01 +2.690844E-01 +5.348104E-01 +2.016075E-01 +-2.609804E-01 +2.138709E-01 +-5.728795E-01 +2.250658E-01 +1.642274E-01 +3.798058E-01 +-4.820927E-02 +1.782804E-01 +-1.217072E+00 +4.264214E-01 +7.722249E-02 +1.569837E-01 +4.762644E-01 +1.486818E-01 +-1.154941E-01 +7.167031E-02 +-5.131035E-01 +3.048317E-01 +4.549123E-01 +1.380401E-01 +-9.594184E-01 +3.120862E-01 +-1.980919E-01 +1.051642E-01 +-1.318185E-01 +4.789331E-02 +1.375739E-01 +3.256717E-02 +-2.735374E-01 +4.667995E-02 +1.737686E-01 +6.821558E-02 +2.284751E-01 +4.756464E-02 +-2.076359E-01 +9.139686E-02 +-1.203908E-02 +1.955920E-02 +5.958901E-01 +1.930652E-01 +4.040923E-01 +1.903999E-01 +-2.109673E-02 +1.771182E-01 +2.678499E-03 +3.839698E-02 +-1.291510E-01 +8.549633E-02 +-5.355609E-02 +8.872016E-02 +-2.460868E-01 +1.952909E-02 +4.739432E-01 +7.754571E-02 +-5.944508E-02 +4.505470E-02 +2.471670E-01 +2.965480E-02 +-4.979193E-01 +1.916606E-01 +-2.393405E-02 +7.535656E-02 +5.654207E-01 +1.170234E-01 +3.792949E-01 +7.850400E-02 +3.867197E+01 +3.002989E+02 +8.242962E-01 +3.328988E-01 +-8.721432E-01 +5.499353E-01 +-2.827445E-02 +3.808892E-01 +3.027429E-01 +3.196213E-01 +-2.668704E-01 +1.775652E-01 +-1.364711E-01 +3.174950E-02 +1.693700E-01 +3.777266E-01 +3.705874E-01 +1.301648E-01 +-6.953870E-01 +2.916220E-01 +-1.232998E-01 +2.594878E-01 +-6.396748E-01 +1.781359E-01 +-4.475990E-01 +8.796908E-02 +-6.138249E-05 +2.637576E-02 +-3.625609E-01 +1.071480E-01 +-3.679703E-01 +9.529230E-02 +1.887420E-01 +9.572426E-02 +4.512924E-01 +8.828619E-02 +-5.301536E-01 +8.836266E-02 +-5.682538E-01 +1.322901E-01 +-5.640185E-01 +2.067701E-01 +3.060929E-04 +1.385243E-01 +-4.297178E-02 +1.481975E-02 +3.664976E-01 +1.014056E-01 +-1.205351E-01 +5.420132E-03 +-2.494761E-01 +2.833111E-02 +-6.076287E-02 +1.659625E-02 +1.750036E-01 +5.761870E-02 +-1.153502E-02 +3.068594E-02 +1.175234E-01 +5.121138E-02 +-1.125845E-01 +2.987398E-02 +1.480797E-01 +1.698824E-02 +-1.492174E-01 +5.382400E-02 +1.748708E-01 +5.434295E-02 +1.527959E-01 +8.471653E-02 +-2.904372E-01 +1.339215E-01 +1.337306E+01 +3.597709E+01 +1.853470E-01 +3.391063E-02 +-3.879575E-01 +4.413806E-02 +1.110756E-01 +4.326819E-02 +1.272136E-02 +3.910615E-02 +-1.708158E-01 +2.766384E-02 +-7.586373E-03 +3.074551E-03 +1.355828E-01 +4.352850E-02 +1.720338E-01 +3.816177E-02 +-1.970199E-01 +3.589295E-02 +-2.050677E-01 +2.927217E-02 +-1.954428E-01 +2.960696E-02 +-9.327484E-02 +1.166779E-02 +-1.277191E-02 +1.154098E-02 +-6.068542E-02 +1.384060E-02 +-1.322775E-01 +1.720248E-02 +5.815605E-02 +1.135370E-02 +5.351119E-02 +1.134169E-02 +-1.303371E-01 +1.231369E-02 +-1.952041E-01 +8.455846E-03 +-1.885871E-01 +1.963780E-02 +4.399038E-02 +2.402598E-02 +2.517287E-02 +5.478896E-04 +1.156188E-01 +1.077148E-02 +-1.591577E-02 +3.461640E-03 +-7.822656E-02 +2.494263E-03 +-3.635209E-02 +1.957925E-03 +1.353515E-01 +7.473575E-03 +-2.381012E-02 +1.967444E-03 +-5.755770E-02 +3.281243E-03 +5.989480E-03 +3.300085E-03 +-6.040364E-02 +4.806497E-03 +2.536904E-02 +4.434017E-03 +3.776001E-02 +1.478745E-02 +3.411739E-02 +7.395006E-03 +-1.019437E-01 +2.119077E-02 +6.603960E+01 +8.773418E+02 +1.757705E+00 +9.979547E-01 +-9.291194E-01 +9.718246E-01 +6.503126E-01 +7.550216E-01 +1.385724E-01 +5.626979E-01 +-6.450707E-01 +5.509520E-01 +1.145305E-01 +2.637490E-01 +6.041981E-01 +1.002567E+00 +7.176775E-01 +8.116244E-01 +-9.631121E-01 +1.095794E+00 +-5.595901E-01 +5.480297E-01 +-5.984362E-01 +5.063932E-01 +-2.603570E-01 +1.151510E-01 +-8.686749E-02 +1.730599E-01 +-2.120689E-01 +3.466831E-01 +-4.031881E-01 +2.600876E-01 +1.789517E-01 +1.259649E-01 +7.595657E-01 +2.791850E-01 +-7.179829E-01 +1.692839E-01 +-6.844193E-01 +1.655593E-01 +-4.480497E-01 +2.887629E-01 +2.385511E-01 +4.148430E-01 +-1.091249E-01 +1.239082E-01 +5.107831E-01 +1.995329E-01 +3.584899E-02 +3.709753E-02 +-5.233802E-01 +9.374099E-02 +-5.070384E-02 +4.140585E-03 +2.924110E-01 +1.951736E-01 +-2.947218E-02 +5.834542E-02 +-2.679404E-03 +1.905263E-02 +1.668483E-01 +1.361507E-01 +-1.784183E-01 +3.440944E-02 +2.258658E-02 +1.802431E-01 +2.878469E-01 +2.782490E-01 +8.842691E-02 +1.041700E-01 +-5.403488E-01 +5.187324E-01 diff --git a/tests/test_score_flux_yn/settings.xml b/tests/test_score_flux_yn/settings.xml new file mode 100644 index 0000000000..517637a59f --- /dev/null +++ b/tests/test_score_flux_yn/settings.xml @@ -0,0 +1,19 @@ + + + + + 10 + 5 + 100 + + + + + + -160 -160 -183 + 160 160 183 + + + + + diff --git a/tests/test_score_flux_yn/tallies.xml b/tests/test_score_flux_yn/tallies.xml new file mode 100644 index 0000000000..71301d0e7f --- /dev/null +++ b/tests/test_score_flux_yn/tallies.xml @@ -0,0 +1,14 @@ + + + + + + flux + + + + + flux-y5 + + + diff --git a/tests/test_score_flux_yn/test_score_flux_yn.py b/tests/test_score_flux_yn/test_score_flux_yn.py new file mode 100644 index 0000000000..6e54389fc1 --- /dev/null +++ b/tests/test_score_flux_yn/test_score_flux_yn.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_output_exists(): + assert os.path.exists(os.path.join(cwd, 'tallies.out')), 'Tally output file does not exist.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'tallies.out')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + test_run() + test_created_statepoint() + test_output_exists() + test_results() + teardown() diff --git a/tests/test_score_nuscatter_n/geometry.xml b/tests/test_score_nuscatter_n/geometry.xml new file mode 100644 index 0000000000..2407bf74b3 --- /dev/null +++ b/tests/test_score_nuscatter_n/geometry.xml @@ -0,0 +1,185 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + + + + \ No newline at end of file diff --git a/tests/test_score_nuscatter_n/materials.xml b/tests/test_score_nuscatter_n/materials.xml new file mode 100644 index 0000000000..9c0b74f3f1 --- /dev/null +++ b/tests/test_score_nuscatter_n/materials.xml @@ -0,0 +1,272 @@ + + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_score_nuscatter_n/results.py b/tests/test_score_nuscatter_n/results.py new file mode 100644 index 0000000000..93371960ae --- /dev/null +++ b/tests/test_score_nuscatter_n/results.py @@ -0,0 +1,37 @@ +#!/usr/bin/env python + +import sys +import numpy as np + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# extract tally results and convert to vector +results = sp.tallies[0].results +shape = results.shape +size = (np.product(shape)) +results = np.reshape(results, size) + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write out tally results +outstr += 'tallies:\n' +for item in results: + outstr += "{0:12.6E}\n".format(item) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_score_nuscatter_n/results_true.dat b/tests/test_score_nuscatter_n/results_true.dat new file mode 100644 index 0000000000..f0bdd471cd --- /dev/null +++ b/tests/test_score_nuscatter_n/results_true.dat @@ -0,0 +1,33 @@ +k-combined: +9.785262E-01 2.895177E-02 +tallies: +9.630000E+00 +1.882770E+01 +9.951694E-01 +2.306142E-01 +4.578746E-01 +5.623761E-02 +4.846244E-01 +6.415248E-02 +1.908107E-01 +1.203246E-02 +1.920000E+00 +7.570000E-01 +2.543372E-01 +1.948219E-02 +1.781533E-01 +7.666342E-03 +2.613666E-02 +1.105500E-03 +-1.067567E-01 +3.627605E-03 +3.114000E+01 +1.988874E+02 +1.607479E+01 +5.265956E+01 +5.942231E+00 +7.358984E+00 +5.014911E-01 +1.018317E-01 +-7.566804E-01 +1.535164E-01 diff --git a/tests/test_score_nuscatter_n/settings.xml b/tests/test_score_nuscatter_n/settings.xml new file mode 100644 index 0000000000..517637a59f --- /dev/null +++ b/tests/test_score_nuscatter_n/settings.xml @@ -0,0 +1,19 @@ + + + + + 10 + 5 + 100 + + + + + + -160 -160 -183 + 160 160 183 + + + + + diff --git a/tests/test_score_nuscatter_n/tallies.xml b/tests/test_score_nuscatter_n/tallies.xml new file mode 100644 index 0000000000..90e28227b6 --- /dev/null +++ b/tests/test_score_nuscatter_n/tallies.xml @@ -0,0 +1,9 @@ + + + + + + nu-scatter nu-scatter-1 nu-scatter-2 nu-scatter-3 nu-scatter-4 + + + diff --git a/tests/test_score_nuscatter_n/test_score_nuscatter_n.py b/tests/test_score_nuscatter_n/test_score_nuscatter_n.py new file mode 100644 index 0000000000..6e54389fc1 --- /dev/null +++ b/tests/test_score_nuscatter_n/test_score_nuscatter_n.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_output_exists(): + assert os.path.exists(os.path.join(cwd, 'tallies.out')), 'Tally output file does not exist.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'tallies.out')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + test_run() + test_created_statepoint() + test_output_exists() + test_results() + teardown() diff --git a/tests/test_score_nuscatter_pn/geometry.xml b/tests/test_score_nuscatter_pn/geometry.xml new file mode 100644 index 0000000000..2407bf74b3 --- /dev/null +++ b/tests/test_score_nuscatter_pn/geometry.xml @@ -0,0 +1,185 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + + + + \ No newline at end of file diff --git a/tests/test_score_nuscatter_pn/materials.xml b/tests/test_score_nuscatter_pn/materials.xml new file mode 100644 index 0000000000..9c0b74f3f1 --- /dev/null +++ b/tests/test_score_nuscatter_pn/materials.xml @@ -0,0 +1,272 @@ + + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_score_nuscatter_pn/results.py b/tests/test_score_nuscatter_pn/results.py new file mode 100644 index 0000000000..f4c32a893f --- /dev/null +++ b/tests/test_score_nuscatter_pn/results.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +import sys +import numpy as np + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# extract tally results and convert to vector +results1 = sp.tallies[0].results +shape1 = results1.shape +size1 = (np.product(shape1)) +results1 = np.reshape(results1, size1) +results2 = sp.tallies[1].results +shape2 = results2.shape +size2 = (np.product(shape2)) +results2 = np.reshape(results2, size2) + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write out tally results +outstr += 'tally 1:\n' +for item in results1: + outstr += "{0:12.6E}\n".format(item) +outstr += 'tally 2:\n' +for item in results2: + outstr += "{0:12.6E}\n".format(item) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_score_nuscatter_pn/results_true.dat b/tests/test_score_nuscatter_pn/results_true.dat new file mode 100644 index 0000000000..61f40f694c --- /dev/null +++ b/tests/test_score_nuscatter_pn/results_true.dat @@ -0,0 +1,24 @@ +k-combined: +9.785262E-01 2.895177E-02 +tally 1: +9.630000E+00 +1.882770E+01 +9.951694E-01 +2.306142E-01 +4.578746E-01 +5.623761E-02 +4.846244E-01 +6.415248E-02 +1.908107E-01 +1.203246E-02 +tally 2: +9.630000E+00 +1.882770E+01 +9.951694E-01 +2.306142E-01 +4.578746E-01 +5.623761E-02 +4.846244E-01 +6.415248E-02 +1.908107E-01 +1.203246E-02 diff --git a/tests/test_score_nuscatter_pn/settings.xml b/tests/test_score_nuscatter_pn/settings.xml new file mode 100644 index 0000000000..517637a59f --- /dev/null +++ b/tests/test_score_nuscatter_pn/settings.xml @@ -0,0 +1,19 @@ + + + + + 10 + 5 + 100 + + + + + + -160 -160 -183 + 160 160 183 + + + + + diff --git a/tests/test_score_nuscatter_pn/tallies.xml b/tests/test_score_nuscatter_pn/tallies.xml new file mode 100644 index 0000000000..dfad48dda2 --- /dev/null +++ b/tests/test_score_nuscatter_pn/tallies.xml @@ -0,0 +1,14 @@ + + + + + + nu-scatter-0 nu-scatter-1 nu-scatter-2 nu-scatter-3 nu-scatter-4 + + + + + nu-scatter-p4 + + + diff --git a/tests/test_score_nuscatter_pn/test_score_nuscatter_pn.py b/tests/test_score_nuscatter_pn/test_score_nuscatter_pn.py new file mode 100644 index 0000000000..6e54389fc1 --- /dev/null +++ b/tests/test_score_nuscatter_pn/test_score_nuscatter_pn.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_output_exists(): + assert os.path.exists(os.path.join(cwd, 'tallies.out')), 'Tally output file does not exist.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'tallies.out')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + test_run() + test_created_statepoint() + test_output_exists() + test_results() + teardown() diff --git a/tests/test_score_nuscatter_yn/geometry.xml b/tests/test_score_nuscatter_yn/geometry.xml new file mode 100644 index 0000000000..2407bf74b3 --- /dev/null +++ b/tests/test_score_nuscatter_yn/geometry.xml @@ -0,0 +1,185 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + + + + \ No newline at end of file diff --git a/tests/test_score_nuscatter_yn/materials.xml b/tests/test_score_nuscatter_yn/materials.xml new file mode 100644 index 0000000000..9c0b74f3f1 --- /dev/null +++ b/tests/test_score_nuscatter_yn/materials.xml @@ -0,0 +1,272 @@ + + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_score_nuscatter_yn/results.py b/tests/test_score_nuscatter_yn/results.py new file mode 100644 index 0000000000..f4c32a893f --- /dev/null +++ b/tests/test_score_nuscatter_yn/results.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +import sys +import numpy as np + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# extract tally results and convert to vector +results1 = sp.tallies[0].results +shape1 = results1.shape +size1 = (np.product(shape1)) +results1 = np.reshape(results1, size1) +results2 = sp.tallies[1].results +shape2 = results2.shape +size2 = (np.product(shape2)) +results2 = np.reshape(results2, size2) + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write out tally results +outstr += 'tally 1:\n' +for item in results1: + outstr += "{0:12.6E}\n".format(item) +outstr += 'tally 2:\n' +for item in results2: + outstr += "{0:12.6E}\n".format(item) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_score_nuscatter_yn/results_true.dat b/tests/test_score_nuscatter_yn/results_true.dat new file mode 100644 index 0000000000..4854c866e8 --- /dev/null +++ b/tests/test_score_nuscatter_yn/results_true.dat @@ -0,0 +1,38 @@ +k-combined: +9.785262E-01 2.895177E-02 +tally 1: +9.630000E+00 +1.882770E+01 +tally 2: +9.630000E+00 +1.882770E+01 +-3.859107E-02 +5.633136E-03 +1.605692E-01 +1.156085E-02 +-3.977190E-02 +6.818184E-03 +3.522369E-02 +3.789031E-03 +-7.407836E-02 +2.699128E-03 +-2.238973E-02 +4.031197E-03 +-8.969310E-02 +2.769455E-03 +4.330738E-02 +1.714798E-03 +-2.173918E-03 +4.710249E-03 +2.499041E-02 +1.104051E-03 +-7.071751E-02 +2.855396E-03 +-5.031274E-02 +1.551567E-03 +-8.277227E-02 +2.137523E-03 +-1.616867E-02 +1.338600E-03 +-1.497058E-02 +2.014934E-03 diff --git a/tests/test_score_nuscatter_yn/settings.xml b/tests/test_score_nuscatter_yn/settings.xml new file mode 100644 index 0000000000..517637a59f --- /dev/null +++ b/tests/test_score_nuscatter_yn/settings.xml @@ -0,0 +1,19 @@ + + + + + 10 + 5 + 100 + + + + + + -160 -160 -183 + 160 160 183 + + + + + diff --git a/tests/test_score_nuscatter_yn/tallies.xml b/tests/test_score_nuscatter_yn/tallies.xml new file mode 100644 index 0000000000..b80b37dc8c --- /dev/null +++ b/tests/test_score_nuscatter_yn/tallies.xml @@ -0,0 +1,14 @@ + + + + + + nu-scatter-0 + + + + + nu-scatter-y3 + + + diff --git a/tests/test_score_nuscatter_yn/test_score_nuscatter_yn.py b/tests/test_score_nuscatter_yn/test_score_nuscatter_yn.py new file mode 100644 index 0000000000..6e54389fc1 --- /dev/null +++ b/tests/test_score_nuscatter_yn/test_score_nuscatter_yn.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_output_exists(): + assert os.path.exists(os.path.join(cwd, 'tallies.out')), 'Tally output file does not exist.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'tallies.out')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + test_run() + test_created_statepoint() + test_output_exists() + test_results() + teardown() diff --git a/tests/test_score_scatter_yn/geometry.xml b/tests/test_score_scatter_yn/geometry.xml new file mode 100644 index 0000000000..2407bf74b3 --- /dev/null +++ b/tests/test_score_scatter_yn/geometry.xml @@ -0,0 +1,185 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + + + + \ No newline at end of file diff --git a/tests/test_score_scatter_yn/materials.xml b/tests/test_score_scatter_yn/materials.xml new file mode 100644 index 0000000000..9c0b74f3f1 --- /dev/null +++ b/tests/test_score_scatter_yn/materials.xml @@ -0,0 +1,272 @@ + + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_score_scatter_yn/results.py b/tests/test_score_scatter_yn/results.py new file mode 100644 index 0000000000..f4c32a893f --- /dev/null +++ b/tests/test_score_scatter_yn/results.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +import sys +import numpy as np + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# extract tally results and convert to vector +results1 = sp.tallies[0].results +shape1 = results1.shape +size1 = (np.product(shape1)) +results1 = np.reshape(results1, size1) +results2 = sp.tallies[1].results +shape2 = results2.shape +size2 = (np.product(shape2)) +results2 = np.reshape(results2, size2) + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write out tally results +outstr += 'tally 1:\n' +for item in results1: + outstr += "{0:12.6E}\n".format(item) +outstr += 'tally 2:\n' +for item in results2: + outstr += "{0:12.6E}\n".format(item) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_score_scatter_yn/results_true.dat b/tests/test_score_scatter_yn/results_true.dat new file mode 100644 index 0000000000..617d117058 --- /dev/null +++ b/tests/test_score_scatter_yn/results_true.dat @@ -0,0 +1,56 @@ +k-combined: +9.785262E-01 2.895177E-02 +tally 1: +9.630000E+00 +1.882770E+01 +tally 2: +9.630000E+00 +1.882770E+01 +-3.859107E-02 +5.633136E-03 +1.605692E-01 +1.156085E-02 +-3.977190E-02 +6.818184E-03 +3.522369E-02 +3.789031E-03 +-7.407836E-02 +2.699128E-03 +-2.238973E-02 +4.031197E-03 +-8.969310E-02 +2.769455E-03 +4.330738E-02 +1.714798E-03 +-2.173918E-03 +4.710249E-03 +2.499041E-02 +1.104051E-03 +-7.071751E-02 +2.855396E-03 +-5.031274E-02 +1.551567E-03 +-8.277227E-02 +2.137523E-03 +-1.616867E-02 +1.338600E-03 +-1.497058E-02 +2.014934E-03 +-2.683538E-02 +1.569513E-03 +-4.613955E-03 +9.640026E-04 +8.947303E-03 +3.299483E-04 +-8.588404E-03 +7.913430E-04 +-1.510489E-02 +2.519720E-03 +3.236471E-03 +1.675379E-03 +-1.237888E-02 +9.699709E-04 +4.744373E-02 +1.066913E-03 +-4.995487E-02 +1.125835E-03 diff --git a/tests/test_score_scatter_yn/settings.xml b/tests/test_score_scatter_yn/settings.xml new file mode 100644 index 0000000000..517637a59f --- /dev/null +++ b/tests/test_score_scatter_yn/settings.xml @@ -0,0 +1,19 @@ + + + + + 10 + 5 + 100 + + + + + + -160 -160 -183 + 160 160 183 + + + + + diff --git a/tests/test_score_scatter_yn/tallies.xml b/tests/test_score_scatter_yn/tallies.xml new file mode 100644 index 0000000000..3c46c13d40 --- /dev/null +++ b/tests/test_score_scatter_yn/tallies.xml @@ -0,0 +1,15 @@ + + + + + + scatter-0 + analog + + + + + scatter-y4 + + + diff --git a/tests/test_score_scatter_yn/test_score_scatter_yn.py b/tests/test_score_scatter_yn/test_score_scatter_yn.py new file mode 100644 index 0000000000..6e54389fc1 --- /dev/null +++ b/tests/test_score_scatter_yn/test_score_scatter_yn.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_output_exists(): + assert os.path.exists(os.path.join(cwd, 'tallies.out')), 'Tally output file does not exist.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'tallies.out')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + test_run() + test_created_statepoint() + test_output_exists() + test_results() + teardown() diff --git a/tests/test_score_total_yn/geometry.xml b/tests/test_score_total_yn/geometry.xml new file mode 100644 index 0000000000..2407bf74b3 --- /dev/null +++ b/tests/test_score_total_yn/geometry.xml @@ -0,0 +1,185 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 2 1 1 2 1 1 2 1 1 2 1 1 2 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 2 1 1 1 1 1 1 1 1 1 2 1 1 1 + 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + rectangular + 17 17 + -10.71 -10.71 + 1.26 1.26 + + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 4 3 3 4 3 3 4 3 3 4 3 3 4 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 4 3 3 3 3 3 3 3 3 3 4 3 3 3 + 3 3 3 3 3 4 3 3 4 3 3 4 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 + 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 + 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 5 5 5 5 5 + 5 5 5 5 5 5 5 6 6 6 6 6 6 6 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 + + + + + + rectangular + 21 21 + -224.91 -224.91 + 21.42 21.42 + + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 + 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 + 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 7 7 7 7 7 + 7 7 7 7 7 7 7 8 8 8 8 8 8 8 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 + + + + \ No newline at end of file diff --git a/tests/test_score_total_yn/materials.xml b/tests/test_score_total_yn/materials.xml new file mode 100644 index 0000000000..9c0b74f3f1 --- /dev/null +++ b/tests/test_score_total_yn/materials.xml @@ -0,0 +1,272 @@ + + + + 71c + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_score_total_yn/results.py b/tests/test_score_total_yn/results.py new file mode 100644 index 0000000000..f4c32a893f --- /dev/null +++ b/tests/test_score_total_yn/results.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +import sys +import numpy as np + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# extract tally results and convert to vector +results1 = sp.tallies[0].results +shape1 = results1.shape +size1 = (np.product(shape1)) +results1 = np.reshape(results1, size1) +results2 = sp.tallies[1].results +shape2 = results2.shape +size2 = (np.product(shape2)) +results2 = np.reshape(results2, size2) + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write out tally results +outstr += 'tally 1:\n' +for item in results1: + outstr += "{0:12.6E}\n".format(item) +outstr += 'tally 2:\n' +for item in results2: + outstr += "{0:12.6E}\n".format(item) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_score_total_yn/results_true.dat b/tests/test_score_total_yn/results_true.dat new file mode 100644 index 0000000000..1722472dfd --- /dev/null +++ b/tests/test_score_total_yn/results_true.dat @@ -0,0 +1,212 @@ +k-combined: +9.785262E-01 2.895177E-02 +tally 1: +0.000000E+00 +0.000000E+00 +1.094646E+01 +2.428429E+01 +2.225169E+00 +1.016830E+00 +3.150587E+01 +2.052037E+02 +tally 2: +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.094646E+01 +2.428429E+01 +5.846652E-02 +4.670734E-02 +1.274633E-01 +5.689843E-03 +-2.769645E-01 +3.469519E-02 +-9.770271E-02 +2.212600E-02 +-1.718138E-03 +6.195671E-03 +5.293044E-02 +1.357406E-02 +-4.557967E-01 +5.262731E-02 +-5.968883E-02 +3.239875E-02 +5.152732E-02 +5.322314E-03 +-4.737115E-02 +2.543009E-03 +-4.137677E-02 +7.737200E-03 +2.349356E-01 +3.675562E-02 +-2.118782E-01 +2.259914E-02 +9.771726E-02 +7.657338E-03 +-1.425459E-01 +6.661165E-03 +5.503163E-02 +3.547136E-03 +-2.043861E-02 +6.938007E-03 +9.716884E-02 +7.890404E-03 +4.580386E-02 +1.770664E-03 +9.386179E-02 +1.633861E-02 +-1.267384E-01 +1.875299E-02 +1.489661E-01 +1.033285E-02 +-1.653923E-02 +4.659923E-03 +-7.241633E-02 +6.087112E-03 +2.225169E+00 +1.016830E+00 +-3.287201E-02 +3.113510E-03 +5.276917E-02 +2.373544E-03 +-3.047290E-02 +2.066181E-03 +-2.215745E-02 +1.448141E-03 +-3.197517E-03 +1.029606E-03 +3.785832E-02 +5.476722E-04 +-7.694649E-02 +2.322476E-03 +-1.318499E-02 +1.294712E-03 +3.252651E-02 +7.647996E-04 +-8.446605E-03 +2.396205E-04 +-2.387179E-02 +5.981366E-04 +4.490141E-02 +9.076461E-04 +-5.073474E-02 +1.479816E-03 +2.337733E-03 +3.629619E-04 +-1.592142E-02 +1.932660E-04 +2.817351E-02 +2.756714E-04 +-7.179521E-03 +2.705344E-04 +5.563912E-03 +5.688500E-04 +2.382506E-02 +2.909462E-04 +1.753681E-02 +6.909020E-04 +-2.592268E-02 +3.202222E-04 +2.160865E-02 +9.134735E-04 +-4.319667E-03 +7.902430E-04 +1.070167E-02 +2.288415E-04 +3.150587E+01 +2.052037E+02 +-1.508479E-01 +3.017029E-01 +1.402160E-01 +7.296606E-02 +-6.067590E-01 +1.366120E-01 +-4.674088E-01 +1.052625E-01 +4.292988E-02 +1.212020E-01 +-8.705007E-02 +7.120546E-02 +-4.832242E-01 +9.504346E-02 +-1.625084E-01 +1.004166E-01 +7.545226E-02 +5.868337E-02 +-6.752172E-03 +3.870347E-02 +-3.514281E-01 +7.616059E-02 +5.581075E-01 +1.419356E-01 +-4.187866E-01 +5.601005E-02 +2.708391E-01 +6.304397E-02 +4.343467E-02 +7.332731E-02 +-3.196333E-02 +1.037211E-02 +-1.138408E-03 +2.321733E-02 +6.915075E-02 +2.899777E-02 +4.666937E-02 +2.139001E-02 +2.022369E-01 +7.382740E-02 +-1.274694E-01 +1.226053E-02 +3.724451E-01 +1.421241E-01 +1.587808E-01 +7.967636E-02 +-8.004346E-03 +4.594356E-02 diff --git a/tests/test_score_total_yn/settings.xml b/tests/test_score_total_yn/settings.xml new file mode 100644 index 0000000000..517637a59f --- /dev/null +++ b/tests/test_score_total_yn/settings.xml @@ -0,0 +1,19 @@ + + + + + 10 + 5 + 100 + + + + + + -160 -160 -183 + 160 160 183 + + + + + diff --git a/tests/test_score_total_yn/tallies.xml b/tests/test_score_total_yn/tallies.xml new file mode 100644 index 0000000000..df6de1c134 --- /dev/null +++ b/tests/test_score_total_yn/tallies.xml @@ -0,0 +1,14 @@ + + + + + + total + + + + + total-y4 + + + diff --git a/tests/test_score_total_yn/test_score_total_yn.py b/tests/test_score_total_yn/test_score_total_yn.py new file mode 100644 index 0000000000..6e54389fc1 --- /dev/null +++ b/tests/test_score_total_yn/test_score_total_yn.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +import glob +from optparse import OptionParser + +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() + +def test_run(): + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0, 'OpenMC did not exit successfully.' + +def test_created_statepoint(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' + +def test_output_exists(): + assert os.path.exists(os.path.join(cwd, 'tallies.out')), 'Tally output file does not exist.' + +def test_results(): + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare, 'Results do not agree.' + +def teardown(): + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'tallies.out')) + output.append(os.path.join(cwd, 'results_test.dat')) + for f in output: + if os.path.exists(f): + os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + test_run() + test_created_statepoint() + test_output_exists() + test_results() + teardown()