diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index ac68f5b00e..34ed206da5 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -5,6 +5,7 @@ #define OPENMC_GEOMETRY_AUX_H #include +#include namespace openmc { @@ -41,7 +42,7 @@ extern "C" void neighbor_lists(); //! Populate all data structures needed for distribcells. //============================================================================== -extern "C" void prepare_distribcell(int32_t* filter_cell_list, int n); +extern "C" void prepare_distribcell(); //============================================================================== //! Recursively search through the geometry and count cell instances. @@ -65,21 +66,6 @@ extern "C" void count_cell_instances(int32_t univ_indx); extern "C" int count_universe_instances(int32_t search_univ, int32_t target_univ_id); -//============================================================================== -//! Find the length necessary for a string to contain a distribcell path. -//! @param target_cell The index of the Cell in the global Cell array. -//! @param map The index of the distribcell mapping corresponding to the target -//! cell. -//! @param target_offset An instance number for a distributed cell. -//! @param root_univ The index of the root Universe in the global Universe -//! array. -//! @return The size of a character array needed to fit the distribcell path. -//============================================================================== - -extern "C" int -distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ); - //============================================================================== //! Build a character array representing the path to a distribcell instance. //! @param target_cell The index of the Cell in the global Cell array. @@ -88,13 +74,12 @@ distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, //! @param target_offset An instance number for a distributed cell. //! @param root_univ The index of the root Universe in the global Universe //! array. -//! @param[out] path The unique traversal through the geometry tree that leads -//! to the desired instance of the target cell. +//! @return The unique traversal through the geometry tree that leads to the +//! desired instance of the target cell. //============================================================================== -extern "C" void -distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ, char *path); +std::string +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset); //============================================================================== //! Determine the maximum number of nested coordinate levels in the geometry. diff --git a/include/openmc/tallies/tally_filter_distribcell.h b/include/openmc/tallies/tally_filter_distribcell.h new file mode 100644 index 0000000000..b29520caa7 --- /dev/null +++ b/include/openmc/tallies/tally_filter_distribcell.h @@ -0,0 +1,93 @@ +#ifndef OPENMC_TALLY_FILTER_DISTRIBCELL_H +#define OPENMC_TALLY_FILTER_DISTRIBCELL_H + +#include + +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/geometry_aux.h" // For distribcell_path +#include "openmc/lattice.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +class DistribcellFilter : public TallyFilter +{ +public: + virtual ~DistribcellFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + auto cells = get_node_array(node, "bins"); + if (cells.size() != 1) { + fatal_error("Only one cell can be specified per distribcell filter."); + } + cell_ = cells[0]; + } + + virtual void + initialize() override + { + auto search = cell_map.find(cell_); + if (search != cell_map.end()) { + cell_ = search->second; + n_bins_ = cells[cell_]->n_instances_; + } else { + std::stringstream err_msg; + err_msg << "Could not find cell " << cell_ + << " specified on tally filter."; + fatal_error(err_msg); + } + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + int offset = 0; + auto distribcell_index = cells[cell_]->distribcell_index_; + for (int i = 0; i < p->n_coord; i++) { + auto& c {*cells[p->coord[i].cell]}; + if (c.type_ == FILL_UNIVERSE) { + offset += c.offset_[distribcell_index]; + } else if (c.type_ == FILL_LATTICE) { + auto& lat {*lattices[p->coord[i+1].lattice-1]}; + int i_xyz[3] {p->coord[i+1].lattice_x, + p->coord[i+1].lattice_y, + p->coord[i+1].lattice_z}; + if (lat.are_valid_indices(i_xyz)) { + offset += lat.offset(distribcell_index, i_xyz); + } + } + if (cell_ == p->coord[i].cell) { + match.bins.push_back(offset + 1); + match.weights.push_back(1); + return; + } + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + write_dataset(filter_group, "type", "distribcell"); + write_dataset(filter_group, "n_bins", n_bins_); + write_dataset(filter_group, "bins", cells[cell_]->id_); + } + + virtual std::string + text_label(int bin) const override + { + auto map = cells[cell_]->distribcell_index_; + auto path = distribcell_path(cell_, map, bin-1); + + return "Distributed Cell " + path; + } + + int32_t cell_; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_DISTRIBCELL_H diff --git a/src/cell.cpp b/src/cell.cpp index 41a15aefdc..545bb7f58f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -803,10 +803,6 @@ extern "C" { int32_t cell_fill(Cell* c) {return c->fill_;} - int32_t cell_n_instances(Cell* c) {return c->n_instances_;} - - int cell_distribcell_index(Cell* c) {return c->distribcell_index_;} - int cell_material_size(Cell* c) {return c->material_.size();} //TODO: off-by-one @@ -821,8 +817,6 @@ extern "C" { double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT_[i];} - int32_t cell_offset(Cell* c, int map) {return c->offset_[map];} - void extend_cells_c(int32_t n) { cells.reserve(cells.size() + n); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 20f64e1876..1da7bd78b7 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -12,6 +12,8 @@ #include "openmc/material.h" #include "openmc/settings.h" #include "openmc/surface.h" +#include "openmc/tallies/tally_filter.h" +#include "openmc/tallies/tally_filter_distribcell.h" namespace openmc { @@ -181,12 +183,16 @@ neighbor_lists() //============================================================================== void -prepare_distribcell(int32_t* filter_cell_list, int n) +prepare_distribcell() { - // Read the list of cells contained in distribcell filters from Fortran. + // Find all cells listed in a DistribcellFilter. std::unordered_set distribcells; - for (int i = 0; i < n; i++) { - distribcells.insert(filter_cell_list[i]); + for (auto* filt : tally_filters) { + //TODO: replace this cast + auto* distrib_filt = dynamic_cast(filt); + if (distrib_filt) { + distribcells.insert(distrib_filt->cell_); + } } // Find all cells with distributed materials or temperatures. Make sure that @@ -393,29 +399,11 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, } } -//============================================================================== - -int -distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ) +std::string +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset) { - Universe& root = *universes[root_univ]; - std::string path_ {distribcell_path_inner(target_cell, map, target_offset, - root, 0)}; - return path_.size() + 1; -} - -//============================================================================== - -void -distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ, char* path) -{ - Universe& root = *universes[root_univ]; - std::string path_ {distribcell_path_inner(target_cell, map, target_offset, - root, 0)}; - path_.copy(path, path_.size()); - path[path_.size()] = '\0'; + auto& root_univ = *universes[openmc_root_universe]; + return distribcell_path_inner(target_cell, map, target_offset, root_univ, 0); } //============================================================================== diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index b58b5415ba..89cfbcf2a7 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -59,20 +59,6 @@ module geometry_header integer(C_INT32_T) :: fill end function cell_fill_c - function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') & - result(n_instances) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT32_T) :: n_instances - end function cell_n_instances_c - - function cell_distribcell_index_c(cell_ptr) & - bind(C, name='cell_distribcell_index') result(distribcell_index) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT) :: distribcell_index - end function cell_distribcell_index_c - function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') & result(n) import C_PTR, C_INT @@ -103,14 +89,6 @@ module geometry_header real(C_DOUBLE) :: sqrtkT end function cell_sqrtkT_c - function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") & - result(offset) - import C_PTR, C_INT, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT), intent(in), value :: map - integer(C_INT32_T) :: offset - end function cell_offset_c - function lattice_pointer(lat_ind) bind(C) result(ptr) import C_PTR, C_INT32_T integer(C_INT32_T), intent(in), value :: lat_ind @@ -123,23 +101,6 @@ module geometry_header integer(C_INT32_T) :: id end function lattice_id_c - function lattice_are_valid_indices_c(lat_ptr, i_xyz) & - bind(C, name='lattice_are_valid_indices') result (is_valid) - import C_PTR, C_INT, C_BOOL - type(C_PTR), intent(in), value :: lat_ptr - integer(C_INT), intent(in) :: i_xyz(3) - logical(C_BOOL) :: is_valid - end function lattice_are_valid_indices_c - - function lattice_offset_c(lat_ptr, map, i_xyz) & - bind(C, name='lattice_offset') result(offset) - import C_PTR, C_INT, C_INT32_T - type(C_PTR), intent(in), value :: lat_ptr - integer(C_INT), intent(in), value :: map - integer(C_INT), intent(in) :: i_xyz(3) - integer(C_INT32_T) :: offset - end function lattice_offset_c - subroutine extend_cells_c(n) bind(C) import C_INT32_t integer(C_INT32_T), intent(in), value :: n @@ -154,8 +115,6 @@ module geometry_header type(C_PTR) :: ptr contains procedure :: id => lattice_id - procedure :: are_valid_indices => lattice_are_valid_indices - procedure :: offset => lattice_offset end type Lattice !=============================================================================== @@ -176,13 +135,10 @@ module geometry_header procedure :: type => cell_type procedure :: universe => cell_universe procedure :: fill => cell_fill - procedure :: n_instances => cell_n_instances - procedure :: distribcell_index => cell_distribcell_index procedure :: material_size => cell_material_size procedure :: material => cell_material procedure :: sqrtkT_size => cell_sqrtkT_size procedure :: sqrtkT => cell_sqrtkT - procedure :: offset => cell_offset end type Cell @@ -208,21 +164,6 @@ contains id = lattice_id_c(this % ptr) end function lattice_id - function lattice_are_valid_indices(this, i_xyz) result (is_valid) - class(Lattice), intent(in) :: this - integer(C_INT), intent(in) :: i_xyz(3) - logical(C_BOOL) :: is_valid - is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz) - end function lattice_are_valid_indices - - function lattice_offset(this, map, i_xyz) result(offset) - class(Lattice), intent(in) :: this - integer(C_INT), intent(in) :: map - integer(C_INT), intent(in) :: i_xyz(3) - integer(C_INT32_T) :: offset - offset = lattice_offset_c(this % ptr, map, i_xyz) - end function lattice_offset - !=============================================================================== function cell_id(this) result(id) @@ -255,18 +196,6 @@ contains fill = cell_fill_c(this % ptr) end function cell_fill - function cell_n_instances(this) result(n_instances) - class(Cell), intent(in) :: this - integer(C_INT32_T) :: n_instances - n_instances = cell_n_instances_c(this % ptr) - end function cell_n_instances - - function cell_distribcell_index(this) result(distribcell_index) - class(Cell), intent(in) :: this - integer(C_INT) :: distribcell_index - distribcell_index = cell_distribcell_index_c(this % ptr) - end function cell_distribcell_index - function cell_material_size(this) result(n) class(Cell), intent(in) :: this integer(C_INT) :: n @@ -293,13 +222,6 @@ contains sqrtkT = cell_sqrtkT_c(this % ptr, i) end function cell_sqrtkT - function cell_offset(this, map) result(offset) - class(Cell), intent(in) :: this - integer(C_INT), intent(in) :: map - integer(C_INT32_T) :: offset - offset = cell_offset_c(this % ptr, map) - end function cell_offset - !=============================================================================== ! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table ! appears at in the model. Later, this list is used to determine the actual diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 49c6352355..e6ebe30587 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -59,12 +59,8 @@ module input_xml integer(C_INT32_T), intent(in), value :: univ_indx end subroutine count_cell_instances - subroutine prepare_distribcell_c(cell_list, n) & - bind(C, name="prepare_distribcell") - import C_INT32_T, C_INT - integer(C_INT), intent(in), value :: n - integer(C_INT32_T), intent(in) :: cell_list(n) - end subroutine prepare_distribcell_c + subroutine prepare_distribcell() bind(C) + end subroutine prepare_distribcell subroutine read_surfaces(node_ptr) bind(C) import C_PTR @@ -3148,33 +3144,4 @@ contains end subroutine read_multipole_data -!=============================================================================== -! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the -! offsets for distribcells -!=============================================================================== - - subroutine prepare_distribcell() - - integer :: i, j - type(SetInt) :: cell_list ! distribcells to track - integer(C_INT32_T), allocatable :: cell_list_c(:) - - ! Find all cells listed in a distribcell filter. - do i = 1, n_tallies - do j = 1, size(tallies(i) % obj % filter) - select type(filt => filters(tallies(i) % obj % filter(j)) % obj) - type is (DistribcellFilter) - call cell_list % add(filt % cell) - end select - end do - end do - - allocate(cell_list_c(cell_list % size())) - do i = 1, cell_list % size() - cell_list_c(i) = cell_list % get_item(i) - 1 - end do - call prepare_distribcell_c(cell_list_c, cell_list % size()) - - end subroutine prepare_distribcell - end module input_xml diff --git a/src/lattice.cpp b/src/lattice.cpp index 21f4768506..14b679903d 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -893,12 +893,6 @@ extern "C" { Lattice* lattice_pointer(int lat_ind) {return lattices[lat_ind];} int32_t lattice_id(Lattice *lat) {return lat->id_;} - - bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3]) - {return lat->are_valid_indices(i_xyz);} - - int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) - {return lat->offset(map, i_xyz);} } } // namespace openmc diff --git a/src/tallies/tally_filter.cpp b/src/tallies/tally_filter.cpp index 944b154033..b0f54dd1ee 100644 --- a/src/tallies/tally_filter.cpp +++ b/src/tallies/tally_filter.cpp @@ -6,6 +6,7 @@ #include "openmc/tallies/tally_filter_cell.h" #include "openmc/tallies/tally_filter_cellborn.h" #include "openmc/tallies/tally_filter_cellfrom.h" +#include "openmc/tallies/tally_filter_distribcell.h" #include "openmc/tallies/tally_filter_mesh.h" #include "openmc/tallies/tally_filter_meshsurface.h" @@ -85,6 +86,8 @@ extern "C" { tally_filters.push_back(new CellbornFilter()); } else if (type_ == "cellfrom") { tally_filters.push_back(new CellFromFilter()); + } else if (type_ == "distribcell") { + tally_filters.push_back(new DistribcellFilter()); } else if (type_ == "mesh") { tally_filters.push_back(new MeshFilter()); } else if (type_ == "meshsurface") { diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index b6975c8a12..4a4b286852 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -1,16 +1,6 @@ module tally_filter_distribcell - use, intrinsic :: ISO_C_BINDING - - use constants - use dict_header, only: EMPTY - use error - use geometry_header - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str use tally_filter_header - use xml_interface implicit none private @@ -20,150 +10,18 @@ module tally_filter_distribcell ! reside in. !=============================================================================== - type, public, extends(TallyFilter) :: DistribcellFilter + type, public, extends(CppTallyFilter) :: DistribcellFilter integer :: cell contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_distribcell - procedure :: to_statepoint => to_statepoint_distribcell - procedure :: text_label => text_label_distribcell procedure :: initialize => initialize_distribcell end type DistribcellFilter contains - subroutine from_xml(this, node) - class(DistribcellFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - n = node_word_count(node, "bins") - if (n /= 1) call fatal_error("Only one cell can be & - &specified per distribcell filter.") - - ! Store bins - call get_node_value(node, "bins", this % cell) - end subroutine from_xml - - subroutine get_all_bins_distribcell(this, p, estimator, match) - class(DistribcellFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: distribcell_index, offset, i - - distribcell_index = cells(this % cell) % distribcell_index() - offset = 0 - do i = 1, p % n_coord - if (cells(p % coord(i) % cell + 1) % type() == FILL_UNIVERSE) then - offset = offset + cells(p % coord(i) % cell + 1) & - % offset(distribcell_index) - elseif (cells(p % coord(i) % cell + 1) % type() == FILL_LATTICE) then - if (lattices(p % coord(i + 1) % lattice) % are_valid_indices([& - p % coord(i + 1) % lattice_x, & - p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z])) then - offset = offset + lattices(p % coord(i + 1) % lattice) & - % offset(distribcell_index, & - [p % coord(i + 1) % lattice_x, & - p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z]) - end if - end if - if (this % cell == p % coord(i) % cell + 1) then - call match % bins_push_back(offset + 1) - call match % weights_push_back(ONE) - return - end if - end do - end subroutine get_all_bins_distribcell - - subroutine to_statepoint_distribcell(this, filter_group) - class(DistribcellFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "distribcell") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", cells(this % cell) % id()) - end subroutine to_statepoint_distribcell - subroutine initialize_distribcell(this) class(DistribcellFilter), intent(inout) :: this - - integer :: id - integer :: val - - ! Convert id to index. - id = this % cell - val = cell_dict % get(id) - if (val /= EMPTY) then - this % cell = val - this % n_bins = cells(this % cell) % n_instances() - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if + call this % initialize_cpp_inner() + this % n_bins = this % n_bins_cpp() end subroutine initialize_distribcell - function text_label_distribcell(this, bin) result(label) - class(DistribcellFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = '' - call find_offset(this % cell, bin-1, label) - label = "Distributed Cell " // label - end function text_label_distribcell - -!=============================================================================== -! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and -! a target offset to build a string which is the path from the base universe to -! the target cell with the given offset -!=============================================================================== - - subroutine find_offset(i_cell, target_offset, path) - integer, intent(in) :: i_cell ! The target cell index - integer, intent(in) :: target_offset ! Target offset - character(*), intent(inout) :: path ! Path to offset - - integer :: map ! Index in maps vector - integer :: i ! Index over cells - - integer(C_INT) :: path_len - character(kind=C_CHAR), allocatable, target :: path_c(:) - - interface - function distribcell_path_len(target_cell, map, target_offset, root_univ)& - bind(C) result(len) - import C_INT32_T, C_INT - integer(C_INT32_T), intent(in), value :: target_cell, map, & - target_offset, root_univ - integer(C_INT) :: len - end function distribcell_path_len - - subroutine distribcell_path(target_cell, map, target_offset, root_univ, & - path) bind(C) - import C_INT32_T, C_CHAR - integer(C_INT32_T), intent(in), value :: target_cell, map, & - target_offset, root_univ - character(kind=C_CHAR), intent(out) :: path(*) - end subroutine distribcell_path - end interface - - ! Get the distribcell index for this cell - map = cells(i_cell) % distribcell_index() - - path_len = distribcell_path_len(i_cell-1, map, target_offset, & - root_universe) - allocate(path_c(path_len)) - call distribcell_path(i_cell-1, map, target_offset, root_universe, & - path_c) - do i = 1, min(path_len, MAX_LINE_LEN) - if (path_c(i) == C_NULL_CHAR) exit - path(i:i) = path_c(i) - end do - end subroutine find_offset - end module tally_filter_distribcell