diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index cd52068fa0..155b0eaf54 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -107,6 +107,14 @@ class TallyView(object): _dll.openmc_tally_results(self._index, data, shape) return as_array(data, tuple(shape[::-1])) + @filters.setter + def filters(self, filters): + # Get filter indices as int32_t[] + n = len(filters) + indices = (c_int32*n)(*(f._index for f in filters)) + + _dll.openmc_tally_set_filters(self._index, n, indices) + @nuclides.setter def nuclides(self, nuclides): nucs = (c_char_p * len(nuclides))() diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 67341ae953..769abf1236 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -261,14 +261,16 @@ contains logical :: energy_filters integer :: i ! loop counter integer :: n ! size of arrays in mesh specification - integer :: ng ! number of energy groups (default 1) + integer(C_INT32_T) :: ng ! number of energy groups (default 1) integer :: n_filter ! number of filters integer :: i_start, i_end integer :: i_filt_start, i_filt_end + integer(C_INT32_T), allocatable :: filter_indices(:) integer(C_INT) :: err integer :: i_filt ! index in filters array integer :: iarray3(3) ! temp integer array real(8) :: rarray3(3) ! temp double array + real(C_DOUBLE), allocatable :: energies(:) type(TallyObject), pointer :: t type(RegularMesh), pointer :: m type(XMLNode) :: node_mesh @@ -385,60 +387,39 @@ contains ! Set up mesh filter i_filt = i_filt_start - allocate(MeshFilter :: filters(i_filt) % obj) - select type (filt => filters(i_filt) % obj) - type is (MeshFilter) - filt % id = i_filt - filt % n_bins = product(m % dimension) - filt % mesh = i_start - ! Add filter to dictionary - call filter_dict % add_key(filt % id, i_filt) - end select + err = openmc_filter_set_type(i_filt, C_CHAR_'mesh' // C_NULL_CHAR) + err = openmc_filter_set_id(i_filt, i_filt) + err = openmc_mesh_filter_set_mesh(i_filt, i_start) if (energy_filters) then ! Read and set incoming energy mesh filter i_filt = i_filt + 1 - allocate(EnergyFilter :: filters(i_filt) % obj) - select type (filt => filters(i_filt) % obj) - type is (EnergyFilter) - filt % id = i_filt - ng = node_word_count(node_mesh, "energy") - filt % n_bins = ng - 1 - allocate(filt % bins(ng)) - call get_node_array(node_mesh, "energy", filt % bins) - ! Add filter to dictionary - call filter_dict % add_key(filt % id, i_filt) - end select + err = openmc_filter_set_type(i_filt, C_CHAR_'energy' // C_NULL_CHAR) + err = openmc_filter_set_id(i_filt, i_filt) + + ! Get energies and set bins + ng = node_word_count(node_mesh, "energy") + allocate(energies(ng)) + call get_node_array(node_mesh, "energy", energies) + err = openmc_energy_filter_set_bins(i_filt, ng, energies) ! Read and set outgoing energy mesh filter i_filt = i_filt + 1 - allocate(EnergyoutFilter :: filters(i_filt) % obj) - select type (filt => filters(i_filt) % obj) - type is (EnergyoutFilter) - filt % id = i_filt - ng = node_word_count(node_mesh, "energy") - filt % n_bins = ng - 1 - allocate(filt % bins(ng)) - call get_node_array(node_mesh, "energy", filt % bins) - ! Add filter to dictionary - call filter_dict % add_key(filt % id, i_filt) - end select + err = openmc_filter_set_type(i_filt, C_CHAR_'energyout' // C_NULL_CHAR) + err = openmc_filter_set_id(i_filt, i_filt) + err = openmc_energy_filter_set_bins(i_filt, ng, energies) end if ! Duplicate the mesh filter for the mesh current tally since other ! tallies use this filter and we need to change the dimension i_filt = i_filt + 1 - allocate(MeshFilter :: filters(i_filt) % obj) - select type (filt => filters(i_filt) % obj) - type is (MeshFilter) - filt % id = i_filt - ! We need to increase the dimension by one since we also need - ! currents coming into and out of the boundary mesh cells. - filt % n_bins = product(m % dimension + 1) - filt % mesh = i_start - ! Add filter to dictionary - call filter_dict % add_key(filt % id, i_filt) - end select + err = openmc_filter_set_type(i_filt, C_CHAR_'mesh' // C_NULL_CHAR) + err = openmc_filter_set_id(i_filt, i_filt) + err = openmc_mesh_filter_set_mesh(i_filt, i_start) + + ! We need to increase the dimension by one since we also need + ! currents coming into and out of the boundary mesh cells. + filters(i_filt) % obj % n_bins = product(m % dimension + 1) ! Set up surface filter i_filt = i_filt + 1 @@ -476,15 +457,11 @@ contains call get_node_value(root, "reset", t % reset) end if - ! Set the mesh filter index in the tally find_filter array - n_filter = 1 - t % find_filter(FILTER_MESH) = n_filter - ! Set the incoming energy mesh filter index in the tally find_filter ! array + n_filter = 1 if (energy_filters) then n_filter = n_filter + 1 - t % find_filter(FILTER_ENERGYIN) = n_filter end if ! Set number of nucilde bins @@ -507,11 +484,13 @@ contains t % type = TALLY_VOLUME ! Allocate and set filters - allocate(t % filter(n_filter)) - t % filter(1) = i_filt_start + allocate(filter_indices(n_filter)) + filter_indices(1) = i_filt_start if (energy_filters) then - t % filter(2) = i_filt_start + 1 + filter_indices(2) = i_filt_start + 1 end if + err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices) + deallocate(filter_indices) ! Allocate scoring bins allocate(t % score_bins(3)) @@ -543,16 +522,17 @@ contains ! array if (energy_filters) then n_filter = n_filter + 1 - t % find_filter(FILTER_ENERGYOUT) = n_filter end if ! Allocate and set indices in filters array - allocate(t % filter(n_filter)) - t % filter(1) = i_filt_start + allocate(filter_indices(n_filter)) + filter_indices(1) = i_filt_start if (energy_filters) then - t % filter(2) = i_filt_start + 1 - t % filter(3) = i_filt_start + 2 + filter_indices(2) = i_filt_start + 1 + filter_indices(3) = i_filt_start + 2 end if + err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices) + deallocate(filter_indices) ! Allocate macro reactions allocate(t % score_bins(2)) @@ -577,15 +557,16 @@ contains ! Set the surface filter index in the tally find_filter array n_filter = n_filter + 1 - t % find_filter(FILTER_SURFACE) = n_filter ! Allocate and set filters - allocate(t % filter(n_filter)) - t % filter(1) = i_filt_end - 1 - t % filter(n_filter) = i_filt_end + allocate(filter_indices(n_filter)) + filter_indices(1) = i_filt_end - 1 + filter_indices(n_filter) = i_filt_end if (energy_filters) then - t % filter(2) = i_filt_start + 1 + filter_indices(2) = i_filt_start + 1 end if + err = openmc_tally_set_filters(i_start + i - 1, n_filter, filter_indices) + deallocate(filter_indices) ! Allocate macro reactions allocate(t % score_bins(1)) diff --git a/src/initialize.F90 b/src/initialize.F90 index 92273a6ca5..42db8c7774 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -29,7 +29,7 @@ module initialize use state_point, only: load_state_point use string, only: to_str, starts_with, ends_with, str_to_int use summary, only: write_summary - use tally_header, only: TallyObject, configure_tallies + use tally_header, only: TallyObject use tally_filter use tally, only: init_tally_routines @@ -116,9 +116,6 @@ contains end if if (run_mode /= MODE_PLOTTING) then - ! Allocate and setup tally stride, filter_matches, and tally maps - call configure_tallies() - ! Set up tally procedure pointers call init_tally_routines() diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3225c51b78..15ee43f1b9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3315,11 +3315,7 @@ contains end select ! Set filter id - f % obj % id = filter_id - - ! Add filter to dictionary - call filter_dict % add_key(filter_id, i) - + err = openmc_filter_set_id(i_start + i - 1, filter_id) end do READ_FILTERS ! ========================================================================== @@ -3394,64 +3390,28 @@ contains allocate(temp_filter(n_filter)) if (n_filter > 0) then call get_node_array(node_tal, "filters", temp_filter) + + do j = 1, n_filter + ! Get pointer to filter + if (filter_dict % has_key(temp_filter(j))) then + i_filt = filter_dict % get_key(temp_filter(j)) + f => filters(i_filt) + else + call fatal_error("Could not find filter " & + // trim(to_str(temp_filter(j))) // " specified on tally " & + // trim(to_str(t % id))) + end if + + ! Store the index of the filter + temp_filter(j) = i_filt + end do + + ! Set the filters + err = openmc_tally_set_filters(i_start + i - 1, n_filter, temp_filter) + else + allocate(t % filter(n_filter)) end if - do j = 1, n_filter - ! Get pointer to filter - if (filter_dict % has_key(temp_filter(j))) then - i_filt = filter_dict % get_key(temp_filter(j)) - f => filters(i_filt) - else - call fatal_error("Could not find filter " & - // trim(to_str(temp_filter(j))) // " specified on tally " & - // trim(to_str(t % id))) - end if - - ! Set the filter index in the tally find_filter array - select type (filt => f % obj) - type is (DistribcellFilter) - t % find_filter(FILTER_DISTRIBCELL) = j - type is (CellFilter) - t % find_filter(FILTER_CELL) = j - type is (CellFromFilter) - t % find_filter(FILTER_CELLFROM) = j - type is (CellbornFilter) - t % find_filter(FILTER_CELLBORN) = j - type is (MaterialFilter) - t % find_filter(FILTER_MATERIAL) = j - type is (UniverseFilter) - t % find_filter(FILTER_UNIVERSE) = j - type is (SurfaceFilter) - t % find_filter(FILTER_SURFACE) = j - type is (MeshFilter) - t % find_filter(FILTER_MESH) = j - type is (EnergyFilter) - t % find_filter(FILTER_ENERGYIN) = j - type is (EnergyoutFilter) - t % find_filter(FILTER_ENERGYOUT) = j - ! Set to analog estimator - t % estimator = ESTIMATOR_ANALOG - type is (DelayedGroupFilter) - t % find_filter(FILTER_DELAYEDGROUP) = j - type is (MuFilter) - t % find_filter(FILTER_MU) = j - ! Set to analog estimator - t % estimator = ESTIMATOR_ANALOG - type is (PolarFilter) - t % find_filter(FILTER_POLAR) = j - type is (AzimuthalFilter) - t % find_filter(FILTER_AZIMUTHAL) = j - type is (EnergyFunctionFilter) - t % find_filter(FILTER_ENERGYFUNCTION) = j - end select - - ! Store the index of the filter - temp_filter(j) = i_filt - end do - - ! Store the filter indices - call move_alloc(FROM=temp_filter, TO=t % filter) - ! ======================================================================= ! READ DATA FOR NUCLIDES diff --git a/src/simulation.F90 b/src/simulation.F90 index 2a063ddeb4..8e977bbdc3 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -21,6 +21,7 @@ module simulation use state_point, only: write_state_point, write_source_point use string, only: to_str use tally, only: accumulate_tallies, setup_active_tallies + use tally_header, only: configure_tallies use trigger, only: check_triggers use tracking, only: transport @@ -371,7 +372,11 @@ contains subroutine initialize_simulation() + ! Allocate tally results arrays if they're not allocated yet + call configure_tallies() + !$omp parallel + ! Allocate array for microscopic cross section cache allocate(micro_xs(n_nuclides)) ! Allocate array for matching filter bins diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index a76d08bf4d..32baad3e90 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -182,6 +182,8 @@ contains if (index >= 1 .and. index <= n_filters) then if (allocated(filters(index) % obj)) then filters(index) % obj % id = id + call filter_dict % add_key(id, index) + err = 0 else err = E_FILTER_NOT_ALLOCATED diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index edfa66f3b9..e9654f6094 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -108,7 +108,7 @@ module tally_header contains procedure :: read_results_hdf5 => tally_read_results_hdf5 procedure :: write_results_hdf5 => tally_write_results_hdf5 - procedure :: setup_arrays => tally_setup_arrays + procedure :: allocate_results => tally_allocate_results end type TallyObject integer(C_INT32_T), public, bind(C) :: n_tallies = 0 ! # of tallies @@ -210,45 +210,32 @@ contains end subroutine tally_read_results_hdf5 !=============================================================================== -! SETUP_ARRAYS allocates and populates several member arrays of the TallyObject -! derived type, including stride, filter_matches, and results. +! ALLOCATE_RESULTS allocates and initializes the results component of the +! TallyObject derived type !=============================================================================== - subroutine tally_setup_arrays(this) + subroutine tally_allocate_results(this) class(TallyObject), intent(inout) :: this - integer :: j ! loop index for filters - integer :: n ! temporary stride - integer :: i_filt ! filter index - - ! Allocate stride - if (allocated(this % filter)) then - if (allocated(this % stride)) deallocate(this % stride) - allocate(this % stride(size(this % filter))) - - ! The filters are traversed in opposite order so that the last filter has - ! the shortest stride in memory and the first filter has the largest - ! stride - n = 1 - STRIDE: do j = size(this % filter), 1, -1 - i_filt = this % filter(j) - this % stride(j) = n - n = n * filters(i_filt) % obj % n_bins - end do STRIDE - this % n_filter_bins = n - else - this % n_filter_bins = 1 - end if - ! Set total number of filter and scoring bins this % total_score_bins = this % n_score_bins * this % n_nuclide_bins - ! Allocate results array - if (allocated(this % results)) deallocate(this % results) - allocate(this % results(3, this % total_score_bins, this % n_filter_bins)) + if (allocated(this % results)) then + ! If results was already allocated but shape is wrong, then reallocate it + ! to the correct shape + if (this % total_score_bins /= size(this % results, 2) .or. & + this % n_filter_bins /= size(this % results, 3)) then + deallocate(this % results) + allocate(this % results(3, this % total_score_bins, this % n_filter_bins)) + end if + else + allocate(this % results(3, this % total_score_bins, this % n_filter_bins)) + end if + + ! Initialize results array to zero this % results(:,:,:) = ZERO - end subroutine tally_setup_arrays + end subroutine tally_allocate_results !=============================================================================== ! CONFIGURE_TALLIES initializes several data structures related to tallies. This @@ -260,12 +247,14 @@ contains integer :: i - ! Allocate global tallies - allocate(global_tallies(3, N_GLOBAL_TALLIES)) + ! Allocate and initialize global tallies + if (.not. allocated(global_tallies)) then + allocate(global_tallies(3, N_GLOBAL_TALLIES)) + end if global_tallies(:,:) = ZERO do i = 1, n_tallies - call tallies(i) % setup_arrays() + call tallies(i) % allocate_results() end do end subroutine configure_tallies @@ -412,24 +401,25 @@ contains integer(C_INT32_T), intent(in) :: filter_indices(n) integer(C_INT) :: err - integer :: i - integer :: j + integer :: i ! index in t % filter/stride + integer :: j ! index in t % find_filter + integer :: k ! index in global filters array + integer :: stride ! filter stride err = 0 if (index >= 1 .and. index <= n_tallies) then associate (t => tallies(index)) - if (allocated(t % filter)) deallocate(t % filter) - allocate(t % filter(n)) t % find_filter(:) = 0 do i = 1, n - if (filter_indices(i) < 1 .or. filter_indices(i) > n_filters) then + k = filter_indices(i) + if (k < 1 .or. k > n_filters) then err = E_OUT_OF_BOUNDS exit end if ! Set the filter index in the tally find_filter array - select type (filt => filters(i) % obj) + select type (filt => filters(k) % obj) type is (DistribcellFilter) j = FILTER_DISTRIBCELL type is (CellFilter) @@ -465,7 +455,28 @@ contains end select t % find_filter(j) = i end do - if (err == 0) t % filter(:) = filter_indices + + if (err == 0) then + if (allocated(t % filter)) deallocate(t % filter) + if (allocated(t % stride)) deallocate(t % stride) + allocate(t % filter(n), t % stride(n)) + + ! Filters are traversed in reverse so that the last filter has the + ! shortest stride in memory and the first filter has the largest stride + stride = 1 + do i = n, 1, -1 + ! Set filter and stride + k = filter_indices(i) + t % filter(i) = k + t % stride(i) = stride + + ! Multiply stride by number of bins in this filter + stride = stride * filters(k) % obj % n_bins + end do + + ! Set total number of filter bins + t % n_filter_bins = stride + end if end associate else err = E_OUT_OF_BOUNDS @@ -509,8 +520,6 @@ contains end select end do - call t % setup_arrays() - err = 0 end associate else @@ -674,8 +683,6 @@ contains end select end do - call t % setup_arrays() - err = 0 end associate else