From 04e48224e6ce8272acd05080fded3326ad774159 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Aug 2018 14:22:55 -0500 Subject: [PATCH] Almost complete mesh conversion --- CMakeLists.txt | 1 + include/openmc/eigenvalue.h | 4 + include/openmc/hdf5_interface.h | 8 + include/openmc/particle.h | 2 +- src/cmfd_data.F90 | 7 +- src/cmfd_header.F90 | 3 +- src/cmfd_input.F90 | 110 +------- src/eigenvalue.F90 | 1 - src/eigenvalue.cpp | 5 + src/input_xml.F90 | 12 +- src/mesh.cpp | 6 +- src/mesh_header.F90 | 42 ++- src/particle.cpp | 2 +- src/physics.F90 | 1 - src/physics_mg.F90 | 1 - src/plot.F90 | 6 +- src/settings.cpp | 1 + src/simulation.F90 | 5 +- src/simulation_header.F90 | 10 +- src/state_point.F90 | 2 +- src/tallies/tally_filter_mesh.F90 | 277 ++++++++++---------- src/tallies/tally_filter_meshsurface.F90 | 309 ++++++++++++----------- src/tallies/trigger.F90 | 11 +- 23 files changed, 380 insertions(+), 446 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 78d1caf4d..713bf7c07 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -385,6 +385,7 @@ add_library(libopenmc SHARED src/distribution_energy.cpp src/distribution_multi.cpp src/distribution_spatial.cpp + src/eigenvalue.cpp src/endf.cpp src/initialize.cpp src/finalize.cpp diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index c2e76e26f..cc879461e 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -2,6 +2,9 @@ #define OPENMC_EIGENVALUE_H #include // for int64_t +#include + +#include "xtensor/xtensor.hpp" #include "openmc/particle.h" @@ -15,6 +18,7 @@ extern std::vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS extern "C" int64_t n_bank; +#pragma omp threadprivate(n_bank) //============================================================================== // Non-member functions diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 944fe75cb..b53703fc9 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -317,6 +317,14 @@ write_attribute(hid_t obj_id, const char* name, const std::array& buffer) write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); } +template inline void +write_attribute(hid_t obj_id, const char* name, const std::vector& buffer) +{ + hsize_t dims[] {buffer.size()}; + write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); +} + + //============================================================================== // Templates/overloads for write_dataset //============================================================================== diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 846d2e5a6..298b2878a 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -161,7 +161,7 @@ extern "C" { {mark_as_lost(message.str());} //! create a particle restart HDF5 file - void write_restart(); + void write_restart() const; }; diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index ce4825426..ec07a4d04 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -80,7 +80,7 @@ contains integer :: i_mesh ! flattend index for mesh logical :: energy_filters! energy filters present real(8) :: flux ! temp variable for flux - type(RegularMesh), pointer :: m ! pointer for mesh object + type(RegularMesh) :: m ! pointer for mesh object ! Extract spatial and energy indices from object nx = cmfd % indices(1) @@ -99,7 +99,7 @@ contains select type(filt => filters(i_filter_mesh) % obj) type is (MeshFilter) - m => meshes(filt % mesh) + m = meshes(filt % mesh) end select ! Set mesh widths @@ -354,9 +354,6 @@ contains ! Normalize openmc source distribution cmfd % openmc_src = cmfd % openmc_src/sum(cmfd % openmc_src)*cmfd%norm - ! Nullify all pointers - if (associated(m)) nullify(m) - end subroutine compute_xs !=============================================================================== diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 2e6162b49..eafba7b8a 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -95,7 +95,8 @@ module cmfd_header ! Main object type(cmfd_type), public :: cmfd - type(RegularMesh), public, pointer :: cmfd_mesh => null() + integer, public :: index_cmfd_mesh + type(RegularMesh), public :: cmfd_mesh ! Pointers for different tallies type(TallyContainer), public, pointer :: cmfd_tallies(:) => null() diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index de73bc03f..3208ad9d0 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -3,7 +3,7 @@ module cmfd_input use, intrinsic :: ISO_C_BINDING use cmfd_header - use mesh_header, only: mesh_dict + use mesh_header use mgxs_interface, only: energy_bins, num_energy_groups use tally use tally_header @@ -241,7 +241,7 @@ contains use constants, only: MAX_LINE_LEN use error, only: fatal_error, warning - use mesh_header, only: RegularMesh, openmc_extend_meshes + use mesh_header use string use tally, only: openmc_tally_allocate use tally_header, only: openmc_extend_tallies @@ -266,112 +266,22 @@ contains integer :: iarray3(3) ! temp integer array real(8) :: rarray3(3) ! temp double array real(C_DOUBLE), allocatable :: energies(:) - type(RegularMesh), pointer :: m type(XMLNode) :: node_mesh - err = openmc_extend_meshes(1, i_start) + ! Read CMFD mesh + call read_meshes(root % ptr) - ! Allocate mesh - cmfd_mesh => meshes(i_start) - m => meshes(i_start) + ! Get index of cmfd mesh and set ID + i_start = n_meshes() - 1 + err = openmc_mesh_set_id(i_start, i_start) - ! Set mesh id - m % id = i_start - - ! Set mesh type to rectangular - m % type = MESH_REGULAR + ! Save reference to CMFD mesh + index_cmfd_mesh = i_start + cmfd_mesh = meshes(i_start) ! Get pointer to mesh XML node node_mesh = root % child("mesh") - ! Determine number of dimensions for mesh - n = node_word_count(node_mesh, "dimension") - if (n /= 2 .and. n /= 3) then - call fatal_error("Mesh must be two or three dimensions.") - end if - m % n_dimension = n - - ! Allocate attribute arrays - allocate(m % dimension(n)) - allocate(m % lower_left(n)) - allocate(m % width(n)) - allocate(m % upper_right(n)) - - ! Check that dimensions are all greater than zero - call get_node_array(node_mesh, "dimension", iarray3(1:n)) - if (any(iarray3(1:n) <= 0)) then - call fatal_error("All entries on the element for a tally mesh& - & must be positive.") - end if - - ! Read dimensions in each direction - m % dimension = iarray3(1:n) - - ! Read mesh lower-left corner location - if (m % n_dimension /= node_word_count(node_mesh, "lower_left")) then - call fatal_error("Number of entries on must be the same as & - &the number of entries on .") - end if - call get_node_array(node_mesh, "lower_left", m % lower_left) - - ! Make sure both upper-right or width were specified - if (check_for_node(node_mesh, "upper_right") .and. & - check_for_node(node_mesh, "width")) then - call fatal_error("Cannot specify both and on a & - &tally mesh.") - end if - - ! Make sure either upper-right or width was specified - if (.not.check_for_node(node_mesh, "upper_right") .and. & - .not.check_for_node(node_mesh, "width")) then - call fatal_error("Must specify either and on a & - &tally mesh.") - end if - - if (check_for_node(node_mesh, "width")) then - ! Check to ensure width has same dimensions - if (node_word_count(node_mesh, "width") /= & - node_word_count(node_mesh, "lower_left")) then - call fatal_error("Number of entries on must be the same as the & - &number of entries on .") - end if - - ! Check for negative widths - call get_node_array(node_mesh, "width", rarray3(1:n)) - if (any(rarray3(1:n) < ZERO)) then - call fatal_error("Cannot have a negative on a tally mesh.") - end if - - ! Set width and upper right coordinate - m % width = rarray3(1:n) - m % upper_right = m % lower_left + m % dimension * m % width - - elseif (check_for_node(node_mesh, "upper_right")) then - ! Check to ensure width has same dimensions - if (node_word_count(node_mesh, "upper_right") /= & - node_word_count(node_mesh, "lower_left")) then - call fatal_error("Number of entries on must be the same & - &as the number of entries on .") - end if - - ! Check that upper-right is above lower-left - call get_node_array(node_mesh, "upper_right", rarray3(1:n)) - if (any(rarray3(1:n) < m % lower_left)) then - call fatal_error("The coordinates must be greater than & - &the coordinates on a tally mesh.") - end if - - ! Set upper right coordinate and width - m % upper_right = rarray3(1:n) - m % width = (m % upper_right - m % lower_left) / real(m % dimension, 8) - end if - - ! Set volume fraction - m % volume_frac = ONE/real(product(m % dimension),8) - - ! Add mesh to dictionary - call mesh_dict % set(m % id, i_start) - ! Determine number of filters energy_filters = check_for_node(node_mesh, "energy") n = merge(5, 3, energy_filters) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 38284c279..14fe59110 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -7,7 +7,6 @@ module eigenvalue use error, only: fatal_error, warning use math, only: t_percentile use mesh, only: count_bank_sites - use mesh_header, only: RegularMesh, meshes use message_passing use random_lcg, only: prn, set_particle_seed, advance_prn_seed use settings diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 9ce5b9ef8..ef4159e78 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -147,5 +147,10 @@ extern "C" double entropy_c(int i) return entropy.at(i - 1); } +extern "C" double entropy_clear() +{ + entropy.clear(); +} + } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index be12870ad..c3fedd881 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4,7 +4,7 @@ module input_xml use algorithm, only: find use cmfd_input, only: configure_cmfd - use cmfd_header, only: cmfd_mesh + use cmfd_header, only: index_cmfd_mesh use constants use dict_header, only: DictIntInt, DictCharInt, DictEntryCI use endf, only: reaction_name @@ -86,11 +86,6 @@ module input_xml type(C_PTR), value :: node_ptr end subroutine read_materials - subroutine read_meshes(node_ptr) bind(C) - import C_PTR - type(C_PTR), value :: node_ptr - end subroutine - function find_root_universe() bind(C) result(root) import C_INT32_T integer(C_INT32_T) :: root @@ -232,7 +227,6 @@ contains if (run_mode == MODE_EIGENVALUE) then ! Preallocate space for keff and entropy by generation call k_generation % reserve(n_max_batches*gen_per_batch) - call entropy % reserve(n_max_batches*gen_per_batch) end if ! Particle tracks @@ -1222,7 +1216,7 @@ contains ! READ MESH DATA ! Check for user meshes and allocate - call read_meshes() + call read_meshes(root % ptr) ! We only need the mesh info for plotting if (run_mode == MODE_PLOTTING) then @@ -2396,7 +2390,7 @@ contains &meshlines on plot " // trim(to_str(pl % id))) end if - pl % meshlines_mesh => cmfd_mesh + pl % index_meshlines_mesh = index_cmfd_mesh case ('entropy') diff --git a/src/mesh.cpp b/src/mesh.cpp index 9c955fae7..a4526e5b6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -532,10 +532,10 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, int n_energy, const double* energies, bool* outside) { // Determine shape of array for counts - int m = xt::prod(shape_)(); + std::size_t m = xt::prod(shape_)(); std::vector shape; if (n_energy > 0) { - shape = {m, n_energy}; + shape = {m, static_cast(n_energy)}; } else { shape = {m}; } @@ -757,6 +757,8 @@ void read_meshes(pugi::xml_node* root) //============================================================================== extern "C" { + int n_meshes() { return meshes.size(); } + RegularMesh* mesh_ptr(int i) { return meshes[i-1].get(); } int32_t mesh_id(RegularMesh* m) { return m->id_; } diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 591263000..2a6abee7d 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -2,35 +2,20 @@ module mesh_header use, intrinsic :: ISO_C_BINDING - use constants - use dict_header, only: DictIntInt - use error - use hdf5_interface - use string, only: to_str, to_lower - use xml_interface - implicit none - private - public :: openmc_extend_meshes - public :: openmc_get_mesh_index - public :: openmc_mesh_get_id - public :: openmc_mesh_get_dimension - public :: openmc_mesh_get_params - public :: openmc_mesh_set_id - public :: openmc_mesh_set_dimension - public :: openmc_mesh_set_params !=============================================================================== ! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by ! congruent squares or cubes !=============================================================================== - type, public :: RegularMesh + type :: RegularMesh type(C_PTR) :: ptr contains procedure :: id => regular_id procedure :: volume_frac => regular_volume_frac procedure :: n_dimension => regular_n_dimension + procedure :: dimension => regular_dimension procedure :: lower_left => regular_lower_left procedure :: upper_right => regular_upper_right procedure :: width => regular_width @@ -41,11 +26,6 @@ module mesh_header procedure :: get_indices_from_bin => regular_get_indices_from_bin end type RegularMesh - integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes - - ! Dictionary that maps user IDs to indices in 'meshes' - type(DictIntInt), public :: mesh_dict - interface function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C) import C_INT32_T, C_INT @@ -158,21 +138,21 @@ module mesh_header real(C_DOUBLE) :: w end function - function mesh_get_bin(ptr, xyz) result(bin) bind(C) + pure function mesh_get_bin(ptr, xyz) result(bin) bind(C) import C_PTR, C_DOUBLE, C_INT type(C_PTR), value :: ptr real(C_DOUBLE), intent(in) :: xyz(*) integer(C_INT) :: bin end function - function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C) + pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C) import C_PTR, C_INT type(C_PTR), value :: ptr integer(C_INT), intent(in) :: ijk(*) integer(C_INT) :: bin end function - subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C) + pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C) import C_PTR, C_DOUBLE, C_INT, C_BOOL type(C_PTR), value :: ptr real(C_DOUBLE), intent(in) :: xyz(*) @@ -180,7 +160,7 @@ module mesh_header logical(C_BOOL), intent(out) :: in_mesh end subroutine - subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C) + pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C) import C_PTR, C_INT type(C_PTR), value :: ptr integer(C_INT), value :: bin @@ -192,6 +172,16 @@ module mesh_header integer(C_INT), value :: i type(C_PTR) :: ptr end function + + subroutine read_meshes(node_ptr) bind(C) + import C_PTR + type(C_PTR), value :: node_ptr + end subroutine + + function n_meshes() result(n) bind(C) + import C_INT + integer(C_INT) :: n + end function end interface contains diff --git a/src/particle.cpp b/src/particle.cpp index 729a860eb..ca57c2faa 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -141,7 +141,7 @@ Particle::mark_as_lost(const char* message) } void -Particle::write_restart() +Particle::write_restart() const { // Dont write another restart file if in particle restart mode if (settings::run_mode == RUN_MODE_PARTICLE) return; diff --git a/src/physics.F90 b/src/physics.F90 index f7d80261a..f1a0d7c57 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -6,7 +6,6 @@ module physics use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math - use mesh_header, only: meshes use message_passing use nuclide_header use particle_header diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 84695153c..bcafc4314 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -7,7 +7,6 @@ module physics_mg use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math, only: rotate_angle - use mesh_header, only: meshes use mgxs_interface use message_passing use nuclide_header, only: material_xs diff --git a/src/plot.F90 b/src/plot.F90 index b8f3db079..bfc8bc225 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -9,7 +9,7 @@ module plot use hdf5_interface use output, only: time_stamp use material_header, only: materials - use mesh_header, only: meshes + use mesh_header, only: meshes, RegularMesh use particle_header use plot_header use progress_header, only: ProgressBar @@ -242,8 +242,8 @@ contains width = xyz_ur_plot - xyz_ll_plot m = meshes(pl % index_meshlines_mesh) - call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh) - call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh) + call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension()), in_mesh) + call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension()), in_mesh) ! sweep through all meshbins on this plane and draw borders do i = ijk_ll(outer), ijk_ur(outer) diff --git a/src/settings.cpp b/src/settings.cpp index 4883d07a3..db6aee7f5 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -16,6 +16,7 @@ #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/mesh.h" +#include "openmc/output.h" #include "openmc/random_lcg.h" #include "openmc/source.h" #include "openmc/string_utils.h" diff --git a/src/simulation.F90 b/src/simulation.F90 index 8b7cfdb27..3ca5d9623 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -10,8 +10,7 @@ module simulation use cmfd_execute, only: cmfd_init_batch, cmfd_tally_init, execute_cmfd use cmfd_header, only: cmfd_on use constants, only: ZERO - use eigenvalue, only: calculate_average_keff, & - calculate_generation_keff, shannon_entropy, & + use eigenvalue, only: calculate_average_keff, calculate_generation_keff, & synchronize_bank, keff_generation, k_sum #ifdef _OPENMP use eigenvalue, only: join_bank_from_threads @@ -479,7 +478,7 @@ contains ! will potentially populate k_generation and entropy) current_batch = 0 call k_generation % clear() - call entropy % clear() + call entropy_clear() need_depletion_rx = .false. ! If this is a restart run, load the state point data and binary source diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 9b7c68cd4..ad5513c61 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -64,6 +64,12 @@ module simulation_header !$omp threadprivate(trace, thread_id, current_work) + interface + subroutine entropy_clear() bind(C) + end subroutine + end interface + + contains !=============================================================================== @@ -80,12 +86,12 @@ contains !=============================================================================== subroutine free_memory_simulation() + if (allocated(work_index)) deallocate(work_index) call k_generation % clear() call k_generation % shrink_to_fit() - call entropy % clear() - call entropy % shrink_to_fit() + call entropy_clear() end subroutine free_memory_simulation end module simulation_header diff --git a/src/state_point.F90 b/src/state_point.F90 index b7f4e1377..d20456d7b 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -173,7 +173,7 @@ contains call write_dataset(file_id, "generations_per_batch", gen_per_batch) k = k_generation % size() call write_dataset(file_id, "k_generation", k_generation % data(1:k)) - call entropy_to_hdf5() + call entropy_to_hdf5(file_id) call write_dataset(file_id, "k_col_abs", k_col_abs) call write_dataset(file_id, "k_col_tra", k_col_tra) call write_dataset(file_id, "k_abs_tra", k_abs_tra) diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 28eafcaaf..3e8eaa1c6 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -92,155 +92,154 @@ contains real(8) :: distance ! distance traveled in mesh cell logical :: start_in_mesh ! starting coordinates inside mesh? logical :: end_in_mesh ! ending coordinates inside mesh? - type(RegularMesh), pointer :: m + type(RegularMesh) :: m ! Get a pointer to the mesh. - m => meshes(this % mesh) - n = m % n_dimension + m = meshes(this % mesh) if (estimator /= ESTIMATOR_TRACKLENGTH) then ! If this is an analog or collision tally, then there can only be one ! valid mesh bin. call m % get_bin(p % coord(1) % xyz, bin) - if (bin /= NO_BIN_FOUND) then + if (bin >= 0) then call match % bins % push_back(bin) call match % weights % push_back(ONE) end if return end if - ! A track can span multiple mesh bins so we need to handle a lot of - ! intersection logic for tracklength tallies. + ! ! A track can span multiple mesh bins so we need to handle a lot of + ! ! intersection logic for tracklength tallies. - ! ======================================================================== - ! Determine if the track intersects the tally mesh. + ! ! ======================================================================== + ! ! Determine if the track intersects the tally mesh. - ! Copy the starting and ending coordinates of the particle. Offset these - ! just a bit for the purposes of determining if there was an intersection - ! in case the mesh surfaces coincide with lattice/geometric surfaces which - ! might produce finite-precision errors. - xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw - xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw + ! ! Copy the starting and ending coordinates of the particle. Offset these + ! ! just a bit for the purposes of determining if there was an intersection + ! ! in case the mesh surfaces coincide with lattice/geometric surfaces which + ! ! might produce finite-precision errors. + ! xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw + ! xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - ! Determine indices for starting and ending location. - call m % get_indices(xyz0, ijk0(:n), start_in_mesh) - call m % get_indices(xyz1, ijk1(:n), end_in_mesh) + ! ! Determine indices for starting and ending location. + ! call m % get_indices(xyz0, ijk0(:n), start_in_mesh) + ! call m % get_indices(xyz1, ijk1(:n), end_in_mesh) - ! If this is the first iteration of the filter loop, check if the track - ! intersects any part of the mesh. - if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then - if (.not. m % intersects(xyz0, xyz1)) return - end if + ! ! If this is the first iteration of the filter loop, check if the track + ! ! intersects any part of the mesh. + ! if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then + ! if (.not. m % intersects(xyz0, xyz1)) return + ! end if - ! ======================================================================== - ! Figure out which mesh cell to tally. + ! ! ======================================================================== + ! ! Figure out which mesh cell to tally. - ! Copy the un-modified coordinates the particle direction. - xyz0 = p % last_xyz - xyz1 = p % coord(1) % xyz - uvw = p % coord(1) % uvw + ! ! Copy the un-modified coordinates the particle direction. + ! xyz0 = p % last_xyz + ! xyz1 = p % coord(1) % xyz + ! uvw = p % coord(1) % uvw - ! Compute the length of the entire track. - total_distance = sqrt(sum((xyz1 - xyz0)**2)) + ! ! Compute the length of the entire track. + ! total_distance = sqrt(sum((xyz1 - xyz0)**2)) - ! We are looking for the first valid mesh bin. Check to see if the - ! particle starts inside the mesh. - if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then - ! The particle does not start in the mesh. Note that we nudged the - ! start and end coordinates by a TINY_BIT each so we will have - ! difficulty resolving tracks that are less than 2*TINY_BIT in length. - ! If the track is that short, it is also insignificant so we can - ! safely ignore it in the tallies. - if (total_distance < 2*TINY_BIT) return + ! ! We are looking for the first valid mesh bin. Check to see if the + ! ! particle starts inside the mesh. + ! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then + ! ! The particle does not start in the mesh. Note that we nudged the + ! ! start and end coordinates by a TINY_BIT each so we will have + ! ! difficulty resolving tracks that are less than 2*TINY_BIT in length. + ! ! If the track is that short, it is also insignificant so we can + ! ! safely ignore it in the tallies. + ! if (total_distance < 2*TINY_BIT) return - ! The particle does not start in the mesh so keep iterating the ijk0 - ! indices to cross the nearest mesh surface until we've found a valid - ! bin. MAX_SEARCH_ITER prevents an infinite loop. - search_iter = 0 - do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) - if (search_iter == MAX_SEARCH_ITER) then - call warning("Failed to find a mesh intersection on a tally mesh & - &filter.") - return - end if + ! ! The particle does not start in the mesh so keep iterating the ijk0 + ! ! indices to cross the nearest mesh surface until we've found a valid + ! ! bin. MAX_SEARCH_ITER prevents an infinite loop. + ! search_iter = 0 + ! do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) + ! if (search_iter == MAX_SEARCH_ITER) then + ! call warning("Failed to find a mesh intersection on a tally mesh & + ! &filter.") + ! return + ! end if - do j = 1, n - if (abs(uvw(j)) < FP_PRECISION) then - d(j) = INFINITY - else if (uvw(j) > 0) then - xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - else - xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - end if - end do - j = minloc(d(:n), 1) - if (uvw(j) > ZERO) then - ijk0(j) = ijk0(j) + 1 - else - ijk0(j) = ijk0(j) - 1 - end if + ! do j = 1, n + ! if (abs(uvw(j)) < FP_PRECISION) then + ! d(j) = INFINITY + ! else if (uvw(j) > 0) then + ! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! else + ! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! end if + ! end do + ! j = minloc(d(:n), 1) + ! if (uvw(j) > ZERO) then + ! ijk0(j) = ijk0(j) + 1 + ! else + ! ijk0(j) = ijk0(j) - 1 + ! end if - search_iter = search_iter + 1 - end do - distance = d(j) - xyz0 = xyz0 + distance * uvw - end if + ! search_iter = search_iter + 1 + ! end do + ! distance = d(j) + ! xyz0 = xyz0 + distance * uvw + ! end if - do - ! ======================================================================== - ! Compute the length of the track segment in the appropiate mesh cell and - ! return. + ! do + ! ! ======================================================================== + ! ! Compute the length of the track segment in the appropiate mesh cell and + ! ! return. - if (all(ijk0(:n) == ijk1(:n))) then - ! The track ends in this cell. Use the particle end location rather - ! than the mesh surface. - distance = sqrt(sum((xyz1 - xyz0)**2)) - else - ! The track exits this cell. Determine the distance to the closest mesh - ! surface. - do j = 1, n - if (abs(uvw(j)) < FP_PRECISION) then - d(j) = INFINITY - else if (uvw(j) > 0) then - xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - else - xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - end if - end do - j = minloc(d(:n), 1) - distance = d(j) - end if + ! if (all(ijk0(:n) == ijk1(:n))) then + ! ! The track ends in this cell. Use the particle end location rather + ! ! than the mesh surface. + ! distance = sqrt(sum((xyz1 - xyz0)**2)) + ! else + ! ! The track exits this cell. Determine the distance to the closest mesh + ! ! surface. + ! do j = 1, n + ! if (abs(uvw(j)) < FP_PRECISION) then + ! d(j) = INFINITY + ! else if (uvw(j) > 0) then + ! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! else + ! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! end if + ! end do + ! j = minloc(d(:n), 1) + ! distance = d(j) + ! end if - ! Assign the next tally bin and the score. - bin = m % get_bin_from_indices(ijk0(:n)) - call match % bins % push_back(bin) - call match % weights % push_back(distance / total_distance) + ! ! Assign the next tally bin and the score. + ! bin = m % get_bin_from_indices(ijk0(:n)) + ! call match % bins % push_back(bin) + ! call match % weights % push_back(distance / total_distance) - ! Find the next mesh cell that the particle enters. + ! ! Find the next mesh cell that the particle enters. - ! If the particle track ends in that bin, then we are done. - if (all(ijk0(:n) == ijk1(:n))) exit + ! ! If the particle track ends in that bin, then we are done. + ! if (all(ijk0(:n) == ijk1(:n))) exit - ! Translate the starting coordintes by the distance to that face. This - ! should be the xyz that we computed the distance to in the last - ! iteration of the filter loop. - xyz0 = xyz0 + distance * uvw + ! ! Translate the starting coordintes by the distance to that face. This + ! ! should be the xyz that we computed the distance to in the last + ! ! iteration of the filter loop. + ! xyz0 = xyz0 + distance * uvw - ! Increment the indices into the next mesh cell. - if (uvw(j) > ZERO) then - ijk0(j) = ijk0(j) + 1 - else - ijk0(j) = ijk0(j) - 1 - end if + ! ! Increment the indices into the next mesh cell. + ! if (uvw(j) > ZERO) then + ! ijk0(j) = ijk0(j) + 1 + ! else + ! ijk0(j) = ijk0(j) - 1 + ! end if - ! If the next indices are invalid, then the track has left the mesh and - ! we are done. - if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit - end do + ! ! If the next indices are invalid, then the track has left the mesh and + ! ! we are done. + ! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit + ! end do end subroutine get_all_bins_mesh @@ -248,9 +247,12 @@ contains class(MeshFilter), intent(in) :: this integer(HID_T), intent(in) :: filter_group + type(RegularMesh) :: m + + m = meshes(this % mesh) call write_dataset(filter_group, "type", "mesh") call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", meshes(this % mesh) % id) + call write_dataset(filter_group, "bins", m % id()) end subroutine to_statepoint_mesh function text_label_mesh(this, bin) result(label) @@ -259,20 +261,20 @@ contains character(MAX_LINE_LEN) :: label integer, allocatable :: ijk(:) + type(RegularMesh) :: m - associate (m => meshes(this % mesh)) - allocate(ijk(m % n_dimension)) - call m % get_indices_from_bin(bin, ijk) - if (m % n_dimension == 1) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" - elseif (m % n_dimension == 2) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" - elseif (m % n_dimension == 3) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" - end if - end associate + m = meshes(this % mesh) + allocate(ijk(m % n_dimension())) + call m % get_indices_from_bin(bin, ijk) + if (m % n_dimension() == 1) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" + elseif (m % n_dimension() == 2) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ")" + elseif (m % n_dimension() == 3) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + end if end function text_label_mesh !=============================================================================== @@ -304,18 +306,25 @@ contains integer(C_INT32_T), value, intent(in) :: index_mesh integer(C_INT) :: err + type(RegularMesh) :: m + integer :: i + err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (MeshFilter) - if (index_mesh >= 1 .and. index_mesh <= n_meshes) then + if (index_mesh >= 0 .and. index_mesh < n_meshes()) then f % mesh = index_mesh - f % n_bins = product(meshes(index_mesh) % dimension) + f % n_bins = 1 + m = meshes(index_mesh) + do i = 1, m % n_dimension() + f % n_bins = f % n_bins * m % dimension(i) + end do else err = E_OUT_OF_BOUNDS call set_errmsg("Index in 'meshes' array is out of bounds.") end if - class default + class default err = E_INVALID_TYPE call set_errmsg("Tried to set mesh on a non-mesh filter.") end select diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index c7363a2db..d1500b8b3 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -42,6 +42,7 @@ contains integer :: id integer :: n integer :: n_dim + integer(C_INT) :: err type(RegularMesh) :: m n = node_word_count(node, "bins") @@ -90,133 +91,130 @@ contains real(8) :: distance ! actual distance traveled logical :: start_in_mesh ! particle's starting xyz in mesh? logical :: end_in_mesh ! particle's ending xyz in mesh? + type(RegularMesh) :: m - ! Copy starting and ending location of particle - xyz0 = p % last_xyz_current - xyz1 = p % coord(1) % xyz + ! ! Copy starting and ending location of particle + ! xyz0 = p % last_xyz_current + ! xyz1 = p % coord(1) % xyz - associate (m => meshes(this % mesh)) - n_dim = m % n_dimension + ! m = meshes(this % mesh) + ! n_dim = m % n_dimension() - ! Determine indices for starting and ending location - call m % get_indices(xyz0, ijk0, start_in_mesh) - call m % get_indices(xyz1, ijk1, end_in_mesh) + ! ! Determine indices for starting and ending location + ! call m % get_indices(xyz0, ijk0, start_in_mesh) + ! call m % get_indices(xyz1, ijk1, end_in_mesh) - ! Check to see 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. m % intersects(xyz0, xyz1)) return - end if + ! ! Check to see 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. m % intersects(xyz0, xyz1)) return + ! end if - ! Calculate number of surface crossings - n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim))) - if (n_cross == 0) return + ! ! Calculate number of surface crossings + ! n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim))) + ! if (n_cross == 0) return - ! Copy particle's direction - uvw = p % coord(1) % uvw + ! ! Copy particle's direction + ! uvw = p % coord(1) % uvw - ! Bounding coordinates - do d1 = 1, n_dim - if (uvw(d1) > 0) then - xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1) - else - xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1) - end if - end do + ! ! Bounding coordinates + ! do d1 = 1, n_dim + ! if (uvw(d1) > 0) then + ! xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1) + ! else + ! xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1) + ! end if + ! end do - do j = 1, n_cross - ! Set the distances to infinity - d = INFINITY + ! do j = 1, n_cross + ! ! Set the distances to infinity + ! d = INFINITY - ! Calculate distance to each bounding surface. We need to treat - ! special case where the cosine of the angle is zero since this would - ! result in a divide-by-zero. - do d1 = 1, n_dim - if (uvw(d1) == 0) then - d(d1) = INFINITY - else - d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1) - end if - end do + ! ! Calculate distance to each bounding surface. We need to treat + ! ! special case where the cosine of the angle is zero since this would + ! ! result in a divide-by-zero. + ! do d1 = 1, n_dim + ! if (uvw(d1) == 0) then + ! d(d1) = INFINITY + ! else + ! d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1) + ! end if + ! end do - ! Determine the closest bounding surface of the mesh cell by - ! calculating the minimum distance. Then use the minimum distance and - ! direction of the particle to determine which surface was crossed. - distance = minval(d) + ! ! Determine the closest bounding surface of the mesh cell by + ! ! calculating the minimum distance. Then use the minimum distance and + ! ! direction of the particle to determine which surface was crossed. + ! distance = minval(d) - ! Loop over the dimensions - do d1 = 1, n_dim + ! ! Loop over the dimensions + ! do d1 = 1, n_dim - ! Check whether distance is the shortest distance - if (distance == d(d1)) then + ! ! Check whether distance is the shortest distance + ! if (distance == d(d1)) then - ! Check whether particle is moving in positive d1 direction - if (uvw(d1) > 0) then + ! ! Check whether particle is moving in positive d1 direction + ! if (uvw(d1) > 0) then - ! Outward current on d1 max surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - 1 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! Outward current on d1 max surface + ! if (all(ijk0(:n_dim) >= 1) .and. & + ! all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 - 1 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if - ! Advance position - ijk0(d1) = ijk0(d1) + 1 - xyz_cross(d1) = xyz_cross(d1) + m % width(d1) + ! ! Advance position + ! ijk0(d1) = ijk0(d1) + 1 + ! xyz_cross(d1) = xyz_cross(d1) + m % width(d1) - ! If the particle crossed the surface, tally the inward current on - ! d1 min surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - 2 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! If the particle crossed the surface, tally the inward current on + ! ! d1 min surface + ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 - 2 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if - else - ! The particle is moving in the negative d1 direction + ! else + ! ! The particle is moving in the negative d1 direction - ! Outward current on d1 min surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - 3 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! Outward current on d1 min surface + ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 - 3 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if - ! Advance position - ijk0(d1) = ijk0(d1) - 1 - xyz_cross(d1) = xyz_cross(d1) - m % width(d1) + ! ! Advance position + ! ijk0(d1) = ijk0(d1) - 1 + ! xyz_cross(d1) = xyz_cross(d1) - m % width(d1) - ! If the particle crossed the surface, tally the inward current on - ! d1 max surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! If the particle crossed the surface, tally the inward current on + ! ! d1 max surface + ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if - end if - end if - end do + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if + ! end if + ! end if + ! end do - ! Calculate new coordinates - xyz0 = xyz0 + distance * uvw - end do - end associate + ! ! Calculate new coordinates + ! xyz0 = xyz0 + distance * uvw + ! end do end subroutine get_all_bins @@ -224,9 +222,12 @@ contains class(MeshSurfaceFilter), intent(in) :: this integer(HID_T), intent(in) :: filter_group + type(RegularMesh) :: m + + m = meshes(this % mesh) call write_dataset(filter_group, "type", "meshsurface") call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", meshes(this % mesh) % id) + call write_dataset(filter_group, "bins", m % id()) end subroutine to_statepoint function text_label(this, bin) result(label) @@ -238,55 +239,55 @@ contains integer :: i_surf integer :: n_dim integer, allocatable :: ijk(:) + type(RegularMesh) :: m - associate (m => meshes(this % mesh)) - n_dim = m % n_dimension - allocate(ijk(n_dim)) + m = meshes(this % mesh) + n_dim = m % n_dimension() + allocate(ijk(n_dim)) - ! Get flattend mesh index and surface index - i_mesh = (bin - 1) / (4*n_dim) + 1 - i_surf = mod(bin - 1, 4*n_dim) + 1 + ! Get flattend mesh index and surface index + i_mesh = (bin - 1) / (4*n_dim) + 1 + i_surf = mod(bin - 1, 4*n_dim) + 1 - ! Get mesh index part of label - call m % get_indices_from_bin(i_mesh, ijk) - if (m % n_dimension == 1) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" - elseif (m % n_dimension == 2) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" - elseif (m % n_dimension == 3) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" - end if + ! Get mesh index part of label + call m % get_indices_from_bin(i_mesh, ijk) + if (m % n_dimension() == 1) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" + elseif (m % n_dimension() == 2) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ")" + elseif (m % n_dimension() == 3) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + end if - ! Get surface part of label - select case (i_surf) - case (OUT_LEFT) - label = trim(label) // " Outgoing, x-min" - case (IN_LEFT) - label = trim(label) // " Incoming, x-min" - case (OUT_RIGHT) - label = trim(label) // " Outgoing, x-max" - case (IN_RIGHT) - label = trim(label) // " Incoming, x-max" - case (OUT_BACK) - label = trim(label) // " Outgoing, y-min" - case (IN_BACK) - label = trim(label) // " Incoming, y-min" - case (OUT_FRONT) - label = trim(label) // " Outgoing, y-max" - case (IN_FRONT) - label = trim(label) // " Incoming, y-max" - case (OUT_BOTTOM) - label = trim(label) // " Outgoing, z-min" - case (IN_BOTTOM) - label = trim(label) // " Incoming, z-min" - case (OUT_TOP) - label = trim(label) // " Outgoing, z-max" - case (IN_TOP) - label = trim(label) // " Incoming, z-max" - end select - end associate + ! Get surface part of label + select case (i_surf) + case (OUT_LEFT) + label = trim(label) // " Outgoing, x-min" + case (IN_LEFT) + label = trim(label) // " Incoming, x-min" + case (OUT_RIGHT) + label = trim(label) // " Outgoing, x-max" + case (IN_RIGHT) + label = trim(label) // " Incoming, x-max" + case (OUT_BACK) + label = trim(label) // " Outgoing, y-min" + case (IN_BACK) + label = trim(label) // " Incoming, y-min" + case (OUT_FRONT) + label = trim(label) // " Outgoing, y-max" + case (IN_FRONT) + label = trim(label) // " Incoming, y-max" + case (OUT_BOTTOM) + label = trim(label) // " Outgoing, z-min" + case (IN_BOTTOM) + label = trim(label) // " Incoming, z-min" + case (OUT_TOP) + label = trim(label) // " Outgoing, z-max" + case (IN_TOP) + label = trim(label) // " Incoming, z-max" + end select end function text_label !=============================================================================== @@ -318,16 +319,22 @@ contains integer(C_INT32_T), value, intent(in) :: index_mesh integer(C_INT) :: err + integer :: i integer :: n_dim + type(RegularMesh) :: m err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (MeshSurfaceFilter) - if (index_mesh >= 1 .and. index_mesh <= n_meshes) then + if (index_mesh >= 0 .and. index_mesh < n_meshes()) then f % mesh = index_mesh - n_dim = meshes(index_mesh) % n_dimension - f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension) + m = meshes(index_mesh) + n_dim = m % n_dimension() + f % n_bins = 4*n_dim + do i = 1, n_dim + f % n_bins = f % n_bins * m % dimension(i) + end do else err = E_OUT_OF_BOUNDS call set_errmsg("Index in 'meshes' array is out of bounds.") diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 2cb3bf819..27aa5ec50 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -257,14 +257,14 @@ contains logical :: print_ebin ! should incoming energy bin be displayed? real(8) :: rel_err = ZERO ! temporary relative error of result real(8) :: std_dev = ZERO ! temporary standard deviration of result - type(RegularMesh), pointer :: m ! surface current mesh + type(RegularMesh) :: m ! surface current mesh ! Get pointer to mesh i_filter_mesh = t % filter(t % find_filter(FILTER_MESH)) i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE)) select type(filt => filters(i_filter_mesh) % obj) type is (MeshFilter) - m => meshes(filt % mesh) + m = meshes(filt % mesh) end select ! initialize bins array @@ -285,8 +285,11 @@ contains end if ! Get the dimensions and number of cells in the mesh - n_dim = m % n_dimension - n_cells = product(m % dimension) + n_dim = m % n_dimension() + n_cells = 1 + do j = 1, n_dim + n_cells = n_cells * m % dimension(j) + end do ! Loop over all the mesh cells do i = 1, n_cells