crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1286 from paulromano/api-improvements

Browse source 
Improve C++ interface for tallies, filters, materials
This commit is contained in:
Amanda Lund 2019-07-16 13:19:28 -05:00 • committed by GitHub
commit 61c911cffd
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
51 changed files with 1527 additions and 755 deletions
Download patch file Download diff file Expand all files Collapse all files

16
include/openmc/capi.h
View file

@ -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);

96
include/openmc/cell.h
View file

@ -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<double, int32_t>
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<double> rotation_;
std::vector<int32_t> 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<double, int32_t>
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() {}
};
//==============================================================================

74
include/openmc/material.h
View file

@ -6,6 +6,7 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#include <hdf5.h>
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
@ -35,6 +36,7 @@ extern std::unordered_map<int32_t, int32_t> 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<std::string>& name,
const std::vector<double>& 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<const int> nuclides() const { return {nuclide_.data(), nuclide_.size()}; }
//! Get densities of each nuclide in material
//! \return Densities in [atom/b-cm]
gsl::span<const double> 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<Bremsstrahlung> 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_;
};
//==============================================================================

54
include/openmc/tallies/filter.h
View file

@ -7,6 +7,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#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);

17
include/openmc/tallies/filter_azimuthal.h
View file

@ -4,6 +4,8 @@
#include <string>
#include <vector>
#include <gsl/gsl>
#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<double> bins);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<double> bins_;
};

21
include/openmc/tallies/filter_cell.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#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<int32_t>& cells() const { return cells_; }
void set_cells(gsl::span<int32_t> cells);
protected:
//----------------------------------------------------------------------------
// Data members
//! The indices of the cells binned by this filter.
std::vector<int32_t> cells_;

3
include/openmc/tallies/filter_cellborn.h
View file

@ -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)

3
include/openmc/tallies/filter_cellfrom.h
View file

@ -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)

19
include/openmc/tallies/filter_delayedgroup.h
View file

@ -3,6 +3,8 @@
#include <vector>
#include <gsl/gsl>
#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<int>& groups() const { return groups_; }
void set_groups(gsl::span<int> groups);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<int> groups_;
};

19
include/openmc/tallies/filter_distribcell.h
View file

@ -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_;
};

23
include/openmc/tallies/filter_energy.h
View file

@ -3,6 +3,8 @@
#include <vector>
#include <gsl/gsl>
#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<double>& bins() const { return bins_; }
void set_bins(gsl::span<const double> bins);
bool matches_transport_groups() const { return matches_transport_groups_; }
protected:
//----------------------------------------------------------------------------
// Data members
std::vector<double> 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)

10
include/openmc/tallies/filter_energyfunc.h
View file

@ -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<double> energy_;

17
include/openmc/tallies/filter_legendre.h
View file

@ -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_;
};

23
include/openmc/tallies/filter_material.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#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<int32_t>& materials() { return materials_; }
const std::vector<int32_t>& materials() const { return materials_; }
void set_materials(gsl::span<const int32_t> materials);
private:
//----------------------------------------------------------------------------
// Data members
//! The indices of the materials binned by this filter.
std::vector<int32_t> materials_;

12
include/openmc/tallies/filter_mesh.h
View file

@ -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_;
};

6
include/openmc/tallies/filter_meshsurface.h
View file

@ -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;
};

17
include/openmc/tallies/filter_mu.h
View file

@ -3,6 +3,8 @@
#include <vector>
#include <gsl/gsl>
#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<double> bins);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<double> bins_;
};

17
include/openmc/tallies/filter_particle.h
View file

@ -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<Particle::Type>& particles() const { return particles_; }
void set_particles(gsl::span<Particle::Type> particles);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<Particle::Type> particles_;
};

17
include/openmc/tallies/filter_polar.h
View file

@ -4,6 +4,8 @@
#include <cmath>
#include <vector>
#include <gsl/gsl>
#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<double> bins);
private:
//----------------------------------------------------------------------------
// Data members
std::vector<double> bins_;
};

23
include/openmc/tallies/filter_sph_harm.h
View file

@ -3,6 +3,8 @@
#include <string>
#include <gsl/gsl>
#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.

23
include/openmc/tallies/filter_sptl_legendre.h
View file

@ -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.

19
include/openmc/tallies/filter_surface.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#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<int32_t> surfaces);
private:
//----------------------------------------------------------------------------
// Data members
//! The indices of the surfaces binned by this filter.
std::vector<int32_t> surfaces_;

19
include/openmc/tallies/filter_universe.h
View file

@ -5,6 +5,8 @@
#include <unordered_map>
#include <vector>
#include <gsl/gsl>
#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<int32_t> universes);
private:
//----------------------------------------------------------------------------
// Data members
//! The indices of the universes binned by this filter.
std::vector<int32_t> universes_;

32
include/openmc/tallies/filter_zernike.h
View file

@ -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;
};

25
include/openmc/tallies/tally.h
View file

@ -2,8 +2,10 @@
#define OPENMC_TALLIES_TALLY_H
#include "openmc/constants.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/trigger.h"
#include <gsl/gsl>
#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<std::string> scores);
void set_scores(const std::vector<std::string>& scores);
void set_nuclides(pugi::xml_node node);
//----------------------------------------------------------------------------
// Methods for getting and setting filter/stride data.
void set_nuclides(const std::vector<std::string>& nuclides);
const std::vector<int32_t>& 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<Filter*> filters);
int32_t strides(int i) const {return strides_[i];}
@ -112,6 +123,8 @@ private:
std::vector<int32_t> strides_;
int32_t n_filter_bins_ {0};
gsl::index index_;
};
//==============================================================================

4
openmc/capi/error.py
View file

@ -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)

13
openmc/capi/filter.py
View file

@ -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)]

6
openmc/capi/tally.py
View file

@ -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)]

80
src/cell.cpp
View file

@ -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;
}

2
src/geometry_aux.cpp
View file

@ -309,7 +309,7 @@ prepare_distribcell()
for (auto& filt : model::tally_filters) {
auto* distrib_filt = dynamic_cast<DistribcellFilter*>(filt.get());
if (distrib_filt) {
distribcells.insert(distrib_filt->cell_);
distribcells.insert(distrib_filt->cell());
}
}

255
src/material.cpp
View file

@ -47,9 +47,10 @@ std::unordered_map<int32_t, int32_t> 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<std::string>& name,
const std::vector<double>& density)
{
auto n = name.size();
Expects(n > 0);
Expects(n == density.size());
if (n != nuclide_.size()) {
nuclide_.resize(n);
atom_density_ = xt::zeros<double>({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<double, 1> atom_density = xt::zeros<double>({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>(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<double, 1> atom_density = xt::zeros<double>({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<double>({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

6
src/state_point.cpp
View file

@ -141,13 +141,13 @@ openmc_statepoint_write(const char* filename, bool* write_source)
std::vector<int32_t> 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<int32_t> 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);
}

89
src/tallies/filter.cpp
View file

@ -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<AzimuthalFilter>());
@ -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;
}

33
src/tallies/filter_azimuthal.cpp
View file

@ -15,22 +15,35 @@ AzimuthalFilter::from_xml(pugi::xml_node node)
{
auto bins = get_node_array<double>(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<double> 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;

50
src/tallies/filter_cell.cpp
View file

@ -12,29 +12,39 @@ namespace openmc {
void
CellFilter::from_xml(pugi::xml_node node)
{
cells_ = get_node_array<int32_t>(node, "bins");
n_bins_ = cells_.size();
// Get cell IDs and convert into indices into the global cells vector
auto cells = get_node_array<int32_t>(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<int32_t> 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<CellFilter*>(filt);
*cells = cell_filt->cells_.data();
*n = cell_filt->cells_.size();
*cells = cell_filt->cells().data();
*n = cell_filt->cells().size();
return 0;
}

25
src/tallies/filter_delayedgroup.cpp
View file

@ -8,21 +8,32 @@ namespace openmc {
void
DelayedGroupFilter::from_xml(pugi::xml_node node)
{
groups_ = get_node_array<int>(node, "bins");
n_bins_ = groups_.size();
auto groups = get_node_array<int>(node, "bins");
this->set_groups(groups);
}
void
DelayedGroupFilter::set_groups(gsl::span<int> 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

27
src/tallies/filter_distribcell.cpp
View file

@ -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

35
src/tallies/filter_energy.cpp
View file

@ -16,7 +16,25 @@ namespace openmc {
void
EnergyFilter::from_xml(pugi::xml_node node)
{
bins_ = get_node_array<double>(node, "bins");
auto bins = get_node_array<double>(node, "bins");
this->set_bins(bins);
}
void
EnergyFilter::set_bins(gsl::span<const double> 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;
}

16
src/tallies/filter_legendre.cpp
View file

@ -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;
}

57
src/tallies/filter_material.cpp
View file

@ -3,7 +3,6 @@
#include <sstream>
#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<int32_t>(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<int32_t>(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<const int32_t> 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;
}

34
src/tallies/filter_mu.cpp
View file

@ -13,22 +13,36 @@ MuFilter::from_xml(pugi::xml_node node)
{
auto bins = get_node_array<double>(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<double> 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;

19
src/tallies/filter_particle.cpp
View file

@ -8,8 +8,25 @@ void
ParticleFilter::from_xml(pugi::xml_node node)
{
auto particles = get_node_array<int>(node, "bins");
// Convert to vector of Particle::Type
std::vector<Particle::Type> types;
for (auto& p : particles) {
particles_.push_back(static_cast<Particle::Type>(p - 1));
types.push_back(static_cast<Particle::Type>(p - 1));
}
this->set_particles(types);
}
void
ParticleFilter::set_particles(gsl::span<Particle::Type> 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();
}

34
src/tallies/filter_polar.cpp
View file

@ -14,22 +14,36 @@ PolarFilter::from_xml(pugi::xml_node node)
{
auto bins = get_node_array<double>(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<double> 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;

58
src/tallies/filter_sph_harm.cpp
View file

@ -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;

67
src/tallies/filter_sptl_legendre.cpp
View file

@ -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<int>(filt->axis_);
*min = filt->min_;
*max = filt->max_;
*axis = static_cast<int>(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<LegendreAxis>(*axis);
if (min) filt->min_ = *min;
if (max) filt->max_ = *max;
if (axis) filt->set_axis(static_cast<LegendreAxis>(*axis));
if (min && max) filt->set_minmax(*min, *max);
return 0;
}

41
src/tallies/filter_surface.cpp
View file

@ -11,30 +11,41 @@ namespace openmc {
void
SurfaceFilter::from_xml(pugi::xml_node node)
{
surfaces_ = get_node_array<int32_t>(node, "bins");
n_bins_ = surfaces_.size();
}
auto surfaces = get_node_array<int32_t>(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<int32_t> 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

41
src/tallies/filter_universe.cpp
View file

@ -11,30 +11,39 @@ namespace openmc {
void
UniverseFilter::from_xml(pugi::xml_node node)
{
universes_ = get_node_array<int32_t>(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<int32_t>(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<int32_t> 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

17
src/tallies/filter_zernike.cpp
View file

@ -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;
}

594
src/tallies/tally.cpp
View file

@ -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 <filter> element. The <tally> element itself should "
"have a list of filters that apply, e.g., <filters>1 2</filters> "
"where 1 and 2 are the IDs of filters specified outside of "
"<tally>."};
}
// Determine number of filters
std::vector<int> filter_ids;
if (check_for_node(node, "filters")) {
filter_ids = get_node_array<int>(node, "filters");
}
// Allocate and store filter user ids
std::vector<Filter*> 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<ParticleFilter*>(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<SphericalHarmonicsFilter*>(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<ParticleFilter*>(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<ParticleFilter*>(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<int>(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<Tally>(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<Filter*> 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<const EnergyoutFilter*>(filt)) {
if (dynamic_cast<const EnergyoutFilter*>(f)) {
energyout_filter_ = i;
} else if (dynamic_cast<const DelayedGroupFilter*>(filt)) {
} else if (dynamic_cast<const DelayedGroupFilter*>(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<std::string> scores)
Tally::set_scores(const std::vector<std::string>& 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<std::string>(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<std::string>& 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<Tally>());
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 <filter> element. The <tally> element itself should "
"have a list of filters that apply, e.g., <filters>1 2</filters> "
"where 1 and 2 are the IDs of filters specified outside of "
"<tally>.");
}
// Determine number of filters
std::vector<int> filters;
if (check_for_node(node_tal, "filters")) {
filters = get_node_array<int>(node_tal, "filters");
}
// Allocate and store filter user ids
if (!filters.empty()) {
std::vector<int> 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<ParticleFilter*>(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<SphericalHarmonicsFilter*>(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<ParticleFilter*>(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<ParticleFilter*>(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<int>(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<Tally>(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<Tally>());
model::tallies.push_back(std::make_unique<Tally>(-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<Filter*> 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;
}

72
src/tallies/tally_scoring.cpp
View file

@ -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<int>(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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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<DelayedGroupFilter*>(
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,

2
tests/unit_tests/test_capi.py
View file

@ -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