From 46ea32139b092d1e164a7a40679ef439b61892f9 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 14 Oct 2018 11:56:56 -0400 Subject: [PATCH] Move MaterialFilter to C++ --- include/openmc/material.h | 2 +- .../openmc/tallies/tally_filter_material.h | 81 +++++++++++ src/cell.cpp | 2 +- src/material.cpp | 10 +- src/tallies/tally_filter.cpp | 21 +++ src/tallies/tally_filter_material.F90 | 132 +++--------------- 6 files changed, 130 insertions(+), 118 deletions(-) create mode 100644 include/openmc/tallies/tally_filter_material.h diff --git a/include/openmc/material.h b/include/openmc/material.h index 1af6f85976..69647dec3e 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -24,7 +24,7 @@ extern std::unordered_map material_map; class Material { public: - int32_t id; //!< Unique ID + int32_t id_; //!< Unique ID double volume_ {-1.0}; //!< Volume in [cm^3] //! \brief Default temperature for cells containing this material. diff --git a/include/openmc/tallies/tally_filter_material.h b/include/openmc/tallies/tally_filter_material.h new file mode 100644 index 0000000000..3fd5893759 --- /dev/null +++ b/include/openmc/tallies/tally_filter_material.h @@ -0,0 +1,81 @@ +#ifndef OPENMC_TALLY_FILTER_MATERIAL_H +#define OPENMC_TALLY_FILTER_MATERIAL_H + +#include +#include +#include + +#include "openmc/error.h" +#include "openmc/material.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +class MaterialFilter : public TallyFilter +{ +public: + virtual std::string type() const override {return "material";} + + virtual ~MaterialFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + materials_ = get_node_array(node, "bins"); + n_bins_ = materials_.size(); + } + + virtual void + initialize() override + { + for (auto& m : materials_) { + auto search = material_map.find(m); + if (search != material_map.end()) { + m = search->second; + } else { + std::stringstream err_msg; + err_msg << "Could not find material " << m + << " specified on tally filter."; + fatal_error(err_msg); + } + } + + for (int i = 0; i < materials_.size(); i++) { + map_[materials_[i]] = i; + } + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + auto search = map_.find(p->material - 1); + if (search != map_.end()) { + // TODO: off-by-one + match.bins.push_back(search->second + 1); + match.weights.push_back(1); + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + TallyFilter::to_statepoint(filter_group); + std::vector material_ids; + for (auto c : materials_) material_ids.push_back(materials[c]->id_); + write_dataset(filter_group, "bins", material_ids); + } + + virtual std::string + text_label(int bin) const override + { + return "Material " + std::to_string(materials[materials_[bin-1]]->id_); + } + + std::vector materials_; + std::unordered_map map_; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_MATERIAL_H diff --git a/src/cell.cpp b/src/cell.cpp index 545bb7f58f..8896cde988 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -476,7 +476,7 @@ CSGCell::to_hdf5(hid_t cell_group) const std::vector mat_ids; for (auto i_mat : material_) { if (i_mat != MATERIAL_VOID) { - mat_ids.push_back(materials[i_mat]->id); + mat_ids.push_back(materials[i_mat]->id_); } else { mat_ids.push_back(MATERIAL_VOID); } diff --git a/src/material.cpp b/src/material.cpp index 1a0c2d2078..fb52d5cd94 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -23,7 +23,7 @@ std::unordered_map material_map; Material::Material(pugi::xml_node node) { if (check_for_node(node, "id")) { - id = std::stoi(get_node_value(node, "id")); + id_ = std::stoi(get_node_value(node, "id")); } else { fatal_error("Must specify id of material in materials XML file."); } @@ -52,7 +52,7 @@ read_materials(pugi::xml_node* node) // Populate the material map. for (int i = 0; i < materials.size(); i++) { - int32_t mid = materials[i]->id; + int32_t mid = materials[i]->id_; auto search = material_map.find(mid); if (search == material_map.end()) { material_map[mid] = i; @@ -78,7 +78,7 @@ openmc_material_get_volume(int32_t index, double* volume) return 0; } else { std::stringstream msg; - msg << "Volume for material with ID=" << m->id << " not set."; + msg << "Volume for material with ID=" << m->id_ << " not set."; set_errmsg(msg); return OPENMC_E_UNASSIGNED; } @@ -113,11 +113,11 @@ openmc_material_set_volume(int32_t index, double volume) extern "C" { Material* material_pointer(int32_t indx) {return materials[indx];} - int32_t material_id(Material* mat) {return mat->id;} + int32_t material_id(Material* mat) {return mat->id_;} void material_set_id(Material* mat, int32_t id, int32_t index) { - mat->id = id; + mat->id_ = id; //TODO: off-by-one material_map[id] = index - 1; } diff --git a/src/tallies/tally_filter.cpp b/src/tallies/tally_filter.cpp index 05de1b46c1..dfe241cd54 100644 --- a/src/tallies/tally_filter.cpp +++ b/src/tallies/tally_filter.cpp @@ -8,6 +8,7 @@ #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_material.h" #include "openmc/tallies/tally_filter_mesh.h" #include "openmc/tallies/tally_filter_meshsurface.h" #include "openmc/tallies/tally_filter_mu.h" @@ -93,6 +94,8 @@ extern "C" { tally_filters.push_back(new CellFromFilter()); } else if (type_ == "distribcell") { tally_filters.push_back(new DistribcellFilter()); + } else if (type_ == "material") { + tally_filters.push_back(new MaterialFilter()); } else if (type_ == "mesh") { tally_filters.push_back(new MeshFilter()); } else if (type_ == "meshsurface") { @@ -140,6 +143,24 @@ extern "C" { *n = filt->cells_.size(); } + void + material_filter_get_bins(MaterialFilter* filt, int32_t** bins, int32_t* n) + { + *bins = filt->materials_.data(); + *n = filt->materials_.size(); + } + + void + material_filter_set_bins(MaterialFilter* filt, int32_t n, int32_t* bins) + { + filt->materials_.clear(); + filt->materials_.resize(n); + for (int i = 0; i < n; i++) filt->materials_[i] = bins[i]; + filt->n_bins_ = filt->materials_.size(); + filt->map_.clear(); + for (int i = 0; i < n; i++) filt->map_[bins[i]] = i; + } + int mesh_filter_get_mesh(MeshFilter* filt) {return filt->mesh_;} void diff --git a/src/tallies/tally_filter_material.F90 b/src/tallies/tally_filter_material.F90 index 2b9a793310..9b1b83d2d9 100644 --- a/src/tallies/tally_filter_material.F90 +++ b/src/tallies/tally_filter_material.F90 @@ -2,117 +2,20 @@ module tally_filter_material use, intrinsic :: ISO_C_BINDING - use constants - use dict_header, only: DictIntInt, EMPTY use error - use hdf5_interface - use material_header, only: materials, material_dict - use particle_header, only: Particle - use string, only: to_str use tally_filter_header - use xml_interface implicit none - private - public :: openmc_material_filter_get_bins - public :: openmc_material_filter_set_bins !=============================================================================== ! MATERIAL specifies which material tally events reside in. !=============================================================================== - type, public, extends(TallyFilter) :: MaterialFilter - integer, allocatable :: materials(:) - type(DictIntInt) :: map - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_material - procedure :: to_statepoint => to_statepoint_material - procedure :: text_label => text_label_material - procedure :: initialize => initialize_material + type, extends(CppTallyFilter) :: MaterialFilter end type MaterialFilter contains - subroutine from_xml(this, node) - class(MaterialFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - allocate(this % materials(n)) - call get_node_array(node, "bins", this % materials) - end subroutine from_xml - - subroutine get_all_bins_material(this, p, estimator, match) - class(MaterialFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: val - - val = this % map % get(p % material) - if (val /= EMPTY) then - call match % bins_push_back(val) - call match % weights_push_back(ONE) - end if - - end subroutine get_all_bins_material - - subroutine to_statepoint_material(this, filter_group) - class(MaterialFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: material_ids(:) - - call write_dataset(filter_group, "type", "material") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(material_ids(size(this % materials))) - do i = 1, size(this % materials) - material_ids(i) = materials(this % materials(i)) % id() - end do - call write_dataset(filter_group, "bins", material_ids) - end subroutine to_statepoint_material - - subroutine initialize_material(this) - class(MaterialFilter), intent(inout) :: this - - integer :: i, id - integer :: val - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % materials(i) - val = material_dict % get(id) - if (val /= EMPTY) then - this % materials(i) = val - else - call fatal_error("Could not find material " // trim(to_str(id)) & - &// " specified on a tally filter.") - end if - end do - - ! Generate mapping from material indices to filter bins. - do i = 1, this % n_bins - call this % map % set(this % materials(i), i) - end do - end subroutine initialize_material - - function text_label_material(this, bin) result(label) - class(MaterialFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Material " // to_str(materials(this % materials(bin)) % id()) - end function text_label_material - !=============================================================================== ! C API FUNCTIONS !=============================================================================== @@ -124,12 +27,20 @@ contains integer(C_INT32_T), intent(out) :: n integer(C_INT) :: err + interface + subroutine material_filter_get_bins(filt, bins, n) bind(C) + import C_PTR, C_INT32_T + type(C_PTR), value :: filt + type(C_PTR) :: bins + integer(C_INT32_T) :: n + end subroutine material_filter_get_bins + end interface + err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (MaterialFilter) - bins = C_LOC(f % materials) - n = size(f % materials) + call material_filter_get_bins(f % ptr, bins, n) err = 0 class default err = E_INVALID_TYPE @@ -147,22 +58,21 @@ contains integer(C_INT32_T), intent(in) :: bins(n) integer(C_INT) :: err - integer :: i + interface + subroutine material_filter_set_bins(filt, n, bins) bind(C) + import C_PTR, C_INT32_T + type(C_PTR), value :: filt + integer(C_INT32_T), value :: n + integer(C_INT32_T) :: bins(n) + end subroutine material_filter_set_bins + end interface err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (MaterialFilter) - f % n_bins = n - if (allocated(f % materials)) deallocate(f % materials) - allocate(f % materials(n)) - f % materials(:) = bins - - ! Generate mapping from material indices to filter bins. - call f % map % clear() - do i = 1, n - call f % map % set(f % materials(i), i) - end do + call material_filter_set_bins(f % ptr, n, bins) + f % n_bins = f % n_bins_cpp() class default err = E_INVALID_TYPE