diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 9f05a56dcc..148a28c888 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -10,15 +10,15 @@ extern "C" { #endif int openmc_calculate_volumes(); - int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n); + int openmc_cell_filter_get_bins(int32_t index, const int32_t** cells, int32_t* n); int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n); int openmc_cell_get_id(int32_t index, int32_t* id); int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T); int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices); int openmc_cell_set_id(int32_t index, int32_t id); int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance); - int openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n); - int openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies); + int openmc_energy_filter_get_bins(int32_t index, const double** energies, size_t* n); + int openmc_energy_filter_set_bins(int32_t index, size_t n, const double* energies); int openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end); int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end); @@ -48,7 +48,7 @@ extern "C" { int openmc_legendre_filter_set_order(int32_t index, int order); int openmc_load_nuclide(const char* name); int openmc_material_add_nuclide(int32_t index, const char name[], double density); - int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n); + int openmc_material_get_densities(int32_t index, const int** nuclides, const double** densities, int* n); int openmc_material_get_id(int32_t index, int32_t* id); int openmc_material_get_fissionable(int32_t index, bool* fissionable); int openmc_material_get_density(int32_t index, double* density); @@ -57,8 +57,8 @@ extern "C" { int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density); int openmc_material_set_id(int32_t index, int32_t id); int openmc_material_set_volume(int32_t index, double volume); - int openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n); - int openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins); + int openmc_material_filter_get_bins(int32_t index, const int32_t** bins, size_t* n); + int openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins); int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh); int openmc_mesh_get_id(int32_t index, int32_t* id); @@ -95,7 +95,7 @@ extern "C" { int openmc_tally_get_active(int32_t index, bool* active); int openmc_tally_get_estimator(int32_t index, int* estimator); int openmc_tally_get_id(int32_t index, int32_t* id); - int openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n); + int openmc_tally_get_filters(int32_t index, const int32_t** indices, size_t* n); int openmc_tally_get_n_realizations(int32_t index, int32_t* n); int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n); int openmc_tally_get_scores(int32_t index, int** scores, int* n); @@ -104,7 +104,7 @@ extern "C" { int openmc_tally_results(int32_t index, double** ptr, size_t shape_[3]); int openmc_tally_set_active(int32_t index, bool active); int openmc_tally_set_estimator(int32_t index, const char* estimator); - int openmc_tally_set_filters(int32_t index, int n, const int32_t* indices); + int openmc_tally_set_filters(int32_t index, size_t n, const int32_t* indices); int openmc_tally_set_id(int32_t index, int32_t id); int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides); int openmc_tally_set_scores(int32_t index, int n, const char** scores); diff --git a/include/openmc/cell.h b/include/openmc/cell.h index ee12305951..afe40fe433 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -72,9 +72,66 @@ public: //! A geometry primitive that links surfaces, universes, and materials //============================================================================== -class Cell -{ +class Cell { public: + //---------------------------------------------------------------------------- + // Constructors, destructors, factory functions + + explicit Cell(pugi::xml_node cell_node); + Cell() {}; + virtual ~Cell() = default; + + //---------------------------------------------------------------------------- + // Methods + + //! \brief Determine if a cell contains the particle at a given location. + //! + //! The bounds of the cell are detemined by a logical expression involving + //! surface half-spaces. At initialization, the expression was converted + //! to RPN notation. + //! + //! The function is split into two cases, one for simple cells (those + //! involving only the intersection of half-spaces) and one for complex cells. + //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon + //! as we know that one half-space is not satisfied, we can exit. This + //! provides a performance benefit for the common case. In + //! contains_complex, we evaluate the RPN expression using a stack, similar to + //! how a RPN calculator would work. + //! \param r The 3D Cartesian coordinate to check. + //! \param u A direction used to "break ties" the coordinates are very + //! close to a surface. + //! \param on_surface The signed index of a surface that the coordinate is + //! known to be on. This index takes precedence over surface sense + //! calculations. + virtual bool + contains(Position r, Direction u, int32_t on_surface) const = 0; + + //! Find the oncoming boundary of this cell. + virtual std::pair + distance(Position r, Direction u, int32_t on_surface) const = 0; + + //! Write all information needed to reconstruct the cell to an HDF5 group. + //! \param group_id An HDF5 group id. + virtual void to_hdf5(hid_t group_id) const = 0; + + //---------------------------------------------------------------------------- + // Accessors + + //! Get the temperature of a cell instance + //! \param[in] instance Instance index. If -1 is given, the temperature for + //! the first instance is returned. + //! \return Temperature in [K] + double temperature(int32_t instance = -1) const; + + //! Set the temperature of a cell instance + //! \param[in] T Temperature in [K] + //! \param[in] instance Instance index. If -1 is given, the temperature for + //! all instances is set. + void set_temperature(double T, int32_t instance = -1); + + //---------------------------------------------------------------------------- + // Data members + int32_t id_; //!< Unique ID std::string name_; //!< User-defined name int type_; //!< Material, universe, or lattice @@ -116,41 +173,6 @@ public: std::vector rotation_; std::vector offset_; //!< Distribcell offset table - - explicit Cell(pugi::xml_node cell_node); - Cell() {}; - - //! \brief Determine if a cell contains the particle at a given location. - //! - //! The bounds of the cell are detemined by a logical expression involving - //! surface half-spaces. At initialization, the expression was converted - //! to RPN notation. - //! - //! The function is split into two cases, one for simple cells (those - //! involving only the intersection of half-spaces) and one for complex cells. - //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon - //! as we know that one half-space is not satisfied, we can exit. This - //! provides a performance benefit for the common case. In - //! contains_complex, we evaluate the RPN expression using a stack, similar to - //! how a RPN calculator would work. - //! \param r The 3D Cartesian coordinate to check. - //! \param u A direction used to "break ties" the coordinates are very - //! close to a surface. - //! \param on_surface The signed index of a surface that the coordinate is - //! known to be on. This index takes precedence over surface sense - //! calculations. - virtual bool - contains(Position r, Direction u, int32_t on_surface) const = 0; - - //! Find the oncoming boundary of this cell. - virtual std::pair - distance(Position r, Direction u, int32_t on_surface) const = 0; - - //! Write all information needed to reconstruct the cell to an HDF5 group. - //! @param group_id An HDF5 group id. - virtual void to_hdf5(hid_t group_id) const = 0; - - virtual ~Cell() {} }; //============================================================================== diff --git a/include/openmc/material.h b/include/openmc/material.h index 3109f6d787..02c5380f35 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -6,6 +6,7 @@ #include #include +#include #include #include "pugixml.hpp" #include "xtensor/xtensor.hpp" @@ -35,6 +36,7 @@ extern std::unordered_map material_map; class Material { public: + //---------------------------------------------------------------------------- // Types struct ThermalTable { int index_table; //!< Index of table in data::thermal_scatt @@ -42,11 +44,15 @@ public: double fraction; //!< How often to use table }; - // Constructors + //---------------------------------------------------------------------------- + // Constructors, destructors, factory functions Material() {}; explicit Material(pugi::xml_node material_node); + ~Material(); + //---------------------------------------------------------------------------- // Methods + void calculate_xs(Particle& p) const; //! Assign thermal scattering tables to specific nuclides within the material @@ -59,12 +65,65 @@ public: //! Finalize the material, assigning tables, normalize density, etc. void finalize(); - //! Set total density of the material - int set_density(double density, std::string units); - //! Write material data to HDF5 void to_hdf5(hid_t group) const; + //! Add nuclide to the material + // + //! \param[in] nuclide Name of the nuclide + //! \param[in] density Density of the nuclide in [atom/b-cm] + void add_nuclide(const std::string& nuclide, double density); + + //! Set atom densities for the material + // + //! \param[in] name Name of each nuclide + //! \param[in] density Density of each nuclide in [atom/b-cm] + void set_densities(const std::vector& name, + const std::vector& density); + + //---------------------------------------------------------------------------- + // Accessors + + //! Get density in [atom/b-cm] + //! \return Density in [atom/b-cm] + double density() const { return density_; } + + //! Get density in [g/cm^3] + //! \return Density in [g/cm^3] + double density_gpcc() const { return density_gpcc_; } + + //! Set total density of the material + // + //! \param[in] density Density value + //! \param[in] units Units of density + void set_density(double density, gsl::cstring_span units); + + //! Get nuclides in material + //! \return Indices into the global nuclides vector + gsl::span nuclides() const { return {nuclide_.data(), nuclide_.size()}; } + + //! Get densities of each nuclide in material + //! \return Densities in [atom/b-cm] + gsl::span densities() const { return {atom_density_.data(), atom_density_.size()}; } + + //! Get ID of material + //! \return ID of material + int32_t id() const { return id_; } + + //! Assign a unique ID to the material + //! \param[in] Unique ID to assign. A value of -1 indicates that an ID + //! should be automatically assigned. + void set_id(int32_t id); + + //! Get whether material is fissionable + //! \return Whether material is fissionable + bool fissionable() const { return fissionable_; } + + //! Get volume of material + //! \return Volume in [cm^3] + double volume() const; + + //---------------------------------------------------------------------------- // Data int32_t id_; //!< Unique ID std::string name_; //!< Name of material @@ -95,6 +154,9 @@ public: std::unique_ptr ttb_; private: + //---------------------------------------------------------------------------- + // Private methods + //! Calculate the collision stopping power void collision_stopping_power(double* s_col, bool positron); @@ -106,6 +168,10 @@ private: void calculate_neutron_xs(Particle& p) const; void calculate_photon_xs(Particle& p) const; + + //---------------------------------------------------------------------------- + // Private data members + gsl::index index_; }; //============================================================================== diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index c2fd4dfb37..ba6ca08830 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -7,6 +7,8 @@ #include #include +#include + #include "openmc/hdf5_interface.h" #include "openmc/particle.h" #include "pugixml.hpp" @@ -43,13 +45,34 @@ namespace openmc { class Filter { public: - virtual ~Filter() = default; + //---------------------------------------------------------------------------- + // Constructors, destructors, factory functions - virtual std::string type() const = 0; + Filter(); + virtual ~Filter(); + + //! Create a new tally filter + // + //! \param[in] type Type of the filter + //! \param[in] id Unique ID for the filter. If none is passed, an ID is + //! automatically assigned + //! \return Pointer to the new filter object + static Filter* create(const std::string& type, int32_t id = -1); + + //! Create a new tally filter from an XML node + // + //! \param[in] node XML node + //! \return Pointer to the new filter object + static Filter* create(pugi::xml_node node); //! Uses an XML input to fill the filter's data fields. virtual void from_xml(pugi::xml_node node) = 0; + //---------------------------------------------------------------------------- + // Methods + + virtual std::string type() const = 0; + //! Matches a tally event to a set of filter bins and weights. //! //! \param[out] match will contain the matching bins and corresponding @@ -71,11 +94,32 @@ public: //! "Incoming Energy [0.625E-6, 20.0)". virtual std::string text_label(int bin) const = 0; - virtual void initialize() {} + //---------------------------------------------------------------------------- + // Accessors - int32_t id_; + //! Get unique ID of filter + //! \return Unique ID + int32_t id() const { return id_; } + //! Assign a unique ID to the filter + //! \param[in] Unique ID to assign. A value of -1 indicates that an ID should + //! be automatically assigned + void set_id(int32_t id); + + //! Get number of bins + //! \return Number of bins + int n_bins() const { return n_bins_; } + + gsl::index index() const { return index_; } + + //---------------------------------------------------------------------------- + // Data members + +protected: int n_bins_; +private: + int32_t id_ {-1}; + gsl::index index_; }; //============================================================================== @@ -99,8 +143,6 @@ namespace model { // Non-member functions //============================================================================== -Filter* allocate_filter(const std::string& type); - //! Make sure index corresponds to a valid filter int verify_filter(int32_t index); diff --git a/include/openmc/tallies/filter_azimuthal.h b/include/openmc/tallies/filter_azimuthal.h index b896400315..84c2a2d8c4 100644 --- a/include/openmc/tallies/filter_azimuthal.h +++ b/include/openmc/tallies/filter_azimuthal.h @@ -4,6 +4,8 @@ #include #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -15,8 +17,14 @@ namespace openmc { class AzimuthalFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~AzimuthalFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "azimuthal";} void from_xml(pugi::xml_node node) override; @@ -28,6 +36,15 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + void set_bins(gsl::span bins); + +private: + //---------------------------------------------------------------------------- + // Data members + std::vector bins_; }; diff --git a/include/openmc/tallies/filter_cell.h b/include/openmc/tallies/filter_cell.h index 9c654e82ed..b570ee0273 100644 --- a/include/openmc/tallies/filter_cell.h +++ b/include/openmc/tallies/filter_cell.h @@ -5,6 +5,8 @@ #include #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -16,14 +18,18 @@ namespace openmc { class CellFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~CellFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "cell";} void from_xml(pugi::xml_node node) override; - void initialize() override; - void get_all_bins(const Particle* p, int estimator, FilterMatch& match) const override; @@ -31,6 +37,17 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + const std::vector& cells() const { return cells_; } + + void set_cells(gsl::span cells); + +protected: + //---------------------------------------------------------------------------- + // Data members + //! The indices of the cells binned by this filter. std::vector cells_; diff --git a/include/openmc/tallies/filter_cellborn.h b/include/openmc/tallies/filter_cellborn.h index 400e82c282..706f6c47a4 100644 --- a/include/openmc/tallies/filter_cellborn.h +++ b/include/openmc/tallies/filter_cellborn.h @@ -14,6 +14,9 @@ namespace openmc { class CellbornFilter : public CellFilter { public: + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "cellborn";} void get_all_bins(const Particle* p, int estimator, FilterMatch& match) diff --git a/include/openmc/tallies/filter_cellfrom.h b/include/openmc/tallies/filter_cellfrom.h index bd3e08dcbb..e86e34854c 100644 --- a/include/openmc/tallies/filter_cellfrom.h +++ b/include/openmc/tallies/filter_cellfrom.h @@ -14,6 +14,9 @@ namespace openmc { class CellFromFilter : public CellFilter { public: + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "cellfrom";} void get_all_bins(const Particle* p, int estimator, FilterMatch& match) diff --git a/include/openmc/tallies/filter_delayedgroup.h b/include/openmc/tallies/filter_delayedgroup.h index 2f65a50593..8a2bbeeaa2 100644 --- a/include/openmc/tallies/filter_delayedgroup.h +++ b/include/openmc/tallies/filter_delayedgroup.h @@ -3,6 +3,8 @@ #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -17,8 +19,14 @@ namespace openmc { class DelayedGroupFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~DelayedGroupFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "delayedgroup";} void from_xml(pugi::xml_node node) override; @@ -30,6 +38,17 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + const std::vector& groups() const { return groups_; } + + void set_groups(gsl::span groups); + +private: + //---------------------------------------------------------------------------- + // Data members + std::vector groups_; }; diff --git a/include/openmc/tallies/filter_distribcell.h b/include/openmc/tallies/filter_distribcell.h index 73857e7cb0..9430a09061 100644 --- a/include/openmc/tallies/filter_distribcell.h +++ b/include/openmc/tallies/filter_distribcell.h @@ -14,14 +14,18 @@ namespace openmc { class DistribcellFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~DistribcellFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "distribcell";} void from_xml(pugi::xml_node node) override; - void initialize() override; - void get_all_bins(const Particle* p, int estimator, FilterMatch& match) const override; @@ -29,6 +33,17 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + int32_t cell() const { return cell_; } + + void set_cell(int32_t cell); + +private: + //---------------------------------------------------------------------------- + // Data members + int32_t cell_; }; diff --git a/include/openmc/tallies/filter_energy.h b/include/openmc/tallies/filter_energy.h index 025a77c624..af7601721e 100644 --- a/include/openmc/tallies/filter_energy.h +++ b/include/openmc/tallies/filter_energy.h @@ -3,6 +3,8 @@ #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -14,8 +16,14 @@ namespace openmc { class EnergyFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~EnergyFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "energy";} void from_xml(pugi::xml_node node) override; @@ -27,6 +35,18 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + const std::vector& bins() const { return bins_; } + void set_bins(gsl::span bins); + + bool matches_transport_groups() const { return matches_transport_groups_; } + +protected: + //---------------------------------------------------------------------------- + // Data members + std::vector bins_; //! True if transport group number can be used directly to get bin number @@ -43,6 +63,9 @@ public: class EnergyoutFilter : public EnergyFilter { public: + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "energyout";} void get_all_bins(const Particle* p, int estimator, FilterMatch& match) diff --git a/include/openmc/tallies/filter_energyfunc.h b/include/openmc/tallies/filter_energyfunc.h index 6f51826267..9b8d839f8b 100644 --- a/include/openmc/tallies/filter_energyfunc.h +++ b/include/openmc/tallies/filter_energyfunc.h @@ -15,6 +15,9 @@ namespace openmc { class EnergyFunctionFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + EnergyFunctionFilter() : Filter {} { @@ -23,6 +26,9 @@ public: ~EnergyFunctionFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "energyfunction";} void from_xml(pugi::xml_node node) override; @@ -34,6 +40,10 @@ public: std::string text_label(int bin) const override; +private: + //---------------------------------------------------------------------------- + // Data members + //! Incident neutron energy interpolation grid. std::vector energy_; diff --git a/include/openmc/tallies/filter_legendre.h b/include/openmc/tallies/filter_legendre.h index 054ba14e79..3a14ec3cf9 100644 --- a/include/openmc/tallies/filter_legendre.h +++ b/include/openmc/tallies/filter_legendre.h @@ -14,8 +14,14 @@ namespace openmc { class LegendreFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~LegendreFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "legendre";} void from_xml(pugi::xml_node node) override; @@ -27,6 +33,17 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + int order() const { return order_; } + + void set_order(int order); + +private: + //---------------------------------------------------------------------------- + // Data members + int order_; }; diff --git a/include/openmc/tallies/filter_material.h b/include/openmc/tallies/filter_material.h index d272451db5..65cf832a3f 100644 --- a/include/openmc/tallies/filter_material.h +++ b/include/openmc/tallies/filter_material.h @@ -5,6 +5,8 @@ #include #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -16,14 +18,18 @@ namespace openmc { class MaterialFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~MaterialFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "material";} void from_xml(pugi::xml_node node) override; - void initialize() override; - void get_all_bins(const Particle* p, int estimator, FilterMatch& match) const override; @@ -31,6 +37,19 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + std::vector& materials() { return materials_; } + + const std::vector& materials() const { return materials_; } + + void set_materials(gsl::span materials); + +private: + //---------------------------------------------------------------------------- + // Data members + //! The indices of the materials binned by this filter. std::vector materials_; diff --git a/include/openmc/tallies/filter_mesh.h b/include/openmc/tallies/filter_mesh.h index a6caa8f565..98dec5d501 100644 --- a/include/openmc/tallies/filter_mesh.h +++ b/include/openmc/tallies/filter_mesh.h @@ -16,8 +16,14 @@ namespace openmc { class MeshFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~MeshFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "mesh";} void from_xml(pugi::xml_node node) override; @@ -29,11 +35,17 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + virtual int32_t mesh() const {return mesh_;} virtual void set_mesh(int32_t mesh); protected: + //---------------------------------------------------------------------------- + // Data members + int32_t mesh_; }; diff --git a/include/openmc/tallies/filter_meshsurface.h b/include/openmc/tallies/filter_meshsurface.h index 32393cfaca..19d178c7ff 100644 --- a/include/openmc/tallies/filter_meshsurface.h +++ b/include/openmc/tallies/filter_meshsurface.h @@ -8,6 +8,9 @@ namespace openmc { class MeshSurfaceFilter : public MeshFilter { public: + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "meshsurface";} void get_all_bins(const Particle* p, int estimator, FilterMatch& match) @@ -15,6 +18,9 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + void set_mesh(int32_t mesh) override; }; diff --git a/include/openmc/tallies/filter_mu.h b/include/openmc/tallies/filter_mu.h index 68b2f7b023..ae9c3e06ac 100644 --- a/include/openmc/tallies/filter_mu.h +++ b/include/openmc/tallies/filter_mu.h @@ -3,6 +3,8 @@ #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -15,8 +17,14 @@ namespace openmc { class MuFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~MuFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "mu";} void from_xml(pugi::xml_node node) override; @@ -28,6 +36,15 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + void set_bins(gsl::span bins); + +private: + //---------------------------------------------------------------------------- + // Data members + std::vector bins_; }; diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h index 268c19c6de..61aa3106b6 100644 --- a/include/openmc/tallies/filter_particle.h +++ b/include/openmc/tallies/filter_particle.h @@ -15,8 +15,14 @@ namespace openmc { class ParticleFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~ParticleFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "particle";} void from_xml(pugi::xml_node node) override; @@ -28,6 +34,17 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + const std::vector& particles() const { return particles_; } + + void set_particles(gsl::span particles); + +private: + //---------------------------------------------------------------------------- + // Data members + std::vector particles_; }; diff --git a/include/openmc/tallies/filter_polar.h b/include/openmc/tallies/filter_polar.h index 2bcdcd07f4..9659504568 100644 --- a/include/openmc/tallies/filter_polar.h +++ b/include/openmc/tallies/filter_polar.h @@ -4,6 +4,8 @@ #include #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -15,8 +17,14 @@ namespace openmc { class PolarFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~PolarFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "polar";} void from_xml(pugi::xml_node node) override; @@ -28,6 +36,15 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + void set_bins(gsl::span bins); + +private: + //---------------------------------------------------------------------------- + // Data members + std::vector bins_; }; diff --git a/include/openmc/tallies/filter_sph_harm.h b/include/openmc/tallies/filter_sph_harm.h index b53b624417..80f998b2f2 100644 --- a/include/openmc/tallies/filter_sph_harm.h +++ b/include/openmc/tallies/filter_sph_harm.h @@ -3,6 +3,8 @@ #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -18,8 +20,14 @@ enum class SphericalHarmonicsCosine { class SphericalHarmonicsFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~SphericalHarmonicsFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "sphericalharmonics";} void from_xml(pugi::xml_node node) override; @@ -31,6 +39,21 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + int order() const { return order_; } + + void set_order(int order); + + SphericalHarmonicsCosine cosine() const { return cosine_; } + + void set_cosine(gsl::cstring_span cosine); + +private: + //---------------------------------------------------------------------------- + // Data members + int order_; //! The type of angle that this filter measures when binning events. diff --git a/include/openmc/tallies/filter_sptl_legendre.h b/include/openmc/tallies/filter_sptl_legendre.h index 995bc73608..fea7275150 100644 --- a/include/openmc/tallies/filter_sptl_legendre.h +++ b/include/openmc/tallies/filter_sptl_legendre.h @@ -18,8 +18,14 @@ enum class LegendreAxis { class SpatialLegendreFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~SpatialLegendreFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "spatiallegendre";} void from_xml(pugi::xml_node node) override; @@ -31,6 +37,23 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + int order() const { return order_; } + void set_order(int order); + + LegendreAxis axis() const { return axis_; } + void set_axis(LegendreAxis axis); + + double min() const { return min_; } + double max() const { return max_; } + void set_minmax(double min, double max); + +private: + //---------------------------------------------------------------------------- + // Data members + int order_; //! The Cartesian coordinate axis that the Legendre expansion is applied to. diff --git a/include/openmc/tallies/filter_surface.h b/include/openmc/tallies/filter_surface.h index ea4dc5b990..23b9be33dd 100644 --- a/include/openmc/tallies/filter_surface.h +++ b/include/openmc/tallies/filter_surface.h @@ -5,6 +5,8 @@ #include #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -16,14 +18,18 @@ namespace openmc { class SurfaceFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~SurfaceFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "surface";} void from_xml(pugi::xml_node node) override; - void initialize() override; - void get_all_bins(const Particle* p, int estimator, FilterMatch& match) const override; @@ -31,6 +37,15 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + void set_surfaces(gsl::span surfaces); + +private: + //---------------------------------------------------------------------------- + // Data members + //! The indices of the surfaces binned by this filter. std::vector surfaces_; diff --git a/include/openmc/tallies/filter_universe.h b/include/openmc/tallies/filter_universe.h index 3ae0093cbd..fc2ace18f8 100644 --- a/include/openmc/tallies/filter_universe.h +++ b/include/openmc/tallies/filter_universe.h @@ -5,6 +5,8 @@ #include #include +#include + #include "openmc/tallies/filter.h" namespace openmc { @@ -16,14 +18,18 @@ namespace openmc { class UniverseFilter : public Filter { public: + //---------------------------------------------------------------------------- + // Constructors, destructors + ~UniverseFilter() = default; + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "universe";} void from_xml(pugi::xml_node node) override; - void initialize() override; - void get_all_bins(const Particle* p, int estimator, FilterMatch& match) const override; @@ -31,6 +37,15 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + + void set_universes(gsl::span universes); + +private: + //---------------------------------------------------------------------------- + // Data members + //! The indices of the universes binned by this filter. std::vector universes_; diff --git a/include/openmc/tallies/filter_zernike.h b/include/openmc/tallies/filter_zernike.h index 19fff13363..e3ae89dec9 100644 --- a/include/openmc/tallies/filter_zernike.h +++ b/include/openmc/tallies/filter_zernike.h @@ -14,10 +14,16 @@ namespace openmc { class ZernikeFilter : public Filter { public: - std::string type() const override {return "zernike";} + //---------------------------------------------------------------------------- + // Constructors, destructors ~ZernikeFilter() = default; + //---------------------------------------------------------------------------- + // Methods + + std::string type() const override {return "zernike";} + void from_xml(pugi::xml_node node) override; void get_all_bins(const Particle* p, int estimator, FilterMatch& match) @@ -27,10 +33,25 @@ public: std::string text_label(int bin) const override; - int order() const {return order_;} + //---------------------------------------------------------------------------- + // Accessors + int order() const { return order_; } virtual void set_order(int order); + double x() const { return x_; } + void set_x(double x) { x_ = x; } + + double y() const { return y_; } + void set_y(double y) { y_ = y; } + + double r() const { return r_; } + void set_r(double r) { r_ = r; } + + //---------------------------------------------------------------------------- + // Data members + +protected: //! Cartesian x coordinate for the origin of this expansion. double x_; @@ -40,7 +61,6 @@ public: //! Maximum radius from the origin covered by this expansion. double r_; -protected: int order_; }; @@ -51,6 +71,9 @@ protected: class ZernikeRadialFilter : public ZernikeFilter { public: + //---------------------------------------------------------------------------- + // Methods + std::string type() const override {return "zernikeradial";} void get_all_bins(const Particle* p, int estimator, FilterMatch& match) @@ -58,6 +81,9 @@ public: std::string text_label(int bin) const override; + //---------------------------------------------------------------------------- + // Accessors + void set_order(int order) override; }; diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 969c064e75..2a166738d8 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -2,8 +2,10 @@ #define OPENMC_TALLIES_TALLY_H #include "openmc/constants.h" +#include "openmc/tallies/filter.h" #include "openmc/tallies/trigger.h" +#include #include "pugixml.hpp" #include "xtensor/xfixed.hpp" #include "xtensor/xtensor.hpp" @@ -21,24 +23,33 @@ namespace openmc { class Tally { public: - Tally(); + //---------------------------------------------------------------------------- + // Constructors, destructors, factory functions + explicit Tally(int32_t id); + explicit Tally(pugi::xml_node node); + ~Tally(); + static Tally* create(int32_t id = -1); - void init_from_xml(pugi::xml_node node); + //---------------------------------------------------------------------------- + // Accessors + + void set_id(int32_t id); + + void set_active(bool active) { active_ = active; } void set_scores(pugi::xml_node node); - void set_scores(std::vector scores); + void set_scores(const std::vector& scores); void set_nuclides(pugi::xml_node node); - //---------------------------------------------------------------------------- - // Methods for getting and setting filter/stride data. + void set_nuclides(const std::vector& nuclides); const std::vector& filters() const {return filters_;} int32_t filters(int i) const {return filters_[i];} - void set_filters(const int32_t filter_indices[], int n); + void set_filters(gsl::span filters); int32_t strides(int i) const {return strides_[i];} @@ -112,6 +123,8 @@ private: std::vector strides_; int32_t n_filter_bins_ {0}; + + gsl::index index_; }; //============================================================================== diff --git a/openmc/capi/error.py b/openmc/capi/error.py index b35de4e60c..89e7b6e390 100644 --- a/openmc/capi/error.py +++ b/openmc/capi/error.py @@ -36,4 +36,6 @@ def _error_handler(err, func, args): elif err == errcode('OPENMC_E_WARNING'): warn(msg) elif err < 0: - raise exc.OpenMCError("Unknown error encountered (code {}).".format(err)) + if not msg: + msg = "Unknown error encountered (code {}).".format(err) + raise exc.OpenMCError(msg) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index c0716bcef6..bda3c3178f 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -28,10 +28,10 @@ _dll.openmc_cell_filter_get_bins.argtypes = [ _dll.openmc_cell_filter_get_bins.restype = c_int _dll.openmc_cell_filter_get_bins.errcheck = _error_handler _dll.openmc_energy_filter_get_bins.argtypes = [ - c_int32, POINTER(POINTER(c_double)), POINTER(c_int32)] + c_int32, POINTER(POINTER(c_double)), POINTER(c_size_t)] _dll.openmc_energy_filter_get_bins.restype = c_int _dll.openmc_energy_filter_get_bins.errcheck = _error_handler -_dll.openmc_energy_filter_set_bins.argtypes = [c_int32, c_int32, POINTER(c_double)] +_dll.openmc_energy_filter_set_bins.argtypes = [c_int32, c_size_t, POINTER(c_double)] _dll.openmc_energy_filter_set_bins.restype = c_int _dll.openmc_energy_filter_set_bins.errcheck = _error_handler _dll.openmc_filter_get_id.argtypes = [c_int32, POINTER(c_int32)] @@ -53,10 +53,10 @@ _dll.openmc_legendre_filter_set_order.argtypes = [c_int32, c_int] _dll.openmc_legendre_filter_set_order.restype = c_int _dll.openmc_legendre_filter_set_order.errcheck = _error_handler _dll.openmc_material_filter_get_bins.argtypes = [ - c_int32, POINTER(POINTER(c_int32)), POINTER(c_int32)] + c_int32, POINTER(POINTER(c_int32)), POINTER(c_size_t)] _dll.openmc_material_filter_get_bins.restype = c_int _dll.openmc_material_filter_get_bins.errcheck = _error_handler -_dll.openmc_material_filter_set_bins.argtypes = [c_int32, c_int32, POINTER(c_int32)] +_dll.openmc_material_filter_set_bins.argtypes = [c_int32, c_size_t, POINTER(c_int32)] _dll.openmc_material_filter_set_bins.restype = c_int _dll.openmc_material_filter_set_bins.errcheck = _error_handler _dll.openmc_mesh_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)] @@ -94,6 +94,7 @@ _dll.openmc_zernike_filter_set_order.restype = c_int _dll.openmc_zernike_filter_set_order.errcheck = _error_handler _dll.tally_filters_size.restype = c_size_t + class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -148,7 +149,7 @@ class EnergyFilter(Filter): @property def bins(self): energies = POINTER(c_double)() - n = c_int32() + n = c_size_t() _dll.openmc_energy_filter_get_bins(self._index, energies, n) return as_array(energies, (n.value,)) @@ -231,7 +232,7 @@ class MaterialFilter(Filter): @property def bins(self): materials = POINTER(c_int32)() - n = c_int32() + n = c_size_t() _dll.openmc_material_filter_get_bins(self._index, materials, n) return [Material(index=materials[i]) for i in range(n.value)] diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 9529a31f21..88c14e4494 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -36,7 +36,7 @@ _dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_tally_get_id.restype = c_int _dll.openmc_tally_get_id.errcheck = _error_handler _dll.openmc_tally_get_filters.argtypes = [ - c_int32, POINTER(POINTER(c_int32)), POINTER(c_int)] + c_int32, POINTER(POINTER(c_int32)), POINTER(c_size_t)] _dll.openmc_tally_get_filters.restype = c_int _dll.openmc_tally_get_filters.errcheck = _error_handler _dll.openmc_tally_get_n_realizations.argtypes = [c_int32, POINTER(c_int32)] @@ -63,7 +63,7 @@ _dll.openmc_tally_results.errcheck = _error_handler _dll.openmc_tally_set_active.argtypes = [c_int32, c_bool] _dll.openmc_tally_set_active.restype = c_int _dll.openmc_tally_set_active.errcheck = _error_handler -_dll.openmc_tally_set_filters.argtypes = [c_int32, c_int, POINTER(c_int32)] +_dll.openmc_tally_set_filters.argtypes = [c_int32, c_size_t, POINTER(c_int32)] _dll.openmc_tally_set_filters.restype = c_int _dll.openmc_tally_set_filters.errcheck = _error_handler _dll.openmc_tally_set_estimator.argtypes = [c_int32, c_char_p] @@ -249,7 +249,7 @@ class Tally(_FortranObjectWithID): @property def filters(self): filt_idx = POINTER(c_int32)() - n = c_int() + n = c_size_t() _dll.openmc_tally_get_filters(self._index, filt_idx, n) return [_get_filter(filt_idx[i]) for i in range(n.value)] diff --git a/src/cell.cpp b/src/cell.cpp index 87b67ef856..04273eded0 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -212,6 +212,39 @@ Universe::to_hdf5(hid_t universes_group) const // Cell implementation //============================================================================== +double +Cell::temperature(int32_t instance) const +{ + if (sqrtkT_.size() < 1) { + throw std::runtime_error{"Cell temperature has not yet been set."}; + } + + if (instance >= 0) { + double sqrtkT = sqrtkT_.size() == 1 ? + sqrtkT_.at(0) : + sqrtkT_.at(instance); + return sqrtkT * sqrtkT / K_BOLTZMANN; + } else { + return sqrtkT_[0] * sqrtkT_[0] / K_BOLTZMANN; + } +} + +void +Cell::set_temperature(double T, int32_t instance) +{ + if (instance >= 0) { + sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T); + } else { + for (auto& T_ : sqrtkT_) { + T_ = std::sqrt(K_BOLTZMANN * T); + } + } +} + +//============================================================================== +// CSGCell implementation +//============================================================================== + CSGCell::CSGCell() {} // empty constructor CSGCell::CSGCell(pugi::xml_node cell_node) @@ -917,27 +950,18 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n, extern "C" int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) { - if (index >= 0 && index < model::cells.size()) { - Cell& c {*model::cells[index]}; - - if (instance) { - if (*instance >= 0 && *instance < c.sqrtkT_.size()) { - c.sqrtkT_[*instance] = std::sqrt(K_BOLTZMANN * T); - } else { - strcpy(openmc_err_msg, "Distribcell instance is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } - } else { - for (auto& T_ : c.sqrtkT_) { - T_ = std::sqrt(K_BOLTZMANN * T); - } - } - - } else { + if (index < 0 || index >= model::cells.size()) { strcpy(openmc_err_msg, "Index in cells array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } + int32_t instance_index = instance ? *instance : -1; + try { + model::cells[index]->set_temperature(T, instance_index); + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; + } return 0; } @@ -949,25 +973,13 @@ openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) return OPENMC_E_OUT_OF_BOUNDS; } - Cell& c {*model::cells[index]}; - - if (c.sqrtkT_.size() < 1) { - strcpy(openmc_err_msg, "Cell temperature has not yet been set."); + int32_t instance_index = instance ? *instance : -1; + try { + *T = model::cells[index]->temperature(instance_index); + } catch (const std::exception& e) { + set_errmsg(e.what()); return OPENMC_E_UNASSIGNED; } - - if (instance) { - if (*instance >= 0 && *instance < c.n_instances_) { - double sqrtkT = c.sqrtkT_.size() == 1 ? c.sqrtkT_[0] : c.sqrtkT_[*instance]; - *T = sqrtkT * sqrtkT / K_BOLTZMANN; - } else { - strcpy(openmc_err_msg, "Distribcell instance is out of bounds."); - return OPENMC_E_OUT_OF_BOUNDS; - } - } else { - *T = c.sqrtkT_[0] * c.sqrtkT_[0] / K_BOLTZMANN; - } - return 0; } diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index f02ee42b14..d195a9de39 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -309,7 +309,7 @@ prepare_distribcell() for (auto& filt : model::tally_filters) { auto* distrib_filt = dynamic_cast(filt.get()); if (distrib_filt) { - distribcells.insert(distrib_filt->cell_); + distribcells.insert(distrib_filt->cell()); } } diff --git a/src/material.cpp b/src/material.cpp index f9116a935b..0a6bebf17a 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -47,9 +47,10 @@ std::unordered_map material_map; //============================================================================== Material::Material(pugi::xml_node node) + : index_{model::materials.size()} { if (check_for_node(node, "id")) { - id_ = std::stoi(get_node_value(node, "id")); + this->set_id(std::stoi(get_node_value(node, "id"))); } else { fatal_error("Must specify id of material in materials XML file."); } @@ -331,6 +332,11 @@ Material::Material(pugi::xml_node node) } } +Material::~Material() +{ + model::material_map.erase(id_); +} + void Material::finalize() { // Set fissionable if any nuclide is fissionable @@ -849,11 +855,41 @@ void Material::calculate_photon_xs(Particle& p) const } } -int Material::set_density(double density, std::string units) +void Material::set_id(int32_t id) { + Expects(id >= -1); + + // Clear entry in material map if an ID was already assigned before + if (id_ != -1) { + model::material_map.erase(id_); + id_ = -1; + } + + // Make sure no other material has same ID + if (model::material_map.find(id) != model::material_map.end()) { + throw std::runtime_error{"Two materials have the same ID: " + std::to_string(id)}; + } + + // If no ID specified, auto-assign next ID in sequence + if (id == -1) { + id = 0; + for (const auto& f : model::materials) { + id = std::max(id, f->id_); + } + ++id; + } + + // Update ID and entry in material map + id_ = id; + model::material_map[id] = index_; +} + +void Material::set_density(double density, gsl::cstring_span units) +{ + Expects(density >= 0.0); + if (nuclide_.empty()) { - set_errmsg("No nuclides exist in material yet."); - return OPENMC_E_ALLOCATE; + throw std::runtime_error{"No nuclides exist in material yet."}; } if (units == "atom/b-cm") { @@ -884,10 +920,51 @@ int Material::set_density(double density, std::string units) density_ *= f; atom_density_ *= f; } else { - set_errmsg("Invalid units '" + units + "' specified."); - return OPENMC_E_INVALID_ARGUMENT; + throw std::invalid_argument{"Invalid units '" + std::string(units.data()) + + "' specified."}; } - return 0; +} + +void Material::set_densities(const std::vector& name, + const std::vector& density) +{ + auto n = name.size(); + Expects(n > 0); + Expects(n == density.size()); + + if (n != nuclide_.size()) { + nuclide_.resize(n); + atom_density_ = xt::zeros({n}); + } + + double sum_density = 0.0; + for (gsl::index i = 0; i < n; ++i) { + const auto& nuc {name[i]}; + if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) { + int err = openmc_load_nuclide(nuc.c_str()); + if (err < 0) throw std::runtime_error{openmc_err_msg}; + } + + nuclide_[i] = data::nuclide_map.at(nuc); + Expects(density[i] > 0.0); + atom_density_(i) = density[i]; + sum_density += density[i]; + } + + // Set total density to the sum of the vector + this->set_density(sum_density, "atom/b-cm"); + + // Assign S(a,b) tables + this->init_thermal(); +} + +double Material::volume() const +{ + if (volume_ < 0.0) { + throw std::runtime_error{"Volume for material with ID=" + + std::to_string(id_) + " not set."}; + } + return volume_; } void Material::to_hdf5(hid_t group) const @@ -948,6 +1025,42 @@ void Material::to_hdf5(hid_t group) const close_group(material_group); } +void Material::add_nuclide(const std::string& name, double density) +{ + // Check if nuclide is already in material + for (int i = 0; i < nuclide_.size(); ++i) { + int i_nuc = nuclide_[i]; + if (data::nuclides[i_nuc]->name_ == name) { + double awr = data::nuclides[i_nuc]->awr_; + density_ += density - atom_density_(i); + density_gpcc_ += (density - atom_density_(i)) + * awr * MASS_NEUTRON / N_AVOGADRO; + atom_density_(i) = density; + return; + } + } + + // If nuclide wasn't found, extend nuclide/density arrays + int err = openmc_load_nuclide(name.c_str()); + if (err < 0) throw std::runtime_error{openmc_err_msg}; + + // Append new nuclide/density + int i_nuc = data::nuclide_map[name]; + nuclide_.push_back(i_nuc); + + auto n = nuclide_.size(); + + // Create copy of atom_density_ array with one extra entry + xt::xtensor atom_density = xt::zeros({n}); + xt::view(atom_density, xt::range(0, n-1)) = atom_density_; + atom_density(n-1) = density; + atom_density_ = atom_density; + + density_ += density; + density_gpcc_ += density * data::nuclides[i_nuc]->awr_ + * MASS_NEUTRON / N_AVOGADRO; +} + //============================================================================== // Non-method functions //============================================================================== @@ -1101,19 +1214,6 @@ void read_materials_xml() model::materials.push_back(std::make_unique(material_node)); } model::materials.shrink_to_fit(); - - // Populate the material map. - for (int i = 0; i < model::materials.size(); i++) { - int32_t mid = model::materials[i]->id_; - auto search = model::material_map.find(mid); - if (search == model::material_map.end()) { - model::material_map[mid] = i; - } else { - std::stringstream err_msg; - err_msg << "Two or more materials use the same unique ID: " << mid; - fatal_error(err_msg); - } - } } void free_memory_material() @@ -1144,40 +1244,10 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density) { int err = 0; if (index >= 0 && index < model::materials.size()) { - auto& m = model::materials[index]; - - // Check if nuclide is already in material - for (int i = 0; i < m->nuclide_.size(); ++i) { - int i_nuc = m->nuclide_[i]; - if (data::nuclides[i_nuc]->name_ == name) { - double awr = data::nuclides[i_nuc]->awr_; - m->density_ += density - m->atom_density_(i); - m->density_gpcc_ += (density - m->atom_density_(i)) - * awr * MASS_NEUTRON / N_AVOGADRO; - m->atom_density_(i) = density; - return 0; - } - } - - // If nuclide wasn't found, extend nuclide/density arrays - err = openmc_load_nuclide(name); - - if (err == 0) { - // Append new nuclide/density - int i_nuc = data::nuclide_map[name]; - m->nuclide_.push_back(i_nuc); - - auto n = m->nuclide_.size(); - - // Create copy of atom_density_ array with one extra entry - xt::xtensor atom_density = xt::zeros({n}); - xt::view(atom_density, xt::range(0, n-1)) = m->atom_density_; - atom_density(n-1) = density; - m->atom_density_ = atom_density; - - m->density_ += density; - m->density_gpcc_ += density * data::nuclides[i_nuc]->awr_ - * MASS_NEUTRON / N_AVOGADRO; + try { + model::materials[index]->add_nuclide(name, density); + } catch (const std::runtime_error& e) { + return OPENMC_E_DATA; } } else { set_errmsg("Index in materials array is out of bounds."); @@ -1187,14 +1257,14 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density) } extern "C" int -openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n) +openmc_material_get_densities(int32_t index, const int** nuclides, const double** densities, int* n) { if (index >= 0 && index < model::materials.size()) { auto& mat = model::materials[index]; - if (!mat->nuclide_.empty()) { - *nuclides = mat->nuclide_.data(); - *densities = mat->atom_density_.data(); - *n = mat->nuclide_.size(); + if (!mat->nuclides().empty()) { + *nuclides = mat->nuclides().data(); + *densities = mat->densities().data(); + *n = mat->nuclides().size(); return 0; } else { set_errmsg("Material atom density array has not been allocated."); @@ -1211,7 +1281,7 @@ openmc_material_get_density(int32_t index, double* density) { if (index >= 0 && index < model::materials.size()) { auto& mat = model::materials[index]; - *density = mat->density_gpcc_; + *density = mat->density_gpcc(); return 0; } else { set_errmsg("Index in materials array is out of bounds."); @@ -1223,7 +1293,7 @@ extern "C" int openmc_material_get_fissionable(int32_t index, bool* fissionable) { if (index >= 0 && index < model::materials.size()) { - *fissionable = model::materials[index]->fissionable_; + *fissionable = model::materials[index]->fissionable(); return 0; } else { set_errmsg("Index in materials array is out of bounds."); @@ -1235,7 +1305,7 @@ extern "C" int openmc_material_get_id(int32_t index, int32_t* id) { if (index >= 0 && index < model::materials.size()) { - *id = model::materials[index]->id_; + *id = model::materials[index]->id(); return 0; } else { set_errmsg("Index in materials array is out of bounds."); @@ -1247,16 +1317,13 @@ extern "C" int openmc_material_get_volume(int32_t index, double* volume) { if (index >= 0 && index < model::materials.size()) { - auto& m = model::materials[index]; - if (m->volume_ >= 0.0) { - *volume = m->volume_; - return 0; - } else { - std::stringstream msg; - msg << "Volume for material with ID=" << m->id_ << " not set."; - set_errmsg(msg); + try { + *volume = model::materials[index]->volume(); + } catch (const std::exception& e) { + set_errmsg(e.what()); return OPENMC_E_UNASSIGNED; } + return 0; } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; @@ -1267,59 +1334,51 @@ extern "C" int openmc_material_set_density(int32_t index, double density, const char* units) { if (index >= 0 && index < model::materials.size()) { - return model::materials[index]->set_density(density, units); + try { + model::materials[index]->set_density(density, units); + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; + } } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } + return 0; } extern "C" int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density) { if (index >= 0 && index < model::materials.size()) { - auto& mat {model::materials[index]}; - if (n != mat->nuclide_.size()) { - mat->nuclide_.resize(n); - mat->atom_density_ = xt::zeros({n}); + try { + model::materials[index]->set_densities({name, name + n}, {density, density + n}); + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; } - - double sum_density = 0.0; - for (int i = 0; i < n; ++i) { - std::string nuc {name[i]}; - if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) { - int err = openmc_load_nuclide(nuc.c_str()); - if (err < 0) return err; - } - - mat->nuclide_[i] = data::nuclide_map[nuc]; - mat->atom_density_(i) = density[i]; - sum_density += density[i]; - } - - // Set total density to the sum of the vector - int err = mat->set_density(sum_density, "atom/b-cm"); - - // Assign S(a,b) tables - mat->init_thermal(); - return err; } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } + return 0; } extern "C" int openmc_material_set_id(int32_t index, int32_t id) { if (index >= 0 && index < model::materials.size()) { - model::materials[index]->id_ = id; - model::material_map[id] = index; - return 0; + try { + model::materials.at(index)->set_id(id); + } catch (const std::exception& e) { + set_errmsg(e.what()); + return OPENMC_E_UNASSIGNED; + } } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } + return 0; } extern "C" int diff --git a/src/state_point.cpp b/src/state_point.cpp index 7516b0d485..6ea666fde2 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -141,13 +141,13 @@ openmc_statepoint_write(const char* filename, bool* write_source) std::vector filter_ids; filter_ids.reserve(model::tally_filters.size()); for (const auto& filt : model::tally_filters) - filter_ids.push_back(filt->id_); + filter_ids.push_back(filt->id()); write_attribute(filters_group, "ids", filter_ids); // Write info for each filter for (const auto& filt : model::tally_filters) { hid_t filter_group = create_group(filters_group, - "filter " + std::to_string(filt->id_)); + "filter " + std::to_string(filt->id())); filt->to_statepoint(filter_group); close_group(filter_group); } @@ -188,7 +188,7 @@ openmc_statepoint_write(const char* filename, bool* write_source) std::vector filter_ids; filter_ids.reserve(tally.filters().size()); for (auto i_filt : tally.filters()) - filter_ids.push_back(model::tally_filters[i_filt]->id_); + filter_ids.push_back(model::tally_filters[i_filt]->id()); write_dataset(tally_group, "filters", filter_ids); } diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index e9aae53a7a..c67eec830d 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -7,6 +7,7 @@ #include "openmc/capi.h" #include "openmc/constants.h" // for MAX_LINE_LEN; #include "openmc/error.h" +#include "openmc/xml_interface.h" #include "openmc/tallies/filter_azimuthal.h" #include "openmc/tallies/filter_cell.h" #include "openmc/tallies/filter_cellborn.h" @@ -50,8 +51,46 @@ namespace model { // Non-member functions //============================================================================== -Filter* -allocate_filter(const std::string& type) +extern "C" size_t tally_filters_size() +{ + return model::tally_filters.size(); +} + +//============================================================================== +// Filter implementation +//============================================================================== + +Filter::Filter() : index_{model::tally_filters.size()} +{ } + +Filter::~Filter() +{ + model::filter_map.erase(id_); +} + +Filter* Filter::create(pugi::xml_node node) +{ + // Copy filter id + if (!check_for_node(node, "id")) { + fatal_error("Must specify id for filter in tally XML file."); + } + int filter_id = std::stoi(get_node_value(node, "id")); + + // Convert filter type to lower case + std::string s; + if (check_for_node(node, "type")) { + s = get_node_value(node, "type", true); + } + + // Allocate according to the filter type + auto f = Filter::create(s, filter_id); + + // Read filter data from XML + f->from_xml(node); + return f; +} + +Filter* Filter::create(const std::string& type, int32_t id) { if (type == "azimuthal") { model::tally_filters.push_back(std::make_unique()); @@ -100,12 +139,40 @@ allocate_filter(const std::string& type) } else { throw std::runtime_error{"Unknown filter type: " + type}; } + + // Assign ID + model::tally_filters.back()->set_id(id); + return model::tally_filters.back().get(); } -extern "C" size_t tally_filters_size() +void Filter::set_id(int32_t id) { - return model::tally_filters.size(); + Expects(id >= -1); + + // Clear entry in filter map if an ID was already assigned before + if (id_ != -1) { + model::filter_map.erase(id_); + id_ = -1; + } + + // Make sure no other filter has same ID + if (model::filter_map.find(id) != model::filter_map.end()) { + throw std::runtime_error{"Two filters have the same ID: " + std::to_string(id)}; + } + + // If no ID specified, auto-assign next ID in sequence + if (id == -1) { + id = 0; + for (const auto& f : model::tally_filters) { + id = std::max(id, f->id_); + } + ++id; + } + + // Update ID and entry in filter map + id_ = id; + model::filter_map[id] = index_; } //============================================================================== @@ -126,7 +193,7 @@ openmc_filter_get_id(int32_t index, int32_t* id) { if (int err = verify_filter(index)) return err; - *id = model::tally_filters[index]->id_; + *id = model::tally_filters[index]->id(); return 0; } @@ -135,13 +202,7 @@ openmc_filter_set_id(int32_t index, int32_t id) { if (int err = verify_filter(index)) return err; - if (model::filter_map.find(id) != model::filter_map.end()) { - set_errmsg("Two filters have the same ID: " + std::to_string(id)); - return OPENMC_E_INVALID_ID; - } - - model::tally_filters[index]->id_ = id; - model::filter_map[id] = index; + model::tally_filters[index]->set_id(id); return 0; } @@ -172,7 +233,7 @@ openmc_get_filter_next_id(int32_t* id) { int32_t largest_filter_id = 0; for (const auto& t : model::tally_filters) { - largest_filter_id = std::max(largest_filter_id, t->id_); + largest_filter_id = std::max(largest_filter_id, t->id()); } *id = largest_filter_id + 1; } @@ -181,7 +242,7 @@ extern "C" int openmc_new_filter(const char* type, int32_t* index) { *index = model::tally_filters.size(); - allocate_filter(type); + Filter::create(type); return 0; } diff --git a/src/tallies/filter_azimuthal.cpp b/src/tallies/filter_azimuthal.cpp index 1a33a6bfe3..eb02319c20 100644 --- a/src/tallies/filter_azimuthal.cpp +++ b/src/tallies/filter_azimuthal.cpp @@ -15,22 +15,35 @@ AzimuthalFilter::from_xml(pugi::xml_node node) { auto bins = get_node_array(node, "bins"); - if (bins.size() > 1) { - bins_ = bins; - - } else { + if (bins.size() == 1) { // Allow a user to input a lone number which will mean that you subdivide // [-pi,pi) evenly with the input being the number of bins int n_angle = bins[0]; - - if (n_angle <= 1) fatal_error("Number of bins for azimuthal filter must " - "be greater than 1."); + if (n_angle <= 1) throw std::runtime_error{ + "Number of bins for azimuthal filter must be greater than 1."}; double d_angle = 2.0 * PI / n_angle; - bins_.resize(n_angle + 1); - for (int i = 0; i < n_angle; i++) bins_[i] = -PI + i * d_angle; - bins_[n_angle] = PI; + bins.resize(n_angle + 1); + for (int i = 0; i < n_angle; i++) bins[i] = -PI + i * d_angle; + bins[n_angle] = PI; + } + + this->set_bins(bins); +} + +void AzimuthalFilter::set_bins(gsl::span bins) +{ + // Clear existing bins + bins_.clear(); + bins_.reserve(bins.size()); + + // Copy bins, ensuring they are valid + for (gsl::index i = 0; i < bins.size(); ++i) { + if (i > 0 && bins[i] <= bins[i-1]) { + throw std::runtime_error{"Azimuthal bins must be monotonically increasing."}; + } + bins_.push_back(bins[i]); } n_bins_ = bins_.size() - 1; diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp index 8286a69760..7cc0077568 100644 --- a/src/tallies/filter_cell.cpp +++ b/src/tallies/filter_cell.cpp @@ -12,29 +12,39 @@ namespace openmc { void CellFilter::from_xml(pugi::xml_node node) { - cells_ = get_node_array(node, "bins"); - n_bins_ = cells_.size(); + // Get cell IDs and convert into indices into the global cells vector + auto cells = get_node_array(node, "bins"); + for (auto& c : cells) { + auto search = model::cell_map.find(c); + if (search == model::cell_map.end()) { + std::stringstream err_msg; + err_msg << "Could not find cell " << c + << " specified on tally filter."; + throw std::runtime_error{err_msg.str()}; + } + c = search->second; + } + + this->set_cells(cells); } void -CellFilter::initialize() +CellFilter::set_cells(gsl::span cells) { - // Convert cell IDs to indices of the global array. - for (auto& c : cells_) { - auto search = model::cell_map.find(c); - if (search != model::cell_map.end()) { - c = search->second; - } else { - std::stringstream err_msg; - err_msg << "Could not find cell " << c << " specified on tally filter."; - fatal_error(err_msg); - } + // Clear existing cells + cells_.clear(); + cells_.reserve(cells.size()); + map_.clear(); + + // Update cells and mapping + for (auto& index : cells) { + Expects(index >= 0); + Expects(index < model::cells.size()); + cells_.push_back(index); + map_[index] = cells_.size() - 1; } - // Populate the index->bin map. - for (int i = 0; i < cells_.size(); i++) { - map_[cells_[i]] = i; - } + n_bins_ = cells_.size(); } void @@ -70,7 +80,7 @@ CellFilter::text_label(int bin) const //============================================================================== extern "C" int -openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n) +openmc_cell_filter_get_bins(int32_t index, const int32_t** cells, int32_t* n) { if (int err = verify_filter(index)) return err; @@ -81,8 +91,8 @@ openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n) } auto cell_filt = static_cast(filt); - *cells = cell_filt->cells_.data(); - *n = cell_filt->cells_.size(); + *cells = cell_filt->cells().data(); + *n = cell_filt->cells().size(); return 0; } diff --git a/src/tallies/filter_delayedgroup.cpp b/src/tallies/filter_delayedgroup.cpp index 4f448358e3..2b69783081 100644 --- a/src/tallies/filter_delayedgroup.cpp +++ b/src/tallies/filter_delayedgroup.cpp @@ -8,21 +8,32 @@ namespace openmc { void DelayedGroupFilter::from_xml(pugi::xml_node node) { - groups_ = get_node_array(node, "bins"); - n_bins_ = groups_.size(); + auto groups = get_node_array(node, "bins"); + this->set_groups(groups); +} + +void +DelayedGroupFilter::set_groups(gsl::span groups) +{ + // Clear existing groups + groups_.clear(); + groups_.reserve(groups.size()); // Make sure all the group index values are valid. // TODO: do these need to be decremented for zero-based indexing? - for (auto group : groups_) { + for (auto group : groups) { if (group < 1) { - fatal_error("Encountered delayedgroup bin with index " - + std::to_string(group) + " which is less than 1"); + throw std::invalid_argument{"Encountered delayedgroup bin with index " + + std::to_string(group) + " which is less than 1"}; } else if (group > MAX_DELAYED_GROUPS) { - fatal_error("Encountered delayedgroup bin with index " + throw std::invalid_argument{"Encountered delayedgroup bin with index " + std::to_string(group) + " which is greater than MAX_DELATED_GROUPS (" - + std::to_string(MAX_DELAYED_GROUPS) + ")"); + + std::to_string(MAX_DELAYED_GROUPS) + ")"}; } + groups_.push_back(group); } + + n_bins_ = groups_.size(); } void diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp index 94f23621b1..f9c7ab4d66 100644 --- a/src/tallies/filter_distribcell.cpp +++ b/src/tallies/filter_distribcell.cpp @@ -15,23 +15,26 @@ DistribcellFilter::from_xml(pugi::xml_node node) if (cells.size() != 1) { fatal_error("Only one cell can be specified per distribcell filter."); } - cell_ = cells[0]; -} -void -DistribcellFilter::initialize() -{ - // Convert the cell ID to an index of the global array. - auto search = model::cell_map.find(cell_); - if (search != model::cell_map.end()) { - cell_ = search->second; - n_bins_ = model::cells[cell_]->n_instances_; - } else { + // Find index in global cells vector corresponding to cell ID + auto search = model::cell_map.find(cells[0]); + if (search == model::cell_map.end()) { std::stringstream err_msg; err_msg << "Could not find cell " << cell_ << " specified on tally filter."; - fatal_error(err_msg); + throw std::runtime_error{err_msg.str()}; } + + this->set_cell(search->second); +} + +void +DistribcellFilter::set_cell(int32_t cell) +{ + Expects(cell >= 0); + Expects(cell < model::cells.size()); + cell_ = cell; + n_bins_ = model::cells[cell]->n_instances_; } void diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index e11ab30803..dde9b692ac 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -16,7 +16,25 @@ namespace openmc { void EnergyFilter::from_xml(pugi::xml_node node) { - bins_ = get_node_array(node, "bins"); + auto bins = get_node_array(node, "bins"); + this->set_bins(bins); +} + +void +EnergyFilter::set_bins(gsl::span bins) +{ + // Clear existing bins + bins_.clear(); + bins_.reserve(bins.size()); + + // Copy bins, ensuring they are valid + for (gsl::index i = 0; i < bins.size(); ++i) { + if (i > 0 && bins[i] <= bins[i-1]) { + throw std::runtime_error{"Energy bins must be monotonically increasing."}; + } + bins_.push_back(bins[i]); + } + n_bins_ = bins_.size() - 1; // In MG mode, check if the filter bins match the transport bins. @@ -27,7 +45,7 @@ EnergyFilter::from_xml(pugi::xml_node node) if (!settings::run_CE) { if (n_bins_ == data::num_energy_groups) { matches_transport_groups_ = true; - for (auto i = 0; i < n_bins_ + 1; i++) { + for (gsl::index i = 0; i < n_bins_ + 1; ++i) { if (data::rev_energy_bins[i] != bins_[i]) { matches_transport_groups_ = false; break; @@ -111,7 +129,7 @@ EnergyoutFilter::text_label(int bin) const //============================================================================== extern"C" int -openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n) +openmc_energy_filter_get_bins(int32_t index, const double** energies, size_t* n) { // Make sure this is a valid index to an allocated filter. if (int err = verify_filter(index)) return err; @@ -127,13 +145,13 @@ openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n) } // Output the bins. - *energies = filt->bins_.data(); - *n = filt->bins_.size(); + *energies = filt->bins().data(); + *n = filt->bins().size(); return 0; } extern "C" int -openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies) +openmc_energy_filter_set_bins(int32_t index, size_t n, const double* energies) { // Make sure this is a valid index to an allocated filter. if (int err = verify_filter(index)) return err; @@ -149,10 +167,7 @@ openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies) } // Update the filter. - filt->bins_.clear(); - filt->bins_.resize(n); - for (int i = 0; i < n; i++) filt->bins_[i] = energies[i]; - filt->n_bins_ = n - 1; + filt->set_bins({energies, n}); return 0; } diff --git a/src/tallies/filter_legendre.cpp b/src/tallies/filter_legendre.cpp index 6fa041c76b..3961e06d40 100644 --- a/src/tallies/filter_legendre.cpp +++ b/src/tallies/filter_legendre.cpp @@ -10,7 +10,16 @@ namespace openmc { void LegendreFilter::from_xml(pugi::xml_node node) { - order_ = std::stoi(get_node_value(node, "order")); + this->set_order(std::stoi(get_node_value(node, "order"))); +} + +void +LegendreFilter::set_order(int order) +{ + if (order < 0) { + throw std::invalid_argument{"Legendre order must be non-negative."}; + } + order_ = order; n_bins_ = order_ + 1; } @@ -60,7 +69,7 @@ openmc_legendre_filter_get_order(int32_t index, int* order) } // Output the order. - *order = filt->order_; + *order = filt->order(); return 0; } @@ -81,8 +90,7 @@ openmc_legendre_filter_set_order(int32_t index, int order) } // Update the filter. - filt->order_ = order; - filt->n_bins_ = order + 1; + filt->set_order(order); return 0; } diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp index fd9e637bde..bfc216244a 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -3,7 +3,6 @@ #include #include "openmc/capi.h" -#include "openmc/error.h" #include "openmc/material.h" #include "openmc/xml_interface.h" @@ -12,30 +11,39 @@ namespace openmc { void MaterialFilter::from_xml(pugi::xml_node node) { - materials_ = get_node_array(node, "bins"); - n_bins_ = materials_.size(); -} - -void -MaterialFilter::initialize() -{ - // Convert material IDs to indices of the global array. - for (auto& m : materials_) { + // Get material IDs and convert to indices in the global materials vector + auto mats = get_node_array(node, "bins"); + for (auto& m : mats) { auto search = model::material_map.find(m); - if (search != model::material_map.end()) { - m = search->second; - } else { + if (search == model::material_map.end()) { std::stringstream err_msg; err_msg << "Could not find material " << m << " specified on tally filter."; - fatal_error(err_msg); + throw std::runtime_error{err_msg.str()}; } + m = search->second; } - // Populate the index->bin map. - for (int i = 0; i < materials_.size(); i++) { - map_[materials_[i]] = i; + this->set_materials(mats); +} + +void +MaterialFilter::set_materials(gsl::span materials) +{ + // Clear existing materials + materials_.clear(); + materials_.reserve(materials.size()); + map_.clear(); + + // Update materials and mapping + for (auto& index : materials) { + Expects(index >= 0); + Expects(index < model::materials.size()); + materials_.push_back(index); + map_[index] = materials_.size() - 1; } + + n_bins_ = materials_.size(); } void @@ -69,7 +77,7 @@ MaterialFilter::text_label(int bin) const //============================================================================== extern "C" int -openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n) +openmc_material_filter_get_bins(int32_t index, const int32_t** bins, size_t* n) { // Make sure this is a valid index to an allocated filter. if (int err = verify_filter(index)) return err; @@ -85,13 +93,13 @@ openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n) } // Output the bins. - *bins = filt->materials_.data(); - *n = filt->materials_.size(); + *bins = filt->materials().data(); + *n = filt->materials().size(); return 0; } extern "C" int -openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins) +openmc_material_filter_set_bins(int32_t index, size_t n, const int32_t* bins) { // Make sure this is a valid index to an allocated filter. if (int err = verify_filter(index)) return err; @@ -107,12 +115,7 @@ openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins) } // Update the filter. - 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_[filt->materials_[i]] = i; + filt->set_materials({bins, n}); return 0; } diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp index aed8371b80..9aaebea8cf 100644 --- a/src/tallies/filter_mu.cpp +++ b/src/tallies/filter_mu.cpp @@ -13,22 +13,36 @@ MuFilter::from_xml(pugi::xml_node node) { auto bins = get_node_array(node, "bins"); - if (bins.size() > 1) { - bins_ = bins; - - } else { + if (bins.size() == 1) { // Allow a user to input a lone number which will mean that you subdivide // [-1,1) evenly with the input being the number of bins int n_angle = bins[0]; - - if (n_angle <= 1) fatal_error("Number of bins for mu filter must " - "be greater than 1."); + if (n_angle <= 1) throw std::runtime_error{ + "Number of bins for mu filter must be greater than 1."}; double d_angle = 2.0 / n_angle; - bins_.resize(n_angle + 1); - for (int i = 0; i < n_angle; i++) bins_[i] = -1 + i * d_angle; - bins_[n_angle] = 1; + bins.resize(n_angle + 1); + for (int i = 0; i < n_angle; i++) bins[i] = -1 + i * d_angle; + bins[n_angle] = 1; + } + + this->set_bins(bins); +} + +void +MuFilter::set_bins(gsl::span bins) +{ + // Clear existing bins + bins_.clear(); + bins_.reserve(bins.size()); + + // Copy bins, ensuring they are valid + for (gsl::index i = 0; i < bins.size(); ++i) { + if (i > 0 && bins[i] <= bins[i-1]) { + throw std::runtime_error{"Mu bins must be monotonically increasing."}; + } + bins_.push_back(bins[i]); } n_bins_ = bins_.size() - 1; diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index 23a305b8dd..de142489ff 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -8,8 +8,25 @@ void ParticleFilter::from_xml(pugi::xml_node node) { auto particles = get_node_array(node, "bins"); + + // Convert to vector of Particle::Type + std::vector types; for (auto& p : particles) { - particles_.push_back(static_cast(p - 1)); + types.push_back(static_cast(p - 1)); + } + this->set_particles(types); +} + +void +ParticleFilter::set_particles(gsl::span particles) +{ + // Clear existing particles + particles_.clear(); + particles_.reserve(particles.size()); + + // Set particles and number of bins + for (auto p : particles) { + particles_.push_back(p); } n_bins_ = particles_.size(); } diff --git a/src/tallies/filter_polar.cpp b/src/tallies/filter_polar.cpp index 9730bbb584..edbbc89229 100644 --- a/src/tallies/filter_polar.cpp +++ b/src/tallies/filter_polar.cpp @@ -14,22 +14,36 @@ PolarFilter::from_xml(pugi::xml_node node) { auto bins = get_node_array(node, "bins"); - if (bins.size() > 1) { - bins_ = bins; - - } else { + if (bins.size() == 1) { // Allow a user to input a lone number which will mean that you subdivide // [0,pi] evenly with the input being the number of bins int n_angle = bins[0]; - - if (n_angle <= 1) fatal_error("Number of bins for polar filter must " - "be greater than 1."); + if (n_angle <= 1) throw std::runtime_error{ + "Number of bins for polar filter must be greater than 1."}; double d_angle = PI / n_angle; - bins_.resize(n_angle + 1); - for (int i = 0; i < n_angle; i++) bins_[i] = i * d_angle; - bins_[n_angle] = PI; + bins.resize(n_angle + 1); + for (int i = 0; i < n_angle; i++) bins[i] = i * d_angle; + bins[n_angle] = PI; + } + + this->set_bins(bins); +} + +void +PolarFilter::set_bins(gsl::span bins) +{ + // Clear existing bins + bins_.clear(); + bins_.reserve(bins.size()); + + // Copy bins, ensuring they are valid + for (gsl::index i = 0; i < bins.size(); ++i) { + if (i > 0 && bins[i] <= bins[i-1]) { + throw std::runtime_error{"Polar bins must be monotonically increasing."}; + } + bins_.push_back(bins[i]); } n_bins_ = bins_.size() - 1; diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp index 0b3a371187..80d6878b70 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -12,21 +12,34 @@ namespace openmc { void SphericalHarmonicsFilter::from_xml(pugi::xml_node node) { - order_ = std::stoi(get_node_value(node, "order")); - n_bins_ = (order_ + 1) * (order_ + 1); - + this->set_order(std::stoi(get_node_value(node, "order"))); if (check_for_node(node, "cosine")) { - auto cos = get_node_value(node, "cosine", true); - if (cos == "scatter") { - cosine_ = SphericalHarmonicsCosine::scatter; - } else if (cos == "particle") { - cosine_ = SphericalHarmonicsCosine::particle; - } else { - std::stringstream err_msg; - err_msg << "Unrecognized cosine type, \"" << cos - << "\" in spherical harmonics filter"; - fatal_error(err_msg); - } + this->set_cosine(get_node_value(node, "cosine", true)); + } +} + +void +SphericalHarmonicsFilter::set_order(int order) +{ + if (order < 0) { + throw std::invalid_argument{"Spherical harmonics order must be non-negative."}; + } + order_ = order; + n_bins_ = (order_ + 1) * (order_ + 1); +} + +void +SphericalHarmonicsFilter::set_cosine(gsl::cstring_span cosine) +{ + if (cosine == "scatter") { + cosine_ = SphericalHarmonicsCosine::scatter; + } else if (cosine == "particle") { + cosine_ = SphericalHarmonicsCosine::particle; + } else { + std::stringstream err_msg; + err_msg << "Unrecognized cosine type, \"" << cosine + << "\" in spherical harmonics filter"; + throw std::invalid_argument{err_msg.str()}; } } @@ -121,7 +134,7 @@ openmc_sphharm_filter_get_order(int32_t index, int* order) if (err) return err; // Output the order. - *order = filt->order_; + *order = filt->order(); return 0; } @@ -135,7 +148,7 @@ openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]) if (err) return err; // Output the cosine. - if (filt->cosine_ == SphericalHarmonicsCosine::scatter) { + if (filt->cosine() == SphericalHarmonicsCosine::scatter) { strcpy(cosine, "scatter"); } else { strcpy(cosine, "particle"); @@ -153,8 +166,7 @@ openmc_sphharm_filter_set_order(int32_t index, int order) if (err) return err; // Update the filter. - filt->order_ = order; - filt->n_bins_ = (order + 1) * (order + 1); + filt->set_order(order); return 0; } @@ -168,12 +180,10 @@ openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]) if (err) return err; // Update the filter. - if (strcmp(cosine, "scatter") == 0) { - filt->cosine_ = SphericalHarmonicsCosine::scatter; - } else if (strcmp(cosine, "particle") == 0) { - filt->cosine_ = SphericalHarmonicsCosine::particle; - } else { - set_errmsg("Invalid spherical harmonics cosine."); + try { + filt->set_cosine(cosine); + } catch (const std::invalid_argument& e) { + set_errmsg(e.what()); return OPENMC_E_INVALID_ARGUMENT; } return 0; diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp index 6562ff01c7..40d82f4c1f 100644 --- a/src/tallies/filter_sptl_legendre.cpp +++ b/src/tallies/filter_sptl_legendre.cpp @@ -12,25 +12,54 @@ namespace openmc { void SpatialLegendreFilter::from_xml(pugi::xml_node node) { - order_ = std::stoi(get_node_value(node, "order")); + this->set_order(std::stoi(get_node_value(node, "order"))); auto axis = get_node_value(node, "axis"); - if (axis == "x") { - axis_ = LegendreAxis::x; - } else if (axis == "y") { - axis_ = LegendreAxis::y; - } else if (axis == "z") { - axis_ = LegendreAxis::z; - } else { - fatal_error("Unrecognized axis on SpatialLegendreFilter"); + switch (axis[0]) { + case 'x': + this->set_axis(LegendreAxis::x); + break; + case 'y': + this->set_axis(LegendreAxis::y); + break; + case 'z': + this->set_axis(LegendreAxis::z); + break; + default: + throw std::runtime_error{"Axis for SpatialLegendreFilter must be 'x', 'y', or 'z'"}; } - min_ = std::stod(get_node_value(node, "min")); - max_ = std::stod(get_node_value(node, "max")); + double min = std::stod(get_node_value(node, "min")); + double max = std::stod(get_node_value(node, "max")); + this->set_minmax(min, max); +} +void +SpatialLegendreFilter::set_order(int order) +{ + if (order < 0) { + throw std::invalid_argument{"Legendre order must be non-negative."}; + } + order_ = order; n_bins_ = order_ + 1; } +void +SpatialLegendreFilter::set_axis(LegendreAxis axis) +{ + axis_ = axis; +} + +void +SpatialLegendreFilter::set_minmax(double min, double max) +{ + if (max < min) { + throw std::invalid_argument{"Maximum value must be greater than minimum value"}; + } + min_ = min; + max_ = max; +} + void SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const @@ -126,7 +155,7 @@ openmc_spatial_legendre_filter_get_order(int32_t index, int* order) if (err) return err; // Output the order. - *order = filt->order_; + *order = filt->order(); return 0; } @@ -141,9 +170,9 @@ openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, if (err) return err; // Output the params. - *axis = static_cast(filt->axis_); - *min = filt->min_; - *max = filt->max_; + *axis = static_cast(filt->axis()); + *min = filt->min(); + *max = filt->max(); return 0; } @@ -157,8 +186,7 @@ openmc_spatial_legendre_filter_set_order(int32_t index, int order) if (err) return err; // Update the filter. - filt->order_ = order; - filt->n_bins_ = order + 1; + filt->set_order(order); return 0; } @@ -173,9 +201,8 @@ openmc_spatial_legendre_filter_set_params(int32_t index, const int* axis, if (err) return err; // Update the filter. - if (axis) filt->axis_ = static_cast(*axis); - if (min) filt->min_ = *min; - if (max) filt->max_ = *max; + if (axis) filt->set_axis(static_cast(*axis)); + if (min && max) filt->set_minmax(*min, *max); return 0; } diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp index faa8aafbf7..72b356f216 100644 --- a/src/tallies/filter_surface.cpp +++ b/src/tallies/filter_surface.cpp @@ -11,30 +11,41 @@ namespace openmc { void SurfaceFilter::from_xml(pugi::xml_node node) { - surfaces_ = get_node_array(node, "bins"); - n_bins_ = surfaces_.size(); -} + auto surfaces = get_node_array(node, "bins"); -void -SurfaceFilter::initialize() -{ - // Convert surface IDs to indices of the global array. - for (auto& s : surfaces_) { + // Convert surface IDs to indices of the global surfaces vector. + for (auto& s : surfaces) { auto search = model::surface_map.find(s); - if (search != model::surface_map.end()) { - s = search->second; - } else { + if (search == model::surface_map.end()) { std::stringstream err_msg; err_msg << "Could not find surface " << s << " specified on tally filter."; - fatal_error(err_msg); + throw std::runtime_error{err_msg.str()}; } + + s = search->second; } - // Populate the index->bin map. - for (int i = 0; i < surfaces_.size(); i++) { - map_[surfaces_[i]] = i; + this->set_surfaces(surfaces); +} + +void +SurfaceFilter::set_surfaces(gsl::span surfaces) +{ + // Clear existing surfaces + surfaces_.clear(); + surfaces_.reserve(surfaces.size()); + map_.clear(); + + // Update surfaces and mapping + for (auto& index : surfaces) { + Expects(index >= 0); + Expects(index < model::surfaces.size()); + surfaces_.push_back(index); + map_[index] = surfaces_.size() - 1; } + + n_bins_ = surfaces_.size(); } void diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp index 50c058c05b..dffdee621f 100644 --- a/src/tallies/filter_universe.cpp +++ b/src/tallies/filter_universe.cpp @@ -11,30 +11,39 @@ namespace openmc { void UniverseFilter::from_xml(pugi::xml_node node) { - universes_ = get_node_array(node, "bins"); - n_bins_ = universes_.size(); -} - -void -UniverseFilter::initialize() -{ - // Convert universe IDs to indices of the global array. - for (auto& u : universes_) { + // Get material IDs and convert to indices in the global materials vector + auto universes = get_node_array(node, "bins"); + for (auto& u : universes) { auto search = model::universe_map.find(u); - if (search != model::universe_map.end()) { - u = search->second; - } else { + if (search == model::universe_map.end()) { std::stringstream err_msg; err_msg << "Could not find universe " << u << " specified on tally filter."; - fatal_error(err_msg); + throw std::runtime_error{err_msg.str()}; } + u = search->second; } - // Populate the index->bin map. - for (int i = 0; i < universes_.size(); i++) { - map_[universes_[i]] = i; + this->set_universes(universes); +} + +void +UniverseFilter::set_universes(gsl::span universes) +{ + // Clear existing universes + universes_.clear(); + universes_.reserve(universes.size()); + map_.clear(); + + // Update universes and mapping + for (auto& index : universes) { + Expects(index >= 0); + Expects(index < model::universes.size()); + universes_.push_back(index); + map_[index] = universes_.size() - 1; } + + n_bins_ = universes_.size(); } void diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp index 125bfada89..b732ef7931 100644 --- a/src/tallies/filter_zernike.cpp +++ b/src/tallies/filter_zernike.cpp @@ -74,6 +74,9 @@ ZernikeFilter::text_label(int bin) const void ZernikeFilter::set_order(int order) { + if (order < 0) { + throw std::invalid_argument{"Zernike order must be non-negative."}; + } order_ = order; n_bins_ = ((order+1) * (order+2)) / 2; } @@ -111,7 +114,7 @@ ZernikeRadialFilter::text_label(int bin) const void ZernikeRadialFilter::set_order(int order) { - order_ = order; + ZernikeFilter::set_order(order); n_bins_ = order / 2 + 1; } @@ -165,9 +168,9 @@ openmc_zernike_filter_get_params(int32_t index, double* x, double* y, if (err) return err; // Output the params. - *x = filt->x_; - *y = filt->y_; - *r = filt->r_; + *x = filt->x(); + *y = filt->y(); + *r = filt->r(); return 0; } @@ -196,9 +199,9 @@ openmc_zernike_filter_set_params(int32_t index, const double* x, if (err) return err; // Update the filter. - if (x) filt->x_ = *x; - if (y) filt->y_ = *y; - if (r) filt->r_ = *r; + if (x) filt->set_x(*x); + if (y) filt->set_y(*y); + if (r) filt->set_r(*r); return 0; } diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 2cb0945ab3..27d6cc58c0 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -240,37 +240,293 @@ score_str_to_int(std::string score_str) // Tally object implementation //============================================================================== -Tally::Tally() +Tally::Tally(int32_t id) + : index_{model::tallies.size()} { - this->set_filters(nullptr, 0); + this->set_id(id); + this->set_filters({}); } -void -Tally::init_from_xml(pugi::xml_node node) +Tally::Tally(pugi::xml_node node) + : index_{model::tallies.size()} { + // Copy and set tally id + if (!check_for_node(node, "id")) { + throw std::runtime_error{"Must specify id for tally in tally XML file."}; + } + int32_t id = std::stoi(get_node_value(node, "id")); + this->set_id(id); + if (check_for_node(node, "name")) name_ = get_node_value(node, "name"); + + // ======================================================================= + // READ DATA FOR FILTERS + + // Check if user is using old XML format and throw an error if so + if (check_for_node(node, "filter")) { + throw std::runtime_error{"Tally filters must be specified independently of " + "tallies in a element. The element itself should " + "have a list of filters that apply, e.g., 1 2 " + "where 1 and 2 are the IDs of filters specified outside of " + "."}; + } + + // Determine number of filters + std::vector filter_ids; + if (check_for_node(node, "filters")) { + filter_ids = get_node_array(node, "filters"); + } + + // Allocate and store filter user ids + std::vector filters; + for (int filter_id : filter_ids) { + // Determine if filter ID is valid + auto it = model::filter_map.find(filter_id); + if (it == model::filter_map.end()) { + throw std::runtime_error{"Could not find filter " + std::to_string(filter_id) + + " specified on tally " + std::to_string(id_)}; + } + + // Store the index of the filter + filters.push_back(model::tally_filters[it->second].get()); + } + + // Set the filters + this->set_filters(filters); + + // Check for the presence of certain filter types + bool has_energyout = energyout_filter_ >= 0; + int particle_filter_index = C_NONE; + for (gsl::index j = 0; j < filters_.size(); ++j) { + int i_filter = filters_[j]; + const auto& f = model::tally_filters[i_filter].get(); + + auto pf = dynamic_cast(f); + if (pf) particle_filter_index = i_filter; + + // Change the tally estimator if a filter demands it + std::string filt_type = f->type(); + if (filt_type == "energyout" || filt_type == "legendre") { + estimator_ = ESTIMATOR_ANALOG; + } else if (filt_type == "sphericalharmonics") { + auto sf = dynamic_cast(f); + if (sf->cosine() == SphericalHarmonicsCosine::scatter) { + estimator_ = ESTIMATOR_ANALOG; + } + } else if (filt_type == "spatiallegendre" || filt_type == "zernike" + || filt_type == "zernikeradial") { + estimator_ = ESTIMATOR_COLLISION; + } + } + + // ======================================================================= + // READ DATA FOR NUCLIDES + + this->set_nuclides(node); + + // ======================================================================= + // READ DATA FOR SCORES + + this->set_scores(node); + + if (!check_for_node(node, "scores")) { + fatal_error("No scores specified on tally " + std::to_string(id_) + + "."); + } + + // Check if tally is compatible with particle type + if (settings::photon_transport) { + if (particle_filter_index == C_NONE) { + for (int score : scores_) { + switch (score) { + case SCORE_INVERSE_VELOCITY: + fatal_error("Particle filter must be used with photon " + "transport on and inverse velocity score"); + break; + case SCORE_FLUX: + case SCORE_TOTAL: + case SCORE_SCATTER: + case SCORE_NU_SCATTER: + case SCORE_ABSORPTION: + case SCORE_FISSION: + case SCORE_NU_FISSION: + case SCORE_CURRENT: + case SCORE_EVENTS: + case SCORE_DELAYED_NU_FISSION: + case SCORE_PROMPT_NU_FISSION: + case SCORE_DECAY_RATE: + warning("Particle filter is not used with photon transport" + " on and " + reaction_name(score) + " score."); + break; + } + } + } else { + const auto& f = model::tally_filters[particle_filter_index].get(); + auto pf = dynamic_cast(f); + for (auto p : pf->particles()) { + if (p == Particle::Type::electron || + p == Particle::Type::positron) { + estimator_ = ESTIMATOR_ANALOG; + } + } + } + } else { + if (particle_filter_index >= 0) { + const auto& f = model::tally_filters[particle_filter_index].get(); + auto pf = dynamic_cast(f); + for (auto p : pf->particles()) { + if (p != Particle::Type::neutron) { + warning("Particle filter other than NEUTRON used with photon " + "transport turned off. All tallies for particle type " + + std::to_string(static_cast(p)) + " will have no scores"); + } + } + } + } + + // Check for a tally derivative. + if (check_for_node(node, "derivative")) { + int deriv_id = std::stoi(get_node_value(node, "derivative")); + + // Find the derivative with the given id, and store it's index. + auto it = model::tally_deriv_map.find(deriv_id); + if (it == model::tally_deriv_map.end()) { + fatal_error("Could not find derivative " + std::to_string(deriv_id) + + " specified on tally " + std::to_string(id_)); + } + + deriv_ = it->second; + + // Only analog or collision estimators are supported for differential + // tallies. + if (estimator_ == ESTIMATOR_TRACKLENGTH) { + estimator_ = ESTIMATOR_COLLISION; + } + + const auto& deriv = model::tally_derivs[deriv_]; + if (deriv.variable == DIFF_NUCLIDE_DENSITY + || deriv.variable == DIFF_TEMPERATURE) { + for (int i_nuc : nuclides_) { + if (has_energyout && i_nuc == -1) { + fatal_error("Error on tally " + std::to_string(id_) + + ": Cannot use a 'nuclide_density' or 'temperature' " + "derivative on a tally with an outgoing energy filter and " + "'total' nuclide rate. Instead, tally each nuclide in the " + "material individually."); + // Note that diff tallies with these characteristics would work + // correctly if no tally events occur in the perturbed material + // (e.g. pertrubing moderator but only tallying fuel), but this + // case would be hard to check for by only reading inputs. + } + } + } + } + + // If settings.xml trigger is turned on, create tally triggers + if (settings::trigger_on) { + this->init_triggers(node); + } + + // ======================================================================= + // SET TALLY ESTIMATOR + + // Check if user specified estimator + if (check_for_node(node, "estimator")) { + std::string est = get_node_value(node, "estimator"); + if (est == "analog") { + estimator_ = ESTIMATOR_ANALOG; + } else if (est == "tracklength" || est == "track-length" + || est == "pathlength" || est == "path-length") { + // If the estimator was set to an analog estimator, this means the + // tally needs post-collision information + if (estimator_ == ESTIMATOR_ANALOG) { + throw std::runtime_error{"Cannot use track-length estimator for tally " + + std::to_string(id_)}; + } + + // Set estimator to track-length estimator + estimator_ = ESTIMATOR_TRACKLENGTH; + + } else if (est == "collision") { + // If the estimator was set to an analog estimator, this means the + // tally needs post-collision information + if (estimator_ == ESTIMATOR_ANALOG) { + throw std::runtime_error{"Cannot use collision estimator for tally " + + std::to_string(id_)}; + } + + // Set estimator to collision estimator + estimator_ = ESTIMATOR_COLLISION; + + } else { + throw std::runtime_error{"Invalid estimator '" + est + "' on tally " + + std::to_string(id_)}; + } + } +} + +Tally::~Tally() +{ + model::tally_map.erase(id_); +} + +Tally* +Tally::create(int32_t id) +{ + model::tallies.push_back(std::make_unique(id)); + return model::tallies.back().get(); } void -Tally::set_filters(const int32_t filter_indices[], int n) +Tally::set_id(int32_t id) +{ + Expects(id >= -1); + + // Clear entry in tally map if an ID was already assigned before + if (id_ != -1) { + model::tally_map.erase(id_); + id_ = -1; + } + + // Make sure no other tally has the same ID + if (model::tally_map.find(id) != model::tally_map.end()) { + throw std::runtime_error{"Two tallies have the same ID: " + std::to_string(id)}; + } + + // If no ID specified, auto-assign next ID in sequence + if (id == -1) { + id = 0; + for (const auto& t : model::tallies) { + id = std::max(id, t->id_); + } + ++id; + } + + // Update ID and entry in tally map + id_ = id; + model::tally_map[id] = index_; +} + +void +Tally::set_filters(gsl::span filters) { // Clear old data. filters_.clear(); strides_.clear(); // Copy in the given filter indices. - filters_.assign(filter_indices, filter_indices + n); + auto n = filters.size(); + filters_.reserve(n); for (int i = 0; i < n; ++i) { - auto i_filt = filters_[i]; - if (i_filt < 0 || i_filt >= model::tally_filters.size()) - throw std::out_of_range("Index in tally filter array out of bounds."); + // Add index to vector of filters + auto& f {filters[i]}; + filters_.push_back(model::filter_map.at(f->id())); // Keep track of indices for special filters. - const auto* filt = model::tally_filters[i_filt].get(); - if (dynamic_cast(filt)) { + if (dynamic_cast(f)) { energyout_filter_ = i; - } else if (dynamic_cast(filt)) { + } else if (dynamic_cast(f)) { delayedgroup_filter_ = i; } } @@ -282,7 +538,7 @@ Tally::set_filters(const int32_t filter_indices[], int n) int stride = 1; for (int i = n-1; i >= 0; --i) { strides_[i] = stride; - stride *= model::tally_filters[filters_[i]]->n_bins_; + stride *= model::tally_filters[filters_[i]]->n_bins(); } n_filter_bins_ = stride; } @@ -298,7 +554,7 @@ Tally::set_scores(pugi::xml_node node) } void -Tally::set_scores(std::vector scores) +Tally::set_scores(const std::vector& scores) { // Reset state and prepare for the new scores. scores_.clear(); @@ -450,17 +706,25 @@ Tally::set_nuclides(pugi::xml_node node) // The user provided specifics nuclides. Parse it as an array with either // "total" or a nuclide name like "U-235" in each position. auto words = get_node_array(node, "nuclides"); - for (auto word : words) { - if (word == "total") { - nuclides_.push_back(-1); - } else { - auto search = data::nuclide_map.find(word); - if (search == data::nuclide_map.end()) - fatal_error("Could not find the nuclide " + word - + " specified in tally " + std::to_string(id_) - + " in any material"); - nuclides_.push_back(search->second); - } + this->set_nuclides(words); + } +} + +void +Tally::set_nuclides(const std::vector& nuclides) +{ + nuclides_.clear(); + + for (const auto& nuc : nuclides) { + if (nuc == "total") { + nuclides_.push_back(-1); + } else { + auto search = data::nuclide_map.find(nuc); + if (search == data::nuclide_map.end()) + fatal_error("Could not find the nuclide " + nuc + + " specified in tally " + std::to_string(id_) + + " in any material"); + nuclides_.push_back(search->second); } } } @@ -614,36 +878,7 @@ void read_tallies_xml() // Check for user filters and allocate for (auto node_filt : root.children("filter")) { - // Copy filter id - if (!check_for_node(node_filt, "id")) { - fatal_error("Must specify id for filter in tally XML file."); - } - int filter_id = std::stoi(get_node_value(node_filt, "id")); - - // Check to make sure 'id' hasn't been used - if (model::filter_map.find(filter_id) != model::filter_map.end()) { - fatal_error("Two or more filters use the same unique ID: " - + std::to_string(filter_id)); - } - - // Convert filter type to lower case - std::string s; - if (check_for_node(node_filt, "type")) { - s = get_node_value(node_filt, "type", true); - } - - // Allocate according to the filter type - Filter* f = allocate_filter(s); - - // Read filter data from XML - f->from_xml(node_filt); - - // Set filter id - f->id_ = filter_id; - model::filter_map[filter_id] = model::tally_filters.size() - 1; - - // Initialize filter - f->initialize(); + auto f = Filter::create(node_filt); } // ========================================================================== @@ -657,229 +892,7 @@ void read_tallies_xml() } for (auto node_tal : root.children("tally")) { - model::tallies.push_back(std::make_unique()); - - auto& t {model::tallies.back()}; - t->init_from_xml(node_tal); - - // Copy and set tally id - if (!check_for_node(node_tal, "id")) { - fatal_error("Must specify id for tally in tally XML file."); - } - t->id_ = std::stoi(get_node_value(node_tal, "id")); - model::tally_map[t->id_] = model::tallies.size() - 1; - - // Copy tally name - if (check_for_node(node_tal, "name")) { - t->name_ = get_node_value(node_tal, "name"); - } - - // ======================================================================= - // READ DATA FOR FILTERS - - // Check if user is using old XML format and throw an error if so - if (check_for_node(node_tal, "filter")) { - fatal_error("Tally filters must be specified independently of " - "tallies in a element. The element itself should " - "have a list of filters that apply, e.g., 1 2 " - "where 1 and 2 are the IDs of filters specified outside of " - "."); - } - - // Determine number of filters - std::vector filters; - if (check_for_node(node_tal, "filters")) { - filters = get_node_array(node_tal, "filters"); - } - - // Allocate and store filter user ids - if (!filters.empty()) { - std::vector filter_indices; - for (int filter_id : filters) { - // Determine if filter ID is valid - auto it = model::filter_map.find(filter_id); - if (it == model::filter_map.end()) { - fatal_error("Could not find filter " + std::to_string(filter_id) - + " specified on tally " + std::to_string(t->id_)); - } - - // Store the index of the filter - filter_indices.push_back(it->second); - } - - // Set the filters - t->set_filters(filter_indices.data(), filter_indices.size()); - } - - // Check for the presence of certain filter types - bool has_energyout = t->energyout_filter_ >= 0; - int particle_filter_index = C_NONE; - for (int j = 0; j < t->filters().size(); ++j) { - int i_filter = t->filters(j); - const auto& f = model::tally_filters[i_filter].get(); - - auto pf = dynamic_cast(f); - if (pf) particle_filter_index = i_filter; - - // Change the tally estimator if a filter demands it - std::string filt_type = f->type(); - if (filt_type == "energyout" || filt_type == "legendre") { - t->estimator_ = ESTIMATOR_ANALOG; - } else if (filt_type == "sphericalharmonics") { - auto sf = dynamic_cast(f); - if (sf->cosine_ == SphericalHarmonicsCosine::scatter) { - t->estimator_ = ESTIMATOR_ANALOG; - } - } else if (filt_type == "spatiallegendre" || filt_type == "zernike" - || filt_type == "zernikeradial") { - t->estimator_ = ESTIMATOR_COLLISION; - } - } - - // ======================================================================= - // READ DATA FOR NUCLIDES - - t->set_nuclides(node_tal); - - // ======================================================================= - // READ DATA FOR SCORES - - t->set_scores(node_tal); - - if (!check_for_node(node_tal, "scores")) { - fatal_error("No scores specified on tally " + std::to_string(t->id_) - + "."); - } - - // Check if tally is compatible with particle type - if (settings::photon_transport) { - if (particle_filter_index == C_NONE) { - for (int score : t->scores_) { - switch (score) { - case SCORE_INVERSE_VELOCITY: - fatal_error("Particle filter must be used with photon " - "transport on and inverse velocity score"); - break; - case SCORE_FLUX: - case SCORE_TOTAL: - case SCORE_SCATTER: - case SCORE_NU_SCATTER: - case SCORE_ABSORPTION: - case SCORE_FISSION: - case SCORE_NU_FISSION: - case SCORE_CURRENT: - case SCORE_EVENTS: - case SCORE_DELAYED_NU_FISSION: - case SCORE_PROMPT_NU_FISSION: - case SCORE_DECAY_RATE: - warning("Particle filter is not used with photon transport" - " on and " + reaction_name(score) + " score."); - break; - } - } - } else { - const auto& f = model::tally_filters[particle_filter_index].get(); - auto pf = dynamic_cast(f); - for (auto p : pf->particles_) { - if (p == Particle::Type::electron || - p == Particle::Type::positron) { - t->estimator_ = ESTIMATOR_ANALOG; - } - } - } - } else { - if (particle_filter_index >= 0) { - const auto& f = model::tally_filters[particle_filter_index].get(); - auto pf = dynamic_cast(f); - for (auto p : pf->particles_) { - if (p != Particle::Type::neutron) { - warning("Particle filter other than NEUTRON used with photon " - "transport turned off. All tallies for particle type " + - std::to_string(static_cast(p)) + " will have no scores"); - } - } - } - } - - // Check for a tally derivative. - if (check_for_node(node_tal, "derivative")) { - int deriv_id = std::stoi(get_node_value(node_tal, "derivative")); - - // Find the derivative with the given id, and store it's index. - auto it = model::tally_deriv_map.find(deriv_id); - if (it == model::tally_deriv_map.end()) { - fatal_error("Could not find derivative " + std::to_string(deriv_id) - + " specified on tally " + std::to_string(t->id_)); - } - - t->deriv_ = it->second; - - // Only analog or collision estimators are supported for differential - // tallies. - if (t->estimator_ == ESTIMATOR_TRACKLENGTH) { - t->estimator_ = ESTIMATOR_COLLISION; - } - - const auto& deriv = model::tally_derivs[t->deriv_]; - if (deriv.variable == DIFF_NUCLIDE_DENSITY - || deriv.variable == DIFF_TEMPERATURE) { - for (int i_nuc : t->nuclides_) { - if (has_energyout && i_nuc == -1) { - fatal_error("Error on tally " + std::to_string(t->id_) - + ": Cannot use a 'nuclide_density' or 'temperature' " - "derivative on a tally with an outgoing energy filter and " - "'total' nuclide rate. Instead, tally each nuclide in the " - "material individually."); - // Note that diff tallies with these characteristics would work - // correctly if no tally events occur in the perturbed material - // (e.g. pertrubing moderator but only tallying fuel), but this - // case would be hard to check for by only reading inputs. - } - } - } - } - - // If settings.xml trigger is turned on, create tally triggers - if (settings::trigger_on) { - t->init_triggers(node_tal); - } - - // ======================================================================= - // SET TALLY ESTIMATOR - - // Check if user specified estimator - if (check_for_node(node_tal, "estimator")) { - std::string est = get_node_value(node_tal, "estimator"); - if (est == "analog") { - t->estimator_ = ESTIMATOR_ANALOG; - } else if (est == "tracklength" || est == "track-length" - || est == "pathlength" || est == "path-length") { - // If the estimator was set to an analog estimator, this means the - // tally needs post-collision information - if (t->estimator_ == ESTIMATOR_ANALOG) { - fatal_error("Cannot use track-length estimator for tally " - + std::to_string(t->id_)); - } - - // Set estimator to track-length estimator - t->estimator_ = ESTIMATOR_TRACKLENGTH; - - } else if (est == "collision") { - // If the estimator was set to an analog estimator, this means the - // tally needs post-collision information - if (t->estimator_ == ESTIMATOR_ANALOG) { - fatal_error("Cannot use collision estimator for tally " + - std::to_string(t->id_)); - } - - // Set estimator to collision estimator - t->estimator_ = ESTIMATOR_COLLISION; - - } else { - fatal_error("Invalid estimator '" + est + "' on tally " + - std::to_string(t->id_)); - } - } + model::tallies.push_back(std::make_unique(node_tal)); } } @@ -1057,7 +1070,7 @@ openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end) if (index_start) *index_start = model::tallies.size(); if (index_end) *index_end = model::tallies.size() + n - 1; for (int i = 0; i < n; ++i) { - model::tallies.push_back(std::make_unique()); + model::tallies.push_back(std::make_unique(-1)); } return 0; } @@ -1141,14 +1154,7 @@ openmc_tally_set_id(int32_t index, int32_t id) return OPENMC_E_OUT_OF_BOUNDS; } - if (model::tally_map.find(id) != model::tally_map.end()) { - set_errmsg("Two or more tallies use the same unique ID: " - + std::to_string(id)); - return OPENMC_E_INVALID_ID; - } - - model::tallies[index]->id_ = id; - model::tally_map[id] = index; + model::tallies[index]->set_id(id); return 0; } @@ -1288,7 +1294,7 @@ openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides) } extern "C" int -openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n) +openmc_tally_get_filters(int32_t index, const int32_t** indices, size_t* n) { if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); @@ -1301,7 +1307,7 @@ openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n) } extern "C" int -openmc_tally_set_filters(int32_t index, int n, const int32_t* indices) +openmc_tally_set_filters(int32_t index, size_t n, const int32_t* indices) { // Make sure the index fits in the array bounds. if (index < 0 || index >= model::tallies.size()) { @@ -1311,9 +1317,15 @@ openmc_tally_set_filters(int32_t index, int n, const int32_t* indices) // Set the filters. try { - model::tallies[index]->set_filters(indices, n); + // Convert indices to filter pointers + std::vector filters; + for (gsl::index i = 0; i < n; ++i) { + int32_t i_filt = indices[i]; + filters.push_back(model::tally_filters.at(i_filt).get()); + } + model::tallies[index]->set_filters(filters); } catch (const std::out_of_range& ex) { - set_errmsg(ex.what()); + set_errmsg("Index in tally filter array out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index a69ab8ca5c..6bf9ca6e9f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -69,7 +69,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, bool end) if (!match.bins_present_) { match.bins_.clear(); match.weights_.clear(); - for (auto i = 0; i < model::tally_filters[i_filt]->n_bins_; ++i) { + for (auto i = 0; i < model::tally_filters[i_filt]->n_bins(); ++i) { match.bins_.push_back(i); match.weights_.push_back(1.0); } @@ -209,14 +209,14 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) if (tally.deriv_ != C_NONE) apply_derivative_to_score(p, i_tally, 0, 0., SCORE_NU_FISSION, score); - if (!settings::run_CE && eo_filt.matches_transport_groups_) { + if (!settings::run_CE && eo_filt.matches_transport_groups()) { // determine outgoing energy group from fission bank auto g_out = static_cast(bank.E); // modify the value so that g_out = 1 corresponds to the highest energy // bin - g_out = eo_filt.n_bins_ - g_out; + g_out = eo_filt.n_bins() - g_out; // change outgoing energy bin simulation::filter_matches[i_eout_filt].bins_[i_bin] = g_out; @@ -231,11 +231,11 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) } // Set EnergyoutFilter bin index - if (E_out < eo_filt.bins_.front() || E_out > eo_filt.bins_.back()) { + if (E_out < eo_filt.bins().front() || E_out > eo_filt.bins().back()) { continue; } else { - auto i_match = lower_bound_index(eo_filt.bins_.begin(), - eo_filt.bins_.end(), E_out); + auto i_match = lower_bound_index(eo_filt.bins().begin(), + eo_filt.bins().end(), E_out); simulation::filter_matches[i_eout_filt].bins_[i_bin] = i_match; } @@ -272,8 +272,8 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) model::tally_filters[i_dg_filt].get())}; // Loop over delayed group bins until the corresponding bin is found - for (auto d_bin = 0; d_bin < dg_filt.n_bins_; ++d_bin) { - if (dg_filt.groups_[d_bin] == g) { + for (auto d_bin = 0; d_bin < dg_filt.n_bins(); ++d_bin) { + if (dg_filt.groups()[d_bin] == g) { // Find the filter index and weight for this filter combination double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { @@ -632,8 +632,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; auto yield = data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->wgt_absorb_ * yield @@ -670,8 +670,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; score = simulation::keff * p->wgt_bank_ / p->n_bank_ * p->n_delayed_bank_[d-1] * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -693,8 +693,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; auto yield = data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->neutron_xs_[i_nuclide].fission * yield @@ -722,8 +722,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; auto yield = data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->neutron_xs_[j_nuclide].fission * yield @@ -770,8 +770,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; @@ -830,8 +830,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Find the corresponding filter bin and then score - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; if (d == g) score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -852,8 +852,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; @@ -892,8 +892,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; @@ -1757,8 +1757,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; score = p->wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= @@ -1807,8 +1807,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; score = simulation::keff * p->wgt_bank_ / p->n_bank_ * p->n_delayed_bank_[d-1] * flux; if (i_nuclide >= 0) { @@ -1839,8 +1839,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; if (i_nuclide >= 0) { score = flux * atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, @@ -1879,8 +1879,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; score = p->wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= @@ -1959,8 +1959,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Find the corresponding filter bin and then score - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; if (d == g) score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1978,8 +1978,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; // Tally each delayed group bin individually - for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { - auto d = filt.groups_[d_bin]; + for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { + auto d = filt.groups()[d_bin]; if (i_nuclide >= 0) { score += atom_density * flux * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 53ad09e906..af86a054d3 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -113,7 +113,7 @@ def test_material(capi_init): m.volume = 10.0 assert m.volume == 10.0 - with pytest.raises(exc.InvalidArgumentError): + with pytest.raises(exc.OpenMCError): m.set_density(1.0, 'goblins') rho = 2.25e-2