diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index c8237effb1..fb41b5b055 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -1,3 +1,18 @@ +ace.o: ace_header.o +ace.o: constants.o +ace.o: datatypes.o +ace.o: datatypes_header.o +ace.o: endf.o +ace.o: error.o +ace.o: fission.o +ace.o: global.o +ace.o: material_header.o +ace.o: output.o +ace.o: string.o + +ace_header.o: constants.o +ace_header.o: endf_header.o + cmfd_execute.o: cmfd_utils.o cmfd_execute.o: global.o cmfd_execute.o: mesh.o @@ -13,21 +28,15 @@ cmfd_utils.o: mesh_header.o cmfd_utils.o: string.o cmfd_utils.o: xml-fortran/templates/cmfd_t.o +cross_section.o: ace_header.o cross_section.o: constants.o -cross_section.o: cross_section_header.o -cross_section.o: datatypes.o -cross_section.o: datatypes_header.o -cross_section.o: endf.o cross_section.o: error.o -cross_section.o: fileio.o cross_section.o: fission.o -cross_section.o: global.o +cross_section.o: global.o cross_section.o: material_header.o -cross_section.o: output.o -cross_section.o: string.o - -cross_section_header.o: constants.o -cross_section_header.o: endf_header.o +cross_section.o: particle_header.o +cross_section.o: random_lcg.o +cross_section.o: search.o datatypes.o: datatypes_header.o @@ -44,12 +53,8 @@ energy_grid.o: datatypes_header.o energy_grid.o: global.o energy_grid.o: output.o -fileio.o: constants.o -fileio.o: global.o -fileio.o: string.o - +fission.o: ace_header.o fission.o: constants.o -fission.o: cross_section_header.o fission.o: error.o fission.o: global.o fission.o: interpolation.o @@ -65,10 +70,10 @@ geometry.o: particle_header.o geometry.o: string.o geometry.o: tally.o +global.o: ace_header.o global.o: bank_header.o global.o: cmfd_header.o global.o: constants.o -global.o: cross_section_header.o global.o: datatypes_header.o global.o: geometry_header.o global.o: material_header.o @@ -78,8 +83,9 @@ global.o: source_header.o global.o: tally_header.o global.o: timing.o +initialize.o: ace.o +initialize.o: bank_header.o initialize.o: constants.o -initialize.o: cross_section.o initialize.o: datatypes.o initialize.o: datatypes_header.o initialize.o: energy_grid.o @@ -88,8 +94,6 @@ initialize.o: geometry.o initialize.o: geometry_header.o initialize.o: global.o initialize.o: input_xml.o -initialize.o: logging.o -initialize.o: mpi_routines.o initialize.o: output.o initialize.o: random_lcg.o initialize.o: source.o @@ -113,6 +117,14 @@ input_xml.o: xml-fortran/templates/materials_t.o input_xml.o: xml-fortran/templates/settings_t.o input_xml.o: xml-fortran/templates/tallies_t.o +intercycle.o: error.o +intercycle.o: global.o +intercycle.o: output.o +intercycle.o: particle_header.o +intercycle.o: random_lcg.o +intercycle.o: tally_header.o +intercycle.o: timing.o + interpolation.o: constants.o interpolation.o: endf_header.o interpolation.o: error.o @@ -127,7 +139,7 @@ main.o: cmfd_execute.o main.o: constants.o main.o: global.o main.o: initialize.o -main.o: mpi_routines.o +main.o: intercycle.o main.o: output.o main.o: particle_header.o main.o: physics.o @@ -140,15 +152,6 @@ main.o: timing.o mesh.o: mesh_header.o -mpi_routines.o: constants.o -mpi_routines.o: error.o -mpi_routines.o: global.o -mpi_routines.o: output.o -mpi_routines.o: particle_header.o -mpi_routines.o: random_lcg.o -mpi_routines.o: tally_header.o -mpi_routines.o: timing.o - output.o: constants.o output.o: datatypes.o output.o: endf.o @@ -162,8 +165,9 @@ output.o: tally_header.o particle_header.o: constants.o particle_header.o: geometry_header.o +physics.o: ace_header.o physics.o: constants.o -physics.o: cross_section_header.o +physics.o: cross_section.o physics.o: endf.o physics.o: error.o physics.o: fission.o @@ -189,9 +193,9 @@ search.o: constants.o search.o: error.o search.o: global.o +source.o: ace_header.o source.o: bank_header.o source.o: constants.o -source.o: cross_section_header.o source.o: error.o source.o: global.o source.o: output.o @@ -207,9 +211,9 @@ string.o: global.o tally.o: constants.o tally.o: error.o tally.o: global.o +tally.o: intercycle.o tally.o: mesh.o tally.o: mesh_header.o -tally.o: mpi_routines.o tally.o: output.o tally.o: search.o tally.o: string.o diff --git a/src/Makefile b/src/Makefile index 7de59c267d..d1a28a27a1 100644 --- a/src/Makefile +++ b/src/Makefile @@ -118,7 +118,7 @@ endif ifeq ($(COMPILER),ibm) F90 = xlf2003 - F90FLAGS := -WF,-DNO_F2008 + F90FLAGS := -WF,-DNO_F2008 -O2 # Debugging ifeq ($(DEBUG),yes) @@ -166,11 +166,19 @@ ifeq ($(USE_MPI),yes) endif #=============================================================================== -# Special options for ORNL Jaguar supercomputer +# Options for IBM Blue Gene/P ANL supercomputer #=============================================================================== -HOSTNAME = $(shell hostname) -ifneq (,$(findstring jaguar,$(HOSTNAME))) +ifeq ($(MACHINE),bluegene) + F90 = /bgsys/drivers/ppcfloor/comm/xl/bin/mpixlf2003 + F90FLAGS = -WF,-DNO_F2008,-DMPI -O3 +endif + +#=============================================================================== +# Options for Cray XK6 ORNL Jaguar supercomputer +#=============================================================================== + +ifeq ($(MACHINE),crayxk6) F90 = ftn F90FLAGS += -DMPI endif @@ -181,7 +189,7 @@ endif all: xml-fortran $(program) xml-fortran: - cd xml-fortran; make F90=$(F90) F90FLAGS="$(F90FLAGS)" + cd xml-fortran; make MACHINE=$(MACHINE) F90=$(F90) F90FLAGS="$(F90FLAGS)" cd xml-fortran/templates; make F90=$(F90) F90FLAGS="$(F90FLAGS)" $(program): $(objects) $(F90) $(objects) $(templates) $(xml_fort) -o $@ $(LDFLAGS) diff --git a/src/OBJECTS b/src/OBJECTS index 543556e8da..cda813fc9b 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -1,9 +1,10 @@ objects = \ +ace.o \ +ace_header.o \ bank_header.o \ cmfd_execute.o \ cmfd_header.o \ cmfd_utils.o \ -cross_section_header.o \ cross_section.o \ datatypes.o \ datatypes_header.o \ @@ -12,20 +13,18 @@ endf.o \ endf_header.o \ energy_grid.o \ error.o \ -fileio.o \ fission.o \ geometry.o \ geometry_header.o \ global.o \ initialize.o \ +intercycle.o \ interpolation.o \ input_xml.o \ -logging.o \ main.o \ material_header.o \ mesh_header.o \ mesh.o \ -mpi_routines.o \ output.o \ particle_header.o \ physics.o \ diff --git a/src/ace.F90 b/src/ace.F90 new file mode 100644 index 0000000000..32894116c3 --- /dev/null +++ b/src/ace.F90 @@ -0,0 +1,1314 @@ +module ace + + use ace_header, only: Nuclide, Reaction, SAB_Table, XsListing, & + DistEnergy + use constants + use datatypes, only: dict_create, dict_add_key, dict_get_key, & + dict_has_key, dict_delete, dict_keys + use datatypes_header, only: DictionaryCI, ListKeyValueCI + use endf, only: reaction_name + use error, only: fatal_error + use fission, only: nu_total + use global + use material_header, only: Material + use output, only: write_message, print_nuclide, header + use string, only: split_string, str_to_int, str_to_real, & + lower_case, to_str + + implicit none + + integer :: NXS(16) ! Descriptors for ACE XSS tables + integer :: JXS(32) ! Pointers into ACE XSS tables + real(8), allocatable :: XSS(:) ! Cross section data + integer :: XSS_index ! current index in XSS data + + type(DictionaryCI), pointer :: already_read => null() + + private :: NXS + private :: JXS + private :: XSS + +contains + +!=============================================================================== +! READ_XS reads all the cross sections for the problem and stores them in +! nuclides and sab_tables arrays +!=============================================================================== + + subroutine read_xs() + + integer :: i ! index in materials array + integer :: j ! index over nuclides in material + integer :: index_list ! index in xs_listings array + integer :: index_nuclides ! index in nuclides + integer :: index_sab ! index in sab_tables + character(12) :: name ! name of isotope, e.g. 92235.03c + character(12) :: alias ! alias of nuclide, e.g. U-235.03c + type(Material), pointer :: mat => null() + type(Nuclide), pointer :: nuc => null() + type(SAB_Table), pointer :: sab => null() + + ! ========================================================================== + ! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES + + ! First we need to determine how many continuous-energy tables and how many + ! S(a,b) thermal scattering tables there are -- this loop doesn't actually + ! read the data, it simply counts the number of nuclides and S(a,b) tables + index_nuclides = 0 + index_sab = 0 + do i = 1, n_materials + ! Get pointer to material + mat => materials(i) + + ! First go through all the nuclide and check if they exist on other + ! materials -- if not, then increment the count for the number of + ! nuclides and add to dictionary + do j = 1, mat % n_nuclides + name = mat % names(j) + + ! First check that this nuclide is listed in the cross_sections.xml + ! file + if (.not. dict_has_key(xs_listing_dict, name)) then + message = "Could not find nuclide " // trim(name) // & + " in cross_sections.xml file!" + call fatal_error() + end if + + ! Find index in xs_listing and set the name and alias according to the + ! listing + index_list = dict_get_key(xs_listing_dict, name) + name = xs_listings(index_list) % name + alias = xs_listings(index_list) % alias + + ! If this nuclide hasn't been encountered yet, we need to add its name + ! and alias to the nuclide_dict + if (.not. dict_has_key(nuclide_dict, name)) then + index_nuclides = index_nuclides + 1 + mat % nuclide(j) = index_nuclides + + call dict_add_key(nuclide_dict, name, index_nuclides) + call dict_add_key(nuclide_dict, alias, index_nuclides) + else + mat % nuclide(j) = dict_get_key(nuclide_dict, name) + end if + end do + + ! Check if S(a,b) exists on other materials + if (mat % has_sab_table) then + name = mat % sab_name + + ! First check that this nuclide is listed in the cross_sections.xml + ! file + if (.not. dict_has_key(xs_listing_dict, name)) then + message = "Could not find S(a,b) table " // trim(name) // & + " in cross_sections.xml file!" + call fatal_error() + end if + + ! Find index in xs_listing and set the name and alias according to the + ! listing + index_list = dict_get_key(xs_listing_dict, name) + name = xs_listings(index_list) % name + alias = xs_listings(index_list) % alias + + ! If this S(a,b) table hasn't been encountered yet, we need to add its + ! name and alias to the sab_dict + if (.not. dict_has_key(sab_dict, name)) then + index_sab = index_sab + 1 + mat % sab_table = index_sab + + call dict_add_key(sab_dict, name, index_sab) + call dict_add_key(sab_dict, alias, index_sab) + else + mat % sab_table = dict_get_key(sab_dict, name) + end if + end if + end do + + n_nuclides_total = index_nuclides + n_sab_tables = index_sab + + ! allocate arrays for ACE table storage and cross section cache + allocate(nuclides(n_nuclides_total)) + allocate(sab_tables(n_sab_tables)) + allocate(micro_xs(n_nuclides_total)) + + ! ========================================================================== + ! READ ALL ACE CROSS SECTION TABLES + + ! display header in summary.out + if (master) call header("CROSS SECTION TABLES", unit=UNIT_SUMMARY) + + call dict_create(already_read) + + ! Loop over all files + do i = 1, n_materials + mat => materials(i) + + do j = 1, mat % n_nuclides + name = mat % names(j) + + if (.not. dict_has_key(already_read, name)) then + index_list = dict_get_key(xs_listing_dict, name) + index_nuclides = dict_get_key(nuclide_dict, name) + name = xs_listings(index_list) % name + alias = xs_listings(index_list) % alias + + call read_ace_table(index_nuclides, index_list) + + ! Print out information on table to summary.out file + nuc => nuclides(index_nuclides) + if (master) call print_nuclide(nuc, unit=UNIT_SUMMARY) + + call dict_add_key(already_read, name, 0) + call dict_add_key(already_read, alias, 0) + end if + end do + + if (mat % has_sab_table) then + name = mat % sab_name + + if (.not. dict_has_key(already_read, name)) then + index_list = dict_get_key(xs_listing_dict, name) + index_sab = dict_get_key(sab_dict, name) + + call read_ace_table(index_sab, index_list) + + call dict_add_key(already_read, name, 0) + end if + end if + end do + + call dict_delete(already_read) + + ! ========================================================================== + ! ASSIGN S(A,B) TABLES TO SPECIFIC NUCLIDES WITHIN MATERIALS + + do i = 1, n_materials + mat => materials(i) + + if (mat % has_sab_table) then + ! In order to know which nuclide the S(a,b) table applies to, we need + ! to search through the list of nuclides for one which has a matching + ! zaid + sab => sab_tables(mat % sab_table) + + do j = 1, mat % n_nuclides + nuc => nuclides(mat % nuclide(j)) + if (nuc % zaid == sab % zaid) then + mat % sab_nuclide = j + end if + end do + + ! Check to make sure S(a,b) table matched a nuclide + if (mat % sab_nuclide == 0) then + message = "S(a,b) table " // trim(mat % sab_name) // " did not match & + &any nuclide on material " // trim(to_str(mat % id)) + call fatal_error() + end if + end if + end do + + end subroutine read_xs + +!=============================================================================== +! READ_ACE_BINARY reads a single cross section table in binary format. This +! routine reads the header data for each table and then calls appropriate +! subroutines to parse the actual data. +!=============================================================================== + + subroutine read_ace_table(index_table, index_list) + + integer, intent(in) :: index_table ! index in nuclides/sab_tables + integer, intent(in) :: index_list ! index in xs_listings + + integer :: i ! loop index for XSS records + integer :: j, j1, j2 ! indices in XSS + integer :: record_length ! Fortran record length + integer :: location ! location of ACE table + integer :: entries ! number of entries on each record + integer :: length ! length of ACE table + integer :: in = 7 ! file unit + integer :: zaids(16) ! list of ZAIDs (only used for S(a,b)) + integer :: filetype ! filetype (ASCII or BINARY) + real(8) :: kT ! temperature of table + real(8) :: awrs(16) ! list of atomic weight ratios (not used) + real(8) :: awr ! atomic weight ratio for table + logical :: file_exists ! does ACE library exist? + character(7) :: readable ! is ACE library readable? + character(10) :: name ! name of ACE table + character(10) :: date_ ! date ACE library was processed + character(10) :: mat ! material identifier + character(70) :: comment ! comment for ACE table + character(MAX_FILE_LEN) :: filename ! path to ACE cross section library + type(Nuclide), pointer :: nuc => null() + type(SAB_Table), pointer :: sab => null() + type(XsListing), pointer :: listing => null() + + ! determine path, record length, and location of table + listing => xs_listings(index_list) + filename = listing % path + record_length = listing % recl + location = listing % location + entries = listing % entries + filetype = listing % filetype + + ! Check if ACE library exists and is readable + inquire(FILE=filename, EXIST=file_exists, READ=readable) + if (.not. file_exists) then + message = "ACE library '" // trim(filename) // "' does not exist!" + call fatal_error() + elseif (readable(1:3) == 'NO') then + message = "ACE library '" // trim(filename) // "' is not readable! & + &Change file permissions with chmod command." + call fatal_error() + end if + + ! display message + message = "Loading ACE cross section table: " // listing % name + call write_message(6) + + if (filetype == ASCII) then + ! ======================================================================= + ! READ ACE TABLE IN ASCII FORMAT + + ! Find location of table + open(UNIT=in, FILE=filename, STATUS='old', ACTION='read') + rewind(UNIT=in) + do i = 1, location - 1 + read(UNIT=in, FMT=*) + end do + + ! Read first line of header + read(UNIT=in, FMT='(A10,2E12.0,1X,A10)') name, awr, kT, date_ + + ! Read more header and NXS and JXS + 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) + + ! determine table length + length = NXS(1) + allocate(XSS(length)) + + ! Read XSS array + read(UNIT=in, FMT='(4E20.0)') XSS + + elseif (filetype == BINARY) then + ! ======================================================================= + ! READ ACE TABLE IN BINARY FORMAT + + open(UNIT=in, FILE=filename, STATUS='old', ACTION='read', & + ACCESS='direct', RECL=record_length) + + ! Read all header information + read(UNIT=in, REC=location) name, awr, kT, date_, & + comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS + + ! determine table length + length = NXS(1) + allocate(XSS(length)) + + ! Read remaining records with XSS + do i = 1, (length + entries - 1)/entries + j1 = 1 + (i-1)*entries + j2 = min(length, j1 + entries - 1) + read(UNIT=IN, REC=location + i) (XSS(j), j=j1,j2) + end do + end if + + ! ========================================================================== + ! PARSE DATA BASED ON NXS, JXS, AND XSS ARRAYS + + select case(listing % type) + case (ACE_NEUTRON) + nuc => nuclides(index_table) + nuc % name = name + nuc % awr = awr + nuc % kT = kT + nuc % zaid = NXS(2) + + ! read all blocks + call read_esz(nuc) + call read_nu_data(nuc) + call read_reactions(nuc) + call read_angular_dist(nuc) + call read_energy_dist(nuc) + call read_unr_res(nuc) + + ! for fissionable nuclides, precalculate microscopic nu-fission cross + ! sections so that we don't need to call the nu_total function during + ! cross section lookups + + if (nuc % fissionable) call generate_nu_fission(nuc) + + case (ACE_THERMAL) + sab => sab_tables(index_table) + sab % name = name + sab % awr = awr + sab % kT = kT + sab % zaid = zaids(1) + + call read_thermal_data(sab) + end select + + deallocate(XSS) + if(associated(nuc)) nullify(nuc) + if(associated(sab)) nullify(sab) + + close(UNIT=in) + + end subroutine read_ace_table + +!=============================================================================== +! READ_ESZ - reads through the ESZ block. This block contains the energy grid, +! total xs, absorption xs, elastic scattering xs, and heating numbers. +!=============================================================================== + + subroutine read_esz(nuc) + + type(Nuclide), pointer :: nuc + + integer :: NE ! number of energy points for total and elastic cross sections + + ! determine number of energy points + NE = NXS(3) + nuc % n_grid = NE + + ! allocate storage for energy grid and cross section arrays + allocate(nuc % energy(NE)) + allocate(nuc % total(NE)) + allocate(nuc % elastic(NE)) + allocate(nuc % fission(NE)) + allocate(nuc % nu_fission(NE)) + allocate(nuc % absorption(NE)) + allocate(nuc % heating(NE)) + + ! initialize cross sections + nuc % total = ZERO + nuc % elastic = ZERO + nuc % fission = ZERO + nuc % nu_fission = ZERO + nuc % absorption = ZERO + nuc % heating = ZERO + + ! Read data from XSS -- only the energy grid, elastic scattering and heating + ! cross section values are actually read from here. The total and absorption + ! cross sections are reconstructed from the partial reaction data. + + XSS_index = 1 + nuc % energy = get_real(NE) + + ! Skip total and absorption + XSS_index = XSS_index + 2*NE + + ! Continue reading elastic scattering and heating + nuc % elastic = get_real(NE) + nuc % heating = get_real(NE) + + end subroutine read_esz + +!=============================================================================== +! 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 read_nu_data(nuc) + + type(Nuclide), pointer :: nuc + + integer :: i ! loop index + integer :: JXS2 ! location for fission nu data + integer :: JXS24 ! location for delayed neutron data + 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 :: NPCR ! number of delayed neutron precursor groups + integer :: LED ! location of energy distribution locators + integer :: LDIS ! location of all energy distributions + integer :: LOCC ! location of energy distributions for given MT + integer :: lc ! locator + integer :: length ! length of data to allocate + type(DistEnergy), pointer :: edist => null() + + JXS2 = JXS(2) + JXS24 = JXS(24) + + if (JXS2 == 0) then + ! ======================================================================= + ! NO PROMPT/TOTAL NU DATA + nuc % nu_t_type = NU_NONE + nuc % nu_p_type = NU_NONE + + elseif (XSS(JXS2) > 0) then + ! ======================================================================= + ! PROMPT OR TOTAL NU DATA + KNU = JXS2 + LNU = int(XSS(KNU)) + if (LNU == 1) then + ! Polynomial data + nuc % nu_t_type = NU_POLYNOMIAL + nuc % 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 + nuc % nu_t_type = NU_TABULAR + nuc % 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(nuc % 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 + nuc % nu_t_data = get_real(length) + + elseif (XSS(JXS2) < 0) then + ! ======================================================================= + ! PROMPT AND TOTAL NU DATA -- read prompt data first + KNU = JXS2 + 1 + LNU = int(XSS(KNU)) + if (LNU == 1) then + ! Polynomial data + nuc % nu_p_type = NU_POLYNOMIAL + + ! allocate determine how many coefficients for polynomial + NC = int(XSS(KNU+1)) + length = NC + 1 + elseif (LNU == 2) then + ! Tabular data + nuc % nu_p_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(nuc % nu_p_data(length)) + + ! read data + XSS_index = KNU + 1 + nuc % nu_p_data = get_real(length) + + ! Now read total nu data + KNU = JXS2 + int(abs(XSS(JXS2))) + 1 + LNU = int(XSS(KNU)) + if (LNU == 1) then + ! Polynomial data + nuc % 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 + nuc % 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(nuc % nu_t_data(length)) + + ! read data + XSS_index = KNU + 1 + nuc % nu_t_data = get_real(length) + end if + + if (JXS24 > 0) then + ! ======================================================================= + ! DELAYED NU DATA + + nuc % 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(nuc % nu_d_data(length)) + + ! read delayed nu data + XSS_index = KNU + 1 + nuc % nu_d_data = get_real(length) + + ! ======================================================================= + ! DELAYED NEUTRON ENERGY DISTRIBUTION + + ! Allocate space for secondary energy distribution + NPCR = NXS(8) + nuc % n_precursor = NPCR + allocate(nuc % nu_d_edist(NPCR)) + + LED = JXS(26) + LDIS = JXS(27) + + ! Loop over all delayed neutron precursor groups + do i = 1, NPCR + ! find location of energy distribution data + LOCC = int(XSS(LED + i - 1)) + + ! read energy distribution data + edist => nuc % nu_d_edist(i) + call get_energy_dist(edist, LOCC, .true.) + end do + + ! ======================================================================= + ! DELAYED NEUTRON PRECUSOR YIELDS AND CONSTANTS + + ! determine length of all precursor constants/yields/interp data + length = 0 + lc = JXS(25) + do i = 1, NPCR + NR = int(XSS(lc + length + 1)) + NE = int(XSS(lc + length + 2 + 2*NR)) + length = length + 3 + 2*NR + 2*NE + end do + + ! allocate space for precusor data + allocate(nuc % nu_d_precursor_data(length)) + + ! read delayed neutron precursor data + XSS_index = lc + nuc % nu_d_precursor_data = get_real(length) + + else + nuc % 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(nuc) + + type(Nuclide), pointer :: nuc + + integer :: i ! loop indices + integer :: i_fission ! index in nuc % index_fission + integer :: LMT ! index of MT list in XSS + integer :: NMT ! Number of reactions + integer :: JXS4 ! index of Q values in XSS + integer :: JXS5 ! index of neutron multiplicities in XSS + integer :: JXS7 ! index of reactions cross-sections in XSS + integer :: LXS ! location of cross-section locators + integer :: LOCA ! location of cross-section for given MT + integer :: IE ! reaction's starting index on energy grid + integer :: NE ! number of energies for reaction + type(Reaction), pointer :: rxn => null() + + LMT = JXS(3) + JXS4 = JXS(4) + JXS5 = JXS(5) + LXS = JXS(6) + JXS7 = JXS(7) + NMT = NXS(4) + + ! allocate array of reactions. Add one since we need to include an elastic + ! scattering channel + nuc % n_reaction = NMT + 1 + allocate(nuc % reactions(NMT+1)) + + ! Store elastic scattering cross-section on reaction one + rxn => nuc % reactions(1) + rxn % MT = 2 + rxn % Q_value = ZERO + rxn % TY = 1 + rxn % IE = 1 + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. + allocate(rxn % sigma(nuc % n_grid)) + rxn % sigma = nuc % elastic + + ! Add contribution of elastic scattering to total cross section + nuc % total = nuc % total + nuc % elastic + + ! By default, set nuclide to not fissionable and then change if fission + ! reactions are encountered + nuc % fissionable = .false. + nuc % has_partial_fission = .false. + nuc % n_fission = 0 + i_fission = 0 + + do i = 1, NMT + rxn => nuc % reactions(i+1) + + ! set defaults + rxn % has_angle_dist = .false. + rxn % has_energy_dist = .false. + + ! read MT number, Q-value, and neutrons produced + rxn % MT = int(XSS(LMT + i - 1)) + rxn % Q_value = XSS(JXS4 + i - 1) + rxn % TY = int(XSS(JXS5 + i - 1)) + + ! read starting energy index + LOCA = int(XSS(LXS + i - 1)) + IE = int(XSS(JXS7 + LOCA - 1)) + rxn % IE = IE + + ! read number of energies cross section values + NE = int(XSS(JXS7 + LOCA)) + allocate(rxn % sigma(NE)) + XSS_index = JXS7 + LOCA + 1 + rxn % sigma = get_real(NE) + + ! Skip redundant reactions -- this includes total inelastic level + ! scattering, gas production cross sections (MT=200+), and (n,p), (n,d), + ! etc. reactions leaving the nucleus in an excited state + if (rxn % MT == N_LEVEL .or. rxn % MT > N_DA) cycle + + ! Add contribution to total cross section + nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma + + ! Add contribution to absorption cross section + if (rxn % MT >= N_GAMMA .and. rxn % MT <= N_DA) then + nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma + end if + + ! Information about fission reactions + if (rxn % MT == N_FISSION) then + allocate(nuc % index_fission(1)) + elseif (rxn % MT == N_F) then + allocate(nuc % index_fission(PARTIAL_FISSION_MAX)) + nuc % has_partial_fission = .true. + end if + + ! Add contribution to fission cross section + if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & + .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) then + nuc % fissionable = .true. + nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma + + ! Also need to add fission cross sections to absorption + nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma + + ! Keep track of this reaction for easy searching later + i_fission = i_fission + 1 + nuc % index_fission(i_fission) = i + 1 + nuc % n_fission = nuc % n_fission + 1 + end if + end do + + end subroutine read_reactions + +!=============================================================================== +! READ_ANGULAR_DIST parses the angular distribution for each reaction with +! secondary neutrons +!=============================================================================== + + subroutine read_angular_dist(nuc) + + type(Nuclide), pointer :: nuc + + integer :: JXS8 ! location of angular distribution locators + integer :: JXS9 ! location of angular distributions + integer :: LOCB ! location of angular distribution for given MT + integer :: NE ! number of incoming energies + integer :: NP ! number of points for cosine distribution + integer :: LC ! locator + integer :: i ! index in reactions array + integer :: j ! index over incoming energies + integer :: length ! length of data array to allocate + type(Reaction), pointer :: rxn => null() + + JXS8 = JXS(8) + JXS9 = JXS(9) + + ! loop over all reactions with secondary neutrons -- NXS(5) does not include + ! elastic scattering + do i = 1, NXS(5) + 1 + rxn => nuc%reactions(i) + + ! find location of angular distribution + LOCB = int(XSS(JXS8 + i - 1)) + if (LOCB == -1) then + ! Angular distribution data are specified through LAWi = 44 in the DLW + ! block + cycle + elseif (LOCB == 0) then + ! No angular distribution data are given for this reaction, isotropic + ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0) + cycle + end if + rxn % has_angle_dist = .true. + + ! allocate space for incoming energies and locations + NE = int(XSS(JXS9 + LOCB - 1)) + 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 nucs + XSS_index = JXS9 + LOCB + rxn % adist % energy = get_real(NE) + rxn % adist % location = get_int(NE) + + ! determine dize of data block + length = 0 + do j = 1, NE + LC = rxn % adist % location(j) + if (LC == 0) then + ! isotropic + rxn % adist % type(j) = ANGLE_ISOTROPIC + elseif (LC > 0) then + ! 32 equiprobable bins + rxn % adist % type(j) = ANGLE_32_EQUI + length = length + 33 + elseif (LC < 0) then + ! tabular distribution + rxn % adist % type(j) = ANGLE_TABULAR + NP = int(XSS(JXS9 + abs(LC))) + length = length + 2 + 3*NP + end if + end do + + ! allocate angular distribution data and read + 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) + XSS_index = JXS9 + abs(LC) - 1 + 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 + + end do + + end subroutine read_angular_dist + +!=============================================================================== +! READ_ENERGY_DIST parses the secondary energy distribution for each reaction +! with seconary neutrons (except elastic scattering) +!=============================================================================== + + subroutine read_energy_dist(nuc) + + type(Nuclide), pointer :: nuc + + integer :: LED ! location of energy distribution locators + integer :: LOCC ! location of energy distributions for given MT + integer :: i ! loop index + type(Reaction), pointer :: rxn => null() + + LED = JXS(10) + + ! Loop over all reactions + do i = 1, NXS(5) + rxn => nuc % reactions(i+1) ! skip over elastic scattering + rxn % has_energy_dist = .true. + + ! find location of energy distribution data + LOCC = int(XSS(LED + i - 1)) + + ! allocate energy distribution + allocate(rxn % edist) + + ! read data for energy distribution + call get_energy_dist(rxn % edist, LOCC) + end do + + end subroutine read_energy_dist + +!=============================================================================== +! GET_ENERGY_DIST reads in data for a single law for an energy distribution and +! calls itself recursively if there are multiple energy distributions for a +! single reaction +!=============================================================================== + + recursive subroutine get_energy_dist(edist, loc_law, delayed_n) + + type(DistEnergy), pointer :: edist ! energy distribution + integer, intent(in) :: loc_law ! locator for data + logical, optional :: delayed_n ! is this for delayed neutrons? + + integer :: LDIS ! location of all energy distributions + integer :: LNW ! location of next energy distribution if multiple + 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 + integer :: lc ! locator + integer :: length ! length of data to allocate + integer :: length_interp_data ! length of interpolation data + + ! determine location of energy distribution + if (present(delayed_n)) then + LDIS = JXS(27) + else + LDIS = JXS(11) + end if + + ! locator for next law and information on this law + LNW = int(XSS(LDIS + loc_law - 1)) + LAW = int(XSS(LDIS + loc_law)) + IDAT = int(XSS(LDIS + loc_law + 1)) + NR = int(XSS(LDIS + loc_law + 2)) + edist % law = LAW + edist % p_valid % n_regions = NR + + ! allocate space for ENDF interpolation parameters + if (NR > 0) then + allocate(edist % p_valid % nbt(NR)) + allocate(edist % p_valid % int(NR)) + end if + + ! read ENDF interpolation parameters + XSS_index = LDIS + loc_law + 3 + if (NR > 0) then + edist % p_valid % nbt = int(get_real(NR)) + edist % p_valid % int = int(get_real(NR)) + end if + + ! allocate space for law validity data + NE = int(XSS(LDIS + loc_law + 3 + 2*NR)) + edist % p_valid % n_pairs = NE + allocate(edist % p_valid % x(NE)) + allocate(edist % p_valid % y(NE)) + + length_interp_data = 5 + 2*(NR + NE) + + ! read law validity data + XSS_index = LDIS + loc_law + 4 + 2*NR + edist % p_valid % x = get_real(NE) + edist % p_valid % y = get_real(NE) + + ! Set index to beginning of IDAT array + lc = LDIS + IDAT - 2 + + ! determine length of energy distribution + length = length_energy_dist(lc, LAW, loc_law, length_interp_data) + + ! allocate secondary energy distribution array + allocate(edist % data(length)) + + ! read secondary energy distribution + XSS_index = lc + 1 + edist % data = get_real(length) + + ! read next energy distribution if present + if (LNW > 0) then + allocate(edist % next) + call get_energy_dist(edist % next, LNW) + end if + + end subroutine get_energy_dist + +!=============================================================================== +! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy +! distribution array based on the secondary energy law and location in XSS +!=============================================================================== + + function length_energy_dist(lc, law, LOCC, lid) result(length) + + integer, intent(in) :: lc ! location in XSS array + integer, intent(in) :: law ! energy distribution law + integer, intent(in) :: LOCC ! location of energy distribution + integer, intent(in) :: lid ! length of interpolation data + integer :: length ! length of energy distribution (LDAT) + + integer :: i,j,k + integer :: NR ! number of interpolation regions + integer :: NE ! number of incoming energies + integer :: NP ! number of points in outgoing energy distribution + integer :: NMU ! number of points in outgoing cosine distribution + integer :: NRa ! number of interpolation regions for Watt 'a' + integer :: NEa ! number of energies for Watt 'a' + integer :: NRb ! number of interpolation regions for Watt 'b' + integer :: NEb ! number of energies for Watt 'b' + + ! initialize length + length = 0 + + select case (law) + case (1) + ! Tabular equiprobable energy bins + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + NP = int(XSS(lc + 3 + 2*NR + NE)) + length = 3 + 2*NR + NE + 3*NP*NE + + case (2) + ! Discrete photon energy + length = 2 + + case (3) + ! Level scattering + length = 2 + + case (4) + ! Continuous tabular distribution + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + length = length + 2 + 2*NR + 2*NE + do i = 1,NE + ! determine length + NP = int(XSS(lc + length + 2)) + length = length + 2 + 3*NP + + ! adjust location for this block + j = lc + 2 + 2*NR + NE + i + XSS(j) = XSS(j) - LOCC - lid + end do + + case (5) + ! General evaporation spectrum + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + NP = int(XSS(lc + 3 + 2*NR + 2*NE)) + length = 3 + 2*NR + 2*NE + NP + + case (7) + ! Maxwell fission spectrum + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + length = 3 + 2*NR + 2*NE + + case (9) + ! Evaporation spectrum + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + length = 3 + 2*NR + 2*NE + + case (11) + ! Watt spectrum + NRa = int(XSS(lc + 1)) + NEa = int(XSS(lc + 2 + 2*NRa)) + NRb = int(XSS(lc + 3 + 2*(NRa+NEa))) + NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb))) + length = 5 + 2*(NRa + NEa + NRb + NEb) + + case (44) + ! Kalbach-Mann correlated scattering + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + length = length + 2 + 2*NR + 2*NE + do i = 1,NE + NP = int(XSS(lc + length + 2)) + length = length + 2 + 5*NP + + ! adjust location for this block + j = lc + 2 + 2*NR + NE + i + XSS(j) = XSS(j) - LOCC - lid + end do + + case (61) + ! Correlated energy and angle distribution + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + length = length + 2 + 2*NR + 2*NE + do i = 1,NE + ! outgoing energy distribution + NP = int(XSS(lc + length + 2)) + + ! adjust locators for angular distribution + do j = 1, NP + k = lc + length + 2 + 3*NP + j + if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid + end do + + length = length + 2 + 4*NP + do j = 1, NP + ! outgoing angle distribution -- NMU here is actually + ! referred to as NP in the MCNP documentation + NMU = int(XSS(lc + length + 2)) + length = length + 2 + 3*NMU + end do + + ! adjust locators for energy distribution + j = lc + 2 + 2*NR + NE + i + XSS(j) = XSS(j) - LOCC - lid + end do + + case (66) + ! N-body phase space distribution + length = 2 + + case (67) + ! Laboratory energy-angle law + NR = int(XSS(lc + 1)) + NE = int(XSS(lc + 2 + 2*NR)) + NMU = int(XSS(lc + 4 + 2*NR + 2*NE)) + length = 4 + 2*(NR + NE + NMU) + + end select + + end function length_energy_dist + +!=============================================================================== +! READ_UNR_RES reads in unresolved resonance probability tables if present. +!=============================================================================== + + subroutine read_unr_res(nuc) + + type(Nuclide), pointer :: nuc + + integer :: JXS23 ! location of URR data + integer :: lc ! locator + integer :: N ! # of incident energies + integer :: M ! # of probabilities + integer :: i ! index over incoming energies + integer :: j ! index over values + integer :: k ! index over probabilities + + ! determine locator for URR data + JXS23 = JXS(23) + + ! check if URR data is present + if (JXS23 /= 0) then + nuc % urr_present = .true. + allocate(nuc % urr_data) + lc = JXS23 + else + nuc % urr_present = .false. + return + end if + + ! read parameters + nuc % urr_data % n_energy = int(XSS(lc)) + nuc % urr_data % n_prob = int(XSS(lc + 1)) + nuc % urr_data % interp = int(XSS(lc + 2)) + nuc % urr_data % inelastic_flag = int(XSS(lc + 3)) + nuc % urr_data % absorption_flag = int(XSS(lc + 4)) + if (int(XSS(lc + 5)) == 0) then + nuc % urr_data % multiply_smooth = .false. + else + nuc % urr_data % multiply_smooth = .true. + end if + + ! if the inelastic competition flag indicates that the inelastic cross + ! section should be determined from a normal reaction cross section, we need + ! to set up a pointer to that reaction + nuc % urr_inelastic = NONE + if (nuc % urr_data % inelastic_flag > 0) then + do i = 1, nuc % n_reaction + if (nuc % reactions(i) % MT == nuc % urr_data % inelastic_flag) then + nuc % urr_inelastic = i + end if + end do + + ! Abort if no corresponding inelastic reaction was found + if (nuc % urr_inelastic == NONE) then + message = "Could not find inelastic reaction specified on " & + // "unresolved resonance probability table." + call fatal_error() + end if + end if + + ! allocate incident energies and probability tables + N = nuc % urr_data % n_energy + M = nuc % urr_data % n_prob + allocate(nuc % urr_data % energy(N)) + allocate(nuc % urr_data % prob(N,6,M)) + + ! read incident energies + XSS_index = lc + 6 + nuc % urr_data % energy = get_real(N) + + ! read probability tables + do i = 1, N + do j = 1, 6 + do k = 1, M + nuc % urr_data % prob(i,j,k) = XSS(XSS_index) + XSS_index = XSS_index + 1 + end do + end do + end do + + end subroutine read_unr_res + +!=============================================================================== +! GENERATE_NU_FISSION precalculates the microscopic nu-fission cross section for +! a given nuclide. This is done so that the nu_total function does not need to +! be called during cross section lookups. +!=============================================================================== + + subroutine generate_nu_fission(nuc) + + type(Nuclide), pointer :: nuc + + integer :: i ! index on nuclide energy grid + real(8) :: E ! energy + real(8) :: nu ! # of neutrons per fission + + do i = 1, nuc % n_grid + ! determine energy + E = nuc % energy(i) + + ! determine total nu at given energy + nu = nu_total(nuc, E) + + ! determine nu-fission microscopic cross section + nuc % nu_fission(i) = nu * nuc % fission(i) + end do + + end subroutine generate_nu_fission + +!=============================================================================== +! READ_THERMAL_DATA reads elastic and inelastic cross sections and corresponding +! secondary energy/angle distributions derived from experimental S(a,b) +! data. Namely, in MCNP language, this routine reads the ITIE, ITCE, ITXE, and +! ITCA blocks. +!=============================================================================== + + subroutine read_thermal_data(table) + + type(SAB_Table), pointer :: table + + integer :: i ! index for incoming energies + integer :: j ! index for outgoing energies + integer :: k ! index for outoging angles + integer :: lc ! location in XSS array + integer :: NE_in ! number of incoming energies + integer :: NE_out ! number of outgoing energies + integer :: NMU ! number of outgoing angles + integer :: JXS4 ! location of elastic energy table + + ! read secondary energy mode for inelastic scattering + table % secondary_mode = NXS(7) + + ! read number of inelastic energies and allocate arrays + NE_in = int(XSS(JXS(1))) + table % n_inelastic_e_in = NE_in + allocate(table % inelastic_e_in(NE_in)) + allocate(table % inelastic_sigma(NE_in)) + + ! read inelastic energies and cross-sections + XSS_index = JXS(1) + 1 + table % inelastic_e_in = get_real(NE_in) + table % inelastic_sigma = get_real(NE_in) + + ! set threshold value + table % threshold_inelastic = table % inelastic_e_in(NE_in) + + ! 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)) + + ! 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) + + ! 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 + end do + + ! read number of elastic energies and allocate arrays + JXS4 = JXS(4) + if (JXS4 /= 0) then + NE_in = int(XSS(JXS4)) + table % n_elastic_e_in = NE_in + allocate(table % elastic_e_in(NE_in)) + allocate(table % elastic_P(NE_in)) + + ! read elastic energies and P + XSS_index = JXS4 + 1 + table % elastic_e_in = get_real(NE_in) + table % elastic_P = get_real(NE_in) + + ! set threshold + table % threshold_elastic = table % elastic_e_in(NE_in) + + ! determine whether sigma=P or sigma = P/E + table % elastic_mode = NXS(5) + else + table % threshold_elastic = ZERO + table % n_elastic_e_in = 0 + end if + + ! allocate space for outgoing energy/angle for elastic scattering + NMU = NXS(6) + 1 + table % n_elastic_mu = NMU + if (NMU > 0) then + allocate(table % elastic_mu(NMU, NE_in)) + end if + + ! read equiprobable outgoing cosines for elastic scattering each + ! incoming energy + if (JXS4 /= 0 .and. NMU /= 0) then + lc = JXS(6) - 1 + do i = 1, NE_in + do j = 1, NMU + table % elastic_mu(j,i) = XSS(lc + j) + end do + lc = lc + NMU + end do + end if + + end subroutine read_thermal_data + +!=============================================================================== +! GET_INT returns an array of integers read from the current position in the XSS +! array +!=============================================================================== + + function get_int(n_values) result(array) + + integer, intent(in) :: n_values ! number of values to read + integer :: array(n_values) ! array of values + + array = int(XSS(XSS_index:XSS_index + n_values - 1)) + XSS_index = XSS_index + n_values + + end function get_int + +!=============================================================================== +! GET_REAL returns an array of real(8)s read from the current position in the +! XSS array +!=============================================================================== + + function get_real(n_values) result(array) + + integer, intent(in) :: n_values ! number of values to read + real(8) :: array(n_values) ! array of values + + array = XSS(XSS_index:XSS_index + n_values - 1) + XSS_index = XSS_index + n_values + + end function get_real + +end module ace diff --git a/src/cross_section_header.F90 b/src/ace_header.F90 similarity index 94% rename from src/cross_section_header.F90 rename to src/ace_header.F90 index f2a64e306f..0c14a1ca5f 100644 --- a/src/cross_section_header.F90 +++ b/src/ace_header.F90 @@ -1,4 +1,4 @@ -module cross_section_header +module ace_header use constants, only: MAX_FILE_LEN use endf_header, only: Tab1 @@ -51,11 +51,16 @@ module cross_section_header end type Reaction !=============================================================================== -! URRDATA contains unresolved resonance data. +! URRDATA contains probability tables for the unresolved resonance range. !=============================================================================== type UrrData - integer, allocatable :: params(:) + integer :: n_energy ! # of incident neutron energies + integer :: n_prob ! # of probabilities + integer :: interp ! inteprolation (2=lin-lin, 5=log-log) + integer :: inelastic_flag ! inelastic competition flag + integer :: absorption_flag ! other absorption flag + logical :: multiply_smooth ! multiply by smooth cross section? real(8), allocatable :: energy(:) real(8), allocatable :: prob(:,:,:) end type UrrData @@ -108,6 +113,7 @@ module cross_section_header ! Unresolved resonance data logical :: urr_present + integer :: urr_inelastic type(UrrData), pointer :: urr_data => null() ! Reactions @@ -203,4 +209,4 @@ module cross_section_header real(8) :: nu_fission ! macroscopic production xs end type MaterialMacroXS -end module cross_section_header +end module ace_header diff --git a/src/cmfd_utils.F90 b/src/cmfd_utils.F90 index d397663984..2130784406 100644 --- a/src/cmfd_utils.F90 +++ b/src/cmfd_utils.F90 @@ -395,9 +395,9 @@ contains res = leakage + interactions - scattering - (ONE/keff)*fission ! write output - label = "MESH (" // trim(int_to_str(i)) // ". " // & - & trim(int_to_str(j)) // ", " // trim(int_to_str(k)) // & - & ") GROUP " // trim(int_to_str(g)) + label = "MESH (" // trim(int4_to_str(i)) // ". " // & + & trim(int4_to_str(j)) // ", " // trim(int4_to_str(k)) // & + & ") GROUP " // trim(int4_to_str(g)) write(UNIT=UNIT_CMFD, FMT='(A,T35,A)') label, & & trim(real_to_str(res)) diff --git a/src/constants.F90 b/src/constants.F90 index 5b7013b1f5..cb791c0d78 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -22,7 +22,7 @@ module constants real(8), parameter :: TINY_BIT = 1e-8_8 ! User for precision in geometry - real(8), parameter :: FP_PRECISION = 1e-7_8 + real(8), parameter :: FP_PRECISION = 1e-5_8 ! Maximum number of collisions/crossings integer, parameter :: MAX_EVENTS = 10000 @@ -54,57 +54,64 @@ module constants ! GEOMETRY-RELATED CONSTANTS ! Boundary conditions - integer, parameter :: & - & BC_TRANSMIT = 0, & ! Transmission boundary condition (default) - & BC_VACUUM = 1, & ! Vacuum boundary condition - & BC_REFLECT = 2, & ! Reflecting boundary condition - & BC_PERIODIC = 3 ! Periodic boundary condition + integer, parameter :: & + BC_TRANSMIT = 0, & ! Transmission boundary condition (default) + BC_VACUUM = 1, & ! Vacuum boundary condition + BC_REFLECT = 2, & ! Reflecting boundary condition + BC_PERIODIC = 3 ! Periodic boundary condition ! Logical operators for cell definitions - integer, parameter :: & - & OP_LEFT_PAREN = huge(0), & ! Left parentheses - & OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses - & OP_UNION = huge(0) - 2, & ! Union operator - & OP_DIFFERENCE = huge(0) - 3 ! Difference operator + integer, parameter :: & + OP_LEFT_PAREN = huge(0), & ! Left parentheses + OP_RIGHT_PAREN = huge(0) - 1, & ! Right parentheses + OP_UNION = huge(0) - 2, & ! Union operator + OP_DIFFERENCE = huge(0) - 3 ! Difference operator ! Cell types - integer, parameter :: & - & CELL_NORMAL = 1, & ! Cell with a specified material - & CELL_FILL = 2, & ! Cell filled by a separate universe - & CELL_LATTICE = 3, & ! Cell filled with a lattice - & CELL_VOID = -1 + integer, parameter :: & + CELL_NORMAL = 1, & ! Cell with a specified material + CELL_FILL = 2, & ! Cell filled by a separate universe + CELL_LATTICE = 3, & ! Cell filled with a lattice + CELL_VOID = -1 ! Lattice types - integer, parameter :: & - & LATTICE_RECT = 1, & - & LATTICE_HEX = 2 + integer, parameter :: & + LATTICE_RECT = 1, & + LATTICE_HEX = 2 + + ! Lattice boundary crossings + integer, parameter :: & + LATTICE_LEFT = 1, & + LATTICE_RIGHT = 2, & + LATTICE_BOTTOM = 3, & + LATTICE_TOP = 4 ! 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_PLANE = 4, & ! Arbitrary plane - & SURF_CYL_X = 5, & ! Cylinder along x-axis - & SURF_CYL_Y = 6, & ! Cylinder along y-axis - & SURF_CYL_Z = 7, & ! Cylinder along z-axis - & SURF_SPHERE = 8, & ! Sphere - & SURF_BOX_X = 9, & ! Box extending infinitely in x-direction - & SURF_BOX_Y = 10, & ! Box extending infinitely in y-direction - & SURF_BOX_Z = 11, & ! Box extending infinitely in z-direction - & SURF_BOX = 12, & ! Rectangular prism - & SURF_GQ = 13 ! General quadratic surface + 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_PLANE = 4, & ! Arbitrary plane + SURF_CYL_X = 5, & ! Cylinder along x-axis + SURF_CYL_Y = 6, & ! Cylinder along y-axis + SURF_CYL_Z = 7, & ! Cylinder along z-axis + SURF_SPHERE = 8, & ! Sphere + SURF_BOX_X = 9, & ! Box extending infinitely in x-direction + SURF_BOX_Y = 10, & ! Box extending infinitely in y-direction + SURF_BOX_Z = 11, & ! Box extending infinitely in z-direction + SURF_BOX = 12, & ! Rectangular prism + SURF_GQ = 13 ! General quadratic surface ! Surface senses - integer, parameter :: & - & SENSE_POSITIVE = 1, & - & SENSE_NEGATIVE = -1 + integer, parameter :: & + SENSE_POSITIVE = 1, & + SENSE_NEGATIVE = -1 ! ============================================================================ ! CROSS SECTION RELATED CONSTANTS ! Interpolation flag - integer, parameter :: & + integer, parameter :: & HISTOGRAM = 1, & ! y is constant in x LINEAR_LINEAR = 2, & ! y is linear in x LINEAR_LOG = 3, & ! y is linear in ln(x) @@ -119,19 +126,19 @@ module constants ! Angular distribution type integer, parameter :: & - & ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution - & ANGLE_32_EQUI = 2, & ! 32 equiprobable bins - & ANGLE_TABULAR = 3 ! Tabular angular distribution + ANGLE_ISOTROPIC = 1, & ! Isotropic angular distribution + ANGLE_32_EQUI = 2, & ! 32 equiprobable bins + ANGLE_TABULAR = 3 ! Tabular angular distribution ! 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_EQUAL = 0, & ! Equally-likely outgoing energy bins + SAB_SECONDARY_SKEWED = 1 ! Skewed outgoing energy bins ! Elastic mode for S(a,b) elastic scattering integer, parameter :: & - & SAB_ELASTIC_DISCRETE = 3, & ! Sample from discrete cosines - & SAB_ELASTIC_EXACT = 4 ! Exact treatment for coherent elastic + SAB_ELASTIC_DISCRETE = 3, & ! Sample from discrete cosines + SAB_ELASTIC_EXACT = 4 ! Exact treatment for coherent elastic ! Reaction types integer, parameter :: & @@ -201,6 +208,15 @@ module constants ASCII = 1, & ! ASCII cross section file BINARY = 2 ! Binary cross section file + ! Probability table parameters + integer, parameter :: & + URR_CUM_PROB = 1, & + URR_TOTAL = 2, & + URR_ELASTIC = 3, & + URR_FISSION = 4, & + URR_N_GAMMA = 5, & + URR_HEATING = 6 + ! Maximum number of partial fission reactions integer, parameter :: PARTIAL_FISSION_MAX = 4 @@ -280,10 +296,10 @@ module constants PROB_CRITICALITY = 2 ! Criticality problem ! Unit numbers - integer, parameter :: UNIT_LOG = 11 ! unit # for writing log file - integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file - integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file - integer, parameter :: UNIT_CMFD = 14 ! unit # for writing cmfd file + integer, parameter :: UNIT_SUMMARY = 11 ! unit # for writing summary file + integer, parameter :: UNIT_TALLY = 12 ! unit # for writing tally file + integer, parameter :: UNIT_PLOT = 13 ! unit # for writing plot file + integer, parameter :: UNIT_CMFD = 14 ! unit # for writing cmfd file end module constants diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 86d15ef1d1..367d277cf5 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -1,1284 +1,408 @@ module cross_section + use ace_header, only: Nuclide, SAB_Table, Reaction, UrrData use constants - use cross_section_header, only: Nuclide, Reaction, SAB_Table, XsListing, & - DistEnergy - use datatypes, only: dict_create, dict_add_key, dict_get_key, & - dict_has_key, dict_delete, dict_keys - use datatypes_header, only: DictionaryCI, ListKeyValueCI - use endf, only: reaction_name - use error, only: fatal_error - use fileio, only: read_line, skip_lines - use fission, only: nu_total + use error, only: fatal_error + use fission, only: nu_total use global - use material_header, only: Material - use output, only: write_message - use string, only: split_string, str_to_int, str_to_real, & - lower_case, int_to_str + use material_header, only: Material + use particle_header, only: Particle + use random_lcg, only: prn + use search, only: binary_search implicit none - integer :: NXS(16) ! Descriptors for ACE XSS tables - integer :: JXS(32) ! Pointers into ACE XSS tables - real(8), allocatable :: XSS(:) ! Cross section data - integer :: XSS_index ! current index in XSS data - - type(DictionaryCI), pointer :: already_read => null() - - private :: NXS - private :: JXS - private :: XSS - contains !=============================================================================== -! READ_XS reads all the cross sections for the problem and stores them in -! nuclides and sab_tables arrays +! CALCULATE_XS determines the macroscopic cross sections for the material the +! particle is currently traveling through. !=============================================================================== - subroutine read_xs() + subroutine calculate_xs(p) - integer :: i ! index in materials array - integer :: j ! index over nuclides in material - integer :: index_list ! index in xs_listings array - integer :: index_nuclides ! index in nuclides - integer :: index_sab ! index in sab_tables - character(12) :: name ! name of isotope, e.g. 92235.03c - character(12) :: alias ! alias of nuclide, e.g. U-235.03c - type(Material), pointer :: mat => null() + type(Particle), pointer :: p + + integer :: i ! loop index over nuclides + integer :: index_nuclide ! index into nuclides array + integer :: index_sab ! index into sab_tables array + real(8) :: atom_density ! atom density of a nuclide + real(8) :: sab_threshold ! threshold for S(a,b) table + type(Material), pointer :: mat => null() ! current material + + ! Set all material macroscopic cross sections to zero + material_xs % total = ZERO + material_xs % elastic = ZERO + material_xs % absorption = ZERO + material_xs % fission = ZERO + material_xs % nu_fission = ZERO + + mat => materials(p % material) + + ! Find energy index on unionized grid + call find_energy_index(p) + + ! Check if there's an S(a,b) table for this material + if (mat % has_sab_table) then + sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic + else + sab_threshold = ZERO + end if + + ! Add contribution from each nuclide in material + do i = 1, mat % n_nuclides + ! Determine microscopic cross sections for this nuclide + index_nuclide = mat % nuclide(i) + + ! Determine whether to use S(a,b) based on energy of particle and whether + ! the nuclide matches the S(a,b) table + if (p % E < sab_threshold .and. i == mat % sab_nuclide) then + index_sab = mat % sab_table + else + index_sab = 0 + end if + + ! Calculate microscopic cross section for this nuclide + if (p % E /= micro_xs(index_nuclide) % last_E) then + call calculate_nuclide_xs(p, index_nuclide, index_sab) + end if + + ! Copy atom density of nuclide in material + atom_density = mat % atom_density(i) + + ! Add contributions to material macroscopic total cross section + material_xs % total = material_xs % total + & + atom_density * micro_xs(index_nuclide) % total + + ! Add contributions to material macroscopic scattering cross section + material_xs % elastic = material_xs % elastic + & + atom_density * micro_xs(index_nuclide) % elastic + + ! Add contributions to material macroscopic absorption cross section + material_xs % absorption = material_xs % absorption + & + atom_density * micro_xs(index_nuclide) % absorption + + ! Add contributions to material macroscopic fission cross section + material_xs % fission = material_xs % fission + & + atom_density * micro_xs(index_nuclide) % fission + + ! Add contributions to material macroscopic nu-fission cross section + material_xs % nu_fission = material_xs % nu_fission + & + atom_density * micro_xs(index_nuclide) % nu_fission + end do + + end subroutine calculate_xs + +!=============================================================================== +! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a +! given index in the nuclides array at the energy of the given particle +!=============================================================================== + + subroutine calculate_nuclide_xs(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + integer, intent(in) :: index_sab ! index into sab_tables array + + integer :: IE ! index on nuclide energy grid + real(8) :: f ! interp factor on nuclide energy grid type(Nuclide), pointer :: nuc => null() + + ! Set pointer to nuclide + nuc => nuclides(index_nuclide) + + ! TODO: If not using unionized energy grid, we need to find the index on the + ! nuclide energy grid using lethargy mapping or whatever other technique + + ! get index on nuclide energy grid + IE = nuc % grid_index(p % IE) + + ! check for rare case where two energy points are the same + if (nuc % energy(IE) == nuc % energy(IE+1)) IE = IE + 1 + + ! calculate interpolation factor + f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) + + micro_xs(index_nuclide) % index_grid = IE + micro_xs(index_nuclide) % interp_factor = f + + ! Initialize sab treatment to false + micro_xs(index_nuclide) % use_sab = .false. + micro_xs(index_nuclide) % elastic_sab = ZERO + + ! Initialize nuclide cross-sections to zero + micro_xs(index_nuclide) % fission = ZERO + micro_xs(index_nuclide) % nu_fission = ZERO + + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % total = & + (ONE-f) * nuc % total(IE) + f * nuc % total(IE+1) + + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % elastic = & + (ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1) + + ! Calculate microscopic nuclide absorption cross section + micro_xs(index_nuclide) % absorption = & + (ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1) + + if (nuc % fissionable) then + ! Calculate microscopic nuclide total cross section + micro_xs(index_nuclide) % fission = & + (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1) + + ! Calculate microscopic nuclide nu-fission cross section + micro_xs(index_nuclide) % nu_fission = & + (ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1) + end if + + ! If there is S(a,b) data for this nuclide, we need to do a few + ! things. Since the total cross section was based on non-S(a,b) data, we + ! need to correct it by subtracting the non-S(a,b) elastic cross section and + ! then add back in the calculated S(a,b) elastic+inelastic cross section. + + if (index_sab > 0) call calculate_sab_xs(p, index_nuclide, index_sab) + + ! if the particle is in the unresolved resonance range and there are + ! probability tables, we need to determine cross sections from the table + + if (urr_ptables_on .and. nuc % urr_present) then + if (p % E > nuc % urr_data % energy(1) .and. & + p % E < nuc % urr_data % energy(nuc % urr_data % n_energy)) then + call calculate_urr_xs(p, index_nuclide) + end if + end if + + ! Set last evaluated energy + micro_xs(index_nuclide) % last_E = p % E + + end subroutine calculate_nuclide_xs + +!=============================================================================== +! CALCULATE_SAB_XS determines the elastic and inelastic scattering +! cross-sections in the thermal energy range. These cross sections replace +! whatever data were taken from the normal Nuclide table. +!=============================================================================== + + subroutine calculate_sab_xs(p, index_nuclide, index_sab) + + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array + integer, intent(in) :: index_sab ! index into sab_tables array + + integer :: IE_sab ! index on S(a,b) energy grid + real(8) :: f_sab ! interp factor on S(a,b) energy grid + real(8) :: inelastic ! S(a,b) inelastic cross section + real(8) :: elastic ! S(a,b) elastic cross section type(SAB_Table), pointer :: sab => null() - ! ========================================================================== - ! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES + ! Set flag that S(a,b) treatment should be used for scattering + micro_xs(index_nuclide) % use_sab = .true. - ! First we need to determine how many continuous-energy tables and how many - ! S(a,b) thermal scattering tables there are -- this loop doesn't actually - ! read the data, it simply counts the number of nuclides and S(a,b) tables - index_nuclides = 0 - index_sab = 0 - do i = 1, n_materials - ! Get pointer to material - mat => materials(i) + ! Get pointer to S(a,b) table + sab => sab_tables(index_sab) - ! First go through all the nuclide and check if they exist on other - ! materials -- if not, then increment the count for the number of - ! nuclides and add to dictionary - do j = 1, mat % n_nuclides - name = mat % names(j) + ! Get index and interpolation factor for inelastic grid + if (p%E < sab % inelastic_e_in(1)) then + IE_sab = 1 + f_sab = ZERO + else + IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) + f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & + (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) + end if - ! First check that this nuclide is listed in the cross_sections.xml - ! file - if (.not. dict_has_key(xs_listing_dict, name)) then - message = "Could not find nuclide " // trim(name) // & - " in cross_sections.xml file!" - call fatal_error() - end if + ! Calculate S(a,b) inelastic scattering cross section + inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & + sab % inelastic_sigma(IE_sab + 1) - ! Find index in xs_listing and set the name and alias according to the - ! listing - index_list = dict_get_key(xs_listing_dict, name) - name = xs_listings(index_list) % name - alias = xs_listings(index_list) % alias + ! Check for elastic data + if (p % E < sab % threshold_elastic) then + ! Determine whether elastic scattering is given in the coherent or + ! incoherent approximation. For coherent, the cross section is + ! represented as P/E whereas for incoherent, it is simply P - ! If this nuclide hasn't been encountered yet, we need to add its name - ! and alias to the nuclide_dict - if (.not. dict_has_key(nuclide_dict, name)) then - index_nuclides = index_nuclides + 1 - mat % nuclide(j) = index_nuclides - - call dict_add_key(nuclide_dict, name, index_nuclides) - call dict_add_key(nuclide_dict, alias, index_nuclides) + if (sab % elastic_mode == SAB_ELASTIC_EXACT) then + if (p % E < sab % elastic_e_in(1)) then + ! If energy is below that of the lowest Bragg peak, the elastic + ! cross section will be zero + elastic = ZERO else - mat % nuclide(j) = dict_get_key(nuclide_dict, name) + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) + elastic = sab % elastic_P(IE_sab) / p % E end if - end do - - ! Check if S(a,b) exists on other materials - if (mat % has_sab_table) then - name = mat % sab_name - - ! First check that this nuclide is listed in the cross_sections.xml - ! file - if (.not. dict_has_key(xs_listing_dict, name)) then - message = "Could not find S(a,b) table " // trim(name) // & - " in cross_sections.xml file!" - call fatal_error() - end if - - ! Find index in xs_listing and set the name and alias according to the - ! listing - index_list = dict_get_key(xs_listing_dict, name) - name = xs_listings(index_list) % name - alias = xs_listings(index_list) % alias - - ! If this S(a,b) table hasn't been encountered yet, we need to add its - ! name and alias to the sab_dict - if (.not. dict_has_key(sab_dict, name)) then - index_sab = index_sab + 1 - mat % sab_table = index_sab - - call dict_add_key(sab_dict, name, index_sab) - call dict_add_key(sab_dict, alias, index_sab) + else + ! Determine index on elastic energy grid + if (p % E < sab % elastic_e_in(1)) then + IE_sab = 1 else - mat % sab_table = dict_get_key(sab_dict, name) + IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) end if + + ! Get interpolation factor for elastic grid + f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & + sab%elastic_e_in(IE_sab)) + + ! Calculate S(a,b) elastic scattering cross section + elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & + sab % elastic_P(IE_sab + 1) end if - end do + else + ! No elastic data + elastic = ZERO + end if - n_nuclides_total = index_nuclides - n_sab_tables = index_sab + ! Correct total and elastic cross sections + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & + micro_xs(index_nuclide) % elastic + inelastic + elastic + micro_xs(index_nuclide) % elastic = inelastic + elastic - ! allocate arrays for ACE table storage and cross section cache - allocate(nuclides(n_nuclides_total)) - allocate(sab_tables(n_sab_tables)) - allocate(micro_xs(n_nuclides_total)) + ! Store S(a,b) elastic cross section for sampling later + micro_xs(index_nuclide) % elastic_sab = elastic - ! ========================================================================== - ! READ ALL ACE CROSS SECTION TABLES - - call dict_create(already_read) - - ! Loop over all files - do i = 1, n_materials - mat => materials(i) - - do j = 1, mat % n_nuclides - name = mat % names(j) - - if (.not. dict_has_key(already_read, name)) then - index_list = dict_get_key(xs_listing_dict, name) - index_nuclides = dict_get_key(nuclide_dict, name) - name = xs_listings(index_list) % name - alias = xs_listings(index_list) % alias - - call read_ace_table(index_nuclides, index_list) - - call dict_add_key(already_read, name, 0) - call dict_add_key(already_read, alias, 0) - end if - end do - - if (mat % has_sab_table) then - name = mat % sab_name - - if (.not. dict_has_key(already_read, name)) then - index_list = dict_get_key(xs_listing_dict, name) - index_sab = dict_get_key(sab_dict, name) - - call read_ace_table(index_sab, index_list) - - call dict_add_key(already_read, name, 0) - end if - end if - end do - - call dict_delete(already_read) - - ! ========================================================================== - ! ASSIGN S(A,B) TABLES TO SPECIFIC NUCLIDES WITHIN MATERIALS - - do i = 1, n_materials - mat => materials(i) - - if (mat % has_sab_table) then - ! In order to know which nuclide the S(a,b) table applies to, we need - ! to search through the list of nuclides for one which has a matching - ! zaid - sab => sab_tables(mat % sab_table) - - do j = 1, mat % n_nuclides - nuc => nuclides(mat % nuclide(j)) - if (nuc % zaid == sab % zaid) then - mat % sab_nuclide = j - end if - end do - - ! Check to make sure S(a,b) table matched a nuclide - if (mat % sab_nuclide == 0) then - message = "S(a,b) table " // trim(mat % sab_name) // " did not match " & - // "any nuclide on material " // trim(int_to_str(mat % id)) - call fatal_error() - end if - end if - end do - - end subroutine read_xs + end subroutine calculate_sab_xs !=============================================================================== -! READ_ACE_BINARY reads a single cross section table in binary format. This -! routine reads the header data for each table and then calls appropriate -! subroutines to parse the actual data. +! CALCULATE_URR_XS determines cross sections in the unresolved resonance range +! from probability tables !=============================================================================== - subroutine read_ace_table(index_table, index_list) + subroutine calculate_urr_xs(p, index_nuclide) - integer, intent(in) :: index_table ! index in nuclides/sab_tables - integer, intent(in) :: index_list ! index in xs_listings + type(Particle), pointer :: p + integer, intent(in) :: index_nuclide ! index into nuclides array - integer :: i ! loop index for XSS records - integer :: j, j1, j2 ! indices in XSS - integer :: record_length ! Fortran record length - integer :: location ! location of ACE table - integer :: entries ! number of entries on each record - integer :: length ! length of ACE table - integer :: in = 7 ! file unit - integer :: zaids(16) ! list of ZAIDs (only used for S(a,b)) - integer :: filetype ! filetype (ASCII or BINARY) - real(8) :: kT ! temperature of table - real(8) :: awrs(16) ! list of atomic weight ratios (not used) - real(8) :: awr ! atomic weight ratio for table - logical :: file_exists ! does ACE library exist? - character(7) :: readable ! is ACE library readable? - character(10) :: name ! name of ACE table - character(10) :: date_ ! date ACE library was processed - character(10) :: mat ! material identifier - character(70) :: comment ! comment for ACE table - character(MAX_FILE_LEN) :: filename ! path to ACE cross section library - type(Nuclide), pointer :: nuc => null() - type(SAB_Table), pointer :: sab => null() - type(XsListing), pointer :: listing => null() + integer :: i_energy ! index for energy + integer :: i_table ! index for table + real(8) :: f ! interpolation factor + real(8) :: r ! pseudo-random number + real(8) :: elastic ! smooth elastic cross section + real(8) :: absorption ! smooth absorption cross section + real(8) :: fission ! smooth fission cross section + real(8) :: inelastic ! smooth inelastic cross section + type(UrrData), pointer :: urr => null() + type(Nuclide), pointer :: nuc => null() + type(Reaction), pointer :: rxn => null() - ! determine path, record length, and location of table - listing => xs_listings(index_list) - filename = listing % path - record_length = listing % recl - location = listing % location - entries = listing % entries - filetype = listing % filetype + ! copy cross-sections already calculated + elastic = micro_xs(index_nuclide) % elastic + absorption = micro_xs(index_nuclide) % absorption + fission = micro_xs(index_nuclide) % fission - ! Check if ACE library exists and is readable - inquire(FILE=filename, EXIST=file_exists, READ=readable) - if (.not. file_exists) then - message = "ACE library '" // trim(filename) // "' does not exist!" - call fatal_error() - elseif (readable(1:3) == 'NO') then - message = "ACE library '" // trim(filename) // "' is not readable! & - &Change file permissions with chmod command." + ! get pointer to probability table + nuc => nuclides(index_nuclide) + urr => nuc % urr_data + + ! determine energy table + i_energy = 1 + do + if (p % E < urr % energy(i_energy + 1)) exit + i_energy = i_energy + 1 + end do + + ! determine interpolation factor on table + f = (p % E - urr % energy(i_energy)) / & + (urr % energy(i_energy + 1) - urr % energy(i_energy)) + + ! sample probability table using the cumulative distribution + r = prn() + i_table = 1 + do + if (urr % prob(i_energy, URR_CUM_PROB, i_table) > r) exit + i_table = i_table + 1 + end do + + ! determine elastic, fission, and capture cross sections from probability + ! table + if (urr % interp == LINEAR_LINEAR) then + micro_xs(index_nuclide) % elastic = & + (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_table) + & + f * urr % prob(i_energy + 1, URR_ELASTIC, i_table) + micro_xs(index_nuclide) % fission = & + (ONE - f) * urr % prob(i_energy, URR_FISSION, i_table) + & + f * urr % prob(i_energy + 1, URR_FISSION, i_table) + micro_xs(index_nuclide) % absorption = & + micro_xs(index_nuclide) % fission + & + (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_table) + & + f * urr % prob(i_energy + 1, URR_N_GAMMA, i_table) + elseif (urr % interp == LOG_LOG) then + message = "Log-log interpolation on probability table not yet supported." call fatal_error() end if - ! display message - message = "Loading ACE cross section table: " // listing % name - call write_message(6) + ! Determine treatment of inelastic scattering + inelastic = ZERO + if (urr % inelastic_flag > 0) then + ! Get pointer to inelastic scattering reaction + rxn => nuc % reactions(nuc % urr_inelastic) - if (filetype == ASCII) then - ! ======================================================================= - ! READ ACE TABLE IN ASCII FORMAT + ! Get index on energy grid and interpolation factor + i_energy = micro_xs(index_nuclide) % index_grid + f = micro_xs(index_nuclide) % interp_factor - ! Find location of table - open(UNIT=in, FILE=filename, STATUS='old', ACTION='read') - rewind(UNIT=in) - do i = 1, location - 1 - read(UNIT=in, FMT=*) - end do - - ! Read first line of header - read(UNIT=in, FMT='(A10,2E12.0,1X,A10)') name, awr, kT, date_ - - ! Read more header and NXS and JXS - 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) - - ! determine table length - length = NXS(1) - allocate(XSS(length)) - - ! Read XSS array - read(UNIT=in, FMT='(4E20.0)') XSS - - elseif (filetype == BINARY) then - ! ======================================================================= - ! READ ACE TABLE IN BINARY FORMAT - - open(UNIT=in, FILE=filename, STATUS='old', ACTION='read', & - ACCESS='direct', RECL=record_length) - - ! Read all header information - read(UNIT=in, REC=location) name, awr, kT, date_, & - comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS - - ! determine table length - length = NXS(1) - allocate(XSS(length)) - - ! Read remaining records with XSS - do i = 1, (length + entries - 1)/entries - j1 = 1 + (i-1)*entries - j2 = min(length, j1 + entries - 1) - read(UNIT=IN, REC=location + i) (XSS(j), j=j1,j2) - end do + ! Determine inelastic scattering cross section + if (i_energy >= rxn % IE) then + inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%IE + 1) + & + f * rxn % sigma(i_energy - rxn%IE + 2) + end if end if - ! ========================================================================== - ! PARSE DATA BASED ON NXS, JXS, AND XSS ARRAYS - - select case(listing % type) - case (ACE_NEUTRON) - nuc => nuclides(index_table) - nuc % name = name - nuc % awr = awr - nuc % kT = kT - nuc % zaid = NXS(2) - - ! read all blocks - call read_esz(nuc) - call read_nu_data(nuc) - call read_reactions(nuc) - call read_angular_dist(nuc) - call read_energy_dist(nuc) - call read_unr_res(nuc) - - ! for fissionable nuclides, precalculate microscopic nu-fission cross - ! sections so that we don't need to call the nu_total function during - ! cross section lookups - - if (nuc % fissionable) call generate_nu_fission(nuc) - - case (ACE_THERMAL) - sab => sab_tables(index_table) - sab % name = name - sab % awr = awr - sab % kT = kT - sab % zaid = zaids(1) - - call read_thermal_data(sab) - end select - - deallocate(XSS) - if(associated(nuc)) nullify(nuc) - if(associated(sab)) nullify(sab) - - close(UNIT=in) - - end subroutine read_ace_table - -!=============================================================================== -! READ_ESZ - reads through the ESZ block. This block contains the energy grid, -! total xs, absorption xs, elastic scattering xs, and heating numbers. -!=============================================================================== - - subroutine read_esz(nuc) - - type(Nuclide), pointer :: nuc - - integer :: NE ! number of energy points for total and elastic cross sections - - ! determine number of energy points - NE = NXS(3) - nuc % n_grid = NE - - ! allocate storage for energy grid and cross section arrays - allocate(nuc % energy(NE)) - allocate(nuc % total(NE)) - allocate(nuc % elastic(NE)) - allocate(nuc % fission(NE)) - allocate(nuc % nu_fission(NE)) - allocate(nuc % absorption(NE)) - allocate(nuc % heating(NE)) - - ! initialize cross sections - nuc % total = ZERO - nuc % elastic = ZERO - nuc % fission = ZERO - nuc % nu_fission = ZERO - nuc % absorption = ZERO - nuc % heating = ZERO - - ! Read data from XSS -- only the energy grid, elastic scattering and heating - ! cross section values are actually read from here. The total and absorption - ! cross sections are reconstructed from the partial reaction data. - - XSS_index = 1 - nuc % energy = get_real(NE) - - ! Skip total and absorption - XSS_index = XSS_index + 2*NE - - ! Continue reading elastic scattering and heating - nuc % elastic = get_real(NE) - nuc % heating = get_real(NE) - - end subroutine read_esz - -!=============================================================================== -! 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 read_nu_data(nuc) - - type(Nuclide), pointer :: nuc - - integer :: i ! loop index - integer :: JXS2 ! location for fission nu data - integer :: JXS24 ! location for delayed neutron data - 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 :: NPCR ! number of delayed neutron precursor groups - integer :: LED ! location of energy distribution locators - integer :: LDIS ! location of all energy distributions - integer :: LOCC ! location of energy distributions for given MT - integer :: lc ! locator - integer :: length ! length of data to allocate - type(DistEnergy), pointer :: edist => null() - - JXS2 = JXS(2) - JXS24 = JXS(24) - - if (JXS2 == 0) then - ! ======================================================================= - ! NO PROMPT/TOTAL NU DATA - nuc % nu_t_type = NU_NONE - nuc % nu_p_type = NU_NONE - - elseif (XSS(JXS2) > 0) then - ! ======================================================================= - ! PROMPT OR TOTAL NU DATA - KNU = JXS2 - LNU = int(XSS(KNU)) - if (LNU == 1) then - ! Polynomial data - nuc % nu_t_type = NU_POLYNOMIAL - nuc % 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 - nuc % nu_t_type = NU_TABULAR - nuc % 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(nuc % 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 - nuc % nu_t_data = get_real(length) - - elseif (XSS(JXS2) < 0) then - ! ======================================================================= - ! PROMPT AND TOTAL NU DATA -- read prompt data first - KNU = JXS2 + 1 - LNU = int(XSS(KNU)) - if (LNU == 1) then - ! Polynomial data - nuc % nu_p_type = NU_POLYNOMIAL - - ! allocate determine how many coefficients for polynomial - NC = int(XSS(KNU+1)) - length = NC + 1 - elseif (LNU == 2) then - ! Tabular data - nuc % nu_p_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(nuc % nu_p_data(length)) - - ! read data - XSS_index = KNU + 1 - nuc % nu_p_data = get_real(length) - - ! Now read total nu data - KNU = JXS2 + int(abs(XSS(JXS2))) + 1 - LNU = int(XSS(KNU)) - if (LNU == 1) then - ! Polynomial data - nuc % 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 - nuc % 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(nuc % nu_t_data(length)) - - ! read data - XSS_index = KNU + 1 - nuc % nu_t_data = get_real(length) + ! Multiply by smooth cross-section if needed + if (urr % multiply_smooth) then + micro_xs(index_nuclide) % elastic = elastic * & + micro_xs(index_nuclide) % elastic + micro_xs(index_nuclide) % absorption = absorption * & + micro_xs(index_nuclide) % absorption + micro_xs(index_nuclide) % fission = fission * & + micro_xs(index_nuclide) % fission end if - if (JXS24 > 0) then - ! ======================================================================= - ! DELAYED NU DATA + ! Determine nu-fission cross section + if (nuc % fissionable) then + micro_xs(index_nuclide) % nu_fission = nu_total(nuc, p % E) * & + micro_xs(index_nuclide) % fission + end if - nuc % nu_d_type = NU_TABULAR - KNU = JXS24 + ! Calculate total cross section as sum of partials + micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % elastic + & + inelastic + micro_xs(index_nuclide) % absorption - ! 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 + end subroutine calculate_urr_xs - ! allocate space for delayed nu data - allocate(nuc % nu_d_data(length)) - - ! read delayed nu data - XSS_index = KNU + 1 - nuc % nu_d_data = get_real(length) - - ! ======================================================================= - ! DELAYED NEUTRON ENERGY DISTRIBUTION +!=============================================================================== +! FIND_ENERGY_INDEX determines the index on the union energy grid and the +! interpolation factor for a particle at a certain energy +!=============================================================================== - ! Allocate space for secondary energy distribution - NPCR = NXS(8) - nuc % n_precursor = NPCR - allocate(nuc % nu_d_edist(NPCR)) + subroutine find_energy_index(p) - LED = JXS(26) - LDIS = JXS(27) + type(Particle), pointer :: p - ! Loop over all delayed neutron precursor groups - do i = 1, NPCR - ! find location of energy distribution data - LOCC = int(XSS(LED + i - 1)) + integer :: IE ! index on union energy grid + real(8) :: E ! energy of particle + real(8) :: interp ! interpolation factor - ! read energy distribution data - edist => nuc % nu_d_edist(i) - call get_energy_dist(edist, LOCC, .true.) - end do - - ! ======================================================================= - ! DELAYED NEUTRON PRECUSOR YIELDS AND CONSTANTS - - ! determine length of all precursor constants/yields/interp data - length = 0 - lc = JXS(25) - do i = 1, NPCR - NR = int(XSS(lc + length + 1)) - NE = int(XSS(lc + length + 2 + 2*NR)) - length = length + 3 + 2*NR + 2*NE - end do - - ! allocate space for precusor data - allocate(nuc % nu_d_precursor_data(length)) - - ! read delayed neutron precursor data - XSS_index = lc - nuc % nu_d_precursor_data = get_real(length) + ! copy particle's energy + E = p % E + ! if particle's energy is outside of energy grid range, set to first or last + ! index. Otherwise, do a binary search through the union energy grid. + if (E < e_grid(1)) then + IE = 1 + elseif (E > e_grid(n_grid)) then + IE = n_grid - 1 else - nuc % 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(nuc) - - type(Nuclide), pointer :: nuc - - integer :: i ! loop indices - integer :: i_fission ! index in nuc % index_fission - integer :: LMT ! index of MT list in XSS - integer :: NMT ! Number of reactions - integer :: JXS4 ! index of Q values in XSS - integer :: JXS5 ! index of neutron multiplicities in XSS - integer :: JXS7 ! index of reactions cross-sections in XSS - integer :: LXS ! location of cross-section locators - integer :: LOCA ! location of cross-section for given MT - integer :: IE ! reaction's starting index on energy grid - integer :: NE ! number of energies for reaction - type(Reaction), pointer :: rxn => null() - - LMT = JXS(3) - JXS4 = JXS(4) - JXS5 = JXS(5) - LXS = JXS(6) - JXS7 = JXS(7) - NMT = NXS(4) - - ! allocate array of reactions. Add one since we need to include an elastic - ! scattering channel - nuc % n_reaction = NMT + 1 - allocate(nuc % reactions(NMT+1)) - - ! Store elastic scattering cross-section on reaction one - rxn => nuc % reactions(1) - rxn % MT = 2 - rxn % Q_value = ZERO - rxn % TY = 1 - rxn % IE = 1 - allocate(rxn % sigma(nuc % n_grid)) - rxn % sigma = nuc % elastic - - ! Add contribution of elastic scattering to total cross section - nuc % total = nuc % total + nuc % elastic - - ! By default, set nuclide to not fissionable and then change if fission - ! reactions are encountered - nuc % fissionable = .false. - nuc % has_partial_fission = .false. - nuc % n_fission = 0 - i_fission = 0 - - do i = 1, NMT - rxn => nuc % reactions(i+1) - - ! read MT number, Q-value, and neutrons produced - rxn % MT = int(XSS(LMT + i - 1)) - rxn % Q_value = XSS(JXS4 + i - 1) - rxn % TY = int(XSS(JXS5 + i - 1)) - - ! read starting energy index - LOCA = int(XSS(LXS + i - 1)) - IE = int(XSS(JXS7 + LOCA - 1)) - rxn % IE = IE - - ! read number of energies cross section values - NE = int(XSS(JXS7 + LOCA)) - allocate(rxn % sigma(NE)) - XSS_index = JXS7 + LOCA + 1 - rxn % sigma = get_real(NE) - - ! Skip redundant reactions -- this includes total inelastic level - ! scattering, gas production cross sections (MT=200+), and (n,p), (n,d), - ! etc. reactions leaving the nucleus in an excited state - if (rxn % MT == N_LEVEL .or. rxn % MT > N_DA) cycle - - ! Add contribution to total cross section - nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma - - ! Add contribution to absorption cross section - if (rxn % MT >= N_GAMMA .and. rxn % MT <= N_DA) then - nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma - end if - - ! Information about fission reactions - if (rxn % MT == N_FISSION) then - allocate(nuc % index_fission(1)) - elseif (rxn % MT == N_F) then - allocate(nuc % index_fission(PARTIAL_FISSION_MAX)) - nuc % has_partial_fission = .true. - end if - - ! Add contribution to fission cross section - if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF & - .or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) then - nuc % fissionable = .true. - nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma - - ! Also need to add fission cross sections to absorption - nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma - - ! Keep track of this reaction for easy searching later - i_fission = i_fission + 1 - nuc % index_fission(i_fission) = i + 1 - nuc % n_fission = nuc % n_fission + 1 - end if - - ! set defaults - rxn % has_angle_dist = .false. - rxn % has_energy_dist = .false. - end do - - end subroutine read_reactions - -!=============================================================================== -! READ_ANGULAR_DIST parses the angular distribution for each reaction with -! secondary neutrons -!=============================================================================== - - subroutine read_angular_dist(nuc) - - type(Nuclide), pointer :: nuc - - integer :: JXS8 ! location of angular distribution locators - integer :: JXS9 ! location of angular distributions - integer :: LOCB ! location of angular distribution for given MT - integer :: NE ! number of incoming energies - integer :: NP ! number of points for cosine distribution - integer :: LC ! locator - integer :: i ! index in reactions array - integer :: j ! index over incoming energies - integer :: length ! length of data array to allocate - type(Reaction), pointer :: rxn => null() - - JXS8 = JXS(8) - JXS9 = JXS(9) - - ! loop over all reactions with secondary neutrons -- NXS(5) does not include - ! elastic scattering - do i = 1, NXS(5) + 1 - rxn => nuc%reactions(i) - - ! find location of angular distribution - LOCB = int(XSS(JXS8 + i - 1)) - if (LOCB == -1) then - ! Angular distribution data are specified through LAWi = 44 in the DLW - ! block - cycle - elseif (LOCB == 0) then - ! No angular distribution data are given for this reaction, isotropic - ! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0) - cycle - end if - rxn % has_angle_dist = .true. - - ! allocate space for incoming energies and locations - NE = int(XSS(JXS9 + LOCB - 1)) - 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 nucs - XSS_index = JXS9 + LOCB - rxn % adist % energy = get_real(NE) - rxn % adist % location = get_int(NE) - - ! determine dize of data block - length = 0 - do j = 1, NE - LC = rxn % adist % location(j) - if (LC == 0) then - ! isotropic - rxn % adist % type(j) = ANGLE_ISOTROPIC - elseif (LC > 0) then - ! 32 equiprobable bins - rxn % adist % type(j) = ANGLE_32_EQUI - length = length + 33 - elseif (LC < 0) then - ! tabular distribution - rxn % adist % type(j) = ANGLE_TABULAR - NP = int(XSS(JXS9 + abs(LC))) - length = length + 2 + 3*NP - end if - end do - - ! allocate angular distribution data and read - 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) - XSS_index = JXS9 + abs(LC) - 1 - 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 - - end do - - end subroutine read_angular_dist - -!=============================================================================== -! READ_ENERGY_DIST parses the secondary energy distribution for each reaction -! with seconary neutrons (except elastic scattering) -!=============================================================================== - - subroutine read_energy_dist(nuc) - - type(Nuclide), pointer :: nuc - - integer :: LED ! location of energy distribution locators - integer :: LOCC ! location of energy distributions for given MT - integer :: i ! loop index - type(Reaction), pointer :: rxn => null() - - LED = JXS(10) - - ! Loop over all reactions - do i = 1, NXS(5) - rxn => nuc % reactions(i+1) ! skip over elastic scattering - rxn % has_energy_dist = .true. - - ! find location of energy distribution data - LOCC = int(XSS(LED + i - 1)) - - ! allocate energy distribution - allocate(rxn % edist) - - ! read data for energy distribution - call get_energy_dist(rxn % edist, LOCC) - end do - - end subroutine read_energy_dist - -!=============================================================================== -! GET_ENERGY_DIST reads in data for a single law for an energy distribution and -! calls itself recursively if there are multiple energy distributions for a -! single reaction -!=============================================================================== - - recursive subroutine get_energy_dist(edist, loc_law, delayed_n) - - type(DistEnergy), pointer :: edist ! energy distribution - integer, intent(in) :: loc_law ! locator for data - logical, optional :: delayed_n ! is this for delayed neutrons? - - integer :: LDIS ! location of all energy distributions - integer :: LNW ! location of next energy distribution if multiple - 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 - integer :: lc ! locator - integer :: length ! length of data to allocate - integer :: length_interp_data ! length of interpolation data - - ! determine location of energy distribution - if (present(delayed_n)) then - LDIS = JXS(27) - else - LDIS = JXS(11) - end if - - ! locator for next law and information on this law - LNW = int(XSS(LDIS + loc_law - 1)) - LAW = int(XSS(LDIS + loc_law)) - IDAT = int(XSS(LDIS + loc_law + 1)) - NR = int(XSS(LDIS + loc_law + 2)) - edist % law = LAW - edist % p_valid % n_regions = NR - - ! allocate space for ENDF interpolation parameters - if (NR > 0) then - allocate(edist % p_valid % nbt(NR)) - allocate(edist % p_valid % int(NR)) - end if - - ! read ENDF interpolation parameters - XSS_index = LDIS + loc_law + 3 - if (NR > 0) then - edist % p_valid % nbt = int(get_real(NR)) - edist % p_valid % int = int(get_real(NR)) - end if - - ! allocate space for law validity data - NE = int(XSS(LDIS + loc_law + 3 + 2*NR)) - edist % p_valid % n_pairs = NE - allocate(edist % p_valid % x(NE)) - allocate(edist % p_valid % y(NE)) - - length_interp_data = 5 + 2*(NR + NE) - - ! read law validity data - XSS_index = LDIS + loc_law + 4 + 2*NR - edist % p_valid % x = get_real(NE) - edist % p_valid % y = get_real(NE) - - ! Set index to beginning of IDAT array - lc = LDIS + IDAT - 2 - - ! determine length of energy distribution - length = length_energy_dist(lc, LAW, loc_law, length_interp_data) - - ! allocate secondary energy distribution array - allocate(edist % data(length)) - - ! read secondary energy distribution - XSS_index = lc + 1 - edist % data = get_real(length) - - ! read next energy distribution if present - if (LNW > 0) then - allocate(edist % next) - call get_energy_dist(edist % next, LNW) - end if - - end subroutine get_energy_dist - -!=============================================================================== -! LENGTH_ENERGY_DIST determines how many values are contained in an LDAT energy -! distribution array based on the secondary energy law and location in XSS -!=============================================================================== - - function length_energy_dist(lc, law, LOCC, lid) result(length) - - integer, intent(in) :: lc ! location in XSS array - integer, intent(in) :: law ! energy distribution law - integer, intent(in) :: LOCC ! location of energy distribution - integer, intent(in) :: lid ! length of interpolation data - integer :: length ! length of energy distribution (LDAT) - - integer :: i,j,k - integer :: NR ! number of interpolation regions - integer :: NE ! number of incoming energies - integer :: NP ! number of points in outgoing energy distribution - integer :: NMU ! number of points in outgoing cosine distribution - integer :: NRa ! number of interpolation regions for Watt 'a' - integer :: NEa ! number of energies for Watt 'a' - integer :: NRb ! number of interpolation regions for Watt 'b' - integer :: NEb ! number of energies for Watt 'b' - - ! initialize length - length = 0 - - select case (law) - case (1) - ! Tabular equiprobable energy bins - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - NP = int(XSS(lc + 3 + 2*NR + NE)) - length = 3 + 2*NR + NE + 3*NP*NE - - case (2) - ! Discrete photon energy - length = 2 - - case (3) - ! Level scattering - length = 2 - - case (4) - ! Continuous tabular distribution - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - length = length + 2 + 2*NR + 2*NE - do i = 1,NE - ! determine length - NP = int(XSS(lc + length + 2)) - length = length + 2 + 3*NP - - ! adjust location for this block - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid - end do - - case (5) - ! General evaporation spectrum - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - NP = int(XSS(lc + 3 + 2*NR + 2*NE)) - length = 3 + 2*NR + 2*NE + NP - - case (7) - ! Maxwell fission spectrum - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - length = 3 + 2*NR + 2*NE - - case (9) - ! Evaporation spectrum - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - length = 3 + 2*NR + 2*NE - - case (11) - ! Watt spectrum - NRa = int(XSS(lc + 1)) - NEa = int(XSS(lc + 2 + 2*NRa)) - NRb = int(XSS(lc + 3 + 2*(NRa+NEa))) - NEb = int(XSS(lc + 4 + 2*(NRa+NEa+NRb))) - length = 5 + 2*(NRa + NEa + NRb + NEb) - - case (44) - ! Kalbach-Mann correlated scattering - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - length = length + 2 + 2*NR + 2*NE - do i = 1,NE - NP = int(XSS(lc + length + 2)) - length = length + 2 + 5*NP - - ! adjust location for this block - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid - end do - - case (61) - ! Correlated energy and angle distribution - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - length = length + 2 + 2*NR + 2*NE - do i = 1,NE - ! outgoing energy distribution - NP = int(XSS(lc + length + 2)) - - ! adjust locators for angular distribution - do j = 1, NP - k = lc + length + 2 + 3*NP + j - if (XSS(k) /= 0) XSS(k) = XSS(k) - LOCC - lid - end do - - length = length + 2 + 4*NP - do j = 1, NP - ! outgoing angle distribution -- NMU here is actually - ! referred to as NP in the MCNP documentation - NMU = int(XSS(lc + length + 2)) - length = length + 2 + 3*NMU - end do - - ! adjust locators for energy distribution - j = lc + 2 + 2*NR + NE + i - XSS(j) = XSS(j) - LOCC - lid - end do - - case (66) - ! N-body phase space distribution - length = 2 - - case (67) - ! Laboratory energy-angle law - NR = int(XSS(lc + 1)) - NE = int(XSS(lc + 2 + 2*NR)) - NMU = int(XSS(lc + 4 + 2*NR + 2*NE)) - length = 4 + 2*(NR + NE + NMU) - - end select - - end function length_energy_dist - -!=============================================================================== -! READ_UNR_RES reads in unresolved resonance probability tables if present. -!=============================================================================== - - subroutine read_unr_res(nuc) - - type(Nuclide), pointer :: nuc - - integer :: JXS23 ! location of URR data - integer :: lc ! locator - integer :: N ! # of incident energies - integer :: M ! # of probabilities - integer :: i ! index over incoming energies - integer :: j ! index over values - integer :: k ! index over probabilities - - ! determine locator for URR data - JXS23 = JXS(23) - - ! check if URR data is present - if (JXS23 /= 0) then - nuc % urr_present = .true. - allocate(nuc % urr_data) - allocate(nuc % urr_data % params(6)) - lc = JXS23 - else - nuc % urr_present = .false. - return - end if - - ! read parameters - nuc % urr_data % params(1) = int(XSS(lc)) ! # of incident energies - nuc % urr_data % params(2) = int(XSS(lc + 1)) ! # of probabilities - nuc % urr_data % params(3) = int(XSS(lc + 2)) ! interpolation parameter - nuc % urr_data % params(4) = int(XSS(lc + 3)) ! inelastic competition flag - nuc % urr_data % params(5) = int(XSS(lc + 4)) ! other absorption flag - nuc % urr_data % params(6) = int(XSS(lc + 5)) ! factors flag - - ! allocate incident energies and probability tables - N = nuc % urr_data % params(1) - M = nuc % urr_data % params(2) - allocate(nuc % urr_data % energy(N)) - allocate(nuc % urr_data % prob(N,6,M)) - - ! read incident energies - XSS_index = lc + 6 - nuc % urr_data % energy = get_real(N) - - ! read probability tables - do i = 1, N - do j = 1, 6 - do k = 1, M - nuc % urr_data % prob(i,j,k) = XSS(XSS_index) - XSS_index = XSS_index + 1 - end do - end do - end do - - end subroutine read_unr_res - -!=============================================================================== -! GENERATE_NU_FISSION precalculates the microscopic nu-fission cross section for -! a given nuclide. This is done so that the nu_total function does not need to -! be called during cross section lookups. -!=============================================================================== - - subroutine generate_nu_fission(nuc) - - type(Nuclide), pointer :: nuc - - integer :: i ! index on nuclide energy grid - real(8) :: E ! energy - real(8) :: nu ! # of neutrons per fission - - do i = 1, nuc % n_grid - ! determine energy - E = nuc % energy(i) - - ! determine total nu at given energy - nu = nu_total(nuc, E) - - ! determine nu-fission microscopic cross section - nuc % nu_fission(i) = nu * nuc % fission(i) - end do - - end subroutine generate_nu_fission - -!=============================================================================== -! READ_THERMAL_DATA reads elastic and inelastic cross sections and corresponding -! secondary energy/angle distributions derived from experimental S(a,b) -! data. Namely, in MCNP language, this routine reads the ITIE, ITCE, ITXE, and -! ITCA blocks. -!=============================================================================== - - subroutine read_thermal_data(table) - - type(SAB_Table), pointer :: table - - integer :: i ! index for incoming energies - integer :: j ! index for outgoing energies - integer :: k ! index for outoging angles - integer :: lc ! location in XSS array - integer :: NE_in ! number of incoming energies - integer :: NE_out ! number of outgoing energies - integer :: NMU ! number of outgoing angles - integer :: JXS4 ! location of elastic energy table - - ! read secondary energy mode for inelastic scattering - table % secondary_mode = NXS(7) - - ! read number of inelastic energies and allocate arrays - NE_in = int(XSS(JXS(1))) - table % n_inelastic_e_in = NE_in - allocate(table % inelastic_e_in(NE_in)) - allocate(table % inelastic_sigma(NE_in)) - - ! read inelastic energies and cross-sections - XSS_index = JXS(1) + 1 - table % inelastic_e_in = get_real(NE_in) - table % inelastic_sigma = get_real(NE_in) - - ! set threshold value - table % threshold_inelastic = table % inelastic_e_in(NE_in) - - ! 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)) - - ! 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) - - ! 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 - end do - - ! read number of elastic energies and allocate arrays - JXS4 = JXS(4) - if (JXS4 /= 0) then - NE_in = int(XSS(JXS4)) - table % n_elastic_e_in = NE_in - allocate(table % elastic_e_in(NE_in)) - allocate(table % elastic_P(NE_in)) - - ! read elastic energies and P - XSS_index = JXS4 + 1 - table % elastic_e_in = get_real(NE_in) - table % elastic_P = get_real(NE_in) - - ! set threshold - table % threshold_elastic = table % elastic_e_in(NE_in) - - ! determine whether sigma=P or sigma = P/E - table % elastic_mode = NXS(5) - else - table % threshold_elastic = ZERO - table % n_elastic_e_in = 0 - end if - - ! allocate space for outgoing energy/angle for elastic scattering - NMU = NXS(6) + 1 - table % n_elastic_mu = NMU - if (NMU > 0) then - allocate(table % elastic_mu(NMU, NE_in)) - end if - - ! read equiprobable outgoing cosines for elastic scattering each - ! incoming energy - if (JXS4 /= 0 .and. NMU /= 0) then - lc = JXS(6) - 1 - do i = 1, NE_in - do j = 1, NMU - table % elastic_mu(j,i) = XSS(lc + j) - end do - lc = lc + NMU - end do + IE = binary_search(e_grid, n_grid, E) end if - end subroutine read_thermal_data + ! calculate the interpolation factor -- note this will be outside of [0,1) + ! for a particle outside the energy range of the union grid + interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) -!=============================================================================== -! GET_INT returns an array of integers read from the current position in the XSS -! array -!=============================================================================== - - function get_int(n_values) result(array) - - integer, intent(in) :: n_values ! number of values to read - integer :: array(n_values) ! array of values - - array = int(XSS(XSS_index:XSS_index + n_values - 1)) - XSS_index = XSS_index + n_values - - end function get_int - -!=============================================================================== -! GET_REAL returns an array of real(8)s read from the current position in the -! XSS array -!=============================================================================== - - function get_real(n_values) result(array) - - integer, intent(in) :: n_values ! number of values to read - real(8) :: array(n_values) ! array of values - - array = XSS(XSS_index:XSS_index + n_values - 1) - XSS_index = XSS_index + n_values - - end function get_real + ! set particle attributes + p % IE = IE + p % interp = interp + + end subroutine find_energy_index end module cross_section diff --git a/src/endf.F90 b/src/endf.F90 index 5f35da5b07..92103b5c7e 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -1,7 +1,7 @@ module endf use constants - use string, only: int_to_str + use string, only: to_str contains @@ -72,7 +72,7 @@ contains case (N_NPA) string = '(n,npa)' case (N_N1 : N_N40) - string = '(n,n' // trim(int_to_str(MT-50)) // ')' + string = '(n,n' // trim(to_str(MT-50)) // ')' case (N_NC) string = '(n,nc)' case (N_GAMMA) @@ -122,27 +122,27 @@ contains case (444) string = '(damage)' case (600 : 648) - string = '(n,p' // trim(int_to_str(MT-600)) // ')' + string = '(n,p' // trim(to_str(MT-600)) // ')' case (649) string = '(n,pc)' case (650 : 698) - string = '(n,d' // trim(int_to_str(MT-650)) // ')' + string = '(n,d' // trim(to_str(MT-650)) // ')' case (699) string = '(n,dc)' case (700 : 748) - string = '(n,t' // trim(int_to_str(MT-700)) // ')' + string = '(n,t' // trim(to_str(MT-700)) // ')' case (749) string = '(n,tc)' case (750 : 798) - string = '(n,3He' // trim(int_to_str(MT-750)) // ')' + string = '(n,3He' // trim(to_str(MT-750)) // ')' case (799) string = '(n,3Hec)' case (800 : 848) - string = '(n,a' // trim(int_to_str(MT-800)) // ')' + string = '(n,a' // trim(to_str(MT-800)) // ')' case (849) string = '(n,tc)' case default - string = 'MT=' // trim(int_to_str(MT)) + string = 'MT=' // trim(to_str(MT)) end select end function reaction_name diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index d2351c94cf..39f5fcd31f 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -99,7 +99,7 @@ contains ! if the original list is empty, we need to allocate the first element and ! store first energy point - if (list_size(list) == 0) then + if (.not. associated(list)) then allocate(list) current => list do i = 1, n diff --git a/src/error.F90 b/src/error.F90 index fab7c8dc84..2eac3f2372 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -2,7 +2,11 @@ module error use ISO_FORTRAN_ENV - use global, only: master, free_memory, message + use global, only: master, free_memory, message, mpi_err + +#ifdef MPI + use mpi +#endif implicit none @@ -25,14 +29,14 @@ contains ! Only allow master to print to screen if (.not. master) return - write(ou, fmt='(1X,A9)', advance='no') 'WARNING: ' + write(OUTPUT_UNIT, fmt='(1X,A9)', advance='no') 'WARNING: ' n_lines = (len_trim(message)-1)/70 + 1 do i = 1, n_lines if (i == 1) then - write(ou, fmt='(A70)') message(70*(i-1)+1:70*i) + write(OUTPUT_UNIT, fmt='(A70)') message(70*(i-1)+1:70*i) else - write(ou, fmt='(10X,A70)') message(70*(i-1)+1:70*i) + write(OUTPUT_UNIT, fmt='(10X,A70)') message(70*(i-1)+1:70*i) end if end do @@ -44,29 +48,42 @@ contains ! the program is aborted. !=============================================================================== - subroutine fatal_error() + subroutine fatal_error(error_code) + integer, optional :: error_code ! error code + + integer :: code ! error code integer :: n_lines ! number of lines integer :: i ! loop index over lines - ! Only allow master to print to screen - if (master) then - write(eu, fmt='(1X,A7)', advance='no') 'ERROR: ' - - n_lines = (len_trim(message)-1)/72 + 1 - do i = 1, n_lines - if (i == 1) then - write(eu, fmt='(A72)') message(72*(i-1)+1:72*i) - else - write(eu, fmt='(7X,A72)') message(72*(i-1)+1:72*i) - end if - end do - write(eu,*) + ! set default error code + if (present(error_code)) then + code = error_code + else + code = -1 end if + ! Only allow master to print to screen + write(ERROR_UNIT, fmt='(1X,A7)', advance='no') 'ERROR: ' + + n_lines = (len_trim(message)-1)/72 + 1 + do i = 1, n_lines + if (i == 1) then + write(ERROR_UNIT, fmt='(A72)') message(72*(i-1)+1:72*i) + else + write(ERROR_UNIT, fmt='(7X,A72)') message(72*(i-1)+1:72*i) + end if + end do + write(ERROR_UNIT,*) + ! Release memory from all allocatable arrays call free_memory() +#ifdef MPI + ! Abort MPI + call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err) +#endif + ! Abort program stop diff --git a/src/fileio.F90 b/src/fileio.F90 deleted file mode 100644 index 84aae2c83a..0000000000 --- a/src/fileio.F90 +++ /dev/null @@ -1,90 +0,0 @@ -module fileio - - use constants - use global, only: MAX_LINE_LEN, MAX_WORD_LEN, MAX_WORDS - use string, only: split_string_wl, lower_case - - implicit none - -contains - -!=============================================================================== -! READ_LINE reads a line from a file open on a unit -!=============================================================================== - - subroutine read_line(unit, line, ioError) - - integer, intent(in) :: unit ! unit to read from - character(*), intent(out) :: line ! line to return - integer, intent(out) :: ioError ! error status - - read(UNIT=unit, FMT='(A)', IOSTAT=ioError) line - - end subroutine read_line - -!=============================================================================== -! GET_NEXT_LINE reads the next line to the file connected on the specified unit -! including any continuation lines. If a line ends in an ampersand, the next -! line is read and its words are appended to the final array -!=============================================================================== - - subroutine get_next_line(unit, words, n, ioError) - - integer, intent(in) :: unit ! unit to read from - character(*), intent(out) :: words(MAX_WORDS) ! words read - integer, intent(out) :: n ! number of words - integer, intent(out) :: ioError ! error status - - character(MAX_LINE_LEN) :: line ! single line - character(MAX_WORD_LEN) :: local_words(MAX_WORDS) ! words on one line - integer :: index_word ! index of words - - index_word = 0 - do - ! read line from file - read(UNIT=unit, FMT='(A100)', IOSTAT=ioError) line - - ! if we're at the end of the file, return - if (ioError /= 0) return - - ! split a single line into words - call split_string_wl(line, local_words, n) - - ! if there are no words, we're done - if (n == 0) exit - - ! Check whether there is a continuation line - if (local_words(n) == '&') then - words(index_word+1:index_word+n-1) = local_words(1:n-1) - index_word = index_word + n - 1 - else - words(index_word+1:index_word+n) = local_words(1:n) - index_word = index_word + n - exit - end if - end do - - ! set total number of words - n = index_word - - end subroutine get_next_line - -!=============================================================================== -! SKIP_LINES skips 'n_lines' lines from a file open on a unit -!=============================================================================== - - subroutine skip_lines(unit, n_lines, ioError) - - integer, intent(in) :: unit ! unit to read from - integer, intent(in) :: n_lines ! number of lines to skip - integer, intent(out) :: ioError ! error status - - integer :: i ! index for number of lines - - do i = 1, n_lines - read(UNIT=unit, FMT=*, IOSTAT=ioError) - end do - - end subroutine skip_lines - -end module fileio diff --git a/src/fission.F90 b/src/fission.F90 index c544385361..f82bcf3540 100644 --- a/src/fission.F90 +++ b/src/fission.F90 @@ -1,11 +1,11 @@ module fission + use ace_header, only: Nuclide use constants - use cross_section_header, only: Nuclide - use error, only: fatal_error - use global, only: message - use interpolation, only: interpolate_tab1 - use search, only: binary_search + use error, only: fatal_error + use global, only: message + use interpolation, only: interpolate_tab1 + use search, only: binary_search implicit none diff --git a/src/geometry.F90 b/src/geometry.F90 index 569755d2e1..35612f23cc 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -7,7 +7,7 @@ module geometry use global use output, only: write_message use particle_header, only: Particle, LocalCoord, deallocate_coord - use string, only: int_to_str + use string, only: to_str use tally, only: score_surface_current implicit none @@ -171,26 +171,42 @@ contains x = ceiling((xyz(1) - lat % x0)/lat % width_x) y = ceiling((xyz(2) - lat % y0)/lat % width_y) - ! Create new level of coordinates - p % in_lower_universe = .true. - allocate(p % coord % next) - - ! adjust local position of particle - p % coord % next % xyz(1) = p % coord % xyz(1) - & - (lat%x0 + (x-0.5_8)*lat%width_x) - p % coord % next % xyz(2) = p % coord % xyz(2) - & - (lat%y0 + (y-0.5_8)*lat%width_y) - p % coord % next % xyz(3) = p % coord % xyz(3) - p % coord % next % uvw = p % coord % uvw + ! Check if lattice coordinates are within bounds + if (x < 1 .or. x > lat % n_x .or. & + y < 1 .or. y > lat % n_y) then - ! Move particle to next level - p % coord => p % coord % next + ! This condition should only get hit in rare circumstances where + ! a neutron hits the corner of a lattice. In this case, the + ! neutron may need to be moved diagonally across the lattice. To + ! do so, we remove all lower coordinate levels and then search + ! from universe 0. + + p % coord => p % coord0 + call deallocate_coord(p % coord % next) + + else + + ! Create new level of coordinates + p % in_lower_universe = .true. + allocate(p % coord % next) - ! set particle lattice indices - p % coord % lattice = c % fill - p % coord % lattice_x = x - p % coord % lattice_y = y - p % coord % universe = lat % element(x,y) + ! adjust local position of particle + p % coord % next % xyz(1) = p % coord % xyz(1) - & + (lat%x0 + (x-0.5_8)*lat%width_x) + p % coord % next % xyz(2) = p % coord % xyz(2) - & + (lat%y0 + (y-0.5_8)*lat%width_y) + p % coord % next % xyz(3) = p % coord % xyz(3) + p % coord % next % uvw = p % coord % uvw + + ! Move particle to next level + p % coord => p % coord % next + + ! set particle lattice indices + p % coord % lattice = c % fill + p % coord % lattice_x = x + p % coord % lattice_y = y + p % coord % universe = lat % element(x,y) + end if call find_cell(p, found) if (.not. found) exit @@ -216,24 +232,24 @@ contains type(Particle), pointer :: p integer, intent(in) :: last_cell ! last cell particle was in - real(8) :: x ! x-x0 for sphere - real(8) :: y ! y-y0 for sphere - real(8) :: z ! z-z0 for sphere - real(8) :: R ! radius of sphere - real(8) :: u ! x-component of direction - real(8) :: v ! y-component of direction - real(8) :: w ! z-component of direction - real(8) :: n1 ! x-component of surface normal - real(8) :: n2 ! y-component of surface normal - real(8) :: n3 ! z-component of surface normal - real(8) :: dot_prod ! dot product of direction and normal - real(8) :: norm ! "norm" of surface normal - logical :: found ! particle found in universe? + real(8) :: x ! x-x0 for sphere + real(8) :: y ! y-y0 for sphere + real(8) :: z ! z-z0 for sphere + real(8) :: R ! radius of sphere + real(8) :: u ! x-component of direction + real(8) :: v ! y-component of direction + real(8) :: w ! z-component of direction + real(8) :: n1 ! x-component of surface normal + real(8) :: n2 ! y-component of surface normal + real(8) :: n3 ! z-component of surface normal + real(8) :: dot_prod ! dot product of direction and normal + real(8) :: norm ! "norm" of surface normal + logical :: found ! particle found in universe? type(Surface), pointer :: surf => null() surf => surfaces(abs(p % surface)) if (verbosity >= 10 .or. trace) then - message = " Crossing surface " // trim(int_to_str(surf % id)) + message = " Crossing surface " // trim(to_str(surf % id)) call write_message() end if @@ -258,7 +274,7 @@ contains ! Display message if (verbosity >= 10 .or. trace) then - message = " Leaked out of surface " // trim(int_to_str(surf % id)) + message = " Leaked out of surface " // trim(to_str(surf % id)) call write_message() end if return @@ -269,7 +285,8 @@ contains ! Do not handle reflective boundary conditions on lower universes if (p % in_lower_universe) then - message = "Cannot reflect particle off surface in a lower universe." + message = "Cannot reflect particle " // trim(to_str(p % id)) // & + " off surface in a lower universe." call fatal_error() end if @@ -367,7 +384,7 @@ contains p % coord0 % uvw = (/ u, v, w /) case default message = "Reflection not supported for surface " // & - trim(int_to_str(surf % id)) + trim(to_str(surf % id)) call fatal_error() end select @@ -380,7 +397,7 @@ contains ! Diagnostic message if (verbosity >= 10 .or. trace) then - message = " Reflected from surface " // trim(int_to_str(surf%id)) + message = " Reflected from surface " // trim(to_str(surf%id)) call write_message() end if return @@ -412,9 +429,9 @@ contains ! Couldn't find next cell anywhere! if ((.not. found) .and. (.not. plotting)) then - message = "After particle crossed surface " // trim(int_to_str( & - surfaces(abs(p%surface)) % id)) // " it could not be located in " & - // "any cell and it did not leak." + message = "After particle " // trim(to_str(p % id)) // " crossed surface " & + // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & + &located in any cell and it did not leak." call fatal_error() end if @@ -424,110 +441,69 @@ contains ! CROSS_LATTICE moves a particle into a new lattice element !=============================================================================== - subroutine cross_lattice(p) + subroutine cross_lattice(p, lattice_crossed) type(Particle), pointer :: p + integer, intent(in) :: lattice_crossed - integer :: i_x ! x index in lattice - integer :: i_y ! y index in lattice - real(8) :: d_left ! distance to left side - real(8) :: d_right ! distance to right side - real(8) :: d_bottom ! distance to bottom side - real(8) :: d_top ! distance to top side - real(8) :: dist ! shortest distance - real(8) :: x ! x coordinate in local lattice element - real(8) :: y ! y coordinate in local lattice element - real(8) :: z ! z coordinate in local lattice element - real(8) :: u ! cosine of angle with x axis - real(8) :: v ! cosine of angle with y axis - real(8) :: x0 ! half the width of lattice element - real(8) :: y0 ! half the height of lattice element - logical :: found ! particle found in cell? + integer :: i_x ! x index in lattice + integer :: i_y ! y index in lattice + real(8) :: x0 ! half the width of lattice element + real(8) :: y0 ! half the height of lattice element + logical :: found ! particle found in cell? type(Lattice), pointer :: lat => null() lat => lattices(p % coord % lattice) if (verbosity >= 10 .or. trace) then - message = " Crossing lattice " // trim(int_to_str(lat % id)) // & - ". Current position (" // trim(int_to_str(p % coord % lattice_x)) & - // "," // trim(int_to_str(p % coord % lattice_y)) // ")" + message = " Crossing lattice " // trim(to_str(lat % id)) // & + ". Current position (" // trim(to_str(p % coord % lattice_x)) & + // "," // trim(to_str(p % coord % lattice_y)) // ")" call write_message() end if - u = p % coord % uvw(1) - v = p % coord % uvw(2) - if (lat % type == LATTICE_RECT) then - x = p % coord % xyz(1) - y = p % coord % xyz(2) - z = p % coord % xyz(3) x0 = lat % width_x * 0.5_8 y0 = lat % width_y * 0.5_8 - - dist = INFINITY - ! left and right sides - if (u == ZERO) then - d_left = INFINITY - d_right = INFINITY - elseif (u > 0) then - d_left = INFINITY - d_right = (x0 - x)/u - else - d_left = -(x0 + x)/u - d_right = INFINITY - end if - - ! top and bottom sides - if (v == ZERO) then - d_bottom = INFINITY - d_top = INFINITY - elseif (v > 0) then - d_bottom = INFINITY - d_top = (y0 - y)/v - else - d_bottom = -(y0 + y)/v - d_top = INFINITY - end if - - dist = min(d_left, d_right, d_top, d_bottom) - if (dist == d_left) then + select case (lattice_crossed) + case (LATTICE_LEFT) ! Move particle to left element p % coord % lattice_x = p % coord % lattice_x - 1 p % coord % xyz(1) = x0 - elseif (dist == d_right) then + case (LATTICE_RIGHT) ! Move particle to right element p % coord % lattice_x = p % coord % lattice_x + 1 p % coord % xyz(1) = -x0 - elseif (dist == d_bottom) then + case (LATTICE_BOTTOM) ! Move particle to bottom element p % coord % lattice_y = p % coord % lattice_y - 1 p % coord % xyz(2) = y0 - elseif (dist == d_top) then + case (LATTICE_TOP) ! Move particle to top element p % coord % lattice_y = p % coord % lattice_y + 1 p % coord % xyz(2) = -y0 - end if + end select elseif (lat % type == LATTICE_HEX) then ! TODO: Add hex lattice support end if - + ! Check to make sure still in lattice i_x = p % coord % lattice_x i_y = p % coord % lattice_y if (i_x < 1 .or. i_x > lat % n_x) then - message = "Reached edge of lattice " // trim(int_to_str(lat % id)) // & - " at position (" // trim(int_to_str(i_x)) // "," // & - trim(int_to_str(i_y)) // ")." + message = "Particle " // trim(to_str(p % id)) // " reached edge of & + &lattice " // trim(to_str(lat % id)) // " at position (" // & + trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")." call fatal_error() elseif (i_y < 1 .or. i_y > lat % n_y) then - message = "Reached edge of lattice " // trim(int_to_str(lat % id)) // & - " at position (" // trim(int_to_str(i_x)) // "," // & - trim(int_to_str(i_y)) // ")." + message = "Particle " // trim(to_str(p % id)) // " reached edge of & + &lattice " // trim(to_str(lat % id)) // " at position (" // & + trim(to_str(i_x)) // "," // trim(to_str(i_y)) // ")." call fatal_error() end if @@ -537,8 +513,8 @@ contains ! Find cell in next lattice element call find_cell(p, found) if (.not. found) then - message = "Could not locate particle in universe: " // & - int_to_str(universes(p % coord % universe) % id) + message = "Could not locate particle " // trim(to_str(p % id)) // & + " in universe " // to_str(universes(p % coord % universe) % id) call fatal_error() end if @@ -555,7 +531,7 @@ contains type(Particle), pointer :: p real(8), intent(out) :: dist integer, intent(out) :: surface_crossed - logical, intent(out) :: lattice_crossed + integer, intent(out) :: lattice_crossed integer :: i ! index for surface in cell integer :: index_surf ! index in surfaces array (with sign) @@ -576,7 +552,7 @@ contains ! inialize distance to infinity (huge) dist = INFINITY - lattice_crossed = .false. + lattice_crossed = NONE nullify(final_coord) ! Get pointer to top-level coordinates @@ -866,7 +842,7 @@ contains if (d < dist) then dist = d surface_crossed = -cl % surfaces(i) - lattice_crossed = .false. + lattice_crossed = NONE final_coord => coord end if @@ -906,7 +882,11 @@ contains if (d < dist) then if (abs(d - dist)/dist >= FP_PRECISION) then dist = d - lattice_crossed = .true. + if (u > 0) then + lattice_crossed = LATTICE_RIGHT + else + lattice_crossed = LATTICE_LEFT + end if final_coord => coord end if end if @@ -923,7 +903,11 @@ contains if (d < dist) then if (abs(d - dist)/dist >= FP_PRECISION) then dist = d - lattice_crossed = .true. + if (v > 0) then + lattice_crossed = LATTICE_TOP + else + lattice_crossed = LATTICE_BOTTOM + end if final_coord => coord end if end if @@ -1072,7 +1056,7 @@ contains y1 = surf % coeffs(5) z1 = surf % coeffs(6) if (x >= x0 .and. x < x1 .and. y >= y0 .and. y < y1 .and. & - & z >= z0 .and. z < z1) then + z >= z0 .and. z < z1) then s = SENSE_NEGATIVE else s = SENSE_POSITIVE @@ -1091,7 +1075,7 @@ contains I = surf % coeffs(9) J = surf % coeffs(10) func = A*x*x + B*y*y + C*z*z + D*x*y + E*y*z + F*x*z + G*x & - & + H*y + I*z + J + + H*y + I*z + J end select diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 5b3858755f..8fbf07353e 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -45,8 +45,8 @@ module geometry_header integer :: type ! Type of surface real(8), allocatable :: coeffs(:) ! Definition of surface integer, allocatable :: & - & neighbor_pos(:), & ! List of cells on positive side - & neighbor_neg(:) ! List of cells on negative side + neighbor_pos(:), & ! List of cells on positive side + neighbor_neg(:) ! List of cells on negative side integer :: bc ! Boundary condition end type Surface diff --git a/src/global.F90 b/src/global.F90 index e1e19702d8..5d4a21f44d 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,18 +1,18 @@ module global - use bank_header, only: Bank + use ace_header, only: Nuclide, SAB_Table, xsListing, NuclideMicroXS, & + MaterialMacroXS + use bank_header, only: Bank use cmfd_header use constants - use cross_section_header, only: Nuclide, SAB_Table, xsListing, & - NuclideMicroXS, MaterialMacroXS - use datatypes_header, only: DictionaryII, DictionaryCI - use geometry_header, only: Cell, Universe, Lattice, Surface - use material_header, only: Material - use mesh_header, only: StructuredMesh - use particle_header, only: Particle - use source_header, only: ExtSource - use tally_header, only: TallyObject, TallyMap - use timing, only: Timer + use datatypes_header, only: DictionaryII, DictionaryCI + use geometry_header, only: Cell, Universe, Lattice, Surface + use material_header, only: Material + use mesh_header, only: StructuredMesh + use particle_header, only: Particle + use source_header, only: ExtSource + use tally_header, only: TallyObject, TallyMap + use timing, only: Timer #ifdef MPI use mpi @@ -74,6 +74,9 @@ module global integer :: n_grid ! number of points on unionized grid real(8), allocatable :: e_grid(:) ! energies on unionized grid + ! Unreoslved resonance probablity tables + logical :: urr_ptables_on = .false. + ! ============================================================================ ! TALLY-RELATED VARIABLES @@ -112,6 +115,13 @@ module global real(8) :: keff = ONE real(8) :: keff_std + ! Shannon entropy + logical :: entropy_on = .false. + real(8) :: entropy ! value of shannon entropy + real(8) :: entropy_lower_left(3) ! lower-left corner for entropy box + real(8) :: entropy_upper_right(3) ! upper-right corner for entropy box + real(8), allocatable :: entropy_p(:,:,:) + ! ============================================================================ ! PARALLEL PROCESSING VARIABLES @@ -120,6 +130,7 @@ module global logical :: master ! master process? logical :: mpi_enabled ! is MPI in use and initialized? integer :: mpi_err ! MPI error code + integer :: MPI_BANK ! MPI datatype for fission bank ! ============================================================================ ! TIMING VARIABLES @@ -212,11 +223,6 @@ contains if (allocated(fission_bank)) deallocate(fission_bank) if (allocated(source_bank)) deallocate(source_bank) -#ifdef MPI - ! If MPI is in use and enabled, terminate it - call MPI_FINALIZE(mpi_err) -#endif - end subroutine free_memory end module global diff --git a/src/initialize.F90 b/src/initialize.F90 index cc5314d125..1bcae2ef36 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1,7 +1,8 @@ module initialize + use ace, only: read_xs + use bank_header, only: Bank use constants - use cross_section, only: read_xs use datatypes, only: dict_create, dict_add_key, dict_get_key, & dict_has_key, dict_keys use datatypes_header, only: ListKeyValueII, DictionaryII @@ -12,16 +13,18 @@ module initialize use global use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict - use logging, only: create_log - use mpi_routines, only: setup_mpi use output, only: title, header, print_summary, print_geometry, & - print_plot + print_plot, create_summary_file use random_lcg, only: initialize_prng use source, only: initialize_source - use string, only: int_to_str, starts_with, ends_with, lower_case + use string, only: to_str, starts_with, ends_with, lower_case use tally, only: create_tally_map, TallyObject use timing, only: timer_start, timer_stop +#ifdef MPI + use mpi +#endif + implicit none type(DictionaryII), pointer :: build_dict => null() @@ -47,13 +50,13 @@ contains ! Read command line arguments call read_command_line() - if (master) call create_log() + if (master) call create_summary_file() ! Print the OpenMC title and version/date/time information if (master) call title() ! Print initialization header block - if (master) call header("INITIALIZATION", 1) + if (master) call header("INITIALIZATION", level=1) ! Initialize random number generator call initialize_prng() @@ -120,6 +123,81 @@ contains end subroutine initialize_run +!=============================================================================== +! 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() + +#ifdef MPI + integer :: i + integer :: bank_blocks(4) ! Count for each datatype + integer :: bank_types(4) ! Datatypes + integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements + integer(MPI_ADDRESS_KIND) :: base + type(Bank) :: b + + mpi_enabled = .true. + + ! Initialize MPI + call MPI_INIT(mpi_err) + if (mpi_err /= MPI_SUCCESS) then + message = "Failed to initialize MPI." + call fatal_error() + end if + + ! Determine number of processors + call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err) + if (mpi_err /= MPI_SUCCESS) then + message = "Could not determine number of processors." + call fatal_error() + end if + + ! Determine rank of each processor + call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err) + if (mpi_err /= MPI_SUCCESS) then + message = "Could not determine MPI rank." + call fatal_error() + end if + + ! Determine master + if (rank == 0) then + master = .true. + else + master = .false. + end if + + ! Determine displacements for MPI_BANK type + call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err) + call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err) + call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) + call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) + + ! Adjust displacements + base = bank_disp(1) + do i = 1, 4 + bank_disp(i) = bank_disp(i) - base + end do + + ! Define MPI_BANK for fission sites + bank_blocks = (/ 1, 3, 3, 1 /) + bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) + call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & + bank_types, MPI_BANK, mpi_err) + call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) + +#else + ! if no MPI, set processor to master + mpi_enabled = .false. + rank = 0 + n_procs = 1 + master = .true. +#endif + + end subroutine setup_mpi + !=============================================================================== ! READ_COMMAND_LINE reads all parameters from the command line !=============================================================================== @@ -297,8 +375,8 @@ contains i_array = dict_get_key(surface_dict, abs(id)) c % surfaces(j) = sign(i_array, id) else - message = "Could not find surface " // trim(int_to_str(abs(id))) // & - & " specified on cell " // trim(int_to_str(c % id)) + message = "Could not find surface " // trim(to_str(abs(id))) // & + " specified on cell " // trim(to_str(c % id)) call fatal_error() end if end if @@ -311,8 +389,8 @@ contains if (dict_has_key(universe_dict, id)) then c % universe = dict_get_key(universe_dict, id) else - message = "Could not find universe " // trim(int_to_str(id)) // & - " specified on cell " // trim(int_to_str(c % id)) + message = "Could not find universe " // trim(to_str(id)) // & + " specified on cell " // trim(to_str(c % id)) call fatal_error() end if @@ -325,8 +403,8 @@ contains c % type = CELL_NORMAL c % material = dict_get_key(material_dict, id) else - message = "Could not find material " // trim(int_to_str(id)) // & - " specified on cell " // trim(int_to_str(c % id)) + message = "Could not find material " // trim(to_str(id)) // & + " specified on cell " // trim(to_str(c % id)) call fatal_error() end if else @@ -338,8 +416,8 @@ contains c % type = CELL_LATTICE c % fill = dict_get_key(lattice_dict, id) else - message = "Specified fill " // trim(int_to_str(id)) // " on cell " // & - trim(int_to_str(c % id)) // " is neither a universe nor a lattice." + message = "Specified fill " // trim(to_str(id)) // " on cell " // & + trim(to_str(c % id)) // " is neither a universe nor a lattice." call fatal_error() end if end if @@ -356,8 +434,8 @@ contains if (dict_has_key(universe_dict, id)) then l % element(j,k) = dict_get_key(universe_dict, id) else - message = "Invalid universe number " // trim(int_to_str(id)) & - // " specified on lattice " // trim(int_to_str(l % id)) + message = "Invalid universe number " // trim(to_str(id)) & + // " specified on lattice " // trim(to_str(l % id)) call fatal_error() end if end do @@ -376,8 +454,8 @@ contains if (dict_has_key(cell_dict, id)) then t % cell_bins(j) % scalar = dict_get_key(cell_dict, id) else - message = "Could not find cell " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find cell " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -392,8 +470,8 @@ contains if (dict_has_key(surface_dict, id)) then t % surface_bins(j) % scalar = dict_get_key(surface_dict, id) else - message = "Could not find surface " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find surface " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -408,8 +486,8 @@ contains if (dict_has_key(universe_dict, id)) then t % universe_bins(j) % scalar = dict_get_key(universe_dict, id) else - message = "Could not find universe " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find universe " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -424,8 +502,8 @@ contains if (dict_has_key(material_dict, id)) then t % material_bins(j) % scalar = dict_get_key(material_dict, id) else - message = "Could not find material " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find material " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -440,8 +518,8 @@ contains if (dict_has_key(cell_dict, id)) then t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, id) else - message = "Could not find material " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find material " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end do @@ -455,8 +533,8 @@ contains if (dict_has_key(mesh_dict, id)) then t % mesh = dict_get_key(mesh_dict, id) else - message = "Could not find mesh " // trim(int_to_str(id)) // & - & " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find mesh " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if end if @@ -546,9 +624,9 @@ contains ! Check to make sure either all atom percents or all weight percents are ! given if (.not. (all(mat%atom_percent > ZERO) .or. & - & all(mat%atom_percent < ZERO))) then + all(mat%atom_percent < ZERO))) then message = "Cannot mix atom and weight percents in material " // & - & int_to_str(mat % id) + to_str(mat % id) call fatal_error() end if @@ -606,7 +684,7 @@ contains end do sum_percent = ONE / sum_percent mat % density = -mat % density * N_AVOGADRO & - & / MASS_NEUTRON * sum_percent + / MASS_NEUTRON * sum_percent end if ! Calculate nuclide atom densities and deallocate atom_percent array diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cc8b2e21d2..1f701440bf 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -9,14 +9,14 @@ module input_xml use global use mesh_header, only: StructuredMesh use output, only: write_message - use string, only: lower_case, int_to_str, str_to_int, str_to_real, & + use string, only: lower_case, to_str, str_to_int, str_to_real, & split_string, starts_with, ends_with use tally_header, only: TallyObject implicit none - type(DictionaryII), pointer :: & ! used to count how many cells each - & cells_in_univ_dict => null() ! universe contains + type(DictionaryII), pointer :: & ! used to count how many cells each + cells_in_univ_dict => null() ! universe contains contains @@ -80,9 +80,8 @@ contains ! variable call get_environment_variable("CROSS_SECTIONS", env_variable) if (len_trim(env_variable) == 0) then - message = "No cross_sections.xml file was specified in " // & - "settings.xml or in the CROSS_SECTIONS environment " // & - "variable." + message = "No cross_sections.xml file was specified in settings.xml & + &or in the CROSS_SECTIONS environment variable." call fatal_error() else path_cross_sections = trim(env_variable) @@ -144,6 +143,30 @@ contains trace_particle = trace_(2) end if + ! Entropy box + if (associated(entropy_box_)) then + ! Check to make sure enough values were supplied + if (size(entropy_box_) /= 6) then + message = "Need to supply lower-left and upper-right coordinates & + &for Shannon entropy box." + call fatal_error() + end if + + ! Copy values + entropy_lower_left = entropy_box_(1:3) + entropy_upper_right = entropy_box_(4:6) + + ! Check on values provided + if (.not. all(entropy_upper_right > entropy_lower_left)) then + message = "Upper-right coordinate must be greater than lower-left & + &coordinate for Shannon entropy box." + call fatal_error() + end if + + ! Turn on Shannon entropy calculation + entropy_on = .true. + end if + end subroutine read_settings_xml !=============================================================================== @@ -203,7 +226,7 @@ contains ! Check to make sure that either material or fill was specified if (c % material == 0 .and. c % fill == 0) then message = "Neither material nor fill was specified for cell " // & - trim(int_to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -217,7 +240,7 @@ contains ! Check to make sure that surfaces were specified if (.not. associated(cell_(i) % surfaces)) then message = "No surfaces specified for cell " // & - trim(int_to_str(c % id)) + trim(to_str(c % id)) call fatal_error() end if @@ -313,12 +336,11 @@ contains n = size(surface_(i) % coeffs) if (n < coeffs_reqd) then message = "Not enough coefficients specified for surface: " // & - trim(int_to_str(s % id)) - print *, n, coeffs_reqd + trim(to_str(s % id)) call fatal_error() elseif (n > coeffs_reqd) then message = "Too many coefficients specified for surface: " // & - trim(int_to_str(s % id)) + trim(to_str(s % id)) call fatal_error() else allocate(s % coeffs(n)) @@ -339,7 +361,7 @@ contains s % bc = BC_PERIODIC case default message = "Unknown boundary condition '" // trim(word) // & - "' specified on surface " // trim(int_to_str(s % id)) + "' specified on surface " // trim(to_str(s % id)) call fatal_error() end select @@ -487,14 +509,14 @@ contains m % density = 1.0e-24 * val case default message = "Unkwown units '" // trim(material_(i) % density % units) & - // "' specified on material " // trim(int_to_str(m % id)) + // "' specified on material " // trim(to_str(m % id)) call fatal_error() end select ! Check to ensure material has at least one nuclide if (.not. associated(material_(i) % nuclides)) then message = "No nuclides specified on material " // & - trim(int_to_str(m % id)) + trim(to_str(m % id)) call fatal_error() end if @@ -514,7 +536,7 @@ contains ! Check for empty name on nuclide if (len_trim(nuc % name) == 0) then message = "No name specified on nuclide in material " // & - trim(int_to_str(m % id)) + trim(to_str(m % id)) call fatal_error() end if @@ -522,7 +544,7 @@ contains if (len_trim(nuc % xs) == 0) then if (default_xs == '') then message = "No cross section specified for nuclide in material " & - // trim(int_to_str(m % id)) + // trim(to_str(m % id)) call fatal_error() else nuc % xs = default_xs @@ -540,8 +562,8 @@ contains trim(name) call fatal_error() elseif (nuc % ao /= ZERO .and. nuc % wo /= ZERO) then - message = "Cannot specify both atom and weight percents for a " & - // "nuclide: " // trim(name) + message = "Cannot specify both atom and weight percents for a & + &nuclide: " // trim(name) call fatal_error() end if @@ -671,16 +693,16 @@ contains ! Read mesh origin location if (m % n_dimension /= size(mesh_(i) % origin)) then - message = "Number of entries on must be the same as " // & - "the number of entries on ." + message = "Number of entries on must be the same as the & + &number of entries on ." call fatal_error() end if m % origin = mesh_(i) % origin ! Read mesh widths if (size(mesh_(i) % width) /= size(mesh_(i) % origin)) then - message = "Number of entries on must be the same as " // & - "the number of entries on ." + message = "Number of entries on must be the same as the & + &number of entries on ." call fatal_error() end if m % width = mesh_(i) % width @@ -713,7 +735,7 @@ contains if (len_trim(tally_(i) % filters % cell) > 0 .and. & len_trim(tally_(i) % filters % surface) > 0) then message = "Cannot specify both cell and surface filters for tally " & - // trim(int_to_str(t % id)) + // trim(to_str(t % id)) call fatal_error() end if @@ -768,8 +790,8 @@ contains i_mesh = dict_get_key(mesh_dict, id) m => meshes(i_mesh) else - message = "Could not find mesh " // trim(int_to_str(id)) // & - " specified on tally " // trim(int_to_str(t % id)) + message = "Could not find mesh " // trim(to_str(id)) // & + " specified on tally " // trim(to_str(t % id)) call fatal_error() end if @@ -823,8 +845,8 @@ contains case ('total') t % macro_bins(j) % scalar = MACRO_TOTAL if (t % n_bins(T_ENERGYOUT) > 0) then - message = "Cannot tally total reaction rate with an " & - // "outgoing energy filter." + message = "Cannot tally total reaction rate with an & + &outgoing energy filter." call fatal_error() end if case ('scatter') @@ -848,15 +870,15 @@ contains case ('absorption') t % macro_bins(j) % scalar = MACRO_ABSORPTION if (t % n_bins(T_ENERGYOUT) > 0) then - message = "Cannot tally absorption rate with an outgoing " & - // "energy filter." + message = "Cannot tally absorption rate with an outgoing & + &energy filter." call fatal_error() end if case ('fission') t % macro_bins(j) % scalar = MACRO_FISSION if (t % n_bins(T_ENERGYOUT) > 0) then - message = "Cannot tally fission rate with an outgoing " & - // "energy filter." + message = "Cannot tally fission rate with an outgoing & + &energy filter." call fatal_error() end if case ('nu-fission') @@ -868,9 +890,9 @@ contains ! Check to make sure that current is the only desired response ! for this tally if (n_words > 1) then - message = "Cannot tally other macro reactions in the same " & - // "tally as surface currents. Separate other macro " & - // "reactions into a distinct tally." + message = "Cannot tally other macro reactions in the same & + &tally as surface currents. Separate other macro & + &reactions into a distinct tally." call fatal_error() end if diff --git a/src/mpi_routines.F90 b/src/intercycle.F90 similarity index 61% rename from src/mpi_routines.F90 rename to src/intercycle.F90 index a523694772..50301927ce 100644 --- a/src/mpi_routines.F90 +++ b/src/intercycle.F90 @@ -1,7 +1,8 @@ -module mpi_routines +module intercycle - use constants, only: MAX_LINE_LEN - use error, only: fatal_error + use ISO_FORTRAN_ENV + + use error, only: fatal_error, warning use global use output, only: write_message use particle_header, only: Particle, initialize_particle @@ -13,88 +14,8 @@ module mpi_routines use mpi #endif - 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() - -#ifdef MPI - integer :: i - integer :: bank_blocks(4) ! Count for each datatype - integer :: bank_types(4) ! Datatypes - integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements - integer(MPI_ADDRESS_KIND) :: base - type(Bank) :: b - - mpi_enabled = .true. - - ! Initialize MPI - call MPI_INIT(mpi_err) - if (mpi_err /= MPI_SUCCESS) then - message = "Failed to initialize MPI." - call fatal_error() - end if - - ! Determine number of processors - call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err) - if (mpi_err /= MPI_SUCCESS) then - message = "Could not determine number of processors." - call fatal_error() - end if - - ! Determine rank of each processor - call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err) - if (mpi_err /= MPI_SUCCESS) then - message = "Could not determine MPI rank." - call fatal_error() - end if - - ! Determine master - if (rank == 0) then - master = .true. - else - master = .false. - end if - - ! Determine displacements for MPI_BANK type - call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err) - call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err) - call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err) - call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err) - - ! Adjust displacements - base = bank_disp(1) - do i = 1, 4 - bank_disp(i) = bank_disp(i) - base - end do - - ! Define MPI_BANK for fission sites - bank_blocks = (/ 1, 3, 3, 1 /) - bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) - call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, & - & bank_types, MPI_BANK, mpi_err) - call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) - -#else - ! if no MPI, set processor to master - mpi_enabled = .false. - rank = 0 - n_procs = 1 - master = .true. -#endif - - end subroutine setup_mpi - !=============================================================================== ! SYNCHRONIZE_BANK samples source sites from the fission sites that were ! accumulated during the cycle. This routine is what allows this Monte Carlo to @@ -115,9 +36,9 @@ contains integer(8) :: sites_needed ! # of sites to be sampled real(8) :: p_sample ! probability of sampling a site type(Bank), allocatable :: & - & temp_sites(:), & ! local array of extra sites on each node - & left_bank(:), & ! bank sites to send/recv to or from left node - & right_bank(:) ! bank sites to send/recv to or fram right node + temp_sites(:), & ! local array of extra sites on each node + left_bank(:), & ! bank sites to send/recv to or from left node + right_bank(:) ! bank sites to send/recv to or fram right node #ifdef MPI integer :: status(MPI_STATUS_SIZE) ! message status @@ -133,11 +54,11 @@ contains ! Determine starting index for fission bank and total sites in fission bank start = 0_8 call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) finish = start + n_bank total = finish call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) #else start = 0_8 @@ -200,11 +121,11 @@ contains #ifdef MPI start = 0_8 call MPI_EXSCAN(index_local, start, 1, MPI_INTEGER8, MPI_SUM, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) finish = start + index_local total = finish call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - & MPI_COMM_WORLD, mpi_err) + MPI_COMM_WORLD, mpi_err) #else start = 0_8 finish = index_local @@ -259,18 +180,18 @@ contains if (send_to_right > 0) then i = index_local - send_to_right + 1 call MPI_ISEND(temp_sites(i), send_to_right, MPI_BANK, rank+1, 0, & - & MPI_COMM_WORLD, request, mpi_err) + MPI_COMM_WORLD, request, mpi_err) else if (send_to_right < 0) then call MPI_IRECV(right_bank, -send_to_right, MPI_BANK, rank+1, 1, & - & MPI_COMM_WORLD, request_right, mpi_err) + MPI_COMM_WORLD, request_right, mpi_err) end if if (send_to_left < 0) then call MPI_IRECV(left_bank, -send_to_left, MPI_BANK, rank-1, 0, & - & MPI_COMM_WORLD, request_left, mpi_err) + MPI_COMM_WORLD, request_left, mpi_err) else if (send_to_left > 0) then call MPI_ISEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, & - & MPI_COMM_WORLD, request, mpi_err) + MPI_COMM_WORLD, request, mpi_err) end if #endif @@ -386,7 +307,8 @@ contains tally_temp = t % scores(:,:) % val_history if (master) then - ! Description of MPI_IN_PLANE + ! The MPI_IN_PLACE specifier allows the master to copy values into a + ! receive buffer without having a temporary variable call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & 0, MPI_COMM_WORLD, mpi_err) @@ -408,4 +330,181 @@ contains end subroutine reduce_tallies #endif -end module mpi_routines +!=============================================================================== +! SHANNON_ENTROPY calculates the Shannon entropy of the fission source +! distribution to assess source convergence +!=============================================================================== + + subroutine shannon_entropy() + + integer :: i ! x-index for entropy mesh + integer :: j ! y-index for entropy mesh + integer :: k ! z-index for entropy mesh + integer :: m ! index for bank sites + integer(8) :: total_bank ! total # of fission bank sites + integer, save :: n_box ! total # of boxes on mesh + integer, save :: n ! # of boxes in each dimension + real(8), save :: width(3) ! width of box in each dimension + logical :: outside_box ! were there sites outside entropy box? + + ! On the first pass through this subroutine, we need to determine how big + ! the entropy mesh should be in each direction and then allocate a + ! three-dimensional array to store the fraction of source sites in each mesh + ! box + + if (.not. allocated(entropy_p)) then + ! determine number of boxes in each direction + n = ceiling((n_particles/20)**(1.0/3.0)) + n_box = n*n*n + + ! determine width + width = (entropy_upper_right - entropy_lower_left)/n + + ! allocate p + allocate(entropy_p(n,n,n)) + end if + + ! initialize p + entropy_p = ZERO + outside_box = .false. + + ! loop over fission sites and count how many are in each mesh box + FISSION_SITES: do m = 1, int(n_bank,4) + ! determine indices for entropy mesh box + i = int((fission_bank(m) % xyz(1) - entropy_lower_left(1))/width(1)) + 1 + j = int((fission_bank(m) % xyz(2) - entropy_lower_left(2))/width(2)) + 1 + k = int((fission_bank(m) % xyz(3) - entropy_lower_left(3))/width(3)) + 1 + + ! if outside mesh, skip particle + if (i < 1 .or. i > n .or. j < 1 .or. & + j > n .or. k < 1 .or. k > n) then + outside_box = .true. + cycle + end if + + ! add to appropriate mesh box + entropy_p(i,j,k) = entropy_p(i,j,k) + 1 + end do FISSION_SITES + + ! display warning message if there were sites outside entropy box + if (outside_box) then + message = "Fission source site(s) outside of entropy box." + call warning() + end if + +#ifdef MPI + ! collect values from all processors + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, entropy_p, n_box, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + else + call MPI_REDUCE(entropy_p, entropy_p, n_box, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + end if + + ! determine total number of bank sites + call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) +#else + total_bank = n_bank +#endif + + ! sum values to obtain shannon entropy + if (master) then + entropy_p = entropy_p / total_bank + entropy = -sum(entropy_p * log(entropy_p)/log(2.0), entropy_p > ZERO) + end if + + end subroutine shannon_entropy + +!=============================================================================== +! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the +! mean and standard deviation of the mean for active cycles and displays them +!=============================================================================== + + subroutine calculate_keff(i_cycle) + + integer, intent(in) :: i_cycle ! index of current cycle + + integer(8) :: total_bank ! total number of source sites + integer :: n ! active cycle number + real(8) :: k_cycle ! single cycle estimate of keff + real(8), save :: k_sum ! accumulated keff + real(8), save :: k_sum_sq ! accumulated keff**2 + + message = "Calculate cycle keff..." + call write_message(8) + + ! initialize sum and square of sum at beginning of run + if (i_cycle == 1) then + k_sum = ZERO + k_sum_sq = ZERO + end if + +#ifdef MPI + ! Collect number bank sites onto master process + call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) +#else + total_bank = n_bank +#endif + + ! Collect statistics and print output + if (master) then + ! Since the creation of bank sites was originally weighted by the last + ! cycle keff, we need to multiply by that keff to get the current cycle's + ! value + + k_cycle = real(total_bank)/real(n_particles)*keff + + if (i_cycle > n_inactive) then + ! Active cycle number + n = i_cycle - n_inactive + + ! Accumulate cycle estimate of k + k_sum = k_sum + k_cycle + k_sum_sq = k_sum_sq + k_cycle*k_cycle + + ! Determine mean and standard deviation of mean + keff = k_sum/n + keff_std = sqrt((k_sum_sq/n - keff*keff)/n) + + ! Display output for this cycle + if (i_cycle > n_inactive+1) then + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=103) i_cycle, k_cycle, entropy, & + keff, keff_std + else + write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std + end if + else + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy + else + write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle + end if + end if + else + ! Display output for inactive cycle + if (entropy_on) then + write(UNIT=OUTPUT_UNIT, FMT=102) i_cycle, k_cycle, entropy + else + write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle + end if + keff = k_cycle + end if + end if + +#ifdef MPI + ! Broadcast new keff value to all processors + call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) +#endif + +100 format (2X,I5,2X,F8.5) +101 format (2X,I5,2X,F8.5,5X,F8.5," +/-",F8.5) +102 format (2X,I5,2X,F8.5,3X,F8.5) +103 format (2X,I5,2X,F8.5,3X,F8.5,3X,F8.5," +/-",F8.5) + + end subroutine calculate_keff + +end module intercycle diff --git a/src/interpolation.F90 b/src/interpolation.F90 index 8e3f083ff5..549b720e0b 100644 --- a/src/interpolation.F90 +++ b/src/interpolation.F90 @@ -5,7 +5,7 @@ module interpolation use error, only: fatal_error use global, only: message use search, only: binary_search - use string, only: int_to_str + use string, only: to_str implicit none @@ -118,7 +118,7 @@ contains r = (log(x) - log(x0))/(log(x1) - log(x0)) y = exp((1-r)*log(y0) + r*log(y1)) case default - message = "Unsupported interpolation scheme: " // int_to_str(interp) + message = "Unsupported interpolation scheme: " // to_str(interp) call fatal_error() end select @@ -204,7 +204,7 @@ contains r = (log(x) - log(x0))/(log(x1) - log(x0)) y = exp((1-r)*log(y0) + r*log(y1)) case default - message = "Unsupported interpolation scheme: " // int_to_str(interp) + message = "Unsupported interpolation scheme: " // to_str(interp) call fatal_error() end select diff --git a/src/logging.F90 b/src/logging.F90 deleted file mode 100644 index 377c1ab827..0000000000 --- a/src/logging.F90 +++ /dev/null @@ -1,43 +0,0 @@ -module logging - - use constants, only: MAX_WORD_LEN, UNIT_LOG - use global, only: path_input - - implicit none - -contains - -!=============================================================================== -! CREATE_LOG creates a new log file (or overwrites the existing log) -!=============================================================================== - - subroutine create_log() - - character(MAX_WORD_LEN) :: path_log ! path of log file - logical :: file_exists ! does log file already exist? - ! integer :: ioError ! error status for file access - - ! Create filename for log file - path_log = trim(path_input) // ".log" - - ! Check if log file already exists - inquire(FILE=path_log, EXIST=file_exists) - if (file_exists) then - ! Possibly copy old log file - end if - - ! Open log file for writing - ! open(FILE=path_log, UNIT=UNIT_LOG, STATUS='replace', & - ! & ACTION='write', IOSTAT=ioError) - - end subroutine create_log - -!=============================================================================== -! LOG_TALLIES -!=============================================================================== - - subroutine log_tallies() - - end subroutine log_tallies - -end module logging diff --git a/src/main.F90 b/src/main.F90 index a08e24eff0..2631f71b15 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -4,16 +4,16 @@ program main use constants use global use initialize, only: initialize_run - use mpi_routines, only: synchronize_bank + use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank use output, only: write_message, header, print_runtime use particle_header, only: Particle use plot, only: run_plot use physics, only: transport use random_lcg, only: set_particle_seed use source, only: get_source_particle - use string, only: int_to_str + use string, only: to_str use tally, only: synchronize_tallies, write_tallies, & - tally_statistics, calculate_keff + tally_statistics use timing, only: timer_start, timer_stop #ifdef MPI @@ -48,6 +48,11 @@ program main ! deallocate arrays call free_memory() + +#ifdef MPI + ! If MPI is in use and enabled, terminate it + call MPI_FINALIZE(mpi_err) +#endif contains @@ -62,11 +67,24 @@ contains integer(8) :: i_particle ! history index type(Particle), pointer :: p => null() - if (master) call header("BEGIN SIMULATION", 1) + if (master) call header("BEGIN SIMULATION", level=1) tallies_on = .false. call timer_start(time_inactive) + ! Display column titles + if (entropy_on) then + message = " Cycle k(cycle) Entropy Average k" + call write_message(1) + message = " ===== ======== ======= ===================" + call write_message(1) + else + message = " Cycle k(cycle) Average k" + call write_message(1) + message = " ===== ======== ===================" + call write_message(1) + end if + ! ========================================================================== ! LOOP OVER CYCLES CYCLE_LOOP: do i_cycle = 1, n_cycles @@ -74,9 +92,9 @@ contains ! Start timer for computation call timer_start(time_compute) - message = "Simulating cycle " // trim(int_to_str(i_cycle)) // "..." + message = "Simulating cycle " // trim(to_str(i_cycle)) // "..." call write_message(8) - + ! Set all tallies to zero n_bank = 0 @@ -121,6 +139,9 @@ contains call timer_stop(time_ic_tallies) end if + ! Calculate shannon entropy + if (entropy_on) call shannon_entropy() + ! Distribute fission bank across processors evenly call synchronize_bank(i_cycle) @@ -146,7 +167,7 @@ contains ! ========================================================================== ! END OF RUN WRAPUP - if (master) call header("SIMULATION FINISHED", 1) + if (master) call header("SIMULATION FINISHED", level=1) end subroutine run_problem diff --git a/src/output.F90 b/src/output.F90 index f0f7d6c9a5..7b2898528c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -2,16 +2,16 @@ module output use ISO_FORTRAN_ENV + use ace_header, only: Nuclide, Reaction, UrrData use constants - use cross_section_header, only: Reaction - use datatypes, only: dict_get_key - use endf, only: reaction_name - use geometry_header, only: Cell, Universe, Surface + use datatypes, only: dict_get_key + use endf, only: reaction_name + use geometry_header, only: Cell, Universe, Surface use global - use mesh_header, only: StructuredMesh - use particle_header, only: Particle, LocalCoord - use string, only: upper_case, int_to_str, real_to_str - use tally_header, only: TallyObject + use mesh_header, only: StructuredMesh + use particle_header, only: Particle, LocalCoord + use string, only: upper_case, to_str + use tally_header, only: TallyObject implicit none @@ -31,31 +31,46 @@ contains character(10) :: today_date character(8) :: today_time - write(ou,*) - write(ou,*) ' .d88888b. 888b d888 .d8888b.' - write(ou,*) ' d88P" "Y88b 8888b d8888 d88P Y88b' - write(ou,*) ' 888 888 88888b.d88888 888 888' - write(ou,*) ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ' - write(ou,*) ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ' - write(ou,*) ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888' - write(ou,*) ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P' - write(ou,*) ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"' - write(ou,*) '__________________888______________________________________________________' - write(ou,*) ' 888' - write(ou,*) ' 888' - write(ou,*) + write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & + ' .d88888b. 888b d888 .d8888b.', & + ' d88P" "Y88b 8888b d8888 d88P Y88b', & + ' 888 888 88888b.d88888 888 888', & + ' 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 ', & + ' 888 888 888 "88b d8P Y8b 888 "88b 888 Y888P 888 888 ', & + ' 888 888 888 888 88888888 888 888 888 Y8P 888 888 888', & + ' Y88b. .d88P 888 d88P Y8b. 888 888 888 " 888 Y88b d88P', & + ' "Y88888P" 88888P" "Y8888 888 888 888 888 "Y8888P"', & + '__________________888______________________________________________________', & + ' 888', & + ' 888' ! Write version information - write(ou,*) ' Developed At: Massachusetts Institute of Technology' - write(ou,*) ' Lead Developer: Paul K. Romano' - write(ou,100) VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE -100 format (6X,"Version:",9X,I1,".",I1,".",I1) + write(UNIT=OUTPUT_UNIT, FMT=*) & + ' Developed At: Massachusetts Institute of Technology' + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Version:",7X,I1,".",I1,".",I1)') & + VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE ! Write the date and time call get_today(today_date, today_time) - write(ou,101) trim(today_date), trim(today_time) -101 format (6X,"Date/Time:",7X,A,1X,A) - write(ou,*) + write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') & + trim(today_date), trim(today_time) + + ! Write information to summary file + call header("OpenMC Monte Carlo Code", unit=UNIT_SUMMARY, level=1) + write(UNIT=UNIT_SUMMARY, FMT=*) & + "Copyright: 2011 Massachusetts Institute of Technology" + write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & + "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') & + trim(today_date), trim(today_time) + + ! Write information on number of processors +#ifdef MPI + write(UNIT=OUTPUT_UNIT, FMT='(1X,A)') ' MPI Processes: ' // & + trim(to_str(n_procs)) + write(UNIT=UNIT_SUMMARY, FMT='(1X,"MPI Processes:",1X,A)') & + trim(to_str(n_procs)) +#endif end subroutine title @@ -64,34 +79,32 @@ contains ! specified, it is assumed to be a minor header block (H3). !=============================================================================== - subroutine header(msg, level, unit_file) + subroutine header(msg, unit, level) character(*), intent(in) :: msg - integer, optional :: level - integer, optional :: unit_file + integer, optional :: unit + integer, optional :: level + integer :: n + integer :: m + integer :: unit_ integer :: header_level - integer :: n, m - integer :: unit - character(75) :: line + character(MAX_LINE_LEN) :: line - ! set default header level - if (.not. present(level)) then - header_level = 3 - else + ! set default level + if (present(level)) then header_level = level + else + header_level = 3 end if ! set default unit - if (present(unit_file)) then - unit = unit_file + if (present(unit)) then + unit_ = unit else - unit = ou + unit_ = OUTPUT_UNIT end if - ! Print first blank line - write(unit,*) - ! determine how many times to repeat '=' character n = (63 - len_trim(msg))/2 m = n @@ -104,23 +117,16 @@ contains ! print header based on level select case (header_level) case (1) - ! determine number of spaces to put in from of header - write(unit,*) repeat('=', 75) - write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // & - & repeat('=', m) - write(unit,*) repeat('=', 75) + write(UNIT=unit_, FMT='(/3(1X,A/))') repeat('=', 75), & + repeat('=', n) // '> ' // trim(line) // ' <' // & + repeat('=', m), repeat('=', 75) case (2) - write(unit,*) trim(line) - write(unit,*) repeat('-', 75) + write(UNIT=unit_, FMT='(/2(1X,A/))') trim(line), repeat('-', 75) case (3) - n = (63 - len_trim(line))/2 - write(unit,*) repeat('=', n) // '> ' // trim(line) // ' <' // & - & repeat('=', m) + write(UNIT=unit_, FMT='(/1X,A/)') repeat('=', n) // '> ' // & + trim(line) // ' <' // repeat('=', m) end select - ! Print trailing blank line - write(unit, *) - end subroutine header !=============================================================================== @@ -207,13 +213,13 @@ contains ! display type of particle select case (p % type) case (NEUTRON) - write(ou,*) 'Neutron ' // int_to_str(p % id) + write(ou,*) 'Neutron ' // to_str(p % id) case (PHOTON) - write(ou,*) 'Photon ' // int_to_str(p % id) + write(ou,*) 'Photon ' // to_str(p % id) case (ELECTRON) - write(ou,*) 'Electron ' // int_to_str(p % id) + write(ou,*) 'Electron ' // to_str(p % id) case default - write(ou,*) 'Unknown Particle ' // int_to_str(p % id) + write(ou,*) 'Unknown Particle ' // to_str(p % id) end select ! loop through each level of universes @@ -221,26 +227,26 @@ contains i = 0 do while(associated(coord)) ! Print level - write(ou,*) ' Level ' // trim(int_to_str(i)) + write(ou,*) ' Level ' // trim(to_str(i)) ! Print cell for this level if (coord % cell /= NONE) then c => cells(coord % cell) - write(ou,*) ' Cell = ' // trim(int_to_str(c % id)) + write(ou,*) ' Cell = ' // trim(to_str(c % id)) end if ! Print universe for this level if (coord % universe /= NONE) then u => universes(coord % universe) - write(ou,*) ' Universe = ' // trim(int_to_str(u % id)) + write(ou,*) ' Universe = ' // trim(to_str(u % id)) end if ! Print information on lattice if (coord % lattice /= NONE) then l => lattices(coord % lattice) - write(ou,*) ' Lattice = ' // trim(int_to_str(l % id)) - write(ou,*) ' Lattice position = (' // trim(int_to_str(& - p % coord % lattice_x)) // ',' // trim(int_to_str(& + write(ou,*) ' Lattice = ' // trim(to_str(l % id)) + write(ou,*) ' Lattice position = (' // trim(to_str(& + p % coord % lattice_x)) // ',' // trim(to_str(& p % coord % lattice_y)) // ')' end if @@ -255,13 +261,13 @@ contains ! Print surface if (p % surface /= NONE) then s => surfaces(p % surface) - write(ou,*) ' Surface = ' // int_to_str(s % id) + write(ou,*) ' Surface = ' // to_str(s % id) end if - write(ou,*) ' Weight = ' // real_to_str(p % wgt) - write(ou,*) ' Energy = ' // real_to_str(p % E) - write(ou,*) ' IE = ' // int_to_str(p % IE) - write(ou,*) ' Interpolation factor = ' // real_to_str(p % interp) + write(ou,*) ' Weight = ' // to_str(p % wgt) + write(ou,*) ' Energy = ' // to_str(p % E) + write(ou,*) ' IE = ' // to_str(p % IE) + write(ou,*) ' Interpolation factor = ' // to_str(p % interp) write(ou,*) end subroutine print_particle @@ -275,12 +281,12 @@ contains type(Reaction), pointer :: rxn write(ou,*) 'Reaction ' // reaction_name(rxn % MT) - write(ou,*) ' MT = ' // int_to_str(rxn % MT) - write(ou,*) ' Q-value = ' // real_to_str(rxn % Q_value) - write(ou,*) ' TY = ' // int_to_str(rxn % TY) - write(ou,*) ' Starting index = ' // int_to_str(rxn % IE) + write(ou,*) ' MT = ' // to_str(rxn % MT) + write(ou,*) ' Q-value = ' // to_str(rxn % Q_value) + write(ou,*) ' TY = ' // to_str(rxn % TY) + write(ou,*) ' Starting index = ' // to_str(rxn % IE) if (rxn % has_energy_dist) then - write(ou,*) ' Energy: Law ' // int_to_str(rxn % edist % law) + write(ou,*) ' Energy: Law ' // to_str(rxn % edist % law) end if write(ou,*) @@ -290,39 +296,47 @@ contains ! PRINT_CELL displays the attributes of a cell !=============================================================================== - subroutine print_cell(c) + subroutine print_cell(c, unit) type(Cell), pointer :: c + integer, optional :: unit - integer :: temp - integer :: i + integer :: temp + integer :: i + integer :: unit_ character(MAX_LINE_LEN) :: string type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() type(Material), pointer :: m => null() - write(ou,*) 'Cell ' // int_to_str(c % id) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Cell ' // to_str(c % id) temp = dict_get_key(cell_dict, c % id) - write(ou,*) ' Array Index = ' // int_to_str(temp) + write(unit_,*) ' Array Index = ' // to_str(temp) u => universes(c % universe) - write(ou,*) ' Universe = ' // int_to_str(u % id) + write(unit_,*) ' Universe = ' // to_str(u % id) select case (c % type) case (CELL_NORMAL) - write(ou,*) ' Fill = NONE' + write(unit_,*) ' Fill = NONE' case (CELL_FILL) u => universes(c % fill) - write(ou,*) ' Fill = Universe ' // int_to_str(u % id) + write(unit_,*) ' Fill = Universe ' // to_str(u % id) case (CELL_LATTICE) l => lattices(c % fill) - write(ou,*) ' Fill = Lattice ' // int_to_str(l % id) + write(unit_,*) ' Fill = Lattice ' // to_str(l % id) end select if (c % material == 0) then - write(ou,*) ' Material = NONE' + write(unit_,*) ' Material = NONE' else m => materials(c % material) - write(ou,*) ' Material = ' // int_to_str(m % id) + write(unit_,*) ' Material = ' // to_str(m % id) end if - write(ou,*) ' Parent Cell = ' // int_to_str(c % parent) + write(unit_,*) ' Parent Cell = ' // to_str(c % parent) string = "" do i = 1, c % n_surfaces select case (c % surfaces(i)) @@ -335,11 +349,11 @@ contains case (OP_DIFFERENCE) string = trim(string) // ' !' case default - string = trim(string) // ' ' // int_to_str(c % surfaces(i)) + string = trim(string) // ' ' // to_str(c % surfaces(i)) end select end do - write(ou,*) ' Surface Specification:' // trim(string) - write(ou,*) + write(unit_,*) ' Surface Specification:' // trim(string) + write(unit_,*) end subroutine print_cell @@ -347,23 +361,31 @@ contains ! PRINT_UNIVERSE displays the attributes of a universe !=============================================================================== - subroutine print_universe(univ) + subroutine print_universe(univ, unit) type(Universe), pointer :: univ + integer, optional :: unit - integer :: i + integer :: i + integer :: unit_ character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() - write(ou,*) 'Universe ' // int_to_str(univ % id) - write(ou,*) ' Level = ' // int_to_str(univ % level) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Universe ' // to_str(univ % id) + write(unit_,*) ' Level = ' // to_str(univ % level) string = "" do i = 1, univ % n_cells c => cells(univ % cells(i)) - string = trim(string) // ' ' // int_to_str(c % id) + string = trim(string) // ' ' // to_str(c % id) end do - write(ou,*) ' Cells =' // trim(string) - write(ou,*) + write(unit_,*) ' Cells =' // trim(string) + write(unit_,*) end subroutine print_universe @@ -371,18 +393,27 @@ contains ! PRINT_LATTICE displays the attributes of a lattice !=============================================================================== - subroutine print_lattice(lat) + subroutine print_lattice(lat, unit) type(Lattice), pointer :: lat + integer, optional :: unit - write(ou,*) 'Lattice ' // int_to_str(lat % id) - write(ou,*) ' n_x = ' // int_to_str(lat % n_x) - write(ou,*) ' n_y = ' // int_to_str(lat % n_y) - write(ou,*) ' x0 = ' // real_to_str(lat % x0) - write(ou,*) ' y0 = ' // real_to_str(lat % y0) - write(ou,*) ' width_x = ' // real_to_str(lat % width_x) - write(ou,*) ' width_y = ' // real_to_str(lat % width_y) - write(ou,*) + integer :: unit_ + + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Lattice ' // to_str(lat % id) + write(unit_,*) ' n_x = ' // to_str(lat % n_x) + write(unit_,*) ' n_y = ' // to_str(lat % n_y) + write(unit_,*) ' x0 = ' // to_str(lat % x0) + write(unit_,*) ' y0 = ' // to_str(lat % y0) + write(unit_,*) ' width_x = ' // to_str(lat % width_x) + write(unit_,*) ' width_y = ' // to_str(lat % width_y) + write(unit_,*) end subroutine print_lattice @@ -390,14 +421,22 @@ contains ! PRINT_SURFACE displays the attributes of a surface !=============================================================================== - subroutine print_surface(surf) + subroutine print_surface(surf, unit) type(Surface), pointer :: surf + integer, optional :: unit integer :: i + integer :: unit_ character(MAX_LINE_LEN) :: string - write(ou,*) 'Surface ' // int_to_str(surf % id) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Surface ' // to_str(surf % id) select case (surf % type) case (SURF_PX) string = "X Plane" @@ -422,40 +461,40 @@ contains case (SURF_GQ) string = "General Quadratic" end select - write(ou,*) ' Type = ' // trim(string) + write(unit_,*) ' Type = ' // trim(string) string = "" do i = 1, size(surf % coeffs) - string = trim(string) // ' ' // real_to_str(surf % coeffs(i), 4) + string = trim(string) // ' ' // to_str(surf % coeffs(i), 4) end do - write(ou,*) ' Coefficients = ' // trim(string) + write(unit_,*) ' Coefficients = ' // trim(string) string = "" if (allocated(surf % neighbor_pos)) then do i = 1, size(surf % neighbor_pos) - string = trim(string) // ' ' // int_to_str(surf % neighbor_pos(i)) + string = trim(string) // ' ' // to_str(surf % neighbor_pos(i)) end do end if - write(ou,*) ' Positive Neighbors = ' // trim(string) + write(unit_,*) ' Positive Neighbors = ' // trim(string) string = "" if (allocated(surf % neighbor_neg)) then do i = 1, size(surf % neighbor_neg) - string = trim(string) // ' ' // int_to_str(surf % neighbor_neg(i)) + string = trim(string) // ' ' // to_str(surf % neighbor_neg(i)) end do end if - write(ou,*) ' Negative Neighbors =' // trim(string) + write(unit_,*) ' Negative Neighbors =' // trim(string) select case (surf % bc) case (BC_TRANSMIT) - write(ou,*) ' Boundary Condition = Transmission' + write(unit_,*) ' Boundary Condition = Transmission' case (BC_VACUUM) - write(ou,*) ' Boundary Condition = Vacuum' + write(unit_,*) ' Boundary Condition = Vacuum' case (BC_REFLECT) - write(ou,*) ' Boundary Condition = Reflective' + write(unit_,*) ' Boundary Condition = Reflective' case (BC_PERIODIC) - write(ou,*) ' Boundary Condition = Periodic' + write(unit_,*) ' Boundary Condition = Periodic' end select - write(ou,*) + write(unit_,*) end subroutine print_surface @@ -463,37 +502,45 @@ contains ! PRINT_MATERIAL displays the attributes of a material !=============================================================================== - subroutine print_material(mat) + subroutine print_material(mat, unit) type(Material), pointer :: mat + integer, optional :: unit - integer :: i - real(8) :: density + integer :: i + integer :: unit_ + real(8) :: density character(MAX_LINE_LEN) :: string type(Nuclide), pointer :: nuc => null() + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + ! Write identifier for material - write(ou,*) 'Material ' // int_to_str(mat % id) + write(unit_,*) 'Material ' // to_str(mat % id) ! Write total atom density in atom/b-cm - write(ou,*) ' Atom Density = ' // trim(real_to_str(mat % density)) & - & // ' atom/b-cm' + write(unit_,*) ' Atom Density = ' // trim(to_str(mat % density)) & + // ' atom/b-cm' ! Write atom density for each nuclide in material - write(ou,*) ' Nuclides:' + write(unit_,*) ' Nuclides:' do i = 1, mat % n_nuclides nuc => nuclides(mat % nuclide(i)) density = mat % atom_density(i) string = ' ' // trim(nuc % name) // ' = ' // & - & trim(real_to_str(density)) // ' atom/b-cm' - write(ou,*) trim(string) + trim(to_str(density)) // ' atom/b-cm' + write(unit_,*) trim(string) end do ! Write information on S(a,b) table if (mat % has_sab_table) then - write(ou,*) ' S(a,b) table = ' // trim(mat % sab_name) + write(unit_,*) ' S(a,b) table = ' // trim(mat % sab_name) end if - write(ou,*) + write(unit_,*) end subroutine print_material @@ -501,29 +548,37 @@ contains ! PRINT_TALLY displays the attributes of a tally !=============================================================================== - subroutine print_tally(t) + subroutine print_tally(t, unit) type(TallyObject), pointer :: t + integer, optional :: unit - integer :: i - integer :: id - character(MAX_LINE_LEN) :: string + integer :: i + integer :: id + integer :: unit_ + character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() type(Material), pointer :: m => null() type(StructuredMesh), pointer :: sm => null() - write(ou,*) 'Tally ' // int_to_str(t % id) + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + write(unit_,*) 'Tally ' // to_str(t % id) if (t % n_bins(T_CELL) > 0) then string = "" do i = 1, t % n_bins(T_CELL) id = t % cell_bins(i) % scalar c => cells(id) - string = trim(string) // ' ' // trim(int_to_str(c % id)) + string = trim(string) // ' ' // trim(to_str(c % id)) end do - write(ou, *) ' Cell Bins:' // trim(string) + write(unit_, *) ' Cell Bins:' // trim(string) end if if (t % n_bins(T_SURFACE) > 0) then @@ -531,9 +586,9 @@ contains do i = 1, t % n_bins(T_SURFACE) id = t % surface_bins(i) % scalar s => surfaces(id) - string = trim(string) // ' ' // trim(int_to_str(s % id)) + string = trim(string) // ' ' // trim(to_str(s % id)) end do - write(ou, *) ' Surface Bins:' // trim(string) + write(unit_, *) ' Surface Bins:' // trim(string) end if if (t % n_bins(T_UNIVERSE) > 0) then @@ -541,9 +596,9 @@ contains do i = 1, t % n_bins(T_UNIVERSE) id = t % universe_bins(i) % scalar u => universes(id) - string = trim(string) // ' ' // trim(int_to_str(u % id)) + string = trim(string) // ' ' // trim(to_str(u % id)) end do - write(ou, *) ' Material Bins:' // trim(string) + write(unit_, *) ' Material Bins:' // trim(string) end if if (t % n_bins(T_MATERIAL) > 0) then @@ -551,20 +606,20 @@ contains do i = 1, t % n_bins(T_MATERIAL) id = t % material_bins(i) % scalar m => materials(id) - string = trim(string) // ' ' // trim(int_to_str(m % id)) + string = trim(string) // ' ' // trim(to_str(m % id)) end do - write(ou, *) ' Material Bins:' // trim(string) + write(unit_, *) ' Material Bins:' // trim(string) end if if (t % n_bins(T_MESH) > 0) then string = "" id = t % mesh sm => meshes(id) - string = trim(string) // ' ' // trim(int_to_str(sm % dimension(1))) + string = trim(string) // ' ' // trim(to_str(sm % dimension(1))) do i = 2, sm % n_dimension - string = trim(string) // ' x ' // trim(int_to_str(sm % dimension(i))) + string = trim(string) // ' x ' // trim(to_str(sm % dimension(i))) end do - write(ou, *) ' Mesh Bins:' // trim(string) + write(unit_, *) ' Mesh Bins:' // trim(string) end if if (t % n_bins(T_CELLBORN) > 0) then @@ -572,27 +627,27 @@ contains do i = 1, t % n_bins(T_CELLBORN) id = t % cellborn_bins(i) % scalar c => cells(id) - string = trim(string) // ' ' // trim(int_to_str(c % id)) + string = trim(string) // ' ' // trim(to_str(c % id)) end do - write(ou, *) ' Birth Region Bins:' // trim(string) + write(unit_, *) ' Birth Region Bins:' // trim(string) end if if (t % n_bins(T_ENERGYIN) > 0) then string = "" do i = 1, t % n_bins(T_ENERGYIN) + 1 - string = trim(string) // ' ' // trim(real_to_str(& + string = trim(string) // ' ' // trim(to_str(& t % energy_in(i))) end do - write(ou,*) ' Incoming Energy Bins:' // trim(string) + write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if if (t % n_bins(T_ENERGYOUT) > 0) then string = "" do i = 1, t % n_bins(T_ENERGYOUT) + 1 - string = trim(string) // ' ' // trim(real_to_str(& + string = trim(string) // ' ' // trim(to_str(& t % energy_out(i))) end do - write(ou,*) ' Outgoing Energy Bins:' // trim(string) + write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if if (t % n_macro_bins > 0) then @@ -613,9 +668,9 @@ contains string = trim(string) // ' nu-fission' end select end do - write(ou,*) ' Macro Reactions:' // trim(string) + write(unit_,*) ' Macro Reactions:' // trim(string) end if - write(ou,*) + write(unit_,*) end subroutine print_tally @@ -632,37 +687,129 @@ contains type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() + ! print summary of surfaces + call header("SURFACE SUMMARY", unit=UNIT_SUMMARY) + do i = 1, n_surfaces + s => surfaces(i) + call print_surface(s, unit=UNIT_SUMMARY) + end do + ! print summary of cells - call header("CELL SUMMARY") + call header("CELL SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_cells c => cells(i) - call print_cell(c) + call print_cell(c, unit=UNIT_SUMMARY) end do ! print summary of universes - call header("UNIVERSE SUMMARY") + call header("UNIVERSE SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_universes u => universes(i) - call print_universe(u) + call print_universe(u, unit=UNIT_SUMMARY) end do ! print summary of lattices if (n_lattices > 0) then - call header("LATTICE SUMMARY") + call header("LATTICE SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_lattices l => lattices(i) - call print_lattice(l) + call print_lattice(l, unit=UNIT_SUMMARY) end do end if - ! print summary of surfaces - call header("SURFACE SUMMARY") - do i = 1, n_surfaces - s => surfaces(i) - call print_surface(s) + end subroutine print_geometry + +!=============================================================================== +! PRINT_NUCLIDE displays information about a continuous-energy neutron +! cross_section table and its reactions and secondary angle/energy distributions +!=============================================================================== + + subroutine print_nuclide(nuc, unit) + + type(Nuclide), pointer :: nuc + integer, optional :: unit + + integer :: i + integer :: unit_ + integer :: size_total + integer :: size_xs + integer :: size_angle + integer :: size_energy + type(Reaction), pointer :: rxn => null() + type(UrrData), pointer :: urr => null() + + ! set default unit for writing information + if (present(unit)) then + unit_ = unit + else + unit_ = OUTPUT_UNIT + end if + + ! Determine size of cross-sections + size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8 + size_total = size_xs + + ! Basic nuclide information + write(unit_,*) 'Nuclide ' // trim(nuc % name) + write(unit_,*) ' zaid = ' // trim(to_str(nuc % zaid)) + write(unit_,*) ' awr = ' // trim(to_str(nuc % awr)) + write(unit_,*) ' kT = ' // trim(to_str(nuc % kT)) + write(unit_,*) ' # of grid points = ' // trim(to_str(nuc % n_grid)) + write(unit_,*) ' Fissionable = ', nuc % fissionable + write(unit_,*) ' # of fission reactions = ' // trim(to_str(nuc % n_fission)) + write(unit_,*) ' # of reactions = ' // trim(to_str(nuc % n_reaction)) + write(unit_,*) ' Size of cross sections = ' // trim(to_str(& + size_xs)) // ' bytes' + + write(unit_,*) ' Reaction Q-value Mult IE size(angle) size(energy)' + do i = 1, nuc % n_reaction + ! Information on each reaction + rxn => nuc % reactions(i) + + ! Determine size of angle distribution + if (rxn % has_angle_dist) then + size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8 + else + size_angle = 0 + end if + + ! Determine size of energy distribution + if (rxn % has_energy_dist) then + size_energy = size(rxn % edist % data) * 8 + else + size_energy = 0 + end if + + write(unit_,'(3X,A11,1X,F8.3,2X,I4,2X,I6,1X,I11,1X,I11)') & + reaction_name(rxn % MT), rxn % Q_value, rxn % TY, rxn % IE, & + size_angle, size_energy + + ! Accumulate data size + size_total = size_total + size_angle + size_energy end do - end subroutine print_geometry + ! Write information about URR probability tables + if (nuc % urr_present) then + urr => nuc % urr_data + write(unit_,*) ' Unresolved resonance probability table:' + write(unit_,*) ' # of energies = ' // trim(to_str(urr % n_energy)) + write(unit_,*) ' # of probabilities = ' // trim(to_str(urr % n_prob)) + write(unit_,*) ' Interpolation = ' // trim(to_str(urr % interp)) + write(unit_,*) ' Inelastic flag = ' // trim(to_str(urr % inelastic_flag)) + write(unit_,*) ' Absorption flag = ' // trim(to_str(urr % absorption_flag)) + write(unit_,*) ' Multiply by smooth? ', urr % multiply_smooth + write(unit_,*) ' Min energy = ', trim(to_str(urr % energy(1))) + write(unit_,*) ' Max energy = ', trim(to_str(urr % energy(urr % n_energy))) + end if + + ! Write total memory used + write(unit_,*) ' Total memory used = ' // trim(to_str(size_total)) & + // ' bytes' + + ! Blank line at end of nuclide + write(unit_,*) + + end subroutine print_nuclide !=============================================================================== ! PRINT_SUMMARY displays summary information about the problem about to be run @@ -677,65 +824,63 @@ contains type(TallyObject), pointer :: t => null() ! Display problem summary - call header("PROBLEM SUMMARY") + call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY) if (problem_type == PROB_CRITICALITY) then - write(ou,100) 'Problem type:', 'Criticality' - write(ou,101) 'Number of Cycles:', n_cycles - write(ou,101) 'Number of Inactive Cycles:', n_inactive + write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality' + write(UNIT_SUMMARY,101) 'Number of Cycles:', n_cycles + write(UNIT_SUMMARY,101) 'Number of Inactive Cycles:', n_inactive elseif (problem_type == PROB_SOURCE) then - write(ou,100) 'Problem type:', 'External Source' + write(UNIT_SUMMARY,100) 'Problem type:', 'External Source' end if - write(ou,101) 'Number of Particles:', n_particles + write(UNIT_SUMMARY,101) 'Number of Particles:', n_particles ! Display geometry summary - call header("GEOMETRY SUMMARY") - write(ou,101) 'Number of Cells:', n_cells - write(ou,101) 'Number of Surfaces:', n_surfaces - write(ou,101) 'Number of Materials:', n_materials + call header("GEOMETRY SUMMARY", unit=UNIT_SUMMARY) + write(UNIT_SUMMARY,101) 'Number of Cells:', n_cells + write(UNIT_SUMMARY,101) 'Number of Surfaces:', n_surfaces + write(UNIT_SUMMARY,101) 'Number of Materials:', n_materials ! print summary of all geometry call print_geometry() ! print summary of materials - call header("MATERIAL SUMMARY") + call header("MATERIAL SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_materials m => materials(i) - call print_material(m) + call print_material(m, unit=UNIT_SUMMARY) end do ! print summary of tallies if (n_tallies > 0) then - call header("TALLY SUMMARY") + call header("TALLY SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_tallies t=> tallies(i) - call print_tally(t) + call print_tally(t, unit=UNIT_SUMMARY) end do end if ! print summary of unionized energy grid - call header("UNIONIZED ENERGY GRID") - write(ou,*) "Points on energy grid: " // trim(int_to_str(n_grid)) - write(ou,*) "Extra storage required: " // trim(int_to_str(& + call header("UNIONIZED ENERGY GRID", unit=UNIT_SUMMARY) + write(UNIT_SUMMARY,*) "Points on energy grid: " // trim(to_str(n_grid)) + write(UNIT_SUMMARY,*) "Extra storage required: " // trim(to_str(& n_grid*n_nuclides_total*4)) // " bytes" ! print summary of variance reduction - call header("VARIANCE REDUCTION") + call header("VARIANCE REDUCTION", unit=UNIT_SUMMARY) if (survival_biasing) then - write(ou,100) "Survival Biasing:", "on" + write(UNIT_SUMMARY,100) "Survival Biasing:", "on" else - write(ou,100) "Survival Biasing:", "off" + write(UNIT_SUMMARY,100) "Survival Biasing:", "off" end if - string = real_to_str(weight_cutoff) - write(ou,100) "Weight Cutoff:", trim(string) - string = real_to_str(weight_survive) - write(ou,100) "Survival weight:", trim(string) - write(ou,*) + string = to_str(weight_cutoff) + write(UNIT_SUMMARY,100) "Weight Cutoff:", trim(string) + string = to_str(weight_survive) + write(UNIT_SUMMARY,100) "Survival weight:", trim(string) ! Format descriptor for columns 100 format (1X,A,T35,A) 101 format (1X,A,T35,I11) - end subroutine print_summary !=============================================================================== @@ -748,16 +893,16 @@ contains call header("PLOTTING SUMMARY") ! Print plotting origin - write(ou,100) "Plotting Origin:", trim(real_to_str(plot_origin(1))) // & - " " // trim(real_to_str(plot_origin(2))) // " " // & - trim(real_to_str(plot_origin(3))) + write(ou,100) "Plotting Origin:", trim(to_str(plot_origin(1))) // & + " " // trim(to_str(plot_origin(2))) // " " // & + trim(to_str(plot_origin(3))) ! Print plotting width - write(ou,100) "Plotting Width:", trim(real_to_str(plot_width(1))) // & - " " // trim(real_to_str(plot_width(2))) + write(ou,100) "Plotting Width:", trim(to_str(plot_width(1))) // & + " " // trim(to_str(plot_width(2))) ! Print pixel width - write(ou,100) "Pixel Width:", trim(real_to_str(pixel)) + write(ou,100) "Pixel Width:", trim(to_str(pixel)) write(ou,*) ! Format descriptor for columns @@ -799,7 +944,7 @@ contains ! display calculate rate and final keff total_particles = n_particles * n_cycles speed = real(total_particles) / time_compute % elapsed - string = real_to_str(speed) + string = to_str(speed) write(ou,101) "Calculation Rate", trim(string) write(ou,102) "Final Keff", keff, keff_std write(ou,*) @@ -811,4 +956,30 @@ contains end subroutine print_runtime +!=============================================================================== +! CREATE_SUMMARY_FILE opens the summary.out file for logging information about +! the simulation +!=============================================================================== + + subroutine create_summary_file() + + integer :: io_error + logical :: file_exists ! does log file already exist? + character(MAX_FILE_LEN) :: path ! path of summary file + + ! Create filename for log file + path = "summary.out" + + ! Check if log file already exists + inquire(FILE=path, EXIST=file_exists) + if (file_exists) then + ! Possibly copy old log file + end if + + ! Open log file for writing + open(UNIT=UNIT_SUMMARY, FILE=path, STATUS='replace', & + ACTION='write', IOSTAT=io_error) + + end subroutine create_summary_file + end module output diff --git a/src/physics.F90 b/src/physics.F90 index 04c829ff86..869acc367a 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,21 +1,22 @@ module physics + use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section_header, only: Nuclide, Reaction, DistEnergy - use endf, only: reaction_name, is_fission, is_scatter - use error, only: fatal_error, warning - use fission, only: nu_total, nu_prompt, nu_delayed - use geometry, only: find_cell, distance_to_boundary, & - cross_surface, cross_lattice - use geometry_header, only: Universe, BASE_UNIVERSE + use cross_section, only: calculate_xs + use endf, only: reaction_name, is_fission, is_scatter + use error, only: fatal_error, warning + use fission, only: nu_total, nu_prompt, nu_delayed + use geometry, only: find_cell, distance_to_boundary, cross_surface, & + cross_lattice + use geometry_header, only: Universe, BASE_UNIVERSE use global - use interpolation, only: interpolate_tab1 - use output, only: write_message - use particle_header, only: Particle, LocalCoord - use random_lcg, only: prn - use search, only: binary_search - use string, only: int_to_str, real_to_str - use tally, only: score_tally, score_surface_current + use interpolation, only: interpolate_tab1 + use output, only: write_message + use particle_header, only: Particle, LocalCoord + use random_lcg, only: prn + use search, only: binary_search + use string, only: to_str + use tally, only: score_tally, score_surface_current implicit none @@ -29,22 +30,21 @@ contains type(Particle), pointer :: p - integer :: surface_crossed ! surface which particle is on - integer :: last_cell ! most recent cell particle was in - integer :: n_event ! number of collisions/crossings - real(8) :: d_boundary ! distance to nearest boundary - real(8) :: d_collision ! sampled distance to collision - real(8) :: distance ! distance particle travels - logical :: found_cell ! found cell which particle is in? - logical :: lattice_crossed ! is surface crossing in lattice? + integer :: surface_crossed ! surface which particle is on + integer :: lattice_crossed ! lattice boundary which particle crossed + integer :: last_cell ! most recent cell particle was in + integer :: n_event ! number of collisions/crossings + real(8) :: d_boundary ! distance to nearest boundary + real(8) :: d_collision ! sampled distance to collision + real(8) :: distance ! distance particle travels + logical :: found_cell ! found cell which particle is in? type(LocalCoord), pointer :: coord => null() if (p % coord % cell == NONE) then call find_cell(p, found_cell) ! Particle couldn't be located if (.not. found_cell) then - write(message, '(A,3ES11.3)') & - "Could not locate cell for particle at: ", p % coord0 % xyz + message = "Could not locate particle " // trim(to_str(p % id)) call fatal_error() end if @@ -53,12 +53,12 @@ contains end if if (verbosity >= 9 .or. trace) then - message = "Simulating Particle " // trim(int_to_str(p % id)) + message = "Simulating Particle " // trim(to_str(p % id)) call write_message() end if if (verbosity >= 10 .or. trace) then - message = " Born in cell " // trim(int_to_str(& + message = " Born in cell " // trim(to_str(& cells(p % coord % cell) % id)) call write_message() end if @@ -94,9 +94,9 @@ contains if (d_collision > d_boundary) then last_cell = p % coord % cell p % coord % cell = NONE - if (lattice_crossed) then + if (lattice_crossed /= NONE) then p % surface = NONE - call cross_lattice(p) + call cross_lattice(p, lattice_crossed) else p % surface = surface_crossed call cross_surface(p, last_cell) @@ -125,8 +125,8 @@ contains ! If particle has too many events, display warning and kill it n_event = n_event + 1 if (n_event == MAX_EVENTS) then - message = "Particle " // trim(int_to_str(p%id)) // " underwent " & - // "maximum number of events." + message = "Particle " // trim(to_str(p%id)) // " underwent maximum & + &number of events." call warning() p % alive = .false. end if @@ -135,259 +135,6 @@ contains end subroutine transport -!=============================================================================== -! CALCULATE_XS determines the macroscopic cross sections for the material the -! particle is currently traveling through. -!=============================================================================== - - subroutine calculate_xs(p) - - type(Particle), pointer :: p - - integer :: i ! loop index over nuclides - integer :: index_nuclide ! index into nuclides array - integer :: index_sab ! index into sab_tables array - real(8) :: atom_density ! atom density of a nuclide - real(8) :: sab_threshold ! threshold for S(a,b) table - type(Material), pointer :: mat => null() ! current material - - ! Set all material macroscopic cross sections to zero - material_xs % total = ZERO - material_xs % elastic = ZERO - material_xs % absorption = ZERO - material_xs % fission = ZERO - material_xs % nu_fission = ZERO - - mat => materials(p % material) - - ! Find energy index on unionized grid - call find_energy_index(p) - - ! Check if there's an S(a,b) table for this material - if (mat % has_sab_table) then - sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic - else - sab_threshold = ZERO - end if - - ! Add contribution from each nuclide in material - do i = 1, mat % n_nuclides - ! Determine microscopic cross sections for this nuclide - index_nuclide = mat % nuclide(i) - - ! Determine whether to use S(a,b) based on energy of particle and whether - ! the nuclide matches the S(a,b) table - if (p % E < sab_threshold .and. i == mat % sab_nuclide) then - index_sab = mat % sab_table - else - index_sab = 0 - end if - - ! Calculate microscopic cross section for this nuclide - if (p % E /= micro_xs(index_nuclide) % last_E) then - call calculate_nuclide_xs(p, index_nuclide, index_sab) - end if - - ! Copy atom density of nuclide in material - atom_density = mat % atom_density(i) - - ! Add contributions to material macroscopic total cross section - material_xs % total = material_xs % total + & - atom_density * micro_xs(index_nuclide) % total - - ! Add contributions to material macroscopic scattering cross section - material_xs % elastic = material_xs % elastic + & - atom_density * micro_xs(index_nuclide) % elastic - - ! Add contributions to material macroscopic absorption cross section - material_xs % absorption = material_xs % absorption + & - atom_density * micro_xs(index_nuclide) % absorption - - ! Add contributions to material macroscopic fission cross section - material_xs % fission = material_xs % fission + & - atom_density * micro_xs(index_nuclide) % fission - - ! Add contributions to material macroscopic nu-fission cross section - material_xs % nu_fission = material_xs % nu_fission + & - atom_density * micro_xs(index_nuclide) % nu_fission - end do - - end subroutine calculate_xs - -!=============================================================================== -! CALCULATE_NUCLIDE_XS determines microscopic cross sections for a nuclide of a -! given index in the nuclides array at the energy of the given particle -!=============================================================================== - - subroutine calculate_nuclide_xs(p, index_nuclide, index_sab) - - type(Particle), pointer :: p - integer, intent(in) :: index_nuclide ! index into nuclides array - integer, intent(in) :: index_sab ! index into sab_tables array - - integer :: IE ! index on nuclide energy grid - integer :: IE_sab ! index on S(a,b) energy grid - real(8) :: f ! interp factor on nuclide energy grid - real(8) :: f_sab ! interp factor on S(a,b) energy grid - real(8) :: inelastic ! S(a,b) inelastic cross section - real(8) :: elastic ! S(a,b) elastic cross section - type(Nuclide), pointer :: nuc => null() - type(SAB_Table), pointer :: sab => null() - - ! Set pointer to nuclide - nuc => nuclides(index_nuclide) - - ! TODO: If not using unionized energy grid, we need to find the index on the - ! nuclide energy grid using lethargy mapping or whatever other technique - - ! search nuclide energy grid - IE = nuc % grid_index(p % IE) - f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE)) - - micro_xs(index_nuclide) % index_grid = IE - micro_xs(index_nuclide) % interp_factor = f - - ! Initialize sab treatment to false - micro_xs(index_nuclide) % use_sab = .false. - micro_xs(index_nuclide) % elastic_sab = ZERO - - ! Initialize nuclide cross-sections to zero - micro_xs(index_nuclide) % fission = ZERO - micro_xs(index_nuclide) % nu_fission = ZERO - - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % total = & - (ONE-f) * nuc % total(IE) + f * nuc % total(IE+1) - - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % elastic = & - (ONE-f) * nuc % elastic(IE) + f * nuc % elastic(IE+1) - - ! Calculate microscopic nuclide absorption cross section - micro_xs(index_nuclide) % absorption = & - (ONE-f) * nuc % absorption(IE) + f * nuc % absorption(IE+1) - - if (nuc % fissionable) then - ! Calculate microscopic nuclide total cross section - micro_xs(index_nuclide) % fission = & - (ONE-f) * nuc % fission(IE) + f * nuc % fission(IE+1) - - ! Calculate microscopic nuclide nu-fission cross section - micro_xs(index_nuclide) % nu_fission = & - (ONE-f) * nuc % nu_fission(IE) + f * nuc % nu_fission(IE+1) - end if - - ! If there is S(a,b) data for this nuclide, we need to do a few - ! things. Since the total cross section was based on non-S(a,b) data, we - ! need to correct it by subtracting the non-S(a,b) elastic cross section and - ! then add back in the calculated S(a,b) elastic+inelastic cross section. - - if (index_sab > 0) then - micro_xs(index_nuclide) % use_sab = .true. - - ! Get pointer to S(a,b) table - sab => sab_tables(index_sab) - - ! Get index and interpolation factor for inelastic grid - if (p%E < sab % inelastic_e_in(1)) then - IE_sab = 1 - f_sab = ZERO - else - IE_sab = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, p%E) - f_sab = (p%E - sab%inelastic_e_in(IE_sab)) / & - (sab%inelastic_e_in(IE_sab+1) - sab%inelastic_e_in(IE_sab)) - end if - - ! Calculate S(a,b) inelastic scattering cross section - inelastic = (ONE-f_sab) * sab % inelastic_sigma(IE_sab) + f_sab * & - sab % inelastic_sigma(IE_sab + 1) - - ! Check for elastic data - if (p % E < sab % threshold_elastic) then - ! Determine whether elastic scattering is given in the coherent or - ! incoherent approximation. For coherent, the cross section is - ! represented as P/E whereas for incoherent, it is simply P - - if (sab % elastic_mode == SAB_ELASTIC_EXACT) then - if (p % E < sab % elastic_e_in(1)) then - ! If energy is below that of the lowest Bragg peak, the elastic - ! cross section will be zero - elastic = ZERO - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - elastic = sab % elastic_P(IE_sab) / p % E - end if - else - ! Determine index on elastic energy grid - if (p % E < sab % elastic_e_in(1)) then - IE_sab = 1 - else - IE_sab = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, p%E) - end if - - ! Get interpolation factor for elastic grid - f_sab = (p%E - sab%elastic_e_in(IE_sab))/(sab%elastic_e_in(IE_sab+1) - & - sab%elastic_e_in(IE_sab)) - - ! Calculate S(a,b) elastic scattering cross section - elastic = (ONE-f_sab) * sab % elastic_P(IE_sab) + f_sab * & - sab % elastic_P(IE_sab + 1) - end if - else - ! No elastic data - elastic = ZERO - end if - - ! Correct total and elastic cross sections - micro_xs(index_nuclide) % total = micro_xs(index_nuclide) % total - & - micro_xs(index_nuclide) % elastic + inelastic + elastic - micro_xs(index_nuclide) % elastic = inelastic + elastic - - ! Store S(a,b) elastic cross section for sampling later - micro_xs(index_nuclide) % elastic_sab = elastic - end if - - ! Set last evaluated energy - micro_xs(index_nuclide) % last_E = p % E - - end subroutine calculate_nuclide_xs - -!=============================================================================== -! FIND_ENERGY_INDEX determines the index on the union energy grid and the -! interpolation factor for a particle at a certain energy -!=============================================================================== - - subroutine find_energy_index(p) - - type(Particle), pointer :: p - - integer :: IE ! index on union energy grid - real(8) :: E ! energy of particle - real(8) :: interp ! interpolation factor - - ! copy particle's energy - E = p % E - - ! if particle's energy is outside of energy grid range, set to first or last - ! index. Otherwise, do a binary search through the union energy grid. - if (E < e_grid(1)) then - IE = 1 - elseif (E > e_grid(n_grid)) then - IE = n_grid - 1 - else - IE = binary_search(e_grid, n_grid, E) - end if - - ! calculate the interpolation factor -- note this will be outside of [0,1) - ! for a particle outside the energy range of the union grid - interp = (E - e_grid(IE))/(e_grid(IE+1) - e_grid(IE)) - - ! set particle attributes - p % IE = IE - p % interp = interp - - end subroutine find_energy_index - !=============================================================================== ! COLLISION samples a nuclide and reaction and then calls the appropriate ! routine for that reaction @@ -420,7 +167,7 @@ contains ! Display information about collision if (verbosity >= 10 .or. trace) then message = " " // trim(reaction_name(MT)) // ". Energy = " // & - trim(real_to_str(p % E * 1e6_8)) // " eV." + trim(to_str(p % E * 1e6_8)) // " eV." call write_message() end if @@ -450,9 +197,6 @@ contains ! Reset number of particles banked during collision p % n_bank = 0 - ! find energy index, interpolation factor - call find_energy_index(p) - end subroutine collision !=============================================================================== @@ -571,7 +315,7 @@ contains ! add to cumulative probability prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & - & + f*(rxn%sigma(IE-rxn%IE+2))) + + f*(rxn%sigma(IE-rxn%IE+2))) end do else ! For nuclides with only total fission reaction, get a pointer to @@ -596,7 +340,7 @@ contains ! add to cumulative probability if (nuc % has_partial_fission) then prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & - & + f*(rxn%sigma(IE-rxn%IE+2))) + + f*(rxn%sigma(IE-rxn%IE+2))) else prob = prob + micro_xs(index_nuclide) % fission end if @@ -672,7 +416,7 @@ contains if (i > nuc % n_reaction) then message = "Did not sample any reaction for nuclide " // & trim(nuc % name) // " on material " // & - trim(int_to_str(mat % id)) + trim(to_str(mat % id)) call fatal_error() end if @@ -692,7 +436,7 @@ contains ! add to cumulative probability prob = prob + ((ONE-f)*rxn%sigma(IE-rxn%IE+1) & - & + f*(rxn%sigma(IE-rxn%IE+2))) + + f*(rxn%sigma(IE-rxn%IE+2))) end do ! Perform collision physics for inelastics scattering @@ -1023,7 +767,7 @@ contains ! determine velocity vector of target nucleus v_target(1) = u*vt v_target(2) = v*vt - v_target(3) = W*vt + v_target(3) = w*vt end subroutine sample_target_velocity @@ -1262,7 +1006,7 @@ contains ! determine outgoing energy in lab E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & - & / ((A+ONE)*(A+ONE)) + / ((A+ONE)*(A+ONE)) ! determine outgoing angle in lab mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) @@ -1338,7 +1082,7 @@ contains else i = binary_search(rxn % adist % energy, n, E) r = (E - rxn % adist % energy(i)) / & - & (rxn % adist % energy(i+1) - rxn % adist % energy(i)) + (rxn % adist % energy(i+1) - rxn % adist % energy(i)) end if ! Sample between the ith and (i+1)th bin @@ -1380,8 +1124,7 @@ contains mu = mu0 + (xi - c_k)/p0 elseif (interp == LINEAR_LINEAR) then - ! Linear-linear interpolation -- not sure how you come about the - ! formula given in the MCNP manual + ! Linear-linear interpolation p1 = rxn % adist % data(lc + NP + k+1) mu1 = rxn % adist % data(lc + k+1) @@ -1392,7 +1135,7 @@ contains mu = mu0 + (sqrt(p0*p0 + 2*frac*(xi - c_k))-p0)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(interp)) + message = "Unknown interpolation type: " // trim(to_str(interp)) call fatal_error() end if @@ -1404,7 +1147,7 @@ contains end if else - message = "Unknown angular distribution type: " // trim(int_to_str(type)) + message = "Unknown angular distribution type: " // trim(to_str(type)) call fatal_error() end if @@ -1559,7 +1302,7 @@ contains lc = 2 + 2*NR i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) ! Sample outgoing energy bin r1 = prn() @@ -1635,7 +1378,7 @@ contains else i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1708,10 +1451,10 @@ contains E_out = E_l_k + (r1 - c_k)/p_l_k else E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - & - & p_l_k)/frac + p_l_k)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(INTT)) + message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -1861,7 +1604,7 @@ contains else i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) end if ! Sample between the ith and (i+1)th bin @@ -1940,7 +1683,7 @@ contains E_out = E_l_k + (r1 - c_k)/p_l_k else E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - & - & p_l_k)/frac + p_l_k)/frac end if ! Determine Kalbach-Mann parameters @@ -1952,7 +1695,7 @@ contains 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 - message = "Unknown interpolation type: " // trim(int_to_str(INTT)) + message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -2004,7 +1747,7 @@ contains else i = binary_search(edist % data(lc+1), NE, E_in) r = (E_in - edist%data(lc+i)) / & - & (edist%data(lc+i+1) - edist%data(lc+i)) + (edist%data(lc+i+1) - edist%data(lc+i)) end if ! Sample between the ith and (i+1)th bin @@ -2079,10 +1822,10 @@ contains E_out = E_l_k + (r1 - c_k)/p_l_k else E_out = E_l_k + (sqrt(p_l_k*p_l_k + 2*frac*(r1 - c_k)) - & - & p_l_k)/frac + p_l_k)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(INTT)) + message = "Unknown interpolation type: " // trim(to_str(INTT)) call fatal_error() end if @@ -2135,7 +1878,7 @@ contains mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac end if else - message = "Unknown interpolation type: " // trim(int_to_str(JJ)) + message = "Unknown interpolation type: " // trim(to_str(JJ)) call fatal_error() end if diff --git a/src/plot.F90 b/src/plot.F90 index a42c2ac4ed..77ae883c03 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -23,6 +23,7 @@ contains integer :: i ! loop index integer :: surface_crossed ! surface which particle is on + integer :: lattice_crossed ! is surface crossing in lattice? integer :: last_cell ! most recent cell particle was in integer :: enter_surface ! entrance surface real(8) :: xyz(3) ! starting coordinates @@ -31,7 +32,6 @@ contains real(8) :: d ! distance to boundary real(8) :: distance ! distance particle travels logical :: found_cell ! found cell which particle is in? - logical :: lattice_crossed ! is surface crossing in lattice? character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file type(Cell), pointer :: c => null() type(Universe), pointer :: univ => null() @@ -161,9 +161,9 @@ contains write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell p % coord % cell = 0 - if (lattice_crossed) then + if (lattice_crossed /= NONE) then p % surface = NONE - call cross_lattice(p) + call cross_lattice(p, lattice_crossed) else p % surface = surface_crossed call cross_surface(p, last_cell) diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 8b7225c7ff..674be52122 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -23,7 +23,7 @@ module random_lcg contains !=============================================================================== -! PRN generates a pseudo-random number using a linean congruential generator +! PRN generates a pseudo-random number using a linear congruential generator !=============================================================================== function prn() result(pseudo_rn) diff --git a/src/source.F90 b/src/source.F90 index aa326be32a..6abee5be01 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,15 +1,15 @@ module source - use bank_header, only: Bank - use constants, only: ONE, MAX_LINE_LEN - use cross_section_header, only: Nuclide - use error, only: fatal_error + use ace_header, only: Nuclide + use bank_header, only: Bank + use constants, only: ONE, MAX_LINE_LEN + use error, only: fatal_error use global - use output, only: write_message - use particle_header, only: Particle, initialize_particle - use physics, only: watt_spectrum - use random_lcg, only: prn, set_particle_seed - use string, only: int_to_str + use output, only: write_message + use particle_header, only: Particle, initialize_particle + use physics, only: watt_spectrum + use random_lcg, only: prn, set_particle_seed + use string, only: to_str implicit none @@ -51,7 +51,7 @@ contains bytes = maxwork * 64 / 8 !#endif message = "Could not allocate source bank. Attempted to allocate " & - // trim(int_to_str(bytes)) // " bytes." + // trim(to_str(bytes)) // " bytes." call fatal_error() end if @@ -64,7 +64,7 @@ contains bytes = 3 * maxwork * 64 / 8 !#endif message = "Could not allocate fission bank. Attempted to allocate " & - // trim(int_to_str(bytes)) // " bytes." + // trim(to_str(bytes)) // " bytes." call fatal_error() end if @@ -72,6 +72,10 @@ contains if (external_source%type == SRC_BOX) then p_min = external_source%values(1:3) p_max = external_source%values(4:6) + else + message = "Unsupported external source type: " // & + to_str(external_source%type) + call fatal_error() end if ! Initialize first cycle source bank diff --git a/src/string.F90 b/src/string.F90 index 839abdea38..2800d84c12 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -6,8 +6,8 @@ module string implicit none - interface int_to_str - module procedure int4_to_str, int8_to_str + interface to_str + module procedure int4_to_str, int8_to_str, real_to_str end interface contains @@ -348,7 +348,7 @@ contains ! Determine number of digits to right of decimal point index_decimal = index(string, '.') index_exponent = max(index(string, 'd'), index(string, 'D'), & - & index(string, 'e'), index(string, 'E')) + index(string, 'e'), index(string, 'E')) if (index_decimal > 0) then if (index_exponent > 0) then d = index_exponent - index_decimal - 1 diff --git a/src/tally.F90 b/src/tally.F90 index 547f212690..a1d66674ec 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,20 +1,18 @@ module tally - use ISO_FORTRAN_ENV - use constants use error, only: fatal_error use global use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices use mesh_header, only: StructuredMesh - use output, only: write_message, header + use output, only: header use search, only: binary_search - use string, only: int_to_str, real_to_str + use string, only: to_str use tally_header, only: TallyScore, TallyMapItem, TallyMapElement #ifdef MPI + use intercycle, only: reduce_tallies use mpi - use mpi_routines, only: reduce_tallies #endif implicit none @@ -23,81 +21,6 @@ module tally contains -!=============================================================================== -! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the -! mean and standard deviation of the mean for active cycles and displays them -!=============================================================================== - - subroutine calculate_keff(i_cycle) - - integer, intent(in) :: i_cycle ! index of current cycle - - integer(8) :: total_bank ! total number of source sites - integer :: n ! active cycle number - real(8) :: k_cycle ! single cycle estimate of keff - real(8), save :: k_sum ! accumulated keff - real(8), save :: k_sum_sq ! accumulated keff**2 - - message = "Calculate cycle keff..." - call write_message(8) - - ! initialize sum and square of sum at beginning of run - if (i_cycle == 1) then - k_sum = ZERO - k_sum_sq = ZERO - end if - -#ifdef MPI - ! Collect number bank sites onto master process - call MPI_REDUCE(n_bank, total_bank, 1, MPI_INTEGER8, MPI_SUM, 0, & - & MPI_COMM_WORLD, mpi_err) -#else - total_bank = n_bank -#endif - - ! Collect statistics and print output - if (master) then - ! Since the creation of bank sites was originally weighted by the last - ! cycle keff, we need to multiply by that keff to get the current cycle's - ! value - - k_cycle = real(total_bank)/real(n_particles)*keff - - if (i_cycle > n_inactive) then - ! Active cycle number - n = i_cycle - n_inactive - - ! Accumulate cycle estimate of k - k_sum = k_sum + k_cycle - k_sum_sq = k_sum_sq + k_cycle*k_cycle - - ! Determine mean and standard deviation of mean - keff = k_sum/n - keff_std = sqrt((k_sum_sq/n - keff*keff)/n) - - ! Display output for this cycle - if (i_cycle > n_inactive+1) then - write(UNIT=OUTPUT_UNIT, FMT=101) i_cycle, k_cycle, keff, keff_std - else - write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle - end if - else - ! Display output for inactive cycle - write(UNIT=OUTPUT_UNIT, FMT=100) i_cycle, k_cycle - keff = k_cycle - end if - end if - -#ifdef MPI - ! Broadcast new keff value to all processors - call MPI_BCAST(keff, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) -#endif - -100 format (2X,I5,2X,F8.5) -101 format (2X,I5,2X,F8.5,9X,F8.5,1X,F8.5) - - end subroutine calculate_keff - !=============================================================================== ! CREATE_TALLY_MAP creates a map that allows a quick determination of which ! tallies and bins need to be scored to when a particle makes a collision. This @@ -1043,7 +966,7 @@ contains t => tallies(i) ! Write header block - call header("TALLY " // trim(int_to_str(t % id)), 3, UNIT_TALLY) + call header("TALLY " // trim(to_str(t % id)), unit=UNIT_TALLY, level=3) ! Handle surface current tallies separately if (t % surface_current) then @@ -1125,8 +1048,8 @@ contains do k = 1, t % n_macro_bins write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & repeat(" ", indent), macro_name(abs(t % macro_bins(k) % scalar)), & - real_to_str(t % scores(score_index,k) % val), & - trim(real_to_str(t % scores(score_index,k) % val_sq)) + to_str(t % scores(score_index,k) % val), & + trim(to_str(t % scores(score_index,k) % val_sq)) end do indent = indent - 2 @@ -1176,14 +1099,14 @@ contains end if do i = 1, m % dimension(1) - string = "Mesh Index (" // trim(int_to_str(i)) // ", " + string = "Mesh Index (" // trim(to_str(i)) // ", " len1 = len_trim(string) do j = 1, m % dimension(2) - string = string(1:len1+1) // trim(int_to_str(j)) // ", " + string = string(1:len1+1) // trim(to_str(j)) // ", " len2 = len_trim(string) do k = 1, m % dimension(3) ! Write mesh cell index - string = string(1:len2+1) // trim(int_to_str(k)) // ")" + string = string(1:len2+1) // trim(to_str(k)) // ")" write(UNIT=UNIT_TALLY, FMT='(1X,A)') trim(string) do l = 1, n @@ -1202,15 +1125,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Left", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Left", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Right Surface bins(1:3) = (/ i, j, k /) + 1 @@ -1218,15 +1141,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Right", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_RIGHT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Right", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Back Surface bins(1:3) = (/ i, j-1, k /) + 1 @@ -1234,15 +1157,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Back", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Back", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Front Surface bins(1:3) = (/ i, j, k /) + 1 @@ -1250,15 +1173,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Front", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_FRONT score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Front", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Bottom Surface bins(1:3) = (/ i, j, k-1 /) + 1 @@ -1266,15 +1189,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Bottom", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Bottom", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) ! Top Surface bins(1:3) = (/ i, j, k /) + 1 @@ -1282,15 +1205,15 @@ contains score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Incoming Current from Top", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) bins(TS_SURFACE) = OUT_TOP score_index = sum((bins - 1) * t % stride) + 1 write(UNIT=UNIT_TALLY, FMT='(5X,A,T35,A,"+/- ",A)') & "Outgoing Current to Top", & - real_to_str(t % scores(score_index,1) % val), & - trim(real_to_str(t % scores(score_index,1) % val_sq)) + to_str(t % scores(score_index,1) % val), & + trim(to_str(t % scores(score_index,1) % val_sq)) end do end do @@ -1320,38 +1243,38 @@ contains select case(filter_type) case (T_UNIVERSE) i = t % universe_bins(bin) % scalar - label = int_to_str(universes(i) % id) + label = to_str(universes(i) % id) case (T_MATERIAL) i = t % material_bins(bin) % scalar - label = int_to_str(materials(i) % id) + label = to_str(materials(i) % id) case (T_CELL) i = t % cell_bins(bin) % scalar - label = int_to_str(cells(i) % id) + label = to_str(cells(i) % id) case (T_CELLBORN) i = t % cellborn_bins(bin) % scalar - label = int_to_str(cells(i) % id) + label = to_str(cells(i) % id) case (T_SURFACE) i = t % surface_bins(bin) % scalar - label = int_to_str(surfaces(i) % id) + label = to_str(surfaces(i) % id) case (T_MESH) m => meshes(t % mesh) allocate(ijk(m % n_dimension)) call bin_to_mesh_indices(m, bin, ijk) if (m % n_dimension == 2) then - label = "Index (" // trim(int_to_str(ijk(1))) // ", " // & - trim(int_to_str(ijk(2))) // ")" + label = "Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ")" elseif (m % n_dimension == 3) then - label = "Index (" // trim(int_to_str(ijk(1))) // ", " // & - trim(int_to_str(ijk(2))) // ", " // trim(int_to_str(ijk(3))) // ")" + label = "Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" end if case (T_ENERGYIN) E0 = t % energy_in(bin) E1 = t % energy_in(bin + 1) - label = "[" // trim(real_to_str(E0)) // ", " // trim(real_to_str(E1)) // ")" + label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")" case (T_ENERGYOUT) E0 = t % energy_out(bin) E1 = t % energy_out(bin + 1) - label = "[" // trim(real_to_str(E0)) // ", " // trim(real_to_str(E1)) // ")" + label = "[" // trim(to_str(E0)) // ", " // trim(to_str(E1)) // ")" end select end function get_label diff --git a/src/xml-fortran/Makefile b/src/xml-fortran/Makefile index 6e963f0595..c0514533a3 100644 --- a/src/xml-fortran/Makefile +++ b/src/xml-fortran/Makefile @@ -8,6 +8,10 @@ objects = $(source:.f90=.o) # Ignore unused variables +ifeq ($(MACHINE),bluegene) + override F90 = xlf2003 +endif + ifeq ($(F90),ifort) override F90FLAGS += -warn nounused endif diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 71ae91bd27..c69439612e 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -26,5 +26,6 @@ +