From a47c7ab91c259fcbc242abaf99f5f8b56b190656 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Mar 2012 10:43:12 -0400 Subject: [PATCH 01/10] Added weight to Bank derived type. --- src/bank_header.F90 | 1 + src/initialize.F90 | 19 ++++++++++--------- src/physics.F90 | 3 +++ src/source.F90 | 7 ++++++- 4 files changed, 20 insertions(+), 10 deletions(-) diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 233ff611d4..10c803d60a 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -15,6 +15,7 @@ module bank_header sequence integer(8) :: id ! Unique ID + real(8) :: wgt ! weight of bank site real(8) :: xyz(3) ! location of bank particle real(8) :: uvw(3) ! diretional cosines real(8) :: E ! energy diff --git a/src/initialize.F90 b/src/initialize.F90 index c119cd9bab..2f50cf93dd 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -149,9 +149,9 @@ contains subroutine setup_mpi() #ifdef MPI - integer :: bank_blocks(4) ! Count for each datatype - integer :: bank_types(4) ! Datatypes - integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements + integer :: bank_blocks(5) ! Count for each datatype + integer :: bank_types(5) ! Datatypes + integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements type(Bank) :: b mpi_enabled = .true. @@ -186,17 +186,18 @@ contains ! 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) + call MPI_GET_ADDRESS(b % wgt, bank_disp(2), mpi_err) + call MPI_GET_ADDRESS(b % xyz, bank_disp(3), mpi_err) + call MPI_GET_ADDRESS(b % uvw, bank_disp(4), mpi_err) + call MPI_GET_ADDRESS(b % E, bank_disp(5), mpi_err) ! Adjust displacements bank_disp = bank_disp - bank_disp(1) ! 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_blocks = (/ 1, 1, 3, 3, 1 /) + bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /) + call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, & bank_types, MPI_BANK, mpi_err) call MPI_TYPE_COMMIT(MPI_BANK, mpi_err) diff --git a/src/physics.F90 b/src/physics.F90 index ca9d389ec1..00c2508ce4 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -906,6 +906,9 @@ contains fission_bank(i) % id = p % id fission_bank(i) % xyz = p % coord0 % xyz + ! Set weight of fission bank site + fission_bank(i) % wgt = ONE + ! sample cosine of angle mu = sample_angle(rxn, E) diff --git a/src/source.F90 b/src/source.F90 index 81d6ab2e6d..654da57af5 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -94,6 +94,9 @@ contains id = bank_first + j - 1 source_bank(j) % id = id + ! Set weight to one + source_bank(j) % wgt = ONE + ! initialize random number seed call set_particle_seed(id) @@ -148,6 +151,8 @@ contains src => source_bank(index_source) ! copy attributes from source bank site + p % wgt = src % wgt + p % last_wgt = src % wgt p % coord % xyz = src % xyz p % coord % uvw = src % uvw p % last_xyz = src % xyz @@ -315,7 +320,7 @@ contains call fatal_error() end if - ! Read position, angle, and energy + ! Read source sites read(UNIT=UNIT_SOURCE) source_bank(1:work) ! Close binary source file From 1fb88c1360ecd6249ffa113e934aa77b2ec5eff1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Mar 2012 12:59:44 -0400 Subject: [PATCH 02/10] Changed count_bank_sites to use weight of bank sites. --- src/intercycle.F90 | 4 ++-- src/mesh.F90 | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index f7026ff858..be5fa749af 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -4,7 +4,7 @@ module intercycle use error, only: fatal_error, warning use global - use mesh, only: count_fission_sites + use mesh, only: count_bank_sites use mesh_header, only: StructuredMesh use output, only: write_message use random_lcg, only: prn, set_particle_seed, prn_skip @@ -320,7 +320,7 @@ contains end if ! count number of fission sites over mesh - call count_fission_sites(m, fission_bank, entropy_p, total, & + call count_bank_sites(m, fission_bank, entropy_p, total, & size_bank=n_bank, sites_outside=sites_outside) ! display warning message if there were sites outside entropy box diff --git a/src/mesh.F90 b/src/mesh.F90 index a9fab038a1..50f916edb0 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -145,13 +145,13 @@ contains end subroutine bin_to_mesh_indices !=============================================================================== -! COUNT_FISSION_SITES determines the number of fission bank sites in each cell -! of a given mesh as well as an optional energy group structure. This can be -! used for a variety of purposes (Shannon entropy, CMFD, uniform fission source +! COUNT_BANK_SITES determines the number of fission bank sites in each cell of a +! given mesh as well as an optional energy group structure. This can be used for +! a variety of purposes (Shannon entropy, CMFD, uniform fission source ! weighting) !=============================================================================== - subroutine count_fission_sites(m, bank_array, cnt, total, & + subroutine count_bank_sites(m, bank_array, cnt, total, & energies, size_bank, sites_outside) type(StructuredMesh), pointer :: m ! mesh to count sites @@ -206,7 +206,7 @@ contains end if ! add weight - weight = weight + ONE + weight = weight + bank_array(i) % wgt ! determine energy bin if (present(energies)) then @@ -222,8 +222,8 @@ contains end if ! add to appropriate mesh box - ! TODO: if tracking weight through bank, add weight instead - cnt(e_bin,ijk(1),ijk(2),ijk(3)) = cnt(e_bin,ijk(1),ijk(2),ijk(3)) + 1 + cnt(e_bin,ijk(1),ijk(2),ijk(3)) = cnt(e_bin,ijk(1),ijk(2),ijk(3)) + & + bank_array(i) % wgt end do FISSION_SITES #ifdef MPI @@ -253,7 +253,7 @@ contains sites_outside = outside #endif - end subroutine count_fission_sites + end subroutine count_bank_sites !=============================================================================== ! MESH_INTERSECT determines if a line between xyz0 and xyz1 intersects the outer From 7819bef5bc514fed390404c07c44631cd5354bb0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Mar 2012 17:22:03 -0400 Subject: [PATCH 03/10] Uniform fission source weighting. Closes #26 on github. --- src/DEPENDENCIES | 1 + src/global.F90 | 5 +++ src/input_xml.F90 | 57 ++++++++++++++++++++++++ src/intercycle.F90 | 51 +++++++++++++++++++++ src/main.F90 | 29 ++++++------ src/mesh_header.F90 | 3 +- src/physics.F90 | 37 +++++++++++++-- src/xml-fortran/templates/settings_t.xml | 5 ++- 8 files changed, 167 insertions(+), 21 deletions(-) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index d6a8591c6c..7d70c783ef 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -179,6 +179,7 @@ physics.o: geometry_header.o physics.o: global.o physics.o: interpolation.o physics.o: material_header.o +physics.o: mesh.o physics.o: output.o physics.o: particle_header.o physics.o: random_lcg.o diff --git a/src/global.F90 b/src/global.F90 index dae61fa6ad..5da8782b5c 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -156,6 +156,11 @@ module global real(8), allocatable :: entropy_p(:,:,:,:) ! % of source sites in each cell type(StructuredMesh), pointer :: entropy_mesh + ! Uniform fission source weighting + logical :: ufs = .false. + type(StructuredMesh), pointer :: ufs_mesh => null() + real(8), allocatable :: source_frac(:,:,:,:) + ! Write source at end of simulation logical :: write_source = .false. diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 039f8cb653..1b7f423c44 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -236,6 +236,60 @@ contains entropy_on = .true. end if + ! Uniform fission source weighting mesh + if (size(uniform_fs_) > 0) then + ! Check to make sure enough values were supplied + if (size(uniform_fs_(1) % lower_left) /= 3) then + message = "Need to specify (x,y,z) coordinates of lower-left corner & + &of UFS mesh." + elseif (size(uniform_fs_(1) % upper_right) /= 3) then + message = "Need to specify (x,y,z) coordinates of upper-right corner & + &of UFS mesh." + elseif (size(uniform_fs_(1) % dimension) /= 3) then + message = "Dimension of UFS mesh must be given as three & + &integers." + call fatal_error() + end if + + ! Allocate mesh object and coordinates on mesh + allocate(ufs_mesh) + allocate(ufs_mesh % lower_left(3)) + allocate(ufs_mesh % upper_right(3)) + allocate(ufs_mesh % width(3)) + + ! Allocate dimensions + ufs_mesh % n_dimension = 3 + allocate(ufs_mesh % dimension(3)) + + ! Copy dimensions + ufs_mesh % dimension = uniform_fs_(1) % dimension + + ! Copy values + ufs_mesh % lower_left = uniform_fs_(1) % lower_left + ufs_mesh % upper_right = uniform_fs_(1) % upper_right + + ! Check on values provided + if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then + message = "Upper-right coordinate must be greater than lower-left & + &coordinate for UFS mesh." + call fatal_error() + end if + + ! Calculate width + ufs_mesh % width = (ufs_mesh % upper_right - & + ufs_mesh % lower_left) / ufs_mesh % dimension + + ! Calculate volume fraction of each cell + ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8) + + ! Turn on uniform fission source weighting + ufs = .true. + + ! Allocate source_frac + allocate(source_frac(1, ufs_mesh % dimension(1), & + ufs_mesh % dimension(2), ufs_mesh % dimension(3))) + end if + ! Check if the user has specified to write binary source file if (trim(write_source_) == 'on') write_source = .true. @@ -812,6 +866,9 @@ contains ! Set upper right coordinate m % upper_right = m % lower_left + m % dimension * m % width + ! Set volume fraction + m % volume_frac = ONE/real(product(m % dimension),8) + ! Add mesh to dictionary call dict_add_key(mesh_dict, m % id, i) end do diff --git a/src/intercycle.F90 b/src/intercycle.F90 index be5fa749af..d44e00c7f3 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -446,4 +446,55 @@ contains end subroutine calculate_keff +!=============================================================================== +! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and +! reweights the source bank so that the sum of the weights is equal to +! n_particles. The 'source_frac' variable is used later to bias the production +! of fission sites +!=============================================================================== + + subroutine count_source_for_ufs() + + real(8) :: total ! total weight in source bank + logical :: sites_outside ! were there sites outside the ufs mesh? +#ifdef MPI + integer :: n ! total number of ufs mesh cells +#endif + + if (current_batch == 1 .and. current_gen == 1) then + ! On the first cycle, just assume that the source is already evenly + ! distributed so that effectively the production of fission sites is not + ! biased + + source_frac = ufs_mesh % volume_frac + + else + ! count number of source sites in each ufs mesh cell + call count_bank_sites(ufs_mesh, source_bank, source_frac, total, & + sites_outside=sites_outside) + + ! Check for sites outside of the mesh + if (master .and. sites_outside) then + message = "Source sites outside of the UFS mesh!" + call fatal_error() + end if + + ! Normalize to total weight to get fraction of source in each cell + if (master) source_frac = source_frac / total + +#ifdef MPI + ! Send source fraction to all processors + n = product(ufs_mesh % dimension) + call MPI_BCAST(source_frac, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) + call MPI_BCAST(total, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) +#endif + + ! Since the total starting weight is not equal to n_particles, we need to + ! renormalize the weight of the source sites + + source_bank % wgt = source_bank % wgt * n_particles / total + end if + + end subroutine count_source_for_ufs + end module intercycle diff --git a/src/main.F90 b/src/main.F90 index bf5b922bc4..9cd90d6d0e 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -2,16 +2,17 @@ program main use constants use global - use finalize, only: finalize_run - use initialize, only: initialize_run - use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank - use output, only: write_message, header - use plotter, only: run_plot - use physics, only: transport - use source, only: get_source_particle - use string, only: to_str - use tally, only: synchronize_tallies, add_to_score - use timing, only: timer_start, timer_stop + use finalize, only: finalize_run + use initialize, only: initialize_run + use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank, & + count_source_for_ufs + use output, only: write_message, header + use plotter, only: run_plot + use physics, only: transport + use source, only: get_source_particle + use string, only: to_str + use tally, only: synchronize_tallies + use timing, only: timer_start, timer_stop #ifdef MPI use mpi @@ -78,6 +79,9 @@ contains ! Set all tallies to zero n_bank = 0 + ! Count source sites if using uniform fission source weighting + if (ufs) call count_source_for_ufs() + ! ==================================================================== ! LOOP OVER HISTORIES @@ -106,11 +110,6 @@ contains ! Distribute fission bank across processors evenly call synchronize_bank() - ! Add to analog estimate of keff -- 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 - call add_to_score(global_tallies(K_ANALOG), n_bank * keff) - ! Stop timer for inter-cycle synchronization call timer_stop(time_intercycle) diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 8088972ce0..a26178b6f3 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -3,7 +3,7 @@ module mesh_header implicit none !=============================================================================== -! STRUCTUREDMESH represents a tesslation of n-dimensional Euclidean space by +! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by ! congruent squares or cubes !=============================================================================== @@ -11,6 +11,7 @@ module mesh_header integer :: id ! user-specified id integer :: type ! rectangular, hexagonal integer :: n_dimension ! rank of mesh + real(8) :: volume_frac ! volume fraction of each cell 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 diff --git a/src/physics.F90 b/src/physics.F90 index 00c2508ce4..382cce60b3 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -12,6 +12,7 @@ module physics use global use interpolation, only: interpolate_tab1 use material_header, only: Material + use mesh, only: get_mesh_indices use output, only: write_message use particle_header, only: LocalCoord use random_lcg, only: prn @@ -839,6 +840,7 @@ contains integer :: nu ! actual number of neutrons produced integer :: law ! energy distribution law integer :: n_sample ! number of times resampling + integer :: ijk(3) ! indices in ufs mesh real(8) :: E ! incoming energy of neutron real(8) :: E_out ! outgoing energy of fission neutron real(8) :: nu_t ! total nu @@ -850,7 +852,9 @@ contains real(8) :: xi ! random number real(8) :: yield ! delayed neutron precursor yield real(8) :: prob ! cumulative probability + real(8) :: weight ! weight adjustment for ufs method logical :: actual_event ! did fission actually occur? (no survival biasing) + logical :: in_mesh ! source site in ufs mesh? type(Nuclide), pointer :: nuc type(DistEnergy), pointer :: edist => null() @@ -886,12 +890,32 @@ contains ! TODO: Heat generation from fission + ! If uniform fission source weighting is turned on, we increase of decrease + ! the expected number of fission sites produced + + if (ufs) then + ! Determine indices on ufs mesh for current location + call get_mesh_indices(ufs_mesh, p % coord0 % xyz, ijk, in_mesh) + if (.not. in_mesh) then + message = "Source site outside UFS mesh!" + call fatal_error() + end if + + if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then + weight = ufs_mesh % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3)) + else + weight = ONE + end if + else + weight = ONE + end if + ! Sample number of neutrons produced if (actual_event) then - nu_t = p % wgt / keff * nu_t + nu_t = p % wgt / keff * nu_t * weight else nu_t = p % last_wgt * micro_xs(index_nuclide) % fission / (keff * & - micro_xs(index_nuclide) % total) * nu_t + micro_xs(index_nuclide) % total) * nu_t * weight end if if (prn() > nu_t - int(nu_t)) then nu = int(nu_t) @@ -899,6 +923,9 @@ contains nu = int(nu_t) + 1 end if + ! Add to analog estimate of keff + call add_to_score(global_tallies(K_ANALOG), nu/weight * keff) + ! Bank source neutrons if (nu == 0 .or. n_bank == 3*work) return do i = int(n_bank,4) + 1, int(min(n_bank + nu, 3*work),4) @@ -907,7 +934,11 @@ contains fission_bank(i) % xyz = p % coord0 % xyz ! Set weight of fission bank site - fission_bank(i) % wgt = ONE + if (ufs) then + fission_bank(i) % wgt = ONE/weight + else + fission_bank(i) % wgt = ONE + end if ! sample cosine of angle mu = sample_angle(rxn, E) diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index f2b387316d..321565955f 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -20,7 +20,7 @@ - + @@ -30,13 +30,14 @@ - + + From cdb2bbbfd1e76fd75a0a5ed094e64211ec9d58be Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Mar 2012 17:31:06 -0400 Subject: [PATCH 04/10] Removed two unnecessary MPI calls relating to count_bank_sites. --- src/intercycle.F90 | 17 ++++++++--------- src/mesh.F90 | 14 ++------------ 2 files changed, 10 insertions(+), 21 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index d44e00c7f3..8957ed4ef2 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -285,8 +285,7 @@ contains integer :: i, j, k ! index for bank sites integer :: n ! # of boxes in each dimension - real(8) :: total ! total weight of fission bank sites - logical :: sites_outside ! were there sites outside entropy box? + logical :: sites_outside ! were there sites outside entropy box? type(StructuredMesh), pointer :: m => null() ! Get pointer to entropy mesh @@ -320,7 +319,7 @@ contains end if ! count number of fission sites over mesh - call count_bank_sites(m, fission_bank, entropy_p, total, & + call count_bank_sites(m, fission_bank, entropy_p, & size_bank=n_bank, sites_outside=sites_outside) ! display warning message if there were sites outside entropy box @@ -332,7 +331,7 @@ contains ! sum values to obtain shannon entropy if (master) then ! Normalize to total weight of bank sites - entropy_p = entropy_p / total + entropy_p = entropy_p / sum(entropy_p) entropy = 0 do i = 1, m % dimension(1) @@ -470,7 +469,7 @@ contains else ! count number of source sites in each ufs mesh cell - call count_bank_sites(ufs_mesh, source_bank, source_frac, total, & + call count_bank_sites(ufs_mesh, source_bank, source_frac, & sites_outside=sites_outside) ! Check for sites outside of the mesh @@ -479,16 +478,16 @@ contains call fatal_error() end if - ! Normalize to total weight to get fraction of source in each cell - if (master) source_frac = source_frac / total - #ifdef MPI ! Send source fraction to all processors n = product(ufs_mesh % dimension) call MPI_BCAST(source_frac, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) - call MPI_BCAST(total, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) #endif + ! Normalize to total weight to get fraction of source in each cell + total = sum(source_frac) + source_frac = source_frac / total + ! Since the total starting weight is not equal to n_particles, we need to ! renormalize the weight of the source sites diff --git a/src/mesh.F90 b/src/mesh.F90 index 50f916edb0..989fec6daa 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -151,14 +151,13 @@ contains ! weighting) !=============================================================================== - subroutine count_bank_sites(m, bank_array, cnt, total, & - energies, size_bank, sites_outside) + subroutine count_bank_sites(m, bank_array, cnt, energies, size_bank, & + sites_outside) type(StructuredMesh), pointer :: m ! mesh to count sites type(Bank), intent(in) :: bank_array(:) ! fission or source bank real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each ! cell and energy group - real(8), intent(out) :: total ! total weight of sites real(8), optional :: energies(:) ! energy grid to search integer(8), optional :: size_bank ! # of bank sites (on each proc) logical, optional :: sites_outside ! were there sites outside mesh? @@ -168,7 +167,6 @@ contains integer :: ijk(3) ! indices on mesh integer :: n_groups ! number of groups in energies integer :: e_bin ! energy_bin - real(8) :: weight ! accumulated weight of sites logical :: in_mesh ! was single site outside mesh? logical :: outside ! was any site outside mesh? #ifdef MPI @@ -177,7 +175,6 @@ contains ! initialize variables cnt = ZERO - weight = ZERO outside = .false. ! Set size of bank @@ -205,9 +202,6 @@ contains cycle end if - ! add weight - weight = weight + bank_array(i) % wgt - ! determine energy bin if (present(energies)) then if (bank_array(i) % E < energies(1)) then @@ -245,11 +239,7 @@ contains MPI_COMM_WORLD, mpi_err) end if - ! determine total weight of bank sites - call MPI_REDUCE(weight, total, 1, MPI_REAL8, MPI_SUM, 0, & - MPI_COMM_WORLD, mpi_err) #else - total = weight sites_outside = outside #endif From 1262759289f61b7b8f541a5ab7d2dbee7d861182 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 1 Apr 2012 12:30:37 -0400 Subject: [PATCH 05/10] Changed names of a few subroutines in tally module to be more consistent. --- src/intercycle.F90 | 4 ++-- src/tally.F90 | 18 +++++++++--------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 8957ed4ef2..c7482c2d00 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -10,7 +10,7 @@ module intercycle use random_lcg, only: prn, set_particle_seed, prn_skip use search, only: binary_search use string, only: to_str - use tally, only: accumulate_batch_estimate + use tally, only: accumulate_score use tally_header, only: TallyObject use timing, only: timer_start, timer_stop @@ -392,7 +392,7 @@ contains n = current_batch - n_inactive ! Accumulate single batch realizations of k - call accumulate_batch_estimate(global_tallies) + call accumulate_score(global_tallies) ! Determine sample mean of keff keff = global_tallies(K_ANALOG) % sum/n diff --git a/src/tally.F90 b/src/tally.F90 index 594d21d353..2c5eb44d8a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1312,11 +1312,11 @@ contains end subroutine add_to_score !=============================================================================== -! ACCUMULATE_BATCH_ESTIMATE accumulates scores from many histories (or many -! generations) into a single realization of a random variable. +! ACCUMULATE_SCORE accumulates scores from many histories (or many generations) +! into a single realization of a random variable. !=============================================================================== - elemental subroutine accumulate_batch_estimate(score) + elemental subroutine accumulate_score(score) type(TallyScore), intent(inout) :: score @@ -1333,7 +1333,7 @@ contains ! Reset the single batch estimate score % value = ZERO - end subroutine accumulate_batch_estimate + end subroutine accumulate_score !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history @@ -1354,7 +1354,7 @@ contains t => tallies(i) ! Loop over all filter and scoring bins - call accumulate_batch_estimate(t % scores) + call accumulate_score(t % scores) end do end subroutine synchronize_tallies @@ -1795,7 +1795,7 @@ contains ! the mean for a TallyScore. !=============================================================================== - elemental subroutine calculate_statistics(score) + elemental subroutine statistics_score(score) type(TallyScore), intent(inout) :: score @@ -1807,7 +1807,7 @@ contains score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / & (n_active - 1)) - end subroutine calculate_statistics + end subroutine statistics_score !=============================================================================== ! TALLY_STATISTICS computes the mean and standard deviation of the mean of each @@ -1824,11 +1824,11 @@ contains do i = 1, n_tallies t => tallies(i) - call calculate_statistics(t % scores) + call statistics_score(t % scores) end do ! Calculate statistics for global tallies - call calculate_statistics(global_tallies) + call statistics_score(global_tallies) end subroutine tally_statistics From 821c672cc7767abc4422ba42d6d976e4a4fbfa99 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 1 Apr 2012 12:41:46 -0400 Subject: [PATCH 06/10] Added reset_score subroutine. Closes #79 on github. --- src/tally.F90 | 134 ++++++++++++++++++++++++++++---------------------- 1 file changed, 75 insertions(+), 59 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 2c5eb44d8a..4ff68ff153 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1295,46 +1295,6 @@ contains end function get_next_bin -!=============================================================================== -! ADD_TO_SCORE accumulates a scoring contribution to a specific tally bin and -! specific response function. Note that we don't need to add the square of the -! contribution since that is done at the cycle level, not the history level -!=============================================================================== - - subroutine add_to_score(score, val) - - type(TallyScore), intent(inout) :: score - real(8), intent(in) :: val - - score % n_events = score % n_events + 1 - score % value = score % value + val - - end subroutine add_to_score - -!=============================================================================== -! ACCUMULATE_SCORE accumulates scores from many histories (or many generations) -! into a single realization of a random variable. -!=============================================================================== - - elemental subroutine accumulate_score(score) - - type(TallyScore), intent(inout) :: score - - real(8) :: val - - ! Add the sum and square of the sum of contributions from each cycle - ! within a cycle to the variables sum and sum_sq. This will later allow us - ! to calculate a variance on the tallies - - val = score % value/(n_particles*gen_per_batch) - score % sum = score % sum + val - score % sum_sq = score % sum_sq + val*val - - ! Reset the single batch estimate - score % value = ZERO - - end subroutine accumulate_score - !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable @@ -1790,25 +1750,6 @@ contains end function get_label -!=============================================================================== -! CALCULATE_STATISTICS determines the sample mean and the standard deviation of -! the mean for a TallyScore. -!=============================================================================== - - elemental subroutine statistics_score(score) - - type(TallyScore), intent(inout) :: score - - ! Calculate sample mean and standard deviation of the mean -- note that we - ! have used Bessel's correction so that the estimator of the variance of the - ! sample mean is unbiased. - - score % sum = score % sum/n_active - score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / & - (n_active - 1)) - - end subroutine statistics_score - !=============================================================================== ! TALLY_STATISTICS computes the mean and standard deviation of the mean of each ! tally and stores them in the val and val_sq attributes of the TallyScores @@ -1832,4 +1773,79 @@ contains end subroutine tally_statistics +!=============================================================================== +! ADD_TO_SCORE accumulates a scoring contribution to a specific tally bin and +! specific response function. Note that we don't need to add the square of the +! contribution since that is done at the cycle level, not the history level +!=============================================================================== + + subroutine add_to_score(score, val) + + type(TallyScore), intent(inout) :: score + real(8), intent(in) :: val + + score % n_events = score % n_events + 1 + score % value = score % value + val + + end subroutine add_to_score + +!=============================================================================== +! ACCUMULATE_SCORE accumulates scores from many histories (or many generations) +! into a single realization of a random variable. +!=============================================================================== + + elemental subroutine accumulate_score(score) + + type(TallyScore), intent(inout) :: score + + real(8) :: val + + ! Add the sum and square of the sum of contributions from each cycle + ! within a cycle to the variables sum and sum_sq. This will later allow us + ! to calculate a variance on the tallies + + val = score % value/(n_particles*gen_per_batch) + score % sum = score % sum + val + score % sum_sq = score % sum_sq + val*val + + ! Reset the single batch estimate + score % value = ZERO + + end subroutine accumulate_score + +!=============================================================================== +! STATISTICS_SCORE determines the sample mean and the standard deviation of the +! mean for a TallyScore. +!=============================================================================== + + elemental subroutine statistics_score(score) + + type(TallyScore), intent(inout) :: score + + ! Calculate sample mean and standard deviation of the mean -- note that we + ! have used Bessel's correction so that the estimator of the variance of the + ! sample mean is unbiased. + + score % sum = score % sum/n_active + score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / & + (n_active - 1)) + + end subroutine statistics_score + +!=============================================================================== +! RESET_SCORE zeroes out the value and accumulated sum and sum-squared for a +! single TallyScore. +!=============================================================================== + + elemental subroutine reset_score(score) + + type(TallyScore), intent(inout) :: score + + score % n_events = 0 + score % value = ZERO + score % sum = ZERO + score % sum_sq = ZERO + + end subroutine reset_score + end module tally From bca4121d02222507994205f34de63138d14dfac8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 1 Apr 2012 14:24:50 -0400 Subject: [PATCH 07/10] Updated command line options. --- src/initialize.F90 | 18 ++++++++++++++---- src/output.F90 | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 2f50cf93dd..677e7ee969 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -14,8 +14,8 @@ module initialize use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict, read_plots_xml use output, only: title, header, print_summary, print_geometry, & - print_plot, create_summary_file, & - create_xs_summary_file + print_plot, create_summary_file, print_usage, & + create_xs_summary_file, print_version use random_lcg, only: initialize_prng use source, only: initialize_source use string, only: to_str, starts_with, ends_with, lower_case @@ -236,9 +236,19 @@ contains ! Check for flags if (starts_with(argv(i), "-")) then - if (argv(i)(2:5) == 'plot') then + select case (argv(i)) + case ('-p', '-plot', '--plot') plotting = .true. - end if + case ('-?', '-help', '--help') + call print_usage() + stop + case ('-v', '-version', '--version') + call print_version() + stop + case default + message = "Unknown command line option: " // argv(i) + call fatal_error() + end select last_flag = i end if diff --git a/src/output.F90 b/src/output.F90 index 2fc6c77e96..d486be1c7f 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -138,6 +138,41 @@ contains end subroutine header +!=============================================================================== +! PRINT_VERSION shows the current version as well as copright and license +! information +!=============================================================================== + + subroutine print_version() + + if (master) then + write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') & + "OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE + write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2012 & + &Massachusetts Institute of Technology" + write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at & + &" + end if + + end subroutine print_version + +!=============================================================================== +! PRINT_USAGE displays information about command line usage of OpenMC +!=============================================================================== + + subroutine print_usage() + + if (master) then + write(OUTPUT_UNIT,*) 'Usage: openmc [options] [directory]' + write(OUTPUT_UNIT,*) + write(OUTPUT_UNIT,*) 'Options:' + write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode' + write(OUTPUT_UNIT,*) ' -v, --version Show version information' + write(OUTPUT_UNIT,*) ' -?, --help Show this message' + end if + + end subroutine print_usage + !=============================================================================== ! WRITE_MESSAGE displays an informational message to the log file and the ! standard output stream. From 08bda29d3727340c5f77f610dcd6da9f066720d4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 1 Apr 2012 14:37:54 -0400 Subject: [PATCH 08/10] More updates to documentation. --- docs/source/index.rst | 1 + docs/source/license.rst | 24 +++++++++++++ docs/source/usersguide/input.rst | 25 +++++++++++++ docs/source/usersguide/setup.rst | 62 +++++++++++++++++++++++++++----- 4 files changed, 104 insertions(+), 8 deletions(-) create mode 100644 docs/source/license.rst diff --git a/docs/source/index.rst b/docs/source/index.rst index 5b8bebff35..eb9775b77c 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -29,4 +29,5 @@ Contents usersguide/index devguide/index publications + license developers diff --git a/docs/source/license.rst b/docs/source/license.rst new file mode 100644 index 0000000000..e8ebbaad1e --- /dev/null +++ b/docs/source/license.rst @@ -0,0 +1,24 @@ +.. _license: + +================= +License Agreement +================= + +Copyright © 2011-2012 Massachusetts Institute of Technology + +Permission is hereby granted, free of charge, to any person obtaining a copy of +this software and associated documentation files (the "Software"), to deal in +the Software without restriction, including without limitation the rights to +use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software is furnished to do so, +subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS +FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR +COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 32f63bbcfa..f6104bba4c 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -210,6 +210,31 @@ integers: the batch number, generation number, and particle number. *Default*: None +```` Element +------------------------ + +The ```` element describes a mesh that is used for re-weighting +source sites at every generation based on the uniform fission site methodology +described in Kelly et al., "MC21 Analysis of the Nuclear Energy Agency Monte +Carlo Performance Benchmark Problem," Proceedings of *Physor 2012*, Knoxville, +TN (2012). This mesh should cover all possible fissionable materials in the +problem. It has the following attributes/sub-elements: + + :dimension: + The number of mesh cells in the x, y, and z directions, respectively. + + *Default*: None + + :lower_left: + The Cartersian coordinates of the lower-left corner of the mesh. + + *Default*: None + + :upper_right: + The Cartersian coordinates of the upper-right corner of the mesh. + + *Default*: None + ```` Element ----------------------- diff --git a/docs/source/usersguide/setup.rst b/docs/source/usersguide/setup.rst index a148000c77..b49eb08d42 100644 --- a/docs/source/usersguide/setup.rst +++ b/docs/source/usersguide/setup.rst @@ -41,9 +41,11 @@ been compiled with the same compiler you intend to compile OpenMC with. Obtaining the Source -------------------- -All OpenMC source code is hosted on GitHub_. This means that you will need to -have git_ installed on your computer in order to get source code and updates -directly from the repository. GitHub has a good set of `instructions +All OpenMC source code is hosted on GitHub_. You can download the source code +directly from GitHub or, if you have the git_ version control software installed +on your computer, you can use git to obtain the source code. The latter method +has the benefit that it is easy to receive updates directly from the GitHub +repository. GitHub has a good set of `instructions `_ for how to set up git to work with GitHub since this involves setting up ssh_ keys. With git installed and setup, the following command will download the full source code from the GitHub @@ -96,12 +98,12 @@ Makefile, you can enter the following from a terminal:: make DEBUG=yes ---------- -Compiling ---------- +------------------------------- +Compiling on Linux and Mac OS X +------------------------------- -To compile the code, run the following commands from within the root directory -for OpenMC: +To compile OpenMC on Linux or Max OS X, run the following commands from within +the root directory of the source code: .. code-block:: sh @@ -110,6 +112,25 @@ for OpenMC: This will build an executable named ``openmc``. +-------------------- +Compiling on Windows +-------------------- + +To compile OpenMC on a Windows operating system, you will need to first install +Cygwin_, a Linux-like environment for Windows. When configuring Cygwin, make +sure you install both the gfortran compiler as well as git. Once you have +obtained the source code, run the following commands from within the source code +root directory: + +.. code-block:: sh + + cd src + make + +This will build an executable named ``openmc``. + +.. _Cygwin: http://www.cygwin.com/ + --------------------------- Cross-Section Configuration --------------------------- @@ -130,3 +151,28 @@ directory of the OpenMC distribution to the location of the Serpent cross-sections. Then, either set the :ref:`cross_sections` in a settings.xml file or the :envvar:`CROSS_SECTIONS` environment variable to the absolute path of the ``cross_sections_serpent.xml`` file. + +-------------- +Running OpenMC +-------------- + +Once you have a model built (see :ref:`usersguide_input`), you can either run +the openmc executable directly from the directory containing your XML input +files, or you can specify as a command-line argument the directory containing +the XML input files. For example, if the path of your OpenMC exectuable is +``/home/username/openmc/src/openmc`` and your XML input files are in the +directory ``/home/username/somemodel/``, one way to run the simulation would be: + +.. code-block:: sh + + cd /home/username/somemodel + openmc + +Alternatively, you could run from any directory: + +.. code-block:: sh + + openmc /home/username/somemodel + +Note that in the latter case, any output files will be placed in the present +working directory which may be different from ``/home/username/somemodel``. From 671f305b3e937718d067aea180aa27f5c05f5a5e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 1 Apr 2012 14:55:25 -0400 Subject: [PATCH 09/10] Added mesh_intersects check in score_surface_current. Closes #58 on github. --- src/tally.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 4ff68ff153..2169af5c14 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1014,8 +1014,11 @@ contains call get_mesh_indices(m, xyz0, ijk0, start_in_mesh) call get_mesh_indices(m, xyz1, ijk1, end_in_mesh) - ! Check to make sure start or end is in mesh - if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) cycle + ! Check to if start or end is in mesh -- if not, check if track still + ! intersects with mesh + if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then + if (.not. mesh_intersects(m, xyz0, xyz1)) cycle + end if ! Calculate number of surface crossings n_cross = sum(abs(ijk1 - ijk0)) From 2f3bbef81cd7199f0f345f1f16f3ac15891f776e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 1 Apr 2012 20:25:26 -0400 Subject: [PATCH 10/10] Fixed subtle bug with S(a,b) cross section calculation. --- src/cross_section.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index aa155eefe7..be8f009a7b 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -187,8 +187,13 @@ contains end if end if - ! Set last evaluated energy - if (index_sab == 0) micro_xs(index_nuclide) % last_E = p % E + ! Set last evaluated energy -- if we're in S(a,b) region, force + ! re-calculation of cross-section + if (index_sab == 0) then + micro_xs(index_nuclide) % last_E = p % E + else + micro_xs(index_nuclide) % last_E = ZERO + end if end subroutine calculate_nuclide_xs