From e6d2648dc19dcc90861e34f0d4547dc4ff57a63d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 5 Nov 2019 08:13:25 -0500 Subject: [PATCH] Implement cell instance filter --- CMakeLists.txt | 1 + include/openmc/tallies/filter_cell_instance.h | 92 ++++++++++++++++++ src/geometry.cpp | 8 +- src/geometry_aux.cpp | 10 +- src/tallies/filter.cpp | 3 + src/tallies/filter_cell_instance.cpp | 94 +++++++++++++++++++ 6 files changed, 202 insertions(+), 6 deletions(-) create mode 100644 include/openmc/tallies/filter_cell_instance.h create mode 100644 src/tallies/filter_cell_instance.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ccc457f7..021b274c2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -253,6 +253,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_cellborn.cpp src/tallies/filter_cellfrom.cpp src/tallies/filter_cell.cpp + src/tallies/filter_cell_instance.cpp src/tallies/filter_delayedgroup.cpp src/tallies/filter_distribcell.cpp src/tallies/filter_energyfunc.cpp diff --git a/include/openmc/tallies/filter_cell_instance.h b/include/openmc/tallies/filter_cell_instance.h new file mode 100644 index 000000000..d9e3785dc --- /dev/null +++ b/include/openmc/tallies/filter_cell_instance.h @@ -0,0 +1,92 @@ +#ifndef OPENMC_TALLIES_FILTER_CELL_INSTANCE_H +#define OPENMC_TALLIES_FILTER_CELL_INSTANCE_H + +#include +#include // for hash +#include +#include + +#include + +#include "openmc/tallies/filter.h" + +//============================================================================== +//! Define an instance of a particular cell +//============================================================================== + +namespace openmc { +struct CellInstance { + //! Check for equality + bool operator==(const CellInstance& other) const + { return index_cell == other.index_cell && instance == other.instance; } + + gsl::index index_cell; + gsl::index instance; +}; +} // namespace openmc + +namespace std { +template<> +struct hash { + // Taken from https://stackoverflow.com/a/17017281 + std::size_t operator()(const openmc::CellInstance& k) const + { + std::size_t res = 17; + res = 31 * res + std::hash()(k.index_cell); + res = 31 * res + std::hash()(k.instance); + return res; + } +}; +} // namespace std + + +namespace openmc { + +//============================================================================== +//! Specifies cell instances that tally events reside in. +//============================================================================== + +class CellInstanceFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + //CellInstanceFilter() = default; + //CellInstanceFilter(gsl::span instances); + ~CellInstanceFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type() const override {return "cellinstance";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //---------------------------------------------------------------------------- + // Accessors + + const std::vector& cell_instances() const { return cell_instances_; } + + void set_cell_instances(gsl::span instances); + +private: + //---------------------------------------------------------------------------- + // Data members + + //! The indices of the cells binned by this filter. + std::vector cell_instances_; + + //! A map from cell/instance indices to filter bin indices. + std::unordered_map map_; +}; + +} // namespace openmc + +#endif // OPENMC_TALLIES_FILTER_CELL_INSTANCE_H diff --git a/src/geometry.cpp b/src/geometry.cpp index 88b48c866..6bcfc9c58 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -132,8 +132,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) //! Found a material cell which means this is the lowest coord level. // Find the distribcell instance number. - if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) { - int offset = 0; + int offset = 0; + if (c.distribcell_index_ >= 0) { for (int i = 0; i < p->n_coord_; i++) { const auto& c_i {*model::cells[p->coord_[i].cell]}; if (c_i.type_ == FILL_UNIVERSE) { @@ -148,10 +148,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) } } } - p->cell_instance_ = offset; - } else { - p->cell_instance_ = 0; } + p->cell_instance_ = offset; // Set the material and temperature. p->material_last_ = p->material_; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 5ebee4aca..965a58aea 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -18,6 +18,7 @@ #include "openmc/settings.h" #include "openmc/surface.h" #include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_cell_instance.h" #include "openmc/tallies/filter_distribcell.h" @@ -319,13 +320,20 @@ find_root_universe() void prepare_distribcell() { - // Find all cells listed in a DistribcellFilter. + // Find all cells listed in a DistribcellFilter or CellInstanceFilter std::unordered_set distribcells; for (auto& filt : model::tally_filters) { auto* distrib_filt = dynamic_cast(filt.get()); if (distrib_filt) { distribcells.insert(distrib_filt->cell()); } + + auto* ci_filt = dynamic_cast(filt.get()); + if (ci_filt) { + for (const auto& i : ci_filt->cell_instances()) { + distribcells.insert(i.index_cell); + } + } } // Find all cells with distributed materials or temperatures. Make sure that diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index d8fe58cc7..b1b6f4d61 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -12,6 +12,7 @@ #include "openmc/tallies/filter_cell.h" #include "openmc/tallies/filter_cellborn.h" #include "openmc/tallies/filter_cellfrom.h" +#include "openmc/tallies/filter_cell_instance.h" #include "openmc/tallies/filter_delayedgroup.h" #include "openmc/tallies/filter_distribcell.h" #include "openmc/tallies/filter_energyfunc.h" @@ -100,6 +101,8 @@ Filter* Filter::create(const std::string& type, int32_t id) model::tally_filters.push_back(std::make_unique()); } else if (type == "cellfrom") { model::tally_filters.push_back(std::make_unique()); + } else if (type == "cellinstance") { + model::tally_filters.push_back(std::make_unique()); } else if (type == "distribcell") { model::tally_filters.push_back(std::make_unique()); } else if (type == "delayedgroup") { diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp new file mode 100644 index 000000000..427ee1d5d --- /dev/null +++ b/src/tallies/filter_cell_instance.cpp @@ -0,0 +1,94 @@ +#include "openmc/tallies/filter_cell_instance.h" + +#include + +#include "openmc/capi.h" +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +CellInstanceFilter::from_xml(pugi::xml_node node) +{ + // Get cell IDs/instances + auto cells = get_node_array(node, "bins"); + Expects(cells.size() % 2 == 0); + + // Convert into vector of CellInstance + std::vector instances; + for (gsl::index i = 0; i < cells.size() / 2; ++i) { + int32_t cell_id = cells[2*i]; + gsl::index instance = cells[2*i + 1]; + auto search = model::cell_map.find(cell_id); + if (search == model::cell_map.end()) { + std::stringstream err_msg; + err_msg << "Could not find cell " << cell_id + << " specified on tally filter."; + throw std::runtime_error{err_msg.str()}; + } + gsl::index index = search->second; + instances.push_back({index, instance}); + } + + this->set_cell_instances(instances); +} + +void +CellInstanceFilter::set_cell_instances(gsl::span instances) +{ + // Clear existing cells + cell_instances_.clear(); + cell_instances_.reserve(instances.size()); + map_.clear(); + + // Update cells and mapping + for (auto& x : instances) { + Expects(x.index_cell >= 0); + Expects(x.index_cell < model::cells.size()); + cell_instances_.push_back(x); + map_[x] = cell_instances_.size() - 1; + } + + n_bins_ = cell_instances_.size(); +} + +void +CellInstanceFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + gsl::index index_cell = p->coord_[p->n_coord_ - 1].cell; + gsl::index instance = p->cell_instance_; + auto search = map_.find({index_cell, instance}); + if (search != map_.end()) { + int index_bin = search->second; + match.bins_.push_back(index_bin); + match.weights_.push_back(1.0); + } +} + +void +CellInstanceFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + size_t n = cell_instances_.size(); + xt::xtensor data({n, 2}); + for (gsl::index i = 0; i < n; ++i) { + const auto& x = cell_instances_[i]; + data(i, 0) = model::cells[x.index_cell]->id_; + data(i, 1) = x.instance; + } + write_dataset(filter_group, "bins", data); +} + +std::string +CellInstanceFilter::text_label(int bin) const +{ + const auto& x = cell_instances_[bin]; + auto cell_id = model::cells[x.index_cell]->id_; + return "Cell " + std::to_string(cell_id) + ", CellInstance " + + std::to_string(x.instance); +} + +} // namespace openmc