From b03786d5363c6b8442b56e8d067d1f9440225bef Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 14 Aug 2013 23:39:07 -0500 Subject: [PATCH 01/31] Simple implementation of random number streams. Setup separate streams for tracking and tallies. --- src/random_lcg.F90 | 43 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 37 insertions(+), 6 deletions(-) diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 3130d42a53..12b85c2c7a 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -5,8 +5,13 @@ module random_lcg private save + ! Random number streams + integer, parameter :: N_STREAMS = 2 + integer, parameter :: STREAM_TRACKING = 1 + integer, parameter :: STREAM_TALLIES = 2 + integer(8) :: prn_seed0 ! original seed - integer(8) :: prn_seed ! current seed + integer(8) :: prn_seed(N_STREAMS) ! current seed integer(8) :: prn_mult ! multiplication factor, g integer(8) :: prn_add ! additive factor, c integer :: prn_bits ! number of bits, M @@ -14,11 +19,14 @@ module random_lcg integer(8) :: prn_mask ! 2^M - 1 integer(8) :: prn_stride ! stride between particles real(8) :: prn_norm ! 2^(-M) + integer :: stream ! current RNG stream public :: prn public :: initialize_prng public :: set_particle_seed public :: prn_skip + public :: prn_set_stream + public :: STREAM_TRACKING, STREAM_TALLIES contains @@ -33,12 +41,12 @@ contains ! This algorithm uses bit-masking to find the next integer(8) value to be ! used to calculate the random number - prn_seed = iand(prn_mult*prn_seed + prn_add, prn_mask) + prn_seed(stream) = iand(prn_mult*prn_seed(stream) + prn_add, prn_mask) ! Once the integer is calculated, we just need to divide by 2**m, ! represented here as multiplying by a pre-calculated factor - pseudo_rn = prn_seed * prn_norm + pseudo_rn = prn_seed(stream) * prn_norm end function prn @@ -51,8 +59,13 @@ contains use global, only: seed + integer :: i + + stream = STREAM_TRACKING prn_seed0 = seed - prn_seed = prn_seed + do i = 1, N_STREAMS + prn_seed(i) = prn_seed0 + i - 1 + end do prn_mult = 2806196910506780709_8 prn_add = 1_8 prn_bits = 63 @@ -72,7 +85,11 @@ contains integer(8), intent(in) :: id - prn_seed = prn_skip_ahead(id*prn_stride, prn_seed0) + integer :: i + + do i = 1, N_STREAMS + prn_seed(i) = prn_skip_ahead(id*prn_stride, prn_seed0 + i - 1) + end do end subroutine set_particle_seed @@ -84,7 +101,7 @@ contains integer(8), intent(in) :: n ! number of seeds to skip - prn_seed = prn_skip_ahead(n, prn_seed) + prn_seed(stream) = prn_skip_ahead(n, prn_seed(stream)) end subroutine prn_skip @@ -149,4 +166,18 @@ contains end function prn_skip_ahead +!=============================================================================== +! PRN_SET_STREAM changes the random number stream. If random numbers are needed +! in routines not used directly for tracking (e.g. physics), this allows the +! numbers to be generated without affecting reproducibility of the physics. +!=============================================================================== + + subroutine prn_set_stream(i) + + integer, intent(in) :: i + + stream = i + + end subroutine prn_set_stream + end module random_lcg From acb16a3e884b387bd2bf1eb75e25876a334ab9fe Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 29 Nov 2013 13:47:46 -0500 Subject: [PATCH 02/31] Prevent unused variable warning in cmfd_input. --- src/cmfd_input.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index ac3d59824e..120db5f649 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -20,7 +20,9 @@ contains use cmfd_header, only: allocate_cmfd +#ifdef PETSC integer :: new_comm ! new mpi communicator +#endif integer :: color ! color group of processor ! Read in cmfd input file @@ -34,17 +36,17 @@ contains end if ! Split up procs -# ifdef PETSC +#ifdef PETSC call MPI_COMM_SPLIT(MPI_COMM_WORLD, color, 0, new_comm, mpi_err) -# endif +#endif ! assign to PETSc -# ifdef PETSC +#ifdef PETSC PETSC_COMM_WORLD = new_comm ! Initialize PETSc on all procs call PetscInitialize(PETSC_NULL_CHARACTER, mpi_err) -# endif +#endif ! Initialize timers call time_cmfd % reset() From 740b60e72ef4c2d8ae9ffff59e252bcadf894282 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 5 Dec 2013 20:07:14 -0500 Subject: [PATCH 03/31] sample_external_source now checks to make sure the sampled xyz exists within the problem domain. If it does not, the source site is resampled. --- src/DEPENDENCIES | 1 + src/source.F90 | 28 ++++++++++++++++++++++------ 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index b4fe997bbc..579bf26fb1 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -329,6 +329,7 @@ solver_interface.o: vector_header.o source.o: bank_header.o source.o: constants.o source.o: error.o +source.o: geometry.o source.o: geometry_header.o source.o: global.o source.o: math.o diff --git a/src/source.F90 b/src/source.F90 index 2082ad6dda..3d4de842a0 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -3,6 +3,7 @@ module source use bank_header, only: Bank use constants use error, only: fatal_error + use geometry, only: find_cell use geometry_header, only: BASE_UNIVERSE use global use math, only: maxwell_spectrum, watt_spectrum @@ -73,6 +74,8 @@ contains real(8) :: p_max(3) ! maximum coordinates of source real(8) :: a ! Arbitrary parameter 'a' real(8) :: b ! Arbitrary parameter 'b' + logical :: found ! Does the source particle exist within geometry? + type(Particle) :: p ! Temporary particle for using find_cell ! Set weight to one by default site % wgt = ONE @@ -80,11 +83,24 @@ contains ! Sample position select case (external_source % type_space) case (SRC_SPACE_BOX) - ! Coordinates sampled uniformly over a box - p_min = external_source % params_space(1:3) - p_max = external_source % params_space(4:6) - r = (/ (prn(), i = 1,3) /) - site % xyz = p_min + r*(p_max - p_min) + ! Set particle defaults + call p % initialize() + ! Repeat sampling source location until a good site has been found + found = .false. + do while (.not.found) + ! Coordinates sampled uniformly over a box + p_min = external_source % params_space(1:3) + p_max = external_source % params_space(4:6) + r = (/ (prn(), i = 1,3) /) + site % xyz = p_min + r*(p_max - p_min) + + ! Fill p with needed data + p % coord0 % xyz = site % xyz + + ! Now search to see if location exists in geometry + call find_cell(p, found) + + end do case (SRC_SPACE_POINT) ! Point source @@ -146,7 +162,7 @@ contains end subroutine sample_external_source !=============================================================================== -! GET_SOURCE_PARTICLE returns the next source particle +! GET_SOURCE_PARTICLE returns the next source particle !=============================================================================== subroutine get_source_particle(p, index_source) From c3cd394f7a27968b2eeab61427cc4840746df139 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 5 Dec 2013 20:13:25 -0500 Subject: [PATCH 04/31] Prettied up the comments. --- src/source.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/source.F90 b/src/source.F90 index 3d4de842a0..57870206de 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -75,7 +75,7 @@ contains real(8) :: a ! Arbitrary parameter 'a' real(8) :: b ! Arbitrary parameter 'b' logical :: found ! Does the source particle exist within geometry? - type(Particle) :: p ! Temporary particle for using find_cell + type(Particle) :: p ! Temporary particle for using find_cell ! Set weight to one by default site % wgt = ONE From 7cc648036b6771ce1d8645d172ddecf56d8a0b10 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 9 Dec 2013 09:42:00 -0500 Subject: [PATCH 05/31] Added limit to number of external source resamples. --- src/constants.F90 | 20 ++++++++++++-------- src/source.F90 | 10 +++++++++- 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index e9999c0bdd..564e0baedd 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -20,7 +20,7 @@ module constants FILETYPE_PARTICLE_RESTART = -2 ! ============================================================================ - ! ADJUSTABLE PARAMETERS + ! ADJUSTABLE PARAMETERS ! NOTE: This is the only section of the constants module that should ever be ! adjusted. Modifying constants in other sections may cause the code to fail. @@ -50,6 +50,10 @@ module constants integer, parameter :: MAX_WORD_LEN = 150 integer, parameter :: MAX_FILE_LEN = 255 + ! Maximum number of external source spatial resamples to encounter before an + ! error is thrown. + integer, parameter :: MAX_EXTSRC_RESAMPLES = 10000 + ! ============================================================================ ! PHYSICAL CONSTANTS @@ -110,9 +114,9 @@ module constants ! Surface types integer, parameter :: & - SURF_PX = 1, & ! Plane parallel to x-plane - SURF_PY = 2, & ! Plane parallel to y-plane - SURF_PZ = 3, & ! Plane parallel to z-plane + SURF_PX = 1, & ! Plane parallel to x-plane + SURF_PY = 2, & ! Plane parallel to y-plane + SURF_PZ = 3, & ! Plane parallel to z-plane SURF_PLANE = 4, & ! Arbitrary plane SURF_CYL_X = 5, & ! Cylinder along x-axis SURF_CYL_Y = 6, & ! Cylinder along y-axis @@ -143,7 +147,7 @@ module constants ELECTRON = 3 ! Angular distribution type - integer, parameter :: & + integer, parameter :: & ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution ANGLE_32_EQUI = 2, & ! 32 equiprobable bins ANGLE_TABULAR = 3 ! Tabular angular distribution @@ -275,7 +279,7 @@ module constants SCORE_KAPPA_FISSION = -12, & ! fission energy production rate SCORE_CURRENT = -13, & ! partial current SCORE_EVENTS = -14 ! number of events - + ! Maximum scattering order supported integer, parameter :: SCATT_ORDER_MAX = 10 character(len=*), parameter :: SCATT_ORDER_MAX_PNSTR = "scatter-p10" @@ -322,7 +326,7 @@ module constants ! Source angular distribution types integer, parameter :: & - SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular + SRC_ANGLE_ISOTROPIC = 1, & ! Isotropic angular SRC_ANGLE_MONO = 2, & ! Monodirectional source SRC_ANGLE_TABULAR = 3 ! Tabular distribution @@ -332,7 +336,7 @@ module constants SRC_ENERGY_MAXWELL = 2, & ! Maxwell fission spectrum SRC_ENERGY_WATT = 3, & ! Watt fission spectrum SRC_ENERGY_TABULAR = 4 ! Tabular distribution - + ! ============================================================================ ! MISCELLANEOUS CONSTANTS diff --git a/src/source.F90 b/src/source.F90 index 57870206de..4fcdd0074a 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -76,6 +76,7 @@ contains real(8) :: b ! Arbitrary parameter 'b' logical :: found ! Does the source particle exist within geometry? type(Particle) :: p ! Temporary particle for using find_cell + integer, save :: num_resamples = 0 ! Number of resamples encountered ! Set weight to one by default site % wgt = ONE @@ -99,7 +100,14 @@ contains ! Now search to see if location exists in geometry call find_cell(p, found) - + if (.not. found) then + num_resamples = num_resamples + 1 + if (num_resamples == MAX_EXTSRC_RESAMPLES) then + message = "Maximum number of external source spatial resamples & + &reached!" + call fatal_error() + end if + end if end do case (SRC_SPACE_POINT) From 8884fb9e3f1657b2b5e0c5b4e8e4a3af5adfbe68 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 10 Dec 2013 05:58:11 -0500 Subject: [PATCH 06/31] Added changes to allow for the entire list of zaids to be checked. --- src/ace.F90 | 30 +++++++----- src/ace_header.F90 | 113 +++++++++++++++++++++++---------------------- src/output.F90 | 66 ++++++++++++++------------ 3 files changed, 112 insertions(+), 97 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index d54aa83c28..d5fa012d36 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -115,9 +115,9 @@ contains ! search through the list of nuclides for one which has a matching zaid sab => sab_tables(mat % i_sab_tables(k)) - ! Loop through nuclides and find match + ! Loop through nuclides and find match FIND_NUCLIDE: do j = 1, mat % n_nuclides - if (nuclides(mat % nuclide(j)) % zaid == sab % zaid) then + if (any(sab % zaid == nuclides(mat % nuclide(j)) % zaid)) then mat % i_sab_nuclides(k) = j exit FIND_NUCLIDE end if @@ -161,16 +161,16 @@ contains mat % i_sab_tables(m) = temp_table end do SORT_SAB end if - + ! Deallocate temporary arrays for names of nuclides and S(a,b) tables if (allocated(mat % names)) deallocate(mat % names) if (allocated(mat % sab_names)) deallocate(mat % sab_names) end do MATERIAL_LOOP2 - + ! Avoid some valgrind leak errors call already_read % clear() - + end subroutine read_xs !=============================================================================== @@ -253,7 +253,7 @@ contains end if ! Read more header and NXS and JXS - read(UNIT=in, FMT=100) comment, mat, & + read(UNIT=in, FMT=100) comment, mat, & (zaids(i), awrs(i), i=1,16), NXS, JXS 100 format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/& ,8I9/8I9/8I9/8I9/8I9/8I9) @@ -277,7 +277,7 @@ contains ACCESS='direct', RECL=record_length) ! Read all header information - read(UNIT=in, REC=location) name, awr, kT, date_, & + read(UNIT=in, REC=location) name, awr, kT, date_, & comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS ! determine table length @@ -333,7 +333,15 @@ contains sab % name = name sab % awr = awr sab % kT = kT - sab % zaid = zaids(1) + ! Find sab % n_zaid + do i = 1, 16 + if (zaids(i) == 0) then + sab % n_zaid = i - 1 + exit + end if + end do + allocate(sab % zaid(sab % n_zaid)) + sab % zaid = zaids(1: sab % n_zaid) call read_thermal_data(sab) end select @@ -406,7 +414,7 @@ contains integer :: LNU ! type of nu data (polynomial or tabular) integer :: NC ! number of polynomial coefficients integer :: NR ! number of interpolation regions - integer :: NE ! number of energies + integer :: NE ! number of energies integer :: NPCR ! number of delayed neutron precursor groups integer :: LED ! location of energy distribution locators integer :: LDIS ! location of all energy distributions @@ -809,7 +817,7 @@ contains LED = JXS(10) - ! Loop over all reactions + ! Loop over all reactions do i = 1, NXS(5) rxn => nuc % reactions(i+1) ! skip over elastic scattering rxn % has_energy_dist = .true. @@ -840,7 +848,7 @@ contains integer :: LDIS ! location of all energy distributions integer :: LNW ! location of next energy distribution if multiple - integer :: LAW ! secondary energy distribution law + integer :: LAW ! secondary energy distribution law integer :: NR ! number of interpolation regions integer :: NE ! number of incoming energies integer :: IDAT ! location of first energy distribution for given MT diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 5e2f39a007..26fbf4cbca 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -16,7 +16,7 @@ module ace_header integer, allocatable :: type(:) ! type of distribution integer, allocatable :: location(:) ! location of each table real(8), allocatable :: data(:) ! angular distribution data - + ! Type-Bound procedures contains procedure :: clear => distangle_clear ! Deallocates DistAngle @@ -35,7 +35,7 @@ module ace_header ! For reactions that may have multiple energy distributions such as (n.2n), ! this pointer allows multiple laws to be stored type(DistEnergy), pointer :: next => null() - + ! Type-Bound procedures contains procedure :: clear => distenergy_clear ! Deallocates DistEnergy @@ -57,7 +57,7 @@ module ace_header logical :: has_energy_dist ! Energy distribution present? type(DistAngle) :: adist ! Secondary angular distribution type(DistEnergy), pointer :: edist => null() ! Secondary energy distribution - + ! Type-Bound procedures contains procedure :: clear => reaction_clear ! Deallocates Reaction @@ -76,7 +76,7 @@ module ace_header logical :: multiply_smooth ! multiply by smooth cross section? real(8), allocatable :: energy(:) ! incident energies real(8), allocatable :: prob(:,:,:) ! actual probabibility tables - + ! Type-Bound procedures contains procedure :: clear => urrdata_clear ! Deallocates UrrData @@ -137,7 +137,7 @@ module ace_header ! Reactions integer :: n_reaction ! # of reactions type(Reaction), pointer :: reactions(:) => null() - + ! Type-Bound procedures contains procedure :: clear => nuclide_clear ! Deallocates Nuclide @@ -147,12 +147,13 @@ module ace_header ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc !=============================================================================== - + type SAlphaBeta - character(10) :: name ! name of table, e.g. lwtr.10t - integer :: zaid ! Z and A identifier, e.g. 6012 for Carbon-12 - real(8) :: awr ! weight of nucleus in neutron masses - real(8) :: kT ! temperature in MeV (k*T) + character(10) :: name ! name of table, e.g. lwtr.10t + real(8) :: awr ! weight of nucleus in neutron masses + real(8) :: kT ! temperature in MeV (k*T) + integer :: n_zaid ! Number of valid zaids + integer, allocatable :: zaid(:) ! List of valid Z and A identifiers, e.g. 6012 ! threshold for S(a,b) treatment (usually ~4 eV) real(8) :: threshold_inelastic @@ -164,7 +165,7 @@ module ace_header integer :: n_inelastic_mu ! # of outgoing angles for inelastic integer :: secondary_mode ! secondary mode (equal/skewed) real(8), allocatable :: inelastic_e_in(:) - real(8), allocatable :: inelastic_sigma(:) + real(8), allocatable :: inelastic_sigma(:) real(8), allocatable :: inelastic_e_out(:,:) real(8), allocatable :: inelastic_mu(:,:,:) @@ -240,125 +241,125 @@ module ace_header !=============================================================================== ! DISTANGLE_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - +!=============================================================================== + subroutine distangle_clear(this) - + class(DistAngle), intent(inout) :: this ! The DistAngle object to clear - + if (allocated(this % energy)) & deallocate(this % energy, this % type, this % location, this % data) - - end subroutine distangle_clear + + end subroutine distangle_clear !=============================================================================== ! DISTENERGY_CLEAR resets and deallocates data in DistEnergy. -!=============================================================================== - +!=============================================================================== + recursive subroutine distenergy_clear(this) - + class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear - + ! Clear p_valid call this % p_valid % clear() - + if (allocated(this % data)) & deallocate(this % data) - + if (associated(this % next)) then ! recursively clear this item call this % next % clear() deallocate(this % next) end if - + end subroutine distenergy_clear - + !=============================================================================== ! REACTION_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - +!=============================================================================== + subroutine reaction_clear(this) - + class(Reaction), intent(inout) :: this ! The Reaction object to clear - + if (allocated(this % sigma)) & deallocate(this % sigma) - + if (associated(this % edist)) then call this % edist % clear() deallocate(this % edist) end if - + call this % adist % clear() - - end subroutine reaction_clear - + + end subroutine reaction_clear + !=============================================================================== ! URRDATA_CLEAR resets and deallocates data in Reaction. -!=============================================================================== - +!=============================================================================== + subroutine urrdata_clear(this) - + class(UrrData), intent(inout) :: this ! The UrrData object to clear - + if (allocated(this % energy)) & deallocate(this % energy, this % prob) - - end subroutine urrdata_clear + + end subroutine urrdata_clear !=============================================================================== ! NUCLIDE_CLEAR resets and deallocates data in Nuclide. -!=============================================================================== - +!=============================================================================== + subroutine nuclide_clear(this) - + class(Nuclide), intent(inout) :: this ! The Nuclide object to clear - + integer :: i ! Loop counter - + if (allocated(this % grid_index)) & deallocate(this % grid_index) - + if (allocated(this % energy)) & deallocate(this % total, this % elastic, this % fission, & this % nu_fission, this % absorption) if (allocated(this % heating)) & deallocate(this % heating) - + if (allocated(this % index_fission)) & deallocate(this % index_fission) - + if (allocated(this % nu_t_data)) & deallocate(this % nu_t_data) - + if (allocated(this % nu_p_data)) & deallocate(this % nu_p_data) - + if (allocated(this % nu_d_data)) & deallocate(this % nu_d_data) - + if (allocated(this % nu_d_precursor_data)) & deallocate(this % nu_d_precursor_data) - + if (associated(this % nu_d_edist)) then do i = 1, size(this % nu_d_edist) call this % nu_d_edist(i) % clear() end do deallocate(this % nu_d_edist) end if - + if (associated(this % urr_data)) then call this % urr_data % clear() deallocate(this % urr_data) end if - + if (associated(this % reactions)) then do i = 1, size(this % reactions) call this % reactions(i) % clear() end do deallocate(this % reactions) end if - - end subroutine nuclide_clear + + end subroutine nuclide_clear end module ace_header diff --git a/src/output.F90 b/src/output.F90 index c60b96bf26..e4fcfecef5 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -126,7 +126,7 @@ contains ! print header based on level select case (header_level) case (1) - write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & + write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & repeat('=', n) // '> ' // trim(line) // ' <' // & repeat('=', m), repeat('=', 75) case (2) @@ -180,7 +180,7 @@ contains end subroutine print_usage !=============================================================================== -! WRITE_MESSAGE displays an informational message to the log file and the +! WRITE_MESSAGE displays an informational message to the log file and the ! standard output stream. !=============================================================================== @@ -215,7 +215,7 @@ contains ! Determine last space in current line last_space = index(message(i_start+1:i_start+line_wrap), & ' ', BACK=.true.) - if (last_space == 0) then + if (last_space == 0) then i_end = min(length + 1, i_start+line_wrap) - 1 write(ou, fmt='(1X,A)') message(i_start+1:i_end) else @@ -1030,6 +1030,7 @@ contains integer :: size_sab ! memory used by S(a,b) table integer :: unit_ ! unit to write to + integer :: i ! Loop counter for parsing through sab % zaid ! set default unit for writing information if (present(unit)) then @@ -1040,7 +1041,12 @@ contains ! Basic S(a,b) table information write(unit_,*) 'S(a,b) Table ' // trim(sab % name) - write(unit_,*) ' zaid = ' // trim(to_str(sab % zaid)) + write(unit_,'(A)',advance="no") ' zaids = ' + do i = 1, sab % n_zaid - 1 + ! This could use some fancy handling for if the line becomes too long + write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(sab % n_zaid))) // ', ' + end do + write(unit_,*) trim(to_str(sab % zaid(sab % n_zaid))) write(unit_,*) ' awr = ' // trim(to_str(sab % awr)) write(unit_,*) ' kT = ' // trim(to_str(sab % kT)) @@ -1242,7 +1248,7 @@ contains end select end if write(UNIT=ou, FMT=*) - + write(UNIT=ou, FMT='(2X,A9,3X)', ADVANCE='NO') "=========" write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========" if (entropy_on) write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========" @@ -1272,7 +1278,7 @@ contains if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & entropy(overall_gen) - if (overall_gen - n_inactive*gen_per_batch > 1) then + if (overall_gen - n_inactive*gen_per_batch > 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & keff, keff_std end if @@ -1300,7 +1306,7 @@ contains entropy(current_batch*gen_per_batch) ! write out accumulated k-effective if after first active batch - if (overall_gen - n_inactive*gen_per_batch > 1) then + if (overall_gen - n_inactive*gen_per_batch > 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & keff, keff_std else @@ -1310,7 +1316,7 @@ contains ! write out cmfd keff if it is active and other display info if (cmfd_on) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % k_cmfd(current_batch) + cmfd % k_cmfd(current_batch) select case(trim(cmfd_display)) case('entropy') write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & @@ -1349,7 +1355,7 @@ contains ! Plot id write(ou,100) "Plot ID:", trim(to_str(pl % id)) - + ! Plot type if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Plot Type:", "Slice" @@ -1387,7 +1393,7 @@ contains trim(to_str(pl % pixels(2))) else if (pl % type == PLOT_TYPE_VOXEL) then write(ou,100) "Voxels:", trim(to_str(pl % pixels(1))) // " " // & - trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3))) + trim(to_str(pl % pixels(2))) // " " // trim(to_str(pl % pixels(3))) end if write(ou,*) @@ -1539,7 +1545,7 @@ contains write(ou,100) 'Cell ID','No. Overlap Checks' - do i = 1, n_cells + do i = 1, n_cells write(ou,101) cells(i) % id, overlap_check_cnt(i) if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1 end do @@ -1573,7 +1579,7 @@ contains integer :: k ! loop index for scoring bins integer :: n ! loop index for nuclides integer :: l ! loop index for user scores - integer :: type ! type of tally filter + integer :: type ! type of tally filter integer :: indent ! number of spaces to preceed output integer :: filter_index ! index in results array for filters integer :: score_index ! scoring bin index @@ -1749,13 +1755,13 @@ contains score_name = 'P' // trim(to_str(t % scatt_order(k))) // & ' Scattering Moment' end if - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) else if (t % score_bins(k) == SCORE_SCATTER_PN) then score_name = "Scattering Rate" - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) @@ -1763,7 +1769,7 @@ contains score_index = score_index + 1 score_name = 'P' // trim(to_str(n_order)) // & ' Scattering Moment' - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) @@ -1775,7 +1781,7 @@ contains else score_name = score_names(abs(t % score_bins(k))) end if - write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), score_name, & to_str(t % results(score_index,filter_index) % sum), & trim(to_str(t % results(score_index,filter_index) % sum_sq)) @@ -1813,8 +1819,8 @@ contains integer :: i_filter_ein ! index for incoming energy filter integer :: i_filter_surf ! index for surface filter integer :: n ! number of incoming energy bins - integer :: len1 ! length of string - integer :: len2 ! length of string + integer :: len1 ! length of string + integer :: len2 ! length of string integer :: filter_index ! index in results array for filters logical :: print_ebin ! should incoming energy bin be displayed? character(MAX_LINE_LEN) :: string @@ -1864,14 +1870,14 @@ contains mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Left", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1881,14 +1887,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_RIGHT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Right", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1898,14 +1904,14 @@ contains mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Back", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1915,14 +1921,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Front", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_FRONT filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Front", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1932,14 +1938,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.) matching_bins(i_filter_surf) = IN_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Bottom", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) @@ -1949,14 +1955,14 @@ contains mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.) matching_bins(i_filter_surf) = IN_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Top", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) matching_bins(i_filter_surf) = OUT_TOP filter_index = sum((matching_bins(1:t%n_filters) - 1) * t % stride) + 1 - write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & + write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Top", & to_str(t % results(1,filter_index) % sum), & trim(to_str(t % results(1,filter_index) % sum_sq)) From ee1edca155ab185f664aa6805ee81c671af0ac53 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 10 Dec 2013 09:06:40 -0500 Subject: [PATCH 07/31] fixed source_write option in input reading --- src/input_xml.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 05f4ae0c56..8cd93f0e10 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -635,8 +635,8 @@ contains if (check_for_node(node_sp, "source_write")) then call get_node_value(node_sp, "source_write", temp_str) call lower_case(temp_str) - if (trim(temp_str) == 'true' .or. & - trim(temp_str) == '1') source_separate = .true. + if (trim(temp_str) == 'false' .or. & + trim(temp_str) == '0') source_write = .false. end if else ! If no tag was present, by default write state point at From b6395729a4ba6e33b321fd08635d8f73baeb307f Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 10 Dec 2013 11:51:40 -0500 Subject: [PATCH 08/31] Added better printing of the new zaid list for cross_sections.out; now it will go to a new line before hitting the 80-char limit. --- src/output.F90 | 35 +++++++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index e4fcfecef5..1507aa4cd5 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1028,9 +1028,10 @@ contains type(SAlphaBeta), pointer :: sab integer, optional :: unit - integer :: size_sab ! memory used by S(a,b) table - integer :: unit_ ! unit to write to - integer :: i ! Loop counter for parsing through sab % zaid + integer :: size_sab ! memory used by S(a,b) table + integer :: unit_ ! unit to write to + integer :: i ! Loop counter for parsing through sab % zaid + integer :: char_count ! Counter for the number of characters on a line ! set default unit for writing information if (present(unit)) then @@ -1041,12 +1042,30 @@ contains ! Basic S(a,b) table information write(unit_,*) 'S(a,b) Table ' // trim(sab % name) - write(unit_,'(A)',advance="no") ' zaids = ' - do i = 1, sab % n_zaid - 1 - ! This could use some fancy handling for if the line becomes too long - write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(sab % n_zaid))) // ', ' + write(unit_,'(A)',advance="no") ' zaids = ' + ! Initialize the counter based on the above string + char_count = 11 + do i = 1, sab % n_zaid + ! Deal with a line thats too long + if (char_count >= 73) then ! 73 = 80 - (5 ZAID chars + 1 space + 1 comma) + ! End the line + write(unit_,*) "" + ! Add 11 leading blanks + write(unit_,'(A)', advance="no") " " + ! reset the counter to 11 + char_count = 11 + end if + if (i < sab % n_zaid) then + ! Include a comma + write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) // ", " + char_count = char_count + len(trim(to_str(sab % zaid(i)))) + 2 + else + ! Don't include a comma, since we are all done + write(unit_,'(A)',advance="no") trim(to_str(sab % zaid(i))) + end if + end do - write(unit_,*) trim(to_str(sab % zaid(sab % n_zaid))) + write(unit_,*) "" ! Move to next line write(unit_,*) ' awr = ' // trim(to_str(sab % awr)) write(unit_,*) ' kT = ' // trim(to_str(sab % kT)) From 0a05d4594aeb57a0179048c3bb4ccdd00a990b31 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Dec 2013 16:53:11 -0500 Subject: [PATCH 09/31] Fix bug in doppler broadening (not used). --- src/doppler.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/doppler.F90 b/src/doppler.F90 index 12fb34bbfe..e888eb2b18 100644 --- a/src/doppler.F90 +++ b/src/doppler.F90 @@ -87,7 +87,7 @@ contains ! ======================================================================= ! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE - if (k == 1 .and. a >= 4.0) then + if (k == 1 .and. a >= -4.0) then ! Since x = 0, this implies that a = -y F_b = F_a a = -y @@ -101,7 +101,7 @@ contains end if ! ======================================================================= - ! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4 + ! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4 k = i b = ZERO From 1bdf693d798b7b4597d65bee5753296bb8829823 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Dec 2013 22:33:59 -0500 Subject: [PATCH 10/31] Fix use statements in several cmfd modules. --- src/cmfd_data.F90 | 2 +- src/cmfd_execute.F90 | 12 +++++++----- src/cmfd_power_solver.F90 | 2 ++ 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index d0e80aaa90..8dd0fa35ff 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -934,7 +934,7 @@ contains subroutine compute_effective_downscatter() use constants, only: ZERO, CMFD_NOACCEL - use global, only: cmfd, cmfd_downscatter + use global, only: cmfd integer :: nx ! number of mesh cells in x direction integer :: ny ! number of mesh cells in y direction diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 7a35a12144..6234d9fcf9 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -117,9 +117,9 @@ contains subroutine process_cmfd_options() +#ifdef PETSC use global, only: cmfd_snes_monitor, cmfd_ksp_monitor, mpi_err -#ifdef PETSC ! Check for snes monitor if (cmfd_snes_monitor) call PetscOptionsSetValue("-snes_monitor", & "stdout", mpi_err) @@ -138,9 +138,10 @@ contains subroutine calc_fission_source() use constants, only: CMFD_NOACCEL, ZERO, TWO - use global, only: cmfd, cmfd_coremap, master, mpi_err, entropy_on, & - current_batch + use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch + #ifdef MPI + use global, only: mpi_err use mpi #endif @@ -261,13 +262,14 @@ contains use constants, only: ZERO, ONE use error, only: warning, fatal_error - use global, only: n_particles, meshes, source_bank, work, & - n_user_meshes, message, cmfd, master, mpi_err + use global, only: meshes, source_bank, work, n_user_meshes, message, & + cmfd, master use mesh_header, only: StructuredMesh use mesh, only: count_bank_sites, get_mesh_indices use search, only: binary_search #ifdef MPI + use global, only: mpi_err use mpi #endif diff --git a/src/cmfd_power_solver.F90 b/src/cmfd_power_solver.F90 index c1309b72cd..f2761d0a1f 100644 --- a/src/cmfd_power_solver.F90 +++ b/src/cmfd_power_solver.F90 @@ -102,7 +102,9 @@ contains subroutine init_data(adjoint) use constants, only: ONE, ZERO +#ifdef PETSC use global, only: cmfd_write_matrices +#endif logical, intent(in) :: adjoint ! adjoint calcualtion From ca4ec3b3cd1e2dd6d1a07e461b2f5ef504262fc3 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 11 Dec 2013 09:47:49 -0500 Subject: [PATCH 11/31] Added ZrH as the top reflector of the test_salphabeta_multiple test in order to include an S(a,b) library which can be applied to multiple problems in our test suite. --- tests/test_salphabeta_multiple/geometry.xml | 6 ++++-- tests/test_salphabeta_multiple/materials.xml | 18 +++++++++++++++++- .../test_salphabeta_multiple/results_true.dat | 2 +- 3 files changed, 22 insertions(+), 4 deletions(-) diff --git a/tests/test_salphabeta_multiple/geometry.xml b/tests/test_salphabeta_multiple/geometry.xml index 2d427d8cb3..63f69f7438 100644 --- a/tests/test_salphabeta_multiple/geometry.xml +++ b/tests/test_salphabeta_multiple/geometry.xml @@ -44,9 +44,11 @@ - + - + + + diff --git a/tests/test_salphabeta_multiple/materials.xml b/tests/test_salphabeta_multiple/materials.xml index 97ff56c05f..6e975bde97 100644 --- a/tests/test_salphabeta_multiple/materials.xml +++ b/tests/test_salphabeta_multiple/materials.xml @@ -1,7 +1,8 @@ - + 70c @@ -111,4 +112,19 @@ + + + + + + + + + + + + + + + diff --git a/tests/test_salphabeta_multiple/results_true.dat b/tests/test_salphabeta_multiple/results_true.dat index f26a309aae..bbd61cdb53 100644 --- a/tests/test_salphabeta_multiple/results_true.dat +++ b/tests/test_salphabeta_multiple/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.013747E+00 3.067630E-02 +9.761880E-01 9.170415E-03 From 549bdad2e371e6afb78353b7dec8914bacbdaee8 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 14 Dec 2013 10:18:35 -0500 Subject: [PATCH 12/31] Added ability to read in continuous secondary energy distributions for inelastic S(a,b) collisions. Still have to implement new collision physics in physics.F90 --- src/ace.F90 | 83 ++++++++++++++++++++++++++++++++-------------- src/ace_header.F90 | 24 ++++++++++++-- src/constants.F90 | 3 +- 3 files changed, 82 insertions(+), 28 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index d5fa012d36..ff8582470d 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -1195,15 +1195,10 @@ contains integer :: NE_out ! number of outgoing energies integer :: NMU ! number of outgoing angles integer :: JXS4 ! location of elastic energy table + real(8), allocatable :: LOCC(:) ! Location of inelastic data - ! read secondary energy mode for inelastic scattering and check + ! read secondary energy mode for inelastic scattering table % secondary_mode = NXS(7) - if (table % secondary_mode /= SAB_SECONDARY_EQUAL .and. & - table % secondary_mode /= SAB_SECONDARY_SKEWED) then - message = "Unsupported secondary mode on S(a,b) table " // & - trim(adjustl(table % name)) // ": " // to_str(table % secondary_mode) - call fatal_error() - end if ! read number of inelastic energies and allocate arrays NE_in = int(XSS(JXS(1))) @@ -1221,29 +1216,67 @@ contains ! allocate space for outgoing energy/angle for inelastic ! scattering - NE_out = NXS(4) - NMU = NXS(3) + 1 - table % n_inelastic_e_out = NE_out - table % n_inelastic_mu = NMU - allocate(table % inelastic_e_out(NE_out, NE_in)) - allocate(table % inelastic_mu(NMU, NE_out, NE_in)) + if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. & + table % secondary_mode == SAB_SECONDARY_SKEWED) then + NMU = NXS(3) + 1 + table % n_inelastic_mu = NMU + NE_out = NXS(4) + table % n_inelastic_e_out = NE_out + allocate(table % inelastic_e_out(NE_out, NE_in)) + allocate(table % inelastic_mu(NMU, NE_out, NE_in)) + else if (table % secondary_mode == SAB_SECONDARY_CONT) then + NMU = NXS(3) - 1 + table % n_inelastic_mu = NMU + allocate(table % inelastic_data(NE_in)) + allocate(LOCC(NE_in)) + ! NE_out will be determined later + end if ! read outgoing energy/angle distribution for inelastic scattering - lc = JXS(3) - 1 - do i = 1, NE_in - do j = 1, NE_out - ! read outgoing energy - table % inelastic_e_out(j,i) = XSS(lc + 1) + if (table % secondary_mode == SAB_SECONDARY_EQUAL .or. & + table % secondary_mode == SAB_SECONDARY_SKEWED) then + lc = JXS(3) - 1 + do i = 1, NE_in + do j = 1, NE_out + ! read outgoing energy + table % inelastic_e_out(j,i) = XSS(lc + 1) - ! read outgoing angles for this outgoing energy - do k = 1, NMU - table % inelastic_mu(k,j,i) = XSS(lc + 1 + k) + ! read outgoing angles for this outgoing energy + do k = 1, NMU + table % inelastic_mu(k,j,i) = XSS(lc + 1 + k) + end do + + ! advance pointer + lc = lc + 1 + NMU end do - - ! advance pointer - lc = lc + 1 + NMU end do - end do + else if (table % secondary_mode == SAB_SECONDARY_CONT) then + ! Get the location pointers to each Ein's DistEnergySAB data + LOCC = get_real(NE_in) + ! Get the number of outgoing energies and allocate space accordingly + do i = 1, NE_in + NE_out = XSS(XSS_index + i - 1) + table % inelastic_data(i) % n_e_out = NE_out + allocate(table % inelastic_data(i) % e_out (NE_out)) + allocate(table % inelastic_data(i) % e_out_pdf (NE_out)) + allocate(table % inelastic_data(i) % e_out_cdf (NE_out)) + allocate(table % inelastic_data(i) % mu (NMU, NE_out)) + end do + + ! Now we can fill the inelastic_data(i) attributes + do i = 1, NE_in + XSS_index = LOCC(i) + NE_out = table % inelastic_data(i) % n_e_out + do j = 1, NE_out + table % inelastic_data(i) % e_out(j) = XSS(XSS_index + 1) + table % inelastic_data(i) % e_out_pdf(j) = XSS(XSS_index + 2) + table % inelastic_data(i) % e_out_cdf(j) = XSS(XSS_index + 3) + table % inelastic_data(i) % mu(:, j) = & + XSS(XSS_index + 4: XSS_index + 4 + NMU - 1) + XSS_index = XSS_index + 4 + NMU - 1 + end do + end do + end if ! read number of elastic energies and allocate arrays JXS4 = JXS(4) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 26fbf4cbca..e5bf1bb3a4 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -32,7 +32,7 @@ module ace_header type(Tab1) :: p_valid ! probability of law validity real(8), allocatable :: data(:) ! energy distribution data - ! For reactions that may have multiple energy distributions such as (n.2n), + ! For reactions that may have multiple energy distributions such as (n,2n), ! this pointer allows multiple laws to be stored type(DistEnergy), pointer :: next => null() @@ -143,6 +143,20 @@ module ace_header procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide +!=============================================================================== +! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic +! thermal scattering collisions which utilize a continuous secondary energy +! representation. +!=============================================================================== + + type DistEnergySab + integer :: n_e_out + real(8), allocatable :: e_out(:) + real(8), allocatable :: e_out_pdf(:) + real(8), allocatable :: e_out_cdf(:) + real(8), allocatable :: mu(:,:) + end type DistEnergySab + !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc @@ -163,11 +177,17 @@ module ace_header integer :: n_inelastic_e_in ! # of incoming E for inelastic integer :: n_inelastic_e_out ! # of outgoing E for inelastic integer :: n_inelastic_mu ! # of outgoing angles for inelastic - integer :: secondary_mode ! secondary mode (equal/skewed) + integer :: secondary_mode ! secondary mode (equal/skewed/continuous) real(8), allocatable :: inelastic_e_in(:) real(8), allocatable :: inelastic_sigma(:) + ! The following are used only if secondary_mode is 0 or 1 real(8), allocatable :: inelastic_e_out(:,:) real(8), allocatable :: inelastic_mu(:,:,:) + ! The following is used only if secondary_mode is 3 + ! The different implementation is necessary because the continuous + ! representation has a variable number of outgoing energy points for each + ! incoming energy + type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein ! Elastic scattering data integer :: elastic_mode ! elastic mode (discrete/exact) diff --git a/src/constants.F90 b/src/constants.F90 index 564e0baedd..6fea8cbc16 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -155,7 +155,8 @@ module constants ! Secondary energy mode for S(a,b) inelastic scattering integer, parameter :: & SAB_SECONDARY_EQUAL = 0, & ! Equally-likely outgoing energy bins - SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins + SAB_SECONDARY_SKEWED = 1, & ! Skewed outgoing energy bins + SAB_SECONDARY_CONT = 2 ! Continuous, linear-linear interpolation ! Elastic mode for S(a,b) elastic scattering integer, parameter :: & From 8f5ce5ed5ca6ff4913c7e8b80def71c673fd98a2 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 15 Dec 2013 13:45:18 -0500 Subject: [PATCH 13/31] Corrected an error message in physics.F90:sample_energy, which referred to Law 44 instead of the correct Law 61. --- src/physics.F90 | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index fb1f891943..5771377c86 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -47,7 +47,7 @@ contains if (verbosity >= 10 .or. trace) then message = " " // trim(reaction_name(p % event_MT)) // " with " // & trim(adjustl(nuclides(p % event_nuclide) % name)) // & - ". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV." + ". Energy = " // trim(to_str(p % E * 1e6_8)) // " eV." call write_message() end if @@ -158,7 +158,7 @@ contains call fatal_error() end if - ! Find atom density + ! Find atom density i_nuclide = mat % nuclide(i) atom_density = mat % atom_density(i) @@ -209,7 +209,7 @@ contains i_reaction = nuc % index_fission(1) return end if - + ! Get grid index and interpolatoin factor and sample fission cdf i_grid = micro_xs(i_nuclide) % index_grid f = micro_xs(i_nuclide) % interp_factor @@ -226,7 +226,7 @@ contains if (i_grid < rxn % threshold) cycle ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + f*(rxn%sigma(i_grid - rxn%threshold + 2))) ! Create fission bank sites if fission occus @@ -382,7 +382,7 @@ contains if (i_grid < rxn % threshold) cycle ! add to cumulative probability - prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) & + f*(rxn%sigma(i_grid - rxn%threshold + 2))) end do @@ -513,7 +513,7 @@ contains f = ZERO else i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E) - f = (E - sab%elastic_e_in(i)) / & + f = (E - sab%elastic_e_in(i)) / & (sab%elastic_e_in(i+1) - sab%elastic_e_in(i)) end if @@ -563,7 +563,7 @@ contains f = ZERO else i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E) - f = (E - sab%inelastic_e_in(i)) / & + f = (E - sab%inelastic_e_in(i)) / & (sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i)) end if @@ -974,7 +974,7 @@ contains E_cm = E ! determine outgoing energy in lab - E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & + E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & / ((A+ONE)*(A+ONE)) ! determine outgoing angle in lab @@ -1037,7 +1037,7 @@ contains r = ONE else i = binary_search(rxn % adist % energy, n, E) - r = (E - rxn % adist % energy(i)) / & + r = (E - rxn % adist % energy(i)) / & (rxn % adist % energy(i+1) - rxn % adist % energy(i)) end if @@ -1166,7 +1166,7 @@ contains end if end function rotate_angle - + !=============================================================================== ! SAMPLE_ENERGY samples an outgoing energy distribution, either for a secondary ! neutron from a collision or for a prompt/delayed fission neutron @@ -1322,7 +1322,7 @@ contains ! ======================================================================= ! CONTINUOUS TABULAR DISTRIBUTION - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR == 1) then @@ -1348,7 +1348,7 @@ contains r = ONE else i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & + r = (E_in - edist%data(lc+i)) / & (edist%data(lc+i+1) - edist%data(lc+i)) end if @@ -1452,7 +1452,7 @@ contains ! ======================================================================= ! MAXWELL FISSION SPECTRUM - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) @@ -1484,7 +1484,7 @@ contains ! ======================================================================= ! EVAPORATION SPECTRUM - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) @@ -1565,7 +1565,7 @@ contains call fatal_error() end if - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR > 0) then @@ -1587,7 +1587,7 @@ contains r = ONE else i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & + r = (E_in - edist%data(lc+i)) / & (edist%data(lc+i+1) - edist%data(lc+i)) end if @@ -1714,11 +1714,11 @@ contains if (.not. present(mu_out)) then ! call write_particle_restart() - message = "Law 44 called without giving mu_out as argument." + message = "Law 61 called without giving mu_out as argument." call fatal_error() end if - ! read number of interpolation regions and incoming energies + ! read number of interpolation regions and incoming energies NR = int(edist % data(1)) NE = int(edist % data(2 + 2*NR)) if (NR > 0) then @@ -1740,7 +1740,7 @@ contains r = ONE else i = binary_search(edist % data(lc+1:lc+NE), NE, E_in) - r = (E_in - edist%data(lc+i)) / & + r = (E_in - edist%data(lc+i)) / & (edist%data(lc+i+1) - edist%data(lc+i)) end if From 1dd13d80641407ca8217635f3a8ecbdb93908e16 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 15 Dec 2013 14:30:29 -0500 Subject: [PATCH 14/31] Added code to physics.F90%sab_scatter() to sample continuous secondary energy distributions; still have to test. --- src/physics.F90 | 190 +++++++++++++++++++++++++++++++++++++----------- 1 file changed, 146 insertions(+), 44 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 5771377c86..372d67bd3e 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -490,13 +490,24 @@ contains integer :: k ! outgoing cosine bin integer :: n_energy_out ! number of outgoing energy bins real(8) :: f ! interpolation factor - real(8) :: r ! used for skewed sampling + real(8) :: r ! used for skewed sampling & continuous real(8) :: E_ij ! outgoing energy j for E_in(i) real(8) :: E_i1j ! outgoing energy j for E_in(i+1) real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j) real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j) real(8) :: prob ! probability for sampling Bragg edge type(SAlphaBeta), pointer, save :: sab => null() + ! Following are needed only for SAB_SECONDARY_CONT scattering + integer :: l ! sampled incoming E bin (is i or i + 1) + real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_J ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_J ! endpoints interpolated between i and i+1 + real(8) :: E_l_j, E_l_j1 ! adjacent E on outgoing grid l + real(8) :: p_l_j, p_l_j1 ! adjacent p on outgoing grid l + real(8) :: c_j, c_j1 ! cumulative probability + real(8) :: frac ! interpolation factor on outgoing energy + real(8) :: r1 ! RNG for outgoing energy + !$omp threadprivate(sab) ! Get pointer to S(a,b) table @@ -554,8 +565,7 @@ contains ! Outgoing energy is same as incoming energy -- no need to do anything else - ! Determine number of outgoing energy and angle bins - n_energy_out = sab % n_inelastic_e_out + ! Perform inelastic calculations ! Get index and interpolation factor for inelastic grid if (E < sab % inelastic_e_in(1)) then @@ -570,54 +580,146 @@ contains ! Now that we have an incoming energy bin, we need to determine the ! outgoing energy bin. This will depend on the "secondary energy ! mode". If the mode is 0, then the outgoing energy bin is chosen from a - ! set of equally-likely bins. However, if the mode is 1, then the first + ! set of equally-likely bins. If the mode is 1, then the first ! two and last two bins are skewed to have lower probabilities than the ! other bins (0.1 for the first and last bins and 0.4 for the second and - ! second to last bins, relative to a normal bin probability of 1) + ! second to last bins, relative to a normal bin probability of 1). + ! Finally, if the mode is 2, then a continuous distribution (with + ! accompanying PDF and CDF is utilized) - if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then - ! All bins equally likely - j = 1 + int(prn() * n_energy_out) - elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then - r = prn() * (n_energy_out - 3) - if (r > ONE) then - ! equally likely N-4 middle bins - j = int(r) + 2 - elseif (r > 0.6) then - ! second to last bin has relative probability of 0.4 - j = n_energy_out - 1 - elseif (r > 0.5) then - ! last bin has relative probability of 0.1 - j = n_energy_out - elseif (r > 0.1) then - ! second bin has relative probability of 0.4 - j = 2 - else - ! first bin has relative probability of 0.1 - j = 1 + if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. & + (sab % secondary_mode == SAB_SECONDARY_SKEWED)) then + if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then + ! All bins equally likely + + j = 1 + int(prn() * sab % n_inelastic_e_out) + elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then + ! Distribution skewed away from edge points + + ! Determine number of outgoing energy and angle bins + n_energy_out = sab % n_inelastic_e_out + + r = prn() * (n_energy_out - 3) + if (r > ONE) then + ! equally likely N-4 middle bins + j = int(r) + 2 + elseif (r > 0.6) then + ! second to last bin has relative probability of 0.4 + j = n_energy_out - 1 + elseif (r > 0.5) then + ! last bin has relative probability of 0.1 + j = n_energy_out + elseif (r > 0.1) then + ! second bin has relative probability of 0.4 + j = 2 + else + ! first bin has relative probability of 0.1 + j = 1 + end if end if + + ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) + E_ij = sab % inelastic_e_out(j,i) + E_i1j = sab % inelastic_e_out(j,i+1) + + ! Outgoing energy + E = (1 - f)*E_ij + f*E_i1j + + ! Sample outgoing cosine bin + k = 1 + int(prn() * sab % n_inelastic_mu) + + ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) + mu_ijk = sab % inelastic_mu(k,j,i) + mu_i1jk = sab % inelastic_mu(k,j,i+1) + + ! Cosine of angle between incoming and outgoing neutron + mu = (1 - f)*mu_ijk + f*mu_i1jk + + else if (sab % secondary_mode == SAB_SECONDARY_CONT) then + ! Continuous secondary energy - this is to be similar to + ! Law 61 interpolation on outgoing energy + + ! Sample between ith and (i+1)th bin + r = prn() + if (f > r) then + l = i + 1 + else + l = i + end if + + ! Determine endpoints on grid i + n_energy_out = sab % inelastic_data(i) % n_e_out + E_i_1 = sab % inelastic_data(i) % e_out(1) + E_i_J = sab % inelastic_data(i) % e_out(n_energy_out) + + ! Determine endpoints on grid i + 1 + n_energy_out = sab % inelastic_data(i + 1) % n_e_out + E_i1_1 = sab % inelastic_data(i + 1) % e_out(1) + E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out) + + E_1 = E_i_1 + f * (E_i1_1 - E_i_1) + E_J = E_i_J + f * (E_i1_J - E_i_J) + + ! Determine outgoing energy bin + ! (First reset n_energy_out to the right value) + n_energy_out = sab % inelastic_data(l) % n_e_out + r1 = prn() + c_j = sab % inelastic_data(l) % e_out_cdf(1) + do j = 1, n_energy_out - 1 + c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1) + if (r1 < c_j1) exit + c_j = c_j1 + end do + + ! check to make sure k is <= n_energy_out - 1 + j = min(j, n_energy_out - 1) + + ! Get the data to interpolate between + E_l_j = sab % inelastic_data(l) % e_out(j) + p_l_j = sab % inelastic_data(l) % e_out_pdf(j) + + ! Next part assumes linear-linear interpolation in standard + E_l_j1 = sab % inelastic_data(l) % e_out(j + 1) + p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1) + + ! Find secondary energy (variable E) + frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j) + if (frac == ZERO) then + E = E_l_j + (r1 - c_j) / p_l_j + else + E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + & + TWO * frac * (r1 - c_j))) - p_l_j) / frac + end if + + ! Now interpolate between incident energy bins i and i + 1 + if (l == i) then + E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1) + else + E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1) + end if + + ! Find correlated angular distribution for closest outgoing energy bin + if (r1 - c_j < c_j1 - r1) then + j = j + else + j = j + 1 + end if + + ! Sample outgoing cosine bin + k = 1 + int(prn() * sab % n_inelastic_mu) + + ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) + mu_ijk = sab % inelastic_data(i) % mu(k, j) + mu_i1jk = sab % inelastic_data(i) % mu(k, j) + + ! Cosine of angle between incoming and outgoing neutron + mu = (ONE - f) * mu_ijk + f * mu_i1jk + else message = "Invalid secondary energy mode on S(a,b) table " // & trim(sab % name) - end if - - ! Determine outgoing energy corresponding to E_in(i) and E_in(i+1) - E_ij = sab % inelastic_e_out(j,i) - E_i1j = sab % inelastic_e_out(j,i+1) - - ! Outgoing energy - E = (1 - f)*E_ij + f*E_i1j - - ! Sample outgoing cosine bin - k = 1 + int(prn() * sab % n_inelastic_mu) - - ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) - mu_ijk = sab % inelastic_mu(k,j,i) - mu_i1jk = sab % inelastic_mu(k,j,i+1) - - ! Cosine of angle between incoming and outgoing neutron - mu = (1 - f)*mu_ijk + f*mu_i1jk - end if + end if ! (inelastic secondary energy treatment) + end if ! (elastic or inelastic) ! change direction of particle uvw = rotate_angle(uvw, mu) From 56d3d382826d9a3df29e60f51d6e5abcdf211f98 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 15 Dec 2013 15:58:55 -0500 Subject: [PATCH 15/31] Modified code after some testing to not perform interpolation on Ein and Eout when finding mu; instead the closest distribution is chosen, as is done for file 4 angular distributions in sample_angle. --- src/physics.F90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 372d67bd3e..2d3cfeb39f 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -698,7 +698,7 @@ contains E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1) end if - ! Find correlated angular distribution for closest outgoing energy bin + ! Find angular distribution for closest outgoing energy bin if (r1 - c_j < c_j1 - r1) then j = j else @@ -708,12 +708,9 @@ contains ! Sample outgoing cosine bin k = 1 + int(prn() * sab % n_inelastic_mu) - ! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1) - mu_ijk = sab % inelastic_data(i) % mu(k, j) - mu_i1jk = sab % inelastic_data(i) % mu(k, j) - - ! Cosine of angle between incoming and outgoing neutron - mu = (ONE - f) * mu_ijk + f * mu_i1jk + ! Will use mu from the randomly chosen incoming and closest outgoing + ! energy bins + mu = sab % inelastic_data(l) % mu(k, j) else message = "Invalid secondary energy mode on S(a,b) table " // & From b64644407e5c04aaefe865ceb3a43be288e2b72c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 17 Dec 2013 21:19:42 -0500 Subject: [PATCH 16/31] Two small fixes in installation instructions. --- docs/source/usersguide/install.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index c05782e71d..164b0151be 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -78,7 +78,7 @@ Prerequisites ./configure --prefix=/opt/hdf5/1.8.11-gnu --enable-fortran \ --enable-fortran2003 --enable-parallel - You may omit '--enable-parallel' if you want to compile HDF5_ in serial. + You may omit ``--enable-parallel`` if you want to compile HDF5_ in serial. * PETSc_ for CMFD acceleration @@ -93,7 +93,7 @@ Prerequisites --with-fortran-datatypes The BLAS/LAPACK library is not required to be downloaded and can be linked - explicitly (e.g., Intel MLK library). + explicitly (e.g., Intel MKL library). * git_ version control software for obtaining source code From 96999e00c6eb2ab1172c590268bc79155fea9d84 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Dec 2013 11:18:58 -0500 Subject: [PATCH 17/31] Make sure uvw is initialized when sampling source position. --- src/source.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/source.F90 b/src/source.F90 index 4fcdd0074a..231a23cc5c 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -97,6 +97,7 @@ contains ! Fill p with needed data p % coord0 % xyz = site % xyz + p % coord0 % uvw = [ ONE, ZERO, ZERO ] ! Now search to see if location exists in geometry call find_cell(p, found) From d28750c34fdec59b3ce955865299eafe980cba16 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Dec 2013 11:48:47 -0500 Subject: [PATCH 18/31] Fix bug in convert_xsdir.py. --- src/utils/convert_xsdir.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/convert_xsdir.py b/src/utils/convert_xsdir.py index fa29f46eca..32315f1ed5 100755 --- a/src/utils/convert_xsdir.py +++ b/src/utils/convert_xsdir.py @@ -71,7 +71,7 @@ class Xsdir(object): # Handle continuation lines while words[-1] == '+': extraWords = self.f.readline().split() - words = words + extraWords + words = words[:-1] + extraWords assert len(words) >= 7 # Create XsdirTable object and add to line From d0028d388dc68b61efa520ac8ddc937efa8eb6a3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Dec 2013 17:10:31 -0500 Subject: [PATCH 19/31] Allow mixed ascii/binary cross_sections.xml. --- src/input_xml.F90 | 14 +++++++++++++- src/relaxng/cross_sections.rnc | 4 +++- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8cd93f0e10..63fb2163e3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2938,7 +2938,19 @@ contains end if ! set filetype, record length, and number of entries - listing % filetype = filetype + if (check_for_node(node_ace, "filetype")) then + temp_str = '' + call get_node_value(node_ace, "filetype", temp_str) + if (temp_str == 'ascii') then + listing % filetype = ASCII + else if (temp_str == 'binary') then + listing % filetype = BINARY + end if + else + listing % filetype = filetype + end if + + ! Set record length and entries for binary files if (filetype == BINARY) then listing % recl = recl listing % entries = entries diff --git a/src/relaxng/cross_sections.rnc b/src/relaxng/cross_sections.rnc index e82bb986fb..75c0ea29c2 100644 --- a/src/relaxng/cross_sections.rnc +++ b/src/relaxng/cross_sections.rnc @@ -10,7 +10,9 @@ element cross_sections { (element temperature { xsd:double } | attribute temperature { xsd:double }) & (element path { xsd:string { maxLength = "255" } } | attribute path { xsd:string { maxLength = "255" } }) & - (element location { xsd:int } | attribute location { xsd:int })? + (element location { xsd:int } | attribute location { xsd:int })? & + (element filetype { ( "ascii" | "binary" ) } | + attribute filetype { ( "ascii" | "binary" ) })? }* & element directory { xsd:string { maxLength = "255" } }? & From 4a8157f28081eb9b470e408bda3db1cb65b71dff Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 20 Dec 2013 07:10:38 -0500 Subject: [PATCH 20/31] Included explicit integer casting in the reading of thermal data. --- src/ace.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index ff8582470d..aa6de5ab64 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -1195,7 +1195,7 @@ contains integer :: NE_out ! number of outgoing energies integer :: NMU ! number of outgoing angles integer :: JXS4 ! location of elastic energy table - real(8), allocatable :: LOCC(:) ! Location of inelastic data + integer(8), allocatable :: LOCC(:) ! Location of inelastic data ! read secondary energy mode for inelastic scattering table % secondary_mode = NXS(7) @@ -1252,10 +1252,10 @@ contains end do else if (table % secondary_mode == SAB_SECONDARY_CONT) then ! Get the location pointers to each Ein's DistEnergySAB data - LOCC = get_real(NE_in) + LOCC = get_int(NE_in) ! Get the number of outgoing energies and allocate space accordingly do i = 1, NE_in - NE_out = XSS(XSS_index + i - 1) + NE_out = int(XSS(XSS_index + i - 1)) table % inelastic_data(i) % n_e_out = NE_out allocate(table % inelastic_data(i) % e_out (NE_out)) allocate(table % inelastic_data(i) % e_out_pdf (NE_out)) From 4bda7f2054fe8174a6eacf72d367499307227a89 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 20 Dec 2013 09:44:09 -0500 Subject: [PATCH 21/31] Updated theory manual to incorporate the new thermal scattering sampling technique --- docs/source/methods/physics.rst | 64 ++++++++++++++++++++++++++------- 1 file changed, 51 insertions(+), 13 deletions(-) diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index c79144dfe3..fc82dc1a98 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -790,6 +790,7 @@ outgoing angle is \mu = \frac{1}{A} \ln \left ( \xi_4 e^A + (1 - \xi_4) e^{-A} \right ). +.. _ace-law-61: ACE Law 61 - Correlated Energy and Angle Distribution +++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -952,7 +953,7 @@ as v_n \bar{\sigma} (v_n, T) = \int d\mathbf{v}_T v_r \sigma(v_r) M (\mathbf{v}_T) - + where :math:`v_n` is the magnitude of the velocity of the neutron, :math:`\bar{\sigma}` is an effective cross section, :math:`T` is the temperature of the target material, :math:`\mathbf{v}_T` is the velocity of the target @@ -1321,7 +1322,7 @@ given analytically by \mu = 1 - \frac{E_i}{E} -where :math:`E_i` is the energy of the Bragg edge that scattered the neutron. +where :math:`E_i` is the energy of the Bragg edge that scattered the neutron. Outgoing Angle for Incoherent Elastic Scattering ------------------------------------------------ @@ -1348,18 +1349,24 @@ where the interpolation factor is defined as Outgoing Energy and Angle for Inelastic Scattering -------------------------------------------------- -On each |sab| table, there is a correlated angle-energy secondary distribution -for neutron thermal inelastic scattering. While the documentation for the ACE -format implies that there are a series of equiprobable outgoing energies, the -outgoing energies may have non-uniform probability distribution. In particular, -if the thermal data were processed with :math:`iwt = 0` in NJOY, then the first -and last outgoing energies have a relative probability of 1, the second and -second to last energies have a relative probability of 4, and all other energies -have a relative probability of 10. The procedure to determine the outgoing +Each |sab| table provides a correlated angle-energy secondary distribution for +neutron thermal inelastic scattering. There are three representations used +in the ACE thermal scattering data: equiprobable discrete outgoing +energies, non-uniform yet still discrete outgoing energies, and continuous +outgoing energies with corresponding probability and cumulative distribution +functions provided in tabular format. These three representations all +represent the angular distribution in a common format, using a series of +discrete equiprobable outgoing cosines. + +Equi-Probable Outgoing Energies ++++++++++++++++++++++++++++++++ + +If the thermal data was processed with :math:`iwt = 1` in NJOY, then the +outgoing energy spectra is represented in the ACE data as a set of discrete and +equiprobable outgoing energies. The procedure to determine the outgoing energy and angle is as such. First, the interpolation factor is determined from -equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is sampled -either from a uniform distribution or from the aforementioned skewed -distribution. The outgoing energy is then interpolated between values +equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is +sampled from ths uniform distribution and then interpolated between values corresponding to neighboring incoming energies: .. math:: @@ -1380,6 +1387,37 @@ uniformly and then the final cosine is interpolated on the incoming energy grid: where :math:`\mu_{i,j,k}` is the k-th outgoing cosine corresponding to the j-th outgoing energy and the i-th incoming energy. +Skewed Equi-Probable Outgoing Energies +++++++++++++++++++++++++++++++++++++++ + +If the thermal data was processed with :math:`iwt=0` in NJOY, then the +outgoing energy spectra is represented in the ACE data according to the +following: the first and last outgoing energies have a relative probability of +1, the second and second-to-last energies have a relative probability of 4, and +all other energies have a relative probability of 10. The procedure to +determine the outgoing energy and angle is similar to the method discussed +above, except that the sampled probability distribution is now skewed +accordingly. + +Continuous Outgoing Energies +++++++++++++++++++++++++++++ + +If the thermal data was processed with :math:`iwt=2` in NJOY, then the +outgoing energy spectra is represented by a continuous outgoing energy spectra +in tabular form with linear-linear interpolation. The sampling of the outgoing +energy portion of this format is very similar to :ref:`ACE Law 61`, +but the sampling of the correlated angle is performed as it was in the other +two representations discussed in this sub-section. In the Law 61 algorithm, +we found an interpolation factor :math:`f`, statistically sampled an incoming +energy bin :math:`\ell`, and sampled an outgoing energy bin :math:`j` based on +the tabulated cumulative distribution function. Once the outgoing energy has +been determined with equation :eq:`ace-law-4-energy`, we then need to decide +which angular distribution data to use. Like the linear-linear interpolation +case in Law 61, the angular distribution closest to the sampled value of the +cumulative distribution function for the outgoing energy is utilized. The +actual algorithm utilized to sample the outgoing angle is shown in equation +:eq:`inelastic-angle`. + .. _probability_tables: ---------------------------------------------- From 402df187428ebf1f985b6ac73cd2252ff0ca434d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 22 Dec 2013 20:41:49 -0500 Subject: [PATCH 22/31] Fixed a typo in theory documentation. --- docs/source/methods/physics.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index fc82dc1a98..54dc913455 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -1363,10 +1363,10 @@ Equi-Probable Outgoing Energies If the thermal data was processed with :math:`iwt = 1` in NJOY, then the outgoing energy spectra is represented in the ACE data as a set of discrete and -equiprobable outgoing energies. The procedure to determine the outgoing -energy and angle is as such. First, the interpolation factor is determined from +equiprobable outgoing energies. The procedure to determine the outgoing energy +and angle is as such. First, the interpolation factor is determined from equation :eq:`sab-interpolation-factor`. Then, an outgoing energy bin is -sampled from ths uniform distribution and then interpolated between values +sampled from a uniform distribution and then interpolated between values corresponding to neighboring incoming energies: .. math:: From 5db982406428e40e09154b530b0de7ee7b7ed00b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 22 Dec 2013 20:57:36 -0500 Subject: [PATCH 23/31] Fix particle track bug when using Intel compiler. --- src/input_xml.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 63fb2163e3..3ce3c9d0d4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -474,7 +474,8 @@ contains allocate(temp_int_array(n_tracks)) call get_node_array(doc, "track", temp_int_array) - ! Reshape into track_identifiers -- note automatic array allocation + ! Reshape into track_identifiers + allocate(track_identifiers(3, n_tracks/3)) track_identifiers = reshape(temp_int_array, [3, n_tracks/3]) end if From 73c542e2016e7dc7290973671bf2151e9518d88c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 26 Dec 2013 22:44:03 -0500 Subject: [PATCH 24/31] Display number of OpenMP threads under title banner. --- src/output.F90 | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/src/output.F90 b/src/output.F90 index 1507aa4cd5..677c35396c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -31,6 +31,10 @@ contains subroutine title() +#ifdef OPENMP + use omp_lib +#endif + write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & ' .d88888b. 888b d888 .d8888b.', & ' d88P" "Y88b 8888b d8888 d88P Y88b', & @@ -46,25 +50,31 @@ contains ! Write version information write(UNIT=OUTPUT_UNIT, FMT=*) & - ' Copyright: 2011-2013 Massachusetts Institute of Technology' + ' Copyright: 2011-2013 Massachusetts Institute of Technology' write(UNIT=OUTPUT_UNIT, FMT=*) & - ' License: http://mit-crpg.github.io/openmc/license.html' - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & + ' License: http://mit-crpg.github.io/openmc/license.html' + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",8X,I1,".",I1,".",I1)') & VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE #ifdef GIT_SHA1 - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",6X,A)') GIT_SHA1 + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Git SHA1:",7X,A)') GIT_SHA1 #endif ! Write the date and time - write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A)') & + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",6X,A)') & time_stamp() #ifdef MPI ! Write number of processors - write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",1X,A)') & + write(UNIT=OUTPUT_UNIT, FMT='(6X,"MPI Processes:",2X,A)') & trim(to_str(n_procs)) #endif +#ifdef OPENMP + ! Write number of OpenMP threads + write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') & + trim(to_str(omp_get_max_threads())) +#endif + end subroutine title !=============================================================================== From f3f024bbc60f499b5ca38f4e4f156f1f2df1a95f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 26 Dec 2013 23:17:03 -0500 Subject: [PATCH 25/31] Avoid passing temporary arrays to write_data. --- src/hdf5_summary.F90 | 5 +++-- src/state_point.F90 | 6 +++--- src/track_output.F90 | 4 +++- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index 915e6af070..c9e367200c 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -100,6 +100,7 @@ contains integer :: i, j, k, m integer :: n_x, n_y, n_z + integer :: length(3) integer, allocatable :: lattice_universes(:,:,:) type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() @@ -312,8 +313,8 @@ contains end do end do end do - call su % write_data(lattice_universes, "universes", & - length=(/n_x, n_y, n_z/), & + length = [n_x, n_y, n_z] + call su % write_data(lattice_universes, "universes", length=length, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) deallocate(lattice_universes) diff --git a/src/state_point.F90 b/src/state_point.F90 index 3793923bb2..e13cd8f486 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -464,6 +464,7 @@ contains character(19) :: current_time integer :: i integer :: j + integer :: length(4) integer :: int_array(3) integer, allocatable :: temp_array(:) real(8) :: real_array(3) @@ -542,10 +543,9 @@ contains call sp % read_data(cmfd % indices, "indicies", length=4, group="cmfd") call sp % read_data(cmfd % k_cmfd, "k_cmfd", length=restart_batch, & group="cmfd") + length = cmfd % indices([4,1,2,3]) call sp % read_data(cmfd % cmfd_src, "cmfd_src", & - length=(/cmfd % indices(4), cmfd % indices(1), & - cmfd % indices(2), cmfd % indices(3)/), & - group="cmfd") + length=length, group="cmfd") call sp % read_data(cmfd % entropy, "cmfd_entropy", & length=restart_batch, group="cmfd") call sp % read_data(cmfd % balance, "cmfd_balance", & diff --git a/src/track_output.F90 b/src/track_output.F90 index 021cbb8ca6..79a2809066 100644 --- a/src/track_output.F90 +++ b/src/track_output.F90 @@ -56,6 +56,7 @@ contains subroutine finalize_particle_track(p) type(Particle), intent(in) :: p + integer :: length(2) character(MAX_FILE_LEN) :: fname type(BinaryOutput) :: binout @@ -69,7 +70,8 @@ contains // '.binary' #endif call binout % file_create(fname) - call binout % write_data(coords, 'coordinates', length=(/3, n_tracks/)) + length = [3, n_tracks] + call binout % write_data(coords, 'coordinates', length=length) call binout % file_close() deallocate(coords) end subroutine finalize_particle_track From 9a0e50574cfb2c8b7220c4ede790e65d61bcf776 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 27 Dec 2013 16:33:25 -0500 Subject: [PATCH 26/31] Fix bug causing CMFD to report incorrect results when using threads. --- src/cmfd_input.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 120db5f649..6f259f9497 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -577,7 +577,9 @@ contains end do ! Put cmfd tallies into active tally array and turn tallies on +!$omp parallel call setup_active_cmfdtallies() +!$omp end parallel tallies_on = .true. end subroutine create_cmfd_tally From 415e5da764cdc6d05dc43402b958c9e74d359048 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Jan 2014 15:44:48 -0500 Subject: [PATCH 27/31] Added description of UFS method in documentation. --- docs/source/methods/eigenvalue.rst | 38 ++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/docs/source/methods/eigenvalue.rst b/docs/source/methods/eigenvalue.rst index a5ba0bf459..fe99ba22ec 100644 --- a/docs/source/methods/eigenvalue.rst +++ b/docs/source/methods/eigenvalue.rst @@ -108,6 +108,40 @@ at plots of :math:`k_{eff}` and the Shannon entropy. A number of methods have been proposed (see e.g. [Romano]_, [Ueki]_), but each of these is not without problems. +--------------------------- +Uniform Fission Site Method +--------------------------- + +Generally speaking, the variance of a Monte Carlo tally will be inversely +proportional to the number of events that score to the tally. In a reactor +problem, this implies that regions with low relative power density will have +higher variance that regions with high relative power density. One method to +circumvent the uneven distribution of relative errors is the uniform fission +site (UFS) method introduced by [Sutton]_. In this method, the portion of the +problem containing fissionable material is subdivided into a number of cells +(typically using a structured mesh). Rather than producing + +.. math:: + + m = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t} + +fission sites at each collision where :math:`w` is the weight of the neutron, +:math:`k` is the previous-generation estimate of the neutron multiplication +factor, :math:`\nu\Sigma_f` is the neutron production cross section, and +:math:`\Sigma_t` is the total cross section, in the UFS method we produce + +.. math:: + + m_{UFS} = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t} \frac{v_i}{s_i} + +fission sites at each collision where :math:`v_i` is the fraction of the total +volume occupied by cell :math:`i` and :math:`s_i` is the fraction of the fission +source contained in cell :math:`i`. To ensure that no bias is introduced, the +weight of each fission site stored in the fission bank is :math:`s_i/v_i` rather +than unity. By ensuring that the expected number of fission sites in each mesh +cell is constant, the collision density across all cells, and hence the variance +of tallies, is more uniform than it would be otherwise. + .. _Shannon entropy: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-06-3737_entropy.pdf .. [Lieberoth] J. Lieberoth, "A Monte Carlo Technique to Solve the Static @@ -119,5 +153,9 @@ problems. *Proc. International Conference on Mathematics, Computational Methods, and Reactor Physics*, Saratoga Springs, New York (2009). +.. [Sutton] Daniel J. Kelly, Thomas M. Sutton, and Stephen C. Wilson, "MC21 + Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark + Problem," *Proc. PHYSOR 2012*, Knoxville, Tennessee, Apr. 15--20 (2012). + .. [Ueki] Taro Ueki, "On-the-Fly Judgments of Monte Carlo Fission Source Convergence," *Trans. Am. Nucl. Soc.*, **98**, 512 (2008). From 0f628689132f503060c06f1419d7c3c487d0f9a3 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 3 Jan 2014 14:20:54 -0500 Subject: [PATCH 28/31] Made some revisions to plot_mesh_tally and statepoint for usability purposes. --- src/utils/plot_mesh_tally.py | 340 ++++++++++++++++++++--------------- src/utils/statepoint.py | 33 ++-- 2 files changed, 208 insertions(+), 165 deletions(-) diff --git a/src/utils/plot_mesh_tally.py b/src/utils/plot_mesh_tally.py index 04b6b4592a..c70b291653 100755 --- a/src/utils/plot_mesh_tally.py +++ b/src/utils/plot_mesh_tally.py @@ -16,123 +16,161 @@ from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as Naviga import numpy as np class AppForm(QMainWindow): - def __init__(self, parent=None): + def __init__(self, argv, parent=None): QMainWindow.__init__(self, parent) # Read data from source or leakage fraction file - self.get_file_data() - self.main_frame = QWidget() - self.setCentralWidget(self.main_frame) - - # Create the Figure, Canvas, and Axes - self.dpi = 100 - self.fig = Figure((5.0, 15.0), dpi=self.dpi) - self.canvas = FigureCanvas(self.fig) - self.canvas.setParent(self.main_frame) - self.axes = self.fig.add_subplot(111) - - # Create the navigation toolbar, tied to the canvas - self.mpl_toolbar = NavigationToolbar(self.canvas, self.main_frame) + self.all_good = False + while not self.all_good: + if len(argv) > 1: + cl_file = str(argv[1]) + else: + cl_file = None + self.get_file_data(cl_file) + # Check that there are any mesh tallies at all + if len(self.tally_ids) != 0: + self.all_good = True + else: + # if there are not, the user will be given the choice to choose + # another file (but only if using interactive chooser) + if cl_file is None: + choice = QMessageBox.critical(None, "Invalid StatePoint File", + "File Does Not Contain Mesh " + + "Tallies!" + + "\nSelect Another File Or Quit", + QMessageBox.Retry, + QMessageBox.Abort) + if choice == QMessageBox.Abort: + self.all_good = False + break + else: + print("Invalid StatePoint File; File Does Not Contain " + + "Mesh Tallies!") + self.all_good = False + break - # Grid layout at bottom - self.grid = QGridLayout() - # Overall layout - self.vbox = QVBoxLayout() - self.vbox.addWidget(self.canvas) - self.vbox.addWidget(self.mpl_toolbar) - self.vbox.addLayout(self.grid) - self.main_frame.setLayout(self.vbox) + if self.all_good: + # Set maximum colorbar value by maximum tally data value + self.maxvalue = self.datafile.tallies[0].results.max() - # Tally selections - label_tally = QLabel("Tally:") - self.tally = QComboBox() - self.tally.addItems([(str(i + 1)) for i in range(self.n_tallies)]) - self.connect(self.tally, SIGNAL('activated(int)'), - self._update) - self.connect(self.tally, SIGNAL('activated(int)'), - self.populate_boxes) - self.connect(self.tally, SIGNAL('activated(int)'), - self.on_draw) + self.main_frame = QWidget() + self.setCentralWidget(self.main_frame) - # Planar basis - label_basis = QLabel("Basis:") - self.basis = QComboBox() - self.basis.addItems(['xy', 'yz', 'xz']) + # Create the Figure, Canvas, and Axes + self.dpi = 100 + self.fig = Figure((5.0, 15.0), dpi=self.dpi) + self.canvas = FigureCanvas(self.fig) + self.canvas.setParent(self.main_frame) + self.axes = self.fig.add_subplot(111) - # Update window when 'Basis' selection is changed - self.connect(self.basis, SIGNAL('activated(int)'), - self._update) - self.connect(self.basis, SIGNAL('activated(int)'), - self.populate_boxes) - self.connect(self.basis, SIGNAL('activated(int)'), - self.on_draw) + # Create the navigation toolbar, tied to the canvas + self.mpl_toolbar = NavigationToolbar(self.canvas, self.main_frame) - # Axial level within selected basis - label_axial_level = QLabel("Axial Level:") - self.axial_level = QComboBox() - self.connect(self.axial_level, SIGNAL('activated(int)'), - self.on_draw) - - # Add Option to plot mean or uncertainty - label_mean = QLabel("Mean or Uncertainty:") - self.mean = QComboBox() - self.mean.addItems(['Mean','Absolute Uncertainty', - 'Relative Uncertainty']) - - # Update window when mean selection is changed - self.connect(self.mean, SIGNAL('activated(int)'), - self.on_draw) - + # Grid layout at bottom + self.grid = QGridLayout() - self.label_filters = QLabel("Filter options:") + # Overall layout + self.vbox = QVBoxLayout() + self.vbox.addWidget(self.canvas) + self.vbox.addWidget(self.mpl_toolbar) + self.vbox.addLayout(self.grid) + self.main_frame.setLayout(self.vbox) - # Labels for all possible filters - self.labels = {'cell': 'Cell: ', 'cellborn': 'Cell born: ', - 'surface': 'Surface: ', 'material': 'Material', - 'universe': 'Universe: ', 'energyin': 'Energy in: ', - 'energyout': 'Energy out: '} - - # Empty reusable labels - self.qlabels = {} - for j in range(8): - self.nextLabel = QLabel - self.qlabels[j] = self.nextLabel - - # Reusable comboboxes labelled with filter names - self.boxes = {} - for key in self.labels.keys(): - self.nextBox = QComboBox() - self.connect(self.nextBox, SIGNAL('activated(int)'), + # Tally selections + label_tally = QLabel("Tally:") + self.tally = QComboBox() + # Only show options for the tallies with meshes + self.tally.addItems([str(i + 1) for i in self.tally_ids]) + self.connect(self.tally, SIGNAL('activated(int)'), + self._update) + self.connect(self.tally, SIGNAL('activated(int)'), + self.populate_boxes) + self.connect(self.tally, SIGNAL('activated(int)'), self.on_draw) - self.boxes[key] = self.nextBox - # Combobox to select among scores - self.score_label = QLabel("Score:") - self.scoreBox = QComboBox() - for item in self.tally_scores[0]: - self.scoreBox.addItems(str(item)) - self.connect(self.scoreBox, SIGNAL('activated(int)'), - self.on_draw) + # Planar basis + label_basis = QLabel("Basis:") + self.basis = QComboBox() + self.basis.addItems(['xy', 'yz', 'xz']) - # Fill layout - self.grid.addWidget(label_tally, 0, 0) - self.grid.addWidget(self.tally, 0, 1) - self.grid.addWidget(label_basis, 1, 0) - self.grid.addWidget(self.basis, 1, 1) - self.grid.addWidget(label_axial_level, 2, 0) - self.grid.addWidget(self.axial_level, 2, 1) - self.grid.addWidget(label_mean, 3, 0) - self.grid.addWidget(self.mean, 3, 1) - self.grid.addWidget(self.label_filters, 4, 0) + # Update window when 'Basis' selection is changed + self.connect(self.basis, SIGNAL('activated(int)'), + self._update) + self.connect(self.basis, SIGNAL('activated(int)'), + self.populate_boxes) + self.connect(self.basis, SIGNAL('activated(int)'), + self.on_draw) - self._update() - self.populate_boxes() - self.on_draw() + # Axial level within selected basis + label_axial_level = QLabel("Axial Level:") + self.axial_level = QComboBox() + self.connect(self.axial_level, SIGNAL('activated(int)'), + self.on_draw) - def get_file_data(self): + # Add Option to plot mean or uncertainty + label_mean = QLabel("Mean or Uncertainty:") + self.mean = QComboBox() + self.mean.addItems(['Mean','Absolute Uncertainty', + 'Relative Uncertainty']) + + # Update window when mean selection is changed + self.connect(self.mean, SIGNAL('activated(int)'), + self.on_draw) + + + self.label_filters = QLabel("Filter options:") + + # Labels for all possible filters + self.labels = {'cell': 'Cell: ', 'cellborn': 'Cell born: ', + 'surface': 'Surface: ', 'material': 'Material', + 'universe': 'Universe: ', 'energyin': 'Energy in: ', + 'energyout': 'Energy out: '} + + # Empty reusable labels + self.qlabels = {} + for j in range(8): + self.nextLabel = QLabel + self.qlabels[j] = self.nextLabel + + # Reusable comboboxes labelled with filter names + self.boxes = {} + for key in self.labels.keys(): + self.nextBox = QComboBox() + self.connect(self.nextBox, SIGNAL('activated(int)'), + self.on_draw) + self.boxes[key] = self.nextBox + + # Combobox to select among scores + self.score_label = QLabel("Score:") + self.scoreBox = QComboBox() + for item in self.tally_scores[0]: + self.scoreBox.addItems(str(item)) + self.connect(self.scoreBox, SIGNAL('activated(int)'), + self.on_draw) + + # Fill layout + self.grid.addWidget(label_tally, 0, 0) + self.grid.addWidget(self.tally, 0, 1) + self.grid.addWidget(label_basis, 1, 0) + self.grid.addWidget(self.basis, 1, 1) + self.grid.addWidget(label_axial_level, 2, 0) + self.grid.addWidget(self.axial_level, 2, 1) + self.grid.addWidget(label_mean, 3, 0) + self.grid.addWidget(self.mean, 3, 1) + self.grid.addWidget(self.label_filters, 4, 0) + + self._update() + self.populate_boxes() + self.on_draw() + + def get_file_data(self, cl_file=None): # Get data file name from "open file" browser - filename = QFileDialog.getOpenFileName(self, 'Select statepoint file', '.') + if cl_file is None: + filename = QFileDialog.getOpenFileName(self, + 'Select statepoint file', '.') + else: + filename = cl_file # Create StatePoint object and read in data self.datafile = StatePoint(str(filename)) @@ -141,32 +179,31 @@ class AppForm(QMainWindow): self.setWindowTitle('Core Map Tool : ' + str(self.datafile.path)) - # Set maximum colorbar value by maximum tally data value - self.maxvalue = self.datafile.tallies[0].results.max() - self.labelList = [] # Read mesh dimensions -# for mesh in self.datafile.meshes: -# self.nx, self.ny, self.nz = mesh.dimension + # for mesh in self.datafile.meshes: + # self.nx, self.ny, self.nz = mesh.dimension - # Read filter types from statepoint file + # Find which tallies have meshes so the rest can be ignored, + # and for these tallies read the filter and score types + self.tally_ids = [] self.n_tallies = len(self.datafile.tallies) self.tally_list = [] - for tally in self.datafile.tallies: - self.filter_types = [] - for f in tally.filters: - self.filter_types.append(f) - self.tally_list.append(self.filter_types) - - # Read score types from statepoint file self.tally_scores = [] - for tally in self.datafile.tallies: - self.score_types = [] - for s in tally.scores: - self.score_types.append(s) - self.tally_scores.append(self.score_types) -# print 'self.tally_scores = ', self.tally_scores + for itally, tally in enumerate(self.datafile.tallies): + if 'mesh' in tally.filters: + # Then we have a good tally, store the ID, filters and + # scores + self.tally_ids.append(itally) + self.filter_types = [] + for f in tally.filters: + self.filter_types.append(f) + self.tally_list.append(self.filter_types) + self.score_types = [] + for s in tally.scores: + self.score_types.append(s) + self.tally_scores.append(self.score_types) def on_draw(self): """ Redraws the figure @@ -178,70 +215,73 @@ class AppForm(QMainWindow): axial_level = self.axial_level.currentIndex() + 1 is_mean = self.mean.currentIndex() + # get current tally index + tally_id = self.tally_ids[self.tally.currentIndex()] + # Create spec_list spec_list = [] - for tally in self.datafile.tallies[self.tally.currentIndex()].filters.values(): + for tally in self.datafile.tallies[tally_id].filters.values(): if tally.type == 'mesh': continue index = self.boxes[tally.type].currentIndex() spec_list.append((tally.type, index)) - + # Take is_mean and convert it to an index of the score score_loc = is_mean if score_loc > 1: score_loc = 1 - + if self.basis.currentText() == 'xy': matrix = np.zeros((self.nx, self.ny)) for i in range(self.nx): for j in range(self.ny): - matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, j, axial_level))], + matrix[i,j] = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, j + 1, axial_level))], self.scoreBox.currentIndex())[score_loc] # Calculate relative uncertainty from absolute, if # requested if is_mean == 2: # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, j, axial_level))], + mean_val = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, j + 1, axial_level))], self.scoreBox.currentIndex())[0] if mean_val > 0.0: matrix[i,j] = matrix[i,j] / mean_val else: matrix[i,j] = 0.0 - + elif self.basis.currentText() == 'yz': matrix = np.zeros((self.ny, self.nz)) for i in range(self.ny): for j in range(self.nz): - matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (axial_level, i, j))], + matrix[i,j] = self.datafile.get_value(tally_id, + spec_list + [('mesh', (axial_level, i + 1, j + 1))], self.scoreBox.currentIndex())[score_loc] # Calculate relative uncertainty from absolute, if # requested if is_mean == 2: # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (axial_level, i, j))], + mean_val = self.datafile.get_value(tally_id, + spec_list + [('mesh', (axial_level, i + 1, j + 1))], self.scoreBox.currentIndex())[0] if mean_val > 0.0: matrix[i,j] = matrix[i,j] / mean_val else: matrix[i,j] = 0.0 - + else: matrix = np.zeros((self.nx, self.nz)) for i in range(self.nx): for j in range(self.nz): - matrix[i,j] = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, axial_level, j))], + matrix[i,j] = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, axial_level, j + 1))], self.scoreBox.currentIndex())[score_loc] # Calculate relative uncertainty from absolute, if # requested if is_mean == 2: # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(self.tally.currentIndex(), - spec_list + [('mesh', (i, axial_level, j))], + mean_val = self.datafile.get_value(tally_id, + spec_list + [('mesh', (i + 1, axial_level, j + 1))], self.scoreBox.currentIndex())[0] if mean_val > 0.0: matrix[i,j] = matrix[i,j] / mean_val @@ -255,8 +295,8 @@ class AppForm(QMainWindow): # Make figure, set up color bar self.axes = self.fig.add_subplot(111) - cax = self.axes.imshow(matrix, vmin=0.0, vmax=matrix.max(), - interpolation="nearest") + cax = self.axes.imshow(matrix.transpose(), vmin=0.0, vmax=matrix.max(), + interpolation="nearest", origin='lower') self.fig.colorbar(cax) self.axes.set_xticks([]) @@ -271,9 +311,11 @@ class AppForm(QMainWindow): ''' # print 'Calling _update...' + # get current tally index + tally_id = self.tally_ids[self.tally.currentIndex()] + self.mesh = self.datafile.meshes[ - self.datafile.tallies[ - self.tally.currentIndex()].filters['mesh'].bins[0] - 1] + self.datafile.tallies[tally_id].filters['mesh'].bins[0] - 1] self.nx, self.ny, self.nz = self.mesh.dimension @@ -289,8 +331,7 @@ class AppForm(QMainWindow): self.axial_level.addItems([str(i+1) for i in range(self.ny)]) # Determine maximum value from current tally data set - self.maxvalue = self.datafile.tallies[ - self.tally.currentIndex()].results.max() + self.maxvalue = self.datafile.tallies[tally_id].results.max() # print self.maxvalue # Clear and hide old filter labels @@ -307,6 +348,9 @@ class AppForm(QMainWindow): def populate_boxes(self): # print 'Calling populate_boxes...' + # get current tally index + tally_id = self.tally_ids[self.tally.currentIndex()] + n = 5 labels = {'cell': 'Cell : ', 'cellborn': 'Cell born: ', @@ -317,8 +361,7 @@ class AppForm(QMainWindow): # For each filter in newly-selected tally, name a label and fill the # relevant combobox with options for element in self.tally_list[self.tally.currentIndex()]: - nextFilter = self.datafile.tallies[ - self.tally.currentIndex()].filters[element] + nextFilter = self.datafile.tallies[tally_id].filters[element] if element == 'mesh': continue @@ -337,7 +380,7 @@ class AppForm(QMainWindow): elif element == 'energyin' or element == 'energyout': for i in range(nextFilter.length): - text = (str(nextFilter.bins[i]) + ' to ' + + text = (str(nextFilter.bins[i]) + ' to ' + str(nextFilter.bins[i+1])) combobox.addItem(text) @@ -351,9 +394,10 @@ class AppForm(QMainWindow): def main(): app = QApplication(sys.argv) - form = AppForm() - form.show() - app.exec_() + form = AppForm(app.arguments()) + if form.all_good: + form.show() + app.exec_() if __name__ == "__main__": diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index ec9d9edf78..430abee8b8 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -1,7 +1,6 @@ #!/usr/bin/env python2 import struct -from math import sqrt from collections import OrderedDict import numpy as np @@ -10,10 +9,10 @@ import scipy.stats filter_types = {1: 'universe', 2: 'material', 3: 'cell', 4: 'cellborn', 5: 'surface', 6: 'mesh', 7: 'energyin', 8: 'energyout'} -score_types = {-1: 'flux', +score_types = {-1: 'flux', -2: 'total', -3: 'scatter', - -4: 'nu-scatter', + -4: 'nu-scatter', -5: 'scatter-n', -6: 'scatter-pn', -7: 'transport', @@ -276,7 +275,7 @@ class StatePoint(object): f.bins = self._get_int(path=base+'bins') else: f.bins = self._get_int(f.length, path=base+'bins') - + base = 'tallies/tally' + str(i+1) + '/' # Read nuclide bins @@ -379,7 +378,7 @@ class StatePoint(object): Calculates the sample mean and standard deviation of the mean for each tally bin. """ - + # Determine number of realizations n = self.n_realizations @@ -387,14 +386,14 @@ class StatePoint(object): for i in range(len(self.global_tallies)): # Get sum and sum of squares s, s2 = self.global_tallies[i] - + # Calculate sample mean and replace value s /= n self.global_tallies[i,0] = s # Calculate standard deviation if s != 0.0: - self.global_tallies[i,1] = t_value*sqrt((s2/n - s*s)/(n-1)) + self.global_tallies[i,1] = t_value*np.sqrt((s2/n - s*s)/(n-1)) # Regular tallies for t in self.tallies: @@ -402,14 +401,14 @@ class StatePoint(object): for j in range(t.results.shape[1]): # Get sum and sum of squares s, s2 = t.results[i,j] - + # Calculate sample mean and replace value s /= n t.results[i,j,0] = s # Calculate standard deviation if s != 0.0: - t.results[i,j,1] = t_value*sqrt((s2/n - s*s)/(n-1)) + t.results[i,j,1] = t_value*np.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. @@ -458,7 +457,7 @@ class StatePoint(object): filter_index += value*t.filters[f_type].stride else: filter_index += f_index*t.filters[f_type].stride - + # 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. @@ -531,7 +530,7 @@ class StatePoint(object): for i in range(n_filters): # compute indices for filter combination - filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % + filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1 # append in dictionary bin with filter @@ -544,14 +543,14 @@ class StatePoint(object): dims.reverse() dims = np.asarray(dims) if score_str == 'current': - dims += 1 - meshmax[1:4] = dims + dims += 1 + meshmax[1:4] = dims mesh_bins = np.zeros((n_bins,3)) - mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % + mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1 - mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % + mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1 - mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % + mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1 data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1], mesh_bins[:,2])}) @@ -576,7 +575,7 @@ class StatePoint(object): def _get_data(self, n, typeCode, size): return list(struct.unpack('={0}{1}'.format(n,typeCode), self._f.read(n*size))) - + def _get_int(self, n=1, path=None): if self._hdf5: return [int(v) for v in self._f[path].value] From a3726724ecf9d287707d87d7863c84dbf5f3ef13 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Jan 2014 23:15:59 -0500 Subject: [PATCH 29/31] Ability to generate PDF documentation. --- .gitignore | 1 + docs/Makefile | 20 +++++++++++++++--- docs/{img => source/_images}/3dcore.png | Bin docs/{img => source/_images}/3dgeomplot.png | Bin docs/{img => source/_images}/Tracks.png | Bin docs/{img => source/_images}/atr.png | Bin docs/{img => source/_images}/fluxplot.png | Bin docs/{img => source/_images}/fork.png | Bin docs/{img => source/_images}/halfspace.svg | 0 docs/{img => source/_images}/master-slave.png | Bin .../_images}/nearest-neighbor-example.png | Bin .../_images}/nearest-neighbor.png | Bin docs/{img => source/_images}/openmc.png | Bin .../{img => source/_images}/plotmeshtally.png | Bin docs/{img => source/_images}/pullrequest.png | Bin docs/{img => source/_images}/union.svg | 0 docs/{img => source/_images}/uniongrid.svg | 0 docs/source/conf.py | 4 +++- docs/source/devguide/workflow.rst | 4 ++-- docs/source/methods/cross_sections.rst | 2 +- docs/source/methods/geometry.rst | 4 ++-- docs/source/methods/parallelization.rst | 6 +++--- docs/source/usersguide/processing.rst | 12 +++++------ 23 files changed, 35 insertions(+), 18 deletions(-) rename docs/{img => source/_images}/3dcore.png (100%) rename docs/{img => source/_images}/3dgeomplot.png (100%) rename docs/{img => source/_images}/Tracks.png (100%) rename docs/{img => source/_images}/atr.png (100%) rename docs/{img => source/_images}/fluxplot.png (100%) rename docs/{img => source/_images}/fork.png (100%) rename docs/{img => source/_images}/halfspace.svg (100%) rename docs/{img => source/_images}/master-slave.png (100%) rename docs/{img => source/_images}/nearest-neighbor-example.png (100%) rename docs/{img => source/_images}/nearest-neighbor.png (100%) rename docs/{img => source/_images}/openmc.png (100%) rename docs/{img => source/_images}/plotmeshtally.png (100%) rename docs/{img => source/_images}/pullrequest.png (100%) rename docs/{img => source/_images}/union.svg (100%) rename docs/{img => source/_images}/uniongrid.svg (100%) diff --git a/.gitignore b/.gitignore index 023c4f2c1a..9f720f9894 100644 --- a/.gitignore +++ b/.gitignore @@ -19,6 +19,7 @@ src/openmc # Documentation builds docs/build +docs/source/_images/*.pdf # xml-fortran reader src/xml-fortran/xmlreader diff --git a/docs/Makefile b/docs/Makefile index 8909256f3d..89c71dc074 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -6,13 +6,19 @@ SPHINXOPTS = SPHINXBUILD = sphinx-build PAPER = BUILDDIR = build +IMAGEDIR = source/_images # Internal variables. PAPEROPT_a4 = -D latex_paper_size=a4 PAPEROPT_letter = -D latex_paper_size=letter ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source -.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest +# SVG to PDF conversion +SVG2PDF = inkscape +PDFS = $(patsubst %.svg,%.pdf,$(wildcard $(IMAGEDIR)/*.svg)) + + +.PHONY: help images clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest help: @echo "Please use \`make ' where is one of" @@ -33,8 +39,16 @@ help: @echo " linkcheck to check all external links for integrity" @echo " doctest to run all doctests embedded in the documentation (if enabled)" +# Pattern rule for converting SVG to PDF +%.pdf: %.svg + $(SVG2PDF) -f $< -A $@ + +# Rule to build PDFs +images: $(PDFS) + clean: -rm -rf $(BUILDDIR)/* + -rm $(PDFS) html: $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html @@ -91,14 +105,14 @@ epub: @echo @echo "Build finished. The epub file is in $(BUILDDIR)/epub." -latex: +latex: images $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex @echo @echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex." @echo "Run \`make' in that directory to run these through (pdf)latex" \ "(use \`make latexpdf' here to do that automatically)." -latexpdf: +latexpdf: images $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex @echo "Running LaTeX files through pdflatex..." make -C $(BUILDDIR)/latex all-pdf diff --git a/docs/img/3dcore.png b/docs/source/_images/3dcore.png similarity index 100% rename from docs/img/3dcore.png rename to docs/source/_images/3dcore.png diff --git a/docs/img/3dgeomplot.png b/docs/source/_images/3dgeomplot.png similarity index 100% rename from docs/img/3dgeomplot.png rename to docs/source/_images/3dgeomplot.png diff --git a/docs/img/Tracks.png b/docs/source/_images/Tracks.png similarity index 100% rename from docs/img/Tracks.png rename to docs/source/_images/Tracks.png diff --git a/docs/img/atr.png b/docs/source/_images/atr.png similarity index 100% rename from docs/img/atr.png rename to docs/source/_images/atr.png diff --git a/docs/img/fluxplot.png b/docs/source/_images/fluxplot.png similarity index 100% rename from docs/img/fluxplot.png rename to docs/source/_images/fluxplot.png diff --git a/docs/img/fork.png b/docs/source/_images/fork.png similarity index 100% rename from docs/img/fork.png rename to docs/source/_images/fork.png diff --git a/docs/img/halfspace.svg b/docs/source/_images/halfspace.svg similarity index 100% rename from docs/img/halfspace.svg rename to docs/source/_images/halfspace.svg diff --git a/docs/img/master-slave.png b/docs/source/_images/master-slave.png similarity index 100% rename from docs/img/master-slave.png rename to docs/source/_images/master-slave.png diff --git a/docs/img/nearest-neighbor-example.png b/docs/source/_images/nearest-neighbor-example.png similarity index 100% rename from docs/img/nearest-neighbor-example.png rename to docs/source/_images/nearest-neighbor-example.png diff --git a/docs/img/nearest-neighbor.png b/docs/source/_images/nearest-neighbor.png similarity index 100% rename from docs/img/nearest-neighbor.png rename to docs/source/_images/nearest-neighbor.png diff --git a/docs/img/openmc.png b/docs/source/_images/openmc.png similarity index 100% rename from docs/img/openmc.png rename to docs/source/_images/openmc.png diff --git a/docs/img/plotmeshtally.png b/docs/source/_images/plotmeshtally.png similarity index 100% rename from docs/img/plotmeshtally.png rename to docs/source/_images/plotmeshtally.png diff --git a/docs/img/pullrequest.png b/docs/source/_images/pullrequest.png similarity index 100% rename from docs/img/pullrequest.png rename to docs/source/_images/pullrequest.png diff --git a/docs/img/union.svg b/docs/source/_images/union.svg similarity index 100% rename from docs/img/union.svg rename to docs/source/_images/union.svg diff --git a/docs/img/uniongrid.svg b/docs/source/_images/uniongrid.svg similarity index 100% rename from docs/img/uniongrid.svg rename to docs/source/_images/uniongrid.svg diff --git a/docs/source/conf.py b/docs/source/conf.py index 8752b48986..24c26bcd87 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -121,7 +121,7 @@ html_title = "OpenMC Documentation" # The name of an image file (relative to this directory) to place at the top # of the sidebar. -html_logo = '../img/openmc.png' +html_logo = '_images/openmc.png' # The name of an image file (within the static path) to use as favicon of the # docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 @@ -188,6 +188,8 @@ latex_documents = [ u'Massachusetts Institute of Technology', 'manual'), ] +latex_elements = {'preamble': '\\usepackage{enumitem}\\setlistdepth{9}'} + # The name of an image file (relative to this directory) to place at the top of # the title page. #latex_logo = None diff --git a/docs/source/devguide/workflow.rst b/docs/source/devguide/workflow.rst index 9562c9584e..f3793628af 100644 --- a/docs/source/devguide/workflow.rst +++ b/docs/source/devguide/workflow.rst @@ -60,7 +60,7 @@ features and bug fixes. The general steps for contributing are as follows: repository with the same name under your personal account. As such, you can commit to it as you please without disrupting other developers. - .. image:: ../../img/fork.png + .. image:: ../_images/fork.png 2. Clone your fork of OpenMC and create a branch that branches off of *develop*: @@ -77,7 +77,7 @@ features and bug fixes. The general steps for contributing are as follows: 4. Issue a pull request from GitHub and select the *develop* branch of mit-crpg/openmc as the target. - .. image:: ../../img/pullrequest.png + .. image:: ../_images/pullrequest.png At a minimum, you should describe what the changes you've made are and why you are making them. If the changes are related to an oustanding issue, make diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index c564896c79..a193dde96d 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -80,7 +80,7 @@ dashed box would need to be stored on a per-nuclide basis, and the union grid would need to be stored once. This method is also referred to as *double indexing* and is available as an option in Serpent (see paper by Leppanen_). -.. figure:: ../../img/uniongrid.svg +.. figure:: ../_images/uniongrid.* :width: 600px :align: center :figclass: align-center diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index d4824f5d1c..1515ffa891 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -45,7 +45,7 @@ surface by a combination of the unique ID of the surface and a positive/negative sign. The following illustration shows an example of an ellipse with unique ID 1 dividing space into two half-spaces. -.. figure:: ../../img/halfspace.svg +.. figure:: ../_images/halfspace.* :align: center :figclass: align-center @@ -60,7 +60,7 @@ half-space references whose intersection defines the region. The region is then assigned a material defined elsewhere. The following illustration shows an example of a cell defined as the intersection of an ellipse and two planes. -.. figure:: ../../img/union.svg +.. figure:: ../_images/union.* :align: center :figclass: align-center diff --git a/docs/source/methods/parallelization.rst b/docs/source/methods/parallelization.rst index 7f0b9bd6eb..66ead52a2a 100644 --- a/docs/source/methods/parallelization.rst +++ b/docs/source/methods/parallelization.rst @@ -65,7 +65,7 @@ in the case of an eigenvalue calculation). This idea is illustrated in .. _figure-master-slave: -.. figure:: ../../img/master-slave.png +.. figure:: ../_images/master-slave.png :align: center :figclass: align-center @@ -122,7 +122,7 @@ needed. This concept is illustrated in :ref:`Figure 2 .. _figure-nearest-neighbor: -.. figure:: ../../img/nearest-neighbor.png +.. figure:: ../_images/nearest-neighbor.png :align: center :figclass: align-center @@ -203,7 +203,7 @@ communicated between adjacent nodes. .. _figure-neighbor-example: -.. figure:: ../../img/nearest-neighbor-example.png +.. figure:: ../_images/nearest-neighbor-example.png :align: center :figclass: align-center diff --git a/docs/source/usersguide/processing.rst b/docs/source/usersguide/processing.rst index 043041ef42..fc426ba2f4 100644 --- a/docs/source/usersguide/processing.rst +++ b/docs/source/usersguide/processing.rst @@ -38,7 +38,7 @@ running OpenMC with the -plot or -p command-line option (See Plotting in 2D -------------- -.. image:: ../../img/atr.png +.. image:: ../_images/atr.png :height: 200px After running OpenMC to obtain PPM files, images should be saved to another @@ -58,7 +58,7 @@ Ubuntu: ``sudo apt-get install imagemagick``). Images are then converted like: Plotting in 3D -------------- -.. image:: ../../img/3dgeomplot.png +.. image:: ../_images/3dgeomplot.png :height: 200px The binary VOXEL files output by OpenMC can not be viewed directly by any @@ -162,7 +162,7 @@ tasks will be described here in the following sections. Plotting in 2D -------------- -.. image:: ../../img/plotmeshtally.png +.. image:: ../_images/plotmeshtally.png :height: 200px For simple viewing of 2D slices of a mesh plot, the utility plot_mesh_tally.py @@ -170,7 +170,7 @@ is provided. This utility provides an interactive GUI to explore and plot mesh tallies for any scores and filter bins. It requires statepoint.py, as well as `PyQt `_. -.. image:: ../../img/fluxplot.png +.. image:: ../_images/fluxplot.png :height: 200px Alternatively, the user can write their own Python script to manipulate the data @@ -249,7 +249,7 @@ two heatmaps in the previous figure. Plotting in 3D -------------- -.. image:: ../../img/3dcore.png +.. image:: ../_images/3dcore.png :height: 200px As with 3D plots of the geometry, meshtally data needs to be put into a standard @@ -357,7 +357,7 @@ and dumped to MATLAB in one step. Particle Track Visualization ---------------------------- -.. image:: ../../img/Tracks.png +.. image:: ../_images/Tracks.png :height: 200px OpenMC can dump particle tracks—the position of particles as they are From c439c3dc5d533c74ca7df25759009eb89d81af9a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 Jan 2014 20:34:29 -0500 Subject: [PATCH 30/31] Change OPENMP macro to _OPENMP which is automatically defined when compiling with OpenMP. --- src/Makefile | 6 +++--- src/eigenvalue.F90 | 4 ++-- src/global.F90 | 6 +++--- src/initialize.F90 | 6 +++--- src/input_xml.F90 | 2 +- src/output.F90 | 4 ++-- 6 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Makefile b/src/Makefile index 03c370e07a..00abca489f 100644 --- a/src/Makefile +++ b/src/Makefile @@ -181,17 +181,17 @@ endif ifeq ($(OPENMP),yes) ifeq ($(COMPILER),intel) - F90FLAGS += -openmp -DOPENMP + F90FLAGS += -openmp LDFLAGS += -openmp endif ifeq ($(COMPILER),gnu) - F90FLAGS += -fopenmp -DOPENMP + F90FLAGS += -fopenmp LDFLAGS += -fopenmp endif ifeq ($(COMPILER),ibm) - F90FLAGS += -qsmp=omp -WF,-DOPENMP + F90FLAGS += -qsmp=omp LDFLAGS += -qsmp=omp endif endif diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index f5375bf3a9..4047f815f5 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -166,7 +166,7 @@ contains subroutine finalize_generation() -#ifdef OPENMP +#ifdef _OPENMP ! Join the fission bank from each thread into one global fission bank call join_bank_from_threads() #endif @@ -821,7 +821,7 @@ contains end subroutine replay_batch_history -#ifdef OPENMP +#ifdef _OPENMP !=============================================================================== ! JOIN_BANK_FROM_THREADS !=============================================================================== diff --git a/src/global.F90 b/src/global.F90 index 637ea315e6..c3dbf0aa24 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -158,7 +158,7 @@ module global ! Source and fission bank type(Bank), allocatable, target :: source_bank(:) type(Bank), allocatable, target :: fission_bank(:) -#ifdef OPENMP +#ifdef _OPENMP type(Bank), allocatable, target :: master_fission_bank(:) #endif integer(8) :: n_bank ! # of sites in fission bank @@ -205,7 +205,7 @@ module global integer :: MPI_BANK ! MPI datatype for fission bank integer :: MPI_TALLYRESULT ! MPI datatype for TallyResult -#ifdef OPENMP +#ifdef _OPENMP integer :: n_threads = NONE ! number of OpenMP threads integer :: thread_id ! ID of a given thread #endif @@ -438,7 +438,7 @@ contains !$omp parallel if (allocated(fission_bank)) deallocate(fission_bank) !$omp end parallel -#ifdef OPENMP +#ifdef _OPENMP if (allocated(master_fission_bank)) deallocate(master_fission_bank) #endif if (allocated(source_bank)) deallocate(source_bank) diff --git a/src/initialize.F90 b/src/initialize.F90 index b75fa20ef8..9f495e4ac9 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -26,7 +26,7 @@ module initialize use mpi #endif -#ifdef OPENMP +#ifdef _OPENMP use omp_lib #endif @@ -372,7 +372,7 @@ contains ! Read number of threads i = i + 1 -#ifdef OPENMP +#ifdef _OPENMP ! Read and set number of OpenMP threads n_threads = str_to_int(argv(i)) if (n_threads < 1) then @@ -830,7 +830,7 @@ contains call fatal_error() end if -#ifdef OPENMP +#ifdef _OPENMP ! If OpenMP is being used, each thread needs its own private fission ! bank. Since the private fission banks need to be combined at the end of a ! generation, there is also a 'master_fission_bank' that is used to collect diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3ce3c9d0d4..ae4e4b1b2e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -228,7 +228,7 @@ contains ! Number of OpenMP threads if (check_for_node(doc, "threads")) then -#ifdef OPENMP +#ifdef _OPENMP if (n_threads == NONE) then call get_node_value(doc, "threads", n_threads) if (n_threads < 1) then diff --git a/src/output.F90 b/src/output.F90 index 677c35396c..3e31ad7183 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -31,7 +31,7 @@ contains subroutine title() -#ifdef OPENMP +#ifdef _OPENMP use omp_lib #endif @@ -69,7 +69,7 @@ contains trim(to_str(n_procs)) #endif -#ifdef OPENMP +#ifdef _OPENMP ! Write number of OpenMP threads write(UNIT=OUTPUT_UNIT, FMT='(6X,"OpenMP Threads:",1X,A)') & trim(to_str(omp_get_max_threads())) From b142ff6458d4515a2e4f9d72ae0b7ee3bfb67da6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 Jan 2014 21:26:40 -0500 Subject: [PATCH 31/31] Move calculation of number of threads outside of .omp parallel. --- src/initialize.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 9f495e4ac9..c1fbf40cda 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -836,8 +836,9 @@ contains ! generation, there is also a 'master_fission_bank' that is used to collect ! the sites from each thread. + n_threads = omp_get_max_threads() + !$omp parallel - n_threads = omp_get_num_threads() thread_id = omp_get_thread_num() if (thread_id == 0) then