From 7819bef5bc514fed390404c07c44631cd5354bb0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 31 Mar 2012 17:22:03 -0400 Subject: [PATCH] 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 @@ - + +