From 569a532af61cc0b2b27306d4e9e73ed04a1ba45c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 21 Aug 2017 11:16:56 -0500 Subject: [PATCH] Store entropy and UFS meshes in global meshes array --- src/eigenvalue.F90 | 43 +++----------- src/global.F90 | 17 ------ src/input_xml.F90 | 137 +++++++++++++++++++++++---------------------- src/mesh.F90 | 2 +- src/physics.F90 | 24 ++++---- src/physics_mg.F90 | 25 +++++---- 6 files changed, 107 insertions(+), 141 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index accf80b310..a1bb356559 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -300,39 +300,9 @@ contains integer :: ent_idx ! entropy index integer :: i, j, k ! index for bank sites - integer :: n ! # of boxes in each dimension logical :: sites_outside ! were there sites outside entropy box? - type(RegularMesh), pointer :: m - ! Get pointer to entropy mesh - m => entropy_mesh - - ! 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 - if (.not. allocated(m % dimension)) then - ! If the user did not specify how many mesh cells are to be used in - ! each direction, we automatically determine an appropriate number of - ! cells - n = ceiling((n_particles/20)**(ONE/THREE)) - - ! copy dimensions - m % n_dimension = 3 - allocate(m % dimension(3)) - m % dimension = n - - ! determine width - m % width = (m % upper_right - m % lower_left) / m % dimension - - end if - - ! allocate p - allocate(entropy_p(1, m % dimension(1), m % dimension(2), & - m % dimension(3))) - end if + associate (m => meshes(index_entropy_mesh)) ! count number of fission sites over mesh call count_bank_sites(m, fission_bank, entropy_p, & @@ -362,6 +332,7 @@ contains end do end if + end associate end subroutine shannon_entropy !=============================================================================== @@ -621,16 +592,18 @@ contains integer :: mpi_err ! MPI error code #endif + associate (m => meshes(index_ufs_mesh)) + if (current_batch == 1 .and. current_gen == 1) then ! On the first generation, 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 + source_frac = m % volume_frac else ! count number of source sites in each ufs mesh cell - call count_bank_sites(ufs_mesh, source_bank, source_frac, & + call count_bank_sites(m, source_bank, source_frac, & sites_outside=sites_outside, size_bank=work) ! Check for sites outside of the mesh @@ -640,7 +613,7 @@ contains #ifdef MPI ! Send source fraction to all processors - n = product(ufs_mesh % dimension) + n = product(m % dimension) call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) #endif @@ -654,6 +627,8 @@ contains source_bank % wgt = source_bank % wgt * n_particles / total end if + end associate + end subroutine count_source_for_ufs #ifdef _OPENMP diff --git a/src/global.F90 b/src/global.F90 index 5940b6a019..dea17a102e 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -202,25 +202,8 @@ contains call statepoint_batch % clear() call sourcepoint_batch % clear() - ! Deallocate entropy mesh - if (associated(entropy_mesh)) then - if (allocated(entropy_mesh % lower_left)) & - deallocate(entropy_mesh % lower_left) - if (allocated(entropy_mesh % upper_right)) & - deallocate(entropy_mesh % upper_right) - if (allocated(entropy_mesh % width)) deallocate(entropy_mesh % width) - deallocate(entropy_mesh) - end if - ! Deallocate ufs if (allocated(source_frac)) deallocate(source_frac) - if (associated(ufs_mesh)) then - if (allocated(ufs_mesh % lower_left)) deallocate(ufs_mesh % lower_left) - if (allocated(ufs_mesh % upper_right)) & - deallocate(ufs_mesh % upper_right) - if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width) - deallocate(ufs_mesh) - end if end subroutine free_memory diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 0cf33c0240..e650f0fdf1 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -144,6 +144,7 @@ contains integer :: n integer :: temp_int integer :: temp_int_array3(3) + integer(C_INT) :: err integer, allocatable :: temp_int_array(:) real(8), allocatable :: temp_real(:) integer :: n_tracks @@ -700,47 +701,54 @@ contains &corner of Shannon entropy mesh.") end if - ! Allocate mesh object and coordinates on mesh - allocate(entropy_mesh) - allocate(entropy_mesh % lower_left(3)) - allocate(entropy_mesh % upper_right(3)) - allocate(entropy_mesh % width(3)) + err = openmc_extend_meshes(1, index_entropy_mesh) - ! Copy values - call get_node_array(node_entropy, "lower_left", & - entropy_mesh % lower_left) - call get_node_array(node_entropy, "upper_right", & - entropy_mesh % upper_right) + associate (m => meshes(index_entropy_mesh)) + ! Allocate mesh object and coordinates on mesh + allocate(m % lower_left(3)) + allocate(m % upper_right(3)) + allocate(m % width(3)) - ! Check on values provided - if (.not. all(entropy_mesh % upper_right > entropy_mesh % lower_left)) & - &then - call fatal_error("Upper-right coordinate must be greater than & - &lower-left coordinate for Shannon entropy mesh.") - end if + ! Copy values + call get_node_array(node_entropy, "lower_left", m % lower_left) + call get_node_array(node_entropy, "upper_right", m % upper_right) - ! Check if dimensions were specified -- if not, they will be calculated - ! automatically upon first entry into shannon_entropy - if (check_for_node(node_entropy, "dimension")) then - - ! If so, make sure proper number of values were given - if (node_word_count(node_entropy, "dimension") /= 3) then - call fatal_error("Dimension of entropy mesh must be given as three & - &integers.") + ! Check on values provided + if (.not. all(m % upper_right > m % lower_left)) & + &then + call fatal_error("Upper-right coordinate must be greater than & + &lower-left coordinate for Shannon entropy mesh.") end if ! Allocate dimensions - entropy_mesh % n_dimension = 3 - allocate(entropy_mesh % dimension(3)) + m % n_dimension = 3 + allocate(m % dimension(3)) - ! Copy dimensions - call get_node_array(node_entropy, "dimension", entropy_mesh % dimension) + ! Check if dimensions were specified -- if not, they will be calculated + ! automatically upon first entry into shannon_entropy + if (check_for_node(node_entropy, "dimension")) then + ! If so, make sure proper number of values were given + if (node_word_count(node_entropy, "dimension") /= 3) then + call fatal_error("Dimension of entropy mesh must be given as three & + &integers.") + end if + + ! Copy dimensions + call get_node_array(node_entropy, "dimension", m % dimension) + else + ! If the user did not specify how many mesh cells are to be used in + ! each direction, we automatically determine an appropriate number of + ! cells + m % dimension = ceiling((n_particles/20)**(ONE/THREE)) + end if ! Calculate width - entropy_mesh % width = (entropy_mesh % upper_right - & - entropy_mesh % lower_left) / entropy_mesh % dimension + m % width = (m % upper_right - m % lower_left) / m % dimension - end if + ! Allocate space for storing number of fission sites in each mesh cell + allocate(entropy_p(1, m % dimension(1), m % dimension(2), & + m % dimension(3))) + end associate ! Turn on Shannon entropy calculation entropy_on = .true. @@ -764,42 +772,44 @@ contains &integers.") end if + err = openmc_extend_meshes(1, index_ufs_mesh) + ! 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)) + associate (m => meshes(index_ufs_mesh)) + allocate(m % lower_left(3)) + allocate(m % upper_right(3)) + allocate(m % width(3)) - ! Allocate dimensions - ufs_mesh % n_dimension = 3 - allocate(ufs_mesh % dimension(3)) + ! Allocate dimensions + m % n_dimension = 3 + allocate(m % dimension(3)) - ! Copy dimensions - call get_node_array(node_ufs, "dimension", ufs_mesh % dimension) + ! Copy dimensions + call get_node_array(node_ufs, "dimension", m % dimension) - ! Copy values - call get_node_array(node_ufs, "lower_left", ufs_mesh % lower_left) - call get_node_array(node_ufs, "upper_right", ufs_mesh % upper_right) + ! Copy values + call get_node_array(node_ufs, "lower_left", m % lower_left) + call get_node_array(node_ufs, "upper_right", m % upper_right) - ! Check on values provided - if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then - call fatal_error("Upper-right coordinate must be greater than & - &lower-left coordinate for UFS mesh.") - end if + ! Check on values provided + if (.not. all(m % upper_right > m % lower_left)) then + call fatal_error("Upper-right coordinate must be greater than & + &lower-left coordinate for UFS mesh.") + end if - ! Calculate width - ufs_mesh % width = (ufs_mesh % upper_right - & - ufs_mesh % lower_left) / ufs_mesh % dimension + ! Calculate width + m % width = (m % upper_right - m % lower_left) / m % dimension - ! Calculate volume fraction of each cell - ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8) + ! Calculate volume fraction of each cell + m % volume_frac = ONE/real(product(m % dimension),8) + + ! Allocate source_frac + allocate(source_frac(1, m % dimension(1), m % dimension(2), & + m % dimension(3))) + end associate ! 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 state points @@ -4566,12 +4576,12 @@ contains select case (trim(meshtype)) case ('ufs') - if (.not. associated(ufs_mesh)) then + if (index_ufs_mesh < 0) then call fatal_error("No UFS mesh for meshlines on plot " & // trim(to_str(pl % id))) end if - pl % meshlines_mesh => ufs_mesh + pl % meshlines_mesh => meshes(index_ufs_mesh) case ('cmfd') @@ -4584,17 +4594,12 @@ contains case ('entropy') - if (.not. associated(entropy_mesh)) then + if (index_entropy_mesh < 0) then call fatal_error("No entropy mesh for meshlines on plot " & // trim(to_str(pl % id))) end if - if (.not. allocated(entropy_mesh % dimension)) then - call fatal_error("No dimension specified on entropy mesh & - &for meshlines on plot " // trim(to_str(pl % id))) - end if - - pl % meshlines_mesh => entropy_mesh + pl % meshlines_mesh => meshes(index_entropy_mesh) case ('tally') diff --git a/src/mesh.F90 b/src/mesh.F90 index 6d94ffba48..9b8943ff25 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -20,7 +20,7 @@ contains subroutine count_bank_sites(m, bank_array, cnt, energies, size_bank, & sites_outside) - type(RegularMesh), pointer :: m ! mesh to count sites + type(RegularMesh), intent(in) :: 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 diff --git a/src/physics.F90 b/src/physics.F90 index 8081aebcbd..774106ddf7 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1088,18 +1088,20 @@ contains ! the expected number of fission sites produced if (ufs) then - ! Determine indices on ufs mesh for current location - call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh) - if (.not. in_mesh) then - call write_particle_restart(p) - call fatal_error("Source site outside UFS mesh!") - end if + associate (m => meshes(index_ufs_mesh)) + ! Determine indices on ufs mesh for current location + call m % get_indices(p % coord(1) % xyz, ijk, in_mesh) + if (.not. in_mesh) then + call write_particle_restart(p) + call fatal_error("Source site outside UFS mesh!") + 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 + if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then + weight = m % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3)) + else + weight = ONE + end if + end associate else weight = ONE end if diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index c412e9e32a..7f1ead66a3 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -193,20 +193,21 @@ contains ! the expected number of fission sites produced if (ufs) then + associate (m => meshes(index_ufs_mesh)) + ! Determine indices on ufs mesh for current location + call m % get_indices(p % coord(1) % xyz, ijk, in_mesh) - ! Determine indices on ufs mesh for current location - call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh) + if (.not. in_mesh) then + call write_particle_restart(p) + call fatal_error("Source site outside UFS mesh!") + end if - if (.not. in_mesh) then - call write_particle_restart(p) - call fatal_error("Source site outside UFS mesh!") - 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 + if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then + weight = m % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3)) + else + weight = ONE + end if + end associate else weight = ONE end if