From 0adb96d042078548bb01de8b915aa9bd919bee0e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 Mar 2012 11:19:39 -0400 Subject: [PATCH] Cleaned up source, added comments, removed unnecessary variables. --- src/DEPENDENCIES | 4 +- src/ace.F90 | 7 +- src/bank_header.F90 | 4 + src/energy_grid.F90 | 2 +- src/initialize.F90 | 14 ++-- src/input_xml.F90 | 3 +- src/main.F90 | 1 - src/mesh_header.F90 | 14 ++-- src/output.F90 | 190 ++++++++++++++++++++++-------------------- src/physics.F90 | 76 +---------------- src/plot.F90 | 111 +++++++++++++----------- src/plot_header.F90 | 2 +- src/ppmlib.F90 | 66 ++++++++++----- src/search.F90 | 1 - src/source.F90 | 2 +- src/source_header.F90 | 2 +- src/tally.F90 | 6 +- src/timing.F90 | 10 +-- 18 files changed, 245 insertions(+), 270 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index bcc3b533ff..d43e0f3539 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -99,6 +99,7 @@ initialize.o: random_lcg.o initialize.o: source.o initialize.o: string.o initialize.o: tally.o +initialize.o: tally_header.o initialize.o: timing.o input_xml.o: constants.o @@ -143,7 +144,6 @@ main.o: intercycle.o main.o: output.o main.o: physics.o main.o: plot.o -main.o: random_lcg.o main.o: source.o main.o: string.o main.o: tally.o @@ -177,6 +177,7 @@ physics.o: geometry.o physics.o: geometry_header.o physics.o: global.o physics.o: interpolation.o +physics.o: material_header.o physics.o: output.o physics.o: particle_header.o physics.o: random_lcg.o @@ -198,7 +199,6 @@ plot.o: ppmlib.o plot.o: source.o plot.o: string.o -search.o: constants.o search.o: error.o search.o: global.o diff --git a/src/ace.F90 b/src/ace.F90 index 0406c7dbca..af256c9e21 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -22,8 +22,6 @@ module ace 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 @@ -47,6 +45,7 @@ contains type(Material), pointer :: mat => null() type(Nuclide), pointer :: nuc => null() type(SAB_Table), pointer :: sab => null() + type(DictionaryCI), pointer :: already_read => null() ! ========================================================================== ! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES @@ -154,6 +153,8 @@ contains name = xs_listings(index_list) % name alias = xs_listings(index_list) % alias + ! Read the ACE table into the appropriate entry on the nuclides + ! array call read_ace_table(index_nuclides, index_list) ! Print out information on table to cross_sections.out file @@ -172,6 +173,8 @@ contains index_list = dict_get_key(xs_listing_dict, name) index_sab = dict_get_key(sab_dict, name) + ! Read the ACE table into the appropriate entry on the sab_tables + ! array call read_ace_table(index_sab, index_list) call dict_add_key(already_read, name, 0) diff --git a/src/bank_header.F90 b/src/bank_header.F90 index dd255c5f0b..233ff611d4 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -9,7 +9,11 @@ module bank_header !=============================================================================== type Bank + ! The 'sequence' attribute is used here to ensure that the data listed + ! appears in the given order. This is important for MPI purposes when bank + ! sites are sent from one processor to another. sequence + integer(8) :: id ! Unique ID real(8) :: xyz(3) ! location of bank particle real(8) :: uvw(3) ! diretional cosines diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 3fb1d940d2..4828338a2e 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -37,7 +37,7 @@ contains ! Create dictionary for keeping track of cross sections already added call dict_create(already_added) - ! Loop over all files + ! Loop over all materials do i = 1, n_materials mat => materials(i) diff --git a/src/initialize.F90 b/src/initialize.F90 index 6bc2e899c8..c119cd9bab 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -3,13 +3,13 @@ module initialize use ace, only: read_xs use bank_header, only: Bank use constants - use datatypes, only: dict_create, dict_add_key, dict_get_key, & - dict_has_key, dict_keys + use datatypes, only: dict_create, dict_get_key, dict_has_key, & + dict_keys use datatypes_header, only: ListKeyValueII, DictionaryII use energy_grid, only: unionized_grid use error, only: fatal_error use geometry, only: neighbor_lists - use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE + use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE use global use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict, read_plots_xml @@ -19,7 +19,8 @@ module initialize use random_lcg, only: initialize_prng use source, only: initialize_source use string, only: to_str, starts_with, ends_with, lower_case - use tally, only: create_tally_map, TallyObject + use tally, only: create_tally_map + use tally_header, only: TallyObject use timing, only: timer_start, timer_stop #ifdef MPI @@ -33,8 +34,6 @@ module initialize implicit none - type(DictionaryII), pointer :: build_dict => null() - contains !=============================================================================== @@ -286,9 +285,6 @@ contains ! Create special dictionary used in input_xml call dict_create(cells_in_univ_dict) - - ! Create special dictionary for building universes - call dict_create(build_dict) end subroutine create_dictionaries diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e334007fda..039f8cb653 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1,8 +1,7 @@ module input_xml use constants - use datatypes, only: dict_create, dict_add_key, dict_has_key, & - dict_get_key + use datatypes, only: dict_add_key, dict_has_key, dict_get_key use error, only: fatal_error, warning use geometry_header, only: Cell, Surface, Lattice use global diff --git a/src/main.F90 b/src/main.F90 index 31092f79f6..bf5b922bc4 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -8,7 +8,6 @@ program main use output, only: write_message, header use plotter, only: run_plot use physics, only: transport - use random_lcg, only: set_particle_seed use source, only: get_source_particle use string, only: to_str use tally, only: synchronize_tallies, add_to_score diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index b70bbd45b4..8088972ce0 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -8,13 +8,13 @@ module mesh_header !=============================================================================== type StructuredMesh - integer :: id - integer :: type - integer :: n_dimension - integer, allocatable :: dimension(:) - real(8), allocatable :: lower_left(:) - real(8), allocatable :: upper_right(:) - real(8), allocatable :: width(:) + integer :: id ! user-specified id + integer :: type ! rectangular, hexagonal + integer :: n_dimension ! rank of mesh + integer, allocatable :: dimension(:) ! number of cells in each direction + real(8), allocatable :: lower_left(:) ! lower-left corner of mesh + real(8), allocatable :: upper_right(:) ! upper-right corner of mesh + real(8), allocatable :: width(:) ! width of each mesh cell end type StructuredMesh end module mesh_header diff --git a/src/output.F90 b/src/output.F90 index c389fab8d1..2fc6c77e96 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -29,8 +29,8 @@ contains subroutine title() - character(10) :: today_date - character(8) :: today_time + character(10) :: date_ + character(10) :: time_ write(UNIT=OUTPUT_UNIT, FMT='(/11(A/))') & ' .d88888b. 888b d888 .d8888b.', & @@ -55,18 +55,23 @@ contains #endif ! Write the date and time - call get_today(today_date, today_time) + call date_and_time(DATE=date_, TIME=time_) + date_ = date_(1:4) // "-" // date_(5:6) // "-" // date_(7:8) + time_ = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) write(UNIT=OUTPUT_UNIT, FMT='(6X,"Date/Time:",5X,A,1X,A)') & - trim(today_date), trim(today_time) + trim(date_), trim(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" + "Copyright: 2011-2012 Massachusetts Institute of Technology" write(UNIT=UNIT_SUMMARY, FMT='(1X,A,7X,2(I1,"."),I1)') & "Version:", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE +#ifdef GIT_SHA1 + write(UNIT=UNIT_SUMMARY, FMT='(1X,"Git SHA1:",6X,A)') GIT_SHA1 +#endif write(UNIT=UNIT_SUMMARY, FMT='(1X,"Date/Time:",5X,A,1X,A)') & - trim(today_date), trim(today_time) + trim(date_), trim(time_) ! Write information on number of processors #ifdef MPI @@ -89,8 +94,8 @@ contains integer, optional :: unit ! unit to write to integer, optional :: level ! specified header level - integer :: n - integer :: m + integer :: n ! number of = signs on left + integer :: m ! number of = signs on right integer :: unit_ ! unit to write to integer :: header_level ! actual header level character(MAX_LINE_LEN) :: line @@ -140,14 +145,15 @@ contains subroutine write_message(level) - integer, optional :: level + integer, optional :: level ! verbosity level - integer :: n_lines - integer :: i + integer :: n_lines ! number of lines needed + integer :: i ! index for lines ! Only allow master to print to screen if (.not. master .and. present(level)) return + ! TODO: Take care of line wrapping so words don't get cut off if (.not. present(level) .or. level <= verbosity) then n_lines = (len_trim(message)-1)/79 + 1 do i = 1, n_lines @@ -157,55 +163,13 @@ contains end subroutine write_message -!=============================================================================== -! GET_TODAY determines the date and time at which the program began execution -! and returns it in a readable format -!=============================================================================== - - subroutine get_today(today_date, today_time) - - character(10), intent(out) :: today_date - character(8), intent(out) :: today_time - - integer :: val(8) - character(8) :: date_ - character(10) :: time_ - character(5) :: zone - - call date_and_time(date_, time_, zone, val) - ! val(1) = year (YYYY) - ! val(2) = month (MM) - ! val(3) = day (DD) - ! val(4) = timezone - ! val(5) = hours (HH) - ! val(6) = minutes (MM) - ! val(7) = seconds (SS) - ! val(8) = milliseconds - - if (val(2) < 10) then - if (val(3) < 10) then - today_date = date_(6:6) // "/" // date_(8:8) // "/" // date_(1:4) - else - today_date = date_(6:6) // "/" // date_(7:8) // "/" // date_(1:4) - end if - else - if (val(3) < 10) then - today_date = date_(5:6) // "/" // date_(8:8) // "/" // date_(1:4) - else - today_date = date_(5:6) // "/" // date_(7:8) // "/" // date_(1:4) - end if - end if - today_time = time_(1:2) // ":" // time_(3:4) // ":" // time_(5:6) - - end subroutine get_today - !=============================================================================== ! PRINT_PARTICLE displays the attributes of a particle !=============================================================================== subroutine print_particle() - integer :: i + integer :: i ! index for coordinate levels type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() @@ -266,6 +230,7 @@ contains write(ou,*) ' Surface = ' // to_str(sign(s % id, p % surface)) end if + ! Display weight, energy, grid index, and interpolation factor write(ou,*) ' Weight = ' // to_str(p % wgt) write(ou,*) ' Energy = ' // to_str(p % E) write(ou,*) ' IE = ' // to_str(p % IE) @@ -301,27 +266,35 @@ contains subroutine print_cell(c, unit) type(Cell), pointer :: c - integer, optional :: unit + integer, optional :: unit ! specified unit to write to - integer :: temp - integer :: i - integer :: unit_ + integer :: index_cell ! index in cells array + integer :: i ! loop index for surfaces + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Universe), pointer :: u => null() type(Lattice), pointer :: l => null() type(Material), pointer :: m => null() + ! Set unit to stdout if not already set if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id for cell write(unit_,*) 'Cell ' // to_str(c % id) - temp = dict_get_key(cell_dict, c % id) - write(unit_,*) ' Array Index = ' // to_str(temp) + + ! Find index in cells array and write + index_cell = dict_get_key(cell_dict, c % id) + write(unit_,*) ' Array Index = ' // to_str(index_cell) + + ! Write what universe this cell is in u => universes(c % universe) write(unit_,*) ' Universe = ' // to_str(u % id) + + ! Write information on fill for cell select case (c % type) case (CELL_NORMAL) write(unit_,*) ' Fill = NONE' @@ -332,12 +305,16 @@ contains l => lattices(c % fill) write(unit_,*) ' Fill = Lattice ' // to_str(l % id) end select + + ! Write information on material if (c % material == 0) then write(unit_,*) ' Material = NONE' else m => materials(c % material) write(unit_,*) ' Material = ' // to_str(m % id) end if + + ! Write surface specification string = "" do i = 1, c % n_surfaces select case (c % surfaces(i)) @@ -367,24 +344,31 @@ contains type(Universe), pointer :: univ integer, optional :: unit - integer :: i - integer :: unit_ + integer :: i ! loop index for cells in this universe + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() - type(Universe), pointer :: base_u + type(Universe), pointer :: base_u => null() + ! Set default unit to stdout if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Get a pointer to the base universe base_u => universes(BASE_UNIVERSE) + ! Write user-specified id for this universe write(unit_,*) 'Universe ' // to_str(univ % id) + + ! If this is the base universe, indicate so if (associated(univ, base_u)) then write(unit_,*) ' Base Universe' end if + + ! Write list of cells in this universe string = "" do i = 1, univ % n_cells c => cells(univ % cells(i)) @@ -404,14 +388,16 @@ contains type(Lattice), pointer :: lat integer, optional :: unit - integer :: unit_ + integer :: unit_ ! unit to write to + ! set default unit if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write information about lattice 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) @@ -430,19 +416,23 @@ contains subroutine print_surface(surf, unit) type(Surface), pointer :: surf - integer, optional :: unit + integer, optional :: unit ! specified unit to write to - integer :: i - integer :: unit_ + integer :: i ! loop index for coefficients + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string + ! set default unit if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id of surface write(unit_,*) 'Surface ' // to_str(surf % id) + + ! Write type of surface select case (surf % type) case (SURF_PX) string = "X Plane" @@ -469,12 +459,14 @@ contains end select write(unit_,*) ' Type = ' // trim(string) + ! Write coefficients for this surface string = "" do i = 1, size(surf % coeffs) string = trim(string) // ' ' // to_str(surf % coeffs(i), 4) end do write(unit_,*) ' Coefficients = ' // trim(string) + ! Write neighboring cells on positive side of this surface string = "" if (allocated(surf % neighbor_pos)) then do i = 1, size(surf % neighbor_pos) @@ -483,6 +475,7 @@ contains end if write(unit_,*) ' Positive Neighbors = ' // trim(string) + ! Write neighboring cells on negative side of this surface string = "" if (allocated(surf % neighbor_neg)) then do i = 1, size(surf % neighbor_neg) @@ -490,6 +483,8 @@ contains end do end if write(unit_,*) ' Negative Neighbors =' // trim(string) + + ! Write boundary condition for this surface select case (surf % bc) case (BC_TRANSMIT) write(unit_,*) ' Boundary Condition = Transmission' @@ -513,12 +508,13 @@ contains type(Material), pointer :: mat integer, optional :: unit - integer :: i - integer :: unit_ - real(8) :: density + integer :: i ! loop index for nuclides + integer :: unit_ ! unit to write to + real(8) :: density ! density in atom/b-cm character(MAX_LINE_LEN) :: string type(Nuclide), pointer :: nuc => null() + ! set default unit to stdout if not specified if (present(unit)) then unit_ = unit else @@ -559,9 +555,9 @@ contains type(TallyObject), pointer :: t integer, optional :: unit - integer :: i - integer :: id - integer :: unit_ + integer :: i ! index for filter or score bins + integer :: id ! user-specified id + integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() @@ -569,14 +565,17 @@ contains type(Material), pointer :: m => null() type(StructuredMesh), pointer :: sm => null() + ! set default unit to stdout if not specified if (present(unit)) then unit_ = unit else unit_ = OUTPUT_UNIT end if + ! Write user-specified id of tally write(unit_,*) 'Tally ' // to_str(t % id) + ! Write any cells bins if present if (t % n_filter_bins(FILTER_CELL) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_CELL) @@ -587,6 +586,7 @@ contains write(unit_, *) ' Cell Bins:' // trim(string) end if + ! Write any surface bins if present if (t % n_filter_bins(FILTER_SURFACE) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_SURFACE) @@ -597,6 +597,7 @@ contains write(unit_, *) ' Surface Bins:' // trim(string) end if + ! Write any universe bins if present if (t % n_filter_bins(FILTER_UNIVERSE) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_UNIVERSE) @@ -607,6 +608,7 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if + ! Write any material bins if present if (t % n_filter_bins(FILTER_MATERIAL) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_MATERIAL) @@ -617,6 +619,7 @@ contains write(unit_, *) ' Material Bins:' // trim(string) end if + ! Write any mesh bins if present if (t % n_filter_bins(FILTER_MESH) > 0) then string = "" id = t % mesh @@ -628,6 +631,7 @@ contains write(unit_, *) ' Mesh Bins:' // trim(string) end if + ! Write any birth region bins if present if (t % n_filter_bins(FILTER_CELLBORN) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_CELLBORN) @@ -638,6 +642,7 @@ contains write(unit_, *) ' Birth Region Bins:' // trim(string) end if + ! Write any incoming energy bins if present if (t % n_filter_bins(FILTER_ENERGYIN) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_ENERGYIN) + 1 @@ -647,6 +652,7 @@ contains write(unit_,*) ' Incoming Energy Bins:' // trim(string) end if + ! Write any outgoing energy bins if present if (t % n_filter_bins(FILTER_ENERGYOUT) > 0) then string = "" do i = 1, t % n_filter_bins(FILTER_ENERGYOUT) + 1 @@ -656,6 +662,7 @@ contains write(unit_,*) ' Outgoing Energy Bins:' // trim(string) end if + ! Write any score bins if present if (t % n_score_bins > 0) then string = "" do i = 1, t % n_score_bins @@ -674,7 +681,7 @@ contains string = trim(string) // ' nu-fission' end select end do - write(unit_,*) ' Macro Reactions:' // trim(string) + write(unit_,*) ' Scores:' // trim(string) end if write(unit_,*) @@ -687,7 +694,7 @@ contains subroutine print_geometry() - integer :: i + integer :: i ! loop index for various arrays type(Surface), pointer :: s => null() type(Cell), pointer :: c => null() type(Universe), pointer :: u => null() @@ -735,12 +742,12 @@ contains type(Nuclide), pointer :: nuc integer, optional :: unit - integer :: i - integer :: unit_ - integer :: size_total - integer :: size_xs - integer :: size_angle - integer :: size_energy + integer :: i ! loop index over nuclides + integer :: unit_ ! unit to write to + integer :: size_total ! memory used by nuclide (bytes) + integer :: size_xs ! memory used for cross-sections (bytes) + integer :: size_angle ! memory used for angle distributions (bytes) + integer :: size_energy ! memory used for energy distributions (bytes) type(Reaction), pointer :: rxn => null() type(UrrData), pointer :: urr => null() @@ -824,7 +831,7 @@ contains subroutine print_summary() - integer :: i + integer :: i ! loop index character(15) :: string type(Material), pointer :: m => null() type(TallyObject), pointer :: t => null() @@ -896,24 +903,24 @@ contains subroutine print_plot() - integer i - type(Plot), pointer :: pl => null() + integer :: i ! loop index for plots + type(Plot), pointer :: pl => null() ! Display header for plotting call header("PLOTTING SUMMARY") - do i=1,n_plots + do i = 1, n_plots pl => plots(i) - ! Print plot id + ! Write plot id write(ou,100) "Plot ID:", trim(to_str(pl % id)) - ! Print plotting origin + ! Write plotting origin write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // & " " // trim(to_str(pl % origin(2))) // " " // & trim(to_str(pl % origin(3))) - ! Print plotting width + ! Write plotting width if (pl % type == PLOT_TYPE_SLICE) then write(ou,100) "Width:", trim(to_str(pl % width(1))) // & @@ -940,8 +947,8 @@ contains subroutine print_runtime() - integer(8) :: total_particles - real(8) :: speed + integer(8) :: total_particles ! total # of particles simulated + real(8) :: speed ! # of neutrons/second character(15) :: string ! display header block @@ -1017,7 +1024,8 @@ contains end subroutine create_summary_file !=============================================================================== -! CREATE_XS_SUMMARY_FILE +! CREATE_XS_SUMMARY_FILE creates an output file to write information about the +! cross section tables used in the simulation. !=============================================================================== subroutine create_xs_summary_file() diff --git a/src/physics.F90 b/src/physics.F90 index 27ec0adc20..a36202ea72 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -3,7 +3,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants use cross_section, only: calculate_xs - use endf, only: reaction_name, is_fission, is_scatter + use endf, only: reaction_name 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, & @@ -11,6 +11,7 @@ module physics use geometry_header, only: Universe, BASE_UNIVERSE use global use interpolation, only: interpolate_tab1 + use material_header, only: Material use output, only: write_message use particle_header, only: LocalCoord use random_lcg, only: prn @@ -2061,77 +2062,4 @@ contains end function watt_spectrum -!=============================================================================== -! WIGNER samples a Wigner distribution of energy level spacings. Note that this -! scheme is C50 in the Monte Carlo Sampler from Los Alamos (LA-9721-MS). -!=============================================================================== - - function wigner(D_avg) result (D) - - real(8), intent(in) :: D_avg ! average level spacing - real(8) :: D ! sampled level spacing - - real(8) :: c - - c = -4.*D_avg*D_avg/PI * log(prn()) - D = sqrt(c) - - end function wigner - -!=============================================================================== -! CHI_SQUARED samples a chi-squared distribution with n degrees of freedom. The -! distribution of resonance widths in the unresolved region is given by a -! chi-squared distribution. For the special case of n=1, this is a Porter-Thomas -! distribution. For cases with n odd, rule C64 is used whereas for cases with n -! even, rule C45 is used. -!=============================================================================== - - function chi_squared(n, G_avg) result(G) - - integer, intent(in) :: n ! number of degrees of freedom - real(8), intent(in), optional :: G_avg ! average resonance width - - integer :: i ! loop index - real(8) :: G ! sampled random variable (or resonance width) - real(8) :: x, y, c ! dummy variables - real(8) :: r1, r2 ! psuedorandom numbers - - select case (mod(n,2)) - case (0) - ! Even number of degrees of freedom can be sampled via rule C45. We can - ! sample x as -2/n*log(product(r_i, i = 1 to n/2)) - x = ONE - do i = 1, n/2 - x = x * prn() - end do - x = -2./n * log(x) - - case (1) - ! Odd number of degrees of freedom can be sampled via rule C64. We can - ! sample x as -2/n*(log(r)*cos^2(pi/2*r) + log(product(r_i, i = 1 to - ! floor(n/2))) - - ! Note that we take advantage of integer division on n/2 - y = ONE - do i = 1, n/2 - y = y * prn() - end do - - r1 = prn() - r2 = prn() - c = cos(PI/2.*r2) - x = -2./n * (log(y) + log(r1)*c*c) - end select - - ! If sampling a chi-squared distribution for a resonance width and the - ! average resonance width has been given, return the sampled resonance - ! width. - if (present(G_avg)) then - G = x * G_avg - else - G = x - end if - - end function chi_squared - end module physics diff --git a/src/plot.F90 b/src/plot.F90 index 61d0ea74dd..7582130454 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -2,12 +2,11 @@ module plotter use constants use error, only: fatal_error - use geometry, only: find_cell, distance_to_boundary, cross_surface, & - cross_lattice, cell_contains + use geometry, only: find_cell use geometry_header, only: Universe, BASE_UNIVERSE use global use output, only: write_message - use particle_header, only: LocalCoord, deallocate_coord + use particle_header, only: deallocate_coord use plot_header use ppmlib, only: Image, init_image, allocate_image, & deallocate_image, set_pixel @@ -18,17 +17,16 @@ module plotter contains - !=============================================================================== -! RUN_PLOT +! RUN_PLOT controls the logic for making one or many plots !=============================================================================== subroutine run_plot() - integer :: i - type(Plot), pointer :: pl => null() + integer :: i ! loop index for plots + type(Plot), pointer :: pl => null() - do i=1,n_plots + do i = 1, n_plots pl => plots(i) ! Display output message @@ -36,31 +34,34 @@ contains call write_message(5) if (pl % type == PLOT_TYPE_SLICE) then - ! create 2d image call create_ppm(pl) - end if end do end subroutine run_plot !=============================================================================== -! create_ppm +! CREATE_PPM creates an image based on user input from a plots.xml +! specification in the portable pixmap format (PPM) !=============================================================================== + subroutine create_ppm(pl) - type(Plot), pointer :: pl + type(Plot), pointer :: pl - type(Image) :: img - integer :: in_i, out_i - integer :: x, y - integer :: r, g, b - real(8) :: in_pixel, out_pixel - real(8) :: xyz(3) - logical :: found_cell - type(Cell), pointer :: c => null() + integer :: in_i + integer :: out_i + integer :: x, y ! pixel location + integer :: r, g, b ! colors (red, green, blue) from 0-255 + real(8) :: in_pixel + real(8) :: out_pixel + real(8) :: xyz(3) + logical :: found_cell + type(Image) :: img + type(Cell), pointer :: c => null() + ! Initialize and allocate space for image call init_image(img) call allocate_image(img, pl % pixels(1), pl % pixels(2)) @@ -97,8 +98,8 @@ contains p % coord % uvw = (/ 1, 0, 0 /) p % coord % universe = BASE_UNIVERSE - do y=1, img%height - do x=1, img%width + do y = 1, img % height + do x = 1, img % width call deallocate_coord(p % coord0 % next) p % coord => p % coord0 @@ -106,63 +107,74 @@ contains call find_cell(found_cell) if (.not. found_cell) then - r = pl % not_found % rgb(1) - g = pl % not_found % rgb(2) - b = pl % not_found % rgb(3) + ! If no cell, revert to default color + r = pl % not_found % rgb(1) + g = pl % not_found % rgb(2) + b = pl % not_found % rgb(3) else - if (pl % color_by == PLOT_COLOR_MATS) then - c => cells(p % coord % cell) - r = pl % colors(c % material) % rgb(1) - g = pl % colors(c % material) % rgb(2) - b = pl % colors(c % material) % rgb(3) - else if (pl % color_by == PLOT_COLOR_CELLS) then - r = pl % colors(p % coord % cell) % rgb(1) - g = pl % colors(p % coord % cell) % rgb(2) - b = pl % colors(p % coord % cell) % rgb(3) - else - r = 0 - g = 0 - b = 0 - end if + if (pl % color_by == PLOT_COLOR_MATS) then + ! Assign color based on material + c => cells(p % coord % cell) + r = pl % colors(c % material) % rgb(1) + g = pl % colors(c % material) % rgb(2) + b = pl % colors(c % material) % rgb(3) + else if (pl % color_by == PLOT_COLOR_CELLS) then + ! Assign color based on cell + r = pl % colors(p % coord % cell) % rgb(1) + g = pl % colors(p % coord % cell) % rgb(2) + b = pl % colors(p % coord % cell) % rgb(3) + else + r = 0 + g = 0 + b = 0 + end if end if + ! Create a pixel at (x,y) with color (r,g,b) call set_pixel(img, x, y, r, g, b) + ! Advance pixel in first direction p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel end do + ! Advance pixel in second direction p % coord0 % xyz(in_i) = xyz(in_i) p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel end do + ! Write out the ppm to a file call output_ppm(pl,img) + ! Free up space call deallocate_image(img) end subroutine create_ppm - !=============================================================================== -! output_ppm +! OUTPUT_PPM writes out a previously generated image to a PPM file !=============================================================================== - subroutine output_ppm(pl,img) - type(Plot), pointer :: pl - type(Image), intent(in) :: img + subroutine output_ppm(pl, img) - integer :: i, j + type(Plot), pointer :: pl + type(Image), intent(in) :: img + integer :: i ! loop index for height + integer :: j ! loop index for width + + ! Open PPM file for writing open(UNIT=UNIT_PLOT, FILE=pl % path_plot) + ! Write header write(UNIT_PLOT, '(A2)') 'P6' write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height write(UNIT_PLOT, '(A)') '255' - do j=1, img%height - do i=1, img%width + ! Write color for each pixel + do j = 1, img % height + do i = 1, img % width write(UNIT_PLOT, '(3A1)', advance='no') achar(img%red(i,j)), & - achar(img%green(i,j)), & - achar(img%blue(i,j)) + achar(img%green(i,j)), achar(img%blue(i,j)) end do end do @@ -171,5 +183,4 @@ contains end subroutine output_ppm - end module plotter diff --git a/src/plot_header.F90 b/src/plot_header.F90 index c0ae470643..6268ffadf9 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -26,7 +26,7 @@ module plot_header integer :: basis ! direction of plot slice integer :: pixels(2) ! pixel width/height of plot slice type(ObjectColor) :: not_found ! color for positions where no cell found - type(ObjectColor),allocatable :: colors(:) ! colors of cells/mats + type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats end type Plot integer :: PLOT_TYPE_SLICE = 1 diff --git a/src/ppmlib.F90 b/src/ppmlib.F90 index 2c9687617a..d67268412e 100644 --- a/src/ppmlib.F90 +++ b/src/ppmlib.F90 @@ -13,6 +13,10 @@ module ppmlib contains +!=============================================================================== +! INIT_IMAGE initializes the Image derived type +!=============================================================================== + subroutine init_image(img) type(Image), intent(out) :: img @@ -26,29 +30,46 @@ contains end subroutine init_image +!=============================================================================== +! ALLOCATE_IMAGE sets the width and height of an image and allocates color +! arrays +!=============================================================================== + subroutine allocate_image(img, w, h) type(Image), intent(inout) :: img - integer, intent(in) :: w, h + integer, intent(in) :: w ! width of image + integer, intent(in) :: h ! height of image + ! allocate red, green, and blue array allocate(img % red(w, h)) allocate(img % green(w, h)) allocate(img % blue(w, h)) + ! set width and height img % width = w img % height = h end subroutine allocate_image + +!=============================================================================== +! DEALLOCATE_IMAGE +!=============================================================================== subroutine deallocate_image(img) type(Image) :: img - if ( associated(img % red) ) deallocate(img % red) - if ( associated(img % green) ) deallocate(img % green) - if ( associated(img % blue) ) deallocate(img % blue) + if (associated(img % red)) deallocate(img % red) + if (associated(img % green)) deallocate(img % green) + if (associated(img % blue)) deallocate(img % blue) + end subroutine deallocate_image +!=============================================================================== +! INSIDE_IMAGE determines whether a point (x,y) is inside the image +!=============================================================================== + function inside_image(img, x, y) result(inside) @@ -58,13 +79,15 @@ contains inside = .false. - if ( ( x < img % width) .and. & - ( y < img % height ) .and. & - ( x >= 0 ) .and. & - ( y >= 0 ) ) inside = .true. + if ((x < img % width) .and. (y < img % height) .and. & + (x >= 0) .and. (y >= 0)) inside = .true. end function inside_image +!=============================================================================== +! VALID_IMAGE checks whether the image has a width and height and if its color +! arrays are allocated +!=============================================================================== function valid_image(img) result(valid) @@ -73,29 +96,32 @@ contains valid = .false. - if ( img % width == 0 ) return - if ( img % height == 0 ) return - if ( .not. associated(img % red) .or. & + if (img % width == 0) return + if (img % height == 0) return + if (.not. associated(img % red) .or. & .not. associated(img % green) .or. & - .not. associated(img % blue) ) return + .not. associated(img % blue)) return valid = .true. end function valid_image - +!=============================================================================== +! SET_PIXEL sets the +!=============================================================================== + subroutine set_pixel(img, x, y, r, g, b) type(Image), intent(inout) :: img - integer, intent(in) :: x, y - integer, intent(in) :: r, g, b + integer, intent(in) :: x, y ! coordinates + integer, intent(in) :: r, g, b ! red, green, and blue - if ( inside_image(img, x, y) .and. valid_image(img)) then - img % red(x+1,y+1) = mod(abs(r), 256) - img % green(x+1, y+1) = mod(abs(g), 256) - img % blue(x+1, y+1) = mod(abs(b), 256) + if (inside_image(img, x, y) .and. valid_image(img)) then + img % red(x+1,y+1) = mod(abs(r), 256) + img % green(x+1, y+1) = mod(abs(g), 256) + img % blue(x+1, y+1) = mod(abs(b), 256) end if - end subroutine set_pixel + end subroutine set_pixel end module ppmlib diff --git a/src/search.F90 b/src/search.F90 index 3cdb51bfcb..5968113fec 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -1,6 +1,5 @@ module search - use constants, only: ONE, MAX_LINE_LEN use error, only: fatal_error use global, only: message diff --git a/src/source.F90 b/src/source.F90 index 45e9d76996..81d6ab2e6d 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -1,7 +1,7 @@ module source use bank_header, only: Bank - use constants, only: ONE + use constants use error, only: fatal_error use geometry_header, only: BASE_UNIVERSE use global diff --git a/src/source_header.F90 b/src/source_header.F90 index 077e6338d8..638a960244 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -8,7 +8,7 @@ module source_header !=============================================================================== type ExtSource - integer :: type ! type of source, e.g. 'box' or 'cell' + integer :: type ! type of source, e.g. 'box' or 'point' real(8), allocatable :: values(:) ! values for particular source type end type ExtSource diff --git a/src/tally.F90 b/src/tally.F90 index db57a56802..50d29e6061 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1062,7 +1062,8 @@ contains end subroutine add_to_score !=============================================================================== -! ACCUMULATE_BATCH_ESTIMATE +! ACCUMULATE_BATCH_ESTIMATE accumulates scores from many histories (or many +! generations) into a single realization of a random variable. !=============================================================================== elemental subroutine accumulate_batch_estimate(score) @@ -1540,7 +1541,8 @@ contains end function get_label !=============================================================================== -! CALCULATE_STATISTICS +! CALCULATE_STATISTICS determines the sample mean and the standard deviation of +! the mean for a TallyScore. !=============================================================================== elemental subroutine calculate_statistics(score) diff --git a/src/timing.F90 b/src/timing.F90 index f6e91c0583..483057125d 100644 --- a/src/timing.F90 +++ b/src/timing.F90 @@ -13,7 +13,7 @@ module timing type Timer logical :: running = .false. ! is timer running? integer :: start_counts = 0 ! counts when started - real(8) :: elapsed = 0. ! total time elapsed in seconds + real(8) :: elapsed = ZERO ! total time elapsed in seconds !!$ contains !!$ procedure :: start => timer_start !!$ procedure :: get_value => timer_get_value @@ -24,7 +24,7 @@ module timing contains !=============================================================================== -! TIMER_START +! TIMER_START starts running a timer and measures the current time !=============================================================================== subroutine timer_start(self) @@ -38,7 +38,7 @@ contains end subroutine timer_start !=============================================================================== -! TIMER_GET_VALUE +! TIMER_GET_VALUE returns the current value of the timer !=============================================================================== function timer_get_value(self) result(elapsed) @@ -61,7 +61,7 @@ contains end function timer_get_value !=============================================================================== -! TIMER_STOP +! TIMER_STOP stops the timer and sets the elapsed time !=============================================================================== subroutine timer_stop(self) @@ -78,7 +78,7 @@ contains end subroutine timer_stop !=============================================================================== -! TIMER_RESET +! TIMER_RESET resets a timer to have a zero value !=============================================================================== subroutine timer_reset(self)