diff --git a/CMakeLists.txt b/CMakeLists.txt index c16d72662..4ac9d9c5d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -373,26 +373,15 @@ add_library(libopenmc SHARED src/tallies/tally_derivative_header.F90 src/tallies/tally_filter.F90 src/tallies/tally_filter_header.F90 - src/tallies/tally_filter_azimuthal.F90 - src/tallies/tally_filter_cell.F90 - src/tallies/tally_filter_cellborn.F90 - src/tallies/tally_filter_cellfrom.F90 + src/tallies/tally_filter_cpp.F90 src/tallies/tally_filter_delayedgroup.F90 src/tallies/tally_filter_distribcell.F90 src/tallies/tally_filter_energy.F90 - src/tallies/tally_filter_energyfunc.F90 src/tallies/tally_filter_legendre.F90 - src/tallies/tally_filter_material.F90 src/tallies/tally_filter_mesh.F90 src/tallies/tally_filter_meshsurface.F90 - src/tallies/tally_filter_mu.F90 src/tallies/tally_filter_particle.F90 - src/tallies/tally_filter_polar.F90 src/tallies/tally_filter_sph_harm.F90 - src/tallies/tally_filter_sptl_legendre.F90 - src/tallies/tally_filter_surface.F90 - src/tallies/tally_filter_universe.F90 - src/tallies/tally_filter_zernike.F90 src/tallies/tally_header.F90 src/tallies/trigger.F90 src/tallies/trigger_header.F90 @@ -442,6 +431,25 @@ add_library(libopenmc SHARED src/string_functions.cpp src/summary.cpp src/surface.cpp + src/tallies/filter.cpp + src/tallies/filter_azimuthal.cpp + src/tallies/filter_cellborn.cpp + src/tallies/filter_cellfrom.cpp + src/tallies/filter_cell.cpp + src/tallies/filter_distribcell.cpp + src/tallies/filter_energyfunc.cpp + src/tallies/filter_legendre.cpp + src/tallies/filter_material.cpp + src/tallies/filter_mesh.cpp + src/tallies/filter_meshsurface.cpp + src/tallies/filter_mu.cpp + src/tallies/filter_polar.cpp + src/tallies/filter_sph_harm.cpp + src/tallies/filter_sptl_legendre.cpp + src/tallies/filter_surface.cpp + src/tallies/filter_universe.cpp + src/tallies/filter_zernike.cpp + src/tallies/tally.cpp src/timer.cpp src/thermal.cpp src/xml_interface.cpp diff --git a/docs/source/io_formats/data_wmp.rst b/docs/source/io_formats/data_wmp.rst index 25916c522..8d7cf6ba0 100644 --- a/docs/source/io_formats/data_wmp.rst +++ b/docs/source/io_formats/data_wmp.rst @@ -4,10 +4,15 @@ Windowed Multipole Library Format ================================= -**/version** (*char[]*) - The format version of the file. The current version is "v1.0" +**/** + +:Attributes: - **filetype** (*char[]*) -- String indicating the type of file + - **version** (*int[2]*) -- Major and minor version of the data + +**//** + +:Datasets: -**/nuclide/** - **broaden_poly** (*int[]*) If 1, Doppler broaden curve fit for window with corresponding index. If 0, do not. diff --git a/examples/jupyter/nuclear-data.ipynb b/examples/jupyter/nuclear-data.ipynb index 91923faa1..f1c4c2b2a 100644 --- a/examples/jupyter/nuclear-data.ipynb +++ b/examples/jupyter/nuclear-data.ipynb @@ -1466,7 +1466,7 @@ }, "outputs": [], "source": [ - "url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/092238.h5'\n", + "url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/092238.h5'\n", "filename, headers = urllib.request.urlretrieve(url, '092238.h5')" ] }, diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 5441c1026..3ead7d833 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -136,7 +136,6 @@ extern "C" { // Global variables extern char openmc_err_msg[256]; extern int32_t n_cells; - extern int32_t n_filters; extern int32_t n_lattices; extern int32_t n_materials; extern int n_nuclides; diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 361cfc158..d249e4cb5 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -62,6 +62,7 @@ constexpr int MAX_SAMPLE {100000}; // Maximum number of words in a single line, length of line, and length of // single word +constexpr int MAX_LINE_LEN {250}; constexpr int MAX_WORD_LEN {150}; // Maximum number of external source spatial resamples to encounter before an diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index ac68f5b00..da55c3d9b 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -5,6 +5,7 @@ #define OPENMC_GEOMETRY_AUX_H #include +#include namespace openmc { @@ -26,7 +27,7 @@ extern "C" void assign_temperatures(); //! //! This function looks for a universe that is not listed in a Cell::fill or in //! a Lattice. -//! @return The index of the root universe. +//! \return The index of the root universe. //============================================================================== extern "C" int32_t find_root_universe(); @@ -41,14 +42,14 @@ extern "C" void neighbor_lists(); //! Populate all data structures needed for distribcells. //============================================================================== -extern "C" void prepare_distribcell(int32_t* filter_cell_list, int n); +extern "C" void prepare_distribcell(); //============================================================================== //! Recursively search through the geometry and count cell instances. //! //! This function will update the Cell::n_instances value for each cell in the //! geometry. -//! @param univ_indx The index of the universe to begin searching from (probably +//! \param univ_indx The index of the universe to begin searching from (probably //! the root universe). //============================================================================== @@ -56,51 +57,33 @@ extern "C" void count_cell_instances(int32_t univ_indx); //============================================================================== //! Recursively search through universes and count universe instances. -//! @param search_univ The index of the universe to begin searching from. -//! @param target_univ_id The ID of the universe to be counted. -//! @return The number of instances of target_univ_id in the geometry tree under +//! \param search_univ The index of the universe to begin searching from. +//! \param target_univ_id The ID of the universe to be counted. +//! \return The number of instances of target_univ_id in the geometry tree under //! search_univ. //============================================================================== extern "C" int count_universe_instances(int32_t search_univ, int32_t target_univ_id); -//============================================================================== -//! Find the length necessary for a string to contain a distribcell path. -//! @param target_cell The index of the Cell in the global Cell array. -//! @param map The index of the distribcell mapping corresponding to the target -//! cell. -//! @param target_offset An instance number for a distributed cell. -//! @param root_univ The index of the root Universe in the global Universe -//! array. -//! @return The size of a character array needed to fit the distribcell path. -//============================================================================== - -extern "C" int -distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ); - //============================================================================== //! Build a character array representing the path to a distribcell instance. -//! @param target_cell The index of the Cell in the global Cell array. -//! @param map The index of the distribcell mapping corresponding to the target +//! \param target_cell The index of the Cell in the global Cell array. +//! \param map The index of the distribcell mapping corresponding to the target //! cell. -//! @param target_offset An instance number for a distributed cell. -//! @param root_univ The index of the root Universe in the global Universe -//! array. -//! @param[out] path The unique traversal through the geometry tree that leads -//! to the desired instance of the target cell. +//! \param target_offset An instance number for a distributed cell. +//! \return The unique traversal through the geometry tree that leads to the +//! desired instance of the target cell. //============================================================================== -extern "C" void -distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ, char *path); +std::string +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset); //============================================================================== //! Determine the maximum number of nested coordinate levels in the geometry. -//! @param univ The index of the universe to begin seraching from (probably the +//! \param univ The index of the universe to begin seraching from (probably the //! root universe). -//! @return The number of coordinate levels. +//! \return The number of coordinate levels. //============================================================================== extern "C" int maximum_levels(int32_t univ); diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 690cb5ae8..8b80f29d5 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -380,5 +380,11 @@ write_dataset(hid_t obj_id, const char* name, Position r) write_dataset(obj_id, name, buffer); } +inline void +write_dataset(hid_t obj_id, const char* name, std::string buffer) +{ + write_string(obj_id, name, buffer.c_str(), false); +} + } // namespace openmc #endif // OPENMC_HDF5_INTERFACE_H diff --git a/include/openmc/material.h b/include/openmc/material.h index 909215ee3..b6d7e19cd 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -24,7 +24,7 @@ extern std::unordered_map material_map; class Material { public: - int32_t id; //!< Unique ID + int32_t id_; //!< Unique ID double volume_ {-1.0}; //!< Volume in [cm^3] bool fissionable {false}; //!< Does this material contain fissionable nuclides diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index 8bfa11899..ca0646961 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -90,7 +90,7 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]); //! evaluated at rho and phi. //============================================================================== -extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); +extern "C" void calc_zn(int n, double rho, double phi, double zn[]); //============================================================================== //! Calculate only the even radial components of n-th order modified Zernike @@ -113,7 +113,7 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); //! evaluated at rho and phi when m = 0. //============================================================================== -extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]); +extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]); //============================================================================== //! Rotate the direction cosines through a polar angle whose cosine is mu and diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h new file mode 100644 index 000000000..7920191aa --- /dev/null +++ b/include/openmc/tallies/filter.h @@ -0,0 +1,91 @@ +#ifndef OPENMC_TALLIES_FILTER_H +#define OPENMC_TALLIES_FILTER_H + +#include +#include +#include + +#include "openmc/hdf5_interface.h" +#include "openmc/particle.h" +#include "pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +extern "C" int32_t n_filters; + +class FilterMatch; +extern std::vector filter_matches; +#pragma omp threadprivate(filter_matches) + +class Filter; +extern std::vector tally_filters; + +//============================================================================== +//! Stores bins and weights for filtered tally events. +//============================================================================== + +class FilterMatch +{ +public: + std::vector bins_; + std::vector weights_; +}; + +//============================================================================== +//! Modifies tally score events. +//============================================================================== + +class Filter +{ +public: + virtual ~Filter() = default; + + virtual std::string type() const = 0; + + //! Uses an XML input to fill the filter's data fields. + virtual void from_xml(pugi::xml_node node) = 0; + + //! Matches a tally event to a set of filter bins and weights. + //! + //! \param[out] match will contain the matching bins and corresponding + //! weights; note that there may be zero matching bins + virtual void + get_all_bins(const Particle* p, int estimator, FilterMatch& match) const = 0; + + //! Writes data describing this filter to an HDF5 statepoint group. + virtual void + to_statepoint(hid_t filter_group) const + { + write_dataset(filter_group, "type", type()); + write_dataset(filter_group, "n_bins", n_bins_); + } + + //! Return a string describing a filter bin for the tallies.out file. + // + //! For example, an `EnergyFilter` might return the string + //! "Incoming Energy [0.625E-6, 20.0)". + virtual std::string text_label(int bin) const = 0; + + virtual void initialize() {} + + int n_bins_; +}; + +//============================================================================== + +extern "C" void free_memory_tally_c(); + +//============================================================================== + +// Filter-related Fortran functions that will be called from C++ +extern "C" int verify_filter(int32_t index); +extern "C" Filter* filter_from_f(int32_t index); +extern "C" void filter_update_n_bins(int32_t index); + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_H diff --git a/include/openmc/tallies/filter_azimuthal.h b/include/openmc/tallies/filter_azimuthal.h new file mode 100644 index 000000000..b89640031 --- /dev/null +++ b/include/openmc/tallies/filter_azimuthal.h @@ -0,0 +1,35 @@ +#ifndef OPENMC_TALLIES_FILTER_AZIMUTHAL_H +#define OPENMC_TALLIES_FILTER_AZIMUTHAL_H + +#include +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Bins the incident neutron azimuthal angle (relative to the global xy-plane). +//============================================================================== + +class AzimuthalFilter : public Filter +{ +public: + ~AzimuthalFilter() = default; + + std::string type() const override {return "azimuthal";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + std::vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_AZIMUTHAL_H diff --git a/include/openmc/tallies/filter_cell.h b/include/openmc/tallies/filter_cell.h new file mode 100644 index 000000000..9c654e82e --- /dev/null +++ b/include/openmc/tallies/filter_cell.h @@ -0,0 +1,42 @@ +#ifndef OPENMC_TALLIES_FILTER_CELL_H +#define OPENMC_TALLIES_FILTER_CELL_H + +#include +#include +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Specifies which geometric cells tally events reside in. +//============================================================================== + +class CellFilter : public Filter +{ +public: + ~CellFilter() = default; + + 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; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //! The indices of the cells binned by this filter. + std::vector cells_; + + //! A map from cell indices to filter bin indices. + std::unordered_map map_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_CELL_H diff --git a/include/openmc/tallies/filter_cellborn.h b/include/openmc/tallies/filter_cellborn.h new file mode 100644 index 000000000..400e82c28 --- /dev/null +++ b/include/openmc/tallies/filter_cellborn.h @@ -0,0 +1,26 @@ +#ifndef OPENMC_TALLIES_FILTER_CELLBORN_H +#define OPENMC_TALLIES_FILTER_CELLBORN_H + +#include + +#include "openmc/tallies/filter_cell.h" + +namespace openmc { + +//============================================================================== +//! Specifies which cell the particle was born in. +//============================================================================== + +class CellbornFilter : public CellFilter +{ +public: + std::string type() const override {return "cellborn";} + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + std::string text_label(int bin) const override; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_CELLBORN_H diff --git a/include/openmc/tallies/filter_cellfrom.h b/include/openmc/tallies/filter_cellfrom.h new file mode 100644 index 000000000..bd3e08dcb --- /dev/null +++ b/include/openmc/tallies/filter_cellfrom.h @@ -0,0 +1,26 @@ +#ifndef OPENMC_TALLIES_FILTER_CELLFROM_H +#define OPENMC_TALLIES_FILTER_CELLFROM_H + +#include + +#include "openmc/tallies/filter_cell.h" + +namespace openmc { + +//============================================================================== +//! Specifies which geometric cells particles exit when crossing a surface. +//============================================================================== + +class CellFromFilter : public CellFilter +{ +public: + std::string type() const override {return "cellfrom";} + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + std::string text_label(int bin) const override; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_CELLFROM_H diff --git a/include/openmc/tallies/filter_distribcell.h b/include/openmc/tallies/filter_distribcell.h new file mode 100644 index 000000000..73857e7cb --- /dev/null +++ b/include/openmc/tallies/filter_distribcell.h @@ -0,0 +1,36 @@ +#ifndef OPENMC_TALLIES_FILTER_DISTRIBCELL_H +#define OPENMC_TALLIES_FILTER_DISTRIBCELL_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Specifies which distributed geometric cells tally events reside in. +//============================================================================== + +class DistribcellFilter : public Filter +{ +public: + ~DistribcellFilter() = default; + + 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; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + int32_t cell_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_DISTRIBCELL_H diff --git a/include/openmc/tallies/filter_energyfunc.h b/include/openmc/tallies/filter_energyfunc.h new file mode 100644 index 000000000..6f5182626 --- /dev/null +++ b/include/openmc/tallies/filter_energyfunc.h @@ -0,0 +1,45 @@ +#ifndef OPENMC_TALLIES_FILTER_ENERGYFUNC_H +#define OPENMC_TALLIES_FILTER_ENERGYFUNC_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Multiplies tally scores by an arbitrary function of incident energy +//! described by a piecewise linear-linear interpolation. +//============================================================================== + +class EnergyFunctionFilter : public Filter +{ +public: + EnergyFunctionFilter() + : Filter {} + { + n_bins_ = 1; + } + + ~EnergyFunctionFilter() = default; + + std::string type() const override {return "energyfunction";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //! Incident neutron energy interpolation grid. + std::vector energy_; + + //! Interpolant values. + std::vector y_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_ENERGYFUNC_H diff --git a/include/openmc/tallies/filter_legendre.h b/include/openmc/tallies/filter_legendre.h new file mode 100644 index 000000000..054ba14e7 --- /dev/null +++ b/include/openmc/tallies/filter_legendre.h @@ -0,0 +1,34 @@ +#ifndef OPENMC_TALLIES_FILTER_LEGENDRE_H +#define OPENMC_TALLIES_FILTER_LEGENDRE_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Gives Legendre moments of the change in scattering angle +//============================================================================== + +class LegendreFilter : public Filter +{ +public: + ~LegendreFilter() = default; + + std::string type() const override {return "legendre";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + int order_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_LEGENDRE_H diff --git a/include/openmc/tallies/filter_material.h b/include/openmc/tallies/filter_material.h new file mode 100644 index 000000000..d272451db --- /dev/null +++ b/include/openmc/tallies/filter_material.h @@ -0,0 +1,42 @@ +#ifndef OPENMC_TALLIES_FILTER_MATERIAL_H +#define OPENMC_TALLIES_FILTER_MATERIAL_H + +#include +#include +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Specifies which material tally events reside in. +//============================================================================== + +class MaterialFilter : public Filter +{ +public: + ~MaterialFilter() = default; + + 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; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //! The indices of the materials binned by this filter. + std::vector materials_; + + //! A map from material indices to filter bin indices. + std::unordered_map map_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_MATERIAL_H diff --git a/include/openmc/tallies/filter_mesh.h b/include/openmc/tallies/filter_mesh.h new file mode 100644 index 000000000..a6caa8f56 --- /dev/null +++ b/include/openmc/tallies/filter_mesh.h @@ -0,0 +1,41 @@ +#ifndef OPENMC_TALLIES_FILTER_MESH_H +#define OPENMC_TALLIES_FILTER_MESH_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Indexes the location of particle events to a regular mesh. For tracklength +//! tallies, it will produce multiple valid bins and the bin weight will +//! correspond to the fraction of the track length that lies in that bin. +//============================================================================== + +class MeshFilter : public Filter +{ +public: + ~MeshFilter() = default; + + std::string type() const override {return "mesh";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + virtual int32_t mesh() const {return mesh_;} + + virtual void set_mesh(int32_t mesh); + +protected: + int32_t mesh_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_MESH_H diff --git a/include/openmc/tallies/filter_meshsurface.h b/include/openmc/tallies/filter_meshsurface.h new file mode 100644 index 000000000..32393cfac --- /dev/null +++ b/include/openmc/tallies/filter_meshsurface.h @@ -0,0 +1,22 @@ +#ifndef OPENMC_TALLIES_FILTER_MESHSURFACE_H +#define OPENMC_TALLIES_FILTER_MESHSURFACE_H + +#include "openmc/tallies/filter_mesh.h" + +namespace openmc { + +class MeshSurfaceFilter : public MeshFilter +{ +public: + std::string type() const override {return "meshsurface";} + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + std::string text_label(int bin) const override; + + void set_mesh(int32_t mesh) override; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_MESHSURFACE_H diff --git a/include/openmc/tallies/filter_mu.h b/include/openmc/tallies/filter_mu.h new file mode 100644 index 000000000..68b2f7b02 --- /dev/null +++ b/include/openmc/tallies/filter_mu.h @@ -0,0 +1,35 @@ +#ifndef OPENMC_TALLIES_FILTER_MU_H +#define OPENMC_TALLIES_FILTER_MU_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Bins the incoming-outgoing direction cosine. This is only used for scatter +//! reactions. +//============================================================================== + +class MuFilter : public Filter +{ +public: + ~MuFilter() = default; + + std::string type() const override {return "mu";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + std::vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_MU_H diff --git a/include/openmc/tallies/filter_polar.h b/include/openmc/tallies/filter_polar.h new file mode 100644 index 000000000..2bcdcd07f --- /dev/null +++ b/include/openmc/tallies/filter_polar.h @@ -0,0 +1,35 @@ +#ifndef OPENMC_TALLIES_FILTER_POLAR_H +#define OPENMC_TALLIES_FILTER_POLAR_H + +#include +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Bins the incident neutron polar angle (relative to the global z-axis). +//============================================================================== + +class PolarFilter : public Filter +{ +public: + ~PolarFilter() = default; + + std::string type() const override {return "polar";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + std::vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_POLAR_H diff --git a/include/openmc/tallies/filter_sph_harm.h b/include/openmc/tallies/filter_sph_harm.h new file mode 100644 index 000000000..786ef1a58 --- /dev/null +++ b/include/openmc/tallies/filter_sph_harm.h @@ -0,0 +1,42 @@ +#ifndef OPENMC_TALLIES_FILTER_SPH_HARM_H +#define OPENMC_TALLIES_FILTER_SPH_HARM_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//TODO: those integer values are not needed when Fortran interop is removed +enum class SphericalHarmonicsCosine { + scatter = 1, particle = 2 +}; + +//============================================================================== +//! Gives spherical harmonics expansion moments of a tally score +//============================================================================== + +class SphericalHarmonicsFilter : public Filter +{ +public: + ~SphericalHarmonicsFilter() = default; + + std::string type() const override {return "sphericalharmonics";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + int order_; + + //! The type of angle that this filter measures when binning events. + SphericalHarmonicsCosine cosine_ {SphericalHarmonicsCosine::particle}; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_SPH_HARM_H diff --git a/include/openmc/tallies/filter_sptl_legendre.h b/include/openmc/tallies/filter_sptl_legendre.h new file mode 100644 index 000000000..1b94ecf3d --- /dev/null +++ b/include/openmc/tallies/filter_sptl_legendre.h @@ -0,0 +1,48 @@ +#ifndef OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H +#define OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//TODO: those integer values are not needed when Fortran interop is removed +enum class LegendreAxis { + x = 1, y = 2, z = 3 +}; + +//============================================================================== +//! Gives Legendre moments of the particle's normalized position along an axis +//============================================================================== + +class SpatialLegendreFilter : public Filter +{ +public: + ~SpatialLegendreFilter() = default; + + std::string type() const override {return "spatiallegendre";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + int order_; + + //! The Cartesian coordinate axis that the Legendre expansion is applied to. + LegendreAxis axis_; + + //! The minimum coordinate along the reference axis that the expansion covers. + double min_; + + //! The maximum coordinate along the reference axis that the expansion covers. + double max_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H diff --git a/include/openmc/tallies/filter_surface.h b/include/openmc/tallies/filter_surface.h new file mode 100644 index 000000000..ea4dc5b99 --- /dev/null +++ b/include/openmc/tallies/filter_surface.h @@ -0,0 +1,42 @@ +#ifndef OPENMC_TALLIES_FILTER_SURFACE_H +#define OPENMC_TALLIES_FILTER_SURFACE_H + +#include +#include +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Specifies which surface particles are crossing +//============================================================================== + +class SurfaceFilter : public Filter +{ +public: + ~SurfaceFilter() = default; + + 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; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //! The indices of the surfaces binned by this filter. + std::vector surfaces_; + + //! A map from surface indices to filter bin indices. + std::unordered_map map_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_SURFACE_H diff --git a/include/openmc/tallies/filter_universe.h b/include/openmc/tallies/filter_universe.h new file mode 100644 index 000000000..3ae0093cb --- /dev/null +++ b/include/openmc/tallies/filter_universe.h @@ -0,0 +1,42 @@ +#ifndef OPENMC_TALLIES_FILTER_UNIVERSE_H +#define OPENMC_TALLIES_FILTER_UNIVERSE_H + +#include +#include +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Specifies which geometric universes tally events reside in. +//============================================================================== + +class UniverseFilter : public Filter +{ +public: + ~UniverseFilter() = default; + + 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; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //! The indices of the universes binned by this filter. + std::vector universes_; + + //! A map from universe indices to filter bin indices. + std::unordered_map map_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_UNIVERSE_H diff --git a/include/openmc/tallies/filter_zernike.h b/include/openmc/tallies/filter_zernike.h new file mode 100644 index 000000000..19fff1336 --- /dev/null +++ b/include/openmc/tallies/filter_zernike.h @@ -0,0 +1,65 @@ +#ifndef OPENMC_TALLIES_FILTER_ZERNIKE_H +#define OPENMC_TALLIES_FILTER_ZERNIKE_H + +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Gives Zernike polynomial moments of a particle's position +//============================================================================== + +class ZernikeFilter : public Filter +{ +public: + std::string type() const override {return "zernike";} + + ~ZernikeFilter() = default; + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + int order() const {return order_;} + + virtual void set_order(int order); + + //! Cartesian x coordinate for the origin of this expansion. + double x_; + + //! Cartesian y coordinate for the origin of this expansion. + double y_; + + //! Maximum radius from the origin covered by this expansion. + double r_; + +protected: + int order_; +}; + +//============================================================================== +//! Gives even order radial Zernike polynomial moments of a particle's position +//============================================================================== + +class ZernikeRadialFilter : public ZernikeFilter +{ +public: + std::string type() const override {return "zernikeradial";} + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + std::string text_label(int bin) const override; + + void set_order(int order) override; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_ZERNIKE_H diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 8d9d794d1..5f6e74bbc 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -238,13 +238,15 @@ class MaterialFilter(Filter): materials = POINTER(c_int32)() n = c_int32() _dll.openmc_material_filter_get_bins(self._index, materials, n) - return [Material(index=materials[i]) for i in range(n.value)] + #TODO: fix this off-by-one when materials become 0-indexed + return [Material(index=materials[i]+1) for i in range(n.value)] @bins.setter def bins(self, materials): # Get material indices as int32_t[] n = len(materials) - bins = (c_int32*n)(*(m._index for m in materials)) + #TODO: fix this off-by-one when materials become 0-indexed + bins = (c_int32*n)(*(m._index-1 for m in materials)) _dll.openmc_material_filter_set_bins(self._index, n, bins) diff --git a/openmc/capi/math.py b/openmc/capi/math.py index a8ea849b3..34beb585f 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -18,11 +18,11 @@ _dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double] _dll.calc_rn_c.restype = None _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)] -_dll.calc_zn_c.restype = None -_dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)] +_dll.calc_zn.restype = None +_dll.calc_zn.argtypes = [c_int, c_double, c_double, ndpointer(c_double)] -_dll.calc_zn_rad_c.restype = None -_dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)] +_dll.calc_zn_rad.restype = None +_dll.calc_zn_rad.argtypes = [c_int, c_double, ndpointer(c_double)] _dll.rotate_angle_c.restype = None _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, @@ -153,7 +153,7 @@ def calc_zn(n, rho, phi): num_bins = ((n + 1) * (n + 2)) // 2 zn = np.zeros(num_bins, dtype=np.float64) - _dll.calc_zn_c(n, rho, phi, zn) + _dll.calc_zn(n, rho, phi, zn) return zn @@ -179,7 +179,7 @@ def calc_zn_rad(n, rho): num_bins = n // 2 + 1 zn_rad = np.zeros(num_bins, dtype=np.float64) - _dll.calc_zn_rad_c(n, rho, zn_rad) + _dll.calc_zn_rad(n, rho, zn_rad) return zn_rad diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index dec11cd3a..9aff59444 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -4,7 +4,9 @@ HDF5_VERSION_MINOR = 0 HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR) # Version of WMP nuclear data format -WMP_VERSION = 'v1.0' +WMP_VERSION_MAJOR = 1 +WMP_VERSION_MINOR = 1 +WMP_VERSION = (WMP_VERSION_MAJOR, WMP_VERSION_MINOR) from .data import * diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index a0d90ac01..5d34c9ff6 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -4,7 +4,7 @@ from math import exp, erf, pi, sqrt import h5py import numpy as np -from . import WMP_VERSION +from . import WMP_VERSION, WMP_VERSION_MAJOR from .data import K_BOLTZMANN import openmc.checkvalue as cv from openmc.mixin import EqualityMixin @@ -133,6 +133,8 @@ class WindowedMultipole(EqualityMixin): Parameters ---------- + name : str + Name of the nuclide using the GND naming convention Attributes ---------- @@ -171,7 +173,8 @@ class WindowedMultipole(EqualityMixin): a/E + b/sqrt(E) + c + d sqrt(E) + ... """ - def __init__(self): + def __init__(self, name): + self.name = name self.spacing = None self.sqrtAWR = None self.E_min = None @@ -181,6 +184,10 @@ class WindowedMultipole(EqualityMixin): self.broaden_poly = None self.curvefit = None + @property + def name(self): + return self._name + @property def fit_order(self): return self.curvefit.shape[1] - 1 @@ -221,6 +228,11 @@ class WindowedMultipole(EqualityMixin): def curvefit(self): return self._curvefit + @name.setter + def name(self, name): + cv.check_type('name', name, str) + self._name = name + @spacing.setter def spacing(self, spacing): if spacing is not None: @@ -322,17 +334,25 @@ class WindowedMultipole(EqualityMixin): group = group_or_filename else: h5file = h5py.File(group_or_filename, 'r') - try: - version = h5file['version'].value.decode() - except AttributeError: - version = h5file['version'].value[0].decode() - if version != WMP_VERSION: - raise ValueError('The given WMP data uses version ' - + version + ' whereas your installation of the OpenMC ' - 'Python API expects version ' + WMP_VERSION) - group = h5file['nuclide'] - out = cls() + # Make sure version matches + if 'version' in h5file.attrs: + major, minor = h5file.attrs['version'] + if major != WMP_VERSION_MAJOR: + raise IOError( + 'WMP data format uses version {}. {} whereas your ' + 'installation of the OpenMC Python API expects version ' + '{}.x.'.format(major, minor, WMP_VERSION_MAJOR)) + else: + raise IOError( + 'WMP data does not indicate a version. Your installation of ' + 'the OpenMC Python API expects version {}.x data.' + .format(WMP_VERSION_MAJOR)) + + group = list(h5file.values())[0] + + name = group.name[1:] + out = cls(name) # Read scalars. @@ -478,13 +498,16 @@ class WindowedMultipole(EqualityMixin): fun = np.vectorize(lambda x: self._evaluate(x, T)) return fun(E) - def export_to_hdf5(self, path, libver='earliest'): + def export_to_hdf5(self, path, mode='a', libver='earliest'): """Export windowed multipole data to an HDF5 file. Parameters ---------- path : str Path to write HDF5 file to + mode : {'r', r+', 'w', 'x', 'a'} + Mode that is used to open the HDF5 file. This is the second argument + to the :class:`h5py.File` constructor. libver : {'earliest', 'latest'} Compatibility mode for the HDF5 file. 'latest' will produce files that are less backwards compatible but have performance benefits. @@ -492,12 +515,11 @@ class WindowedMultipole(EqualityMixin): """ # Open file and write version. - with h5py.File(path, 'w', libver=libver) as f: - f.create_dataset('version', (1, ), dtype='S10') - f['version'][:] = WMP_VERSION.encode('ASCII') + with h5py.File(path, mode, libver=libver) as f: + f.attrs['filetype'] = np.string_('data_wmp') + f.attrs['version'] = np.array(WMP_VERSION) - # Make a nuclide group. - g = f.create_group('nuclide') + g = f.create_group(self.name) # Write scalars. g.create_dataset('spacing', data=np.array(self.spacing)) diff --git a/openmc/polynomial.py b/openmc/polynomial.py index 8bcb74f68..ac1ad675d 100644 --- a/openmc/polynomial.py +++ b/openmc/polynomial.py @@ -1,6 +1,5 @@ import numpy as np import openmc -import openmc.capi as capi from collections.abc import Iterable @@ -74,6 +73,7 @@ class ZernikeRadial(Polynomial): return self._order def __call__(self, r): + import openmc.capi as capi if isinstance(r, Iterable): return [np.sum(self._norm_coef * capi.calc_zn_rad(self.order, r_i)) for r_i in r] diff --git a/scripts/openmc-get-multipole-data b/scripts/openmc-get-multipole-data index c8a2d1e9b..554a18428 100755 --- a/scripts/openmc-get-multipole-data +++ b/scripts/openmc-get-multipole-data @@ -29,9 +29,9 @@ parser.add_argument('-b', '--batch', action='store_true', args = parser.parse_args() -baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/' -files = ['WMP_Library_v1.0.tar.gz'] -checksums = ['22cb675734cfccb278dffd40dcfbf26a'] +baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/' +files = ['WMP_Library_v1.1.tar.gz'] +checksums = ['8523895928dd6ba63fba803e3a45d4f3'] block_size = 16384 # ============================================================================== diff --git a/scripts/openmc-plot-mesh-tally b/scripts/openmc-plot-mesh-tally index 0dfb29b9b..61ada852c 100755 --- a/scripts/openmc-plot-mesh-tally +++ b/scripts/openmc-plot-mesh-tally @@ -11,8 +11,10 @@ import tkinter.font as font import tkinter.messagebox as messagebox import tkinter.ttk as ttk +import matplotlib +matplotlib.use("TkAgg") from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg -from matplotlib.backends.backend_tkagg import NavigationToolbar2TkAgg +from matplotlib.backends.backend_tkagg import NavigationToolbar2Tk from matplotlib.figure import Figure import matplotlib.pyplot as plt import numpy as np @@ -56,7 +58,7 @@ class MeshPlotter(tk.Frame): self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1) # Create the navigation toolbar, tied to the canvas - self.mpl_toolbar = NavigationToolbar2TkAgg(self.canvas, figureFrame) + self.mpl_toolbar = NavigationToolbar2Tk(self.canvas, figureFrame) self.mpl_toolbar.update() self.canvas._tkcanvas.pack(side=tk.TOP, fill=tk.BOTH, expand=1) diff --git a/src/cell.cpp b/src/cell.cpp index f851922e1..65582d780 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -476,7 +476,7 @@ CSGCell::to_hdf5(hid_t cell_group) const std::vector mat_ids; for (auto i_mat : material_) { if (i_mat != MATERIAL_VOID) { - mat_ids.push_back(materials[i_mat]->id); + mat_ids.push_back(materials[i_mat]->id_); } else { mat_ids.push_back(MATERIAL_VOID); } @@ -649,6 +649,19 @@ read_cells(pugi::xml_node* node) cells.push_back(new CSGCell(cell_node)); } + // Fill the cell map. + for (int i = 0; i < cells.size(); i++) { + int32_t id = cells[i]->id_; + auto search = cell_map.find(id); + if (search == cell_map.end()) { + cell_map[id] = i; + } else { + std::stringstream err_msg; + err_msg << "Two or more cells use the same unique ID: " << id; + fatal_error(err_msg); + } + } + // Populate the Universe vector and map. for (int i = 0; i < cells.size(); i++) { int32_t uid = cells[i]->universe_; @@ -766,7 +779,19 @@ extern "C" { int32_t cell_id(Cell* c) {return c->id_;} - void cell_set_id(Cell* c, int32_t id) {c->id_ = id;} + void + cell_set_id(Cell* c, int32_t id) + { + c->id_ = id; + + // Find the index of this cell and update the cell map. + for (int i = 0; i < cells.size(); i++) { + if (cells[i] == c) { + cell_map[id] = i; + break; + } + } + } int cell_type(Cell* c) {return c->type_;} @@ -789,10 +814,6 @@ extern "C" { int32_t cell_fill(Cell* c) {return c->fill_;} - int32_t cell_n_instances(Cell* c) {return c->n_instances_;} - - int cell_distribcell_index(Cell* c) {return c->distribcell_index_;} - int cell_material_size(Cell* c) {return c->material_.size();} //TODO: off-by-one @@ -807,8 +828,6 @@ extern "C" { double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT_[i];} - int32_t cell_offset(Cell* c, int map) {return c->offset_[map];} - void extend_cells_c(int32_t n) { cells.reserve(cells.size() + n); diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index ef85949b9..f2fb4dc26 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -99,7 +99,7 @@ contains select type(filt => filters(i_filter_mesh) % obj) type is (MeshFilter) - m = meshes(filt % mesh) + m = meshes(filt % mesh()) end select ! Set mesh widths @@ -154,27 +154,27 @@ contains ! Reset all bins to 1 do l = 1, size(t % filter) - call filter_matches(t % filter(l)) % bins % clear() - call filter_matches(t % filter(l)) % bins % push_back(1) + call filter_matches(t % filter(l)) % bins_clear() + call filter_matches(t % filter(l)) % bins_push_back(1) end do ! Set ijk as mesh indices ijk = (/ i, j, k /) ! Get bin number for mesh indices - filter_matches(i_filter_mesh) % bins % data(1) = & - m % get_bin_from_indices(ijk) + call filter_matches(i_filter_mesh) & + % bins_set_data(1, m % get_bin_from_indices(ijk)) ! Apply energy in filter if (energy_filters) then - filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1 + call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1) end if ! Calculate score index from bins score_index = 1 do l = 1, size(t % filter) score_index = score_index + (filter_matches(t % filter(l)) & - % bins % data(1) - 1) * t % stride(l) + % bins_data(1) - 1) * t % stride(l) end do ! Get flux @@ -198,30 +198,30 @@ contains ! Reset all bins to 1 do l = 1, size(t % filter) - call filter_matches(t % filter(l)) % bins % clear() - call filter_matches(t % filter(l)) % bins % push_back(1) + call filter_matches(t % filter(l)) % bins_clear() + call filter_matches(t % filter(l)) % bins_push_back(1) end do ! Set ijk as mesh indices ijk = (/ i, j, k /) ! Get bin number for mesh indices - filter_matches(i_filter_mesh) % bins % data(1) = & - m % get_bin_from_indices(ijk) + call filter_matches(i_filter_mesh) & + % bins_set_data(1, m % get_bin_from_indices(ijk)) if (energy_filters) then ! Apply energy in filter - filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1 + call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1) ! Set energy out bin - filter_matches(i_filter_eout) % bins % data(1) = ng - g + 1 + call filter_matches(i_filter_eout) % bins_set_data(1, ng - g + 1) end if ! Calculate score index from bins score_index = 1 do l = 1, size(t % filter) score_index = score_index + (filter_matches(t % filter(l)) & - % bins % data(1) - 1) * t % stride(l) + % bins_data(1) - 1) * t % stride(l) end do ! Get scattering @@ -244,23 +244,23 @@ contains ! Initialize and filter for energy do l = 1, size(t % filter) - call filter_matches(t % filter(l)) % bins % clear() - call filter_matches(t % filter(l)) % bins % push_back(1) + call filter_matches(t % filter(l)) % bins_clear() + call filter_matches(t % filter(l)) % bins_push_back(1) end do ! Set the bin for this mesh cell i_mesh = m % get_bin_from_indices([ i, j, k ]) - filter_matches(i_filter_mesh) % bins % data(1) = 12*(i_mesh - 1) + 1 + call filter_matches(i_filter_mesh) % bins_set_data(1, 12*(i_mesh - 1) + 1) ! Set the energy bin if needed if (energy_filters) then - filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1 + call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1) end if score_index = 0 do l = 1, size(t % filter) score_index = score_index + (filter_matches(t % filter(l)) & - % bins % data(1) - 1) * t % stride(l) + % bins_data(1) - 1) * t % stride(l) end do ! Left surface @@ -303,30 +303,30 @@ contains ! Reset all bins to 1 do l = 1, size(t % filter) - call filter_matches(t % filter(l)) % bins % clear() - call filter_matches(t % filter(l)) % bins % push_back(1) + call filter_matches(t % filter(l)) % bins_clear() + call filter_matches(t % filter(l)) % bins_push_back(1) end do ! Set ijk as mesh indices ijk = (/ i, j, k /) ! Get bin number for mesh indices - filter_matches(i_filter_mesh) % bins % data(1) = & - m % get_bin_from_indices(ijk) + call filter_matches(i_filter_mesh) & + % bins_set_data(1, m % get_bin_from_indices(ijk)) ! Apply energy in filter if (energy_filters) then - filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1 + call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1) end if ! Apply Legendre filter - filter_matches(i_filter_legendre) % bins % data(1) = 2 + call filter_matches(i_filter_legendre) % bins_set_data(1, 2) ! Calculate score index from bins score_index = 1 do l = 1, size(t % filter) score_index = score_index + (filter_matches(t % filter(l)) & - % bins % data(1) - 1) * t % stride(l) + % bins_data(1) - 1) * t % stride(l) end do ! Get p1 scatter rr and convert to p1 scatter xs diff --git a/src/constants.F90 b/src/constants.F90 index a3118f84d..250b1ef37 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -17,6 +17,9 @@ module constants ! HDF5 data format integer, parameter :: HDF5_VERSION(2) = [1, 0] + ! WMP data format + integer, parameter :: WMP_VERSION(2) = [1, 1] + ! Version numbers for binary files integer, parameter :: VERSION_STATEPOINT(2) = [17, 0] integer, parameter :: VERSION_PARTICLE_RESTART(2) = [2, 0] @@ -25,7 +28,6 @@ module constants integer, parameter :: VERSION_VOLUME(2) = [1, 0] integer, parameter :: VERSION_VOXEL(2) = [1, 0] integer, parameter :: VERSION_MGXS_LIBRARY(2) = [1, 0] - character(10), parameter :: VERSION_MULTIPOLE = "v1.0" ! ============================================================================ ! ADJUSTABLE PARAMETERS diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 959dabd12..d4e113977 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -54,7 +54,7 @@ void load_dagmc_geometry() c->fill_ = C_NONE; // no fill, single universe cells.push_back(c); - cell_map[c->id_] = c->id_; + cell_map[c->id_] = i; // Populate the Universe vector and dict auto it = universe_map.find(dagmc_univ_id); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 20f64e187..c0d9066a0 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -12,6 +12,8 @@ #include "openmc/material.h" #include "openmc/settings.h" #include "openmc/surface.h" +#include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_distribcell.h" namespace openmc { @@ -181,12 +183,15 @@ neighbor_lists() //============================================================================== void -prepare_distribcell(int32_t* filter_cell_list, int n) +prepare_distribcell() { - // Read the list of cells contained in distribcell filters from Fortran. + // Find all cells listed in a DistribcellFilter. std::unordered_set distribcells; - for (int i = 0; i < n; i++) { - distribcells.insert(filter_cell_list[i]); + for (auto* filt : tally_filters) { + if (filt->type() == "distribcell") { + auto* distrib_filt = static_cast(filt); + distribcells.insert(distrib_filt->cell_); + } } // Find all cells with distributed materials or temperatures. Make sure that @@ -393,29 +398,11 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, } } -//============================================================================== - -int -distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ) +std::string +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset) { - Universe& root = *universes[root_univ]; - std::string path_ {distribcell_path_inner(target_cell, map, target_offset, - root, 0)}; - return path_.size() + 1; -} - -//============================================================================== - -void -distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, - int32_t root_univ, char* path) -{ - Universe& root = *universes[root_univ]; - std::string path_ {distribcell_path_inner(target_cell, map, target_offset, - root, 0)}; - path_.copy(path, path_.size()); - path[path_.size()] = '\0'; + auto& root_univ = *universes[openmc_root_universe]; + return distribcell_path_inner(target_cell, map, target_offset, root_univ, 0); } //============================================================================== diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index b58b5415b..9af1fcbaa 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -59,20 +59,6 @@ module geometry_header integer(C_INT32_T) :: fill end function cell_fill_c - function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') & - result(n_instances) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT32_T) :: n_instances - end function cell_n_instances_c - - function cell_distribcell_index_c(cell_ptr) & - bind(C, name='cell_distribcell_index') result(distribcell_index) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT) :: distribcell_index - end function cell_distribcell_index_c - function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') & result(n) import C_PTR, C_INT @@ -103,14 +89,6 @@ module geometry_header real(C_DOUBLE) :: sqrtkT end function cell_sqrtkT_c - function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") & - result(offset) - import C_PTR, C_INT, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT), intent(in), value :: map - integer(C_INT32_T) :: offset - end function cell_offset_c - function lattice_pointer(lat_ind) bind(C) result(ptr) import C_PTR, C_INT32_T integer(C_INT32_T), intent(in), value :: lat_ind @@ -123,23 +101,6 @@ module geometry_header integer(C_INT32_T) :: id end function lattice_id_c - function lattice_are_valid_indices_c(lat_ptr, i_xyz) & - bind(C, name='lattice_are_valid_indices') result (is_valid) - import C_PTR, C_INT, C_BOOL - type(C_PTR), intent(in), value :: lat_ptr - integer(C_INT), intent(in) :: i_xyz(3) - logical(C_BOOL) :: is_valid - end function lattice_are_valid_indices_c - - function lattice_offset_c(lat_ptr, map, i_xyz) & - bind(C, name='lattice_offset') result(offset) - import C_PTR, C_INT, C_INT32_T - type(C_PTR), intent(in), value :: lat_ptr - integer(C_INT), intent(in), value :: map - integer(C_INT), intent(in) :: i_xyz(3) - integer(C_INT32_T) :: offset - end function lattice_offset_c - subroutine extend_cells_c(n) bind(C) import C_INT32_t integer(C_INT32_T), intent(in), value :: n @@ -154,8 +115,6 @@ module geometry_header type(C_PTR) :: ptr contains procedure :: id => lattice_id - procedure :: are_valid_indices => lattice_are_valid_indices - procedure :: offset => lattice_offset end type Lattice !=============================================================================== @@ -176,13 +135,10 @@ module geometry_header procedure :: type => cell_type procedure :: universe => cell_universe procedure :: fill => cell_fill - procedure :: n_instances => cell_n_instances - procedure :: distribcell_index => cell_distribcell_index procedure :: material_size => cell_material_size procedure :: material => cell_material procedure :: sqrtkT_size => cell_sqrtkT_size procedure :: sqrtkT => cell_sqrtkT - procedure :: offset => cell_offset end type Cell @@ -197,7 +153,6 @@ module geometry_header ! Dictionaries which map user IDs to indices in the global arrays type(DictIntInt) :: cell_dict - type(DictIntInt) :: universe_dict type(DictIntInt) :: lattice_dict contains @@ -208,21 +163,6 @@ contains id = lattice_id_c(this % ptr) end function lattice_id - function lattice_are_valid_indices(this, i_xyz) result (is_valid) - class(Lattice), intent(in) :: this - integer(C_INT), intent(in) :: i_xyz(3) - logical(C_BOOL) :: is_valid - is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz) - end function lattice_are_valid_indices - - function lattice_offset(this, map, i_xyz) result(offset) - class(Lattice), intent(in) :: this - integer(C_INT), intent(in) :: map - integer(C_INT), intent(in) :: i_xyz(3) - integer(C_INT32_T) :: offset - offset = lattice_offset_c(this % ptr, map, i_xyz) - end function lattice_offset - !=============================================================================== function cell_id(this) result(id) @@ -255,18 +195,6 @@ contains fill = cell_fill_c(this % ptr) end function cell_fill - function cell_n_instances(this) result(n_instances) - class(Cell), intent(in) :: this - integer(C_INT32_T) :: n_instances - n_instances = cell_n_instances_c(this % ptr) - end function cell_n_instances - - function cell_distribcell_index(this) result(distribcell_index) - class(Cell), intent(in) :: this - integer(C_INT) :: distribcell_index - distribcell_index = cell_distribcell_index_c(this % ptr) - end function cell_distribcell_index - function cell_material_size(this) result(n) class(Cell), intent(in) :: this integer(C_INT) :: n @@ -293,13 +221,6 @@ contains sqrtkT = cell_sqrtkT_c(this % ptr, i) end function cell_sqrtkT - function cell_offset(this, map) result(offset) - class(Cell), intent(in) :: this - integer(C_INT), intent(in) :: map - integer(C_INT32_T) :: offset - offset = cell_offset_c(this % ptr, map) - end function cell_offset - !=============================================================================== ! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table ! appears at in the model. Later, this list is used to determine the actual @@ -383,7 +304,6 @@ contains if (allocated(lattices)) deallocate(lattices) call cell_dict % clear() - call universe_dict % clear() call lattice_dict % clear() end subroutine free_memory_geometry diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8db0b1ee4..466bff9e7 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -22,6 +22,7 @@ module input_xml use mgxs_data, only: create_macro_xs, read_mgxs use mgxs_interface use nuclide_header + use multipole_header use output, only: title, header, print_plot use photon_header use plot_header @@ -59,12 +60,8 @@ module input_xml integer(C_INT32_T), intent(in), value :: univ_indx end subroutine count_cell_instances - subroutine prepare_distribcell_c(cell_list, n) & - bind(C, name="prepare_distribcell") - import C_INT32_T, C_INT - integer(C_INT), intent(in), value :: n - integer(C_INT32_T), intent(in) :: cell_list(n) - end subroutine prepare_distribcell_c + subroutine prepare_distribcell() bind(C) + end subroutine prepare_distribcell subroutine read_surfaces(node_ptr) bind(C) import C_PTR @@ -253,7 +250,7 @@ contains subroutine read_geometry_dagmc() - integer :: i, j + integer :: i integer :: univ_id integer :: n_cells_in_univ logical :: file_exists @@ -285,7 +282,6 @@ contains if (.not. cells_in_univ_dict % has(univ_id)) then n_universes = n_universes + 1 n_cells_in_univ = 1 - call universe_dict % set(univ_id, n_universes) call univ_ids % push_back(univ_id) else n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id) @@ -363,9 +359,6 @@ contains surfaces(i) % ptr = surface_pointer(i - 1); if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true. - - ! Add surface to dictionary - call surface_dict % set(surfaces(i) % id(), i) end do ! Check to make sure a boundary condition was applied to at least one @@ -455,7 +448,6 @@ contains if (.not. cells_in_univ_dict % has(univ_id)) then n_universes = n_universes + 1 n_cells_in_univ = 1 - call universe_dict % set(univ_id, n_universes - 1) call univ_ids % push_back(univ_id) else n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id) @@ -524,8 +516,6 @@ contains do i = 1, n_surfaces surfaces(i) % ptr = surface_pointer(i - 1); - ! Add surface to dictionary - call surface_dict % set(surfaces(i) % id(), i) end do end subroutine allocate_surfaces @@ -3014,7 +3004,9 @@ contains logical :: file_exists ! Does multipole library exist? character(7) :: readable ! Is multipole library readable? - character(MAX_FILE_LEN) :: filename ! Path to multipole xs library + character(MAX_FILE_LEN) :: filename ! Path to multipole xs library + integer(HID_T) :: file_id + integer(HID_T) :: group_id ! For the time being, and I know this is a bit hacky, we just assume ! that the file will be ZZZAAAmM.h5. @@ -3040,45 +3032,24 @@ contains end if ! Display message - call write_message("Loading Multipole XS table: " // filename, 6) + call write_message("Loading Windowed Multipole XS from " // filename, 6) + ! Open file and make sure version is sufficient + file_id = file_open(filename, 'r') + call check_wmp_version(file_id) + + ! Read nuclide data from HDF5 + group_id = open_group(file_id, nuc % name) allocate(nuc % multipole) - - ! Call the read routine - call nuc % multipole % from_hdf5(filename) + call nuc % multipole % from_hdf5(group_id) nuc % mp_present = .true. + ! Close the group and file. + call close_group(group_id) + call file_close(file_id) + end associate end subroutine read_multipole_data -!=============================================================================== -! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the -! offsets for distribcells -!=============================================================================== - - subroutine prepare_distribcell() - - integer :: i, j - type(SetInt) :: cell_list ! distribcells to track - integer(C_INT32_T), allocatable :: cell_list_c(:) - - ! Find all cells listed in a distribcell filter. - do i = 1, n_tallies - do j = 1, size(tallies(i) % obj % filter) - select type(filt => filters(tallies(i) % obj % filter(j)) % obj) - type is (DistribcellFilter) - call cell_list % add(filt % cell) - end select - end do - end do - - allocate(cell_list_c(cell_list % size())) - do i = 1, cell_list % size() - cell_list_c(i) = cell_list % get_item(i) - 1 - end do - call prepare_distribcell_c(cell_list_c, cell_list % size()) - - end subroutine prepare_distribcell - end module input_xml diff --git a/src/lattice.cpp b/src/lattice.cpp index 21f476850..14b679903 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -893,12 +893,6 @@ extern "C" { Lattice* lattice_pointer(int lat_ind) {return lattices[lat_ind];} int32_t lattice_id(Lattice *lat) {return lat->id_;} - - bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3]) - {return lat->are_valid_indices(i_xyz);} - - int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) - {return lat->offset(map, i_xyz);} } } // namespace openmc diff --git a/src/material.cpp b/src/material.cpp index ecfc7f55e..ecc6eed2b 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -23,7 +23,7 @@ std::unordered_map material_map; Material::Material(pugi::xml_node node) { if (check_for_node(node, "id")) { - id = std::stoi(get_node_value(node, "id")); + id_ = std::stoi(get_node_value(node, "id")); } else { fatal_error("Must specify id of material in materials XML file."); } @@ -52,7 +52,7 @@ read_materials(pugi::xml_node* node) // Populate the material map. for (int i = 0; i < materials.size(); i++) { - int32_t mid = materials[i]->id; + int32_t mid = materials[i]->id_; auto search = material_map.find(mid); if (search == material_map.end()) { material_map[mid] = i; @@ -78,7 +78,7 @@ openmc_material_get_volume(int32_t index, double* volume) return 0; } else { std::stringstream msg; - msg << "Volume for material with ID=" << m->id << " not set."; + msg << "Volume for material with ID=" << m->id_ << " not set."; set_errmsg(msg); return OPENMC_E_UNASSIGNED; } @@ -113,11 +113,11 @@ openmc_material_set_volume(int32_t index, double volume) extern "C" { Material* material_pointer(int32_t indx) {return materials[indx];} - int32_t material_id(Material* mat) {return mat->id;} + int32_t material_id(Material* mat) {return mat->id_;} void material_set_id(Material* mat, int32_t id, int32_t index) { - mat->id = id; + mat->id_ = id; //TODO: off-by-one material_map[id] = index - 1; } diff --git a/src/math.F90 b/src/math.F90 index 143e21e03..fe4156b9d 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -10,8 +10,6 @@ module math public :: t_percentile public :: calc_pn public :: calc_rn - public :: calc_zn - public :: calc_zn_rad public :: rotate_angle public :: faddeeva public :: w_derivative @@ -47,23 +45,6 @@ module math real(C_DOUBLE), intent(out) :: rn(2 * n + 1) end subroutine calc_rn - pure subroutine calc_zn(n, rho, phi, zn) bind(C, name='calc_zn_c') - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), value, intent(in) :: rho - real(C_DOUBLE), value, intent(in) :: phi - real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) - end subroutine calc_zn - - pure subroutine calc_zn_rad(n, rho, zn_rad) bind(C, name='calc_zn_rad_c') - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), value, intent(in) :: rho - real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1) - end subroutine calc_zn_rad - subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c') use ISO_C_BINDING implicit none diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 8d30ea84d..a8646f09d 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -510,7 +510,7 @@ void calc_rn_c(int n, const double uvw[3], double rn[]){ } -void calc_zn_c(int n, double rho, double phi, double zn[]) { +void calc_zn(int n, double rho, double phi, double zn[]) { // =========================================================================== // Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the // following recurrence relations so that only a single sin/cos have to be @@ -587,7 +587,7 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) { } -void calc_zn_rad_c(int n, double rho, double zn_rad[]) { +void calc_zn_rad(int n, double rho, double zn_rad[]) { // Calculate R_p0(rho) as Zn_p0(rho) // Set up the array of the coefficients diff --git a/src/mesh.cpp b/src/mesh.cpp index af088ffd0..a9b716ecd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -20,6 +20,7 @@ #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" #include "openmc/search.h" +#include "openmc/tallies/filter.h" #include "openmc/xml_interface.h" namespace openmc { @@ -899,10 +900,6 @@ void meshes_to_hdf5(hid_t group) //============================================================================== extern "C" { - // Declaration of Fortran procedures - void vector_int_push_back(void* ptr, int value); - void vector_real_push_back(void* ptr, double value); - int n_meshes() { return meshes.size(); } RegularMesh* mesh_ptr(int i) { return meshes.at(i).get(); } @@ -941,35 +938,6 @@ extern "C" { m->get_indices_from_bin(bin, ijk); } - void mesh_bins_crossed(RegularMesh* m, const Particle* p, void* match_bins, - void* match_weights) - { - // Get bins crossed - std::vector bins; - std::vector lengths; - m->bins_crossed(p, bins, lengths); - - // Call bindings for VectorInt and VectorReal on Fortran side - for (int i = 0; i < bins.size(); ++i) { - vector_int_push_back(match_bins, bins[i]); - vector_real_push_back(match_weights, lengths[i]); - } - } - - void mesh_surface_bins_crossed(RegularMesh* m, const Particle* p, - void* match_bins, void* match_weights) - { - // Get surface bins crossed - std::vector bins; - m->surface_bins_crossed(p, bins); - - // Call bindings for VectorInt and VectorReal - for (auto b : bins) { - vector_int_push_back(match_bins, b); - vector_real_push_back(match_weights, 1.0); - } - } - void free_memory_mesh() { meshes.clear(); diff --git a/src/multipole_header.F90 b/src/multipole_header.F90 index 64610e471..359bb3ac1 100644 --- a/src/multipole_header.F90 +++ b/src/multipole_header.F90 @@ -1,8 +1,9 @@ module multipole_header use constants - use error, only: fatal_error use hdf5_interface + use error, only: fatal_error + use string, only: to_str implicit none @@ -30,17 +31,15 @@ module multipole_header type MultipoleArray - !========================================================================= + character(20) :: name ! name of nuclide, e.g. U235 + ! Isotope Properties + logical :: fissionable ! Is this isotope fissionable? + complex(8), allocatable :: data(:,:) ! Poles and residues + real(8) :: sqrtAWR ! Square root of the atomic + ! weight ratio - logical :: fissionable ! Is this isotope fissionable? - complex(8), allocatable :: data(:,:) ! Poles and residues - real(8) :: sqrtAWR ! Square root of the atomic - ! weight ratio - - !========================================================================= ! Windows - integer :: fit_order ! Order of the fit. 1 linear, ! 2 quadratic, etc. real(8) :: E_min ! Start energy for the windows @@ -53,11 +52,9 @@ module multipole_header real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function. ! (reaction type, coeff index, ! window index) - integer, allocatable :: broaden_poly(:) ! if 1, broaden, if 0, don't. contains - procedure :: from_hdf5 => multipole_from_hdf5 end type MultipoleArray @@ -65,29 +62,22 @@ module multipole_header contains !=============================================================================== -! FROM_HDF5 loads multipole data from an HDF5 file. +! FROM_HDF5 loads multipole data from a HDF5 file. !=============================================================================== - subroutine multipole_from_hdf5(this, filename) + subroutine multipole_from_hdf5(this, group_id) class(MultipoleArray), intent(inout) :: this - character(len=*), intent(in) :: filename + integer(HID_T), intent(in) :: group_id - character(len=10) :: version integer :: n_poles, n_residues, n_windows integer(HSIZE_T) :: dims_1d(1), dims_2d(2), dims_3d(3) - integer(HID_T) :: file_id - integer(HID_T) :: group_id integer(HID_T) :: dset - ! Open file for reading and move into the /isotope group - file_id = file_open(filename, 'r', parallel=.true.) - group_id = open_group(file_id, "/nuclide") + ! Get name of nuclide from group + this % name = get_name(group_id) - ! Check the file version number. - call read_dataset(version, file_id, "version") - if (version /= VERSION_MULTIPOLE) call fatal_error("The current multipole& - & format version is " // trim(VERSION_MULTIPOLE) // " but the file "& - // trim(filename) // " uses version " // trim(version) // ".") + ! Get rid of leading '/' + this % name = trim(this % name(2:)) ! Read scalar values. call read_dataset(this % spacing, group_id, "spacing") @@ -121,8 +111,8 @@ contains dset = open_dataset(group_id, "broaden_poly") call get_shape(dset, dims_1d) if (dims_1d(1) /= n_windows) call fatal_error("broaden_poly array shape is& - ¬ consistent with the windows array shape in multipole library"& - // trim(filename) // ".") + ¬ consistent with the windows array shape in WMP library for"& + // trim(this % name) // ".") allocate(this % broaden_poly(n_windows)) call read_dataset(this % broaden_poly, dset) call close_dataset(dset) @@ -131,15 +121,38 @@ contains dset = open_dataset(group_id, "curvefit") call get_shape(dset, dims_3d) if (dims_3d(3) /= n_windows) call fatal_error("curvefit array shape is not& - &consistent with the windows array shape in multipole library"& - // trim(filename) // ".") + &consistent with the windows array shape in WMP library for"& + // trim(this % name) // ".") allocate(this % curvefit(dims_3d(1), dims_3d(2), dims_3d(3))) call read_dataset(this % curvefit, dset) call close_dataset(dset) this % fit_order = int(dims_3d(2), 4) - 1 - ! Close the group and file. - call close_group(group_id) - call file_close(file_id) end subroutine multipole_from_hdf5 + +!=============================================================================== +! CHECK_WMP_VERSION checks for the right version of WMP data within HDF5 +! files +!=============================================================================== + + subroutine check_wmp_version(file_id) + integer(HID_T), intent(in) :: file_id + + integer, allocatable :: version(:) + + if (attribute_exists(file_id, 'version')) then + call read_attribute(version, file_id, 'version') + if (version(1) /= WMP_VERSION(1)) then + call fatal_error("WMP data format uses version " // trim(to_str(& + version(1))) // "." // trim(to_str(version(2))) // " whereas & + &your installation of OpenMC expects version " // trim(to_str(& + WMP_VERSION(1))) // ".x data.") + end if + else + call fatal_error("WMP data does not indicate a version. Your & + &installation of OpenMC expects version " // trim(to_str(& + WMP_VERSION(1))) // ".x data.") + end if + end subroutine check_wmp_version + end module multipole_header diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 341ff6511..da66aea4a 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -53,7 +53,7 @@ module nuclide_header type :: Nuclide ! Nuclide meta-data - character(20) :: name ! name of nuclide, e.g. U235.71c + character(20) :: name ! name of nuclide, e.g. U235 integer :: Z ! atomic number integer :: A ! mass number integer :: metastable ! metastable state diff --git a/src/output.F90 b/src/output.F90 index ab37ceb50..ef3d3259c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -653,17 +653,13 @@ contains character(36) :: score_names(N_SCORE_TYPES) ! names of scoring function character(36) :: score_name ! names of scoring function ! to be applied at write-time - type(TallyFilterMatch), allocatable :: matches(:) + integer, allocatable :: filter_bins(:) character(MAX_WORD_LEN) :: temp_name ! Skip if there are no tallies if (n_tallies == 0) return - allocate(matches(n_filters)) - do i = 1, n_filters - allocate(matches(i) % bins) - allocate(matches(i) % weights) - end do + allocate(filter_bins(n_filters)) ! Initialize names for scores score_names(abs(SCORE_FLUX)) = "Flux" @@ -747,8 +743,7 @@ contains ! Initialize bins, filter level, and indentation do h = 1, size(t % filter) - call matches(t % filter(h)) % bins % clear() - call matches(t % filter(h)) % bins % push_back(0) + filter_bins(t % filter(h)) = 0 end do j = 1 indent = 0 @@ -759,18 +754,17 @@ contains if (size(t % filter) == 0) exit find_bin ! Increment bin combination - matches(t % filter(j)) % bins % data(1) = & - matches(t % filter(j)) % bins % data(1) + 1 + filter_bins(t % filter(j)) = filter_bins(t % filter(j)) + 1 ! ================================================================= ! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER - if (matches(t % filter(j)) % bins % data(1) > & + if (filter_bins(t % filter(j)) > & filters(t % filter(j)) % obj % n_bins) then ! If this is the first filter, then exit if (j == 1) exit print_bin - matches(t % filter(j)) % bins % data(1) = 0 + filter_bins(t % filter(j)) = 0 j = j - 1 indent = indent - 2 @@ -784,7 +778,7 @@ contains ! Print current filter information write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), & trim(filters(t % filter(j)) % obj % & - text_label(matches(t % filter(j)) % bins % data(1))) + text_label(filter_bins(t % filter(j)))) indent = indent + 2 j = j + 1 end if @@ -795,7 +789,7 @@ contains if (size(t % filter) > 0) then write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), & trim(filters(t % filter(j)) % obj % & - text_label(matches(t % filter(j)) % bins % data(1))) + text_label(filter_bins(t % filter(j)))) end if ! Determine scoring index for this bin combination -- note that unlike @@ -804,8 +798,8 @@ contains filter_index = 1 do h = 1, size(t % filter) - filter_index = filter_index + (max(matches(t % filter(h)) & - % bins % data(1),1) - 1) * t % stride(h) + filter_index = filter_index & + + (max(filter_bins(t % filter(h)) ,1) - 1) * t % stride(h) end do ! Write results for this filter bin combination @@ -859,11 +853,6 @@ contains close(UNIT=unit_tally) - do i = 1, n_filters - deallocate(matches(i) % bins) - deallocate(matches(i) % weights) - end do - end subroutine write_tallies !=============================================================================== diff --git a/src/simulation.F90 b/src/simulation.F90 index 01119b4e7..3eabbfecb 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -16,7 +16,7 @@ module simulation use settings use simulation_header use tally_header - use tally_filter_header, only: filter_matches, n_filters + use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer implicit none private @@ -44,8 +44,7 @@ contains ! Allocate array for matching filter bins allocate(filter_matches(n_filters)) do i = 1, n_filters - allocate(filter_matches(i) % bins) - allocate(filter_matches(i) % weights) + filter_matches(i) % ptr = filter_match_pointer(i - 1) end do !$omp end parallel @@ -65,10 +64,6 @@ contains deallocate(materials(i) % mat_nuclide_index) end do !$omp parallel - do i = 1, size(filter_matches) - deallocate(filter_matches(i) % bins) - deallocate(filter_matches(i) % weights) - end do deallocate(micro_xs, micro_photon_xs, filter_matches) !$omp end parallel diff --git a/src/simulation.cpp b/src/simulation.cpp index 60b22af9d..d132aab16 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -12,6 +12,7 @@ #include "openmc/source.h" #include "openmc/state_point.h" #include "openmc/timer.h" +#include "openmc/tallies/filter.h" #include "openmc/tallies/tally.h" #include @@ -53,7 +54,8 @@ extern "C" void write_tallies(); // over the batches, generations, and histories in a fixed source or k-eigenvalue // calculation. -int openmc_run() { +int openmc_run() +{ openmc_simulation_init(); int err = 0; @@ -76,6 +78,12 @@ int openmc_simulation_init() // Determine how much work each process should do calculate_work(); + // Allocate array for matching filter bins +#pragma omp parallel + { + filter_matches.resize(n_filters); + } + // Set up tally procedure pointers init_tally_routines(); @@ -134,6 +142,11 @@ int openmc_simulation_finalize() time_active.stop(); time_finalize.start(); +#pragma omp parallel + { + filter_matches.clear(); + } + // Deallocate Fortran variables, set tallies to inactive simulation_finalize_f(); diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 9d7e11461..a2ef9c056 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -2,7 +2,6 @@ module surface_header use, intrinsic :: ISO_C_BINDING - use dict_header, only: DictIntInt use hdf5_interface implicit none @@ -83,9 +82,6 @@ module surface_header type(Surface), allocatable, target :: surfaces(:) - ! Dictionary that maps user IDs to indices in 'surfaces' - type(DictIntInt) :: surface_dict - contains pure function surface_id(this) result(id) @@ -128,7 +124,6 @@ contains subroutine free_memory_surfaces() if (allocated(surfaces)) deallocate(surfaces) - call surface_dict % clear() call free_memory_surfaces_c() end subroutine free_memory_surfaces diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp new file mode 100644 index 000000000..44758f379 --- /dev/null +++ b/src/tallies/filter.cpp @@ -0,0 +1,170 @@ +#include "openmc/tallies/filter.h" + +#include + +#include "openmc/capi.h" +#include "openmc/constants.h" // for MAX_LINE_LEN; +#include "openmc/error.h" +#include "openmc/tallies/filter_azimuthal.h" +#include "openmc/tallies/filter_cell.h" +#include "openmc/tallies/filter_cellborn.h" +#include "openmc/tallies/filter_cellfrom.h" +#include "openmc/tallies/filter_distribcell.h" +#include "openmc/tallies/filter_energyfunc.h" +#include "openmc/tallies/filter_legendre.h" +#include "openmc/tallies/filter_material.h" +#include "openmc/tallies/filter_mesh.h" +#include "openmc/tallies/filter_meshsurface.h" +#include "openmc/tallies/filter_mu.h" +#include "openmc/tallies/filter_polar.h" +#include "openmc/tallies/filter_sph_harm.h" +#include "openmc/tallies/filter_sptl_legendre.h" +#include "openmc/tallies/filter_surface.h" +#include "openmc/tallies/filter_universe.h" +#include "openmc/tallies/filter_zernike.h" + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +std::vector filter_matches; +std::vector tally_filters; + +//============================================================================== +// Non-member functions +//============================================================================== + +void +free_memory_tally_c() +{ + #pragma omp parallel + { + filter_matches.clear(); + } + + for (Filter* filt : tally_filters) {delete filt;} + tally_filters.clear(); +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + FilterMatch* filter_match_pointer(int indx) + {return &filter_matches[indx];} + + void + filter_match_bins_push_back(FilterMatch* match, int val) + {match->bins_.push_back(val);} + + void + filter_match_weights_push_back(FilterMatch* match, double val) + {match->weights_.push_back(val);} + + void + filter_match_bins_clear(FilterMatch* match) + {match->bins_.clear();} + + void + filter_match_weights_clear(FilterMatch* match) + {match->weights_.clear();} + + int + filter_match_bins_size(FilterMatch* match) + {return match->bins_.size();} + + int + filter_match_bins_data(FilterMatch* match, int indx) + {return match->bins_[indx-1];} + + double + filter_match_weights_data(FilterMatch* match, int indx) + {return match->weights_[indx-1];} + + void + filter_match_bins_set_data(FilterMatch* match, int indx, int val) + {match->bins_[indx-1] = val;} + + Filter* + allocate_filter(const char* type) + { + std::string type_ {type}; + if (type_ == "azimuthal") { + tally_filters.push_back(new AzimuthalFilter()); + } else if (type_ == "cell") { + tally_filters.push_back(new CellFilter()); + } else if (type_ == "cellborn") { + tally_filters.push_back(new CellbornFilter()); + } else if (type_ == "cellfrom") { + tally_filters.push_back(new CellFromFilter()); + } else if (type_ == "distribcell") { + tally_filters.push_back(new DistribcellFilter()); + } else if (type_ == "energyfunction") { + tally_filters.push_back(new EnergyFunctionFilter()); + } else if (type_ == "legendre") { + tally_filters.push_back(new LegendreFilter()); + } else if (type_ == "material") { + tally_filters.push_back(new MaterialFilter()); + } else if (type_ == "mesh") { + tally_filters.push_back(new MeshFilter()); + } else if (type_ == "meshsurface") { + tally_filters.push_back(new MeshSurfaceFilter()); + } else if (type_ == "mu") { + tally_filters.push_back(new MuFilter()); + } else if (type_ == "polar") { + tally_filters.push_back(new PolarFilter()); + } else if (type_ == "surface") { + tally_filters.push_back(new SurfaceFilter()); + } else if (type_ == "spatiallegendre") { + tally_filters.push_back(new SpatialLegendreFilter()); + } else if (type_ == "sphericalharmonics") { + tally_filters.push_back(new SphericalHarmonicsFilter()); + } else if (type_ == "universe") { + tally_filters.push_back(new UniverseFilter()); + } else if (type_ == "zernike") { + tally_filters.push_back(new ZernikeFilter()); + } else if (type_ == "zernikeradial") { + tally_filters.push_back(new ZernikeRadialFilter()); + } else { + return nullptr; + } + return tally_filters.back(); + } + + void filter_from_xml(Filter* filt, pugi::xml_node* node) + {filt->from_xml(*node);} + + void + filter_get_all_bins(Filter* filt, Particle* p, int estimator, + FilterMatch* match) + { + filt->get_all_bins(p, estimator, *match); + } + + void filter_to_statepoint(Filter* filt, hid_t group) + {filt->to_statepoint(group);} + + void filter_text_label(Filter* filt, int bin, char* label) + { + std::string label_str = filt->text_label(bin); + int i = 0; + for (; i < label_str.size() && i < MAX_LINE_LEN; i++) + label[i] = label_str[i]; + label[i] = '\0'; + } + + void filter_initialize(Filter* filt) {filt->initialize();} + + int filter_n_bins(Filter* filt) {return filt->n_bins_;} + + int mesh_filter_get_mesh(MeshFilter* filt) {return filt->mesh();} + + int sphharm_filter_get_cosine(SphericalHarmonicsFilter* filt) + {return static_cast(filt->cosine_);} +} + +} // namespace openmc diff --git a/src/tallies/filter_azimuthal.cpp b/src/tallies/filter_azimuthal.cpp new file mode 100644 index 000000000..99a649cd9 --- /dev/null +++ b/src/tallies/filter_azimuthal.cpp @@ -0,0 +1,73 @@ +#include "openmc/tallies/filter_azimuthal.h" + +#include +#include + +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +AzimuthalFilter::from_xml(pugi::xml_node node) +{ + auto bins = get_node_array(node, "bins"); + + if (bins.size() > 1) { + bins_ = bins; + + } else { + // 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."); + + 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; + } + + n_bins_ = bins_.size() - 1; +} + +void +AzimuthalFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + double phi; + if (estimator == ESTIMATOR_TRACKLENGTH) { + phi = std::atan2(p->coord[0].uvw[1], p->coord[0].uvw[0]); + } else { + phi = std::atan2(p->last_uvw[1], p->last_uvw[0]); + } + + if (phi >= bins_.front() && phi <= bins_.back()) { + //TODO: off-by-one + auto bin = lower_bound_index(bins_.begin(), bins_.end(), phi) + 1; + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } +} + +void +AzimuthalFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); +} + +std::string +AzimuthalFilter::text_label(int bin) const +{ + std::stringstream out; + out << "Azimuthal Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); +} + +} // namespace openmc diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp new file mode 100644 index 000000000..a09c331bf --- /dev/null +++ b/src/tallies/filter_cell.cpp @@ -0,0 +1,91 @@ +#include "openmc/tallies/filter_cell.h" + +#include + +#include "openmc/capi.h" +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +CellFilter::from_xml(pugi::xml_node node) +{ + cells_ = get_node_array(node, "bins"); + n_bins_ = cells_.size(); +} + +void +CellFilter::initialize() +{ + // Convert cell IDs to indices of the global array. + for (auto& c : cells_) { + auto search = cell_map.find(c); + if (search != 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); + } + } + + // Populate the index->bin map. + for (int i = 0; i < cells_.size(); i++) { + map_[cells_[i]] = i; + } +} + +void +CellFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + for (int i = 0; i < p->n_coord; i++) { + auto search = map_.find(p->coord[i].cell); + if (search != map_.end()) { + //TODO: off-by-one + match.bins_.push_back(search->second + 1); + match.weights_.push_back(1.0); + } + } +} + +void +CellFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + std::vector cell_ids; + for (auto c : cells_) cell_ids.push_back(cells[c]->id_); + write_dataset(filter_group, "bins", cell_ids); +} + +std::string +CellFilter::text_label(int bin) const +{ + return "Cell " + std::to_string(cells[cells_[bin-1]]->id_); +} + +//============================================================================== +// C-API functions +//============================================================================== + +extern "C" int +openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n) +{ + int err = verify_filter(index); + if (err) return err; + + auto filt = filter_from_f(index); + if (filt->type() != "cell") { + set_errmsg("Tried to get cells from a non-cell filter."); + return OPENMC_E_INVALID_TYPE; + } + + auto cell_filt = static_cast(filt); + *cells = cell_filt->cells_.data(); + *n = cell_filt->cells_.size(); + return 0; +} + +} // namespace openmc diff --git a/src/tallies/filter_cellborn.cpp b/src/tallies/filter_cellborn.cpp new file mode 100644 index 000000000..e77c2e42f --- /dev/null +++ b/src/tallies/filter_cellborn.cpp @@ -0,0 +1,25 @@ +#include "openmc/tallies/filter_cellborn.h" + +#include "openmc/cell.h" + +namespace openmc { + +void +CellbornFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + auto search = map_.find(p->cell_born); + if (search != map_.end()) { + //TODO: off-by-one + match.bins_.push_back(search->second + 1); + match.weights_.push_back(1.0); + } +} + +std::string +CellbornFilter::text_label(int bin) const +{ + return "Birth Cell " + std::to_string(cells[cells_[bin-1]]->id_); +} + +} // namespace openmc diff --git a/src/tallies/filter_cellfrom.cpp b/src/tallies/filter_cellfrom.cpp new file mode 100644 index 000000000..4ce24e890 --- /dev/null +++ b/src/tallies/filter_cellfrom.cpp @@ -0,0 +1,27 @@ +#include "openmc/tallies/filter_cellfrom.h" + +#include "openmc/cell.h" + +namespace openmc { + +void +CellFromFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + for (int i = 0; i < p->last_n_coord; i++) { + auto search = map_.find(p->last_cell[i]); + if (search != map_.end()) { + //TODO: off-by-one + match.bins_.push_back(search->second + 1); + match.weights_.push_back(1.0); + } + } +} + +std::string +CellFromFilter::text_label(int bin) const +{ + return "Cell from " + std::to_string(cells[cells_[bin-1]]->id_); +} + +} // namespace openmc diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp new file mode 100644 index 000000000..f1e69b62a --- /dev/null +++ b/src/tallies/filter_distribcell.cpp @@ -0,0 +1,80 @@ +#include "openmc/tallies/filter_distribcell.h" + +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/geometry_aux.h" // For distribcell_path +#include "openmc/lattice.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +DistribcellFilter::from_xml(pugi::xml_node node) +{ + auto cells = get_node_array(node, "bins"); + if (cells.size() != 1) { + fatal_error("Only one cell can be specified per distribcell filter."); + } + cell_ = cells[0]; +} + +void +DistribcellFilter::initialize() +{ + // Convert the cell ID to an index of the global array. + auto search = cell_map.find(cell_); + if (search != cell_map.end()) { + cell_ = search->second; + n_bins_ = cells[cell_]->n_instances_; + } else { + std::stringstream err_msg; + err_msg << "Could not find cell " << cell_ + << " specified on tally filter."; + fatal_error(err_msg); + } +} + +void +DistribcellFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + int offset = 0; + auto distribcell_index = cells[cell_]->distribcell_index_; + for (int i = 0; i < p->n_coord; i++) { + auto& c {*cells[p->coord[i].cell]}; + if (c.type_ == FILL_UNIVERSE) { + offset += c.offset_[distribcell_index]; + } else if (c.type_ == FILL_LATTICE) { + auto& lat {*lattices[p->coord[i+1].lattice-1]}; + int i_xyz[3] {p->coord[i+1].lattice_x, + p->coord[i+1].lattice_y, + p->coord[i+1].lattice_z}; + if (lat.are_valid_indices(i_xyz)) { + offset += lat.offset(distribcell_index, i_xyz); + } + } + if (cell_ == p->coord[i].cell) { + //TODO: off-by-one + match.bins_.push_back(offset + 1); + match.weights_.push_back(1.0); + return; + } + } +} + +void +DistribcellFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", cells[cell_]->id_); +} + +std::string +DistribcellFilter::text_label(int bin) const +{ + auto map = cells[cell_]->distribcell_index_; + auto path = distribcell_path(cell_, map, bin-1); + return "Distributed Cell " + path; +} + +} // namespace openmc diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp new file mode 100644 index 000000000..b47e7b7b6 --- /dev/null +++ b/src/tallies/filter_energyfunc.cpp @@ -0,0 +1,69 @@ +#include "openmc/tallies/filter_energyfunc.h" + +#include // for setprecision +#include // for scientific +#include + +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/settings.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +EnergyFunctionFilter::from_xml(pugi::xml_node node) +{ + if (!settings::run_CE) + fatal_error("EnergyFunction filters are only supported for " + "continuous-energy transport calculations"); + + if (!check_for_node(node, "energy")) + fatal_error("Energy grid not specified for EnergyFunction filter."); + + energy_ = get_node_array(node, "energy"); + + if (!check_for_node(node, "y")) + fatal_error("y values not specified for EnergyFunction filter."); + + y_ = get_node_array(node, "y"); +} + +void +EnergyFunctionFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + if (p->last_E >= energy_.front() && p->last_E <= energy_.back()) { + // Search for the incoming energy bin. + auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E); + + // Compute the interpolation factor between the nearest bins. + double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]); + + // Interpolate on the lin-lin grid. + //TODO: off-by-one + match.bins_.push_back(1); + match.weights_.push_back((1-f) * y_[i] + f * y_[i+1]); + } +} + +void +EnergyFunctionFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "energy", energy_); + write_dataset(filter_group, "y", y_); +} + +std::string +EnergyFunctionFilter::text_label(int bin) const +{ + std::stringstream out; + out << std::scientific << std::setprecision(1) + << "Energy Function f" + << "([ " << energy_.front() << ", ..., " << energy_.back() << "]) = " + << "[" << y_.front() << ", ..., " << y_.back() << "]"; + return out.str(); +} + +} // namespace openmc diff --git a/src/tallies/filter_legendre.cpp b/src/tallies/filter_legendre.cpp new file mode 100644 index 000000000..3b59fa763 --- /dev/null +++ b/src/tallies/filter_legendre.cpp @@ -0,0 +1,93 @@ +#include "openmc/tallies/filter_legendre.h" + +#include "openmc/capi.h" +#include "openmc/error.h" +#include "openmc/math_functions.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +LegendreFilter::from_xml(pugi::xml_node node) +{ + order_ = std::stoi(get_node_value(node, "order")); + n_bins_ = order_ + 1; +} + +void +LegendreFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + double wgt[n_bins_]; + calc_pn_c(order_, p->mu, wgt); + for (int i = 0; i < n_bins_; i++) { + //TODO: off-by-one + match.bins_.push_back(i + 1); + match.weights_.push_back(wgt[i]); + } +} + +void +LegendreFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "order", order_); +} + +std::string +LegendreFilter::text_label(int bin) const +{ + return "Legendre expansion, P" + std::to_string(bin - 1); +} + +//============================================================================== +// C-API functions +//============================================================================== + +extern "C" int +openmc_legendre_filter_get_order(int32_t index, int* order) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) return err; + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Not a legendre filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Output the order. + *order = filt->order_; + return 0; +} + +extern "C" int +openmc_legendre_filter_set_order(int32_t index, int order) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) return err; + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Not a legendre filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Update the filter. + filt->order_ = order; + filt->n_bins_ = order + 1; + filter_update_n_bins(index); + return 0; +} + +} // namespace openmc diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp new file mode 100644 index 000000000..0334efae9 --- /dev/null +++ b/src/tallies/filter_material.cpp @@ -0,0 +1,123 @@ +#include "openmc/tallies/filter_material.h" + +#include + +#include "openmc/capi.h" +#include "openmc/error.h" +#include "openmc/material.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +MaterialFilter::from_xml(pugi::xml_node node) +{ + materials_ = get_node_array(node, "bins"); + n_bins_ = materials_.size(); +} + +void +MaterialFilter::initialize() +{ + // Convert material IDs to indices of the global array. + for (auto& m : materials_) { + auto search = material_map.find(m); + if (search != material_map.end()) { + m = search->second; + } else { + std::stringstream err_msg; + err_msg << "Could not find material " << m + << " specified on tally filter."; + fatal_error(err_msg); + } + } + + // Populate the index->bin map. + for (int i = 0; i < materials_.size(); i++) { + map_[materials_[i]] = i; + } +} + +void +MaterialFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + auto search = map_.find(p->material - 1); + if (search != map_.end()) { + //TODO: off-by-one + match.bins_.push_back(search->second + 1); + match.weights_.push_back(1.0); + } +} + +void +MaterialFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + std::vector material_ids; + for (auto c : materials_) material_ids.push_back(materials[c]->id_); + write_dataset(filter_group, "bins", material_ids); +} + +std::string +MaterialFilter::text_label(int bin) const +{ + return "Material " + std::to_string(materials[materials_[bin-1]]->id_); +} + +//============================================================================== +// C-API functions +//============================================================================== + +extern "C" int +openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) return err; + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Tried to get material filter bins on a non-material filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Output the bins. + *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) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) return err; + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Tried to set material filter bins on a non-material filter."); + return OPENMC_E_INVALID_TYPE; + } + + // 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; + filter_update_n_bins(index); + return 0; +} + +} // namespace openmc diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp new file mode 100644 index 000000000..5ff54bcf5 --- /dev/null +++ b/src/tallies/filter_mesh.cpp @@ -0,0 +1,137 @@ +#include "openmc/tallies/filter_mesh.h" + +#include + +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/mesh.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +MeshFilter::from_xml(pugi::xml_node node) +{ + auto bins_ = get_node_array(node, "bins"); + if (bins_.size() != 1) { + fatal_error("Only one mesh can be specified per " + type() + + " mesh filter."); + } + + auto id = bins_[0]; + auto search = mesh_map.find(id); + if (search != mesh_map.end()) { + set_mesh(search->second); + } else{ + std::stringstream err_msg; + err_msg << "Could not find cell " << id << " specified on tally filter."; + fatal_error(err_msg); + } +} + +void +MeshFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) +const +{ + if (estimator != ESTIMATOR_TRACKLENGTH) { + auto bin = meshes[mesh_]->get_bin(p->coord[0].xyz); + if (bin >= 0) { + //TODO: off-by-one + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } + } else { + meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_); + } +} + +void +MeshFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", meshes[mesh_]->id_); +} + +std::string +MeshFilter::text_label(int bin) const +{ + auto& mesh = *meshes[mesh_]; + int n_dim = mesh.n_dimension_; + + int ijk[n_dim]; + mesh.get_indices_from_bin(bin, ijk); + + std::stringstream out; + out << "Mesh Index (" << ijk[0]; + if (n_dim > 1) out << ", " << ijk[1]; + if (n_dim > 2) out << ", " << ijk[2]; + out << ")"; + + return out.str(); +} + +void +MeshFilter::set_mesh(int32_t mesh) +{ + mesh_ = mesh; + n_bins_ = 1; + for (auto dim : meshes[mesh_]->shape_) n_bins_ *= dim; +} + +//============================================================================== +// C-API functions +//============================================================================== + +extern "C" int +openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) return err; + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Tried to get mesh on a non-mesh filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Output the mesh. + *index_mesh = filt->mesh(); + return 0; +} + +extern "C" int +openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) return err; + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Tried to set mesh on a non-mesh filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Check the mesh index. + if (index_mesh < 0 || index_mesh >= meshes.size()) { + set_errmsg("Index in 'meshes' array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + // Update the filter. + filt->set_mesh(index_mesh); + filter_update_n_bins(index); + return 0; +} + +} // namespace openmc diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp new file mode 100644 index 000000000..8074b4a4c --- /dev/null +++ b/src/tallies/filter_meshsurface.cpp @@ -0,0 +1,94 @@ +#include "openmc/tallies/filter_meshsurface.h" + +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/mesh.h" + +namespace openmc { + +void +MeshSurfaceFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + meshes[mesh_]->surface_bins_crossed(p, match.bins_); + for (auto b : match.bins_) match.weights_.push_back(1.0); +} + +std::string +MeshSurfaceFilter::text_label(int bin) const +{ + auto& mesh = *meshes[mesh_]; + int n_dim = mesh.n_dimension_; + + // Get flattend mesh index and surface index. + int i_mesh = (bin - 1) / (4 * n_dim) + 1; + int i_surf = ((bin - 1) % (4 * n_dim)) + 1; + + // Get mesh index part of label. + std::string out = MeshFilter::text_label(i_mesh); + + // Get surface part of label. + switch (i_surf) { + case OUT_LEFT: + out += " Outgoing, x-min"; + break; + case IN_LEFT: + out += " Incoming, x-min"; + break; + case OUT_RIGHT: + out += " Outgoing, x-max"; + break; + case IN_RIGHT: + out += " Incoming, x-max"; + break; + case OUT_BACK: + out += " Outgoing, y-min"; + break; + case IN_BACK: + out += " Incoming, y-min"; + break; + case OUT_FRONT: + out += " Outgoing, y-max"; + break; + case IN_FRONT: + out += " Incoming, y-max"; + break; + case OUT_BOTTOM: + out += " Outgoing, z-min"; + break; + case IN_BOTTOM: + out += " Incoming, z-min"; + break; + case OUT_TOP: + out += " Outgoing, z-max"; + break; + case IN_TOP: + out += " Incoming, z-max"; + break; + } + + return out; +} + +void +MeshSurfaceFilter::set_mesh(int32_t mesh) +{ + mesh_ = mesh; + n_bins_ = 4 * meshes[mesh_]->n_dimension_; + for (auto dim : meshes[mesh_]->shape_) n_bins_ *= dim; +} + +//============================================================================== +// C-API functions +//============================================================================== + +extern"C" int +openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh) +{return openmc_mesh_filter_get_mesh(index, index_mesh);} + +extern"C" int +openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh) +{return openmc_mesh_filter_set_mesh(index, index_mesh);} + +} // namespace openmc diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp new file mode 100644 index 000000000..e57a8aa2e --- /dev/null +++ b/src/tallies/filter_mu.cpp @@ -0,0 +1,64 @@ +#include "openmc/tallies/filter_mu.h" + +#include + +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +MuFilter::from_xml(pugi::xml_node node) +{ + auto bins = get_node_array(node, "bins"); + + if (bins.size() > 1) { + bins_ = bins; + + } else { + // 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."); + + 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; + } + + n_bins_ = bins_.size() - 1; +} + +void +MuFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) +const +{ + if (p->mu >= bins_.front() && p->mu <= bins_.back()) { + //TODO: off-by-one + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu) + 1; + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } +} + +void +MuFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); +} + +std::string +MuFilter::text_label(int bin) const +{ + std::stringstream out; + out << "Change-in-Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); +} + +} // namespace openmc diff --git a/src/tallies/filter_polar.cpp b/src/tallies/filter_polar.cpp new file mode 100644 index 000000000..bfa87a341 --- /dev/null +++ b/src/tallies/filter_polar.cpp @@ -0,0 +1,72 @@ +#include "openmc/tallies/filter_polar.h" + +#include + +#include "openmc/constants.h" +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +PolarFilter::from_xml(pugi::xml_node node) +{ + auto bins = get_node_array(node, "bins"); + + if (bins.size() > 1) { + bins_ = bins; + + } else { + // 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."); + + 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; + } + + n_bins_ = bins_.size() - 1; +} + +void +PolarFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) +const +{ + double theta; + if (estimator == ESTIMATOR_TRACKLENGTH) { + theta = std::acos(p->coord[0].uvw[2]); + } else { + theta = std::acos(p->last_uvw[2]); + } + + if (theta >= bins_.front() && theta <= bins_.back()) { + //TODO: off-by-one + auto bin = lower_bound_index(bins_.begin(), bins_.end(), theta) + 1; + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } +} + +void +PolarFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); +} + +std::string +PolarFilter::text_label(int bin) const +{ + std::stringstream out; + out << "Polar Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); +} + +} // namespace openmc diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp new file mode 100644 index 000000000..764c56223 --- /dev/null +++ b/src/tallies/filter_sph_harm.cpp @@ -0,0 +1,182 @@ +#include "openmc/tallies/filter_sph_harm.h" + +#include // For pair + +#include "openmc/capi.h" +#include "openmc/error.h" +#include "openmc/math_functions.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +SphericalHarmonicsFilter::from_xml(pugi::xml_node node) +{ + order_ = std::stoi(get_node_value(node, "order")); + n_bins_ = (order_ + 1) * (order_ + 1); + + 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); + } + } +} + +void +SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + // Determine cosine term for scatter expansion if necessary + double wgt[order_ + 1]; + if (cosine_ == SphericalHarmonicsCosine::scatter) { + calc_pn_c(order_, p->mu, wgt); + } else { + for (int i = 0; i < order_ + 1; i++) wgt[i] = 1; + } + + // Find the Rn,m values + double rn[n_bins_]; + calc_rn_c(order_, p->last_uvw, rn); + + int j = 0; + for (int n = 0; n < order_ + 1; n++) { + // Calculate n-th order spherical harmonics for (u,v,w) + int num_nm = 2*n + 1; + + // Append the matching (bin,weight) for each moment + for (int i = 0; i < num_nm; i++) { + match.weights_.push_back(wgt[n] * rn[j]); + //TODO: off-by-one + match.bins_.push_back(++j); + } + } +} + +void +SphericalHarmonicsFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "order", order_); + if (cosine_ == SphericalHarmonicsCosine::scatter) { + write_dataset(filter_group, "cosine", "scatter"); + } else { + write_dataset(filter_group, "cosine", "particle"); + } +} + +std::string +SphericalHarmonicsFilter::text_label(int bin) const +{ + std::stringstream out; + for (int n = 0; n < order_ + 1; n++) { + if (bin <= (n + 1) * (n + 1)) { + int m = (bin - n*n - 1) - n; + out << "Spherical harmonic expansion, Y" << n << "," << m; + return out.str(); + } + } +} + +//============================================================================== +// C-API functions +//============================================================================== + +std::pair +check_sphharm_filter(int32_t index) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) { + return {err, nullptr}; + } + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Not a spherical harmonics filter."); + err = OPENMC_E_INVALID_TYPE; + } + return {err, filt}; +} + +extern "C" int +openmc_sphharm_filter_get_order(int32_t index, int* order) +{ + // Check the filter. + auto check_result = check_sphharm_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Output the order. + *order = filt->order_; + return 0; +} + +extern "C" int +openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]) +{ + // Check the filter. + auto check_result = check_sphharm_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Output the cosine. + if (filt->cosine_ == SphericalHarmonicsCosine::scatter) { + strcpy(cosine, "scatter"); + } else { + strcpy(cosine, "particle"); + } + return 0; +} + +extern "C" int +openmc_sphharm_filter_set_order(int32_t index, int order) +{ + // Check the filter. + auto check_result = check_sphharm_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Update the filter. + filt->order_ = order; + filt->n_bins_ = (order + 1) * (order + 1); + filter_update_n_bins(index); + return 0; +} + +extern "C" int +openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]) +{ + // Check the filter. + auto check_result = check_sphharm_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + 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."); + return OPENMC_E_INVALID_ARGUMENT; + } + return 0; +} + +} // namespace openmc diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp new file mode 100644 index 000000000..3f93be1cb --- /dev/null +++ b/src/tallies/filter_sptl_legendre.cpp @@ -0,0 +1,184 @@ +#include "openmc/tallies/filter_sptl_legendre.h" + +#include // For pair + +#include "openmc/capi.h" +#include "openmc/error.h" +#include "openmc/math_functions.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +SpatialLegendreFilter::from_xml(pugi::xml_node node) +{ + 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"); + } + + min_ = std::stod(get_node_value(node, "min")); + max_ = std::stod(get_node_value(node, "max")); + + n_bins_ = order_ + 1; +} + +void +SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + // Get the coordinate along the axis of interest. + double x; + if (axis_ == LegendreAxis::x) { + x = p->coord[0].xyz[0]; + } else if (axis_ == LegendreAxis::y) { + x = p->coord[0].xyz[1]; + } else { + x = p->coord[0].xyz[2]; + } + + if (x >= min_ && x <= max_) { + // Compute the normalized coordinate value. + double x_norm = 2.0*(x - min_) / (max_ - min_) - 1.0; + + // Compute and return the Legendre weights. + double wgt[order_ + 1]; + calc_pn_c(order_, x_norm, wgt); + for (int i = 0; i < order_ + 1; i++) { + //TODO: off-by-one + match.bins_.push_back(i + 1); + match.weights_.push_back(wgt[i]); + } + } +} + +void +SpatialLegendreFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "order", order_); + if (axis_ == LegendreAxis::x) { + write_dataset(filter_group, "axis", "x"); + } else if (axis_ == LegendreAxis::y) { + write_dataset(filter_group, "axis", "y"); + } else { + write_dataset(filter_group, "axis", "z"); + } + write_dataset(filter_group, "min", min_); + write_dataset(filter_group, "max", max_); +} + +std::string +SpatialLegendreFilter::text_label(int bin) const +{ + std::stringstream out; + out << "Legendre expansion, "; + if (axis_ == LegendreAxis::x) { + out << "x"; + } else if (axis_ == LegendreAxis::y) { + out << "y"; + } else { + out << "z"; + } + out << " axis, P" << std::to_string(bin - 1); + return out.str(); +} + +//============================================================================== +// C-API functions +//============================================================================== + +std::pair +check_sptl_legendre_filter(int32_t index) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) { + return {err, nullptr}; + } + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Not a spatial Legendre filter."); + err = OPENMC_E_INVALID_TYPE; + } + return {err, filt}; +} + +extern "C" int +openmc_spatial_legendre_filter_get_order(int32_t index, int* order) +{ + // Check the filter. + auto check_result = check_sptl_legendre_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Output the order. + *order = filt->order_; + return 0; +} + +extern "C" int +openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, + double* min, double* max) +{ + // Check the filter. + auto check_result = check_sptl_legendre_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Output the params. + *axis = static_cast(filt->axis_); + *min = filt->min_; + *max = filt->max_; + return 0; +} + +extern "C" int +openmc_spatial_legendre_filter_set_order(int32_t index, int order) +{ + // Check the filter. + auto check_result = check_sptl_legendre_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Update the filter. + filt->order_ = order; + filt->n_bins_ = order + 1; + filter_update_n_bins(index); + return 0; +} + +extern "C" int +openmc_spatial_legendre_filter_set_params(int32_t index, const int* axis, + const double* min, const double* max) +{ + // Check the filter. + auto check_result = check_sptl_legendre_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Update the filter. + if (axis) filt->axis_ = static_cast(*axis); + if (min) filt->min_ = *min; + if (max) filt->max_ = *max; + return 0; +} + +} // namespace openmc diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp new file mode 100644 index 000000000..7fd2c6e24 --- /dev/null +++ b/src/tallies/filter_surface.cpp @@ -0,0 +1,71 @@ +#include "openmc/tallies/filter_surface.h" + +#include + +#include "openmc/error.h" +#include "openmc/surface.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +SurfaceFilter::from_xml(pugi::xml_node node) +{ + surfaces_ = get_node_array(node, "bins"); + n_bins_ = surfaces_.size(); +} + +void +SurfaceFilter::initialize() +{ + // Convert surface IDs to indices of the global array. + for (auto& s : surfaces_) { + auto search = surface_map.find(s); + if (search != surface_map.end()) { + s = search->second; + } else { + std::stringstream err_msg; + err_msg << "Could not find surface " << s + << " specified on tally filter."; + fatal_error(err_msg); + } + } + + // Populate the index->bin map. + for (int i = 0; i < surfaces_.size(); i++) { + map_[surfaces_[i]] = i; + } +} + +void +SurfaceFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + auto search = map_.find(std::abs(p->surface)-1); + if (search != map_.end()) { + //TODO: off-by-one + match.bins_.push_back(search->second + 1); + if (p->surface < 0) { + match.weights_.push_back(-1.0); + } else { + match.weights_.push_back(1.0); + } + } +} + +void +SurfaceFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + std::vector surface_ids; + for (auto c : surfaces_) surface_ids.push_back(surfaces[c]->id_); + write_dataset(filter_group, "bins", surface_ids); +} + +std::string +SurfaceFilter::text_label(int bin) const +{ + return "Surface " + std::to_string(surfaces[surfaces_[bin-1]]->id_); +} + +} // namespace openmc diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp new file mode 100644 index 000000000..de0d33d29 --- /dev/null +++ b/src/tallies/filter_universe.cpp @@ -0,0 +1,69 @@ +#include "openmc/tallies/filter_universe.h" + +#include + +#include "openmc/cell.h" +#include "openmc/error.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +void +UniverseFilter::from_xml(pugi::xml_node node) +{ + universes_ = get_node_array(node, "bins"); + n_bins_ = universes_.size(); +} + +void +UniverseFilter::initialize() +{ + // Convert universe IDs to indices of the global array. + for (auto& u : universes_) { + auto search = universe_map.find(u); + if (search != universe_map.end()) { + u = search->second; + } else { + std::stringstream err_msg; + err_msg << "Could not find universe " << u + << " specified on tally filter."; + fatal_error(err_msg); + } + } + + // Populate the index->bin map. + for (int i = 0; i < universes_.size(); i++) { + map_[universes_[i]] = i; + } +} + +void +UniverseFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + for (int i = 0; i < p->n_coord; i++) { + auto search = map_.find(p->coord[i].universe); + if (search != map_.end()) { + //TODO: off-by-one + match.bins_.push_back(search->second + 1); + match.weights_.push_back(1.0); + } + } +} + +void +UniverseFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + std::vector universe_ids; + for (auto u : universes_) universe_ids.push_back(universes[u]->id_); + write_dataset(filter_group, "bins", universe_ids); +} + +std::string +UniverseFilter::text_label(int bin) const +{ + return "Universe " + std::to_string(universes[universes_[bin-1]]->id_); +} + +} // namespace openmc diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp new file mode 100644 index 000000000..198bd8cb7 --- /dev/null +++ b/src/tallies/filter_zernike.cpp @@ -0,0 +1,207 @@ +#include "openmc/tallies/filter_zernike.h" + +#include +#include +#include // For pair + +#include "openmc/capi.h" +#include "openmc/error.h" +#include "openmc/math_functions.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +//============================================================================== +// ZernikeFilter implementation +//============================================================================== + +void +ZernikeFilter::from_xml(pugi::xml_node node) +{ + set_order(std::stoi(get_node_value(node, "order"))); + x_ = std::stod(get_node_value(node, "x")); + y_ = std::stod(get_node_value(node, "y")); + r_ = std::stod(get_node_value(node, "r")); +} + +void +ZernikeFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + // Determine the normalized (r,theta) coordinates. + double x = p->coord[0].xyz[0] - x_; + double y = p->coord[0].xyz[1] - y_; + double r = std::sqrt(x*x + y*y) / r_; + double theta = std::atan2(y, x); + + if (r <= 1.0) { + // Compute and return the Zernike weights. + double zn[n_bins_]; + calc_zn(order_, r, theta, zn); + for (int i = 0; i < n_bins_; i++) { + //TODO: off-by-one + match.bins_.push_back(i+1); + match.weights_.push_back(zn[i]); + } + } +} + +void +ZernikeFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "order", order_); + write_dataset(filter_group, "x", x_); + write_dataset(filter_group, "y", y_); + write_dataset(filter_group, "r", r_); +} + +std::string +ZernikeFilter::text_label(int bin) const +{ + std::stringstream out; + for (int n = 0; n < order_+1; n++) { + int last = (n + 1) * (n + 2) / 2; + if (bin <= last) { + int first = last - n; + int m = -n + (bin - first) * 2; + out << "Zernike expansion, Z" << n << "," << m; + return out.str(); + } + } +} + +void +ZernikeFilter::set_order(int order) +{ + order_ = order; + n_bins_ = ((order+1) * (order+2)) / 2; +} + +//============================================================================== +// ZernikeRadialFilter implementation +//============================================================================== + +void +ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + // Determine the normalized radius coordinate. + double x = p->coord[0].xyz[0] - x_; + double y = p->coord[0].xyz[1] - y_; + double r = std::sqrt(x*x + y*y) / r_; + + if (r <= 1.0) { + // Compute and return the Zernike weights. + double zn[n_bins_]; + calc_zn_rad(order_, r, zn); + for (int i = 0; i < n_bins_; i++) { + //TODO: off-by-one + match.bins_.push_back(i+1); + match.weights_.push_back(zn[i]); + } + } +} + +std::string +ZernikeRadialFilter::text_label(int bin) const +{ + return "Zernike expansion, Z" + std::to_string(2*(bin-1)) + ",0"; +} + +void +ZernikeRadialFilter::set_order(int order) +{ + order_ = order; + n_bins_ = order / 2 + 1; +} + +//============================================================================== +// C-API functions +//============================================================================== + +std::pair +check_zernike_filter(int32_t index) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) { + return {err, nullptr}; + } + + // Get a pointer to the filter and downcast. + auto* filt_base = filter_from_f(index); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Not a Zernike filter."); + err = OPENMC_E_INVALID_TYPE; + } + return {err, filt}; +} + +extern "C" int +openmc_zernike_filter_get_order(int32_t index, int* order) +{ + // Check the filter. + auto check_result = check_zernike_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Output the order. + *order = filt->order(); + return 0; +} + +extern "C" int +openmc_zernike_filter_get_params(int32_t index, double* x, double* y, + double* r) +{ + // Check the filter. + auto check_result = check_zernike_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Output the params. + *x = filt->x_; + *y = filt->y_; + *r = filt->r_; + return 0; +} + +extern "C" int +openmc_zernike_filter_set_order(int32_t index, int order) +{ + // Check the filter. + auto check_result = check_zernike_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Update the filter. + filt->set_order(order); + filter_update_n_bins(index); + return 0; +} + +extern "C" int +openmc_zernike_filter_set_params(int32_t index, const double* x, + const double* y, const double* r) +{ + // Check the filter. + auto check_result = check_zernike_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) return err; + + // Update the filter. + if (x) filt->x_ = *x; + if (y) filt->y_ = *y; + if (r) filt->r_ = *r; + return 0; +} + +} // namespace openmc diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 19d99a014..ba152353e 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -2079,15 +2079,15 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) if (.not. filter_matches(i_filt) % bins_present) then - call filter_matches(i_filt) % bins % clear() - call filter_matches(i_filt) % weights % clear() + call filter_matches(i_filt) % bins_clear() + call filter_matches(i_filt) % weights_clear() call filters(i_filt) % obj % get_all_bins(p, t % estimator, & filter_matches(i_filt)) filter_matches(i_filt) % bins_present = .true. end if ! If there are no valid bins for this filter, then there is nothing to ! score and we can move on to the next tally. - if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP + if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP ! Set the index of the bin used in the first filter combination filter_matches(i_filt) % i_bin = 1 @@ -2106,10 +2106,10 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) i_bin = filter_matches(i_filt) % i_bin - filter_index = filter_index + (filter_matches(i_filt) % bins % & - data(i_bin) - 1) * t % stride(j) - filter_weight = filter_weight * filter_matches(i_filt) % weights % & - data(i_bin) + filter_index = filter_index + (filter_matches(i_filt) & + % bins_data(i_bin) - 1) * t % stride(j) + filter_weight = filter_weight * filter_matches(i_filt) & + % weights_data(i_bin) end do ! ====================================================================== @@ -2166,7 +2166,7 @@ contains do j = size(t % filter), 1, -1 i_filt = t % filter(j) if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % & - bins % size()) then + bins_size()) then filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1 finished = .false. exit @@ -2225,15 +2225,15 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) if (.not. filter_matches(i_filt) % bins_present) then - call filter_matches(i_filt) % bins % clear() - call filter_matches(i_filt) % weights % clear() + call filter_matches(i_filt) % bins_clear() + call filter_matches(i_filt) % weights_clear() call filters(i_filt) % obj % get_all_bins(p, t % estimator, & filter_matches(i_filt)) filter_matches(i_filt) % bins_present = .true. end if ! If there are no valid bins for this filter, then there is nothing to ! score and we can move on to the next tally. - if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP + if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP ! Set the index of the bin used in the first filter combination filter_matches(i_filt) % i_bin = 1 @@ -2252,10 +2252,10 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) i_bin = filter_matches(i_filt) % i_bin - filter_index = filter_index + (filter_matches(i_filt) % bins % & - data(i_bin) - 1) * t % stride(j) - filter_weight = filter_weight * filter_matches(i_filt) % weights % & - data(i_bin) + filter_index = filter_index + (filter_matches(i_filt) & + % bins_data(i_bin) - 1) * t % stride(j) + filter_weight = filter_weight * filter_matches(i_filt) & + % weights_data(i_bin) end do ! ====================================================================== @@ -2295,7 +2295,7 @@ contains do j = size(t % filter), 1, -1 i_filt = t % filter(j) if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % & - bins % size()) then + bins_size()) then filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1 finished = .false. exit @@ -2360,7 +2360,7 @@ contains ! save original outgoing energy bin and score index i = t % filter(t % find_filter(FILTER_ENERGYOUT)) i_bin = filter_matches(i) % i_bin - bin_energyout = filter_matches(i) % bins % data(i_bin) + bin_energyout = filter_matches(i) % bins_data(i_bin) ! declare the energyout filter type select type(eo_filt => filters(i) % obj) @@ -2400,7 +2400,7 @@ contains g_out = size(eo_filt % bins) - g_out ! change outgoing energy bin - filter_matches(i) % bins % data(i_bin) = g_out + call filter_matches(i) % bins_set_data(i_bin, g_out) else @@ -2416,8 +2416,8 @@ contains if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle ! change outgoing energy bin - filter_matches(i) % bins % data(i_bin) = & - binary_search(eo_filt % bins, n, E_out) + call filter_matches(i) % bins_set_data(i_bin, & + binary_search(eo_filt % bins, n, E_out)) end if @@ -2428,8 +2428,8 @@ contains ! determine scoring index and weight for this filter combination i_filter = 1 do l = 1, size(t % filter) - i_filter = i_filter + (filter_matches(t % filter(l)) % bins % & - data(filter_matches(t % filter(l)) % i_bin) - 1) * & + i_filter = i_filter + (filter_matches(t % filter(l)) & + % bins_data(filter_matches(t % filter(l)) % i_bin) - 1) * & t % stride(l) end do @@ -2470,10 +2470,10 @@ contains do l = 1, size(t % filter) f = t % filter(l) b = filter_matches(f) % i_bin - i_filter = i_filter + (filter_matches(f) % bins % & - data(b) - 1) * t % stride(l) - filter_weight = filter_weight * filter_matches(f) % & - weights % data(b) + i_filter = i_filter + (filter_matches(f) & + % bins_data(b) - 1) * t % stride(l) + filter_weight = filter_weight * filter_matches(f) & + % weights_data(b) end do call score_fission_delayed_dg(t, d_bin, & @@ -2493,10 +2493,10 @@ contains do l = 1, size(t % filter) f = t % filter(l) b = filter_matches(f) % i_bin - i_filter = i_filter + (filter_matches(f) % bins % data(b) - 1) & + i_filter = i_filter + (filter_matches(f) % bins_data(b) - 1) & * t % stride(l) - filter_weight = filter_weight * filter_matches(f) % weights % & - data(b) + filter_weight = filter_weight * filter_matches(f) & + % weights_data(b) end do ! Add score to tally @@ -2509,7 +2509,7 @@ contains end select ! reset outgoing energy bin and score index - filter_matches(i) % bins % data(i_bin) = bin_energyout + call filter_matches(i) % bins_set_data(i_bin, bin_energyout) end subroutine score_fission_eout @@ -2534,14 +2534,14 @@ contains ! save original delayed group bin i_filt = t % filter(t % find_filter(FILTER_DELAYEDGROUP)) i_bin = filter_matches(i_filt) % i_bin - bin_original = filter_matches(i_filt) % bins % data(i_bin) - filter_matches(i_filt) % bins % data(i_bin) = d_bin + bin_original = filter_matches(i_filt) % bins_data(i_bin) + call filter_matches(i_filt) % bins_set_data(i_bin, d_bin) ! determine scoring index and weight on the modified matching bins filter_index = 1 do i = 1, size(t % filter) - filter_index = filter_index + (filter_matches(t % filter(i)) % bins % & - data(filter_matches(t % filter(i)) % i_bin) - 1) * t % stride(i) + filter_index = filter_index + (filter_matches(t % filter(i)) % & + bins_data(filter_matches(t % filter(i)) % i_bin) - 1) * t % stride(i) end do !$omp atomic @@ -2549,7 +2549,7 @@ contains t % results(RESULT_VALUE, score_index, filter_index) + score ! reset original delayed group bin - filter_matches(i_filt) % bins % data(i_bin) = bin_original + call filter_matches(i_filt) % bins_set_data(i_bin, bin_original) end subroutine score_fission_delayed_dg @@ -2595,15 +2595,15 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) if (.not. filter_matches(i_filt) % bins_present) then - call filter_matches(i_filt) % bins % clear() - call filter_matches(i_filt) % weights % clear() + call filter_matches(i_filt) % bins_clear() + call filter_matches(i_filt) % weights_clear() call filters(i_filt) % obj % get_all_bins(p, t % estimator, & filter_matches(i_filt)) filter_matches(i_filt) % bins_present = .true. end if ! If there are no valid bins for this filter, then there is nothing to ! score and we can move on to the next tally. - if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP + if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP ! Set the index of the bin used in the first filter combination filter_matches(i_filt) % i_bin = 1 @@ -2622,10 +2622,10 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) i_bin = filter_matches(i_filt) % i_bin - filter_index = filter_index + (filter_matches(i_filt) % bins % & - data(i_bin) - 1) * t % stride(j) - filter_weight = filter_weight * filter_matches(i_filt) % weights % & - data(i_bin) + filter_index = filter_index + (filter_matches(i_filt) & + % bins_data(i_bin) - 1) * t % stride(j) + filter_weight = filter_weight * filter_matches(i_filt) & + % weights_data(i_bin) end do ! ====================================================================== @@ -2674,7 +2674,7 @@ contains do j = size(t % filter), 1, -1 i_filt = t % filter(j) if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % & - bins % size()) then + bins_size()) then filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1 finished = .false. exit @@ -2752,15 +2752,15 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) if (.not. filter_matches(i_filt) % bins_present) then - call filter_matches(i_filt) % bins % clear() - call filter_matches(i_filt) % weights % clear() + call filter_matches(i_filt) % bins_clear() + call filter_matches(i_filt) % weights_clear() call filters(i_filt) % obj % get_all_bins(p, t % estimator, & filter_matches(i_filt)) filter_matches(i_filt) % bins_present = .true. end if ! If there are no valid bins for this filter, then there is nothing to ! score and we can move on to the next tally. - if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP + if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP ! Set the index of the bin used in the first filter combination filter_matches(i_filt) % i_bin = 1 @@ -2779,10 +2779,10 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) i_bin = filter_matches(i_filt) % i_bin - filter_index = filter_index + (filter_matches(i_filt) % bins % & - data(i_bin) - 1) * t % stride(j) - filter_weight = filter_weight * filter_matches(i_filt) % weights % & - data(i_bin) + filter_index = filter_index + (filter_matches(i_filt) & + % bins_data(i_bin) - 1) * t % stride(j) + filter_weight = filter_weight * filter_matches(i_filt) & + % weights_data(i_bin) end do ! ====================================================================== @@ -2831,7 +2831,7 @@ contains do j = size(t % filter), 1, -1 i_filt = t % filter(j) if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % & - bins % size()) then + bins_size()) then filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1 finished = .false. exit @@ -2898,15 +2898,15 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) if (.not. filter_matches(i_filt) % bins_present) then - call filter_matches(i_filt) % bins % clear() - call filter_matches(i_filt) % weights % clear() + call filter_matches(i_filt) % bins_clear() + call filter_matches(i_filt) % weights_clear() call filters(i_filt) % obj % get_all_bins(p, t % estimator, & filter_matches(i_filt)) filter_matches(i_filt) % bins_present = .true. end if ! If there are no valid bins for this filter, then there is nothing to ! score and we can move on to the next tally. - if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP + if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP ! Set the index of the bin used in the first filter combination filter_matches(i_filt) % i_bin = 1 @@ -2925,10 +2925,10 @@ contains do j = 1, size(t % filter) i_filt = t % filter(j) i_bin = filter_matches(i_filt) % i_bin - filter_index = filter_index + (filter_matches(i_filt) % bins % & - data(i_bin) - 1) * t % stride(j) - filter_weight = filter_weight * filter_matches(i_filt) % weights % & - data(i_bin) + filter_index = filter_index + (filter_matches(i_filt) & + % bins_data(i_bin) - 1) * t % stride(j) + filter_weight = filter_weight * filter_matches(i_filt) & + % weights_data(i_bin) end do ! Determine score @@ -2961,7 +2961,7 @@ contains do j = size(t % filter), 1, -1 i_filt = t % filter(j) if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % & - bins % size()) then + bins_size()) then filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1 finished = .false. exit diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index 81b195816..98b7bf9a3 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -7,26 +7,15 @@ module tally_filter use tally_filter_header ! Inherit other filters - use tally_filter_azimuthal - use tally_filter_cell - use tally_filter_cellborn - use tally_filter_cellfrom + use tally_filter_cpp use tally_filter_delayedgroup use tally_filter_distribcell use tally_filter_energy - use tally_filter_energyfunc use tally_filter_legendre - use tally_filter_material use tally_filter_mesh use tally_filter_meshsurface - use tally_filter_mu use tally_filter_particle - use tally_filter_polar use tally_filter_sph_harm - use tally_filter_sptl_legendre - use tally_filter_surface - use tally_filter_universe - use tally_filter_zernike implicit none @@ -114,6 +103,14 @@ contains character(:), allocatable :: type_ + interface + function allocate_filter(type) result(ptr) bind(C) + import C_CHAR, C_PTR + character(kind=C_CHAR), intent(in) :: type(*) + type(C_PTR) :: ptr + end function allocate_filter + end interface + ! Convert C string to Fortran string type_ = to_f_string(type) @@ -131,7 +128,7 @@ contains case ('cellborn') allocate(CellbornFilter :: filters(index) % obj) case ('cellfrom') - allocate(CellfromFilter :: filters(index) % obj) + allocate(CellFromFilter :: filters(index) % obj) case ('delayedgroup') allocate(DelayedGroupFilter :: filters(index) % obj) case ('distribcell') @@ -172,6 +169,16 @@ contains err = E_UNASSIGNED call set_errmsg("Unknown filter type: " // trim(type_)) end select + + select type(filt => filters(index) % obj) + class is (CppTallyFilter) + filt % ptr = allocate_filter(type) + if (.not. c_associated(filt % ptr)) then + err = E_UNASSIGNED + call set_errmsg("Could not allocate C++ tally filter") + end if + end select + end if else err = E_OUT_OF_BOUNDS diff --git a/src/tallies/tally_filter_azimuthal.F90 b/src/tallies/tally_filter_azimuthal.F90 deleted file mode 100644 index fc7f15a17..000000000 --- a/src/tallies/tally_filter_azimuthal.F90 +++ /dev/null @@ -1,120 +0,0 @@ -module tally_filter_azimuthal - - use, intrinsic :: ISO_C_BINDING - - use algorithm, only: binary_search - use constants - use error, only: fatal_error - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! AZIMUTHALFILTER bins the incident neutron azimuthal angle (relative to the -! global xy-plane). -!=============================================================================== - - type, public, extends(TallyFilter) :: AzimuthalFilter - real(8), allocatable :: bins(:) - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_azimuthal - procedure :: to_statepoint => to_statepoint_azimuthal - procedure :: text_label => text_label_azimuthal - end type AzimuthalFilter - -contains - - subroutine from_xml(this, node) - class(AzimuthalFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: n - integer :: n_angle - real(8) :: d_angle - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - 1 - allocate(this % bins(n)) - call get_node_array(node, "bins", this % bins) - - ! 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 - if (n == 1) then - n_angle = int(this % bins(1)) - if (n_angle > 1) then - this % n_bins = n_angle - d_angle = TWO * PI / n_angle - deallocate(this % bins) - allocate(this % bins(n_angle + 1)) - do i = 1, n_angle - this % bins(i) = -PI + (i - 1) * d_angle - end do - this % bins(n_angle + 1) = PI - else - call fatal_error("Number of bins for azimuthal filter must be& - & greater than 1.") - end if - end if - end subroutine from_xml - - subroutine get_all_bins_azimuthal(this, p, estimator, match) - class(AzimuthalFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: phi - - n = this % n_bins - - ! Make sure the correct direction vector is used. - if (estimator == ESTIMATOR_TRACKLENGTH) then - phi = atan2(p % coord(1) % uvw(2), p % coord(1) % uvw(1)) - else - phi = atan2(p % last_uvw(2), p % last_uvw(1)) - end if - - ! Search to find azimuthal angle bin. - bin = binary_search(this % bins, n + 1, phi) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - - end subroutine get_all_bins_azimuthal - - subroutine to_statepoint_azimuthal(this, filter_group) - class(AzimuthalFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "azimuthal") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_azimuthal - - function text_label_azimuthal(this, bin) result(label) - class(AzimuthalFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - real(8) :: E0, E1 - - E0 = this % bins(bin) - E1 = this % bins(bin + 1) - label = "Azimuthal Angle [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_azimuthal - -end module tally_filter_azimuthal diff --git a/src/tallies/tally_filter_cell.F90 b/src/tallies/tally_filter_cell.F90 deleted file mode 100644 index a394ed122..000000000 --- a/src/tallies/tally_filter_cell.F90 +++ /dev/null @@ -1,144 +0,0 @@ -module tally_filter_cell - - use, intrinsic :: ISO_C_BINDING - - use constants, only: ONE, MAX_LINE_LEN - use dict_header, only: EMPTY - use error, only: fatal_error - use hdf5_interface - use geometry_header - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - public :: openmc_cell_filter_get_bins - -!=============================================================================== -! CELLFILTER specifies which geometric cells tally events reside in. -!=============================================================================== - - type, public, extends(TallyFilter) :: CellFilter - integer(C_INT32_T), allocatable :: cells(:) - type(DictIntInt) :: map - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_cell - procedure :: to_statepoint => to_statepoint_cell - procedure :: text_label => text_label_cell - procedure :: initialize => initialize_cell - end type CellFilter - -contains - - subroutine from_xml(this, node) - class(CellFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - ! Determine how many bins were given - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - allocate(this % cells(n)) - call get_node_array(node, "bins", this % cells) - end subroutine from_xml - - subroutine get_all_bins_cell(this, p, estimator, match) - class(CellFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - integer :: val - - ! Iterate over coordinate levels to see with cells match - do i = 1, p % n_coord - val = this % map % get(p % coord(i) % cell + 1) - if (val /= EMPTY) then - call match % bins % push_back(val) - call match % weights % push_back(ONE) - end if - end do - - end subroutine get_all_bins_cell - - subroutine to_statepoint_cell(this, filter_group) - class(CellFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cell") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id() - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cell - - subroutine initialize_cell(this) - class(CellFilter), intent(inout) :: this - - integer :: i, id - integer :: val - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) - val = cell_dict % get(id) - if (val /= EMPTY) then - this % cells(i) = val - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % set(this % cells(i), i) - end do - end subroutine initialize_cell - - function text_label_cell(this, bin) result(label) - class(CellFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Cell " // to_str(cells(this % cells(bin)) % id()) - end function text_label_cell - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_cell_filter_get_bins(index, cells, n) result(err) bind(C) - ! Return the cells associated with a cell filter - integer(C_INT32_T), value :: index - type(C_PTR), intent(out) :: cells - integer(C_INT32_T), intent(out) :: n - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (CellFilter) - cells = C_LOC(f % cells) - n = size(f % cells) - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to get cells from a non-cell filter.") - end select - end if - end function openmc_cell_filter_get_bins - -end module tally_filter_cell diff --git a/src/tallies/tally_filter_cellborn.F90 b/src/tallies/tally_filter_cellborn.F90 deleted file mode 100644 index 37bf04202..000000000 --- a/src/tallies/tally_filter_cellborn.F90 +++ /dev/null @@ -1,114 +0,0 @@ -module tally_filter_cellborn - - use, intrinsic :: ISO_C_BINDING - - use constants, only: ONE, MAX_LINE_LEN - use dict_header, only: EMPTY - use error, only: fatal_error - use hdf5_interface - use geometry_header - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! CELLBORNFILTER specifies which cell the particle was born in. -!=============================================================================== - - type, public, extends(TallyFilter) :: CellbornFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_cellborn - procedure :: to_statepoint => to_statepoint_cellborn - procedure :: text_label => text_label_cellborn - procedure :: initialize => initialize_cellborn - end type CellbornFilter - -contains - - subroutine from_xml(this, node) - class(CellbornFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - ! Determine number of bins - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - allocate(this % cells(n)) - call get_node_array(node, "bins", this % cells) - end subroutine from_xml - - subroutine get_all_bins_cellborn(this, p, estimator, match) - class(CellbornFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: val - - val = this % map % get(p % cell_born + 1) - if (val /= EMPTY) then - call match % bins % push_back(val) - call match % weights % push_back(ONE) - end if - - end subroutine get_all_bins_cellborn - - subroutine to_statepoint_cellborn(this, filter_group) - class(CellbornFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cellborn") - call write_dataset(filter_group, "n_bins", this % n_bins) - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id() - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cellborn - - subroutine initialize_cellborn(this) - class(CellbornFilter), intent(inout) :: this - - integer :: i, id - integer :: val - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) - val = cell_dict % get(id) - if (val /= EMPTY) then - this % cells(i) = val - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % set(this % cells(i), i) - end do - end subroutine initialize_cellborn - - function text_label_cellborn(this, bin) result(label) - class(CellbornFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Birth Cell " // to_str(cells(this % cells(bin)) % id()) - end function text_label_cellborn - -end module tally_filter_cellborn diff --git a/src/tallies/tally_filter_cellfrom.F90 b/src/tallies/tally_filter_cellfrom.F90 deleted file mode 100644 index a99aa6883..000000000 --- a/src/tallies/tally_filter_cellfrom.F90 +++ /dev/null @@ -1,80 +0,0 @@ -module tally_filter_cellfrom - - use, intrinsic :: ISO_C_BINDING - - use constants, only: ONE, MAX_LINE_LEN - use dict_header, only: EMPTY - use error, only: fatal_error - use hdf5_interface - use geometry_header - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use tally_filter_cell - use xml_interface - - implicit none - private - -!=============================================================================== -! CELLFROMFILTER specifies which geometric cells particles exit when crossing a -! surface. -!=============================================================================== - - type, public, extends(CellFilter) :: CellFromFilter - contains - ! Inherit from_xml from CellFilter - procedure :: get_all_bins => get_all_bins_cell_from - procedure :: to_statepoint => to_statepoint_cell_from - procedure :: text_label => text_label_cell_from - end type CellFromFilter - -contains - - subroutine get_all_bins_cell_from(this, p, estimator, match) - class(CellFromFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - integer :: val - - ! Starting one coordinate level deeper, find the next bin. - do i = 1, p % last_n_coord - val = this % map % get(p % last_cell(i) + 1) - if (val /= EMPTY) then - call match % bins % push_back(val) - call match % weights % push_back(ONE) - exit - end if - end do - - end subroutine get_all_bins_cell_from - - subroutine to_statepoint_cell_from(this, filter_group) - class(CellFromFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cellfrom") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id() - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cell_from - - function text_label_cell_from(this, bin) result(label) - class(CellFromFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Cell from " // to_str(cells(this % cells(bin)) % id()) - end function text_label_cell_from - -end module tally_filter_cellfrom diff --git a/src/tallies/tally_filter_cpp.F90 b/src/tallies/tally_filter_cpp.F90 new file mode 100644 index 000000000..d0e8abfb4 --- /dev/null +++ b/src/tallies/tally_filter_cpp.F90 @@ -0,0 +1,235 @@ +module tally_filter_cpp + use, intrinsic :: ISO_C_BINDING + + use constants, only: MAX_LINE_LEN + use hdf5_interface, only: HID_T + use particle_header, only: Particle + use tally_filter_header + use xml_interface, only: XMLNode + + implicit none + +!=============================================================================== +! CPPTALLYFILTER is a TallyFilter that is at least partially implemented in C++. +! +! Moving a tally to C++ is easier if done piece-by-piece, so this filter has +! *_cpp_inner procedures that allows the Fortran code to override e.g. +! get_all_bins and still call the C++ implementation for side-by-side comparison +!=============================================================================== + + type, abstract, extends(TallyFilter) :: CppTallyFilter + type(C_PTR) :: ptr + contains + procedure :: n_bins_cpp + procedure :: from_xml_cpp_inner + procedure :: get_all_bins_cpp_inner + procedure :: to_statepoint_cpp_inner + procedure :: text_label_cpp_inner + procedure :: initialize_cpp_inner + procedure :: from_xml => from_xml_cpp_default + procedure :: get_all_bins => get_all_bins_cpp_default + procedure :: to_statepoint => to_statepoint_cpp_default + procedure :: text_label => text_label_cpp_default + procedure :: initialize => initialize_cpp_default + end type CppTallyFilter + +!=============================================================================== +! Pure C++ filters +!=============================================================================== + + type, extends(CppTallyFilter) :: AzimuthalFilter + end type + + type, extends(CppTallyFilter) :: CellFilter + end type + + type, extends(CppTallyFilter) :: CellbornFilter + end type + + type, extends(CellFilter) :: CellFromFilter + end type + + type, extends(CppTallyFilter) :: EnergyFunctionFilter + end type + + type, extends(CppTallyFilter) :: MaterialFilter + end type + + type, extends(CppTallyFilter) :: MuFilter + end type + + type, extends(CppTallyFilter) :: PolarFilter + end type + + type, extends(CppTallyFilter) :: SpatialLegendreFilter + end type + + type, extends(CppTallyFilter) :: SurfaceFilter + ! True if this filter is used for surface currents + logical :: current = .false. + end type + + type, extends(CppTallyFilter) :: UniverseFilter + end type + + type, extends(CppTallyFilter) :: ZernikeFilter + end type + + type, extends(ZernikeFilter) :: ZernikeRadialFilter + end type + +contains + +!=============================================================================== +! CppTallyFilter implementation +!=============================================================================== + + function n_bins_cpp(this) result(n_bins) + class(CppTallyFilter), intent(inout) :: this + integer :: n_bins + interface + function filter_n_bins(filt) result(n_bins) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: n_bins + end function filter_n_bins + end interface + n_bins = filter_n_bins(this % ptr) + end function n_bins_cpp + + subroutine from_xml_cpp_inner(this, node) + class(CppTallyFilter), intent(inout) :: this + class(XMLNode), intent(in) :: node + interface + subroutine filter_from_xml(filt, node) bind(C) + import C_PTR + type(C_PTR), value :: filt + type(C_PTR) :: node + end subroutine filter_from_xml + end interface + call filter_from_xml(this % ptr, node % ptr) + end subroutine from_xml_cpp_inner + + subroutine get_all_bins_cpp_inner(this, p, estimator, match) + class(CppTallyFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + interface + subroutine filter_get_all_bins(filt, p, estimator, match) bind(C) + import C_PTR, Particle, C_INT + type(C_PTR), value :: filt + type(Particle), intent(in) :: p + integer(C_INT), intent(in), value :: estimator + type(C_PTR), value :: match + end subroutine filter_get_all_bins + end interface + call filter_get_all_bins(this % ptr, p, estimator, match % ptr) + end subroutine get_all_bins_cpp_inner + + subroutine to_statepoint_cpp_inner(this, filter_group) + class(CppTallyFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + interface + subroutine filter_to_statepoint(filt, filter_group) bind(C) + import C_PTR, HID_T + type(C_PTR), value :: filt + integer(HID_T), intent(in), value :: filter_group + end subroutine filter_to_statepoint + end interface + call filter_to_statepoint(this % ptr, filter_group) + end subroutine to_statepoint_cpp_inner + + function text_label_cpp_inner(this, bin) result(label) + class(CppTallyFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + character(kind=C_CHAR) :: label_(MAX_LINE_LEN+1) + integer :: i + interface + subroutine filter_text_label(filt, bin, label) bind(C) + import C_PTR, C_INT, C_CHAR + type(C_PTR), value :: filt + integer(C_INT), value :: bin + character(kind=C_CHAR) :: label(*) + end subroutine filter_text_label + end interface + call filter_text_label(this % ptr, bin, label_) + label = " " + do i = 1, MAX_LINE_LEN + if (label_(i) == C_NULL_CHAR) exit + label(i:i) = label_(i) + end do + end function text_label_cpp_inner + + subroutine initialize_cpp_inner(this) + class(CppTallyFilter), intent(inout) :: this + interface + subroutine filter_initialize(filt) bind(C) + import C_PTR + type(C_PTR), value :: filt + end subroutine filter_initialize + end interface + call filter_initialize(this % ptr) + end subroutine initialize_cpp_inner + + subroutine from_xml_cpp_default(this, node) + class(CppTallyFilter), intent(inout) :: this + type(XMLNode), intent(in) :: node + call this % from_xml_cpp_inner(node) + this % n_bins = this % n_bins_cpp() + end subroutine from_xml_cpp_default + + subroutine get_all_bins_cpp_default(this, p, estimator, match) + class(CppTallyFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + call this % get_all_bins_cpp_inner(p, estimator, match) + end subroutine get_all_bins_cpp_default + + subroutine to_statepoint_cpp_default(this, filter_group) + class(CppTallyFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + call this % to_statepoint_cpp_inner(filter_group) + end subroutine to_statepoint_cpp_default + + function text_label_cpp_default(this, bin) result(label) + class(CppTallyFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + label = this % text_label_cpp_inner(bin) + end function text_label_cpp_default + + subroutine initialize_cpp_default(this) + class(CppTallyFilter), intent(inout) :: this + call this % initialize_cpp_inner() + end subroutine initialize_cpp_default + +!=============================================================================== +! FILTER_FROM_F given a Fortran index, return a pointer to a C++ filter. +!=============================================================================== + + function filter_from_f(index) result(filt) bind(C) + integer(C_INT32_T), intent(in), value :: index + type(C_PTR) :: filt + + filt = C_NULL_PTR + select type(f => filters(index) % obj) + class is (CppTallyFilter) + filt = f % ptr + end select + end function + +!=============================================================================== +! FILTER_UPDATE_N_BINS given a Fortran index, updates filt % n_bins using C++. +!=============================================================================== + + subroutine filter_update_n_bins(index) bind(C) + integer(C_INT32_T), intent(in), value :: index + select type(f => filters(index) % obj) + class is (CppTallyFilter) + f % n_bins = f % n_bins_cpp() + end select + end subroutine +end module tally_filter_cpp diff --git a/src/tallies/tally_filter_delayedgroup.F90 b/src/tallies/tally_filter_delayedgroup.F90 index 1a662005b..130469a71 100644 --- a/src/tallies/tally_filter_delayedgroup.F90 +++ b/src/tallies/tally_filter_delayedgroup.F90 @@ -62,8 +62,8 @@ contains integer, intent(in) :: estimator type(TallyFilterMatch), intent(inout) :: match - call match % bins % push_back(1) - call match % weights % push_back(ONE) + call match % bins_push_back(1) + call match % weights_push_back(ONE) end subroutine get_all_bins_dg subroutine to_statepoint_dg(this, filter_group) diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index f35ffc3ec..9033ae6a0 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -1,169 +1,21 @@ module tally_filter_distribcell - use, intrinsic :: ISO_C_BINDING - - use constants - use dict_header, only: EMPTY - use error - use geometry_header - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface + use tally_filter_cpp implicit none private -!=============================================================================== -! DISTRIBCELLFILTER specifies which distributed geometric cells tally events -! reside in. -!=============================================================================== - - type, public, extends(TallyFilter) :: DistribcellFilter - integer :: cell + type, public, extends(CppTallyFilter) :: DistribcellFilter contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_distribcell - procedure :: to_statepoint => to_statepoint_distribcell - procedure :: text_label => text_label_distribcell procedure :: initialize => initialize_distribcell end type DistribcellFilter contains - subroutine from_xml(this, node) - class(DistribcellFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - n = node_word_count(node, "bins") - if (n /= 1) call fatal_error("Only one cell can be & - &specified per distribcell filter.") - - ! Store bins - call get_node_value(node, "bins", this % cell) - end subroutine from_xml - - subroutine get_all_bins_distribcell(this, p, estimator, match) - class(DistribcellFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: distribcell_index, offset, i - - distribcell_index = cells(this % cell) % distribcell_index() - offset = 0 - do i = 1, p % n_coord - if (cells(p % coord(i) % cell + 1) % type() == FILL_UNIVERSE) then - offset = offset + cells(p % coord(i) % cell + 1) & - % offset(distribcell_index) - elseif (cells(p % coord(i) % cell + 1) % type() == FILL_LATTICE) then - if (lattices(p % coord(i + 1) % lattice) % are_valid_indices([& - p % coord(i + 1) % lattice_x, & - p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z])) then - offset = offset + lattices(p % coord(i + 1) % lattice) & - % offset(distribcell_index, & - [p % coord(i + 1) % lattice_x, & - p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z]) - end if - end if - if (this % cell == p % coord(i) % cell + 1) then - call match % bins % push_back(offset + 1) - call match % weights % push_back(ONE) - return - end if - end do - end subroutine get_all_bins_distribcell - - subroutine to_statepoint_distribcell(this, filter_group) - class(DistribcellFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "distribcell") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", cells(this % cell) % id()) - end subroutine to_statepoint_distribcell - subroutine initialize_distribcell(this) class(DistribcellFilter), intent(inout) :: this - - integer :: id - integer :: val - - ! Convert id to index. - id = this % cell - val = cell_dict % get(id) - if (val /= EMPTY) then - this % cell = val - this % n_bins = cells(this % cell) % n_instances() - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if + call this % initialize_cpp_inner() + this % n_bins = this % n_bins_cpp() end subroutine initialize_distribcell - function text_label_distribcell(this, bin) result(label) - class(DistribcellFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = '' - call find_offset(this % cell, bin-1, label) - label = "Distributed Cell " // label - end function text_label_distribcell - -!=============================================================================== -! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and -! a target offset to build a string which is the path from the base universe to -! the target cell with the given offset -!=============================================================================== - - subroutine find_offset(i_cell, target_offset, path) - integer, intent(in) :: i_cell ! The target cell index - integer, intent(in) :: target_offset ! Target offset - character(*), intent(inout) :: path ! Path to offset - - integer :: map ! Index in maps vector - integer :: i ! Index over cells - - integer(C_INT) :: path_len - character(kind=C_CHAR), allocatable, target :: path_c(:) - - interface - function distribcell_path_len(target_cell, map, target_offset, root_univ)& - bind(C) result(len) - import C_INT32_T, C_INT - integer(C_INT32_T), intent(in), value :: target_cell, map, & - target_offset, root_univ - integer(C_INT) :: len - end function distribcell_path_len - - subroutine distribcell_path(target_cell, map, target_offset, root_univ, & - path) bind(C) - import C_INT32_T, C_CHAR - integer(C_INT32_T), intent(in), value :: target_cell, map, & - target_offset, root_univ - character(kind=C_CHAR), intent(out) :: path(*) - end subroutine distribcell_path - end interface - - ! Get the distribcell index for this cell - map = cells(i_cell) % distribcell_index() - - path_len = distribcell_path_len(i_cell-1, map, target_offset, & - root_universe) - allocate(path_c(path_len)) - call distribcell_path(i_cell-1, map, target_offset, root_universe, & - path_c) - do i = 1, min(path_len, MAX_LINE_LEN) - if (path_c(i) == C_NULL_CHAR) exit - path(i:i) = path_c(i) - end do - end subroutine find_offset - end module tally_filter_distribcell diff --git a/src/tallies/tally_filter_energy.F90 b/src/tallies/tally_filter_energy.F90 index f230d438e..70142156e 100644 --- a/src/tallies/tally_filter_energy.F90 +++ b/src/tallies/tally_filter_energy.F90 @@ -95,11 +95,11 @@ contains if (p % g /= NONE .and. this % matches_transport_groups) then if (estimator == ESTIMATOR_TRACKLENGTH) then - call match % bins % push_back(num_energy_groups - p % g + 1) - call match % weights % push_back(ONE) + call match % bins_push_back(num_energy_groups - p % g + 1) + call match % weights_push_back(ONE) else - call match % bins % push_back(num_energy_groups - p % last_g + 1) - call match % weights % push_back(ONE) + call match % bins_push_back(num_energy_groups - p % last_g + 1) + call match % weights_push_back(ONE) end if else @@ -109,8 +109,8 @@ contains ! Search to find incoming energy bin. bin = binary_search(this % bins, n + 1, E) if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) + call match % bins_push_back(bin) + call match % weights_push_back(ONE) end if end if end subroutine get_all_bins_energy @@ -155,16 +155,16 @@ contains if (p % g /= NONE .and. this % matches_transport_groups) then ! Tallies are ordered in increasing groups, group indices ! however are the opposite, so switch - call match % bins % push_back(num_energy_groups - p % g + 1) - call match % weights % push_back(ONE) + call match % bins_push_back(num_energy_groups - p % g + 1) + call match % weights_push_back(ONE) else ! Search to find incoming energy bin. bin = binary_search(this % bins, n + 1, p % E) if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) + call match % bins_push_back(bin) + call match % weights_push_back(ONE) end if end if end subroutine get_all_bins_energyout diff --git a/src/tallies/tally_filter_energyfunc.F90 b/src/tallies/tally_filter_energyfunc.F90 deleted file mode 100644 index 1402707b4..000000000 --- a/src/tallies/tally_filter_energyfunc.F90 +++ /dev/null @@ -1,124 +0,0 @@ -module tally_filter_energyfunc - - use, intrinsic :: ISO_C_BINDING - - use algorithm, only: binary_search - use constants - use error, only: fatal_error - use hdf5_interface - use particle_header, only: Particle - use settings, only: run_CE - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! EnergyFunctionFilter multiplies tally scores by an arbitrary function of -! incident energy described by a piecewise linear-linear interpolation. -!=============================================================================== - - type, public, extends(TallyFilter) :: EnergyFunctionFilter - real(8), allocatable :: energy(:) - real(8), allocatable :: y(:) - - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_energyfunction - procedure :: to_statepoint => to_statepoint_energyfunction - procedure :: text_label => text_label_energyfunction - end type EnergyFunctionFilter - -contains - - subroutine from_xml(this, node) - class(EnergyFunctionFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - this % n_bins = 1 - ! Make sure this is continuous-energy mode. - if (.not. run_CE) then - call fatal_error("EnergyFunction filters are only supported for & - &continuous-energy transport calculations") - end if - - ! Allocate and store energy grid. - if (.not. check_for_node(node, "energy")) then - call fatal_error("Energy grid not specified for EnergyFunction & - &filter.") - end if - n = node_word_count(node, "energy") - allocate(this % energy(n)) - call get_node_array(node, "energy", this % energy) - - ! Allocate and store interpolant values. - if (.not. check_for_node(node, "y")) then - call fatal_error("y values not specified for EnergyFunction & - &filter.") - end if - n = node_word_count(node, "y") - allocate(this % y(n)) - call get_node_array(node, "y", this % y) - end subroutine from_xml - - subroutine get_all_bins_energyfunction(this, p, estimator, match) - class(EnergyFunctionFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n, indx - real(8) :: E, f, weight - - select type(this) - type is (EnergyFunctionFilter) - n = size(this % energy) - - ! Get pre-collision energy of particle - E = p % last_E - - ! Search to find incoming energy bin. - indx = binary_search(this % energy, n, E) - - ! Compute an interpolation factor between nearest bins. - f = (E - this % energy(indx)) & - / (this % energy(indx+1) - this % energy(indx)) - - ! Interpolate on the lin-lin grid. - call match % bins % push_back(1) - weight = (ONE - f) * this % y(indx) + f * this % y(indx+1) - call match % weights % push_back(weight) - end select - end subroutine get_all_bins_energyfunction - - subroutine to_statepoint_energyfunction(this, filter_group) - class(EnergyFunctionFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - select type(this) - type is (EnergyFunctionFilter) - call write_dataset(filter_group, "type", "energyfunction") - call write_dataset(filter_group, "energy", this % energy) - call write_dataset(filter_group, "y", this % y) - end select - end subroutine to_statepoint_energyfunction - - function text_label_energyfunction(this, bin) result(label) - class(EnergyFunctionFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - select type(this) - type is (EnergyFunctionFilter) - write(label, FMT="(A, ES8.1, A, ES8.1, A, ES8.1, A, ES8.1, A)") & - "Energy Function f([", this % energy(1), ", ..., ", & - this % energy(size(this % energy)), "]) = [", this % y(1), & - ", ..., ", this % y(size(this % y)), "]" - end select - end function text_label_energyfunction - -end module tally_filter_energyfunc diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index d6aef4abd..4116934b3 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -7,7 +7,6 @@ module tally_filter_header use error use hdf5_interface, only: HID_T use particle_header, only: Particle - use stl_vector, only: VectorInt, VectorReal use string, only: to_str use xml_interface, only: XMLNode @@ -20,19 +19,38 @@ module tally_filter_header public :: openmc_filter_set_id public :: openmc_get_filter_index public :: openmc_get_filter_next_id + public :: filter_match_pointer + + interface + function filter_match_pointer(indx) bind(C) result(ptr) + import C_PTR, C_INT + integer(C_INT), intent(in), value :: indx + type(C_PTR) :: ptr + end function filter_match_pointer + end interface !=============================================================================== ! TALLYFILTERMATCH stores every valid bin and weight for a filter !=============================================================================== type, public :: TallyFilterMatch + type(C_PTR) :: ptr + ! Index of the bin and weight being used in the current filter combination integer :: i_bin - type(VectorInt), pointer :: bins - type(VectorReal), pointer :: weights ! Indicates whether all valid bins for this filter have been found logical :: bins_present = .false. + + contains + procedure :: bins_push_back + procedure :: weights_push_back + procedure :: bins_clear + procedure :: weights_clear + procedure :: bins_size + procedure :: bins_data + procedure :: weights_data + procedure :: bins_set_data end type TallyFilterMatch !=============================================================================== @@ -129,6 +147,116 @@ module tally_filter_header contains +!=============================================================================== +! TallyFilterMatch implementation +!=============================================================================== + + subroutine bins_push_back(this, val) + class(TallyFilterMatch), intent(inout) :: this + integer, intent(in) :: val + interface + subroutine filter_match_bins_push_back(ptr, val) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), intent(in), value :: val + end subroutine + end interface + call filter_match_bins_push_back(this % ptr, val) + end subroutine bins_push_back + + subroutine weights_push_back(this, val) + class(TallyFilterMatch), intent(inout) :: this + real(8), intent(in) :: val + interface + subroutine filter_match_weights_push_back(ptr, val) bind(C) + import C_PTR, C_DOUBLE + type(C_PTR), value :: ptr + real(C_DOUBLE), intent(in), value :: val + end subroutine + end interface + call filter_match_weights_push_back(this % ptr, val) + end subroutine weights_push_back + + subroutine bins_clear(this) + class(TallyFilterMatch), intent(inout) :: this + interface + subroutine filter_match_bins_clear(ptr) bind(C) + import C_PTR + type(C_PTR), value :: ptr + end subroutine + end interface + call filter_match_bins_clear(this % ptr) + end subroutine bins_clear + + subroutine weights_clear(this) + class(TallyFilterMatch), intent(inout) :: this + interface + subroutine filter_match_weights_clear(ptr) bind(C) + import C_PTR + type(C_PTR), value :: ptr + end subroutine + end interface + call filter_match_weights_clear(this % ptr) + end subroutine weights_clear + + function bins_size(this) result(len) + class(TallyFilterMatch), intent(inout) :: this + integer :: len + interface + function filter_match_bins_size(ptr) bind(C) result(len) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT) :: len + end function + end interface + len = filter_match_bins_size(this % ptr) + end function bins_size + + function bins_data(this, indx) result(val) + class(TallyFilterMatch), intent(inout) :: this + integer, intent(in) :: indx + integer :: val + interface + function filter_match_bins_data(ptr, indx) bind(C) result(val) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), intent(in), value :: indx + integer(C_INT) :: val + end function + end interface + val = filter_match_bins_data(this % ptr, indx) + end function bins_data + + function weights_data(this, indx) result(val) + class(TallyFilterMatch), intent(inout) :: this + integer, intent(in) :: indx + real(8) :: val + interface + function filter_match_weights_data(ptr, indx) bind(C) result(val) + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), value :: ptr + integer(C_INT), intent(in), value :: indx + real(C_DOUBLE) :: val + end function + end interface + val = filter_match_weights_data(this % ptr, indx) + end function weights_data + + subroutine bins_set_data(this, indx, val) + class(TallyFilterMatch), intent(inout) :: this + integer, intent(in) :: indx + integer, intent(in) :: val + interface + subroutine filter_match_bins_set_data(ptr, indx, val) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), value, intent(in) :: indx + integer(C_INT), value, intent(in) :: val + end subroutine + end interface + call filter_match_bins_set_data(this % ptr, indx, val) + end subroutine bins_set_data + !=============================================================================== ! INITIALIZE sets up any internal data, as necessary. If this procedure is not ! overriden by the derived class, then it will do nothing by default. @@ -153,8 +281,8 @@ contains ! array is sufficient and a filter object has already been allocated. !=============================================================================== - function verify_filter(index) result(err) - integer(C_INT32_T), intent(in) :: index + function verify_filter(index) result(err) bind(C) + integer(C_INT32_T), intent(in), value :: index integer(C_INT) :: err err = 0 diff --git a/src/tallies/tally_filter_legendre.F90 b/src/tallies/tally_filter_legendre.F90 index 4663df651..9c0f9e41a 100644 --- a/src/tallies/tally_filter_legendre.F90 +++ b/src/tallies/tally_filter_legendre.F90 @@ -2,121 +2,20 @@ module tally_filter_legendre use, intrinsic :: ISO_C_BINDING - use constants - use error - use hdf5_interface - use math, only: calc_pn - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface + use tally_filter_cpp implicit none - private - public :: openmc_legendre_filter_get_order - public :: openmc_legendre_filter_set_order -!=============================================================================== -! LEGENDREFILTER gives Legendre moments of the change in scattering angle -!=============================================================================== + interface + function openmc_legendre_filter_set_order(index, order) result(err) bind(C) + import C_INT32_T, C_INT + integer(C_INT32_T), value :: index + integer(C_INT), value :: order + integer(C_INT) :: err + end function + end interface - type, public, extends(TallyFilter) :: LegendreFilter - integer(C_INT) :: order - contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label + type, extends(CppTallyFilter) :: LegendreFilter end type LegendreFilter -contains - -!=============================================================================== -! LegendreFilter methods -!=============================================================================== - - subroutine from_xml(this, node) - class(LegendreFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - ! Get specified order - call get_node_value(node, "order", this % order) - this % n_bins = this % order + 1 - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(LegendreFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - real(C_DOUBLE) :: wgt(this % n_bins) - - call calc_pn(this % order, p % mu, wgt) - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(wgt(i)) - end do - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(LegendreFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "legendre") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(LegendreFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Legendre expansion, P" // trim(to_str(bin - 1)) - end function text_label - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_legendre_filter_get_order(index, order) result(err) bind(C) - ! Get the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(out) :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (LegendreFilter) - order = f % order - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to get order on a non-expansion filter.") - end select - end if - end function openmc_legendre_filter_get_order - - - function openmc_legendre_filter_set_order(index, order) result(err) bind(C) - ! Set the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), value :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (LegendreFilter) - f % order = order - f % n_bins = order + 1 - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to set order on a non-expansion filter.") - end select - end if - end function openmc_legendre_filter_set_order - end module tally_filter_legendre diff --git a/src/tallies/tally_filter_material.F90 b/src/tallies/tally_filter_material.F90 deleted file mode 100644 index 51ab26594..000000000 --- a/src/tallies/tally_filter_material.F90 +++ /dev/null @@ -1,175 +0,0 @@ -module tally_filter_material - - use, intrinsic :: ISO_C_BINDING - - use constants - use dict_header, only: DictIntInt, EMPTY - use error - use hdf5_interface - use material_header, only: materials, material_dict - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - public :: openmc_material_filter_get_bins - public :: openmc_material_filter_set_bins - -!=============================================================================== -! MATERIAL specifies which material tally events reside in. -!=============================================================================== - - type, public, extends(TallyFilter) :: MaterialFilter - integer, allocatable :: materials(:) - type(DictIntInt) :: map - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_material - procedure :: to_statepoint => to_statepoint_material - procedure :: text_label => text_label_material - procedure :: initialize => initialize_material - end type MaterialFilter - -contains - - subroutine from_xml(this, node) - class(MaterialFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - allocate(this % materials(n)) - call get_node_array(node, "bins", this % materials) - end subroutine from_xml - - subroutine get_all_bins_material(this, p, estimator, match) - class(MaterialFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: val - - val = this % map % get(p % material) - if (val /= EMPTY) then - call match % bins % push_back(val) - call match % weights % push_back(ONE) - end if - - end subroutine get_all_bins_material - - subroutine to_statepoint_material(this, filter_group) - class(MaterialFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: material_ids(:) - - call write_dataset(filter_group, "type", "material") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(material_ids(size(this % materials))) - do i = 1, size(this % materials) - material_ids(i) = materials(this % materials(i)) % id() - end do - call write_dataset(filter_group, "bins", material_ids) - end subroutine to_statepoint_material - - subroutine initialize_material(this) - class(MaterialFilter), intent(inout) :: this - - integer :: i, id - integer :: val - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % materials(i) - val = material_dict % get(id) - if (val /= EMPTY) then - this % materials(i) = val - else - call fatal_error("Could not find material " // trim(to_str(id)) & - &// " specified on a tally filter.") - end if - end do - - ! Generate mapping from material indices to filter bins. - do i = 1, this % n_bins - call this % map % set(this % materials(i), i) - end do - end subroutine initialize_material - - function text_label_material(this, bin) result(label) - class(MaterialFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Material " // to_str(materials(this % materials(bin)) % id()) - end function text_label_material - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_material_filter_get_bins(index, bins, n) result(err) bind(C) - ! Return the bins for a material filter - integer(C_INT32_T), value :: index - type(C_PTR), intent(out) :: bins - integer(C_INT32_T), intent(out) :: n - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (MaterialFilter) - bins = C_LOC(f % materials) - n = size(f % materials) - err = 0 - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to get material filter bins on a & - &non-material filter.") - end select - end if - end function openmc_material_filter_get_bins - - - function openmc_material_filter_set_bins(index, n, bins) result(err) bind(C) - ! Set the materials for the filter - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: n - integer(C_INT32_T), intent(in) :: bins(n) - integer(C_INT) :: err - - integer :: i - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (MaterialFilter) - f % n_bins = n - if (allocated(f % materials)) deallocate(f % materials) - allocate(f % materials(n)) - f % materials(:) = bins - - ! Generate mapping from material indices to filter bins. - call f % map % clear() - do i = 1, n - call f % map % set(f % materials(i), i) - end do - - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to set material filter bins on a & - &non-material filter.") - end select - end if - end function openmc_material_filter_set_bins - -end module tally_filter_material diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index dda0bc0ce..8f6dd185e 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -2,197 +2,37 @@ module tally_filter_mesh use, intrinsic :: ISO_C_BINDING - use constants - use dict_header, only: EMPTY - use error - use mesh_header - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface + use tally_filter_cpp implicit none - private - public :: openmc_mesh_filter_get_mesh - public :: openmc_mesh_filter_set_mesh -!=============================================================================== -! MESHFILTER indexes the location of particle events to a regular mesh. For -! tracklength tallies, it will produce multiple valid bins and the bin weight -! will correspond to the fraction of the track length that lies in that bin. -!=============================================================================== + interface + function openmc_mesh_filter_set_mesh(index, index_mesh) result(err) bind(C) + import C_INT32_T, C_INT + integer(C_INT32_T), value, intent(in) :: index + integer(C_INT32_T), value, intent(in) :: index_mesh + integer(C_INT) :: err + end function + end interface - type, public, extends(TallyFilter) :: MeshFilter - integer :: mesh + type, public, extends(CppTallyFilter) :: MeshFilter contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_mesh - procedure :: to_statepoint => to_statepoint_mesh - procedure :: text_label => text_label_mesh + procedure :: mesh => get_mesh end type MeshFilter contains - subroutine from_xml(this, node) - class(MeshFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: id - integer :: n - integer(C_INT) :: err - type(RegularMesh) :: m - - n = node_word_count(node, "bins") - - if (n /= 1) call fatal_error("Only one mesh can be & - &specified per mesh filter.") - - ! Determine id of mesh - call get_node_value(node, "bins", id) - - ! Get pointer to mesh - err = openmc_get_mesh_index(id, this % mesh) - if (err /= 0) then - call fatal_error("Could not find mesh " // trim(to_str(id)) & - // " specified on filter.") - end if - - ! Determine number of bins - m = meshes(this % mesh) - this % n_bins = 1 - do i = 1, m % n_dimension() - this % n_bins = this % n_bins * m % dimension(i) - end do - end subroutine from_xml - - subroutine get_all_bins_mesh(this, p, estimator, match) - class(MeshFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: bin - type(RegularMesh) :: m - type(C_PTR) :: ptr_bins, ptr_weights - + function get_mesh(this) result(mesh) + class(MeshFilter) :: this + integer :: mesh interface - subroutine mesh_bins_crossed(m, p, bins, weights) bind(C) - import C_PTR, Particle - type(C_PTR), value :: m - type(Particle), intent(in) :: p - type(C_PTR), value :: bins - type(C_PTR), value :: weights - end subroutine + function mesh_filter_get_mesh(filt) result(index_mesh) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: index_mesh + end function mesh_filter_get_mesh end interface - - ! Get a pointer to the mesh. - m = meshes(this % mesh) - - if (estimator /= ESTIMATOR_TRACKLENGTH) then - ! If this is an analog or collision tally, then there can only be one - ! valid mesh bin. - call m % get_bin(p % coord(1) % xyz, bin) - if (bin >= 0) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - return - else - ptr_bins = C_LOC(match % bins) - ptr_weights = C_LOC(match % weights) - call mesh_bins_crossed(m % ptr, p, ptr_bins, ptr_weights) - end if - - end subroutine get_all_bins_mesh - - subroutine to_statepoint_mesh(this, filter_group) - class(MeshFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - type(RegularMesh) :: m - - m = meshes(this % mesh) - call write_dataset(filter_group, "type", "mesh") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", m % id()) - end subroutine to_statepoint_mesh - - function text_label_mesh(this, bin) result(label) - class(MeshFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer, allocatable :: ijk(:) - type(RegularMesh) :: m - - m = meshes(this % mesh) - allocate(ijk(m % n_dimension())) - call m % get_indices_from_bin(bin, ijk) - if (m % n_dimension() == 1) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" - elseif (m % n_dimension() == 2) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" - elseif (m % n_dimension() == 3) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" - end if - end function text_label_mesh - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_mesh_filter_get_mesh(index, index_mesh) result(err) bind(C) - ! Get the mesh for a mesh filter - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), intent(out) :: index_mesh - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (MeshFilter) - index_mesh = f % mesh - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to set mesh on a non-mesh filter.") - end select - end if - end function openmc_mesh_filter_get_mesh - - - function openmc_mesh_filter_set_mesh(index, index_mesh) result(err) bind(C) - ! Set the mesh for a mesh filter - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: index_mesh - integer(C_INT) :: err - - type(RegularMesh) :: m - integer :: i - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (MeshFilter) - if (index_mesh >= 0 .and. index_mesh < n_meshes()) then - f % mesh = index_mesh - f % n_bins = 1 - m = meshes(index_mesh) - do i = 1, m % n_dimension() - f % n_bins = f % n_bins * m % dimension(i) - end do - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in 'meshes' array is out of bounds.") - end if - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to set mesh on a non-mesh filter.") - end select - end if - end function openmc_mesh_filter_set_mesh + mesh = mesh_filter_get_mesh(this % ptr) + end function get_mesh end module tally_filter_mesh diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index f08bc35ef..f70f5425f 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -2,224 +2,21 @@ module tally_filter_meshsurface use, intrinsic :: ISO_C_BINDING - use constants - use error - use mesh_header - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface + use tally_filter_cpp implicit none - private - public :: openmc_meshsurface_filter_get_mesh - public :: openmc_meshsurface_filter_set_mesh -!=============================================================================== -! MESHFILTER indexes the location of particle events to a regular mesh. For -! tracklength tallies, it will produce multiple valid bins and the bin weight -! will correspond to the fraction of the track length that lies in that bin. -!=============================================================================== + interface + function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) & + bind(C) + import C_INT32_T, C_INT + integer(C_INT32_T), value, intent(in) :: index + integer(C_INT32_T), value, intent(in) :: index_mesh + integer(C_INT) :: err + end function + end interface - type, public, extends(TallyFilter) :: MeshSurfaceFilter - integer :: mesh - contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label + type, extends(CppTallyFilter) :: MeshSurfaceFilter end type MeshSurfaceFilter -contains - - subroutine from_xml(this, node) - class(MeshSurfaceFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: id - integer :: n - integer(C_INT) :: err - type(RegularMesh) :: m - - n = node_word_count(node, "bins") - - if (n /= 1) call fatal_error("Only one mesh can be & - &specified per meshsurface filter.") - - ! Determine id of mesh - call get_node_value(node, "bins", id) - - ! Get pointer to mesh - err = openmc_get_mesh_index(id, this % mesh) - if (err /= 0) then - call fatal_error("Could not find mesh " // trim(to_str(id)) & - // " specified on filter.") - end if - - ! Determine number of bins - m = meshes(this % mesh) - this % n_bins = 4 * m % n_dimension() - do i = 1, m % n_dimension() - this % n_bins = this % n_bins * m % dimension(i) - end do - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(MeshSurfaceFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - type(RegularMesh) :: m - type(C_PTR) :: ptr_bins, ptr_weights - - interface - subroutine mesh_surface_bins_crossed(m, p, bins, weights) bind(C) - import C_PTR, Particle - type(C_PTR), value :: m - type(Particle), intent(in) :: p - type(C_PTR), value :: bins - type(C_PTR), value :: weights - end subroutine - end interface - - ! Get a pointer to the mesh. - m = meshes(this % mesh) - - ptr_bins = C_LOC(match % bins) - ptr_weights = C_LOC(match % weights) - call mesh_surface_bins_crossed(m % ptr, p, ptr_bins, ptr_weights) - - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(MeshSurfaceFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - type(RegularMesh) :: m - - m = meshes(this % mesh) - call write_dataset(filter_group, "type", "meshsurface") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", m % id()) - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(MeshSurfaceFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer :: i_mesh - integer :: i_surf - integer :: n_dim - integer, allocatable :: ijk(:) - type(RegularMesh) :: m - - m = meshes(this % mesh) - n_dim = m % n_dimension() - allocate(ijk(n_dim)) - - ! Get flattend mesh index and surface index - i_mesh = (bin - 1) / (4*n_dim) + 1 - i_surf = mod(bin - 1, 4*n_dim) + 1 - - ! Get mesh index part of label - call m % get_indices_from_bin(i_mesh, ijk) - if (m % n_dimension() == 1) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" - elseif (m % n_dimension() == 2) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" - elseif (m % n_dimension() == 3) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" - end if - - ! Get surface part of label - select case (i_surf) - case (OUT_LEFT) - label = trim(label) // " Outgoing, x-min" - case (IN_LEFT) - label = trim(label) // " Incoming, x-min" - case (OUT_RIGHT) - label = trim(label) // " Outgoing, x-max" - case (IN_RIGHT) - label = trim(label) // " Incoming, x-max" - case (OUT_BACK) - label = trim(label) // " Outgoing, y-min" - case (IN_BACK) - label = trim(label) // " Incoming, y-min" - case (OUT_FRONT) - label = trim(label) // " Outgoing, y-max" - case (IN_FRONT) - label = trim(label) // " Incoming, y-max" - case (OUT_BOTTOM) - label = trim(label) // " Outgoing, z-min" - case (IN_BOTTOM) - label = trim(label) // " Incoming, z-min" - case (OUT_TOP) - label = trim(label) // " Outgoing, z-max" - case (IN_TOP) - label = trim(label) // " Incoming, z-max" - end select - end function text_label - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_meshsurface_filter_get_mesh(index, index_mesh) result(err) bind(C) - ! Get the mesh for a mesh surface filter - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), intent(out) :: index_mesh - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (MeshSurfaceFilter) - index_mesh = f % mesh - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to set mesh on a non-mesh filter.") - end select - end if - end function openmc_meshsurface_filter_get_mesh - - - function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) bind(C) - ! Set the mesh for a mesh surface filter - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: index_mesh - integer(C_INT) :: err - - integer :: i - integer :: n_dim - type(RegularMesh) :: m - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (MeshSurfaceFilter) - if (index_mesh >= 0 .and. index_mesh < n_meshes()) then - f % mesh = index_mesh - m = meshes(index_mesh) - n_dim = m % n_dimension() - f % n_bins = 4*n_dim - do i = 1, n_dim - f % n_bins = f % n_bins * m % dimension(i) - end do - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in 'meshes' array is out of bounds.") - end if - class default - err = E_INVALID_TYPE - call set_errmsg("Tried to set mesh on a non-mesh filter.") - end select - end if - end function openmc_meshsurface_filter_set_mesh - end module tally_filter_meshsurface diff --git a/src/tallies/tally_filter_mu.F90 b/src/tallies/tally_filter_mu.F90 deleted file mode 100644 index 7f8bde923..000000000 --- a/src/tallies/tally_filter_mu.F90 +++ /dev/null @@ -1,110 +0,0 @@ -module tally_filter_mu - - use, intrinsic :: ISO_C_BINDING - - use algorithm, only: binary_search - use constants, only: ONE, TWO, MAX_LINE_LEN, NO_BIN_FOUND - use error, only: fatal_error - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! MUFILTER bins the incoming-outgoing direction cosine. This is only used for -! scatter reactions. -!=============================================================================== - - type, public, extends(TallyFilter) :: MuFilter - real(8), allocatable :: bins(:) - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_mu - procedure :: to_statepoint => to_statepoint_mu - procedure :: text_label => text_label_mu - end type MuFilter - -contains - - subroutine from_xml(this, node) - class(MuFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: n_angle - integer :: n - real(8) :: d_angle - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - 1 - allocate(this % bins(n)) - call get_node_array(node, "bins", this % bins) - - ! 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 - if (n == 1) then - n_angle = int(this % bins(1)) - if (n_angle > 1) then - this % n_bins = n_angle - d_angle = TWO / n_angle - deallocate(this % bins) - allocate(this % bins(n_angle + 1)) - do i = 1, n_angle - this % bins(i) = -ONE + (i - 1) * d_angle - end do - this % bins(n_angle + 1) = ONE - else - call fatal_error("Number of bins for mu filter must be& - & greater than 1.") - end if - end if - end subroutine from_xml - - subroutine get_all_bins_mu(this, p, estimator, match) - class(MuFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - - n = this % n_bins - - ! Search to find incoming energy bin. - bin = binary_search(this % bins, n + 1, p % mu) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - end subroutine get_all_bins_mu - - subroutine to_statepoint_mu(this, filter_group) - class(MuFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "mu") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_mu - - function text_label_mu(this, bin) result(label) - class(MuFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - real(8) :: E0, E1 - - E0 = this % bins(bin) - E1 = this % bins(bin + 1) - label = "Change-in-Angle [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_mu - -end module tally_filter_mu diff --git a/src/tallies/tally_filter_particle.F90 b/src/tallies/tally_filter_particle.F90 index 5ac846310..665ea5cf3 100644 --- a/src/tallies/tally_filter_particle.F90 +++ b/src/tallies/tally_filter_particle.F90 @@ -52,8 +52,8 @@ contains do i = 1, this % n_bins if (this % particles(i) == p % type) then - call match % bins % push_back(i) - call match % weights % push_back(ONE) + call match % bins_push_back(i) + call match % weights_push_back(ONE) end if end do end subroutine get_all_bins diff --git a/src/tallies/tally_filter_polar.F90 b/src/tallies/tally_filter_polar.F90 deleted file mode 100644 index 7cbb55f8b..000000000 --- a/src/tallies/tally_filter_polar.F90 +++ /dev/null @@ -1,118 +0,0 @@ -module tally_filter_polar - - use, intrinsic :: ISO_C_BINDING - - use algorithm, only: binary_search - use constants - use error, only: fatal_error - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! POLARFILTER bins the incident neutron polar angle (relative to the global -! z-axis). -!=============================================================================== - - type, public, extends(TallyFilter) :: PolarFilter - real(8), allocatable :: bins(:) - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_polar - procedure :: to_statepoint => to_statepoint_polar - procedure :: text_label => text_label_polar - end type PolarFilter - -contains - - subroutine from_xml(this, node) - class(PolarFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: n_angle - integer :: n - real(8) :: d_angle - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - 1 - allocate(this % bins(n)) - call get_node_array(node, "bins", this % bins) - - ! 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 - if (n == 1) then - n_angle = int(this % bins(1)) - if (n_angle > 1) then - this % n_bins = n_angle - d_angle = PI / real(n_angle,8) - deallocate(this % bins) - allocate(this % bins(n_angle + 1)) - do i = 1, n_angle - this % bins(i) = (i - 1) * d_angle - end do - this % bins(n_angle + 1) = PI - else - call fatal_error("Number of bins for polar filter must be& - & greater than 1.") - end if - end if - end subroutine from_xml - - subroutine get_all_bins_polar(this, p, estimator, match) - class(PolarFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: theta - - n = this % n_bins - - ! Make sure the correct direction vector is used. - if (estimator == ESTIMATOR_TRACKLENGTH) then - theta = acos(p % coord(1) % uvw(3)) - else - theta = acos(p % last_uvw(3)) - end if - - ! Search to find polar angle bin. - bin = binary_search(this % bins, n + 1, theta) - if (bin /= NO_BIN_FOUND) then - call match % bins % push_back(bin) - call match % weights % push_back(ONE) - end if - end subroutine get_all_bins_polar - - subroutine to_statepoint_polar(this, filter_group) - class(PolarFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "polar") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_polar - - function text_label_polar(this, bin) result(label) - class(PolarFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - real(8) :: E0, E1 - - E0 = this % bins(bin) - E1 = this % bins(bin + 1) - label = "Polar Angle [" // trim(to_str(E0)) // ", " // trim(to_str(E1)) & - // ")" - end function text_label_polar - -end module tally_filter_polar diff --git a/src/tallies/tally_filter_sph_harm.F90 b/src/tallies/tally_filter_sph_harm.F90 index 4a1f43273..455da042c 100644 --- a/src/tallies/tally_filter_sph_harm.F90 +++ b/src/tallies/tally_filter_sph_harm.F90 @@ -2,241 +2,32 @@ module tally_filter_sph_harm use, intrinsic :: ISO_C_BINDING - use constants - use error - use hdf5_interface - use math, only: calc_pn, calc_rn - use particle_header, only: Particle - use string, only: to_str, to_lower, to_f_string - use tally_filter_header - use xml_interface + use tally_filter_cpp implicit none private - public :: openmc_sphharm_filter_get_order - public :: openmc_sphharm_filter_get_cosine - public :: openmc_sphharm_filter_set_order - public :: openmc_sphharm_filter_set_cosine integer, public, parameter :: COSINE_SCATTER = 1 integer, public, parameter :: COSINE_PARTICLE = 2 -!=============================================================================== -! SPHERICALHARMONICSFILTER gives spherical harmonics expansion moments of a -! tally score -!=============================================================================== - - type, public, extends(TallyFilter) :: SphericalHarmonicsFilter - integer :: order - integer :: cosine + type, public, extends(CppTallyFilter) :: SphericalHarmonicsFilter contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label + procedure :: cosine end type SphericalHarmonicsFilter contains -!=============================================================================== -! SphericalHarmonicsFilter methods -!=============================================================================== - - subroutine from_xml(this, node) - class(SphericalHarmonicsFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - character(MAX_WORD_LEN) :: temp_str - - ! Get specified order - call get_node_value(node, "order", this % order) - this % n_bins = (this % order + 1)**2 - - ! Determine how cosine term is to be treated - if (check_for_node(node, "cosine")) then - call get_node_value(node, "cosine", temp_str) - select case (to_lower(temp_str)) - case ('scatter') - this % cosine = COSINE_SCATTER - case ('particle') - this % cosine = COSINE_PARTICLE - end select - else - this % cosine = COSINE_PARTICLE - end if - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(SphericalHarmonicsFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i, j, n - integer :: num_nm - real(C_DOUBLE) :: wgt(this % order + 1) - real(C_DOUBLE) :: rn(this % n_bins) - - ! Determine cosine term for scatter expansion if necessary - if (this % cosine == COSINE_SCATTER) then - call calc_pn(this % order, p % mu, wgt) - else - wgt = ONE - end if - - ! Find the Rn,m values - call calc_rn(this % order, p % last_uvw, rn) - - j = 0 - do n = 0, this % order - ! Calculate n-th order spherical harmonics for (u,v,w) - num_nm = 2*n + 1 - - ! Append matching (bin,weight) for each moment - do i = 1, num_nm - j = j + 1 - call match % bins % push_back(j) - call match % weights % push_back(wgt(n + 1) * rn(j)) - end do - end do - - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(SphericalHarmonicsFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "sphericalharmonics") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - if (this % cosine == COSINE_SCATTER) then - call write_dataset(filter_group, "cosine", "scatter") - else - call write_dataset(filter_group, "cosine", "particle") - end if - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(SphericalHarmonicsFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer :: n, m - - do n = 0, this % order - if (bin <= (n + 1)**2) then - m = (bin - n**2 - 1) - n - label = "Spherical harmonic expansion, Y" // trim(to_str(n)) // & - "," // trim(to_str(m)) - exit - end if - end do - end function text_label - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_sphharm_filter_get_order(index, order) result(err) bind(C) - ! Get the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(out) :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (SphericalHarmonicsFilter) - order = f % order - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spherical harmonics filter.") - end select - end if - end function openmc_sphharm_filter_get_order - - - function openmc_sphharm_filter_get_cosine(index, cosine) result(err) bind(C) - ! Get the order of an expansion filter - integer(C_INT32_T), value :: index - character(kind=C_CHAR), intent(out) :: cosine(*) - integer(C_INT) :: err - - integer :: i - character(10) :: cosine_ - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (SphericalHarmonicsFilter) - select case (f % cosine) - case (COSINE_SCATTER) - cosine_ = 'scatter' - case (COSINE_PARTICLE) - cosine_ = 'particle' - end select - - ! Convert to C string - do i = 1, len_trim(cosine_) - cosine(i) = cosine_(i:i) - end do - cosine(len_trim(cosine_) + 1) = C_NULL_CHAR - - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spherical harmonics filter.") - end select - end if - end function openmc_sphharm_filter_get_cosine - - - function openmc_sphharm_filter_set_order(index, order) result(err) bind(C) - ! Set the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), value :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (SphericalHarmonicsFilter) - f % order = order - f % n_bins = (order + 1)**2 - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spherical harmonics filter.") - end select - end if - end function openmc_sphharm_filter_set_order - - - function openmc_sphharm_filter_set_cosine(index, cosine) result(err) bind(C) - ! Set the cosine parameter - integer(C_INT32_T), value :: index - character(kind=C_CHAR), intent(in) :: cosine(*) - integer(C_INT) :: err - - character(:), allocatable :: cosine_ - - ! Convert C string to Fortran string - cosine_ = to_f_string(cosine) - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (SphericalHarmonicsFilter) - select case (cosine_) - case ('scatter') - f % cosine = COSINE_SCATTER - case ('particle') - f % cosine = COSINE_PARTICLE - end select - - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spherical harmonics filter.") - end select - end if - end function openmc_sphharm_filter_set_cosine + function cosine(this) result(val) + class(SphericalHarmonicsFilter) :: this + integer :: val + interface + function sphharm_filter_get_cosine(filt) result(val) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: val + end function + end interface + val = sphharm_filter_get_cosine(this % ptr) + end function cosine end module tally_filter_sph_harm diff --git a/src/tallies/tally_filter_sptl_legendre.F90 b/src/tallies/tally_filter_sptl_legendre.F90 deleted file mode 100644 index 1bfbd0e3b..000000000 --- a/src/tallies/tally_filter_sptl_legendre.F90 +++ /dev/null @@ -1,225 +0,0 @@ -module tally_filter_sptl_legendre - - use, intrinsic :: ISO_C_BINDING - - use constants - use error - use hdf5_interface - use math, only: calc_pn - use particle_header, only: Particle - use string, only: to_str, to_lower - use tally_filter_header - use xml_interface - - implicit none - private - public :: openmc_spatial_legendre_filter_get_order - public :: openmc_spatial_legendre_filter_get_params - public :: openmc_spatial_legendre_filter_set_order - public :: openmc_spatial_legendre_filter_set_params - - integer, parameter :: AXIS_X = 1 - integer, parameter :: AXIS_Y = 2 - integer, parameter :: AXIS_Z = 3 - -!=============================================================================== -! SPATIALLEGENDREFILTER gives Legendre moments of the particle's normalized -! position along an axis -!=============================================================================== - - type, public, extends(TallyFilter) :: SpatialLegendreFilter - integer(C_INT) :: order - integer(C_INT) :: axis - real(C_DOUBLE) :: min - real(C_DOUBLE) :: max - contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label - end type SpatialLegendreFilter - -contains - -!=============================================================================== -! SpatialLegendreFilter methods -!=============================================================================== - - subroutine from_xml(this, node) - class(SpatialLegendreFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - character(MAX_WORD_LEN) :: axis - - ! Get attributes from XML - call get_node_value(node, "order", this % order) - call get_node_value(node, "axis", axis) - select case (to_lower(axis)) - case ('x') - this % axis = AXIS_X - case ('y') - this % axis = AXIS_Y - case ('z') - this % axis = AXIS_Z - end select - call get_node_value(node, "min", this % min) - call get_node_value(node, "max", this % max) - - this % n_bins = this % order + 1 - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(SpatialLegendreFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - real(C_DOUBLE) :: wgt(this % n_bins) - real(C_DOUBLE) :: x ! Position on specified axis - real(C_DOUBLE) :: x_norm ! Normalized position - - x = p % coord(1) % xyz(this % axis) - if (this % min <= x .and. x <= this % max) then - ! Calculate normalized position between min and max - x_norm = TWO*(x - this % min)/(this % max - this % min) - ONE - - call calc_pn(this % order, x_norm, wgt) - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(wgt(i)) - end do - end if - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(SpatialLegendreFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - character(kind=C_CHAR) :: axis - - call write_dataset(filter_group, "type", "spatiallegendre") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - select case (this % axis) - case (AXIS_X) - axis = 'x' - case (AXIS_Y) - axis = 'y' - case (AXIS_Z) - axis = 'z' - end select - call write_dataset(filter_group, 'axis', axis) - call write_dataset(filter_group, 'min', this % min) - call write_dataset(filter_group, 'max', this % max) - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(SpatialLegendreFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - character(1) :: axis - - select case (this % axis) - case (AXIS_X) - axis = 'x' - case (AXIS_Y) - axis = 'y' - case (AXIS_Z) - axis = 'z' - end select - label = "Legendre expansion, " // axis // " axis, P" // & - trim(to_str(bin - 1)) - end function text_label - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_spatial_legendre_filter_get_order(index, order) result(err) bind(C) - ! Get the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(out) :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (SpatialLegendreFilter) - order = f % order - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spatial Legendre filter.") - end select - end if - end function openmc_spatial_legendre_filter_get_order - - - function openmc_spatial_legendre_filter_set_order(index, order) result(err) bind(C) - ! Set the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), value :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (SpatialLegendreFilter) - f % order = order - f % n_bins = order + 1 - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spatial Legendre filter.") - end select - end if - end function openmc_spatial_legendre_filter_set_order - - - function openmc_spatial_legendre_filter_get_params(index, axis, min, max) & - result(err) bind(C) - ! Get the parameters for a spatial Legendre filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(out) :: axis - real(C_DOUBLE), intent(out) :: min - real(C_DOUBLE), intent(out) :: max - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type(f => filters(index) % obj) - type is (SpatialLegendreFilter) - axis = f % axis - min = f % min - max = f % max - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spatial Legendre filter.") - end select - end if - end function openmc_spatial_legendre_filter_get_params - - - function openmc_spatial_legendre_filter_set_params(index, axis, min, max) & - result(err) bind(C) - ! Set the parameters for a spatial Legendre filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(in), optional :: axis - real(C_DOUBLE), intent(in), optional :: min - real(C_DOUBLE), intent(in), optional :: max - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type(f => filters(index) % obj) - type is (SpatialLegendreFilter) - if (present(axis)) f % axis = axis - if (present(min)) f % min = min - if (present(max)) f % max = max - class default - err = E_INVALID_TYPE - call set_errmsg("Not a spatial Legendre filter.") - end select - end if - end function openmc_spatial_legendre_filter_set_params - -end module tally_filter_sptl_legendre diff --git a/src/tallies/tally_filter_surface.F90 b/src/tallies/tally_filter_surface.F90 deleted file mode 100644 index ec642bc8d..000000000 --- a/src/tallies/tally_filter_surface.F90 +++ /dev/null @@ -1,109 +0,0 @@ -module tally_filter_surface - - use, intrinsic :: ISO_C_BINDING - - use constants, only: ONE, MAX_LINE_LEN - use dict_header, only: EMPTY - use error, only: fatal_error - use hdf5_interface - use surface_header - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! SURFACEFILTER specifies which surface particles are crossing -!=============================================================================== - - type, public, extends(TallyFilter) :: SurfaceFilter - integer, allocatable :: surfaces(:) - - ! True if this filter is used for surface currents - logical :: current = .false. - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_surface - procedure :: to_statepoint => to_statepoint_surface - procedure :: text_label => text_label_surface - procedure :: initialize => initialize_surface - end type SurfaceFilter - -contains - - subroutine from_xml(this, node) - class(SurfaceFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - allocate(this % surfaces(n)) - call get_node_array(node, "bins", this % surfaces) - end subroutine from_xml - - subroutine get_all_bins_surface(this, p, estimator, match) - class(SurfaceFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - - do i = 1, this % n_bins - if (abs(p % surface) == this % surfaces(i)) then - call match % bins % push_back(i) - if (p % surface < 0) then - call match % weights % push_back(-ONE) - else - call match % weights % push_back(ONE) - end if - exit - end if - end do - - end subroutine get_all_bins_surface - - subroutine to_statepoint_surface(this, filter_group) - class(SurfaceFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "surface") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % surfaces) - end subroutine to_statepoint_surface - - subroutine initialize_surface(this) - class(SurfaceFilter), intent(inout) :: this - - integer :: i, id - integer :: val - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % surfaces(i) - val = surface_dict % get(id) - if (val /= EMPTY) then - this % surfaces(i) = val - else - call fatal_error("Could not find surface " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - end subroutine initialize_surface - - function text_label_surface(this, bin) result(label) - class(SurfaceFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Surface " // to_str(surfaces(this % surfaces(bin)) % id()) - end function text_label_surface - -end module tally_filter_surface diff --git a/src/tallies/tally_filter_universe.F90 b/src/tallies/tally_filter_universe.F90 deleted file mode 100644 index 9e0e09fcf..000000000 --- a/src/tallies/tally_filter_universe.F90 +++ /dev/null @@ -1,118 +0,0 @@ -module tally_filter_universe - - use, intrinsic :: ISO_C_BINDING - - use constants, only: ONE, MAX_LINE_LEN - use dict_header, only: EMPTY - use error, only: fatal_error - use hdf5_interface - use geometry_header - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! UNIVERSEFILTER specifies which geometric universes tally events reside in. -!=============================================================================== - - type, public, extends(TallyFilter) :: UniverseFilter - integer, allocatable :: universes(:) - type(DictIntInt) :: map - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_universe - procedure :: to_statepoint => to_statepoint_universe - procedure :: text_label => text_label_universe - procedure :: initialize => initialize_universe - end type UniverseFilter - -contains - - subroutine from_xml(this, node) - class(UniverseFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - allocate(this % universes(n)) - call get_node_array(node, "bins", this % universes) - end subroutine from_xml - - subroutine get_all_bins_universe(this, p, estimator, match) - class(UniverseFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - integer :: val - - ! Iterate over coordinate levels to see which universes match - do i = 1, p % n_coord - val = this % map % get(p % coord(i) % universe) - if (val /= EMPTY) then - call match % bins % push_back(val) - call match % weights % push_back(ONE) - end if - end do - - end subroutine get_all_bins_universe - - subroutine to_statepoint_universe(this, filter_group) - class(UniverseFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: universe_ids(:) - - call write_dataset(filter_group, "type", "universe") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(universe_ids(size(this % universes))) - do i = 1, size(this % universes) - universe_ids(i) = universe_id(this % universes(i)) - end do - call write_dataset(filter_group, "bins", universe_ids) - end subroutine to_statepoint_universe - - subroutine initialize_universe(this) - class(UniverseFilter), intent(inout) :: this - - integer :: i, id - integer :: val - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % universes(i) - val = universe_dict % get(id) - if (val /= EMPTY) then - this % universes(i) = val - else - call fatal_error("Could not find universe " // trim(to_str(id)) & - &// " specified on a tally filter.") - end if - end do - - ! Generate mapping from universe indices to filter bins. - do i = 1, this % n_bins - call this % map % set(this % universes(i), i) - end do - end subroutine initialize_universe - - function text_label_universe(this, bin) result(label) - class(UniverseFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Universe " // to_str(universe_id(this % universes(bin))) - end function text_label_universe - -end module tally_filter_universe diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 deleted file mode 100644 index 04515cf31..000000000 --- a/src/tallies/tally_filter_zernike.F90 +++ /dev/null @@ -1,295 +0,0 @@ -module tally_filter_zernike - - use, intrinsic :: ISO_C_BINDING - - use constants - use error - use hdf5_interface - use math, only: calc_zn, calc_zn_rad - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - public :: openmc_zernike_filter_get_order - public :: openmc_zernike_filter_get_params - public :: openmc_zernike_filter_set_order - public :: openmc_zernike_filter_set_params - -!=============================================================================== -! ZERNIKEFILTER gives Zernike polynomial moments of a particle's position -!=============================================================================== - - type, public, extends(TallyFilter) :: ZernikeFilter - integer :: order - real(8) :: x - real(8) :: y - real(8) :: r - contains - procedure :: calc_n_bins => calc_n_bins_zn - procedure :: from_xml => from_xml_zn - procedure :: get_all_bins => get_all_bins_zn - procedure :: to_statepoint => to_statepoint_zn - procedure :: text_label => text_label_zn - end type ZernikeFilter - -!=============================================================================== -! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a -! particle's position -!=============================================================================== - - type, public, extends(ZernikeFilter) :: ZernikeRadialFilter - contains - procedure :: calc_n_bins => calc_n_bins_zn_rad - ! Inherit from_xml from ZernikeFilter - procedure :: get_all_bins => get_all_bins_zn_rad - procedure :: to_statepoint => to_statepoint_zn_rad - procedure :: text_label => text_label_zn_rad - end type ZernikeRadialFilter - -contains - -!=============================================================================== -! ZernikeFilter methods -!=============================================================================== - - function calc_n_bins_zn(this) result(n_bins) - class(ZernikeFilter), intent(in) :: this - integer :: n - integer :: n_bins - - n = this % order - n_bins = ((n+1) * (n+2))/2 - end function calc_n_bins_zn - - subroutine from_xml_zn(this, node) - class(ZernikeFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - ! Get center of cylinder and radius - call get_node_value(node, "x", this % x) - call get_node_value(node, "y", this % y) - call get_node_value(node, "r", this % r) - - ! Get specified order - call get_node_value(node, "order", n) - this % order = n - this % n_bins = this % calc_n_bins() - end subroutine from_xml_zn - - subroutine get_all_bins_zn(this, p, estimator, match) - class(ZernikeFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - real(8) :: x, y, r, theta - real(C_DOUBLE) :: zn(this % n_bins) - - ! Determine normalized (r,theta) positions - x = p % coord(1) % xyz(1) - this % x - y = p % coord(1) % xyz(2) - this % y - r = sqrt(x*x + y*y)/this % r - if (r <= ONE) then - theta = atan2(y, x) - - ! Get moments for Zernike polynomial orders 0..n - call calc_zn(this % order, r, theta, zn) - - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(zn(i)) - end do - endif - end subroutine get_all_bins_zn - - subroutine to_statepoint_zn(this, filter_group) - class(ZernikeFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "zernike") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - call write_dataset(filter_group, "x", this % x) - call write_dataset(filter_group, "y", this % y) - call write_dataset(filter_group, "r", this % r) - end subroutine to_statepoint_zn - - function text_label_zn(this, bin) result(label) - class(ZernikeFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer :: n, m - integer :: first, last - - do n = 0, this % order - last = (n + 1)*(n + 2)/2 - if (bin <= last) then - first = last - n - m = -n + (bin - first)*2 - label = "Zernike expansion, Z" // trim(to_str(n)) // "," & - // trim(to_str(m)) - exit - end if - end do - end function text_label_zn - -!=============================================================================== -! ZernikeRadialFilter methods -!=============================================================================== - - function calc_n_bins_zn_rad(this) result(num_n_bins) - class(ZernikeRadialFilter), intent(in) :: this - integer :: n - integer :: num_n_bins - - n = this % order - num_n_bins = n/2 + 1 - end function calc_n_bins_zn_rad - - subroutine get_all_bins_zn_rad(this, p, estimator, match) - class(ZernikeRadialFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - real(8) :: x, y, r - real(C_DOUBLE) :: zn_rad(this % n_bins) - - ! Determine normalized (r,theta) positions - x = p % coord(1) % xyz(1) - this % x - y = p % coord(1) % xyz(2) - this % y - r = sqrt(x*x + y*y)/this % r - if (r <= ONE) then - - ! Get moments for even order Zernike polynomial orders 0..n - call calc_zn_rad(this % order, r, zn_rad) - - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(zn_rad(i)) - end do - endif - end subroutine get_all_bins_zn_rad - - subroutine to_statepoint_zn_rad(this, filter_group) - class(ZernikeRadialFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "zernikeradial") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - call write_dataset(filter_group, "x", this % x) - call write_dataset(filter_group, "y", this % y) - call write_dataset(filter_group, "r", this % r) - end subroutine to_statepoint_zn_rad - - function text_label_zn_rad(this, bin) result(label) - class(ZernikeRadialFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Zernike expansion, Z" // trim(to_str(2*(bin-1))) // ",0" - end function text_label_zn_rad - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_zernike_filter_get_order(index, order) result(err) bind(C) - ! Get the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(out) :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeFilter) - order = f % order - type is (ZernikeRadialFilter) - order = f % order - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike filter.") - end select - end if - end function openmc_zernike_filter_get_order - - - function openmc_zernike_filter_get_params(index, x, y, r) result(err) bind(C) - ! Get the Zernike filter parameters - integer(C_INT32_T), value :: index - real(C_DOUBLE), intent(out) :: x - real(C_DOUBLE), intent(out) :: y - real(C_DOUBLE), intent(out) :: r - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - class is (ZernikeFilter) - x = f % x - y = f % y - r = f % r - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike filter.") - end select - end if - end function openmc_zernike_filter_get_params - - - function openmc_zernike_filter_set_order(index, order) result(err) bind(C) - ! Set the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), value :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeFilter) - f % order = order - f % n_bins = ((order + 1)*(order + 2))/2 - type is (ZernikeRadialFilter) - f % order = order - f % n_bins = order/2 + 1 - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike filter.") - end select - end if - end function openmc_zernike_filter_set_order - - - function openmc_zernike_filter_set_params(index, x, y, r) result(err) bind(C) - ! Set the Zernike filter parameters - integer(C_INT32_T), value :: index - real(C_DOUBLE), intent(in), optional :: x - real(C_DOUBLE), intent(in), optional :: y - real(C_DOUBLE), intent(in), optional :: r - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - class is (ZernikeFilter) - if (present(x)) f % x = x - if (present(y)) f % y = y - if (present(r)) f % r = r - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike filter.") - end select - end if - end function openmc_zernike_filter_set_params - -end module tally_filter_zernike diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index d468ac7b7..0b80255d5 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -344,7 +344,7 @@ contains this % estimator = ESTIMATOR_ANALOG type is (SphericalHarmonicsFilter) j = FILTER_SPH_HARMONICS - if (filt % cosine == COSINE_SCATTER) then + if (filt % cosine() == COSINE_SCATTER) then this % estimator = ESTIMATOR_ANALOG end if type is (SpatialLegendreFilter) @@ -413,6 +413,13 @@ contains !=============================================================================== subroutine free_memory_tally() + interface + subroutine free_memory_tally_c() bind(C) + end subroutine free_memory_tally_c + end interface + + call free_memory_tally_c() + n_tallies = 0 if (allocated(tallies)) deallocate(tallies) call tally_dict % clear() diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 715b3c206..f5efd0437 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -173,8 +173,8 @@ contains ! Initialize bins, filter level do j = 1, size(t % filter) - call filter_matches(t % filter(j)) % bins % clear() - call filter_matches(t % filter(j)) % bins % push_back(0) + call filter_matches(t % filter(j)) % bins_clear() + call filter_matches(t % filter(j)) % bins_push_back(0) end do FILTER_LOOP: do filter_index = 1, t % n_filter_bins @@ -264,13 +264,13 @@ contains i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE)) select type(filt => filters(i_filter_mesh) % obj) type is (MeshFilter) - m = meshes(filt % mesh) + m = meshes(filt % mesh()) end select ! initialize bins array do j = 1, size(t % filter) - call filter_matches(t % filter(j)) % bins % clear() - call filter_matches(t % filter(j)) % bins % push_back(1) + call filter_matches(t % filter(j)) % bins_clear() + call filter_matches(t % filter(j)) % bins_push_back(1) end do ! determine how many energyin bins there are @@ -296,20 +296,20 @@ contains ! Get the indices for this cell call m % get_indices_from_bin(i, ijk) - filter_matches(i_filter_mesh) % bins % data(1) = i + call filter_matches(i_filter_mesh) % bins_set_data(1, i) do l = 1, n if (print_ebin) then - filter_matches(i_filter_ein) % bins % data(1) = l + call filter_matches(i_filter_ein) % bins_set_data(1, l) end if ! Left Surface - filter_matches(i_filter_surf) % bins % data(1) = OUT_LEFT + call filter_matches(i_filter_surf) % bins_set_data(1, OUT_LEFT) filter_index = 1 do j = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins % data(1) - 1) * t % stride(j) + bins_data(1) - 1) * t % stride(j) end do call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) if (trigger % std_dev < std_dev) then @@ -321,11 +321,11 @@ contains trigger % variance = std_dev**2 ! Right Surface - filter_matches(i_filter_surf) % bins % data(1) = OUT_RIGHT + call filter_matches(i_filter_surf) % bins_set_data(1, OUT_RIGHT) filter_index = 1 do j = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins % data(1) - 1) * t % stride(j) + bins_data(1) - 1) * t % stride(j) end do call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) if (trigger % std_dev < std_dev) then @@ -337,11 +337,11 @@ contains trigger % variance = trigger % std_dev**2 ! Back Surface - filter_matches(i_filter_surf) % bins % data(1) = OUT_BACK + call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BACK) filter_index = 1 do j = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins % data(1) - 1) * t % stride(j) + bins_data(1) - 1) * t % stride(j) end do call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) if (trigger % std_dev < std_dev) then @@ -353,11 +353,11 @@ contains trigger % variance = trigger % std_dev**2 ! Front Surface - filter_matches(i_filter_surf) % bins % data(1) = OUT_FRONT + call filter_matches(i_filter_surf) % bins_set_data(1, OUT_FRONT) filter_index = 1 do j = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins % data(1) - 1) * t % stride(j) + bins_data(1) - 1) * t % stride(j) end do call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) if (trigger % std_dev < std_dev) then @@ -369,11 +369,11 @@ contains trigger % variance = trigger % std_dev**2 ! Bottom Surface - filter_matches(i_filter_surf) % bins % data(1) = OUT_BOTTOM + call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BOTTOM) filter_index = 1 do j = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins % data(1) - 1) * t % stride(j) + bins_data(1) - 1) * t % stride(j) end do call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) if (trigger % std_dev < std_dev) then @@ -385,11 +385,11 @@ contains trigger % variance = trigger % std_dev**2 ! Top Surface - filter_matches(i_filter_surf) % bins % data(1) = OUT_TOP + call filter_matches(i_filter_surf) % bins_set_data(1, OUT_TOP) filter_index = 1 do j = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % filter(j)) % & - bins % data(1) - 1) * t % stride(j) + bins_data(1) - 1) * t % stride(j) end do call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t) if (trigger % std_dev < std_dev) then diff --git a/tests/regression_tests/dagmc/test.py b/tests/regression_tests/dagmc/test.py index 943245d27..9817cff95 100644 --- a/tests/regression_tests/dagmc/test.py +++ b/tests/regression_tests/dagmc/test.py @@ -1,7 +1,7 @@ from tests.testing_harness import TestHarness import os import pytest -import openmc +import openmc.capi pytestmark = pytest.mark.skipif( not openmc.capi.dagmc_enabled, @@ -9,4 +9,4 @@ pytestmark = pytest.mark.skipif( def test_dagmc(): harness = TestHarness('statepoint.5.h5') - harness.main() + harness.main() diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index ae2d92719..059557dbf 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -168,6 +168,8 @@ def test_tally(capi_init): openmc.capi.MaterialFilter(uid=1) mats = openmc.capi.materials f = openmc.capi.MaterialFilter([mats[2], mats[1]]) + assert f.bins[0] == mats[2] + assert f.bins[1] == mats[1] t.filters = [f] assert t.filters == [f] diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index ff430ab2c..eed9e23c4 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -14,6 +14,6 @@ fi # Download multipole library if [[ ! -e $HOME/WMP_Library/092235.h5 ]]; then - wget -q https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/WMP_Library_v1.0.tar.gz - tar -C $HOME -xzf WMP_Library_v1.0.tar.gz + wget -q https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/WMP_Library_v1.1.tar.gz + tar -C $HOME -xzf WMP_Library_v1.1.tar.gz fi