diff --git a/src/Makefile b/src/Makefile index 359e34c36..84b0ea254 100644 --- a/src/Makefile +++ b/src/Makefile @@ -11,6 +11,7 @@ modules = ace.f90 \ geometry.f90 \ global.f90 \ mcnp_random.f90 \ + mpi_routines.f90 \ output.f90 \ physics.f90 \ search.f90 \ @@ -23,20 +24,24 @@ test_objects = $(modules:.f90=.o) $(test:.f90=.o) #-------------------------------------------------------------------- # Compiler Options -F90 = gfortran -F90FLAGS = -g +F90 = ifort +F90FLAGS = -g -fpp +MPI = /opt/mpich2/1.3.3-gcc +MPIF90 = $(MPI)/bin/mpif90 +MPI_COMPILE_FLAGS = $(shell mpif90 -f90='' -compile_info) +MPI_LINK_FLAGS = $(shell mpif90 -f90='' -link_info) #-------------------------------------------------------------------- # Targets all: $(program) $(program): $(main_objects) - $(F90) $(F90FLAGS) $(main_objects) -o $@ + $(MPIF90) -f90=$(F90) $(MPI_LINK_FLAGS) $(F90FLAGS) $(main_objects) -o $@ clean: @rm -f *.o *.mod $(program) neat: @rm -f *.o *.mod unittest: $(test_objects) - $(F90) $(F90FLAGS) $(test_objects) -o $@ + $(F90) $(F90FLAGS) $(MPI_COMPILE_FLAGS) $(test_objects) -o $@ #-------------------------------------------------------------------- # Rules & misc @@ -59,7 +64,8 @@ global.o: types.o main.o: global.o fileio.o output.o geometry.o mcnp_random.o \ source.o physics.o cross_section.o data_structures.o \ ace.o energy_grid.o -output.o: global.o types.o data_structures.o +mpi_routines.o: global.o output.o types.o +output.o: global.o types.o data_structures.o endf.o physics.o: types.o global.o mcnp_random.o geometry.o output.o \ search.o endf.o search.o: output.o diff --git a/src/ace.f90 b/src/ace.f90 index bdd61a677..fe858e22d 100644 --- a/src/ace.f90 +++ b/src/ace.f90 @@ -38,7 +38,7 @@ contains integer :: n integer :: index_continuous integer :: index_thermal - type(DictionaryCI), pointer :: temp_dict + type(DictionaryCI), pointer :: temp_dict => null() call dict_create(ace_dict) @@ -206,8 +206,9 @@ contains end do - call parse_ESZ(table) - call parse_MTR(table) + call read_esz(table) + call read_nu_data(table) + call read_reactions(table) call read_angular_dist(table) call read_energy_dist(table) @@ -220,14 +221,14 @@ contains end subroutine read_ACE_continuous !===================================================================== -! PARSE_ESZ - reads through the ESZ block. This block contains the +! READ_ESZ - reads through the ESZ block. This block contains the ! energy grid, total xs, absorption xs, elastic scattering xs, and ! heating numbers. !===================================================================== - subroutine parse_ESZ(table) + subroutine read_esz(table) - type(AceContinuous), pointer, intent(inout) :: table + type(AceContinuous), pointer :: table integer :: NE @@ -254,19 +255,170 @@ contains table%sigma_el = get_real(NE) table%heating = get_real(NE) - end subroutine parse_ESZ + end subroutine read_esz !===================================================================== -! PARSE_MTR - Get the list of reaction MTs for this cross-section -! table. The MT values are somewhat arbitrary. Also read in Q-values, -! neutron multiplicities, and cross-sections. -! ===================================================================== +! READ_NU_DATA reads data given on the number of neutrons emitted from +! fission as a function of the incoming energy of a neutron. This data +! may be broken down into prompt and delayed neutrons emitted as well. +!===================================================================== - subroutine parse_MTR(table) + subroutine read_nu_data(table) type(AceContinuous), pointer :: table - type(AceReaction), pointer :: rxn + integer :: JXS2 + integer :: KNU ! Location for nu data + 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 :: length ! Length of data block to allocate + + integer :: JXS24 + integer :: DEC1 + integer :: DEC2 + + JXS2 = JXS(2) + JXS24 = JXS(24) + + if (JXS2 == 0) then + ! ============================================================= + ! No prompt/total nu data is present + table % nu_t_type = NU_NONE + table % nu_p_type = NU_NONE + + elseif (XSS(JXS2) > 0) then + ! ============================================================= + ! Prompt or total nu data is present + KNU = JXS2 + LNU = int(XSS(KNU)) + if (LNU == 1) then + ! Polynomial data + table % nu_t_type = NU_POLYNOMIAL + table % nu_p_type = NU_NONE + + ! allocate determine how many coefficients for polynomial + NC = int(XSS(KNU+1)) + length = NC + 1 + elseif (LNU == 2) then + ! Tabular data + table % nu_t_type = NU_TABULAR + table % nu_p_type = NU_NONE + + ! determine number of interpolation regions and number of + ! energies + NR = int(XSS(KNU+1)) + NE = int(XSS(KNU+2+2*NR)) + length = 2 + 2*NR + 2*NE + end if + + ! allocate space for nu data storage + allocate(table % nu_t_data(length)) + + ! read data -- for polynomial, this is the number of + ! coefficients and the coefficients themselves, and for + ! tabular, this is interpolation data and tabular E/nu + XSS_index = KNU + 1 + table % nu_t_data = get_real(length) + + elseif (XSS(JXS2) < 0) then + ! ============================================================= + ! Prompt and total nu data is present -- read prompt data first + KNU = JXS2 + 1 + LNU = XSS(KNU) + if (LNU == 1) then + ! Polynomial data + table % nu_p_type = NU_POLYNOMIAL + + ! allocate determine how many coefficients for polynomial + NC = XSS(KNU+1) + length = NC + 1 + elseif (LNU == 2) then + ! Tabular data + table % nu_p_type = NU_TABULAR + + ! determine number of interpolation regions and number of + ! energies + NR = XSS(KNU+1) + NE = XSS(KNU+2+2*NR) + length = 2 + 2*NR + 2*NE + end if + + ! allocate space for nu data storage + allocate(table % nu_p_data(length)) + + ! read data + XSS_index = KNU + 1 + table % nu_p_data = get_real(length) + + ! Now read total nu data + KNU = JXS2 + abs(XSS(JXS2)) + 1 + LNU = XSS(KNU) + if (LNU == 1) then + ! Polynomial data + table % nu_t_type = NU_POLYNOMIAL + + ! allocate determine how many coefficients for polynomial + NC = int(XSS(KNU+1)) + length = NC + 1 + elseif (LNU == 2) then + ! Tabular data + table % nu_t_type = NU_TABULAR + + ! determine number of interpolation regions and number of + ! energies + NR = int(XSS(KNU+1)) + NE = int(XSS(KNU+2+2*NR)) + length = 2 + 2*NR + 2*NE + end if + + ! allocate space for nu data storage + allocate(table % nu_t_data(length)) + + ! read data + XSS_index = KNU + 1 + table % nu_t_data = get_real(length) + end if + + if (JXS24 > 0) then + ! ============================================================= + ! Delayed nu data is present + + table % nu_d_type = NU_TABULAR + KNU = JXS24 + + ! determine size of tabular delayed nu data + NR = int(XSS(KNU+1)) + NE = int(XSS(KNU+2+2*NR)) + length = 2 + 2*NR + 2*NE + + ! allocate space for delayed nu data + allocate(table % nu_d_data(length)) + + ! read delayed nu data + XSS_index = KNU + 1 + table % nu_d_data = get_real(length) + + ! TODO: Read secondary energy distribution + ! TODO: Read precursor data + else + table % nu_d_type = NU_NONE + end if + + end subroutine read_nu_data + +!===================================================================== +! READ_REACTIONS - Get the list of reaction MTs for this cross-section +! table. The MT values are somewhat arbitrary. Also read in Q-values, +! neutron multiplicities, and cross-sections. +!===================================================================== + + subroutine read_reactions(table) + + type(AceContinuous), pointer :: table + + type(AceReaction), pointer :: rxn => null() integer :: LMT ! index of MT list in XSS integer :: NMT ! Number of reactions integer :: JXS4 ! index of Q values in XSS @@ -291,7 +443,7 @@ contains ! Store elastic scattering cross-section on reaction one rxn => table%reactions(1) rxn%MT = 2 - rxn%Q_value = 0.0_8 + rxn%Q_value = ZERO rxn%TY = 1 rxn%IE = 1 allocate(rxn%sigma(table%n_grid)) @@ -318,7 +470,7 @@ contains rxn % has_energy_dist = .false. end do - end subroutine parse_MTR + end subroutine read_reactions !===================================================================== ! READ_ANGULAR_DIST parses the angular distribution for each reaction @@ -328,7 +480,7 @@ contains type(AceContinuous), pointer :: table - type(AceReaction), pointer :: rxn + type(AceReaction), pointer :: rxn => null() integer :: JXS8 integer :: JXS9 integer :: NMT @@ -337,7 +489,7 @@ contains integer :: NP integer :: LC integer :: i, j - integer :: size + integer :: length JXS8 = JXS(8) JXS9 = JXS(9) @@ -363,49 +515,49 @@ contains ! allocate space for incoming energies and locations NE = XSS(JXS9 + LOCB - 1) - rxn % adist_n_energy = NE - allocate(rxn % adist_energy(NE)) - allocate(rxn % adist_type(NE)) - allocate(rxn % adist_location(NE)) + rxn % adist % n_energy = NE + allocate(rxn % adist % energy(NE)) + allocate(rxn % adist % type(NE)) + allocate(rxn % adist % location(NE)) ! read incoming energy grid and location of tables XSS_index = JXS9 + LOCB - rxn % adist_energy = get_real(NE) - rxn % adist_location = get_int(NE) + rxn % adist % energy = get_real(NE) + rxn % adist % location = get_int(NE) ! determine dize of data block - size = 0 + length = 0 do j = 1, NE - LC = rxn % adist_location(j) + LC = rxn % adist % location(j) if (LC == 0) then ! isotropic - rxn % adist_type(j) = ANGLE_ISOTROPIC + rxn % adist % type(j) = ANGLE_ISOTROPIC elseif (LC > 0) then ! 32 equiprobable bins - rxn % adist_type(j) = ANGLE_32_EQUI - size = size + 33 + rxn % adist % type(j) = ANGLE_32_EQUI + length = length + 33 elseif (LC < 0) then ! tabular distribution - rxn % adist_type(j) = ANGLE_TABULAR + rxn % adist % type(j) = ANGLE_TABULAR NP = XSS(JXS9 + abs(LC)) - size = size + 2 + 3*NP + length = length + 2 + 3*NP end if end do ! allocate angular distribution data and read - allocate(rxn % adist_data(size)) + allocate(rxn % adist % data(length)) ! read angular distribution -- currently this does not actually ! parse the angular distribution tables for each incoming ! energy, that must be done on-the-fly - LC = rxn % adist_location(1) + LC = rxn % adist % location(1) XSS_index = JXS9 + abs(LC) - 1 - rxn % adist_data = get_real(size) + rxn % adist % data = get_real(length) ! change location pointers since they are currently relative to ! JXS(9) - LC = abs(rxn % adist_location(1)) - rxn % adist_location = abs(rxn % adist_location) - LC + 1 + LC = abs(rxn % adist % location(1)) + rxn % adist % location = abs(rxn % adist % location) - LC end do @@ -419,7 +571,7 @@ contains type(AceContinuous), pointer :: table - type(AceReaction), pointer :: rxn + type(AceReaction), pointer :: rxn => null() integer :: LOCC integer :: LED ! location of LDLW block integer :: LDIS ! location of DLW block @@ -431,7 +583,7 @@ contains integer :: NP2 integer :: NRa, NEa, NRb, NEb integer :: IDAT - integer :: start, size, size_interp_data + integer :: start, length, length_interp_data integer :: i, j, k LED = JXS(10) @@ -440,6 +592,7 @@ contains ! Loop over all reactions do i = 1, NXS(5) rxn => table % reactions(i+1) ! skip over elastic scattering + rxn % has_energy_dist = .true. ! find location of energy distribution data LOCC = XSS(LED + i - 1) @@ -448,80 +601,87 @@ contains LAW = XSS(LDIS + LOCC) IDAT = XSS(LDIS + LOCC + 1) NR = XSS(LDIS + LOCC + 2) + rxn % edist % law = LAW ! allocate space for ENDF interpolation parameters if (NR > 0) then - allocate(rxn % edist_nbt(NR)) - allocate(rxn % edist_int(NR)) + allocate(rxn % edist % nbt(NR)) + allocate(rxn % edist % int(NR)) end if ! read ENDF interpolation parameters XSS_index = LDIS + LOCC + 3 - rxn % edist_nbt = get_real(NR) - rxn % edist_int = get_real(NR) + rxn % edist % nbt = get_real(NR) + rxn % edist % int = get_real(NR) ! allocate space for law validity data NE = XSS(LDIS + LOCC + 3 + 2*NR) - allocate(rxn % edist_energy(NE)) - allocate(rxn % edist_pvalid(NE)) + allocate(rxn % edist % energy(NE)) + allocate(rxn % edist % pvalid(NE)) - size_interp_data = 5 + 2*(NR + NE) + length_interp_data = 5 + 2*(NR + NE) ! read law validity data XSS_index = LDIS + LOCC + 4 + 2*NR - rxn % edist_energy = get_real(NE) - rxn % edist_pvalid = get_real(NE) + rxn % edist % energy = get_real(NE) + rxn % edist % pvalid = get_real(NE) ! Set index to beginning of IDAT array start = LDIS + IDAT - 2 ! Determine size of LDAT array based on which secondary energy ! law it is - size = 0 + length = 0 select case (LAW) case (1) ! Tabular equiprobable energy bins NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) NET = XSS(start + 3 + 2*NR + NE) - size = 3 + 2*NR + NE + 3*NET*NE + length = 3 + 2*NR + NE + 3*NET*NE case (2) ! Discrete photon energy - size = 2 + length = 2 case (3) ! Level scattering - size = 2 + length = 2 case (4) ! Continuous tabular distribution NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) - size = size + 2 + 2*NR + 2*NE + length = length + 2 + 2*NR + 2*NE do j = 1,NE - NP = XSS(start + size + 2) - size = size + 2 + 3*NP + ! determine length + NP = XSS(start + length + 2) + length = length + 2 + 3*NP + + ! adjust location for this block + k = start + 2 + 2*NR + NE + j + XSS(k) = XSS(k) - LOCC - length_interp_data end do + k = start + 2 + 2*NR + NE case (5) ! General evaporation spectrum NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) NET = XSS(start + 3 + 2*NR + 2*NE) - size = 3 + 2*NR + 2*NE + NET + length = 3 + 2*NR + 2*NE + NET case (7) ! Maxwell fission spectrum NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) - size = 3 + 2*NR + 2*NE + length = 3 + 2*NR + 2*NE case (9) ! Evaporation spectrum NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) - size = 3 + 2*NR + 2*NE + length = 3 + 2*NR + 2*NE case (11) ! Watt spectrum @@ -529,53 +689,57 @@ contains NEa = XSS(start + 2 + 2*NRa) NRb = XSS(start + 3 + 2*(NRa+NEa)) NEb = XSS(start + 4 + 2*(NRa+NEa+NRb)) - size = 5 + 2*(NRa + NEa + NRb + NEb) + length = 5 + 2*(NRa + NEa + NRb + NEb) case (44) ! Kalbach-Mann correlated scattering NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) - size = size + 2 + 2*NR + 2*NE + length = length + 2 + 2*NR + 2*NE do j = 1,NE - NP = XSS(start + size + 2) - size = size + 2 + 5*NP + NP = XSS(start + length + 2) + length = length + 2 + 5*NP + + ! adjust location for this block + k = start + 2 + 2*NR + NE + j + XSS(k) = XSS(k) - LOCC - length_interp_data end do case (61) ! Correlated energy and angle distribution NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) - size = size + 2 + 2*NR + 2*NE + length = length + 2 + 2*NR + 2*NE do j = 1,NE ! outgoing energy distribution - NP = XSS(start + size + 2) - size = size + 2 + 4*NP + NP = XSS(start + length + 2) + length = length + 2 + 4*NP do k = 1, NP ! outgoing angle distribution - NP2 = XSS(start + size + 2) - size = size + 2 + 3*NP2 + NP2 = XSS(start + length + 2) + length = length + 2 + 3*NP2 end do end do case (66) ! N-body phase space distribution - size = 2 + length = 2 case (67) ! Laboratory energy-angle law NR = XSS(start + 1) NE = XSS(start + 2 + 2*NR) NMU = XSS(start + 4 + 2*NR + 2*NE) - size = 4 + 2*(NR + NE + NMU) + length = 4 + 2*(NR + NE + NMU) end select ! allocate secondary energy distribution array - allocate(rxn % edist_data(size)) + allocate(rxn % edist % data(length)) ! read secondary energy distribution XSS_index = start + 1 - rxn % edist_data = get_real(size) + rxn % edist % data = get_real(length) end do diff --git a/src/fileio.f90 b/src/fileio.f90 index 73dfbe070..0b1f51054 100644 --- a/src/fileio.f90 +++ b/src/fileio.f90 @@ -271,6 +271,13 @@ contains case ('xs_data') path_xsdata = words(2) + case ('verbosity') + verbosity = str_to_int(words(2)) + if (verbosity == ERROR_CODE) then + msg = "Invalid verbosity: " // words(2) + call error(msg) + end if + case default ! do nothing end select diff --git a/src/geometry.f90 b/src/geometry.f90 index 483d61c04..5873b4b55 100644 --- a/src/geometry.f90 +++ b/src/geometry.f90 @@ -107,8 +107,6 @@ contains c => cells(univ % cells(i)) if (cell_contains(c, p)) then - p%cell = dict_get_key(cell_dict, c % uid) - ! If this cell contains a universe of lattice, search for ! the particle in that universe/lattice if (c % fill > 0) then @@ -125,6 +123,11 @@ contains found = .true. cCell => c cMaterial => materials(cCell%material) + cUniverse => univ + + ! set particle attributes + p%cell = dict_get_key(cell_dict, c % uid) + p%universe = dict_get_key(universe_dict, univ % uid) exit end if end if diff --git a/src/global.f90 b/src/global.f90 index de41e590e..21332ddab 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -56,6 +56,18 @@ module global ! Source and fission bank type(Particle), allocatable, target :: source_bank(:) type(Bank), allocatable, target :: fission_bank(:) + integer :: n_bank ! # of sites in fission bank + integer :: bank_first ! index of first particle in bank + integer :: bank_last ! index of last particle in bank + + ! cycle keff + real(8) :: keff + + ! Parallel processing variables + integer :: n_procs ! number of processes + integer :: rank ! rank of process + logical :: master ! master process? + logical :: mpi_enabled ! is MPI in use and initialized? ! Paths to input file, cross section data, etc character(100) :: & @@ -72,7 +84,8 @@ module global & K_BOLTZMANN = 8.617342e-5, & ! Boltzmann constant in eV/K & INFINITY = huge(0.0_8), & ! positive infinity & ZERO = 0.0_8, & - & ONE = 1.0_8 + & ONE = 1.0_8, & + & TWO = 2.0_8 ! Boundary conditions integer, parameter :: & @@ -146,19 +159,28 @@ module global & PHOTON = 2, & & ELECTRON = 3 + ! Fission neutron emission (nu) type + integer, parameter :: & + & NU_NONE = 0, & ! No nu values (non-fissionable) + & NU_POLYNOMIAL = 1, & ! Nu values given by polynomial + & NU_TABULAR = 2 ! Nu values given by tabular distribution + ! Integer code for read error -- better hope this number is never ! used in an input file! integer, parameter :: & & ERROR_CODE = -huge(0) - integer :: verbosity = 5 + ! The verbosity controls how much information will be printed to the + ! screen and in logs + integer :: verbosity + integer, parameter :: max_words = 500 integer, parameter :: max_line = 250 ! Versioning numbers integer, parameter :: VERSION_MAJOR = 0 integer, parameter :: VERSION_MINOR = 2 - integer, parameter :: VERSION_RELEASE = 0 + integer, parameter :: VERSION_RELEASE = 1 contains @@ -174,6 +196,9 @@ contains ! Default number of particles n_particles = 10000 + ! Default verbosity + verbosity = 5 + end subroutine set_defaults !===================================================================== @@ -183,6 +208,8 @@ contains subroutine free_memory() + integer :: ierr + ! Deallocate cells, surfaces, materials if (allocated(cells)) deallocate(cells) if (allocated(surfaces)) deallocate(surfaces) @@ -195,6 +222,9 @@ contains if (allocated(fission_bank)) deallocate(fission_bank) if (allocated(source_bank)) deallocate(source_bank) + ! If MPI is in use and enabled, terminate it + call MPI_FINALIZE(ierr) + ! End program stop diff --git a/src/input_sample b/src/input_sample index b4319adf7..a90cdaa0b 100644 --- a/src/input_sample +++ b/src/input_sample @@ -20,4 +20,5 @@ source box -4 -4 -4 4 4 4 xs_library endfb7 xs_data /opt/serpent/xsdata/endfb7 -criticality 1 1 15 +criticality 1 1 100 +verbosity 8 diff --git a/src/main.f90 b/src/main.f90 index baedd5cd6..0e7855b5c 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -12,6 +12,7 @@ program main use cross_section, only: read_xsdata, material_total_xs use ace, only: read_xs use energy_grid, only: unionized_grid, original_indices + use mpi_routines implicit none @@ -19,12 +20,16 @@ program main character(250) :: msg type(Universe), pointer :: univ - ! Print the OpenMC title and version/date/time information - call title() - verbosity = 10 + ! Setup MPI + call setup_mpi() - ! Initialize random number generator - call RN_init_problem( 3, 0_8, 0_8, 0_8, 0 ) + ! Print the OpenMC title and version/date/time information + if (master) call title() + + ! Initialize random number generator. The first argument corresponds + ! to which random number generator to use- in this case one of the + ! L'Ecuyer 63-bit RNGs. + call RN_init_problem(3, 0_8, 0_8, 0_8, 0) ! Set default values for settings call set_defaults() @@ -65,8 +70,10 @@ program main ! calculate total material cross-sections for sampling path lenghts call material_total_xs() - call echo_input() - call print_summary() + if (master) then + ! call echo_input() + ! call print_summary() + end if ! create source particles call init_source() @@ -95,6 +102,7 @@ contains CYCLE_LOOP: do i_cycle = 1, n_cycles ! Set all tallies to zero + n_bank = 0 HISTORY_LOOP: do j = 1, n_particles @@ -115,6 +123,7 @@ contains end do HISTORY_LOOP ! Collect results and statistics + print *, n_bank ! print cycle information diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90 new file mode 100644 index 000000000..7eae9e6e3 --- /dev/null +++ b/src/mpi_routines.f90 @@ -0,0 +1,56 @@ +module mpi_routines + + use mpi + use global + use output, only: message, error + + implicit none + + integer :: MPI_bank ! MPI datatype for fission bank + integer(8) :: bank_index ! Fission bank site unique identifier + +contains + +!===================================================================== +! SETUP_MPI initilizes the Message Passing Interface (MPI) and +! determines the number of processors the problem is being run with as +! well as the rank of each processor. +!===================================================================== + + subroutine setup_mpi() + + integer :: ierr + character(250) :: msg + + ! Initialize MPI + call MPI_INIT(ierr) + if (ierr /= MPI_SUCCESS) then + msg = "Failed to initialize MPI." + call error(msg) + end if + mpi_enabled = .true. + + ! Determine number of processors + call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, ierr) + if (ierr /= MPI_SUCCESS) then + msg = "Could not determine number of processors." + call error(msg) + end if + + ! Determine rank of each processor + call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr) + if (ierr /= MPI_SUCCESS) then + msg = "Could not determine MPI rank." + call error(msg) + end if + + ! Determine master + if (rank == 0) then + master = .true. + else + master = .false. + end if + + end subroutine setup_mpi + +end module mpi_routines diff --git a/src/output.f90 b/src/output.f90 index 4ec9f6b2a..4628c1897 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -4,9 +4,13 @@ module output use global use types, only: Cell, Universe, Surface use data_structures, only: dict_get_key + use endf, only: reaction_name implicit none + integer :: ou = OUTPUT_UNIT + integer :: eu = ERROR_UNIT + contains !===================================================================== @@ -19,9 +23,6 @@ module output character(10) :: date character(8) :: time - integer :: ou - - ou = OUTPUT_UNIT write(ou,*) write(ou,*) ' .d88888b. 888b d888 .d8888b. ' @@ -61,9 +62,6 @@ module output subroutine echo_input() character(32) :: string - integer :: ou - - ou = OUTPUT_UNIT ! Display problem summary write(ou,*) '=============================================' @@ -104,18 +102,17 @@ module output character(*), intent(in) :: msg integer, intent(in) :: level - integer :: ou integer :: n_lines integer :: i + ! Only allow master to print to screen + if (.not. master) return + if ( level <= verbosity ) then - ou = OUTPUT_UNIT - n_lines = (len_trim(msg)-1)/79 + 1 do i = 1, n_lines write(ou, fmt='(1X,A79)') msg(79*(i-1)+1:79*i) end do - end if end subroutine message @@ -130,9 +127,10 @@ module output character(*), intent(in) :: msg integer :: n_lines - integer :: i, ou + integer :: i - ou = OUTPUT_UNIT + ! Only allow master to print to screen + if (.not. master) return write(ou, fmt='(1X,A9)', advance='no') 'WARNING: ' @@ -158,9 +156,10 @@ module output character(*), intent(in) :: msg integer :: n_lines - integer :: i, eu + integer :: i - eu = ERROR_UNIT + ! Only allow master to print to screen + if (.not. master) return write(eu, fmt='(1X,A7)', advance='no') 'ERROR: ' @@ -183,7 +182,7 @@ module output ! execution and returns it in a readable format !===================================================================== - subroutine get_today( today_date, today_time ) + subroutine get_today(today_date, today_time) character(10), intent(out) :: today_date character(8), intent(out) :: today_time @@ -220,6 +219,91 @@ module output end subroutine get_today +!===================================================================== +! PRINT_PARTICLE displays the attributes of a particle +!===================================================================== + + subroutine print_particle(p) + + type(Particle), pointer :: p + + integer :: i + type(Cell), pointer :: c => null() + type(Surface), pointer :: s => null() + type(Universe), pointer :: u => null() + character(250) :: string + + select case (p % type) + case (NEUTRON) + write(ou,*) 'Neutron ' // int_to_str(p % uid) + case (PHOTON) + write(ou,*) 'Photon ' // int_to_str(p % uid) + case (ELECTRON) + write(ou,*) 'Electron ' // int_to_str(p % uid) + case default + write(ou,*) 'Unknown Particle ' // int_to_str(p % uid) + end select + write(ou,100) 'x = ', p % xyz(1) + write(ou,100) 'y = ', p % xyz(2) + write(ou,100) 'z = ', p % xyz(3) + write(ou,100) 'u = ', p % uvw(1) + write(ou,100) 'v = ', p % uvw(2) + write(ou,100) 'w = ', p % uvw(3) + write(ou,100) 'Weight = ', p % wgt + write(ou,100) 'Energy = ', p % E + write(ou,*) ' IE = ' // int_to_str(p % IE) + write(ou,100) 'Interpolation factor = ', p % interp + + if (p % cell > 0) then + c => cells(p % cell) + write(ou,*) ' Cell = ' // int_to_str(c % uid) + else + write(ou,*) ' Cell not determined' + end if + + if (p % surface > 0) then + s => surfaces(p % surface) + write(ou,*) ' Surface = ' // int_to_str(s % uid) + else + write(ou,*) ' Surface = None' + end if + + u => universes(p % universe) + write(ou,*) ' Universe = ' // int_to_str(u % uid) + write(ou,*) + + ! Format for a single real + 100 format (5X,A,G10.3) + + nullify(c) + nullify(s) + nullify(u) + + end subroutine print_particle + +!===================================================================== +! PRINT_REACTION displays the attributes of a reaction +!===================================================================== + + subroutine print_reaction(rxn) + + type(AceReaction), pointer :: rxn + + write(ou,*) 'Reaction ' // reaction_name(rxn % MT) + write(ou,*) ' MT = ' // int_to_str(rxn % MT) + write(ou,100) 'Q-value = ', rxn % Q_value + write(ou,*) ' TY = ' // int_to_str(rxn % TY) + write(ou,*) ' Starting index = ' // int_to_str(rxn % IE) + if (rxn % has_energy_dist) then + write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law) + end if + write(ou,*) + + ! Format for a single real + 100 format (5X,A,G10.3) + + end subroutine print_reaction + !===================================================================== ! PRINT_CELL displays the attributes of a cell !===================================================================== @@ -228,15 +312,12 @@ module output type(Cell), pointer :: c - integer :: ou integer :: temp integer :: i type(Universe), pointer :: u => null() type(Material), pointer :: m => null() character(250) :: string - ou = OUTPUT_UNIT - write(ou,*) 'Cell ' // int_to_str(c % uid) temp = dict_get_key(cell_dict, c % uid) write(ou,*) ' Array Index = ' // int_to_str(temp) @@ -287,12 +368,9 @@ module output type(Universe), pointer :: univ - integer :: ou integer :: i character(250) :: string - type(Cell), pointer :: c - - ou = OUTPUT_UNIT + type(Cell), pointer :: c => null() write(ou,*) 'Universe ' // int_to_str(univ % uid) write(ou,*) ' Level = ' // int_to_str(univ % level) @@ -304,6 +382,8 @@ module output write(ou,*) ' Cells =' // trim(string) write(ou,*) + nullify(c) + end subroutine print_universe !===================================================================== @@ -314,12 +394,9 @@ module output type(Surface), pointer :: surf - integer :: ou integer :: i character(80) :: string - ou = OUTPUT_UNIT - write(ou,*) 'Surface ' // int_to_str(surf % uid) select case (surf % type) case (SURF_PX) @@ -354,7 +431,7 @@ module output string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i)) end do end if - write(ou,*) ' Positive Neighbors = ', string + write(ou,*) ' Positive Neighbors = ', trim(string) string = "" if (allocated(surf % neighbor_neg)) then @@ -362,7 +439,7 @@ module output string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i)) end do end if - write(ou,*) ' Negative Neighbors =', string + write(ou,*) ' Negative Neighbors =', trim(string) end subroutine print_surface @@ -374,14 +451,11 @@ module output type(Material), pointer :: mat - integer :: ou integer :: i integer :: n_lines type(AceContinuous), pointer :: table character(250) :: string - ou = OUTPUT_UNIT - write(ou,*) 'Material ' // int_to_str(mat % uid) ! Make string of all isotopes string = "" @@ -403,6 +477,8 @@ module output & ' atom/b-cm' write(ou,*) + nullify(table) + end subroutine print_material !===================================================================== @@ -418,9 +494,6 @@ module output type(Universe), pointer :: u type(Material), pointer :: m integer :: i - integer :: ou - - ou = OUTPUT_UNIT ! print summary of cells write(ou,*) '=============================================' @@ -462,6 +535,11 @@ module output call print_material(m) end do + nullify(s) + nullify(c) + nullify(u) + nullify(m) + end subroutine print_summary end module output diff --git a/src/physics.f90 b/src/physics.f90 index aed1a3cf0..0d9c3edd6 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -175,7 +175,7 @@ contains ! sample nuclide r1 = rang() - prob = 0.0_8 + prob = ZERO total = sum(Sigma_t) do i = 1, n_isotopes prob = prob + Sigma_t(i) / total @@ -194,7 +194,7 @@ contains ! sample reaction channel r1 = rang()*Sigma - prob = 0.0_8 + prob = ZERO do i = 1, table%n_reaction rxn => table%reactions(i) if (rxn%MT >= 200) cycle @@ -213,6 +213,10 @@ contains select case (rxn%MT) case (2) call elastic_scatter(p, table%awr) + case (18:21, 38) + call n_fission(p, table, rxn) + case (51:90) + call level_scatter(p, table, rxn) case (102) call n_gamma(p) case default @@ -246,9 +250,9 @@ contains vz = vel*p % uvw(3) - vcx = vx/(awr + 1.0_8) - vcy = vy/(awr + 1.0_8) - vcz = vz/(awr + 1.0_8) + vcx = vx/(awr + ONE) + vcy = vy/(awr + ONE) + vcz = vz/(awr + ONE) ! Transform to CM frame vx = vx - vcx @@ -259,11 +263,11 @@ contains ! Select isotropic direcion -- this is only valid for s-wave ! scattering - phi = 2.0_8*PI*rang() - mu = 2.0_8*rang() - 1.0_8 + phi = TWO*PI*rang() + mu = TWO*rang() - ONE u = mu - v = sqrt(1.0_8 - mu**2) * cos(phi) - w = sqrt(1.0_8 - mu**2) * sin(phi) + v = sqrt(ONE - mu**2) * cos(phi) + w = sqrt(ONE - mu**2) * sin(phi) vx = u*vel vy = v*vel @@ -288,12 +292,253 @@ contains end subroutine elastic_scatter !===================================================================== -! LEVEL_INELASTIC +! N_FISSION determines the average total, prompt, and delayed neutrons +! produced from fission and creates appropriate bank sites. This +! routine will not work with implicit absorption, namely sampling of +! the number of neutrons! !===================================================================== - subroutine level_inelastic + subroutine n_fission(p, table, rxn) - end subroutine level_inelastic + type(Particle), pointer :: p + type(AceContinuous), pointer :: table + type(AceReaction), pointer :: rxn + + integer :: i ! loop index + integer :: j ! index on nu energy grid + integer :: loc ! index before start of energies/nu values + integer :: NC ! number of polynomial coefficients + integer :: NR ! number of interpolation regions + integer :: NE ! number of energies tabulated + real(8) :: E ! incoming energy of neutron + real(8) :: E_out ! outgoing energy of fission neutron + real(8) :: mu_out ! outgoing angle of fission neutron (if needed) + real(8) :: c ! polynomial coefficient + real(8) :: f ! interpolation factor + real(8) :: nu_total ! total nu + real(8) :: nu_prompt ! prompt nu + real(8) :: nu_delay ! delayed nu + integer :: nu ! actual number of neutrons produced + real(8) :: mu ! fission neutron angular cosine + real(8) :: phi ! fission neutron azimuthal angle + real(8) :: beta ! delayed neutron fraction + character(250) :: msg ! error message + + ! copy energy of neutron + E = p % E + + ! ================================================================ + ! DETERMINE TOTAL NU + if (table % nu_t_type == NU_NONE) then + msg = "No neutron emission data for table: " // table % name + call error(msg) + elseif (table % nu_t_type == NU_POLYNOMIAL) then + ! determine number of coefficients + NC = int(table % nu_t_data(1)) + + ! sum up polynomial in energy + nu_total = ZERO + do i = 0, NC - 1 + c = table % nu_t_data(i+2) + nu_total = nu_total + c * E**i + end do + elseif (table % nu_t_type == NU_TABULAR) then + ! determine number of interpolation regions -- as far as I can + ! tell, no nu data has multiple interpolation + ! regions. Furthermore, it seems all are lin-lin. + NR = int(table % nu_t_data(1)) + if (NR /= 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to determine total nu." + call error(msg) + end if + + ! determine number of energies + loc = 2 + 2*NR + NE = int(table % nu_t_data(loc)) + + ! do binary search over tabuled energies to determine + ! appropriate index and interpolation factor + j = binary_search(table % nu_t_data(loc+1), NE, E) + f = (E - table % nu_t_data(loc+j)) / & + & (table % nu_t_data(loc+j+1) - table % nu_t_data(loc+j)) + + ! determine nu total + loc = loc + NE + nu_total = table % nu_t_data(loc+j) + f * & + & (table % nu_t_data(loc+j+1) - table % nu_t_data(loc+j)) + end if + + ! ================================================================ + ! DETERMINE PROMPT NU + if (table % nu_p_type == NU_NONE) then + ! since no prompt or delayed data is present, this means all + ! neutron emission is prompt + nu_prompt = nu_total + elseif (table % nu_p_type == NU_POLYNOMIAL) then + ! determine number of coefficients + NC = int(table % nu_p_data(1)) + + ! sum up polynomial in energy + nu_prompt = ZERO + do i = 0, NC - 1 + c = table % nu_p_data(i+2) + nu_prompt = nu_prompt + c * E**i + end do + elseif (table % nu_p_type == NU_TABULAR) then + ! determine number of interpolation regions + NR = int(table % nu_p_data(1)) + if (NR /= 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to determine prompt nu." + call error(msg) + end if + + ! determine number of energies + loc = 2 + 2*NR + NE = int(table % nu_p_data(loc)) + + ! do binary search over tabuled energies to determine + ! appropriate index and interpolation factor + j = binary_search(table % nu_p_data(loc+1), NE, E) + f = (E - table % nu_p_data(loc+j)) / & + & (table % nu_p_data(loc+j+1) - table % nu_p_data(loc+j)) + + ! determine nu total + loc = loc + NE + nu_prompt = table % nu_p_data(loc+j) + f * & + & (table % nu_p_data(loc+j+1) - table % nu_p_data(loc+j)) + end if + + ! ================================================================ + ! DETERMINE DELAYED NU + if (table % nu_d_type == NU_NONE) then + nu_delay = ZERO + elseif (table % nu_d_type == NU_TABULAR) then + ! determine number of interpolation regions + NR = int(table % nu_d_data(1)) + if (NR /= 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to determine delayed nu." + call error(msg) + end if + + ! determine number of energies + loc = 2 + 2*NR + NE = int(table % nu_d_data(loc)) + + ! do binary search over tabuled energies to determine + ! appropriate index and interpolation factor + j = binary_search(table % nu_d_data(loc+1), NE, E) + f = (E - table % nu_d_data(loc+j)) / & + & (table % nu_d_data(loc+j+1) - table % nu_d_data(loc+j)) + + ! determine nu total + loc = loc + NE + nu_delay = table % nu_d_data(loc+j) + f * & + & (table % nu_d_data(loc+j+1) - table % nu_d_data(loc+j)) + end if + + beta = nu_delay / nu_total + + ! TODO: Heat generation from fission + + ! Sample number of neutrons produced + nu_total = p % wgt * nu_total + if (rang() > nu_total - int(nu_total)) then + nu = int(nu_total) + else + nu = int(nu_total) + 1 + end if + + ! Bank source neutrons + if (nu == 0 .or. n_bank == 3*n_particles) return + do i = n_bank + 1, min(n_bank + nu, 3*n_particles) + ! Bank source neutrons by copying particle data + fission_bank(i) % uid = p % uid + fission_bank(i) % xyz = p % xyz + + ! sample cosine of angle + if (rxn % has_angle_dist) then + mu = sample_angle(rxn, E) + else + mu = TWO * rang() - ONE + end if + + ! Sample azimuthal angle uniformly in [0,2*pi) + phi = TWO*PI*rang() + fission_bank(i) % uvw(1) = mu + fission_bank(i) % uvw(2) = sqrt(1. - mu**2) * cos(phi) + fission_bank(i) % uvw(3) = sqrt(1. - mu**2) * sin(phi) + + ! determine energy of fission neutron + call sample_energy(rxn, E, E_out, mu_out) + fission_bank(i) % E = E_out + end do + + ! increment number of bank sites + n_bank = min(n_bank + nu, 3*n_particles) + + ! kill original neutron + p % alive = .false. + + end subroutine n_fission + +!===================================================================== +! LEVEL_SCATTER +!===================================================================== + + subroutine level_scatter(p, table, rxn) + + type(Particle), pointer :: p + type(AceContinuous), pointer :: table + type(AceReaction), pointer :: rxn + + real(8) :: A ! atomic weight ratio of nuclide + real(8) :: E_in ! incoming energy + real(8) :: mu_cm ! cosine of scattering angle in center-of-mass + real(8) :: mu_lab ! cosine of scattering angle in laboratory + real(8) :: E_cm ! outgoing energy in center-of-mass + real(8) :: E_lab ! outgoing energy in laboratory + character(250) :: msg ! error message + + + ! copy energy of neutron + E_in = p % E + + ! determine A + A = table % awr + + ! determine if scattering is in CM (it should be!) + if (rxn % TY < 0) then + ! scattering angle in center-of-mass + mu_cm = sample_angle(rxn, E_in) + else + msg = "Level inelastic scattering should not sample angle & + &in laboratory system!" + call error(msg) + end if + + ! sample outgoing energy in center-of-mass + call sample_energy(rxn, E_in, E_cm) + + ! determine outgoing energy in lab + E_lab = E_cm + (E_in + TWO * mu_cm * (A+ONE) * sqrt(E_in * E_cm)) & + & / ((A+ONE)*(A+ONE)) + + ! determine outgoing angle in lab + mu_lab = mu_cm * sqrt(E_cm/E_lab) + ONE/(A+ONE) * sqrt(E_in/E_lab) + + ! change direction of particle + call rotate_angle(p, mu_lab) + + ! change energy of particle + p % E = E_lab + + ! find energy index, interpolation factor + call find_energy_index(p) + + end subroutine level_scatter !===================================================================== ! N_GAMMA @@ -321,7 +566,7 @@ contains ! a tabular distribution. !===================================================================== - subroutine sample_angle(rxn, E) + function sample_angle(rxn, E) result(mu) type(AceReaction), pointer :: rxn ! reaction real(8), intent(in) :: E ! incoming energy @@ -332,81 +577,86 @@ contains integer :: i ! incoming energy bin integer :: n ! number of incoming energy bins integer :: loc ! location in data array - integer :: np ! number of points in cos distribution - integer :: bin ! cosine bin - real(8) :: f ! interpolation factor - real(8) :: mu0 ! cosine in bin b - real(8) :: mu1 ! cosine in bin b+1 + integer :: NP ! number of points in cos distribution + integer :: k ! index on cosine grid + real(8) :: r ! interpolation factor on incoming energy + real(8) :: frac ! interpolation fraction on cosine + real(8) :: mu0 ! cosine in bin k + real(8) :: mu1 ! cosine in bin k+1 real(8) :: mu ! final cosine sampled - real(8) :: c ! cumulative distribution frequency + real(8) :: c_k ! cumulative frequency at k + real(8) :: c_k1 ! cumulative frequency at k+1 real(8) :: p0,p1 ! probability distribution character(250) :: msg ! error message ! determine number of incoming energies - n = rxn % adist_n_energy + n = rxn % adist % n_energy ! find energy bin and calculate interpolation factor -- if the ! energy is outside the range of the tabulated energies, choose ! the first or last bins - if (E < rxn%adist_energy(1)) then + if (E < rxn % adist % energy(1)) then i = 1 - f = 0.0 - elseif (E > rxn%adist_energy(n)) then + r = 0.0 + elseif (E > rxn % adist % energy(n)) then i = n - 1 - f = 1.0 + r = 1.0 else - i = binary_search(rxn%adist_energy, n, E) - f = (E - rxn % adist_energy(i)) / & - & (rxn % adist_energy(i+1) - rxn % adist_energy(i)) + i = binary_search(rxn % adist % energy, n, E) + r = (E - rxn % adist % energy(i)) / & + & (rxn % adist % energy(i+1) - rxn % adist % energy(i)) end if ! Sample between the ith and (i+1)th bin - if (f > rang()) i = i + 1 + if (r > rang()) i = i + 1 ! check whether this is a 32-equiprobable bin or a tabular ! distribution - loc = rxn % adist_location(i) - type = rxn % adist_type(i) + loc = rxn % adist % location(i) + type = rxn % adist % type(i) if (type == ANGLE_ISOTROPIC) then - mu = 2.0_8 * rang() - 1 + mu = TWO * rang() - ONE elseif (type == ANGLE_32_EQUI) then ! sample cosine bin xi = rang() - bin = 1 + int(32.0_8*xi) + k = 1 + int(32.0_8*xi) ! calculate cosine - mu0 = rxn % adist_data(loc + bin - 1) - mu1 = rxn % adist_data(loc + bin) - mu = mu0 + (32.0_8 * xi - bin) * (mu1 - mu0) + mu0 = rxn % adist % data(loc + k) + mu1 = rxn % adist % data(loc + k+1) + mu = mu0 + (32.0_8 * xi - k) * (mu1 - mu0) elseif (type == ANGLE_TABULAR) then - interp = rxn % adist_data(loc) - np = rxn % adist_data(loc+1) + interp = rxn % adist % data(loc + 1) + NP = rxn % adist % data(loc + 2) ! determine outgoing cosine bin xi = rang() - do bin = loc+2, loc+1+NP - c = rxn % adist_data(bin+2*np) - if (xi > c) exit + loc = loc + 2 + c_k = rxn % adist % data(loc + 2*NP + 1) + do k = 1, NP-1 + c_k1 = rxn % adist % data(loc + 2*NP + k+1) + if (xi < c_k1) exit + c_k = c_k1 end do - p0 = rxn % adist_data(bin+np) - mu0 = rxn % adist_data(bin) + p0 = rxn % adist % data(loc + NP + k) + mu0 = rxn % adist % data(loc + k) if (interp == HISTOGRAM) then ! Histogram interpolation - mu = mu0 + (xi - c)/p0 + mu = mu0 + (xi - c_k)/p0 elseif (interp == LINEARLINEAR) then ! Linear-linear interpolation -- not sure how you come about ! the formula given in the MCNP manual - p1 = rxn % adist_data(bin+np+1) - mu1 = rxn % adist_data(bin+1) - - f = (p1 - p0)/(mu1 - mu0) - if (f == ZERO) then - mu = mu0 + (xi - c)/p0 + p1 = rxn % adist % data(loc + NP + k+1) + mu1 = rxn % adist % data(loc + k+1) + + frac = (p1 - p0)/(mu1 - mu0) + if (frac == ZERO) then + mu = mu0 + (xi - c_k)/p0 else - mu = mu0 + (sqrt(p0*p0 + 2*f*(xi - c))-p0)/f + mu = mu0 + (sqrt(p0*p0 + 2*frac*(xi - c_k))-p0)/frac end if else msg = "Unknown interpolation type: " // trim(int_to_str(interp)) @@ -418,7 +668,7 @@ contains call error(msg) end if - end subroutine sample_angle + end function sample_angle !===================================================================== ! ROTATE_ANGLE rotates direction cosines through a polar angle whose @@ -427,46 +677,572 @@ contains ! done in MCNP and SERPENT. !===================================================================== - subroutine rotate_angle(mu, u, v, w) + subroutine rotate_angle(p, mu) - real(8), intent(in) :: mu ! cosine of angle - real(8), intent(inout) :: u - real(8), intent(inout) :: v - real(8), intent(inout) :: w + type(Particle), pointer :: p + real(8), intent(in) :: mu ! cosine of angle in lab real(8) :: phi, sinphi, cosphi real(8) :: a,b real(8) :: u0, v0, w0 ! Copy original directional cosines - u0 = u - v0 = v - w0 = w + u0 = p % uvw(1) + v0 = p % uvw(2) + w0 = p % uvw(3) ! Sample azimuthal angle in [0,2pi) - phi = 2.0_8 * PI * rang() + phi = TWO * PI * rang() ! Precompute factors to save flops sinphi = sin(phi) cosphi = cos(phi) - a = sqrt(1.0_8 - mu*mu) - b = sqrt(1.0_8 - w*w) + a = sqrt(ONE - mu*mu) + b = sqrt(ONE - w0*w0) ! Need to treat special case where sqrt(1 - w**2) is close to zero ! by expanding about the v component rather than the w component if (b > 1e-10) then - u = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b - v = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b - w = mu*w0 - a*b*cosphi + p % uvw(1) = mu*u0 + a*(u0*w0*cosphi - v0*sinphi)/b + p % uvw(2) = mu*v0 + a*(v0*w0*cosphi + u0*sinphi)/b + p % uvw(3) = mu*w0 - a*b*cosphi else - b = sqrt(1.0_8 - v*v) - u = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b - v = mu*v0 - a*b*cosphi - w = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b + b = sqrt(ONE - v0*v0) + p % uvw(1) = mu*u0 + a*(u0*v0*cosphi + w0*sinphi)/b + p % uvw(2) = mu*v0 - a*b*cosphi + p % uvw(3) = mu*w0 + a*(v0*w0*cosphi - u0*sinphi)/b end if end subroutine rotate_angle +!===================================================================== +! SAMPLE_ENERGY +!===================================================================== + + subroutine sample_energy(rxn, E_in, E_out, mu_out) + + type(AceReaction), pointer :: rxn + real(8), intent(in) :: E_in + real(8), intent(out) :: E_out + real(8), intent(inout), optional :: mu_out + + integer :: i ! index on incoming energy grid + integer :: k ! sampled index on outgoing grid + integer :: l ! sampled index on incoming grid + integer :: loc ! dummy index + integer :: NR ! number of interpolation regions + integer :: NE ! number of energies + integer :: NET ! number of outgoing energies + integer :: INTTp ! combination of INTT and ND + integer :: INTT ! 1 = histogram, 2 = linear-linear + integer :: ND ! number of discrete lines + integer :: NP ! number of points in distribution + + real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i + real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1 + real(8) :: E_1, E_K ! endpoints interpolated between i and i+1 + + real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l + real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l + real(8) :: c_k, c_k1 ! cumulative probability + + real(8) :: KM_A ! Kalbach-Mann parameter R + real(8) :: KM_R ! Kalbach-Mann parameter R + real(8) :: A_k, A_k1 ! Kalbach-Mann A on outgoing grid l + real(8) :: R_k, R_k1 ! Kalbach-Mann R on outgoing grid l + + real(8) :: Watt_a, Watt_b ! Watt spectrum parameters + + real(8) :: E_cm + real(8) :: xi1, xi2, xi3, xi4 + real(8) :: r ! interpolation factor on incoming energy + real(8) :: frac ! interpolation factor on outgoing energy + real(8) :: U ! restriction energy + real(8) :: T ! nuclear temperature + character(250) :: msg ! error message + + ! TODO: If there are multiple scattering laws, sample scattering + ! law + + ! Check for multiple interpolation regions + if (rxn % edist % n_interp > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sampling secondary energy distribution." + call error(msg) + end if + + ! Determine which secondary energy distribution law to use + select case (rxn % edist % law) + case (1) + ! ============================================================= + ! TABULAR EQUIPROBABLE ENERGY BINS + + ! read number of interpolation regions, incoming energies, and + ! outgoing energies + NR = rxn % edist % data(1) + NE = rxn % edist % data(2 + 2*NR) + NET = rxn % edist % data(3 + 2*NR + NE) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample equiprobable energy bins." + call error(msg) + end if + + ! determine index on incoming energy grid and interpolation + ! factor + loc = 2 + 2*NR + i = binary_search(rxn % edist % data(loc+1), NE, E_in) + r = (E_in - rxn%edist%data(loc+i)) / & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + + ! Sample outgoing energy bin + xi1 = rang() + k = 1 + int(NET * xi1) + + ! Randomly select between the outgoing table for incoming + ! energy E_i and E_(i+1) + if (rang() < r) then + l = i + 1 + else + l = i + end if + + loc = 3 + 2*NR + NE + (l-1)*NET + E_l_k = rxn % edist % data(loc+k) + E_l_k1 = rxn % edist % data(loc+k+1) + xi2 = rang() + E_out = E_l_k + xi2*(E_l_k1 - E_l_k) + + ! TODO: Add scaled interpolation + + case (3) + ! ============================================================= + ! INELASTIC LEVEL SCATTERING + + E_cm = rxn%edist%data(2) * (E_in - rxn%edist%data(1)) + + E_out = E_cm + + case (4) + ! ============================================================= + ! CONTINUOUS TABULAR DISTRIBUTION + + ! read number of interpolation regions and incoming energies + NR = rxn % edist % data(1) + NE = rxn % edist % data(2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample continuous tabular distribution." + call error(msg) + end if + + ! find energy bin and calculate interpolation factor -- if the + ! energy is outside the range of the tabulated energies, choose + ! the first or last bins + loc = 2 + 2*NR + if (E_in < rxn % edist % data(loc+1)) then + i = 1 + r = 0.0 + elseif (E_in > rxn % edist % energy(loc+NE)) then + i = NE - 1 + r = 1.0 + else + i = binary_search(rxn % edist % energy(loc+1), NE, E_in) + r = (E_in - rxn%edist%energy(loc+i)) / & + & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + end if + + ! Sample between the ith and (i+1)th bin + xi2 = rang() + if (r > xi2) then + l = i + 1 + else + l = i + end if + + ! interpolation on grid i for energy E1 + loc = rxn % edist % data(2 + 2*NR + NE + i) + 2 ! start of EOUT(i) + E_i_1 = rxn%edist%data(loc + 1) + E_i_K = rxn%edist%data(loc + NE) + + loc = rxn % edist % data(2 + 2*NR + NE + i + 1) + 2 ! start of EOUT(i+1) + E_i1_1 = rxn%edist%data(loc + 1) + E_i1_K = rxn%edist%data(loc + NE) + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + + ! determine location of outgoing energies, pdf, cdf for E(l) + loc = rxn % edist % data(2 + 2*NR + NE + l) + + ! determine type of interpolation and number of discrete lines + INTTp = rxn % edist % data(loc + 1) + NP = rxn % edist % data(loc + 2) + if (INTTp > 10) then + INTT = mod(INTTp,10) + ND = (INTTp - INTT)/10 + else + INTT = INTTp + ND = 0 + end if + + if (ND > 0) then + ! discrete lines present + msg = "Discrete lines in continuous tabular distributed not & + &yet supported" + call error(msg) + end if + + ! determine outgoing energy bin + xi1 = rang() + loc = loc + 2 ! start of EOUT + c_k = rxn % edist % data(loc + 2*NP + 1) + do k = 1, NP-1 + c_k1 = rxn % edist % data(loc + 2*NP + k+1) + if (xi1 < c_k1) exit + c_k = c_k1 + end do + + E_l_k = rxn % edist % data(loc+k) + p_l_k = rxn % edist % data(loc+NP+k) + if (INTT == HISTOGRAM) then + ! Histogram interpolation + E_out = E_l_k + (xi1 - c_k)/p_l_k + + elseif (INTT == LINEARLINEAR) then + ! Linear-linear interpolation -- not sure how you come about + ! the formula given in the MCNP manual + E_l_k1 = rxn % edist % data(loc+k+1) + p_l_k1 = rxn % edist % data(loc+NP+k+1) + + frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) + if (frac == ZERO) then + E_out = E_l_k + (xi1 - c_k)/p_l_k + else + E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(xi1 - c_k)) - & + & p_l_k)/frac + end if + else + msg = "Unknown interpolation type: " // trim(int_to_str(INTT)) + call error(msg) + end if + + ! Now interpolate between incident enregy bins i and i + 1 + if (l == i) then + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + end if + + case (5) + ! ============================================================= + ! GENERAL EVAPORATION SPECTRUM + + case (7) + ! ============================================================= + ! MAXWELL FISSION SPECTRUM + + ! read number of interpolation regions and incoming energies + NR = rxn % edist % data(1) + NE = rxn % edist % data(2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample Maxwell fission spectrum." + call error(msg) + end if + + ! find incident energy bin and calculate interpolation factor + loc = 2 + 2*NR + if (E_in < rxn % edist % data(loc+1)) then + i = 1 + r = 0.0 + elseif (E_in > rxn % edist % energy(loc+NE)) then + i = NE - 1 + r = 1.0 + else + i = binary_search(rxn % edist % energy(loc+1), NE, E_in) + r = (E_in - rxn%edist%energy(loc+i)) / & + & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + end if + + ! determine nuclear temperature from tabulated function + loc = loc + NE + T = rxn%edist%data(loc+i) + r * & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + + ! sample maxwell fission spectrum + E_out = maxwell_spectrum(T) + + case (9) + ! ============================================================= + ! EVAPORATION SPECTRUM + + ! read number of interpolation regions and incoming energies + NR = rxn % edist % data(1) + NE = rxn % edist % data(2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample evaporation spectrum." + call error(msg) + end if + + ! find energy bin and calculate interpolation factor -- if the + ! energy is outside the range of the tabulated energies, choose + ! the first or last bins + loc = 2 + 2*NR + if (E_in < rxn % edist % data(loc+1)) then + i = 1 + r = 0.0 + elseif (E_in > rxn % edist % energy(loc+NE)) then + i = NE - 1 + r = 1.0 + else + i = binary_search(rxn % edist % energy(loc+1), NE, E_in) + r = (E_in - rxn%edist%energy(loc+i)) / & + & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + end if + + ! determine nuclear temperature from tabulated function + loc = loc + NE + T = rxn%edist%data(loc+i) + r * & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + + ! sample outgoing energy based on evaporation spectrum + ! probability density function + do + xi1 = rang() + xi2 = rang() + E_out = -T * log(xi1*xi2) + if (E_out <= E_in - U) exit + end do + + case (11) + ! ============================================================= + ! ENERGY-DEPENDENT WATT SPECTRUM + + ! read number of interpolation regions and incoming energies + ! for parameter 'a' + NR = rxn % edist % data(1) + NE = rxn % edist % data(2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample Watt fission spectrum." + call error(msg) + end if + + ! find incident energy bin and calculate interpolation factor + loc = 2 + 2*NR + if (E_in < rxn % edist % data(loc+1)) then + i = 1 + r = 0.0 + elseif (E_in > rxn % edist % energy(loc+NE)) then + i = NE - 1 + r = 1.0 + else + i = binary_search(rxn % edist % energy(loc+1), NE, E_in) + r = (E_in - rxn%edist%energy(loc+i)) / & + & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + end if + + ! determine Watt parameter 'a' from tabulated function + loc = loc + NE + Watt_a = rxn%edist%data(loc+i) + r * & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + + ! read number of interpolation regions and incoming energies + ! for parameter 'b' + loc = loc + NE + NR = rxn % edist % data(loc + 1) + NE = rxn % edist % data(loc + 2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample Watt fission spectrum." + call error(msg) + end if + + ! find incident energy bin and calculate interpolation factor + loc = loc + 2 + 2*NR + if (E_in < rxn % edist % data(loc+1)) then + i = 1 + r = 0.0 + elseif (E_in > rxn % edist % energy(loc+NE)) then + i = NE - 1 + r = 1.0 + else + i = binary_search(rxn % edist % energy(loc+1), NE, E_in) + r = (E_in - rxn%edist%energy(loc+i)) / & + & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + end if + + ! determine Watt parameter 'b' from tabulated function + loc = loc + NE + Watt_b = rxn%edist%data(loc+i) + r * & + & (rxn%edist%data(loc+i+1) - rxn%edist%data(loc+i)) + + ! Sample energy-dependent Watt fission spectrum + E_out = watt_spectrum(Watt_a, Watt_b) + + case (44) + ! ============================================================= + ! KALBACH-MANN CORRELATED SCATTERING + + if (.not. present(mu_out)) then + msg = "Law 44 called without giving mu_out as argument." + call error(msg) + end if + + ! read number of interpolation regions and incoming energies + NR = rxn % edist % data(1) + NE = rxn % edist % data(2 + 2*NR) + if (NR > 0) then + msg = "Multiple interpolation regions not supported while & + &attempting to sample Kalbach-Mann distribution." + call error(msg) + end if + + ! find energy bin and calculate interpolation factor -- if the + ! energy is outside the range of the tabulated energies, choose + ! the first or last bins + loc = 2 + 2*NR + if (E_in < rxn % edist % data(loc+1)) then + i = 1 + r = 0.0 + elseif (E_in > rxn % edist % energy(loc+NE)) then + i = NE - 1 + r = 1.0 + else + i = binary_search(rxn % edist % energy(loc+1), NE, E_in) + r = (E_in - rxn%edist%energy(loc+i)) / & + & (rxn%edist%energy(loc+i+1) - rxn%edist%energy(loc+i)) + end if + + ! Sample between the ith and (i+1)th bin + xi2 = rang() + if (r > xi2) then + l = i + 1 + else + l = i + end if + + ! interpolation on grid i for energy E1 + loc = rxn % edist % data(2 + 2*NR + NE + i) + 2 ! start of EOUT(i) + E_i_1 = rxn%edist%data(loc + 1) + E_i_K = rxn%edist%data(loc + NE) + + loc = rxn % edist % data(2 + 2*NR + NE + i + 1) + 2 ! start of EOUT(i+1) + E_i1_1 = rxn%edist%data(loc + 1) + E_i1_K = rxn%edist%data(loc + NE) + + E_1 = E_i_1 + r*(E_i1_1 - E_i_1) + E_K = E_i_K + r*(E_i1_K - E_i_K) + + ! determine location of outgoing energies, pdf, cdf for E(l) + loc = rxn % edist % data(2 + 2*NR + NE + l) + + ! determine type of interpolation and number of discrete lines + INTTp = rxn % edist % data(loc + 1) + NP = rxn % edist % data(loc + 2) + if (INTTp > 10) then + INTT = mod(INTTp,10) + ND = (INTTp - INTT)/10 + else + INTT = INTTp + ND = 0 + end if + + if (ND > 0) then + ! discrete lines present + msg = "Discrete lines in continuous tabular distributed not & + &yet supported" + call error(msg) + end if + + ! determine outgoing energy bin + xi1 = rang() + loc = loc + 2 ! start of EOUT + c_k = rxn % edist % data(loc + 2*NP + 1) + do k = 1, NP-1 + c_k1 = rxn % edist % data(loc + 2*NP + k+1) + if (xi1 < c_k1) exit + c_k = c_k1 + end do + + E_l_k = rxn % edist % data(loc+k) + p_l_k = rxn % edist % data(loc+NP+k) + if (INTT == HISTOGRAM) then + ! Histogram interpolation + E_out = E_l_k + (xi1 - c_k)/p_l_k + + ! Determine Kalbach-Mann parameters + KM_R = rxn % edist % data(loc + 3*NP + k) + KM_A = rxn % edist % data(loc + 4*NP + k) + + elseif (INTT == LINEARLINEAR) then + ! Linear-linear interpolation -- not sure how you come about + ! the formula given in the MCNP manual + E_l_k1 = rxn % edist % data(loc+k+1) + p_l_k1 = rxn % edist % data(loc+NP+k+1) + + ! Find E prime + frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k) + if (frac == ZERO) then + E_out = E_l_k + (xi1 - c_k)/p_l_k + else + E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(xi1 - c_k)) - & + & p_l_k)/frac + end if + + ! Determine Kalbach-Mann parameters + R_k = rxn % edist % data(loc + 3*NP + k) + R_k1 = rxn % edist % data(loc + 3*NP + k+1) + A_k = rxn % edist % data(loc + 4*NP + k) + A_k1 = rxn % edist % data(loc + 4*NP + k+1) + + KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) + KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k) + else + msg = "Unknown interpolation type: " // trim(int_to_str(INTT)) + call error(msg) + end if + + ! Now interpolate between incident enregy bins i and i + 1 + if (l == i) then + E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1) + else + E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1) + end if + + ! Sampled correlated angle from Kalbach-Mann parameters + xi3 = rang() + xi4 = rang() + T = (TWO*xi4 - ONE) * sinh(KM_A) + if (xi3 > KM_R) then + mu_out = log(T + sqrt(T*T + ONE))/KM_A + else + mu_out = log(xi4*exp(KM_A) + (ONE - xi4)*exp(-KM_A))/KM_A + end if + + case (61) + ! ============================================================= + ! CORRELATED ENERGY AND ANGLE DISTRIBUTION + + if (.not. present(mu_out)) then + msg = "Law 44 called without giving mu_out as argument." + call error(msg) + end if + + case (66) + ! ============================================================= + ! N-BODY PHASE SPACE DISTRIBUTION + + case (67) + ! ============================================================= + ! LABORATORY ENERGY-ANGLE LAW + + end select + + end subroutine sample_energy + !===================================================================== ! MAXWELL_SPECTRUM samples an energy from the Maxwell fission ! distribution based on a rejection sampling scheme. This is described @@ -485,7 +1261,7 @@ contains do r2 = rang() d = r1*r1 + r2*r2 - if (d < 1) exit + if (d < ONE) exit r1 = r2 end do diff --git a/src/source.f90 b/src/source.f90 index 81b382f56..15c6c3386 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -20,6 +20,7 @@ contains real(8) :: phi ! azimuthal angle real(8) :: mu ! cosine of polar angle real(8) :: p_min(3), p_max(3) + type(Particle), pointer :: p => null() ! Allocate fission and source banks allocate( source_bank(n_particles) ) @@ -34,24 +35,28 @@ contains ! Initialize first cycle source bank do i = 1, n_particles call RN_init_particle(int(i,8)) + p => source_bank(i) ! position r = (/ (rang(), j = 1,3) /) - source_bank(i)%uid = i - source_bank(i)%xyz = p_min + r*(p_max - p_min) + p % uid = i + p % xyz = p_min + r*(p_max - p_min) ! angle - phi = 2.0_8*PI*rang() - mu = 2.0_8*rang() - 1.0_8 - source_bank(i)%uvw(1) = mu - source_bank(i)%uvw(2) = sqrt(1. - mu**2) * cos(phi) - source_bank(i)%uvw(3) = sqrt(1. - mu**2) * sin(phi) + phi = TWO*PI*rang() + mu = TWO*rang() - ONE + p % uvw(1) = mu + p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) + p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) - ! cell - source_bank(i)%cell = 0 + ! set defaults + p % type = NEUTRON + p % cell = 0 + p % surface = 0 + p % universe = 0 + p % wgt = ONE + p % alive = .true. - ! set particle to be alive - source_bank(i)%alive = .true. end do ! Reset source index @@ -81,4 +86,28 @@ contains end function get_source_particle +!===================================================================== +! ADD_BANK_SITES +!===================================================================== + + subroutine add_bank_sites(p, table, n) + + type(Particle), pointer :: p + type(AceContinuous), pointer :: table + integer, intent(in) :: n + + integer :: i + + if (n == 0 .or. n_bank == 3*n_particles) return + do i = n_bank + 1, min(n_bank + n, 3*n_particles) + ! Copy particle data + fission_bank(i)%uid = p%uid + fission_bank(i)%xyz = p%xyz + + ! TODO: Sample angle and energy from secondary distribution + end do + n_bank = min(n_bank + n, 3*n_particles) + + end subroutine add_bank_sites + end module source diff --git a/src/types.f90 b/src/types.f90 index 14707d824..bedeb27cf 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -71,14 +71,15 @@ module types integer :: type ! Particle type (n, p, e, etc) real(8) :: xyz(3) ! location real(8) :: uvw(3) ! directional cosines + real(8) :: wgt ! particle weight real(8) :: E ! energy integer :: IE ! index on energy grid real(8) :: interp ! interpolation factor for energy grid integer :: cell ! current cell integer :: universe ! current universe integer :: surface ! current surface - real(8) :: wgt ! particle weight logical :: alive ! is particle alive? + integer :: n_coll ! # of collisions end type Particle !===================================================================== @@ -89,7 +90,9 @@ module types type Bank integer :: uid ! Unique ID - real(8) :: xyz(3) ! Location of bank particle + real(8) :: xyz(3) ! location of bank particle + real(8) :: uvw(3) ! diretional cosines + real(8) :: E ! energy end type Bank !===================================================================== @@ -133,6 +136,35 @@ module types real(8), allocatable :: values(:) ! values for particular source type end type ExtSource +!===================================================================== +! ACEDISTANGLE contains data for a tabular secondary angle +! distribution whether it be tabular or 32 equiprobable cosine bins +!===================================================================== + + type AceDistAngle + integer :: n_energy ! # of incoming energies + real(8), allocatable :: energy(:) ! incoming energy grid + integer, allocatable :: type(:) ! type of distribution + integer, allocatable :: location(:) ! location of each table + real(8), allocatable :: data(:) ! angular distribution data + end type AceDistAngle + +!===================================================================== +! ACEDISTENERGY contains data for a secondary energy distribution for +! all scattering laws +!===================================================================== + + type AceDistEnergy + integer :: law ! secondary distribution law + integer :: n_interp ! # of interpolation regions + integer, allocatable :: nbt(:) ! ENDF interpolation parameters + integer, allocatable :: int(:) ! '' + integer :: n_energy ! # of incoming energies + real(8), allocatable :: energy(:) ! energy grid for law validity + real(8), allocatable :: pvalid(:) ! probability of law validity + real(8), allocatable :: data(:) ! energy distribution data + end type AceDistEnergy + !===================================================================== ! ACEREACTION contains the cross-section and secondary energy and ! angle distributions for a single reaction in a continuous-energy @@ -142,28 +174,13 @@ module types type AceReaction integer :: MT ! ENDF MT value real(8) :: Q_value ! Reaction Q value - real(8) :: TY ! Number of neutrons released + integer :: TY ! Number of neutrons released integer :: IE ! Starting energy grid index real(8), allocatable :: sigma(:) ! Cross section values - logical :: has_angle_dist - logical :: has_energy_dist - - ! Secondary angle distribution - integer :: adist_n_energy ! # of incoming energies - real(8), allocatable :: adist_energy(:) ! incoming energy grid - integer, allocatable :: adist_type(:) ! type of distribution - integer, allocatable :: adist_location(:) ! location of each table - real(8), allocatable :: adist_data(:) ! angular distribution data - - ! Secondary energy distribution - integer :: edist_law ! secondary distribution law - integer :: edist_n_interp ! # of interpolation regions - integer, allocatable :: edist_nbt(:) ! ENDF interpolation parameters - integer, allocatable :: edist_int(:) ! '' - integer :: edist_n_energy ! # of incoming energies - real(8), allocatable :: edist_energy(:) ! energy grid for law validity - real(8), allocatable :: edist_pvalid(:) ! probability of law validity - real(8), allocatable :: edist_data(:) ! energy distribution data + logical :: has_angle_dist ! Angle distribution present? + logical :: has_energy_dist ! Energy distribution present? + type(AceDistAngle) :: adist ! Secondary angular distribution + type(AceDistEnergy) :: edist ! Secondary energy distribution end type AceReaction !===================================================================== @@ -184,19 +201,23 @@ module types real(8), allocatable :: sigma_el(:) real(8), allocatable :: heating(:) - integer :: nu_p_type - real(8), allocatable :: nu_p_energy(:) - real(8), allocatable :: nu_p_value(:) - + ! Total fission neutron emission integer :: nu_t_type - real(8), allocatable :: nu_t_energy(:) - real(8), allocatable :: nu_t_value(:) + real(8), allocatable :: nu_t_data(:) - real(8), allocatable :: nu_d_energy(:) - real(8), allocatable :: nu_d_value(:) + ! Prompt fission neutron emission + integer :: nu_p_type + real(8), allocatable :: nu_p_data(:) + + ! Delayed fission neutron emission + integer :: nu_d_type + real(8), allocatable :: nu_d_data(:) + type(AceDistEnergy) :: nu_d_edist real(8), allocatable :: nu_d_precursor_const(:,:) real(8), allocatable :: nu_d_precursor_energy(:,:) real(8), allocatable :: nu_d_precursor_prob(:,:) + + ! Reactions integer :: n_reaction type(AceReaction), pointer :: reactions(:) => null() @@ -220,27 +241,6 @@ module types real(8), allocatable :: elastic_mu_out(:) end type AceThermal -!===================================================================== -! ACEDISTANGLE contains data for a tabular secondary angle -! distribution whether it be tabular or 32 equiprobable cosine bins -!===================================================================== - -!!$ type AceDistAngle -!!$ integer :: n_energy -!!$ real(8), allocatable :: energy -!!$ integer, allocatable :: location -!!$ real(8), allocatable :: data -!!$ end type AceDistAngle - -!===================================================================== -! ACEDISTTAB contains data for a tabular secondary energy distribution -! according to MCNP Law 4 (ENDF Law 1) -!===================================================================== - - type AceDistEnergy - integer :: law - end type AceDistEnergy - !===================================================================== ! XSDATA contains data read in from a SERPENT xsdata file !=====================================================================