diff --git a/docs/source/releasenotes/notes_0.5.0.rst b/docs/source/releasenotes/notes_0.5.0.rst index ddff604ab8..edf4b4b6ce 100644 --- a/docs/source/releasenotes/notes_0.5.0.rst +++ b/docs/source/releasenotes/notes_0.5.0.rst @@ -32,6 +32,9 @@ New Features Bug Fixes --------- +- 737b90_: Coincident surfaces from separate universes / particle traveling + tangent to a surface. +- a819b4_: Output of surface neighbors in summary.out file. - b11696_: Reading long attribute lists in XML input. - 2bd46a_: Search for tallying nuclides when no default_xs specified. - 7a1f08_: Fix word wrapping when writing messages. @@ -39,6 +42,8 @@ Bug Fixes - 6f8d9d_: Set default tally labels. - 6a3a5e_: Fix problem with corner-crossing in lattices. +.. _737b90: https://github.com/mit-crpg/openmc/commit/737b90 +.. _a819b4: https://github.com/mit-crpg/openmc/commit/a819b4 .. _b11696: https://github.com/mit-crpg/openmc/commit/b11696 .. _2bd46a: https://github.com/mit-crpg/openmc/commit/2bd46a .. _7a1f08: https://github.com/mit-crpg/openmc/commit/7a1f08 diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index b568aaa02d..9f89cb5042 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -69,37 +69,6 @@ deviation. *Default*: off -```` Element -------------------------- - -The ```` element indicates that a criticality calculation should be -performed. It has the following attributes/sub-elements: - - :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. - - *Default*: None - - :generations_per_batch: - The number of total fission source iterations per batch. - - *Default*: 1 - - :inactive: - The number of inactive batches. In general, the starting cycles in a - criticality calculation can not be used to contribute to tallies since the - fission source distribution and eigenvalue are generally not converged - immediately. - - *Default*: None - - :particles: - The number of neutrons to simulate per fission source iteration. - - *Default*: None - .. _cross_sections: ```` Element @@ -130,6 +99,37 @@ default. This element has the following attributes/sub-elements: *Default*: 1.0 +```` Element +------------------------ + +The ```` element indicates that a :math:`k`-eigenvalue calculation +should be performed. It has the following attributes/sub-elements: + + :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. + + *Default*: None + + :generations_per_batch: + The number of total fission source iterations per batch. + + *Default*: 1 + + :inactive: + The number of inactive batches. In general, the starting cycles in a + criticality calculation can not be used to contribute to tallies since the + fission source distribution and eigenvalue are generally not converged + immediately. + + *Default*: None + + :particles: + The number of neutrons to simulate per fission source iteration. + + *Default*: None + ```` Element ------------------------- diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 9b7f715b53..0d2e4e2def 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -90,23 +90,23 @@ cmfd_snes_solver.o: constants.o cmfd_snes_solver.o: global.o cmfd_snes_solver.o: string.o -criticality.o: cmfd_execute.o -criticality.o: constants.o -criticality.o: error.o -criticality.o: global.o -criticality.o: hdf5_interface.o -criticality.o: math.o -criticality.o: mesh.o -criticality.o: mesh_header.o -criticality.o: output.o -criticality.o: physics.o -criticality.o: random_lcg.o -criticality.o: search.o -criticality.o: source.o -criticality.o: state_point.o -criticality.o: string.o -criticality.o: tally.o -criticality.o: timing.o +eigenvalue.o: cmfd_execute.o +eigenvalue.o: constants.o +eigenvalue.o: error.o +eigenvalue.o: global.o +eigenvalue.o: hdf5_interface.o +eigenvalue.o: math.o +eigenvalue.o: mesh.o +eigenvalue.o: mesh_header.o +eigenvalue.o: output.o +eigenvalue.o: physics.o +eigenvalue.o: random_lcg.o +eigenvalue.o: search.o +eigenvalue.o: source.o +eigenvalue.o: state_point.o +eigenvalue.o: string.o +eigenvalue.o: tally.o +eigenvalue.o: timing.o cross_section.o: ace_header.o cross_section.o: constants.o @@ -243,7 +243,7 @@ interpolation.o: search.o interpolation.o: string.o main.o: constants.o -main.o: criticality.o +main.o: eigenvalue.o main.o: finalize.o main.o: fixed_source.o main.o: global.o diff --git a/src/OBJECTS b/src/OBJECTS index 324f16d761..1a0aac3828 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -14,10 +14,10 @@ cmfd_power_solver.o \ cmfd_prod_operator.o \ cmfd_snes_solver.o \ cross_section.o \ -criticality.o \ datatypes.o \ datatypes_header.o \ doppler.o \ +eigenvalue.o \ endf.o \ endf_header.o \ energy_grid.o \ diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 2d49dead7e..1a91f86f7a 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -3,7 +3,7 @@ module bank_header implicit none !=============================================================================== -! BANK is used for storing fission sites in criticality calculations. Since all +! BANK is used for storing fission sites in eigenvalue calculations. Since all ! the state information of a neutron is not needed, this type allows sites to be ! stored with less memory !=============================================================================== diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index e6e3665feb..211d0e1489 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -61,7 +61,7 @@ contains meshes use mesh, only: mesh_indices_to_bin use mesh_header, only: StructuredMesh - use tally_header, only: TallyObject, TallyScore + use tally_header, only: TallyObject integer :: nx ! number of mesh cells in x direction integer :: ny ! number of mesh cells in y direction @@ -155,7 +155,7 @@ contains score_index = sum((t % matching_bins - 1) * t%stride) + 1 ! get flux - flux = t % scores(1,score_index) % sum + flux = t % results(1,score_index) % sum cmfd % flux(h,i,j,k) = flux ! detect zero flux @@ -166,13 +166,13 @@ contains end if ! get total rr and convert to total xs - cmfd % totalxs(h,i,j,k) = t % scores(2,score_index) % sum / flux + cmfd % totalxs(h,i,j,k) = t % results(2,score_index) % sum / flux ! get p1 scatter rr and convert to p1 scatter xs - cmfd % p1scattxs(h,i,j,k) = t % scores(3,score_index) % sum / flux + cmfd % p1scattxs(h,i,j,k) = t % results(3,score_index) % sum / flux ! extract diffusion coefficient tally - cmfd % diffusion(h,i,j,k) = t % scores(4,score_index) % sum / flux + cmfd % diffusion(h,i,j,k) = t % results(4,score_index) % sum / flux ! calculate diffusion coefficient ! cmfd % diffcof(h,i,j,k) = ONE/(3.0_8*cmfd%totalxs(h,i,j,k)) @@ -206,16 +206,16 @@ contains score_index = sum((t % matching_bins - 1) * t%stride) + 1 ! get scattering - cmfd % scattxs(h,g,i,j,k) = t % scores(1,score_index) % sum /& + cmfd % scattxs(h,g,i,j,k) = t % results(1,score_index) % sum /& cmfd % flux(h,i,j,k) ! get nu-fission - cmfd % nfissxs(h,g,i,j,k) = t % scores(2,score_index) % sum /& + cmfd % nfissxs(h,g,i,j,k) = t % results(2,score_index) % sum /& cmfd % flux(h,i,j,k) ! bank source cmfd % openmc_src(g,i,j,k) = cmfd % openmc_src(g,i,j,k) + & - t % scores(2,score_index) % sum + t % results(2,score_index) % sum end do INGROUP @@ -232,60 +232,60 @@ contains (/ i-1, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(1,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(1,h,i,j,k) = t % results(1,score_index) % sum t % matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming - cmfd % current(2,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(2,h,i,j,k) = t % results(1,score_index) % sum ! right surface t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming - cmfd % current(3,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(3,h,i,j,k) = t % results(1,score_index) % sum t % matching_bins(i_filter_surf) = OUT_RIGHT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(4,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(4,h,i,j,k) = t % results(1,score_index) % sum ! back surface t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j-1, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(5,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(5,h,i,j,k) = t % results(1,score_index) % sum t % matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming - cmfd % current(6,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(6,h,i,j,k) = t % results(1,score_index) % sum ! front surface t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming - cmfd % current(7,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(7,h,i,j,k) = t % results(1,score_index) % sum t % matching_bins(i_filter_surf) = OUT_FRONT score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(8,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(8,h,i,j,k) = t % results(1,score_index) % sum ! bottom surface t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k-1 /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(9,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(9,h,i,j,k) = t % results(1,score_index) % sum t % matching_bins(i_filter_surf) = OUT_TOP score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming - cmfd % current(10,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(10,h,i,j,k) = t % results(1,score_index) % sum ! top surface t % matching_bins(i_filter_mesh) = mesh_indices_to_bin(m, & (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! incoming - cmfd % current(11,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(11,h,i,j,k) = t % results(1,score_index) % sum t % matching_bins(i_filter_surf) = OUT_TOP score_index = sum((t % matching_bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(12,h,i,j,k) = t % scores(1,score_index) % sum + cmfd % current(12,h,i,j,k) = t % results(1,score_index) % sum end if TALLY diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index fed2425822..c962fd997d 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -478,7 +478,7 @@ contains use global, only: n_user_tallies, n_tallies, tallies, message use output, only: write_message - use tally, only: reset_score + use tally, only: reset_result integer :: i ! loop counter @@ -491,7 +491,7 @@ contains ! reset that tally tallies(i) % n_realizations = 0 - call reset_score(tallies(i) % scores) + call reset_result(tallies(i) % results) end do diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index cd1d95fd17..f96dda3868 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -173,7 +173,7 @@ contains use mesh_header, only: StructuredMesh use string use tally, only: setup_active_cmfdtallies - use tally_header, only: TallyObject, TallyScore, TallyFilter + use tally_header, only: TallyObject, TallyFilter use xml_data_cmfd_t integer :: i ! loop counter diff --git a/src/constants.F90 b/src/constants.F90 index e0988df2db..fa73e86250 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -32,6 +32,7 @@ module constants ! User for precision in geometry real(8), parameter :: FP_PRECISION = 1e-14_8 real(8), parameter :: FP_REL_PRECISION = 1e-5_8 + real(8), parameter :: FP_COINCIDENT = 1e-12_8 ! Maximum number of collisions/crossings integer, parameter :: MAX_EVENTS = 10000 @@ -348,7 +349,7 @@ module constants ! Running modes integer, parameter :: & MODE_FIXEDSOURCE = 1, & ! Fixed source mode - MODE_CRITICALITY = 2, & ! Criticality mode + MODE_EIGENVALUE = 2, & ! K eigenvalue mode MODE_PLOTTING = 3, & ! Plotting mode MODE_TALLIES = 4 ! Tally results mode diff --git a/src/criticality.F90 b/src/eigenvalue.F90 similarity index 96% rename from src/criticality.F90 rename to src/eigenvalue.F90 index a19acdc6df..6647603c46 100644 --- a/src/criticality.F90 +++ b/src/eigenvalue.F90 @@ -1,4 +1,4 @@ -module criticality +module eigenvalue #ifdef MPI use mpi @@ -26,18 +26,21 @@ module criticality use hdf5_interface, only: hdf5_write_state_point #endif + private + public :: run_eigenvalue + contains !=============================================================================== -! RUN_CRITICALITY encompasses all the main logic where iterations are performed -! over the batches, generations, and histories. +! RUN_EIGENVALUE encompasses all the main logic where iterations are performed +! over the batches, generations, and histories in a k-eigenvalue calculation. !=============================================================================== - subroutine run_criticality() + subroutine run_eigenvalue() integer(8) :: i ! index over individual particles - if (master) call header("CRITICALITY TRANSPORT SIMULATION", level=1) + if (master) call header("K EIGENVALUE SIMULATION", level=1) ! Allocate particle allocate(p) @@ -45,6 +48,9 @@ contains ! Display column titles call print_columns() + ! Turn on inactive timer + call timer_start(time_inactive) + ! ========================================================================== ! LOOP OVER BATCHES BATCH_LOOP: do current_batch = 1, n_batches @@ -99,7 +105,7 @@ contains if (master) call header("SIMULATION FINISHED", level=1) - end subroutine run_criticality + end subroutine run_eigenvalue !=============================================================================== ! INITIALIZE_BATCH @@ -113,19 +119,22 @@ contains ! Reset total starting particle weight used for normalizing tallies total_weight = ZERO + if (current_batch == n_inactive + 1) then + ! Switch from inactive batch timer to active batch timer + call timer_stop(time_inactive) + call timer_start(time_active) + + ! Enable active batches (and tallies_on if it hasn't been enabled) + active_batches = .true. + tallies_on = .true. + + ! Add user tallies to active tallies list + call setup_active_usertallies() + end if + ! check CMFD initialize batch if (cmfd_run) call cmfd_init_batch() - if (current_batch == n_inactive + 1) then - ! This will start the active timer at the first non-inactive batch - ! (including batch 1 if there are no inactive batches). - call timer_start(time_active) - elseif (current_batch == 1) then - ! If there are inactive batches, start the inactive timer on the first - ! batch. - call timer_start(time_inactive) - end if - end subroutine initialize_batch !=============================================================================== @@ -184,15 +193,6 @@ contains end if end do - ! Turn tallies on once inactive batches are complete - if (current_batch == n_inactive) then - tallies_on = .true. - active_batches = .true. - call timer_stop(time_inactive) - call timer_start(time_active) - call setup_active_usertallies() - end if - end subroutine finalize_batch !=============================================================================== @@ -692,4 +692,4 @@ contains end subroutine count_source_for_ufs -end module criticality +end module eigenvalue diff --git a/src/geometry.F90 b/src/geometry.F90 index 00af1c7fbd..fb15f6f3ce 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -8,7 +8,7 @@ module geometry use output, only: write_message use particle_header, only: LocalCoord, deallocate_coord use string, only: to_str - use tally, only: score_surface_current, add_to_score + use tally, only: score_surface_current implicit none @@ -47,7 +47,7 @@ contains ! Determine the specified sense of the surface in the cell and the actual ! sense of the particle with respect to the surface s => surfaces(abs(i_surface)) - actual_sense = sense(s, p % coord % xyz) + actual_sense = sense(s) specified_sense = (c % surfaces(i) > 0) ! Compare sense of point to specified sense @@ -272,10 +272,9 @@ contains call score_surface_current() end if - if (tallies_on) then - ! Score to global leakage tally - call add_to_score(global_tallies(LEAKAGE), p % wgt) - end if + ! Score to global leakage tally + if (tallies_on) global_tallies(LEAKAGE) % value = & + global_tallies(LEAKAGE) % value + p % wgt ! Display message if (verbosity >= 10 .or. trace) then @@ -481,14 +480,32 @@ contains ! ========================================================================== ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS + ! Remove lower coordinate levels and assignment of surface + p % surface = NONE + p % coord => p % coord0 + call deallocate_coord(p % coord % next) call find_cell(found) - ! Couldn't find next cell anywhere! - if ((.not. found) .and. (run_mode /= MODE_PLOTTING)) then - message = "After particle " // trim(to_str(p % id)) // " crossed surface " & - // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & - &located in any cell and it did not leak." - call fatal_error() + if (run_mode /= MODE_PLOTTING .and. (.not. found)) then + ! If a cell is still not found, there are two possible causes: 1) there is + ! a void in the model, and 2) the particle hit a surface at a tangent. If + ! the particle is really traveling tangent to a surface, if we move it + ! forward a tiny bit it should fix the problem. + + p % coord => p % coord0 + call deallocate_coord(p % coord % next) + p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw + call find_cell(found) + + ! Couldn't find next cell anywhere! This probably means there is an actual + ! undefined region in the geometry. + + if (.not. found) then + message = "After particle " // trim(to_str(p % id)) // " crossed surface " & + // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & + &located in any cell and it did not leak." + call fatal_error() + end if end if end subroutine cross_surface @@ -1142,10 +1159,9 @@ contains ! is in. !=============================================================================== - function sense(surf, xyz) result(s) + recursive function sense(surf) result(s) type(Surface), pointer :: surf ! surface - real(8), intent(in) :: xyz(3) ! coordinates of particle logical :: s ! sense of particle real(8) :: x,y,z ! coordinates of particle @@ -1164,9 +1180,9 @@ contains real(8) :: r ! radius for quadratic surfaces real(8) :: x1,y1,z1 ! upper-right corner of box - x = xyz(1) - y = xyz(2) - z = xyz(3) + x = p % coord % xyz(1) + y = p % coord % xyz(2) + z = p % coord % xyz(3) select case (surf % type) case (SURF_PX) @@ -1320,7 +1336,12 @@ contains end select ! Check which side of surface the point is on - if (func > 0) then + if (abs(func) < FP_COINCIDENT) then + ! Particle may be coincident with this surface. Artifically move the + ! particle forward a tiny bit. + p % coord % xyz = p % coord % xyz + TINY_BIT * p % coord % uvw + s = sense(surf) + elseif (func > 0) then s = .true. else s = .false. diff --git a/src/global.F90 b/src/global.F90 index 47d778bac2..813de56730 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -13,7 +13,7 @@ module global use particle_header, only: Particle use plot_header, only: Plot use source_header, only: ExtSource - use tally_header, only: TallyObject, TallyMap, TallyScore + use tally_header, only: TallyObject, TallyMap, TallyResult use timing, only: Timer #ifdef MPI @@ -111,7 +111,7 @@ module global ! 3) track-length estimate of k-eff ! 4) leakage fraction - type(TallyScore), target :: global_tallies(N_GLOBAL_TALLIES) + type(TallyResult), target :: global_tallies(N_GLOBAL_TALLIES) ! Tally map structure type(TallyMap), allocatable :: tally_maps(:) @@ -154,7 +154,7 @@ module global type(ListInt), pointer :: active_tallies => null() ! ============================================================================ - ! CRITICALITY SIMULATION VARIABLES + ! EIGENVALUE SIMULATION VARIABLES integer(8) :: n_particles = 0 ! # of particles per generation integer :: n_batches ! # of batches @@ -208,7 +208,7 @@ module global logical :: mpi_enabled = .false. ! is MPI in use and initialized? integer :: mpi_err ! MPI error code integer :: MPI_BANK ! MPI datatype for fission bank - integer :: MPI_TALLYSCORE ! MPI datatype for TallyScore + integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult ! No reduction at end of batch logical :: reduce_tallies = .true. @@ -240,17 +240,17 @@ module global ! HDF5 VARIABLES #ifdef HDF5 - integer(HID_T) :: hdf5_output_file ! identifier for output file - integer(HID_T) :: hdf5_tallyscore_t ! Compound type for TallyScore - integer(HID_T) :: hdf5_bank_t ! Compound type for Bank - integer(HID_T) :: hdf5_integer8_t ! type for integer(8) - integer :: hdf5_err ! error flag + integer(HID_T) :: hdf5_output_file ! identifier for output file + integer(HID_T) :: hdf5_tallyresult_t ! Compound type for TallyResult + integer(HID_T) :: hdf5_bank_t ! Compound type for Bank + integer(HID_T) :: hdf5_integer8_t ! type for integer(8) + integer :: hdf5_err ! error flag #endif ! ============================================================================ ! MISCELLANEOUS VARIABLES - ! Mode to run in (fixed source, criticality, plotting, etc) + ! Mode to run in (fixed source, eigenvalue, plotting, etc) integer :: run_mode = NONE ! Restart run diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index c6ec49c4d5..e13be7ca96 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -30,20 +30,20 @@ contains subroutine hdf5_initialize() - type(TallyScore), target :: tmp(2) ! temporary TallyScore - type(Bank), target :: tmpb(2) ! temporary Bank - integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals - integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates + type(TallyResult), target :: tmp(2) ! temporary TallyResult + type(Bank), target :: tmpb(2) ! temporary Bank + integer(HID_T) :: coordinates_t ! HDF5 type for 3 reals + integer(HSIZE_T) :: dims(1) = (/3/) ! size of coordinates ! Initialize FORTRAN interface. call h5open_f(hdf5_err) - ! Create the compound datatype for TallyScore + ! Create the compound datatype for TallyResult call h5tcreate_f(H5T_COMPOUND_F, h5offsetof(c_loc(tmp(1)), & - c_loc(tmp(2))), hdf5_tallyscore_t, hdf5_err) - call h5tinsert_f(hdf5_tallyscore_t, "sum", h5offsetof(c_loc(tmp(1)), & + c_loc(tmp(2))), hdf5_tallyresult_t, hdf5_err) + call h5tinsert_f(hdf5_tallyresult_t, "sum", h5offsetof(c_loc(tmp(1)), & c_loc(tmp(1)%sum)), H5T_NATIVE_DOUBLE, hdf5_err) - call h5tinsert_f(hdf5_tallyscore_t, "sum_sq", h5offsetof(c_loc(tmp(1)), & + call h5tinsert_f(hdf5_tallyresult_t, "sum_sq", h5offsetof(c_loc(tmp(1)), & c_loc(tmp(1)%sum_sq)), H5T_NATIVE_DOUBLE, hdf5_err) ! Create compound type for xyz and uvw @@ -73,7 +73,7 @@ contains subroutine hdf5_finalize() ! Release compound datatypes - call h5tclose_f(hdf5_tallyscore_t, hdf5_err) + call h5tclose_f(hdf5_tallyresult_t, hdf5_err) ! Close FORTRAN interface. call h5close_f(hdf5_err) @@ -94,8 +94,8 @@ contains ! Write header information call hdf5_write_header() - ! Write criticality information - if (run_mode == MODE_CRITICALITY) then + ! Write eigenvalue information + if (run_mode == MODE_EIGENVALUE) then ! Need to write integer(8)'s using double instead since there is no H5LT ! call for making a dataset of type long call hdf5_write_double(hdf5_output_file, "n_particles", real(n_particles,8)) @@ -865,8 +865,8 @@ contains ! Write out current batch number call hdf5_write_integer(hdf5_state_point, "current_batch", current_batch) - ! Write out information for criticality run - if (run_mode == MODE_CRITICALITY) then + ! Write out information for eigenvalue run + if (run_mode == MODE_EIGENVALUE) then call hdf5_write_integer(hdf5_state_point, "n_inactive", n_inactive) call hdf5_write_integer(hdf5_state_point, "gen_per_batch", gen_per_batch) @@ -1031,10 +1031,10 @@ contains ! Write global tallies dims(1) = N_GLOBAL_TALLIES call h5screate_simple_f(1, dims, dspace, hdf5_err) - call h5dcreate_f(hdf5_state_point, "global_tallies", hdf5_tallyscore_t, & + call h5dcreate_f(hdf5_state_point, "global_tallies", hdf5_tallyresult_t, & dspace, dset, hdf5_err) f_ptr = c_loc(global_tallies(1)) - CALL h5dwrite_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err) + CALL h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) @@ -1043,7 +1043,7 @@ contains call hdf5_write_integer(tallies_group, "tallies_present", 1) ! Write tally sum and sum_sq - TALLY_SCORES: do i = 1, n_tallies + TALLY_RESULTS: do i = 1, n_tallies ! Get pointer to tally t => tallies(i) @@ -1052,18 +1052,18 @@ contains temp_group, hdf5_err) ! Write sum and sum_sq for each bin - dims2 = shape(t % scores) + dims2 = shape(t % results) call h5screate_simple_f(2, dims2, dspace, hdf5_err) - call h5dcreate_f(temp_group, "values", hdf5_tallyscore_t, & + call h5dcreate_f(temp_group, "results", hdf5_tallyresult_t, & dspace, dset, hdf5_err) - f_ptr = c_loc(t % scores(1, 1)) - CALL h5dwrite_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err) + f_ptr = c_loc(t % results(1, 1)) + CALL h5dwrite_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) call h5dclose_f(dset, hdf5_err) call h5sclose_f(dspace, hdf5_err) ! Close group for the i-th tally call h5gclose_f(temp_group, hdf5_err) - end do TALLY_SCORES + end do TALLY_RESULTS else ! Indicate that tallies are off call hdf5_write_integer(tallies_group, "tallies_present", 0) @@ -1132,8 +1132,8 @@ contains ! Read batch number to restart at call hdf5_read_integer(hdf5_state_point, "current_batch", restart_batch) - ! Read information specific to criticality run - if (mode == MODE_CRITICALITY) then + ! Read information specific to eigenvalue run + if (mode == MODE_EIGENVALUE) then call hdf5_read_integer(hdf5_state_point, "n_inactive", n_inactive) call hdf5_read_integer(hdf5_state_point, "gen_per_batch", gen_per_batch) @@ -1153,7 +1153,7 @@ contains ! Read global tallies f_ptr = c_loc(global_tallies(1)) - call h5dread_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err) + call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) ! Close global tallies dataset call h5dclose_f(dset, hdf5_err) @@ -1167,14 +1167,14 @@ contains call hdf5_read_integer(tally_group, "n_realizations", & tallies(i) % n_realizations) - ! Open dataset for tally values - call h5dopen_f(tally_group, "values", dset, hdf5_err) + ! Open dataset for tally results + call h5dopen_f(tally_group, "results", dset, hdf5_err) ! Read sum and sum_sq for each tally bin - f_ptr = c_loc(tallies(i) % scores(1,1)) - call h5dread_f(dset, hdf5_tallyscore_t, f_ptr, hdf5_err) + f_ptr = c_loc(tallies(i) % results(1,1)) + call h5dread_f(dset, hdf5_tallyresult_t, f_ptr, hdf5_err) - ! Close dataset for tally values + ! Close dataset for tally results call h5dclose_f(dset, hdf5_err) ! Close tally group @@ -1191,7 +1191,7 @@ contains f_ptr = c_loc(source_bank(1)) call h5dread_f(dset, hdf5_bank_t, f_ptr, hdf5_err) - ! Close dataset for tally values + ! Close dataset for source bank call h5dclose_f(dset, hdf5_err) ! Close HDF5 state point file diff --git a/src/initialize.F90 b/src/initialize.F90 index cdd85dfdb8..3eb2a156e8 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -120,7 +120,7 @@ contains ! Allocate banks and create source particles -- for a fixed source ! calculation, the external source distribution is sampled during the ! run, not at initialization - if (run_mode == MODE_CRITICALITY) then + if (run_mode == MODE_EIGENVALUE) then call allocate_banks() if (.not. restart_run) call initialize_source() end if @@ -169,14 +169,14 @@ contains integer :: bank_types(4) ! Datatypes integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements integer :: temp_type ! temporary derived type - integer :: score_blocks(1) ! Count for each datatype - integer :: score_types(1) ! Datatypes - integer(MPI_ADDRESS_KIND) :: score_disp(1) ! Displacements - integer(MPI_ADDRESS_KIND) :: score_base_disp ! Base displacement - integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyScore - integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyScore + integer :: result_blocks(1) ! Count for each datatype + integer :: result_types(1) ! Datatypes + integer(MPI_ADDRESS_KIND) :: result_disp(1) ! Displacements + integer(MPI_ADDRESS_KIND) :: result_base_disp ! Base displacement + integer(MPI_ADDRESS_KIND) :: lower_bound ! Lower bound for TallyResult + integer(MPI_ADDRESS_KIND) :: extent ! Extent for TallyResult type(Bank) :: b - type(TallyScore) :: ts + type(TallyResult) :: tr ! Indicate that MPI is turned on mpi_enabled = .true. @@ -215,29 +215,29 @@ contains call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) ! ========================================================================== - ! CREATE MPI_TALLYSCORE TYPE + ! CREATE MPI_TALLYRESULT TYPE ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(ts % value, score_base_disp, mpi_err) - call MPI_GET_ADDRESS(ts % sum, score_disp(1), mpi_err) + call MPI_GET_ADDRESS(tr % value, result_base_disp, mpi_err) + call MPI_GET_ADDRESS(tr % sum, result_disp(1), mpi_err) ! Adjust displacements - score_disp = score_disp - score_base_disp + result_disp = result_disp - result_base_disp - ! Define temporary type for tallyscore - score_blocks = (/ 2 /) - score_types = (/ MPI_REAL8 /) - call MPI_TYPE_CREATE_STRUCT(1, score_blocks, score_disp, score_types, & + ! Define temporary type for TallyResult + result_blocks = (/ 2 /) + result_types = (/ MPI_REAL8 /) + call MPI_TYPE_CREATE_STRUCT(1, result_blocks, result_disp, result_types, & temp_type, mpi_err) ! Adjust lower-bound and extent of type for tally score lower_bound = 0 - extent = score_disp(1) + 16 + extent = result_disp(1) + 16 call MPI_TYPE_CREATE_RESIZED(temp_type, lower_bound, extent, & - MPI_TALLYSCORE, mpi_err) + MPI_TALLYRESULT, mpi_err) ! Commit derived type for tally scores - call MPI_TYPE_COMMIT(MPI_TALLYSCORE, mpi_err) + call MPI_TYPE_COMMIT(MPI_TALLYRESULT, mpi_err) end subroutine initialize_mpi #endif diff --git a/src/input_xml.F90 b/src/input_xml.F90 index b511d171e9..b833e124c0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -19,7 +19,7 @@ module input_xml implicit none type(DictionaryII), pointer :: & ! used to count how many cells each - cells_in_univ_dict => null() ! universe contains + cells_in_univ_dict => null() ! universe contains contains @@ -92,7 +92,7 @@ contains call get_environment_variable("CROSS_SECTIONS", env_variable) if (len_trim(env_variable) == 0) then message = "No cross_sections.xml file was specified in settings.xml & - &or in the CROSS_SECTIONS environment variable." + &or in the CROSS_SECTIONS environment variable." call fatal_error() else path_cross_sections = trim(env_variable) @@ -102,32 +102,41 @@ contains end if ! Make sure that either criticality or fixed source was specified - if (criticality_ % batches == 0 .and. fixed_source_ % batches == 0) then - message = "Number of batches on or " & - // "tag was zero." + if (eigenvalue_ % batches == 0 .and. fixed_source_ % batches == 0 & + .and. criticality_ % batches == 0) then + message = "Number of batches on or & + &tag was zero." call fatal_error() end if - ! Criticality information + ! Check for old tag if (criticality_ % batches > 0) then + eigenvalue_ = criticality_ + message = "The element has been deprecated and & + &replaced by ." + call warning() + end if + + ! Eigenvalue information + if (eigenvalue_ % batches > 0) then ! Set run mode - if (run_mode == NONE) run_mode = MODE_CRITICALITY + if (run_mode == NONE) run_mode = MODE_EIGENVALUE ! Check number of particles - if (len_trim(criticality_ % particles) == 0) then + if (len_trim(eigenvalue_ % particles) == 0) then message = "Need to specify number of particles per cycles." call fatal_error() end if ! If the number of particles was specified as a command-line argument, we ! don't set it here - if (n_particles == 0) n_particles = str_to_int(criticality_ % particles) + if (n_particles == 0) n_particles = str_to_int(eigenvalue_ % particles) ! Copy batch and generation information - n_batches = criticality_ % batches - n_inactive = criticality_ % inactive + n_batches = eigenvalue_ % batches + n_inactive = eigenvalue_ % inactive n_active = n_batches - n_inactive - gen_per_batch = criticality_ % generations_per_batch + gen_per_batch = eigenvalue_ % generations_per_batch ! Allocate array for batch keff and entropy allocate(k_batch(n_batches)) @@ -168,9 +177,6 @@ contains call fatal_error() end if - ! Turn on tallies if no inactive batches - if (n_inactive == 0) tallies_on = .true. - ! Copy random number seed if specified if (seed_ > 0) seed = seed_ @@ -202,7 +208,7 @@ contains inquire(FILE=path_source, EXIST=file_exists) if (.not. file_exists) then message = "Binary source file '" // trim(path_source) // & - "' does not exist!" + "' does not exist!" call fatal_error() end if @@ -221,7 +227,7 @@ contains coeffs_reqd = 3 case default message = "Invalid spatial distribution for external source: " & - // trim(source_ % space % type) + // trim(source_ % space % type) call fatal_error() end select @@ -235,11 +241,11 @@ contains ! Read parameters for spatial distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for spatial " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for spatial " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_space(n)) @@ -266,7 +272,7 @@ contains external_source % type_angle = SRC_ANGLE_TABULAR case default message = "Invalid angular distribution for external source: " & - // trim(source_ % angle % type) + // trim(source_ % angle % type) call fatal_error() end select @@ -280,11 +286,11 @@ contains ! Read parameters for angle distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for angle " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for angle " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_angle(n)) @@ -314,7 +320,7 @@ contains external_source % type_energy = SRC_ENERGY_TABULAR case default message = "Invalid energy distribution for external source: " & - // trim(source_ % energy % type) + // trim(source_ % energy % type) call fatal_error() end select @@ -328,11 +334,11 @@ contains ! Read parameters for energy distribution if (n < coeffs_reqd) then message = "Not enough parameters specified for energy " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > coeffs_reqd) then message = "Too many parameters specified for energy " & - // "distribution of external source." + // "distribution of external source." call fatal_error() elseif (n > 0) then allocate(external_source % params_energy(n)) @@ -370,10 +376,10 @@ contains ! Check to make sure enough values were supplied if (size(entropy_(1) % lower_left) /= 3) then message = "Need to specify (x,y,z) coordinates of lower-left corner & - &of Shannon entropy mesh." + &of Shannon entropy mesh." elseif (size(entropy_(1) % upper_right) /= 3) then message = "Need to specify (x,y,z) coordinates of upper-right corner & - &of Shannon entropy mesh." + &of Shannon entropy mesh." end if ! Allocate mesh object and coordinates on mesh @@ -388,7 +394,7 @@ contains ! Check on values provided if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) then message = "Upper-right coordinate must be greater than lower-left & - &coordinate for Shannon entropy mesh." + &coordinate for Shannon entropy mesh." call fatal_error() end if @@ -399,7 +405,7 @@ contains ! If so, make sure proper number of values were given if (size(entropy_(1) % dimension) /= 3) then message = "Dimension of entropy mesh must be given as three & - &integers." + &integers." call fatal_error() end if @@ -420,13 +426,13 @@ contains ! Check to make sure enough values were supplied if (size(uniform_fs_(1) % lower_left) /= 3) then message = "Need to specify (x,y,z) coordinates of lower-left corner & - &of UFS mesh." + &of UFS mesh." elseif (size(uniform_fs_(1) % upper_right) /= 3) then message = "Need to specify (x,y,z) coordinates of upper-right corner & - &of UFS mesh." + &of UFS mesh." elseif (size(uniform_fs_(1) % dimension) /= 3) then message = "Dimension of UFS mesh must be given as three & - &integers." + &integers." call fatal_error() end if @@ -450,13 +456,13 @@ contains ! Check on values provided if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then message = "Upper-right coordinate must be greater than lower-left & - &coordinate for UFS mesh." + &coordinate for UFS mesh." call fatal_error() end if ! Calculate width ufs_mesh % width = (ufs_mesh % upper_right - & - ufs_mesh % lower_left) / ufs_mesh % dimension + ufs_mesh % lower_left) / ufs_mesh % dimension ! Calculate volume fraction of each cell ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8) @@ -466,7 +472,7 @@ contains ! Allocate source_frac allocate(source_frac(1, ufs_mesh % dimension(1), & - ufs_mesh % dimension(2), ufs_mesh % dimension(3))) + ufs_mesh % dimension(2), ufs_mesh % dimension(3))) end if ! Check if the user has specified to write state points @@ -488,7 +494,7 @@ contains n_state_points = n_batches / state_point_(1) % interval allocate(statepoint_batch(n_state_points)) statepoint_batch = (/ (state_point_(1) % interval * i, i = 1, & - n_state_points) /) + n_state_points) /) else ! If neither were specified, write state point at last batch n_state_points = 1 @@ -498,7 +504,7 @@ contains ! Check if the user has specified to write binary source file if (trim(state_point_(1) % source_separate) == 'on') & - source_separate = .true. + source_separate = .true. else ! If no tag was present, by default write state point at ! last batch only @@ -629,7 +635,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 " // & - trim(to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -643,7 +649,7 @@ contains ! Check to make sure that surfaces were specified if (.not. associated(cell_(i) % surfaces)) then message = "No surfaces specified for cell " // & - trim(to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -659,7 +665,7 @@ contains ! another universe if (c % fill == NONE) then message = "Cannot apply a rotation to cell " // trim(to_str(& - c % id)) // " because it is not filled with another universe" + c % id)) // " because it is not filled with another universe" call fatal_error() end if @@ -667,7 +673,7 @@ contains n = size(cell_(i) % rotation) if (n /= 3) then message = "Incorrect number of rotation parameters on cell " // & - to_str(c % id) + to_str(c % id) call fatal_error() end if @@ -679,13 +685,13 @@ contains ! Calculate rotation matrix based on angles given allocate(c % rotation(3,3)) c % rotation = reshape((/ & - cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & - -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & - cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & - sin(phi)*cos(theta), & - sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), & - -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), & - cos(phi)*cos(theta) /), (/ 3,3 /)) + cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & + -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & + cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & + sin(phi)*cos(theta), & + sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), & + -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), & + cos(phi)*cos(theta) /), (/ 3,3 /)) end if ! Translation vector @@ -694,7 +700,7 @@ contains ! another universe if (c % fill == NONE) then message = "Cannot apply a translation to cell " // trim(to_str(& - c % id)) // " because it is not filled with another universe" + c % id)) // " because it is not filled with another universe" call fatal_error() end if @@ -702,7 +708,7 @@ contains n = size(cell_(i) % translation) if (n /= 3) then message = "Incorrect number of translation parameters on cell " & - // to_str(c % id) + // to_str(c % id) call fatal_error() end if @@ -814,11 +820,11 @@ contains n = size(surface_(i) % coeffs) if (n < coeffs_reqd) then message = "Not enough coefficients specified for surface: " // & - trim(to_str(s % id)) + trim(to_str(s % id)) call fatal_error() elseif (n > coeffs_reqd) then message = "Too many coefficients specified for surface: " // & - trim(to_str(s % id)) + trim(to_str(s % id)) call fatal_error() else allocate(s % coeffs(n)) @@ -839,7 +845,7 @@ contains s % bc = BC_PERIODIC case default message = "Unknown boundary condition '" // trim(word) // & - "' specified on surface " // trim(to_str(s % id)) + "' specified on surface " // trim(to_str(s % id)) call fatal_error() end select @@ -888,7 +894,7 @@ contains ! Read lattice lower-left location if (size(lattice_(i) % dimension) /= size(lattice_(i) % lower_left)) then message = "Number of entries on must be the same as " // & - "the number of entries on ." + "the number of entries on ." call fatal_error() end if l % x0 = lattice_(i) % lower_left(1) @@ -897,7 +903,7 @@ contains ! Read lattice widths if (size(lattice_(i) % width) /= size(lattice_(i) % lower_left)) then message = "Number of entries on must be the same as " // & - "the number of entries on ." + "the number of entries on ." call fatal_error() end if l % width_x = lattice_(i) % width(1) @@ -998,7 +1004,7 @@ contains sum_density = .false. if (val <= ZERO) then message = "Need to specify a positive density on material " // & - trim(to_str(mat % id)) // "." + trim(to_str(mat % id)) // "." call fatal_error() end if @@ -1015,7 +1021,7 @@ contains mat % density = 1.0e-24 * val case default message = "Unkwown units '" // trim(material_(i) % density % units) & - // "' specified on material " // trim(to_str(mat % id)) + // "' specified on material " // trim(to_str(mat % id)) call fatal_error() end select end if @@ -1026,7 +1032,7 @@ contains ! Check to ensure material has at least one nuclide if (.not. associated(material_(i) % nuclides)) then message = "No nuclides specified on material " // & - trim(to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -1044,7 +1050,7 @@ contains ! Check for empty name on nuclide if (len_trim(nuc % name) == 0) then message = "No name specified on nuclide in material " // & - trim(to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -1052,7 +1058,7 @@ contains if (len_trim(nuc % xs) == 0) then if (default_xs == '') then message = "No cross section specified for nuclide in material " & - // trim(to_str(mat % id)) + // trim(to_str(mat % id)) call fatal_error() else nuc % xs = default_xs @@ -1066,7 +1072,7 @@ contains ! Check that this nuclide is listed in the cross_sections.xml file if (.not. dict_has_key(xs_listing_dict, name)) then message = "Could not find nuclide " // trim(name) // & - " in cross_sections.xml file!" + " in cross_sections.xml file!" call fatal_error() end if @@ -1074,7 +1080,7 @@ contains n = len_trim(name) if (name(n:n) /= 'c') then message = "Cross-section table " // trim(name) // & - " is not a continuous-energy neutron table." + " is not a continuous-energy neutron table." call fatal_error() end if @@ -1099,11 +1105,11 @@ contains ! weight percents were specified if (nuc % ao == ZERO .and. nuc % wo == ZERO) then message = "No atom or weight percent specified for nuclide " // & - trim(name) + trim(name) call fatal_error() elseif (nuc % ao /= ZERO .and. nuc % wo /= ZERO) then message = "Cannot specify both atom and weight percents for a & - &nuclide: " // trim(name) + &nuclide: " // trim(name) call fatal_error() end if @@ -1118,9 +1124,9 @@ contains ! Check to make sure either all atom percents or all weight percents are ! given if (.not. (all(mat % atom_density > ZERO) .or. & - all(mat % atom_density < ZERO))) then + all(mat % atom_density < ZERO))) then message = "Cannot mix atom and weight percents in material " // & - to_str(mat % id) + to_str(mat % id) call fatal_error() end if @@ -1141,7 +1147,7 @@ contains ! Check that this nuclide is listed in the cross_sections.xml file if (.not. dict_has_key(xs_listing_dict, name)) then message = "Could not find S(a,b) table " // trim(name) // & - " in cross_sections.xml file!" + " in cross_sections.xml file!" call fatal_error() end if mat % has_sab_table = .true. @@ -1304,7 +1310,7 @@ contains ! Check that dimensions are all greater than zero if (any(mesh_(i) % dimension <= 0)) then message = "All entries on the element for a tally mesh & - &must be positive." + &must be positive." call fatal_error() end if @@ -1314,24 +1320,24 @@ contains ! Read mesh lower-left corner location if (m % n_dimension /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as & - &the number of entries on ." + &the number of entries on ." call fatal_error() end if m % lower_left = mesh_(i) % lower_left ! Make sure either upper-right or width was specified if (associated(mesh_(i) % upper_right) .and. & - associated(mesh_(i) % width)) then + associated(mesh_(i) % width)) then message = "Cannot specify both and on a & - &tally mesh." + &tally mesh." call fatal_error() end if ! Make sure either upper-right or width was specified if (.not. associated(mesh_(i) % upper_right) .and. & - .not. associated(mesh_(i) % width)) then + .not. associated(mesh_(i) % width)) then message = "Must specify either and on a & - &tally mesh." + &tally mesh." call fatal_error() end if @@ -1339,7 +1345,7 @@ contains ! Check to ensure width has same dimensions if (size(mesh_(i) % width) /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as the & - &number of entries on ." + &number of entries on ." call fatal_error() end if @@ -1357,14 +1363,14 @@ contains ! Check to ensure width has same dimensions if (size(mesh_(i) % upper_right) /= size(mesh_(i) % lower_left)) then message = "Number of entries on must be the same as & - &the number of entries on ." + &the number of entries on ." call fatal_error() end if ! Check that upper-right is above lower-left if (any(mesh_(i) % upper_right < mesh_(i) % lower_left)) then message = "The coordinates must be greater than the & - & coordinates on a tally mesh." + & coordinates on a tally mesh." call fatal_error() end if @@ -1413,7 +1419,7 @@ contains if (size(tally_(i) % filters) > 0) then message = "Tally filters should be specified with multiple & - &elements. Did you forget to change your element?" + &elements. Did you forget to change your element?" call fatal_error() end if @@ -1445,7 +1451,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('cellborn') @@ -1459,7 +1465,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('material') @@ -1473,7 +1479,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('universe') @@ -1487,7 +1493,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('surface') @@ -1501,7 +1507,7 @@ contains allocate(t % filters(j) % int_bins(n_words)) do k = 1, n_words t % filters(j) % int_bins(k) = int(str_to_int(& - tally_(i) % filter(j) % bins(k)),4) + tally_(i) % filter(j) % bins(k)),4) end do case ('mesh') @@ -1523,7 +1529,7 @@ contains m => meshes(i_mesh) else message = "Could not find mesh " // trim(to_str(id)) // & - " specified on tally " // trim(to_str(t % id)) + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if @@ -1547,7 +1553,7 @@ contains allocate(t % filters(j) % real_bins(n_words)) do k = 1, n_words t % filters(j) % real_bins(k) = str_to_real(& - tally_(i) % filter(j) % bins(k)) + tally_(i) % filter(j) % bins(k)) end do case ('energyout') @@ -1561,7 +1567,7 @@ contains allocate(t % filters(j) % real_bins(n_words)) do k = 1, n_words t % filters(j) % real_bins(k) = str_to_real(& - tally_(i) % filter(j) % bins(k)) + tally_(i) % filter(j) % bins(k)) end do end select @@ -1575,9 +1581,9 @@ contains ! Check that both cell and surface weren't specified if (t % find_filter(FILTER_CELL) > 0 .and. & - t % find_filter(FILTER_SURFACE) > 0) then + t % find_filter(FILTER_SURFACE) > 0) then message = "Cannot specify both cell and surface filters for tally " & - // trim(to_str(t % id)) + // trim(to_str(t % id)) call fatal_error() end if @@ -1596,7 +1602,7 @@ contains ! Set bins to 1, 2, 3, ..., n_nuclides_total, -1 t % nuclide_bins(1:n_nuclides_total) = & - (/ (j, j=1, n_nuclides_total) /) + (/ (j, j=1, n_nuclides_total) /) t % nuclide_bins(n_nuclides_total + 1) = -1 ! Set number of nuclide bins @@ -1624,7 +1630,7 @@ contains key_list => dict_keys(nuclide_dict) do while (associated(key_list)) if (starts_with(key_list % data % key, & - tally_(i) % nuclides(j))) then + tally_(i) % nuclides(j))) then word = key_list % data % key exit end if @@ -1636,8 +1642,8 @@ contains ! Check if no nuclide was found if (.not. associated(key_list)) then message = "Could not find the nuclide " // trim(& - tally_(i) % nuclides(j)) // " specified in tally " & - // trim(to_str(t % id)) // " in any material." + tally_(i) % nuclides(j)) // " specified in tally " & + // trim(to_str(t % id)) // " in any material." call fatal_error() end if deallocate(key_list) @@ -1650,7 +1656,7 @@ contains ! Check to make sure nuclide specified is in problem if (.not. dict_has_key(nuclide_dict, word)) then message = "The nuclide " // trim(word) // " from tally " // & - trim(to_str(t % id)) // " is not present in any material." + trim(to_str(t % id)) // " is not present in any material." call fatal_error() end if @@ -1749,7 +1755,7 @@ contains case ('flux') ! Prohibit user from tallying flux for an individual nuclide if (.not. (t % n_nuclide_bins == 1 .and. & - t % nuclide_bins(1) == -1)) then + t % nuclide_bins(1) == -1)) then message = "Cannot tally flux for an individual nuclide." call fatal_error() end if @@ -1763,7 +1769,7 @@ contains t % score_bins(j) = SCORE_TOTAL if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally total reaction rate with an & - &outgoing energy filter." + &outgoing energy filter." call fatal_error() end if case ('scatter') @@ -1824,14 +1830,14 @@ contains t % score_bins(j) = SCORE_ABSORPTION if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally absorption rate with an outgoing & - &energy filter." + &energy filter." call fatal_error() end if case ('fission') t % score_bins(j) = SCORE_FISSION if (t % find_filter(FILTER_ENERGYOUT) > 0) then message = "Cannot tally fission rate with an outgoing & - &energy filter." + &energy filter." call fatal_error() end if case ('nu-fission') @@ -1844,8 +1850,8 @@ contains ! for this tally if (n_words > 1) then message = "Cannot tally other scoring functions in the same & - &tally as surface currents. Separate other scoring & - &functions into a distinct tally." + &tally as surface currents. Separate other scoring & + &functions into a distinct tally." call fatal_error() end if @@ -1877,13 +1883,13 @@ contains t % filters(t % n_filters) % type = FILTER_SURFACE t % filters(t % n_filters) % n_bins = 2 * m % n_dimension allocate(t % filters(t % n_filters) % int_bins(& - 2 * m % n_dimension)) + 2 * m % n_dimension)) if (m % n_dimension == 2) then t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT /) + OUT_RIGHT, IN_FRONT, OUT_FRONT /) elseif (m % n_dimension == 3) then t % filters(t % n_filters) % int_bins = (/ IN_RIGHT, & - OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /) + OUT_RIGHT, IN_FRONT, OUT_FRONT, IN_TOP, OUT_TOP /) end if t % find_filter(FILTER_SURFACE) = t % n_filters @@ -1892,15 +1898,14 @@ contains case default message = "Unknown scoring function: " // & - trim(tally_(i) % scores(l)) + trim(tally_(i) % scores(j)) call fatal_error() end select end do - t % n_score_bins = n_scores - t % n_nonPN_score_bins = n_words + t % n_score_bins = n_words else message = "No specified on tally " // trim(to_str(t % id)) & - // "." + // "." call fatal_error() end if @@ -1918,7 +1923,7 @@ contains ! tally needs post-collision information if (t % estimator == ESTIMATOR_ANALOG) then message = "Cannot use track-length estimator for tally " & - // to_str(t % id) + // to_str(t % id) call fatal_error() end if @@ -1927,7 +1932,7 @@ contains case default message = "Invalid estimator '" // trim(tally_(i) % estimator) & - // "' on tally " // to_str(t % id) + // "' on tally " // to_str(t % id) call fatal_error() end select end if @@ -2031,7 +2036,7 @@ contains ! Set output file path pl % path_plot = trim(path_input) // trim(to_str(pl % id)) // & - "_" // trim(plot_(i) % filename) // ".ppm" + "_" // trim(plot_(i) % filename) // ".ppm" ! Copy plot pixel size if (size(plot_(i) % pixels) == 2) then @@ -2047,7 +2052,7 @@ contains pl % not_found % rgb = plot_(i) % background else message = "Bad background RGB " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if else @@ -2062,7 +2067,7 @@ contains ! pl % type = PLOT_TYPE_POINTS case default message = "Unsupported plot type '" // plot_(i) % type & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2076,7 +2081,7 @@ contains pl % basis = PLOT_BASIS_YZ case default message = "Unsupported plot basis '" // plot_(i) % basis & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2085,7 +2090,7 @@ contains pl % origin = plot_(i) % origin else message = "Origin must be length 3 " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2097,7 +2102,7 @@ contains pl % width(2) = plot_(i) % width(2) else message = "Bad plot width " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2125,7 +2130,7 @@ contains case default message = "Unsupported plot color type '" // plot_(i) % color & - // "' in plot " // trim(to_str(pl % id)) + // "' in plot " // trim(to_str(pl % id)) call fatal_error() end select @@ -2135,7 +2140,7 @@ contains do j = 1, n_cols if (size(plot_(i) % col_spec_(j) % rgb) /= 3) then message = "Bad RGB " & - // "in plot " // trim(to_str(pl % id)) + // "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2147,7 +2152,7 @@ contains pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find cell " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2157,7 +2162,7 @@ contains pl % colors(col_id) % rgb = plot_(i) % col_spec_(j) % rgb else message = "Could not find material " // trim(to_str(col_id)) // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -2169,7 +2174,7 @@ contains if (associated(plot_(i) % mask_)) then if (size(plot_(i) % mask_) > 1) then message = "Mutliple masks" // & - " specified in plot " // trim(to_str(pl % id)) + " specified in plot " // trim(to_str(pl % id)) call fatal_error() else if (.not. size(plot_(i) % mask_) == 0) then do j=1,size(pl % colors) @@ -2204,10 +2209,10 @@ contains ! Check if cross_sections.xml exists inquire(FILE=path_cross_sections, EXIST=file_exists) if (.not. file_exists) then - ! Could not find cross_sections.xml file - message = "Cross sections XML file '" // trim(path_cross_sections) // & + ! Could not find cross_sections.xml file + message = "Cross sections XML file '" // trim(path_cross_sections) // & "' does not exist!" - call fatal_error() + call fatal_error() end if message = "Reading cross sections XML file..." @@ -2223,25 +2228,25 @@ contains call read_xml_file_cross_sections_t(path_cross_sections) if (len_trim(directory_) > 0) then - ! Copy directory information if present - directory = trim(directory_) + ! Copy directory information if present + directory = trim(directory_) else - ! If no directory is listed in cross_sections.xml, by default select the - ! directory in which the cross_sections.xml file resides - i = index(path_cross_sections, "/", BACK=.true.) - directory = path_cross_sections(1:i) + ! If no directory is listed in cross_sections.xml, by default select the + ! directory in which the cross_sections.xml file resides + i = index(path_cross_sections, "/", BACK=.true.) + directory = path_cross_sections(1:i) end if ! determine whether binary/ascii if (filetype_ == 'ascii') then - filetype = ASCII + filetype = ASCII elseif (filetype_ == 'binary') then - filetype = BINARY + filetype = BINARY elseif (len_trim(filetype_) == 0) then - filetype = ASCII + filetype = ASCII else - message = "Unknown filetype in cross_sections.xml: " // trim(filetype_) - call fatal_error() + message = "Unknown filetype in cross_sections.xml: " // trim(filetype_) + call fatal_error() end if ! copy default record length and entries for binary files @@ -2250,57 +2255,57 @@ contains ! Allocate xs_listings array if (.not. associated(ace_tables_)) then - message = "No ACE table listings present in cross_sections.xml file!" - call fatal_error() + message = "No ACE table listings present in cross_sections.xml file!" + call fatal_error() else - n_listings = size(ace_tables_) - allocate(xs_listings(n_listings)) + n_listings = size(ace_tables_) + allocate(xs_listings(n_listings)) end if do i = 1, n_listings - listing => xs_listings(i) + listing => xs_listings(i) - ! copy a number of attributes - listing % name = trim(ace_tables_(i) % name) - listing % alias = trim(ace_tables_(i) % alias) - listing % zaid = ace_tables_(i) % zaid - listing % awr = ace_tables_(i) % awr - listing % kT = ace_tables_(i) % temperature - listing % location = ace_tables_(i) % location + ! copy a number of attributes + listing % name = trim(ace_tables_(i) % name) + listing % alias = trim(ace_tables_(i) % alias) + listing % zaid = ace_tables_(i) % zaid + listing % awr = ace_tables_(i) % awr + listing % kT = ace_tables_(i) % temperature + listing % location = ace_tables_(i) % location - ! determine type of cross section - if (ends_with(listing % name, 'c')) then + ! determine type of cross section + if (ends_with(listing % name, 'c')) then listing % type = ACE_NEUTRON - elseif (ends_with(listing % name, 't')) then + elseif (ends_with(listing % name, 't')) then listing % type = ACE_THERMAL - end if + end if - ! set filetype, record length, and number of entries - listing % filetype = filetype - listing % recl = recl - listing % entries = entries + ! set filetype, record length, and number of entries + listing % filetype = filetype + listing % recl = recl + listing % entries = entries - ! determine metastable state - if (ace_tables_(i) % metastable == 0) then + ! determine metastable state + if (ace_tables_(i) % metastable == 0) then listing % metastable = .false. - else + else listing % metastable = .true. - end if + end if - ! determine path of cross section table - if (starts_with(ace_tables_(i) % path, '/')) then + ! determine path of cross section table + if (starts_with(ace_tables_(i) % path, '/')) then listing % path = ace_tables_(i) % path - else + else if (ends_with(directory,'/')) then - listing % path = trim(directory) // trim(ace_tables_(i) % path) + listing % path = trim(directory) // trim(ace_tables_(i) % path) else - listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path) + listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path) end if - end if + end if - ! create dictionary entry for both name and alias - call dict_add_key(xs_listing_dict, listing % name, i) - call dict_add_key(xs_listing_dict, listing % alias, i) + ! create dictionary entry for both name and alias + call dict_add_key(xs_listing_dict, listing % name, i) + call dict_add_key(xs_listing_dict, listing % alias, i) end do end subroutine read_cross_sections_xml diff --git a/src/main.F90 b/src/main.F90 index 89cc7e8570..c5875d5d0d 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -1,7 +1,7 @@ program main use constants - use criticality, only: run_criticality + use eigenvalue, only: run_eigenvalue use finalize, only: finalize_run use fixed_source, only: run_fixedsource use global @@ -17,8 +17,8 @@ program main select case (run_mode) case (MODE_FIXEDSOURCE) call run_fixedsource() - case (MODE_CRITICALITY) - call run_criticality() + case (MODE_EIGENVALUE) + call run_eigenvalue() case (MODE_PLOTTING) call run_plot() case (MODE_TALLIES) diff --git a/src/mesh.F90 b/src/mesh.F90 index 1c80a94e08..cc3734d430 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -93,7 +93,7 @@ contains !=============================================================================== ! MESH_INDICES_TO_BIN maps (i,j) or (i,j,k) indices to a single bin number for -! use in a TallyObject scores array +! use in a TallyObject results array !=============================================================================== function mesh_indices_to_bin(m, ijk, surface_current) result(bin) @@ -126,7 +126,7 @@ contains end function mesh_indices_to_bin !=============================================================================== -! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject scores array to +! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to ! (i,j) or (i,j,k) indices !=============================================================================== diff --git a/src/output.F90 b/src/output.F90 index 530cd42fb3..c7f2b3ff2d 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -32,35 +32,35 @@ contains subroutine title() write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & - ' .d88888b. 888b d888 .d8888b.', & - ' d88P" "Y88b 8888b d8888 d88P Y88b', & - ' 888 888 88888b.d88888 888 888', & - ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & - ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & - ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & - ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & - ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & - '__________________888______________________________________________________', & - ' 888', & - ' 888' + ' .d88888b. 888b d888 .d8888b.', & + ' d88P" "Y88b 8888b d8888 d88P Y88b', & + ' 888 888 88888b.d88888 888 888', & + ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & + ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & + ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & + ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & + ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & + '__________________888______________________________________________________', & + ' 888', & + ' 888' ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Developed At: Massachusetts Institute of Technology' + ' Developed At: Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & - VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1 #endif ! Write the date and time write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') & - time_stamp() + time_stamp() #ifdef MPI ! Write number of processors write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') & - trim(to_str(n_procs)) + trim(to_str(n_procs)) #endif end subroutine title @@ -77,7 +77,7 @@ contains call date_and_time(DATE=date_, TIME=time_) current_time = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) // & - " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) + " " // time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) end function time_stamp @@ -125,13 +125,13 @@ contains select case (header_level) case (1) write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & - repeat('=', n) // '> ' // trim(line) // ' <' // & - repeat('=', m), repeat('=', 75) + repeat('=', n) // '> ' // trim(line) // ' <' // & + repeat('=', m), repeat('=', 75) case (2) write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75) case (3) write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // & - trim(line) // ' <' // repeat('=', m) + trim(line) // ' <' // repeat('=', m) end select end subroutine header @@ -145,11 +145,11 @@ contains if (master) then write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') & - "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2012 & - &Massachusetts Institute of Technology" + &Massachusetts Institute of Technology" write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at & - &" + &" end if end subroutine print_version @@ -207,7 +207,7 @@ contains else ! Determine last space in current line i_end = i_start + index(message(i_start+1:i_start+line_wrap), & - ' ', BACK=.true.) + ' ', BACK=.true.) ! Write up to last space write(ou, fmt='(1X,A)') message(i_start+1:i_end-1) @@ -270,8 +270,8 @@ contains l => lattices(coord % lattice) write(ou,*) ' Lattice = ' // trim(to_str(l % id)) write(ou,*) ' Lattice position = (' // trim(to_str(& - p % coord % lattice_x)) // ',' // trim(to_str(& - p % coord % lattice_y)) // ')' + p % coord % lattice_x)) // ',' // trim(to_str(& + p % coord % lattice_y)) // ')' end if ! Print local coordinates @@ -390,7 +390,7 @@ contains case default index_surf = abs(c % surfaces(i)) string = trim(string) // ' ' // to_str(sign(& - surfaces(index_surf) % id, c % surfaces(i))) + surfaces(index_surf) % id, c % surfaces(i))) end select end do write(unit_,*) ' Surface Specification:' // trim(string) @@ -484,6 +484,7 @@ contains integer :: i ! loop index for coefficients integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string + type(Cell), pointer :: c => null() ! set default unit if not specified if (present(unit)) then @@ -539,7 +540,9 @@ contains string = "" if (allocated(surf % neighbor_pos)) then do i = 1, size(surf % neighbor_pos) - string = trim(string) // ' ' // to_str(surf % neighbor_pos(i)) + c => cells(abs(surf % neighbor_pos(i))) + string = trim(string) // ' ' // to_str(& + sign(c % id, surf % neighbor_pos(i))) end do end if write(unit_,*) ' Positive Neighbors = ' // trim(string) @@ -548,7 +551,9 @@ contains string = "" if (allocated(surf % neighbor_neg)) then do i = 1, size(surf % neighbor_neg) - string = trim(string) // ' ' // to_str(surf % neighbor_neg(i)) + c => cells(abs(surf % neighbor_neg(i))) + string = trim(string) // ' ' // to_str(& + sign(c % id, surf % neighbor_neg(i))) end do end if write(unit_,*) ' Negative Neighbors =' // trim(string) @@ -595,7 +600,7 @@ contains ! Write total atom density in atom/b-cm write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) & - // ' atom/b-cm' + // ' atom/b-cm' ! Write atom density for each nuclide in material write(unit_,*) ' Nuclides:' @@ -603,7 +608,7 @@ contains nuc => nuclides(mat % nuclide(i)) density = mat % atom_density(i) string = ' ' // trim(nuc % name) // ' = ' // & - trim(to_str(density)) // ' atom/b-cm' + trim(to_str(density)) // ' atom/b-cm' write(unit_,*) trim(string) end do @@ -742,7 +747,7 @@ contains string = "" do i = 1, t % filters(j) % n_bins + 1 string = trim(string) // ' ' // trim(to_str(& - t % filters(j) % real_bins(i))) + t % filters(j) % real_bins(i))) end do write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if @@ -753,7 +758,7 @@ contains string = "" do i = 1, t % filters(j) % n_bins + 1 string = trim(string) // ' ' // trim(to_str(& - t % filters(j) % real_bins(i))) + t % filters(j) % real_bins(i))) end do write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if @@ -765,10 +770,10 @@ contains write(unit_,fmt='(A)',advance='no') ' total' else write(unit_,fmt='(A)',advance='no') ' ' // trim(adjustl(& - nuclides(t % nuclide_bins(i)) % name)) + nuclides(t % nuclide_bins(i)) % name)) end if if (mod(i,4) == 0 .and. i /= t % n_nuclide_bins) & - write(unit_,'(/18X)',advance='no') + write(unit_,'(/18X)',advance='no') end do write(unit_,*) @@ -938,8 +943,8 @@ contains end if write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') & - reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & - law(1:4), rxn % threshold, size_angle, size_energy + reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & + law(1:4), rxn % threshold, size_angle, size_energy ! Accumulate data size size_xs = size_xs + (nuc % n_grid - rxn%threshold + 1) * 8 @@ -976,11 +981,11 @@ contains write(unit_,*) ' Memory Requirements' write(unit_,*) ' Cross sections = ' // trim(to_str(size_xs)) // ' bytes' write(unit_,*) ' Secondary angle distributions = ' // & - trim(to_str(size_angle_total)) // ' bytes' + trim(to_str(size_angle_total)) // ' bytes' write(unit_,*) ' Secondary energy distributions = ' // & - trim(to_str(size_energy_total)) // ' bytes' + trim(to_str(size_energy_total)) // ' bytes' write(unit_,*) ' Probability Tables = ' // & - trim(to_str(size_urr)) // ' bytes' + trim(to_str(size_urr)) // ' bytes' write(unit_,*) ' Total = ' // trim(to_str(size_total)) // ' bytes' ! Blank line at end of nuclide @@ -1016,28 +1021,28 @@ contains ! Inelastic data write(unit_,*) ' # of Incoming Energies (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_e_in)) + trim(to_str(sab % n_inelastic_e_in)) write(unit_,*) ' # of Outgoing Energies (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_e_out)) + trim(to_str(sab % n_inelastic_e_out)) write(unit_,*) ' # of Outgoing Angles (Inelastic) = ' // & - trim(to_str(sab % n_inelastic_mu)) + trim(to_str(sab % n_inelastic_mu)) write(unit_,*) ' Threshold for Inelastic = ' // & - trim(to_str(sab % threshold_inelastic)) + trim(to_str(sab % threshold_inelastic)) ! Elastic data if (sab % n_elastic_e_in > 0) then write(unit_,*) ' # of Incoming Energies (Elastic) = ' // & - trim(to_str(sab % n_elastic_e_in)) + trim(to_str(sab % n_elastic_e_in)) write(unit_,*) ' # of Outgoing Angles (Elastic) = ' // & - trim(to_str(sab % n_elastic_mu)) + trim(to_str(sab % n_elastic_mu)) write(unit_,*) ' Threshold for Elastic = ' // & - trim(to_str(sab % threshold_elastic)) + trim(to_str(sab % threshold_elastic)) end if ! Determine memory used by S(a,b) table and write out size_sab = 8 * (sab % n_inelastic_e_in * (2 + sab % n_inelastic_e_out * & - (1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * & - (2 + sab % n_elastic_mu)) + (1 + sab % n_inelastic_mu)) + sab % n_elastic_e_in * & + (2 + sab % n_elastic_mu)) write(unit_,*) ' Memory Used = ' // trim(to_str(size_sab)) // ' bytes' ! Blank line at end @@ -1065,31 +1070,31 @@ contains call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) write(UNIT=UNIT_SUMMARY, FMT=*) & - "Copyright: 2011-2012 Massachusetts Institute of Technology" + "Copyright: 2011-2012 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & - "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1 #endif write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A)') & - time_stamp() + time_stamp() ! Write information on number of processors #ifdef MPI write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') & - trim(to_str(n_procs)) + trim(to_str(n_procs)) #endif ! Display problem summary call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY) select case(run_mode) - case (MODE_CRITICALITY) - write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality' + case (MODE_EIGENVALUE) + write(UNIT_SUMMARY,100) 'Problem type:', 'k eigenvalue' write(UNIT_SUMMARY,101) 'Number of Batches:', n_batches write(UNIT_SUMMARY,101) 'Number of Inactive Batches:', n_inactive write(UNIT_SUMMARY,101) 'Generations per Batch:', gen_per_batch case (MODE_FIXEDSOURCE) - write(UNIT_SUMMARY,100) 'Problem type:', 'External Source' + write(UNIT_SUMMARY,100) 'Problem type:', 'fixed source' end select write(UNIT_SUMMARY,101) 'Number of Particles:', n_particles @@ -1122,7 +1127,7 @@ contains call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY) write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid)) write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(& - n_grid*n_nuclides_total*4)) // " bytes" + n_grid*n_nuclides_total*4)) // " bytes" ! print summary of variance reduction call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY) @@ -1232,16 +1237,16 @@ contains ! write out information batch and option independent output write(UNIT=OUTPUT_UNIT, FMT='(2X,I5)', ADVANCE='NO') current_batch write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & - k_batch(current_batch) + k_batch(current_batch) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy(current_batch) + entropy(current_batch) ! write out accumulated k-effective if after first active batch if (current_batch > n_inactive + 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & - keff, keff_std + keff, keff_std else write(UNIT=OUTPUT_UNIT, FMT='(23X)', ADVANCE='NO') end if @@ -1249,11 +1254,11 @@ contains ! write out cmfd keff if it is active if (cmfd_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % keff + cmfd % keff ! write out cmfd entopy if (cmfd_on .and. entropy_on) write(UNIT=OUTPUT_UNIT, & - FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy + FMT='(3X, F8.5)', ADVANCE='NO') cmfd % entropy ! next line write(UNIT=OUTPUT_UNIT, FMT=*) @@ -1280,18 +1285,18 @@ contains ! Write plotting origin write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // & - " " // trim(to_str(pl % origin(2))) // " " // & - trim(to_str(pl % origin(3))) + " " // trim(to_str(pl % origin(2))) // " " // & + trim(to_str(pl % origin(3))) ! Write plotting width if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Width:", trim(to_str(pl % width(1))) // & - " " // trim(to_str(pl % width(2))) + " " // trim(to_str(pl % width(2))) write(ou,100) "Coloring:", trim(to_str(pl % color_by)) write(ou,100) "Basis:", trim(to_str(pl % basis)) write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // & - trim(to_str(pl % pixels(2))) + trim(to_str(pl % pixels(2))) end if write(ou,*) @@ -1324,7 +1329,7 @@ contains write(ou,100) " Reading cross sections", time_read_xs % elapsed write(ou,100) " Unionizing energy grid", time_unionize % elapsed write(ou,100) "Total time in simulation", time_inactive % elapsed + & - time_active % elapsed + time_active % elapsed write(ou,100) " Time in transport only", time_transport % elapsed if(cmfd_run) write(ou,100) "Total CMFD time", time_cmfd % elapsed write(ou,100) " Time in inactive batches", time_inactive % elapsed @@ -1338,14 +1343,14 @@ contains if (restart_run) then total_particles = n_particles * (n_batches - & - restart_batch) * gen_per_batch + restart_batch) * gen_per_batch else total_particles = n_particles * n_batches * gen_per_batch end if ! display calculate rate speed = real(total_particles) / (time_inactive % elapsed + & - time_active % elapsed) + time_active % elapsed) string = to_str(speed) write(ou,101) "Calculation Rate", trim(string) @@ -1363,13 +1368,13 @@ contains ! write global tallies write(ou,102) "k-effective (Analog)", global_tallies(K_ANALOG) % sum, & - global_tallies(K_ANALOG) % sum_sq + global_tallies(K_ANALOG) % sum_sq write(ou,102) "k-effective (Collision)", global_tallies(K_COLLISION) % sum, & - global_tallies(K_COLLISION) % sum_sq + global_tallies(K_COLLISION) % sum_sq write(ou,102) "k-effective (Track-length)", global_tallies(K_TRACKLENGTH) % sum, & - global_tallies(K_TRACKLENGTH) % sum_sq + global_tallies(K_TRACKLENGTH) % sum_sq write(ou,102) "Leakage Fraction", global_tallies(LEAKAGE) % sum, & - global_tallies(LEAKAGE) % sum_sq + global_tallies(LEAKAGE) % sum_sq write(ou,*) ! format for write statements @@ -1393,11 +1398,11 @@ contains integer :: l ! loop index for nonPN scores integer :: type ! type of tally filter integer :: indent ! number of spaces to preceed output - integer :: filter_index ! index in scores array for filters + integer :: filter_index ! index in results array for filters 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 scattering 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 @@ -1459,14 +1464,14 @@ contains ! Multiply uncertainty by t-value if (confidence_intervals) then - do k = 1, size(t % scores, 2) - do j = 1, size(t % scores, 1) + do k = 1, size(t % results, 2) + do j = 1, size(t % results, 1) ! Calculate t-value for confidence intervals if (confidence_intervals) then alpha = ONE - CONFIDENCE_LEVEL t_value = t_percentile(ONE - alpha/TWO, t % n_realizations - 1) end if - t % scores(j,k) % sum_sq = t_value * t % scores(j,k) % sum_sq + t % results(j,k) % sum_sq = t_value * t % results(j,k) % sum_sq end do end do end if @@ -1474,10 +1479,10 @@ contains ! Write header block if (t % label == "") then call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, & - level=3) + level=3) else call header("TALLY " // trim(to_str(t % id)) // ": " & - // trim(t % label), unit=UNIT_TALLY, level=3) + // trim(t % label), unit=UNIT_TALLY, level=3) endif ! Handle surface current tallies separately @@ -1486,7 +1491,7 @@ contains cycle end if - ! WARNING: Admittedly, the logic for moving for printing scores is + ! WARNING: Admittedly, the logic for moving for printing results is ! extremely confusing and took quite a bit of time to get correct. The ! logic is structured this way since it is not practical to have a do ! loop for each filter variable (given that only a few filters are likely @@ -1517,7 +1522,7 @@ contains indent = indent - 2 ! ================================================================= - ! VALID BIN -- WRITE FILTER INFORMATION OR EXIT TO WRITE SCORES + ! VALID BIN -- WRITE FILTER INFORMATION OR EXIT TO WRITE RESULTS else ! Check if this is last filter @@ -1526,7 +1531,7 @@ contains ! Print current filter information type = t % filters(j) % type write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(filter_name(type)), trim(get_label(t, j)) + trim(filter_name(type)), trim(get_label(t, j)) indent = indent + 2 j = j + 1 end if @@ -1537,7 +1542,7 @@ contains if (t % n_filters > 0) then type = t % filters(j) % type write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(filter_name(type)), trim(get_label(t, j)) + trim(filter_name(type)), trim(get_label(t, j)) end if ! Determine scoring index for this bin combination -- note that unlike @@ -1549,8 +1554,8 @@ contains else filter_index = 1 end if - - ! Write scores for this filter bin combination + + ! Write results for this filter bin combination score_index = 0 if (t % n_filters > 0) indent = indent + 2 do n = 1, t % n_nuclide_bins @@ -1558,11 +1563,11 @@ contains i_nuclide = t % nuclide_bins(n) if (i_nuclide == -1) then write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - "Total Material" + "Total Material" else i_listing = nuclides(i_nuclide) % listing write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(xs_listings(i_listing) % alias) + trim(xs_listings(i_listing) % alias) end if indent = indent + 2 @@ -1639,7 +1644,7 @@ contains integer :: n ! number of incoming energy bins integer :: len1 ! length of string integer :: len2 ! length of string - integer :: filter_index ! index in scores array for filters + integer :: filter_index ! index in results array for filters logical :: print_ebin ! should incoming energy bin be displayed? character(MAX_LINE_LEN) :: string type(StructuredMesh), pointer :: m => null() @@ -1680,110 +1685,110 @@ contains ! Write incoming energy bin write(UNIT=UNIT_TALLY, FMT='(3X,A,1X,A)') & - "Incoming Energy", trim(get_label(t, i_filter_ein)) + "Incoming Energy", trim(get_label(t, i_filter_ein)) end if ! Left Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Left", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Left", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Left", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Left", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) ! Right Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Right", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Right", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Right", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Right", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) ! Back Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Back", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Back", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Back", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Back", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) ! Front Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Front", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Front", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Front", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Front", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) ! Bottom Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Bottom", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Bottom", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Bottom", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Bottom", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) ! Top Surface t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) + mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) t % matching_bins(i_filter_surf) = IN_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Incoming Current from Top", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Incoming Current from Top", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) t % matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((t % matching_bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & - "Outgoing Current to Top", & - to_str(t % scores(1,filter_index) % sum), & - trim(to_str(t % scores(1,filter_index) % sum_sq)) + "Outgoing Current to Top", & + to_str(t % results(1,filter_index) % sum), & + trim(to_str(t % results(1,filter_index) % sum_sq)) end do end do @@ -1831,10 +1836,10 @@ contains call bin_to_mesh_indices(m, bin, ijk) if (m % n_dimension == 2) then label = "Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" + trim(to_str(ijk(2))) // ")" elseif (m % n_dimension == 3) then label = "Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" end if case (FILTER_ENERGYIN, FILTER_ENERGYOUT) E0 = t % filters(i_filter) % real_bins(bin) diff --git a/src/physics.F90 b/src/physics.F90 index 58676d7152..23f83d5ecd 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -19,7 +19,7 @@ module physics use search, only: binary_search use string, only: to_str use tally, only: score_analog_tally, score_tracklength_tally, & - score_surface_current, add_to_score + score_surface_current implicit none @@ -103,12 +103,10 @@ contains if (associated(active_tracklength_tallies)) & call score_tracklength_tally(distance) - ! Score global tracklength tallies - if (tallies_on) then - ! Score track-length estimate of k-eff - call add_to_score(global_tallies(K_TRACKLENGTH), & - p % wgt * distance * material_xs % nu_fission) - end if + ! Score track-length estimate of k-eff + if (tallies_on) global_tallies(K_TRACKLENGTH) % value = & + global_tallies(K_TRACKLENGTH) % value + p % wgt * distance * & + material_xs % nu_fission if (d_collision > d_boundary) then ! ==================================================================== @@ -132,10 +130,9 @@ contains ! PARTICLE HAS COLLISION ! Score collision estimate of keff - if (tallies_on) then - call add_to_score(global_tallies(K_COLLISION), & - p % wgt * material_xs % nu_fission / material_xs % total) - end if + if (tallies_on) global_tallies(K_COLLISION) % value = & + global_tallies(K_COLLISION) % value + p % wgt * & + material_xs % nu_fission / material_xs % total p % surface = NONE call collision() @@ -908,7 +905,8 @@ contains end if ! Add to analog estimate of keff - call add_to_score(global_tallies(K_ANALOG), nu/weight * keff) + global_tallies(K_ANALOG) % value = & + global_tallies(K_ANALOG) % value + nu/weight * keff ! Bank source neutrons if (nu == 0 .or. n_bank == 3*work) return diff --git a/src/plot.F90 b/src/plot.F90 index a0ccbabe67..4e9d48c257 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -95,7 +95,7 @@ contains allocate(p) call initialize_particle() p % coord % xyz = xyz - p % coord % uvw = (/ 1, 0, 0 /) + p % coord % uvw = (/ 0.5, 0.5, 0.5 /) p % coord % universe = BASE_UNIVERSE do y = 1, img % height diff --git a/src/source_header.F90 b/src/source_header.F90 index cb0eb290c3..11c74d110c 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -4,7 +4,7 @@ module source_header !=============================================================================== ! EXTSOURCE describes an external source of neutrons for a fixed-source problem -! or for the starting source in a criticality problem +! or for the starting source in a k eigenvalue problem !=============================================================================== type ExtSource diff --git a/src/state_point.F90 b/src/state_point.F90 index 72ec8463ab..e328334768 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -82,7 +82,7 @@ contains ! If using the no-tally-reduction method, we need to collect tally ! results before writing them to the state point file. - call write_tally_scores_nr(fh) + call write_tally_results_nr(fh) elseif (master) then ! Write number of realizations @@ -93,7 +93,7 @@ contains call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) call MPI_FILE_WRITE(fh, global_tallies, N_GLOBAL_TALLIES, & - MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err) + MPI_TALLYRESULT, MPI_STATUS_IGNORE, mpi_err) if (tallies_on) then ! Indicate that tallies are on @@ -101,14 +101,14 @@ contains call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) - ! Write all tally scores - TALLY_SCORES: do i = 1, n_tallies + ! Write all tally results + TALLY_RESULTS: do i = 1, n_tallies t => tallies(i) - n = size(t % scores, 1) * size(t % scores, 2) - call MPI_FILE_WRITE(fh, t % scores, n, MPI_TALLYSCORE, & + n = size(t % results, 1) * size(t % results, 2) + call MPI_FILE_WRITE(fh, t % results, n, MPI_TALLYRESULT, & MPI_STATUS_IGNORE, mpi_err) - end do TALLY_SCORES + end do TALLY_RESULTS else ! Indicate that tallies are off temp = 0 @@ -120,7 +120,7 @@ contains ! ========================================================================== ! SOURCE BANK - if (run_mode == MODE_CRITICALITY) then + if (run_mode == MODE_EIGENVALUE) then if (source_separate) then ! If the user has specified that the source sites should be written in ! a separate file, we make a call to the appropriate subroutine to @@ -179,8 +179,8 @@ contains ! Write out current batch number write(UNIT_STATE) current_batch - ! Write out information for criticality run - if (run_mode == MODE_CRITICALITY) then + ! Write out information for eigenvalue run + if (run_mode == MODE_EIGENVALUE) then write(UNIT_STATE) n_inactive, gen_per_batch write(UNIT_STATE) k_batch(1:current_batch) write(UNIT_STATE) entropy(1:current_batch) @@ -209,7 +209,7 @@ contains ! Number of realizations write(UNIT_STATE) t % n_realizations - ! Write size of each dimension of tally scores array + ! Write size of each dimension of tally results array write(UNIT_STATE) t % total_score_bins write(UNIT_STATE) t % total_filter_bins @@ -263,25 +263,25 @@ contains ! Indicate that tallies are on write(UNIT_STATE) 1 - TALLY_SCORES: do i = 1, n_tallies + TALLY_RESULTS: do i = 1, n_tallies ! Get pointer to tally t => tallies(i) ! Write tally sum and sum_sq for each bin - do k = 1, size(t % scores, 2) - do j = 1, size(t % scores, 1) - write(UNIT_STATE) t % scores(j,k) % sum - write(UNIT_STATE) t % scores(j,k) % sum_sq + do k = 1, size(t % results, 2) + do j = 1, size(t % results, 1) + write(UNIT_STATE) t % results(j,k) % sum + write(UNIT_STATE) t % results(j,k) % sum_sq end do end do - end do TALLY_SCORES + end do TALLY_RESULTS else ! Indicate that tallies are off write(UNIT_STATE) 0 end if ! Write out source bank - if (run_mode == MODE_CRITICALITY) then + if (run_mode == MODE_EIGENVALUE) then if (source_separate) then ! If the user has specified that the source sites should be written in ! a separate file, we make a call to the appropriate subroutine to @@ -349,8 +349,8 @@ contains call MPI_FILE_WRITE(fh, current_batch, 1, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) - ! Write out information for criticality run - if (run_mode == MODE_CRITICALITY) then + ! Write out information for eigenvalue run + if (run_mode == MODE_EIGENVALUE) then call MPI_FILE_WRITE(fh, n_inactive, 1, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) call MPI_FILE_WRITE(fh, gen_per_batch, 1, MPI_INTEGER, & @@ -454,10 +454,10 @@ contains #ifdef MPI !=============================================================================== -! WRITE_TALLY_SCORES_NR +! WRITE_TALLY_RESULTS_NR !=============================================================================== - subroutine write_tally_scores_nr(fh) + subroutine write_tally_results_nr(fh) integer, intent(in) :: fh ! file handle @@ -466,7 +466,7 @@ contains integer :: m ! number of score bins integer :: temp ! temporary variable integer :: n_bins ! total number of bins - real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of scores + real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of results real(8) :: global_temp(2,N_GLOBAL_TALLIES) real(8) :: dummy ! temporary receive buffer for non-root reduces type(TallyObject), pointer :: t => null() @@ -518,20 +518,20 @@ contains MPI_STATUS_IGNORE, mpi_err) end if - ! Write all tally scores - TALLY_SCORES: do i = 1, n_tallies + ! Write all tally results + TALLY_RESULTS: do i = 1, n_tallies t => tallies(i) - ! Determine size of tally scores array - m = size(t % scores, 1) - n = size(t % scores, 2) + ! Determine size of tally results array + m = size(t % results, 1) + n = size(t % results, 2) n_bins = m*n*2 ! Allocate array for storing sums and sums of squares, but ! contiguously in memory for each allocate(tally_temp(2,m,n)) - tally_temp(1,:,:) = t % scores(:,:) % sum - tally_temp(2,:,:) = t % scores(:,:) % sum_sq + tally_temp(1,:,:) = t % results(:,:) % sum + tally_temp(2,:,:) = t % results(:,:) % sum_sq if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into @@ -544,10 +544,10 @@ contains MPI_STATUS_IGNORE, mpi_err) ! At the end of the simulation, store the results back in the - ! regular TallyScores array + ! regular TallyResults array if (current_batch == n_batches) then - t % scores(:,:) % sum = tally_temp(1,:,:) - t % scores(:,:) % sum_sq = tally_temp(2,:,:) + t % results(:,:) % sum = tally_temp(1,:,:) + t % results(:,:) % sum_sq = tally_temp(2,:,:) end if else ! Receive buffer not significant at other processors @@ -557,7 +557,7 @@ contains ! Deallocate temporary copy of tally results deallocate(tally_temp) - end do TALLY_SCORES + end do TALLY_RESULTS else if (master) then ! Indicate that tallies are off @@ -567,7 +567,7 @@ contains end if end if - end subroutine write_tally_scores_nr + end subroutine write_tally_results_nr #endif !=============================================================================== @@ -644,8 +644,8 @@ contains call MPI_FILE_READ_ALL(fh, restart_batch, 1, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) - ! Read information specific to criticality run - if (mode == MODE_CRITICALITY) then + ! Read information specific to eigenvalue run + if (mode == MODE_EIGENVALUE) then call MPI_FILE_READ_ALL(fh, n_inactive, 1, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) call MPI_FILE_READ_ALL(fh, gen_per_batch, 1, MPI_INTEGER, & @@ -687,12 +687,12 @@ contains call MPI_FILE_READ(fh, tallies(i) % n_realizations, 1, & MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) - ! Read dimensions of tally filters and scores and make sure they + ! Read dimensions of tally filters and results and make sure they ! match call MPI_FILE_READ(fh, temp, 2, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) - if (temp(1) /= size(tallies(i) % scores, 1) .or. & - temp(2) /= size(tallies(i) % scores, 2)) then + if (temp(1) /= size(tallies(i) % results, 1) .or. & + temp(2) /= size(tallies(i) % results, 2)) then message = "Tally dimensions do not match in state point." call fatal_error() end if @@ -736,7 +736,7 @@ contains MPI_STATUS_IGNORE, mpi_err) deallocate(int_array) - ! Read number of scores + ! Read number of results call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, & MPI_STATUS_IGNORE, mpi_err) @@ -760,7 +760,7 @@ contains ! Read global tally data call MPI_FILE_READ(fh, global_tallies, N_GLOBAL_TALLIES, & - MPI_TALLYSCORE, MPI_STATUS_IGNORE, mpi_err) + MPI_TALLYRESULT, MPI_STATUS_IGNORE, mpi_err) ! Check if tally results are present call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err) @@ -770,18 +770,18 @@ contains ! Read sum and sum squared if (temp(1) == 1) then - TALLY_SCORES: do i = 1, n_tallies - n = size(tallies(i) % scores, 1) * size(tallies(i) % scores, 2) - call MPI_FILE_READ(fh, tallies(i) % scores, n, MPI_TALLYSCORE, & + TALLY_RESULTS: do i = 1, n_tallies + n = size(tallies(i) % results, 1) * size(tallies(i) % results, 2) + call MPI_FILE_READ(fh, tallies(i) % results, n, MPI_TALLYRESULT, & MPI_STATUS_IGNORE, mpi_err) - end do TALLY_SCORES + end do TALLY_RESULTS end if end if ! ========================================================================== ! SOURCE BANK - if (run_mode == MODE_CRITICALITY) then + if (run_mode == MODE_EIGENVALUE) then ! Get current offset for master if (master) call MPI_FILE_GET_POSITION(fh, offset, mpi_err) @@ -841,8 +841,8 @@ contains ! Read batch number to restart at read(UNIT_STATE) restart_batch - ! Read information specific to criticality run - if (mode == MODE_CRITICALITY) then + ! Read information specific to eigenvalue run + if (mode == MODE_EIGENVALUE) then read(UNIT_STATE) n_inactive, gen_per_batch read(UNIT_STATE) k_batch(1:restart_batch) read(UNIT_STATE) entropy(1:restart_batch) @@ -886,11 +886,11 @@ contains ! Read number of realizations read(UNIT_STATE) tallies(i) % n_realizations - ! Read dimensions of tally filters and scores and make sure they + ! Read dimensions of tally filters and results and make sure they ! match read(UNIT_STATE) temp(1:2) - if (temp(1) /= size(tallies(i) % scores, 1) .or. & - temp(2) /= size(tallies(i) % scores, 2)) then + if (temp(1) /= size(tallies(i) % results, 1) .or. & + temp(2) /= size(tallies(i) % results, 2)) then message = "Tally dimensions do not match in state point." call fatal_error() end if @@ -927,7 +927,7 @@ contains read(UNIT_STATE) int_array deallocate(int_array) - ! Read number of scores + ! Read number of results read(UNIT_STATE) temp(1) ! Read nuclide bins @@ -955,19 +955,19 @@ contains ! Read sum and sum squared read(UNIT_STATE) temp(1) if (temp(1) == 1) then - TALLY_SCORES: do i = 1, n_tallies - do k = 1, size(tallies(i) % scores, 2) - do j = 1, size(tallies(i) % scores, 1) - read(UNIT_STATE) tallies(i) % scores(j,k) % sum - read(UNIT_STATE) tallies(i) % scores(j,k) % sum_sq + TALLY_RESULTS: do i = 1, n_tallies + do k = 1, size(tallies(i) % results, 2) + do j = 1, size(tallies(i) % results, 1) + read(UNIT_STATE) tallies(i) % results(j,k) % sum + read(UNIT_STATE) tallies(i) % results(j,k) % sum_sq end do end do - end do TALLY_SCORES + end do TALLY_RESULTS end if end if - ! Read source bank for criticality run - if (mode == MODE_CRITICALITY .and. run_mode /= MODE_TALLIES) then + ! Read source bank for eigenvalue run + if (mode == MODE_EIGENVALUE .and. run_mode /= MODE_TALLIES) then read(UNIT_STATE) source_bank end if @@ -995,7 +995,7 @@ contains call write_message(1) end if - ! For criticality calculations, turn on tallies if we've reached active + ! For eigenvalue calculations, turn on tallies if we've reached active ! batches if (current_batch == n_inactive) then tallies_on = .true. diff --git a/src/tally.F90 b/src/tally.F90 index 2f81d8cfee..ba20f6aa53 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -6,12 +6,12 @@ module tally use global use math, only: t_percentile, calc_pn use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices, & - mesh_indices_to_bin, mesh_intersects + mesh_indices_to_bin, mesh_intersects use mesh_header, only: StructuredMesh use output, only: header use search, only: binary_search use string, only: to_str - use tally_header, only: TallyScore, TallyMapItem, TallyMapElement + use tally_header, only: TallyResult, TallyMapItem, TallyMapElement #ifdef MPI use mpi @@ -76,10 +76,10 @@ contains end if ! ======================================================================= - ! CALCULATE SCORES AND ACCUMULATE TALLY + ! CALCULATE RESULTS AND ACCUMULATE TALLY ! If we have made it here, we have a scoring combination of bins for this - ! tally -- now we need to determine where in the scores array we should + ! tally -- now we need to determine where in the results array we should ! be accumulating the tally values ! Determine scoring index for this filter combination @@ -197,7 +197,7 @@ contains ! the end of this select case will handle that one. if (n < t % scatt_order(j)) then t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) + score + t % results(score_index, filter_index) % value + score end if end do j = j + t % scatt_order(j) @@ -229,7 +229,7 @@ contains ! mu*Sigma_s)). score = last_wgt / (3.0_8 * score * (material_xs % total - & - mu * score)) + mu * score)) case (SCORE_N_1N) ! Skip any event where the particle didn't scatter @@ -299,7 +299,7 @@ contains ! bin. However, in the case of fission, since multiple ! fission neutrons were emitted with different energies, ! multiple outgoing energy bins may have been scored to. The - ! following logic treats this special case and scores to + ! following logic treats this special case and results to ! multiple bins call score_fission_eout(t, score_index) @@ -322,13 +322,12 @@ contains case default message = "Invalid score type on tally " // to_str(t % id) // "." - write(*,*) score_index, score_bin, l,j, t % n_score_bins call fatal_error() end select ! Add score to tally t % results(score_index, filter_index) % value = & - t % results(score_index, filter_index) + score + t % results(score_index, filter_index) % value + score end do SCORE_LOOP @@ -391,13 +390,14 @@ contains ! change outgoing energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - size(t % filters(i) % real_bins), E_out) + size(t % filters(i) % real_bins), E_out) ! determine scoring index i_filter = sum((t % matching_bins - 1) * t % stride) + 1 ! Add score to tally - call add_to_score(t % scores(i_score, i_filter), score) + t % results(i_score, i_filter) % value = & + t % results(i_score, i_filter) % value + score end do ! reset outgoing energy bin and score index @@ -461,10 +461,10 @@ contains end if ! ======================================================================= - ! CALCULATE SCORES AND ACCUMULATE TALLY + ! CALCULATE RESULTS AND ACCUMULATE TALLY ! If we have made it here, we have a scoring combination of bins for this - ! tally -- now we need to determine where in the scores array we should + ! tally -- now we need to determine where in the results array we should ! be accumulating the tally values ! Determine scoring index for this filter combination @@ -472,7 +472,7 @@ contains if (t % all_nuclides) then call score_all_nuclides(tracklength_tallies(curr_ptr % data), flux, & - filter_index) + filter_index) else NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins @@ -509,20 +509,20 @@ contains select case(score_bin) case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * & - atom_density * flux + atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux + micro_xs(i_nuclide) % absorption) * & + atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * & - atom_density * flux + atom_density * flux case (SCORE_FISSION) score = micro_xs(i_nuclide) % fission * & - atom_density * flux + atom_density * flux case (SCORE_NU_FISSION) score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux + atom_density * flux case (SCORE_EVENTS) score = ONE case default @@ -557,7 +557,8 @@ contains score_index = (k - 1)*t % n_score_bins + j ! Add score to tally - call add_to_score(t % scores(score_index, filter_index), score) + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score end do SCORE_LOOP @@ -632,7 +633,7 @@ contains score = micro_xs(i_nuclide) % total * atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * atom_density * flux + micro_xs(i_nuclide) % absorption) * atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * atom_density * flux case (SCORE_FISSION) @@ -651,7 +652,8 @@ contains score_index = (i_nuclide - 1)*t % n_score_bins + j ! Add score to tally - call add_to_score(t % scores(score_index, filter_index), score) + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score end do SCORE_LOOP @@ -691,7 +693,8 @@ contains score_index = n_nuclides_total*t % n_score_bins + j ! Add score to tally - call add_to_score(t % scores(score_index, filter_index), score) + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score end do MATERIAL_SCORE_LOOP @@ -775,27 +778,27 @@ contains ! determine next universe bin ! TODO: Account for multiple universes when performing this filter t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord % universe, i_tally) case (FILTER_MATERIAL) t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & - p % material, i_tally) + p % material, i_tally) case (FILTER_CELL) ! determine next cell bin ! TODO: Account for cells in multiple levels when performing this filter t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + p % coord % cell, i_tally) case (FILTER_CELLBORN) ! determine next cellborn bin t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & - p % cell_born, i_tally) + p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & - p % surface, i_tally) + p % surface, i_tally) case (FILTER_ENERGYIN) ! determine incoming energy bin @@ -803,12 +806,12 @@ contains ! check if energy of the particle is within energy bins if (p % E < t % filters(i) % real_bins(1) .or. & - p % E > t % filters(i) % real_bins(k + 1)) then + p % E > t % filters(i) % real_bins(k + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - k + 1, p % E) + k + 1, p % E) end if end select @@ -929,25 +932,25 @@ contains select case(score_bin) case (SCORE_TOTAL) score = micro_xs(i_nuclide) % total * & - atom_density * flux + atom_density * flux case (SCORE_SCATTER) score = (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) * & - atom_density * flux + micro_xs(i_nuclide) % absorption) * & + atom_density * flux case (SCORE_ABSORPTION) score = micro_xs(i_nuclide) % absorption * & - atom_density * flux + atom_density * flux case (SCORE_FISSION) score = micro_xs(i_nuclide) % fission * & - atom_density * flux + atom_density * flux case (SCORE_NU_FISSION) score = micro_xs(i_nuclide) % nu_fission * & - atom_density * flux + atom_density * flux case (SCORE_EVENTS) score = ONE case default message = "Invalid score type on tally " // & - to_str(t % id) // "." + to_str(t % id) // "." call fatal_error() end select @@ -970,7 +973,7 @@ contains score = ONE case default message = "Invalid score type on tally " // & - to_str(t % id) // "." + to_str(t % id) // "." call fatal_error() end select end if @@ -979,7 +982,8 @@ contains score_index = (b - 1)*t % n_score_bins + j ! Add score to tally - call add_to_score(t % scores(score_index, filter_index), score) + t % results(score_index, filter_index) % value = & + t % results(score_index, filter_index) % value + score end do SCORE_LOOP @@ -1028,27 +1032,27 @@ contains ! determine next universe bin ! TODO: Account for multiple universes when performing this filter t % matching_bins(i) = get_next_bin(FILTER_UNIVERSE, & - p % coord % universe, i_tally) + p % coord % universe, i_tally) case (FILTER_MATERIAL) t % matching_bins(i) = get_next_bin(FILTER_MATERIAL, & - p % material, i_tally) + p % material, i_tally) case (FILTER_CELL) ! determine next cell bin ! TODO: Account for cells in multiple levels when performing this filter t % matching_bins(i) = get_next_bin(FILTER_CELL, & - p % coord % cell, i_tally) + p % coord % cell, i_tally) case (FILTER_CELLBORN) ! determine next cellborn bin t % matching_bins(i) = get_next_bin(FILTER_CELLBORN, & - p % cell_born, i_tally) + p % cell_born, i_tally) case (FILTER_SURFACE) ! determine next surface bin t % matching_bins(i) = get_next_bin(FILTER_SURFACE, & - p % surface, i_tally) + p % surface, i_tally) case (FILTER_ENERGYIN) ! determine incoming energy bin @@ -1063,12 +1067,12 @@ contains ! check if energy of the particle is within energy bins if (E < t % filters(i) % real_bins(1) .or. & - E > t % filters(i) % real_bins(n + 1)) then + E > t % filters(i) % real_bins(n + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - n + 1, E) + n + 1, E) end if case (FILTER_ENERGYOUT) @@ -1077,12 +1081,12 @@ contains ! check if energy of the particle is within energy bins if (E < t % filters(i) % real_bins(1) .or. & - E > t % filters(i) % real_bins(n + 1)) then + E > t % filters(i) % real_bins(n + 1)) then t % matching_bins(i) = NO_BIN_FOUND else ! search to find incoming energy bin t % matching_bins(i) = binary_search(t % filters(i) % real_bins, & - n + 1, p % E) + n + 1, p % E) end if end select @@ -1173,14 +1177,14 @@ contains n = t % filters(j) % n_bins ! check if energy of the particle is within energy bins if (p % E < t % filters(j) % real_bins(1) .or. & - p % E > t % filters(j) % real_bins(n + 1)) then + p % E > t % filters(j) % real_bins(n + 1)) then curr_ptr => curr_ptr % next ! select next tally cycle end if ! search to find incoming energy bin t % matching_bins(j) = binary_search(t % filters(j) % real_bins, & - n + 1, p % E) + n + 1, p % E) end if ! ======================================================================= @@ -1198,9 +1202,10 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 - call add_to_score(t % scores(1, filter_index), p % wgt) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if end do else @@ -1209,9 +1214,10 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 - call add_to_score(t % scores(1, filter_index), p % wgt) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if end do end if @@ -1225,9 +1231,10 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 - call add_to_score(t % scores(1, filter_index), p % wgt) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if end do else @@ -1236,9 +1243,10 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 - call add_to_score(t % scores(1, filter_index), p % wgt) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if end do end if @@ -1252,9 +1260,10 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 - call add_to_score(t % scores(1, filter_index), p % wgt) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if end do else @@ -1263,9 +1272,10 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) filter_index = sum((t % matching_bins - 1) * t % stride) + 1 - call add_to_score(t % scores(1, filter_index), p % wgt) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if end do end if @@ -1316,7 +1326,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(1) = ijk0(1) + 1 xyz_cross(1) = xyz_cross(1) + m % width(1) @@ -1328,7 +1338,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_RIGHT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if elseif (distance == d(2)) then @@ -1338,7 +1348,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(2) = ijk0(2) + 1 xyz_cross(2) = xyz_cross(2) + m % width(2) @@ -1350,7 +1360,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_FRONT t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if else if (distance == d(3)) then @@ -1360,7 +1370,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = OUT_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if ijk0(3) = ijk0(3) + 1 xyz_cross(3) = xyz_cross(3) + m % width(3) @@ -1372,7 +1382,7 @@ contains if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then t % matching_bins(i_filter_surf) = IN_TOP t % matching_bins(i_filter_mesh) = & - mesh_indices_to_bin(m, ijk0 + 1, .true.) + mesh_indices_to_bin(m, ijk0 + 1, .true.) end if end if end if @@ -1383,13 +1393,14 @@ contains ! Check for errors if (filter_index <= 0 .or. filter_index > & - t % total_filter_bins) then + t % total_filter_bins) then message = "Score index outside range." call fatal_error() end if ! Add to surface current tally - call add_to_score(t % scores(1, filter_index), p % wgt) + t % results(1, filter_index) % value = & + t % results(1, filter_index) % value + p % wgt end if ! Calculate new coordinates @@ -1440,7 +1451,7 @@ contains end if i_tally_check = tally_maps(filter_type) % items(filter_value) % & - elements(position(filter_type)) % index_tally + elements(position(filter_type)) % index_tally if (i_tally_check > i_tally) then ! Since the index being checked against is greater than the index we @@ -1452,7 +1463,7 @@ contains elseif (i_tally_check == i_tally) then ! Found a match bin = tally_maps(filter_type) % items(filter_value) % & - elements(position(filter_type)) % index_bin + elements(position(filter_type)) % index_bin return end if @@ -1484,21 +1495,23 @@ contains ! Accumulate on master only unless run is not reduced then do it on all if (master .or. (.not. reduce_tallies)) then - ! Accumulate scores for each tally + ! Accumulate results for each tally curr_ptr => active_tallies do while(associated(curr_ptr)) call accumulate_tally(tallies(curr_ptr % data)) curr_ptr => curr_ptr % next end do - if (run_mode == MODE_CRITICALITY) then - ! Before accumulating scores for global_tallies, we need to get the - ! current batch estimate of k_analog for displaying to output - if (active_batches) k_batch(current_batch) = global_tallies(K_ANALOG) % value + if (run_mode == MODE_EIGENVALUE) then + ! Get the current batch estimate of k_analog for displaying to output + ! --- this has to be performed after reduce_tally_values and before + ! accumulate_result + + k_batch(current_batch) = global_tallies(K_ANALOG) % value end if - ! Accumulate scores for global tallies - if (active_batches) call accumulate_score(global_tallies) + ! Accumulate results for global tallies + if (active_batches) call accumulate_result(global_tallies) end if if (associated(curr_ptr)) nullify(curr_ptr) @@ -1515,7 +1528,7 @@ contains integer :: n ! number of filter bins integer :: m ! number of score bins integer :: n_bins ! total number of bins - real(8), allocatable :: tally_temp(:,:) ! contiguous array of scores + real(8), allocatable :: tally_temp(:,:) ! contiguous array of results real(8) :: global_temp(N_GLOBAL_TALLIES) real(8) :: dummy ! temporary receive buffer for non-root reduces type(TallyObject), pointer :: t => null() @@ -1531,23 +1544,23 @@ contains allocate(tally_temp(m,n)) - tally_temp = t % scores(:,:) % value + tally_temp = t % results(:,:) % value if (master) then ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values to value on master - t % scores(:,:) % value = tally_temp + t % results(:,:) % value = tally_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors - t % scores(:,:) % value = 0 + t % results(:,:) % value = 0 end if deallocate(tally_temp) @@ -1560,14 +1573,14 @@ contains if (master) then call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values back to global_tallies on master global_tallies(:) % value = global_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors global_tallies(:) % value = ZERO @@ -1578,11 +1591,11 @@ contains ! last realization if (master) then call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, MPI_COMM_WORLD, mpi_err) else ! Receive buffer not significant at other processors call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, MPI_COMM_WORLD, mpi_err) end if if (associated(curr_ptr)) nullify(curr_ptr) @@ -1605,14 +1618,14 @@ contains t % n_realizations = t % n_realizations + n_procs end if - ! Accumulate each TallyScore - call accumulate_score(t % scores) + ! Accumulate each TallyResult + call accumulate_result(t % results) end subroutine accumulate_tally !=============================================================================== ! TALLY_STATISTICS computes the mean and standard deviation of the mean of each -! tally and stores them in the val and val_sq attributes of the TallyScores +! tally and stores them in the val and val_sq attributes of the TallyResults ! respectively !=============================================================================== @@ -1625,87 +1638,72 @@ contains do i = 1, n_tallies t => tallies(i) - call statistics_score(t % scores, t % n_realizations) + call statistics_result(t % results, t % n_realizations) end do ! Calculate statistics for global tallies - call statistics_score(global_tallies, n_realizations) + call statistics_result(global_tallies, n_realizations) end subroutine tally_statistics !=============================================================================== -! ADD_TO_SCORE accumulates a scoring contribution to a specific tally bin and -! specific response function. Note that we don't need to add the square of the -! contribution since that is done at the batch level, not the history level -!=============================================================================== - - subroutine add_to_score(score, val) - - type(TallyScore), intent(inout) :: score - real(8), intent(in) :: val - - score % value = score % value + val - - end subroutine add_to_score - -!=============================================================================== -! ACCUMULATE_SCORE accumulates scores from many histories (or many generations) +! ACCUMULATE_RESULT accumulates results from many histories (or many generations) ! into a single realization of a random variable. !=============================================================================== - elemental subroutine accumulate_score(score) + elemental subroutine accumulate_result(this) - type(TallyScore), intent(inout) :: score + type(TallyResult), intent(inout) :: this real(8) :: val - ! Add the sum and square of the sum of contributions from a tally score to + ! Add the sum and square of the sum of contributions from a tally result to ! the variables sum and sum_sq. This will later allow us to calculate a ! variance on the tallies. - val = score % value/total_weight - score % sum = score % sum + val - score % sum_sq = score % sum_sq + val*val + val = this % value/total_weight + this % sum = this % sum + val + this % sum_sq = this % sum_sq + val*val ! Reset the single batch estimate - score % value = ZERO + this % value = ZERO - end subroutine accumulate_score + end subroutine accumulate_result !=============================================================================== -! STATISTICS_SCORE determines the sample mean and the standard deviation of the -! mean for a TallyScore. +! STATISTICS_RESULT determines the sample mean and the standard deviation of the +! mean for a TallyResult. !=============================================================================== - elemental subroutine statistics_score(score, n) + elemental subroutine statistics_result(this, n) - type(TallyScore), intent(inout) :: score - integer, intent(in) :: n + type(TallyResult), intent(inout) :: this + integer, intent(in) :: n ! Calculate sample mean and standard deviation of the mean -- note that we ! have used Bessel's correction so that the estimator of the variance of the ! sample mean is unbiased. - score % sum = score % sum/n - score % sum_sq = sqrt((score % sum_sq/n - score % sum * & - score % sum) / (n - 1)) + this % sum = this % sum/n + this % sum_sq = sqrt((this % sum_sq/n - this % sum * & + this % sum) / (n - 1)) - end subroutine statistics_score + end subroutine statistics_result !=============================================================================== -! RESET_SCORE zeroes out the value and accumulated sum and sum-squared for a -! single TallyScore. +! RESET_RESULT zeroes out the value and accumulated sum and sum-squared for a +! single TallyResult. !=============================================================================== - elemental subroutine reset_score(score) + elemental subroutine reset_result(this) - type(TallyScore), intent(inout) :: score + type(TallyResult), intent(inout) :: this - score % value = ZERO - score % sum = ZERO - score % sum_sq = ZERO + this % value = ZERO + this % sum = ZERO + this % sum_sq = ZERO - end subroutine reset_score + end subroutine reset_result !=============================================================================== ! SETUP_ACTIVE_USERTALLIES @@ -1909,7 +1907,7 @@ contains end if do i = n_cmfd_current_tallies + n_user_current_tallies, & - n_user_current_tallies + 1, -1 + n_user_current_tallies + 1, -1 ! allocate node allocate(curr_ptr) diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 47d992f4dd..a9b096f718 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -35,14 +35,14 @@ module tally_header end type TallyMap !=============================================================================== -! TALLYSCORE provides accumulation of scores in a particular tally bin +! TALLYRESULT provides accumulation of results in a particular tally bin !=============================================================================== - type TallyScore + type TallyResult real(8) :: value = 0. real(8) :: sum = 0. real(8) :: sum_sq = 0. - end type TallyScore + end type TallyResult !=============================================================================== ! TALLYFILTER describes a filter that limits what events score to a tally. For @@ -59,8 +59,8 @@ module tally_header !=============================================================================== ! TALLYOBJECT describes a user-specified tally. The region of phase space to -! tally in is given by the TallyBins and the scores are stored in a TallyScore -! array. +! tally in is given by the TallyFilters and the results are stored in a +! TallyResult array. !=============================================================================== type TallyObject @@ -77,10 +77,10 @@ module tally_header integer :: n_filters ! Number of filters type(TallyFilter), allocatable :: filters(:) ! Filter data (type/bins) - ! The stride attribute is used for determining the index in the scores + ! The stride attribute is used for determining the index in the results ! array for a matching_bin combination. Since multiple dimensions are - ! mapped onto one dimension in the scores array, the stride attribute gives - ! the stride for a given filter type within the scores array + ! mapped onto one dimension in the results array, the stride attribute gives + ! the stride for a given filter type within the results array integer, allocatable :: matching_bins(:) integer, allocatable :: stride(:) @@ -105,14 +105,14 @@ module tally_header integer, allocatable :: scatt_order(:) integer :: n_nonPN_score_bins = 0 - ! Scores for each bin -- the first dimenion of the array is for scores + ! Results for each bin -- the first dimension of the array is for scores ! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the ! second dimension of the array is for the combination of filters ! (e.g. specific cell, specific energy group, etc.) integer :: total_filter_bins integer :: total_score_bins - type(TallyScore), allocatable :: scores(:,:) + type(TallyResult), allocatable :: results(:,:) ! reset property - allows a tally to be reset after every batch logical :: reset = .false. diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 916bcfde25..d8e503a85c 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -25,7 +25,7 @@ contains !=============================================================================== ! SETUP_TALLY_ARRAYS allocates and populates several member arrays of the -! TallyObject derived type, including stride, matching_bins, and scores. +! TallyObject derived type, including stride, matching_bins, and results. !=============================================================================== subroutine setup_tally_arrays() @@ -57,8 +57,8 @@ contains t % total_filter_bins = n t % total_score_bins = t % n_score_bins * t % n_nuclide_bins - ! Allocate scores array - allocate(t % scores(t % total_score_bins, t % total_filter_bins)) + ! Allocate results array + allocate(t % results(t % total_score_bins, t % total_filter_bins)) end do TALLY_LOOP @@ -68,7 +68,7 @@ contains ! SETUP_TALLY_MAPS creates a map that allows a quick determination of which ! tallies and bins need to be scored to when a particle makes a collision. This ! subroutine also sets the stride attribute for each tally as well as allocating -! storage for the scores array. +! storage for the results array. !=============================================================================== subroutine setup_tally_maps() diff --git a/src/templates/settings.rnc b/src/templates/settings.rnc index 772a182e96..a2d75254fd 100644 --- a/src/templates/settings.rnc +++ b/src/templates/settings.rnc @@ -2,7 +2,7 @@ element settings { element confidence_intervals { ( "on" | "off" ) }? & ( - element criticality { + element eigenvalue { (element batches { xsd:positiveInteger } | attribute batches { xsd:positiveInteger }) & (element inactive { xsd:nonNegativeInteger } | diff --git a/src/templates/settings_t.xml b/src/templates/settings_t.xml index 2a220b7269..c897c99f83 100644 --- a/src/templates/settings_t.xml +++ b/src/templates/settings_t.xml @@ -42,9 +42,9 @@ - + @@ -60,4 +60,7 @@ + + + diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 3e1cdab3a4..0f400f8fea 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -12,9 +12,10 @@ filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', score_types = {-1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter', -5: 'scatter-1', -6: 'scatter-2', -7: 'scatter-3', - -8: 'transport', -9: 'diffusion', -10: 'n1n', -11: 'n2n', - -12: 'n3n', -13: 'n4n', -14: 'absorption', -15: 'fission', - -16: 'nu-fission', -17: 'current', -18: 'events'} + -8: 'scatter-4', -9: 'scatter-5', -10: 'transport', + -11: 'diffusion', -12: 'n1n', -13: 'n2n', -14: 'n3n', + -15: 'n4n', -16: 'absorption', -17: 'fission', + -18: 'nu-fission', -19: 'current', -20: 'events'} class BinaryFile(object): def __init__(self, filename): @@ -80,7 +81,7 @@ class StatePoint(BinaryFile): # Set flags for what data was read self._metadata = False - self._values = False + self._results = False self._source = False # Initialize arrays for meshes and tallies @@ -193,7 +194,7 @@ class StatePoint(BinaryFile): # Set flag indicating metadata has already been read self._metadata = True - def read_values(self): + def read_results(self): # Check whether metadata has been read if not self._metadata: self._read_metadata() @@ -214,16 +215,16 @@ class StatePoint(BinaryFile): for t in self.tallies: n = t.total_score_bins * t.total_filter_bins - t.values = np.array(self._get_double(2*n)) - t.values.shape = (t.total_filter_bins, t.total_score_bins, 2) + t.results = np.array(self._get_double(2*n)) + t.results.shape = (t.total_filter_bins, t.total_score_bins, 2) - # Indicate that tally values have been read - self._values = True + # Indicate that tally results have been read + self._results = True def read_source(self): - # Check whether tally values have been read - if not self._values: - self.read_values() + # Check whether tally results have been read + if not self._results: + self.read_results() for i in range(self.n_particles): s = SourceSite() @@ -273,18 +274,18 @@ class StatePoint(BinaryFile): # Regular tallies for t in self.tallies: - for i in range(t.values.shape[0]): - for j in range(t.values.shape[1]): + for i in range(t.results.shape[0]): + for j in range(t.results.shape[1]): # Get sum and sum of squares - s, s2 = t.values[i,j] + s, s2 = t.results[i,j] # Calculate sample mean and replace value s /= n - t.values[i,j,0] = s + t.results[i,j,0] = s # Calculate standard deviation if s != 0.0: - t.values[i,j,1] = t_value*sqrt((s2/n - s*s)/(n-1)) + t.results[i,j,1] = t_value*sqrt((s2/n - s*s)/(n-1)) def get_value(self, tally_index, spec_list, score_index): """Returns a tally score given a list of filters to satisfy. @@ -304,14 +305,14 @@ class StatePoint(BinaryFile): score_index : int Index corresponding to score for tally, i.e. the second index in - Tally.values[:,:,:]. + Tally.results[:,:,:]. """ # Get Tally object given the index t = self.tallies[tally_index] - # Initialize index for filter in Tally.values[:,:,:] + # Initialize index for filter in Tally.results[:,:,:] filter_index = 0 # Loop over specified filters in spec_list @@ -334,7 +335,7 @@ class StatePoint(BinaryFile): else: filter_index += f_index*t.filters[f_type].stride - # Return the desired result from Tally.values. This could be the sum and + # Return the desired result from Tally.results. This could be the sum and # sum of squares, or it could be mean and stdev if self.generate_stdev() # has been called already. - return t.values[filter_index, score_index] + return t.results[filter_index, score_index] diff --git a/src/utils/statepoint_cmp.py b/src/utils/statepoint_cmp.py index 52ac000dab..c4ef18a810 100755 --- a/src/utils/statepoint_cmp.py +++ b/src/utils/statepoint_cmp.py @@ -17,8 +17,8 @@ sp1 = StatePoint(path1) sp2 = StatePoint(path2) # Read tally results -sp1.read_values() -sp2.read_values() +sp1.read_results() +sp2.read_results() # Compare header information assert sp1.revision == sp2.revision @@ -71,7 +71,7 @@ for t1, t2 in zip(sp1.tallies, sp2.tallies): assert t1.scores == t2.scores # Compare tally results - assert_allclose(t1.values, t2.values) + assert_allclose(t1.results, t2.results) # If criticality, compare source sites if sp1.run_mode == 2: