diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6fc130947b..580d409c5d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,4 +1,4 @@ -name: CI +name: Tests and Coverage on: # allows us to run workflows manually @@ -21,7 +21,25 @@ env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} jobs: + filter-changes: + runs-on: ubuntu-latest + outputs: + source_changed: ${{ steps.filter.outputs.source_changed }} + steps: + - name: Check out the repository + uses: actions/checkout@v4 + - name: Examine changed files + id: filter + uses: dorny/paths-filter@668c092af3649c4b664c54e4b704aa46782f6f7c # latest master commit, not released yet + with: + filters: | + source_changed: + - '!docs/**' + - '!**/*.md' + predicate-quantifier: 'every' main: + needs: filter-changes + if: ${{ needs.filter-changes.outputs.source_changed == 'true' }} runs-on: ubuntu-22.04 strategy: matrix: @@ -205,12 +223,33 @@ jobs: flag-name: C++ and Python path-to-lcov: coverage.lcov - finish: - needs: main + coverage: + needs: [filter-changes, main] + if: ${{ always() }} runs-on: ubuntu-latest steps: - name: Coveralls Finished + if: ${{ needs.filter-changes.outputs.source_changed == 'true' }} uses: coverallsapp/github-action@v2 with: github-token: ${{ secrets.GITHUB_TOKEN }} parallel-finished: true + + ci-pass: + needs: [filter-changes, main, coverage] + name: Check CI status + if: ${{ always() }} + runs-on: ubuntu-latest + steps: + - name: Check CI status + run: | + if [[ "${{ needs.filter-changes.outputs.source_changed }}" == "false" ]]; then + echo "Documentation-only change - CI skipped successfully" + exit 0 + fi + if [[ "${{ needs.main.result }}" == "success" && "${{ needs.coverage.result }}" == "success" ]]; then + echo "CI passed" + exit 0 + fi + echo "CI failed" + exit 1 diff --git a/.gitignore b/.gitignore index 3b2a24a4db..780059f307 100644 --- a/.gitignore +++ b/.gitignore @@ -30,7 +30,7 @@ docs/source/_images/*.aux docs/source/pythonapi/generated/ # Source build -build +build*/ # build from src/utils/setup.py src/utils/build diff --git a/AGENTS.md b/AGENTS.md index 32e30213ae..44962d1ade 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -31,6 +31,24 @@ C++17/Python codebase where: - **Depletion**: `openmc.deplete` implements burnup via operator-splitting with various integrators (Predictor, CECM, etc.) - **Nuclear Data**: `openmc.data` provides programmatic access to nuclear data files (ENDF, ACE, HDF5) +## Git Branching Workflow + +OpenMC uses a git flow branching model with two primary branches: + +- **`develop` branch**: The main development branch where all ongoing development takes place. This is the **primary branch against which pull requests are submitted and merged**. This branch is not guaranteed to be stable and may contain work-in-progress features. +- **`master` branch**: The stable release branch containing the latest stable release of OpenMC. This branch only receives merges from `develop` when the development team decides a release should occur. + +### Instructions for Code Review + +When analyzing code changes on a feature or bugfix branch (e.g., when a user asks "what do you think of these changes?"), **compare the branch changes against `develop`, not `master`**. Pull requests are submitted to merge into `develop`, so differences relative to `develop` represent the actual proposed changes. Comparing against `master` will include unrelated changes from other features that have already been merged to `develop`. + +### Workflow for contributors + +1. Create a feature/bugfix branch off `develop` +2. Make changes and commit to the feature branch +3. Open a pull request to merge the feature branch into `develop` +4. A committer reviews and merges the PR into `develop` + ## Critical Build & Test Workflows ### Build Dependencies diff --git a/CMakeLists.txt b/CMakeLists.txt index 87b8789d10..d3119fb875 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -372,6 +372,7 @@ list(APPEND libopenmc_SOURCES src/particle.cpp src/particle_data.cpp src/particle_restart.cpp + src/particle_type.cpp src/photon.cpp src/physics.cpp src/physics_common.cpp diff --git a/docs/source/io_formats/collision_track.rst b/docs/source/io_formats/collision_track.rst index a364597544..8e1e00ffb8 100644 --- a/docs/source/io_formats/collision_track.rst +++ b/docs/source/io_formats/collision_track.rst @@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested (``collision_track.N.h5``) or when the run is performed in parallel. The file contains the information needed to reconstruct each recorded collision. -The current revision of the collision track file format is 1.0. +The current revision of the collision track file format is 1.1. **/** @@ -37,9 +37,9 @@ The current revision of the collision track file format is 1.0. - ``material_id`` (*int*) -- ID of the material containing the collision site. - ``universe_id`` (*int*) -- ID of the universe containing the collision site. - ``n_collision`` (*int*) -- Collision counter for the particle history. - - ``particle`` (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, 3=positron). - - ``parent_id`` (*int64*) -- Unique ID of the parent particle. - - ``progeny_id`` (*int64*) -- Progeny ID of the particle. + - ``particle`` (*int32_t*) -- Particle type (PDG number). + - ``parent_id`` (*int64_t*) -- Unique ID of the parent particle. + - ``progeny_id`` (*int64_t*) -- Progeny ID of the particle. In an MPI run, OpenMC writes the combined dataset by gathering collision-track entries from all ranks before flushing them to disk, so the final file appears diff --git a/docs/source/io_formats/depletion_results.rst b/docs/source/io_formats/depletion_results.rst index b2a726dad0..7fb088268c 100644 --- a/docs/source/io_formats/depletion_results.rst +++ b/docs/source/io_formats/depletion_results.rst @@ -34,6 +34,7 @@ The current version of the depletion results file format is 1.2. :Attributes: - **index** (*int*) -- Index used in results for this material - **volume** (*double*) -- Volume of this material in [cm^3] + - **name** (*char[]*) -- Name of this material **/nuclides//** diff --git a/docs/source/io_formats/particle_restart.rst b/docs/source/io_formats/particle_restart.rst index 2734f0470e..e0fe76b9d8 100644 --- a/docs/source/io_formats/particle_restart.rst +++ b/docs/source/io_formats/particle_restart.rst @@ -4,7 +4,7 @@ Particle Restart File Format ============================ -The current version of the particle restart file format is 2.0. +The current version of the particle restart file format is 2.1. **/** @@ -26,8 +26,7 @@ The current version of the particle restart file format is 2.0. - **run_mode** (*char[]*) -- Run mode used, either 'fixed source', 'eigenvalue', or 'particle restart'. - **id** (*int8_t*) -- Unique identifier of the particle. - - **type** (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, - 3=positron) + - **type** (*int32_t*) -- Particle type (PDG number) - **weight** (*double*) -- Weight of the particle. - **energy** (*double*) -- Energy of the particle in eV for continuous-energy mode, or the energy group of the particle for diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 8fc26ae070..ed7c6273aa 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -721,7 +721,10 @@ attributes/sub-elements: is present. :particle: - The source particle type, either ``neutron`` or ``photon``. + The source particle type, specified as a PDG number or a string alias (e.g., + ``neutron``/``n``, ``photon``/``gamma``, ``electron``, ``positron``, + ``proton``/``p``, ``deuteron``/``d``, ``triton``/``t``, ``alpha``, or GNDS + nuclide names like ``Fe57``). *Default*: neutron @@ -1537,7 +1540,8 @@ sub-elements/attributes: *Default*: None :particle_type: - The particle that the weight windows will apply to (e.g., 'neutron') + The particle that the weight windows will apply to, specified as a PDG + code or string (e.g., ``neutron``). *Default*: 'neutron' @@ -1597,7 +1601,8 @@ mesh-based weight windows. *Default*: None :particle_type: - The particle that the weight windows will apply to (e.g., 'neutron') + The particle that the weight windows will apply to, specified as a PDG + code or string (e.g., ``neutron``). *Default*: neutron diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index 4ce2229ed8..2e81e0a3f2 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -15,6 +15,8 @@ following the same format. **/** :Attributes: - **filetype** (*char[]*) -- String indicating the type of file. + - **version** (*int[2]*) -- Major and minor version of the source + file format. :Datasets: @@ -22,5 +24,5 @@ following the same format. particle. The compound type has fields ``r``, ``u``, ``E``, ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``, which represent the position, direction, energy, time, weight, - delayed group, surface ID, and particle type (0=neutron, 1=photon, - 2=electron, 3=positron), respectively. + delayed group, surface ID, and particle type (PDG number), + respectively. diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 2309643dc8..7d7765b849 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -4,7 +4,7 @@ State Point File Format ======================= -The current version of the statepoint file format is 18.1. +The current version of the statepoint file format is 18.2. **/** @@ -56,8 +56,8 @@ The current version of the statepoint file format is 18.1. ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and ``particle``, which represent the position, direction, energy, time, weight, delayed group, surface ID, and particle type - (0=neutron, 1=photon, 2=electron, 3=positron), respectively. Only - present when `run_mode` is 'eigenvalue'. + (PDG number), respectively. Only present when `run_mode` is + 'eigenvalue'. **/tallies/** diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 9f29a949ac..0ba0e061fd 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -318,8 +318,8 @@ should be set to: they use ``energy`` and ``y``. :particle: - A list of integers indicating the type of particles to tally ('neutron' = 1, - 'photon' = 2, 'electron' = 3, 'positron' = 4). + A list of particle identifiers to tally, specified as strings (e.g., + ``neutron``, ``photon``, ``He4``) or as integer PDG numbers. ------------------ ```` Element diff --git a/docs/source/io_formats/track.rst b/docs/source/io_formats/track.rst index a97d75e58e..9c05ac4c37 100644 --- a/docs/source/io_formats/track.rst +++ b/docs/source/io_formats/track.rst @@ -4,7 +4,7 @@ Track File Format ================= -The current revision of the particle track file format is 3.0. +The current revision of the particle track file format is 3.1. **/** @@ -32,6 +32,5 @@ The current revision of the particle track file format is 3.0. the array for each primary/secondary particle. The last offset should match the total size of the array. - - **particles** (*int[]*) -- Particle type for each - primary/secondary particle (0=neutron, 1=photon, - 2=electron, 3=positron). + - **particles** (*int32_t[]*) -- Particle type for + each primary/secondary particle (PDG number). diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 382381a9ee..64ef7de6d5 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -513,6 +513,7 @@ Supported scores: - total - fission - nu-fission + - kappa-fission - events Supported Estimators: diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index b9c693ba1c..5a04fedd70 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -400,7 +400,7 @@ below. { openmc::SourceSite particle; // weight - particle.particle = openmc::ParticleType::neutron; + particle.particle = openmc::ParticleType::neutron(); particle.wgt = 1.0; // position double angle = 2.0 * M_PI * openmc::prn(seed); @@ -477,7 +477,7 @@ parameters to the source class when it is created: { openmc::SourceSite particle; // weight - particle.particle = openmc::ParticleType::neutron; + particle.particle = openmc::ParticleType::neutron(); particle.wgt = 1.0; // position particle.r.x = 0.0; diff --git a/examples/custom_source/source_ring.cpp b/examples/custom_source/source_ring.cpp index b87afbd176..5ab531a4a7 100644 --- a/examples/custom_source/source_ring.cpp +++ b/examples/custom_source/source_ring.cpp @@ -1,17 +1,16 @@ -#include // for M_PI +#include // for M_PI #include // for unique_ptr +#include "openmc/particle.h" #include "openmc/random_lcg.h" #include "openmc/source.h" -#include "openmc/particle.h" -class RingSource : public openmc::Source -{ +class RingSource : public openmc::Source { openmc::SourceSite sample(uint64_t* seed) const { openmc::SourceSite particle; // particle type - particle.particle = openmc::ParticleType::neutron; + particle.particle = openmc::ParticleType::neutron(); // position double angle = 2.0 * M_PI * openmc::prn(seed); double radius = 3.0; @@ -25,10 +24,11 @@ class RingSource : public openmc::Source } }; -// A function to create a unique pointer to an instance of this class when generated -// via a plugin call using dlopen/dlsym. -// You must have external C linkage here otherwise dlopen will not find the file -extern "C" std::unique_ptr openmc_create_source(std::string parameters) +// A function to create a unique pointer to an instance of this class when +// generated via a plugin call using dlopen/dlsym. You must have external C +// linkage here otherwise dlopen will not find the file +extern "C" std::unique_ptr openmc_create_source( + std::string parameters) { return std::make_unique(); } diff --git a/examples/parameterized_custom_source/parameterized_source_ring.cpp b/examples/parameterized_custom_source/parameterized_source_ring.cpp index e70dc8bb19..a983414f1c 100644 --- a/examples/parameterized_custom_source/parameterized_source_ring.cpp +++ b/examples/parameterized_custom_source/parameterized_source_ring.cpp @@ -1,63 +1,65 @@ -#include // for M_PI +#include // for M_PI #include // for unique_ptr #include +#include "openmc/particle.h" #include "openmc/random_lcg.h" #include "openmc/source.h" -#include "openmc/particle.h" class RingSource : public openmc::Source { - public: - RingSource(double radius, double energy) : radius_(radius), energy_(energy) { } +public: + RingSource(double radius, double energy) : radius_(radius), energy_(energy) {} - // Defines a function that can create a unique pointer to a new instance of this class - // by extracting the parameters from the provided string. - static std::unique_ptr from_string(std::string parameters) - { - std::unordered_map parameter_mapping; + // Defines a function that can create a unique pointer to a new instance of + // this class by extracting the parameters from the provided string. + static std::unique_ptr from_string(std::string parameters) + { + std::unordered_map parameter_mapping; - std::stringstream ss(parameters); - std::string parameter; - while (std::getline(ss, parameter, ',')) { - parameter.erase(0, parameter.find_first_not_of(' ')); - std::string key = parameter.substr(0, parameter.find_first_of('=')); - std::string value = parameter.substr(parameter.find_first_of('=') + 1, parameter.length()); - parameter_mapping[key] = value; - } - - double radius = std::stod(parameter_mapping["radius"]); - double energy = std::stod(parameter_mapping["energy"]); - return std::make_unique(radius, energy); + std::stringstream ss(parameters); + std::string parameter; + while (std::getline(ss, parameter, ',')) { + parameter.erase(0, parameter.find_first_not_of(' ')); + std::string key = parameter.substr(0, parameter.find_first_of('=')); + std::string value = + parameter.substr(parameter.find_first_of('=') + 1, parameter.length()); + parameter_mapping[key] = value; } - // Samples from an instance of this class. - openmc::SourceSite sample(uint64_t* seed) const - { - openmc::SourceSite particle; - // particle type - particle.particle = openmc::ParticleType::neutron; - // position - double angle = 2.0 * M_PI * openmc::prn(seed); - double radius = this->radius_; - particle.r.x = radius * std::cos(angle); - particle.r.y = radius * std::sin(angle); - particle.r.z = 0.0; - // angle - particle.u = {1.0, 0.0, 0.0}; - particle.E = this->energy_; + double radius = std::stod(parameter_mapping["radius"]); + double energy = std::stod(parameter_mapping["energy"]); + return std::make_unique(radius, energy); + } - return particle; - } + // Samples from an instance of this class. + openmc::SourceSite sample(uint64_t* seed) const + { + openmc::SourceSite particle; + // particle type + particle.particle = openmc::ParticleType::neutron(); + // position + double angle = 2.0 * M_PI * openmc::prn(seed); + double radius = this->radius_; + particle.r.x = radius * std::cos(angle); + particle.r.y = radius * std::sin(angle); + particle.r.z = 0.0; + // angle + particle.u = {1.0, 0.0, 0.0}; + particle.E = this->energy_; - private: - double radius_; - double energy_; + return particle; + } + +private: + double radius_; + double energy_; }; -// A function to create a unique pointer to an instance of this class when generated -// via a plugin call using dlopen/dlsym. -// You must have external C linkage here otherwise dlopen will not find the file -extern "C" std::unique_ptr openmc_create_source(std::string parameters) +// A function to create a unique pointer to an instance of this class when +// generated via a plugin call using dlopen/dlsym. You must have external C +// linkage here otherwise dlopen will not find the file +extern "C" std::unique_ptr openmc_create_source( + std::string parameters) { return RingSource::from_string(parameters); } diff --git a/include/openmc/capi.h b/include/openmc/capi.h index d8041ef414..019a418a12 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -116,7 +116,7 @@ int openmc_mesh_set_id(int32_t index, int32_t id); int openmc_mesh_get_n_elements(int32_t index, size_t* n); int openmc_mesh_get_volumes(int32_t index, double* volumes); int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz, - int max_mats, int32_t* materials, double* volumes); + int max_mats, int32_t* materials, double* volumes, double* bboxes); int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh); int openmc_new_filter(const char* type, int32_t* index); @@ -201,8 +201,8 @@ int openmc_weight_windows_set_energy_bounds( int32_t index, double* e_bounds, size_t e_bounds_size); int openmc_weight_windows_get_energy_bounds( int32_t index, const double** e_bounds, size_t* e_bounds_size); -int openmc_weight_windows_set_particle(int32_t index, int particle); -int openmc_weight_windows_get_particle(int32_t index, int* particle); +int openmc_weight_windows_set_particle(int32_t index, int32_t particle); +int openmc_weight_windows_get_particle(int32_t index, int32_t* particle); int openmc_weight_windows_get_bounds(int32_t index, const double** lower_bounds, const double** upper_bounds, size_t* size); int openmc_weight_windows_set_bounds(int32_t index, const double* lower_bounds, @@ -227,7 +227,7 @@ int openmc_zernike_filter_set_order(int32_t index, int order); int openmc_zernike_filter_set_params( int32_t index, const double* x, const double* y, const double* r); -int openmc_particle_filter_get_bins(int32_t idx, int bins[]); +int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]); //! Sets the mesh and energy grid for CMFD reweight //! \param[in] meshtyally_id id of CMFD Mesh Tally diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 9dfe2339ab..26c34ebd41 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -123,11 +123,11 @@ private: //! BoundingBox if the particle is in a complex cell. BoundingBox bounding_box_complex(vector postfix) const; - //! Enfource precedence: Parenthases, Complement, Intersection, Union - void add_precedence(); + //! Enforce precedence between intersections and unions + void enforce_precedence(); //! Add parenthesis to enforce precedence - int64_t add_parentheses(int64_t start); + void add_parentheses(int64_t start); //! Remove complement operators from the expression void remove_complement_ops(); diff --git a/include/openmc/constants.h b/include/openmc/constants.h index bba260c3a4..980cff447e 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -25,16 +25,16 @@ using double_4dvec = vector>>>; constexpr int HDF5_VERSION[] {3, 0}; // Version numbers for binary files -constexpr array VERSION_STATEPOINT {18, 1}; -constexpr array VERSION_PARTICLE_RESTART {2, 0}; -constexpr array VERSION_TRACK {3, 0}; +constexpr array VERSION_STATEPOINT {18, 2}; +constexpr array VERSION_PARTICLE_RESTART {2, 1}; +constexpr array VERSION_TRACK {3, 1}; constexpr array VERSION_SUMMARY {6, 1}; constexpr array VERSION_VOLUME {1, 0}; constexpr array VERSION_VOXEL {2, 0}; constexpr array VERSION_MGXS_LIBRARY {1, 0}; constexpr array VERSION_PROPERTIES {1, 1}; constexpr array VERSION_WEIGHT_WINDOWS {1, 0}; -constexpr array VERSION_COLLISION_TRACK {1, 0}; +constexpr array VERSION_COLLISION_TRACK {1, 1}; // ============================================================================ // ADJUSTABLE PARAMETERS diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 26f106fca1..3e31b17a94 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -167,6 +167,7 @@ private: class SpatialBox : public SpatialDistribution { public: explicit SpatialBox(pugi::xml_node node, bool fission = false); + SpatialBox(Position lower_left, Position upper_right, bool fission = false); //! Sample a position from the distribution //! \param seed Pseudorandom number seed pointer diff --git a/include/openmc/material.h b/include/openmc/material.h index e36946c71b..c10f25551e 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -14,6 +14,7 @@ #include "openmc/memory.h" // for unique_ptr #include "openmc/ncrystal_interface.h" #include "openmc/particle.h" +#include "openmc/settings.h" #include "openmc/vector.h" namespace openmc { @@ -110,9 +111,12 @@ public: //! \return Density in [atom/b-cm] double density() const { return density_; } - //! Get density in [g/cm^3] + //! Get density in [g/cm^3]. //! \return Density in [g/cm^3] - double density_gpcc() const { return density_gpcc_; } + double density_gpcc() const + { + return settings::run_CE ? density_gpcc_ : density(); + } //! Get charge density in [e/b-cm] //! \return Charge density in [e/b-cm] diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 7ceb623dae..5dc327fe35 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -4,6 +4,7 @@ #ifndef OPENMC_MESH_H #define OPENMC_MESH_H +#include #include #include "hdf5.h" @@ -87,8 +88,12 @@ namespace detail { class MaterialVolumes { public: + MaterialVolumes(int32_t* mats, double* vols, double* bboxes, int table_size) + : materials_(mats), volumes_(vols), bboxes_(bboxes), table_size_(table_size) + {} + MaterialVolumes(int32_t* mats, double* vols, int table_size) - : materials_(mats), volumes_(vols), table_size_(table_size) + : MaterialVolumes(mats, vols, nullptr, table_size) {} //! Add volume for a given material in a mesh element @@ -96,8 +101,11 @@ public: //! \param[in] index_elem Index of the mesh element //! \param[in] index_material Index of the material within the model //! \param[in] volume Volume to add - void add_volume(int index_elem, int index_material, double volume); - void add_volume_unsafe(int index_elem, int index_material, double volume); + //! \param[in] bbox Bounding box to union into the result (optional) + void add_volume(int index_elem, int index_material, double volume, + const BoundingBox* bbox = nullptr); + void add_volume_unsafe(int index_elem, int index_material, double volume, + const BoundingBox* bbox = nullptr); // Accessors int32_t& materials(int i, int j) { return materials_[i * table_size_ + j]; } @@ -112,11 +120,23 @@ public: return volumes_[i * table_size_ + j]; } + double& bboxes(int i, int j, int k) + { + return bboxes_[(i * table_size_ + j) * 6 + k]; + } + const double& bboxes(int i, int j, int k) const + { + return bboxes_[(i * table_size_ + j) * 6 + k]; + } + + bool has_bboxes() const { return bboxes_ != nullptr; } + bool table_full() const { return table_full_; } private: int32_t* materials_; //!< material index (bins, table_size) double* volumes_; //!< volume in [cm^3] (bins, table_size) + double* bboxes_; //!< bounding boxes (bins, table_size, 6) int table_size_; //!< Size of hash table for each mesh element bool table_full_ {false}; //!< Whether the hash table is full }; @@ -239,6 +259,19 @@ public: void material_volumes(int nx, int ny, int nz, int max_materials, int32_t* materials, double* volumes) const; + //! Determine volume and bounding boxes of materials within each mesh element + // + //! \param[in] nx Number of samples in x direction + //! \param[in] ny Number of samples in y direction + //! \param[in] nz Number of samples in z direction + //! \param[in] max_materials Maximum number of materials in a single mesh + //! element + //! \param[inout] materials Array storing material indices + //! \param[inout] volumes Array storing volumes + //! \param[inout] bboxes Array storing bounding boxes (n_elems, table_size, 6) + void material_volumes(int nx, int ny, int nz, int max_materials, + int32_t* materials, double* volumes, double* bboxes) const; + //! Determine bounding box of mesh // //! \return Bounding box of mesh @@ -969,7 +1002,7 @@ public: Position sample_element(int32_t bin, uint64_t* seed) const override; - int get_bin(Position r) const override; + virtual int get_bin(Position r) const override; int n_bins() const override; @@ -1007,16 +1040,21 @@ public: protected: // Methods - //! Translate a bin value to an element reference virtual const libMesh::Elem& get_element_from_bin(int bin) const; //! Translate an element pointer to a bin index virtual int get_bin_from_element(const libMesh::Elem* elem) const; + // Data members libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set //!< during intialization + vector> + pl_; //!< per-thread point locators + libMesh::BoundingBox bbox_; //!< bounding box of the mesh + private: + // Methods void initialize() override; void set_mesh_pointer_from_filename(const std::string& filename); void build_eqn_sys(); @@ -1025,8 +1063,6 @@ private: unique_ptr unique_m_ = nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is //!< created inside OpenMC - vector> - pl_; //!< per-thread point locators unique_ptr equation_systems_; //!< pointer to the libMesh EquationSystems //!< instance @@ -1035,7 +1071,6 @@ private: std::unordered_map variable_map_; //!< mapping of variable names (tally scores) to libMesh //!< variable numbers - libMesh::BoundingBox bbox_; //!< bounding box of the mesh libMesh::dof_id_type first_element_id_; //!< id of the first element in the mesh }; @@ -1043,8 +1078,9 @@ private: class AdaptiveLibMesh : public LibMesh { public: // Constructor - AdaptiveLibMesh( - libMesh::MeshBase& input_mesh, double length_multiplier = 1.0); + AdaptiveLibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0, + const std::set& block_ids = + std::set()); // Overridden methods int n_bins() const override; @@ -1056,6 +1092,8 @@ public: void write(const std::string& filename) const override; + int get_bin(Position r) const override; + protected: // Overridden methods int get_bin_from_element(const libMesh::Elem* elem) const override; @@ -1064,6 +1102,9 @@ protected: private: // Data members + const std::set + block_ids_; //!< subdomains of the mesh to tally on + const bool block_restrict_; //!< whether a subset of the mesh is being used const libMesh::dof_id_type num_active_; //!< cached number of active elements std::vector diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 60b88a153b..58d8339396 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -163,7 +163,7 @@ bool multipole_in_range(const Nuclide& nuc, double E); namespace data { // Minimum/maximum transport energy for each particle type. Order corresponds to -// that of the ParticleType enum +// transport_index() for supported transport particles. extern array energy_min; extern array energy_max; diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 0f37719b94..2f6e6196bf 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -126,10 +126,6 @@ public: //! Functions //============================================================================ -std::string particle_type_to_str(ParticleType type); - -ParticleType str_to_particle_type(std::string str); - void add_surf_source_to_bank(Particle& p, const Surface& surf); } // namespace openmc diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index fdacfa765b..75166b36af 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -3,6 +3,7 @@ #include "openmc/array.h" #include "openmc/constants.h" +#include "openmc/particle_type.h" #include "openmc/position.h" #include "openmc/random_lcg.h" #include "openmc/tallies/filter_match.h" @@ -30,9 +31,6 @@ constexpr double CACHE_INVALID {-1.0}; //========================================================================== // Aliases and type definitions -//! Particle types -enum class ParticleType { neutron, photon, electron, positron }; - //! Saved ("banked") state of a particle //! NOTE: This structure's MPI type is built in initialize_mpi() of //! initialize.cpp. Any changes made to the struct here must also be @@ -496,7 +494,7 @@ private: MacroXS macro_xs_; CacheDataMG mg_xs_cache_; - ParticleType type_ {ParticleType::neutron}; + ParticleType type_; double E_; double E_last_; diff --git a/include/openmc/particle_type.h b/include/openmc/particle_type.h new file mode 100644 index 0000000000..fed3c92846 --- /dev/null +++ b/include/openmc/particle_type.h @@ -0,0 +1,172 @@ +//============================================================================== +// ParticleType class definition +//============================================================================== + +#ifndef OPENMC_PARTICLE_TYPE_H +#define OPENMC_PARTICLE_TYPE_H + +#include +#include +#include +#include +#include + +#include "openmc/constants.h" + +namespace openmc { + +//------------------------------------------------------------------------------ +// PDG constants (canonical particle identity as simple integers) +//------------------------------------------------------------------------------ + +inline constexpr int32_t PDG_NEUTRON = 2112; +inline constexpr int32_t PDG_PHOTON = 22; +inline constexpr int32_t PDG_ELECTRON = 11; +inline constexpr int32_t PDG_POSITRON = -11; +inline constexpr int32_t PDG_PROTON = 2212; +inline constexpr int32_t PDG_DEUTERON = 1000010020; +inline constexpr int32_t PDG_TRITON = 1000010030; +inline constexpr int32_t PDG_ALPHA = 1000020040; + +//------------------------------------------------------------------------------ +// ParticleType class (standard-layout, trivially copyable) +//------------------------------------------------------------------------------ + +class ParticleType { +public: + //---------------------------------------------------------------------------- + // Constructors + + // Default constructor: defaults to neutron + constexpr ParticleType() : pdg_number_(PDG_NEUTRON) {} + + // Constructor from PDG number + constexpr explicit ParticleType(int32_t pdg_number) : pdg_number_(pdg_number) + {} + + // Constructor from particle name string (e.g., "neutron", "photon", "Fe56") + explicit ParticleType(std::string_view str); + + // Constructor from Z, A, and metastable state for nuclear particles + constexpr ParticleType(int Z, int A, int m = 0) + : pdg_number_(1000000000 + Z * 10000 + A * 10 + m) + {} + + //---------------------------------------------------------------------------- + // Accessors + + // Accessor for the underlying PDG number + constexpr int32_t pdg_number() const { return pdg_number_; } + + //---------------------------------------------------------------------------- + // Methods + + // Convert to string representation + std::string str() const; + + // Check if this represents a nucleus (vs elementary particle) + constexpr bool is_nucleus() const + { + // PDG nuclear codes are >= 1000000000 (100ZZZAAAI format) + return pdg_number_ >= 1000000000; + } + + // Get transport index (0-3 for transportable particles, C_NONE otherwise) + constexpr int transport_index() const; + + // Check if this is a neutron + constexpr bool is_neutron() const { return pdg_number_ == PDG_NEUTRON; } + + // Check if this is a photon + constexpr bool is_photon() const { return pdg_number_ == PDG_PHOTON; } + + constexpr bool is_transportable() const + { + return this->transport_index() != C_NONE; + } + + //---------------------------------------------------------------------------- + // Static factory methods + + static constexpr ParticleType neutron() { return ParticleType {PDG_NEUTRON}; } + static constexpr ParticleType photon() { return ParticleType {PDG_PHOTON}; } + static constexpr ParticleType electron() + { + return ParticleType {PDG_ELECTRON}; + } + static constexpr ParticleType positron() + { + return ParticleType {PDG_POSITRON}; + } + static constexpr ParticleType proton() { return ParticleType {PDG_PROTON}; } + static constexpr ParticleType deuteron() + { + return ParticleType {PDG_DEUTERON}; + } + static constexpr ParticleType triton() { return ParticleType {PDG_TRITON}; } + static constexpr ParticleType alpha() { return ParticleType {PDG_ALPHA}; } + +private: + int32_t pdg_number_; +}; + +//------------------------------------------------------------------------------ +// Static assertions to ensure standard-layout and trivially copyable +//------------------------------------------------------------------------------ + +static_assert(std::is_standard_layout_v, + "ParticleType must be standard-layout"); +static_assert(std::is_trivially_copyable_v, + "ParticleType must be trivially copyable"); +static_assert(sizeof(ParticleType) == sizeof(int32_t), + "ParticleType must be same size as int32_t"); + +//------------------------------------------------------------------------------ +// Comparison operators (free functions for symmetry) +//------------------------------------------------------------------------------ + +constexpr bool operator==(ParticleType lhs, ParticleType rhs) +{ + return lhs.pdg_number() == rhs.pdg_number(); +} + +constexpr bool operator!=(ParticleType lhs, ParticleType rhs) +{ + return lhs.pdg_number() != rhs.pdg_number(); +} + +constexpr bool operator<(ParticleType lhs, ParticleType rhs) +{ + return lhs.pdg_number() < rhs.pdg_number(); +} + +//------------------------------------------------------------------------------ +// ParticleType member function implementations (inline) +//------------------------------------------------------------------------------ + +constexpr int ParticleType::transport_index() const +{ + switch (pdg_number_) { + case PDG_NEUTRON: + return 0; + case PDG_PHOTON: + return 1; + case PDG_ELECTRON: + return 2; + case PDG_POSITRON: + return 3; + default: + return C_NONE; + } +} + +//------------------------------------------------------------------------------ +// Legacy conversion helpers +//------------------------------------------------------------------------------ + +// Legacy enum code (0..3) to ParticleType conversion +ParticleType legacy_particle_index_to_type(int code); + +} // namespace openmc + +#endif // OPENMC_PARTICLE_TYPE_H diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 4df4e5d8d3..0d4086966d 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -107,6 +107,7 @@ public: vector nu_sigma_f_; vector sigma_f_; vector chi_; + vector kappa_fission_; // 3D arrays stored in 1D representing values for all materials x energy // groups x energy groups diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index 3dec48bf26..3db6510697 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -19,9 +19,9 @@ public: //---------------------------------------------------------------------------- // Methods - void compute_segment_correction_factors(); void apply_fixed_sources_and_mesh_domains(); void prepare_fixed_sources_adjoint(); + void prepare_adjoint_simulation(); void simulate(); void output_simulation_results() const; void instability_check( @@ -34,9 +34,15 @@ public: // Accessors FlatSourceDomain* domain() const { return domain_.get(); } + //---------------------------------------------------------------------------- + // Public data members + + // Flag for adjoint simulation; + bool adjoint_needed_; + private: //---------------------------------------------------------------------------- - // Data members + // Private data members // Contains all flat source region data unique_ptr domain_; @@ -51,6 +57,9 @@ private: // Number of energy groups int negroups_; + // Toggle for first simulation + bool is_first_simulation_; + }; // class RandomRaySimulation //============================================================================ @@ -60,6 +69,7 @@ private: void openmc_run_random_ray(); void validate_random_ray_inputs(); void openmc_reset_random_ray(); +void print_adjoint_header(); } // namespace openmc diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 0f5a747fff..c20d46abe9 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -146,6 +146,7 @@ public: // Scalar fields int* material_; + double* density_mult_; int* is_small_; int* n_hits_; int* birthday_; @@ -195,6 +196,9 @@ public: int& material() { return *material_; } const int material() const { return *material_; } + double& density_mult() { return *density_mult_; } + const double density_mult() const { return *density_mult_; } + int& is_small() { return *is_small_; } const int is_small() const { return *is_small_; } @@ -316,7 +320,9 @@ public: //--------------------------------------- // Scalar fields - int material_ {0}; //!< Index in openmc::model::materials array + int material_ {0}; //!< Index in openmc::model::materials array + double density_mult_ {1.0}; //!< A density multiplier queried from the cell + //!< corresponding to the source region. OpenMPMutex lock_; double volume_ { 0.0}; //!< Volume (computed from the sum of ray crossing lengths) @@ -394,6 +400,9 @@ public: int& material(int64_t sr) { return material_[sr]; } const int material(int64_t sr) const { return material_[sr]; } + double& density_mult(int64_t sr) { return density_mult_[sr]; } + const double density_mult(int64_t sr) const { return density_mult_[sr]; } + int& is_small(int64_t sr) { return is_small_[sr]; } const int is_small(int64_t sr) const { return is_small_[sr]; } @@ -625,6 +634,7 @@ private: // SoA storage for scalar fields (one item per source region) vector material_; + vector density_mult_; vector is_small_; vector n_hits_; vector mesh_; diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index 4fbbc1b626..9a8eab7d91 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -10,7 +10,7 @@ #include "openmc/chain.h" #include "openmc/endf.h" #include "openmc/memory.h" // for unique_ptr -#include "openmc/particle.h" +#include "openmc/particle_type.h" #include "openmc/vector.h" // for vector namespace openmc { diff --git a/include/openmc/source.h b/include/openmc/source.h index 36c821fc07..2f32aa2a05 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -12,7 +12,7 @@ #include "openmc/distribution_multi.h" #include "openmc/distribution_spatial.h" #include "openmc/memory.h" -#include "openmc/particle.h" +#include "openmc/particle_type.h" #include "openmc/vector.h" namespace openmc { @@ -148,11 +148,11 @@ protected: private: // Data members - ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted - UPtrSpace space_; //!< Spatial distribution - UPtrAngle angle_; //!< Angular distribution - UPtrDist energy_; //!< Energy distribution - UPtrDist time_; //!< Time distribution + ParticleType particle_; //!< Type of particle emitted + UPtrSpace space_; //!< Spatial distribution + UPtrAngle angle_; //!< Angular distribution + UPtrDist energy_; //!< Energy distribution + UPtrDist time_; //!< Time distribution }; //============================================================================== diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h index 863a6d282f..9932a6037e 100644 --- a/include/openmc/tallies/filter_particle.h +++ b/include/openmc/tallies/filter_particle.h @@ -1,7 +1,7 @@ #ifndef OPENMC_TALLIES_FILTER_PARTICLE_H #define OPENMC_TALLIES_FILTER_PARTICLE_H -#include "openmc/particle.h" +#include "openmc/particle_type.h" #include "openmc/span.h" #include "openmc/tallies/filter.h" #include "openmc/vector.h" diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 7638155228..0d7435ca2a 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -10,7 +10,7 @@ #include "openmc/constants.h" #include "openmc/memory.h" #include "openmc/mesh.h" -#include "openmc/particle.h" +#include "openmc/particle_type.h" #include "openmc/span.h" #include "openmc/tallies/tally.h" #include "openmc/vector.h" @@ -193,10 +193,9 @@ public: private: //---------------------------------------------------------------------------- // Data members - int32_t id_; //!< Unique ID - int64_t index_; //!< Index into weight windows vector - ParticleType particle_type_ { - ParticleType::neutron}; //!< Particle type to apply weight windows to + int32_t id_; //!< Unique ID + int64_t index_; //!< Index into weight windows vector + ParticleType particle_type_; //!< Particle type to apply weight windows to vector energy_bounds_; //!< Energy boundaries [eV] xt::xtensor lower_ww_; //!< Lower weight window bounds (shape: //!< energy_bins, mesh_bins (k, j, i)) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 25ded24cbe..bc21b2e56a 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -203,7 +203,9 @@ class AtomicRelaxation(EqualityMixin): for subshell, df in transitions.items(): cv.check_value('subshell', subshell, _SUBSHELLS) cv.check_type('transitions', df, pd.DataFrame) - self._transitions = transitions + self._transitions = { + subshell: df.convert_dtypes() for subshell, df in transitions.items() + } @classmethod def from_ace(cls, ace): diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ee742f22bc..056f7c2737 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -209,6 +209,8 @@ class TransportOperator(ABC): simulation. full_burn_list : list of int All burnable materials in the geometry. + name_list : list of str + Material names corresponding to materials in burn_list """ def finalize(self): diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 2d7694093b..7928e89ee3 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -211,7 +211,7 @@ class OpenMCOperator(TransportOperator): "section data.") warn(msg) if mat.depletable: - burnable_mats.add(str(mat.id)) + burnable_mats.add((str(mat.id), mat.name)) if mat.volume is None: if mat.name is None: msg = ("Volume not specified for depletable material " @@ -229,9 +229,12 @@ class OpenMCOperator(TransportOperator): "No depletable materials were found in the model.") # Sort the sets - burnable_mats = sorted(burnable_mats, key=int) + burnable_mats = sorted(burnable_mats, key=lambda x: int(x[0])) model_nuclides = sorted(model_nuclides) + # Store material names for later use + burnable_mats, self.name_list = zip(*burnable_mats) + # Construct a global nuclide dictionary, burned first nuclides = list(self.chain.nuclide_dict) for nuc in model_nuclides: @@ -541,6 +544,8 @@ class OpenMCOperator(TransportOperator): A list of all material IDs to be burned. Used for sorting the simulation. full_burn_list : list List of all burnable material IDs + name_list : list of str + Material names corresponding to materials in burn_list """ nuc_list = self.number.burnable_nuclides @@ -554,4 +559,4 @@ class OpenMCOperator(TransportOperator): volume_list = comm.allgather(volume) volume = {k: v for d in volume_list for k, v in d.items()} - return volume, nuc_list, burn_list, self.burnable_mats + return volume, nuc_list, burn_list, self.burnable_mats, self.name_list diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 97277cbda8..57bbe437ff 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -269,6 +269,7 @@ class R2SManager: # Compute material volume fractions on the mesh if mat_vol_kwargs is None: mat_vol_kwargs = {} + mat_vol_kwargs.setdefault('bounding_boxes', True) self.results['mesh_material_volumes'] = mmv = comm.bcast( self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs)) @@ -539,14 +540,15 @@ class R2SManager: def get_decay_photon_source_mesh( self, time_index: int = -1 - ) -> list[openmc.MeshSource]: + ) -> list[openmc.IndependentSource]: """Create decay photon source for a mesh-based calculation. - This function creates N :class:`MeshSource` objects where N is the - maximum number of unique materials that appears in a single mesh - element. For each mesh element-material combination, and - IndependentSource instance is created with a spatial constraint limited - the sampled decay photons to the correct region. + For each mesh element-material combination, an + :class:`~openmc.IndependentSource` is created with a + :class:`~openmc.stats.Box` spatial distribution based on the bounding + box of the material within the mesh element. A material constraint is + also applied so that sampled source sites are limited to the correct + region. When the photon transport model is different from the neutron model, the photon MeshMaterialVolumes is used to determine whether an (element, @@ -559,19 +561,15 @@ class R2SManager: Returns ------- - list of openmc.MeshSource - A list of MeshSource objects, each containing IndependentSource - instances for the decay photons in the corresponding mesh element. + list of openmc.IndependentSource + A list of IndependentSource objects for the decay photons, one for + each mesh element-material combination with non-zero source strength. """ mat_dict = self.neutron_model._get_all_materials() - # Some MeshSource objects will have empty positions; create a "null source" - # that is used for this case - null_source = openmc.IndependentSource(particle='photon', strength=0.0) - - # List to hold sources for each MeshSource (length = N) - source_lists = [] + # List to hold all sources + sources = [] # Index in the overall list of activated materials index_mat = 0 @@ -594,7 +592,7 @@ class R2SManager: if mat_id is not None } - for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): + for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True): # Skip void volume if mat_id is None: continue @@ -604,30 +602,27 @@ class R2SManager: index_mat += 1 continue - # Check whether a new MeshSource object is needed - if j >= len(source_lists): - source_lists.append([null_source]*n_elements) - # Get activated material composition original_mat = materials[index_mat] activated_mat = results[time_index].get_material(str(original_mat.id)) - # Create decay photon source source + # Create decay photon source energy = activated_mat.get_decay_photon_energy() if energy is not None: strength = energy.integral() - source_lists[j][index_elem] = openmc.IndependentSource( + space = openmc.stats.Box(*bbox) + sources.append(openmc.IndependentSource( + space=space, energy=energy, particle='photon', strength=strength, constraints={'domains': [mat_dict[mat_id]]} - ) + )) # Increment index of activated material index_mat += 1 - # Return list of mesh sources - return [openmc.MeshSource(self.domains, sources) for sources in source_lists] + return sources def load_results(self, path: PathLike): """Load results from a previous R2S calculation. diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py index ff39e9acb6..27420246f6 100644 --- a/openmc/deplete/stepresult.py +++ b/openmc/deplete/stepresult.py @@ -70,6 +70,7 @@ class StepResult: self.index_mat = None self.index_nuc = None self.mat_to_hdf5_ind = None + self.name_list = None self.data = None @@ -138,7 +139,7 @@ class StepResult: def n_hdf5_mats(self): return len(self.mat_to_hdf5_ind) - def allocate(self, volume, nuc_list, burn_list, full_burn_list): + def allocate(self, volume, nuc_list, burn_list, full_burn_list, name_list=None): """Allocate memory for depletion step data Parameters @@ -151,12 +152,15 @@ class StepResult: A list of all mat IDs to be burned. Used for sorting the simulation. full_burn_list : list of str List of all burnable material IDs + name_list : list of str, optional + Material names corresponding to materials in burn_list """ self.volume = copy.deepcopy(volume) self.index_nuc = {nuc: i for i, nuc in enumerate(nuc_list)} self.index_mat = {mat: i for i, mat in enumerate(burn_list)} self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)} + self.mat_to_name = dict(zip(burn_list, name_list)) if name_list is not None else {} # Create storage array self.data = np.zeros((self.n_mat, self.n_nuc)) @@ -184,7 +188,7 @@ class StepResult: # Direct transfer direct_attrs = ("time", "k", "source_rate", "index_nuc", - "mat_to_hdf5_ind", "proc_time") + "mat_to_hdf5_ind", "mat_to_name", "proc_time") for attr in direct_attrs: setattr(new, attr, getattr(self, attr)) # Get applicable slice of data @@ -223,6 +227,8 @@ class StepResult: f'mat_id {mat_id} not found in StepResult. Available mat_id ' f'values are {list(self.volume.keys())}' ) from e + if mat_id in self.mat_to_name: + material.name = self.mat_to_name[mat_id] for nuc, _ in sorted(self.index_nuc.items(), key=lambda x: x[1]): atoms = self[mat_id, nuc] if atoms <= 0.0: @@ -313,6 +319,8 @@ class StepResult: mat_single_group = mat_group.create_group(mat) mat_single_group.attrs["index"] = self.mat_to_hdf5_ind[mat] mat_single_group.attrs["volume"] = self.volume[mat] + if mat in self.mat_to_name: + mat_single_group.attrs["name"] = self.mat_to_name[mat] nuc_group = handle.create_group("nuclides") @@ -495,6 +503,7 @@ class StepResult: results.volume = {} results.index_mat = {} results.index_nuc = {} + results.mat_to_name = {} rxn_nuc_to_ind = {} rxn_to_ind = {} @@ -504,6 +513,8 @@ class StepResult: results.volume[mat] = vol results.index_mat[mat] = ind + if "name" in mat_handle.attrs: + results.mat_to_name[mat] = mat_handle.attrs["name"] for nuc, nuc_handle in handle["/nuclides"].items(): ind_atom = nuc_handle.attrs["atom number index"] @@ -569,11 +580,11 @@ class StepResult: .. versionadded:: 0.14.0 """ # Get indexing terms - vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info() + vol_dict, nuc_list, burn_list, full_burn_list, name_list = op.get_results_info() # Create results results = StepResult() - results.allocate(vol_dict, nuc_list, burn_list, full_burn_list) + results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, name_list) n_mat = len(burn_list) diff --git a/openmc/element.py b/openmc/element.py index f1d8f5f24d..5deede9710 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -184,8 +184,6 @@ class Element(str): for nuclide in absent_nuclides: if nuclide in ['O17', 'O18'] and 'O16' in mutual_nuclides: abundances['O16'] += NATURAL_ABUNDANCE[nuclide] - elif nuclide == 'Ta180_m1' and 'Ta180' in library_nuclides: - abundances['Ta180'] = NATURAL_ABUNDANCE[nuclide] elif nuclide == 'Ta180_m1' and 'Ta181' in mutual_nuclides: abundances['Ta181'] += NATURAL_ABUNDANCE[nuclide] elif nuclide == 'W180' and 'W182' in mutual_nuclides: diff --git a/openmc/examples.py b/openmc/examples.py index 01dd9d01f9..350a4d24d5 100644 --- a/openmc/examples.py +++ b/openmc/examples.py @@ -4,7 +4,7 @@ import numpy as np import openmc - +PINCELL_PITCH = 1.26 # cm def pwr_pin_cell() -> openmc.Model: """Create a PWR pin-cell model. @@ -51,7 +51,7 @@ def pwr_pin_cell() -> openmc.Model: model.materials = (fuel, clad, hot_water) # Instantiate ZCylinder surfaces - pitch = 1.26 + pitch = PINCELL_PITCH fuel_or = openmc.ZCylinder(x0=0, y0=0, r=0.39218, name='Fuel OR') clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR') left = openmc.XPlane(x0=-pitch/2, name='left', boundary_type='reflective') @@ -319,14 +319,14 @@ def pwr_core() -> openmc.Model: l100 = openmc.RectLattice( name='Fuel assembly (lower half)', lattice_id=100) l100.lower_left = (-10.71, -10.71) - l100.pitch = (1.26, 1.26) + l100.pitch = (PINCELL_PITCH, PINCELL_PITCH) l100.universes = np.tile(fuel_cold, (17, 17)) l100.universes[tube_x, tube_y] = tube_cold l101 = openmc.RectLattice( name='Fuel assembly (upper half)', lattice_id=101) l101.lower_left = (-10.71, -10.71) - l101.pitch = (1.26, 1.26) + l101.pitch = (PINCELL_PITCH, PINCELL_PITCH) l101.universes = np.tile(fuel_hot, (17, 17)) l101.universes[tube_x, tube_y] = tube_hot @@ -350,7 +350,7 @@ def pwr_core() -> openmc.Model: # Define core lattices l200 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=200) l200.lower_left = (-224.91, -224.91) - l200.pitch = (21.42, 21.42) + l200.pitch = (17 * PINCELL_PITCH, 17 * PINCELL_PITCH) l200.universes = [ [fa_cw]*21, [fa_cw]*21, @@ -376,7 +376,7 @@ def pwr_core() -> openmc.Model: l201 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=201) l201.lower_left = (-224.91, -224.91) - l201.pitch = (21.42, 21.42) + l201.pitch = (17 * PINCELL_PITCH, 17 * PINCELL_PITCH) l201.universes = [ [fa_hw]*21, [fa_hw]*21, @@ -488,7 +488,7 @@ def pwr_assembly() -> openmc.Model: clad_or = openmc.ZCylinder(x0=0, y0=0, r=0.45720, name='Clad OR') # Create boundary planes to surround the geometry - pitch = 21.42 + pitch = 17 * PINCELL_PITCH min_x = openmc.XPlane(x0=-pitch/2, boundary_type='reflective') max_x = openmc.XPlane(x0=+pitch/2, boundary_type='reflective') min_y = openmc.YPlane(y0=-pitch/2, boundary_type='reflective') @@ -514,7 +514,7 @@ def pwr_assembly() -> openmc.Model: # Create fuel assembly Lattice assembly = openmc.RectLattice(name='Fuel Assembly') - assembly.pitch = (pitch/17, pitch/17) + assembly.pitch = (PINCELL_PITCH, PINCELL_PITCH) assembly.lower_left = (-pitch/2, -pitch/2) # Create array indices for guide tube locations in lattice @@ -656,24 +656,21 @@ def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5') -> openmc.Model return model -def random_ray_lattice() -> openmc.Model: - """Create a 2x2 PWR pincell asymmetrical lattic eexample. - - This model is a 2x2 reflective lattice of fuel pins with one of the lattice - locations having just moderator instead of a fuel pin. It uses 7 group - cross section data. +def _generate_c5g7_materials() -> openmc.Materials: + """Generate materials utilizing multi-group cross sections based on the + the C5G7 Benchmark. Returns ------- - model : openmc.Model - A PWR 2x2 lattice model + materials : openmc.Materials + Materials object containing UO2 and water materials. + Data Sources + ------------ + All cross section data are from: + Lewis et al., "Benchmark specification for determinisitc 2D/3D MOX fuel + assembly transport calculations without spatial homogenization" """ - model = openmc.Model() - - ########################################################################### - # Create MGXS data for the problem - # Instantiate the energy group data group_edges = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6] groups = openmc.mgxs.EnergyGroups(group_edges) @@ -704,9 +701,10 @@ def random_ray_lattice() -> openmc.Model: uo2_xsdata.set_fission([7.21206e-03, 8.19301e-04, 6.45320e-03, 1.85648e-02, 1.78084e-02, 8.30348e-02, 2.16004e-01]) - uo2_xsdata.set_nu_fission([2.005998e-02, 2.027303e-03, 1.570599e-02, - 4.518301e-02, 4.334208e-02, 2.020901e-01, - 5.257105e-01]) + nu_fission = np.array([2.005998e-02, 2.027303e-03, 1.570599e-02, + 4.518301e-02, 4.334208e-02, 2.020901e-01, + 5.257105e-01]) + uo2_xsdata.set_nu_fission(nu_fission) uo2_xsdata.set_chi([5.8791e-01, 4.1176e-01, 3.3906e-04, 1.1761e-07, 0.0000e+00, 0.0000e+00, 0.0000e+00]) @@ -752,14 +750,28 @@ def random_ray_lattice() -> openmc.Model: # Instantiate a Materials collection and export to XML materials = openmc.Materials([uo2, water]) materials.cross_sections = "mgxs.h5" + return materials - ########################################################################### - # Define problem geometry +def _generate_subdivided_pin_cell(uo2, water) -> openmc.Universe: + """Create a radially and azimuthally subdivided pin cell universe. Helper + function for random_ray_pin_cell() and random_ray_lattice() + + Parameters + ---------- + uo2 : openmc.Material + UO2 material + water : openmc.Material + Water material + + Returns + ------- + pincell : openmc.Universe + Universe containing an unbounded pin cell + + """ ######################################## - # Define an unbounded pincell universe - - pitch = 1.26 + # Define an unbounded pin cell universe # Create a surface for the fuel outer radius fuel_or = openmc.ZCylinder(r=0.54, name='Fuel OR') @@ -781,7 +793,7 @@ def random_ray_lattice() -> openmc.Model: moderator_c = openmc.Cell( fill=water, region=+outer_ring_b, name='moderator outer c') - # Create pincell universe + # Create pin cell universe pincell_base = openmc.Universe() # Register Cells with Universe @@ -801,19 +813,125 @@ def random_ray_lattice() -> openmc.Model: for i in range(8): azimuthal_cell = openmc.Cell(name=f'azimuthal_cell_{i}') azimuthal_cell.fill = pincell_base - azimuthal_cell.region = +azimuthal_planes[i] & -azimuthal_planes[(i+1) % 8] + azimuthal_cell.region = + \ + azimuthal_planes[i] & -azimuthal_planes[(i+1) % 8] azimuthal_cells.append(azimuthal_cell) # Create a geometry with the azimuthal universes pincell = openmc.Universe(cells=azimuthal_cells, name='pincell') + return pincell + + +def random_ray_pin_cell() -> openmc.Model: + """Create a PWR pin cell example using C5G7 cross section data. + cross section data. + + Returns + ------- + model : openmc.Model + A PWR pin cell model + + """ + model = openmc.Model() + + ########################################################################### + # Create Materials for the problem + materials = _generate_c5g7_materials() + uo2 = materials[0] + water = materials[1] + + ########################################################################### + # Define problem geometry + pincell = _generate_subdivided_pin_cell(uo2, water) + + ######################################## + # Define cell containing lattice and other stuff + pitch = PINCELL_PITCH + box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective') + + pincell = openmc.Cell(fill=pincell, region=-box, name='pincell') + + # Create a geometry with the top-level cell + geometry = openmc.Geometry([pincell]) + + ########################################################################### + # Define problem settings + + # Instantiate a Settings object, set all runtime parameters, and export to XML + settings = openmc.Settings() + settings.energy_mode = "multi-group" + settings.batches = 400 + settings.inactive = 200 + settings.particles = 100 + + # Create an initial uniform spatial source distribution over fissionable zones + lower_left = (-pitch / 2, -pitch / 2, -1) + upper_right = (pitch / 2, pitch / 2, 1) + uniform_dist = openmc.stats.Box(lower_left, upper_right) + rr_source = openmc.IndependentSource(space=uniform_dist) + + settings.random_ray['distance_active'] = 100.0 + settings.random_ray['distance_inactive'] = 20.0 + settings.random_ray['ray_source'] = rr_source + settings.random_ray['volume_normalized_flux_tallies'] = True + + ########################################################################### + # Define tallies + # Now use the mesh filter in a tally and indicate what scores are desired + tally = openmc.Tally(name="Pin tally") + tally.scores = ['flux', 'fission', 'nu-fission'] + tally.estimator = 'analog' + + # Instantiate a Tallies collection and export to XML + tallies = openmc.Tallies([tally]) + + ########################################################################### + # Exporting to OpenMC model + ########################################################################### + + model.geometry = geometry + model.materials = materials + model.settings = settings + model.tallies = tallies + return model + + +def random_ray_lattice() -> openmc.Model: + """Create a 2x2 PWR pin cell asymmetrical lattice example. + + This model is a 2x2 reflective lattice of fuel pins with one of the lattice + locations having just moderator instead of a fuel pin. It uses C5G7 + cross section data. + + Returns + ------- + model : openmc.Model + A PWR 2x2 lattice model + + """ + model = openmc.Model() + + ########################################################################### + # Create Materials for the problem + materials = _generate_c5g7_materials() + uo2 = materials[0] + water = materials[1] + + ########################################################################### + # Define problem geometry + pincell = _generate_subdivided_pin_cell(uo2, water) + ######################################## # Define a moderator lattice universe - moderator_infinite = openmc.Cell(fill=water, name='moderator infinite') + moderator_infinite = openmc.Cell(name='moderator infinite') + moderator_infinite.fill = water + mu = openmc.Universe() mu.add_cells([moderator_infinite]) + pitch = PINCELL_PITCH lattice = openmc.RectLattice() lattice.lower_left = [-pitch/2.0, -pitch/2.0] lattice.pitch = [pitch/10.0, pitch/10.0] @@ -837,8 +955,7 @@ def random_ray_lattice() -> openmc.Model: ######################################## # Define cell containing lattice and other stuff - box = openmc.model.RectangularPrism( - pitch*2, pitch*2, boundary_type='reflective') + box = openmc.model.RectangularPrism(pitch*2, pitch*2, boundary_type='reflective') assembly = openmc.Cell(fill=lattice2x2, region=-box, name='assembly') diff --git a/openmc/filter.py b/openmc/filter.py index e7c54f2f0e..53ec93e21d 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -35,7 +35,6 @@ _CURRENT_NAMES = ( 'z-min out', 'z-min in', 'z-max out', 'z-max in' ) -_PARTICLES = {'neutron', 'photon', 'electron', 'positron'} class FilterMeta(ABCMeta): @@ -735,9 +734,8 @@ class ParticleFilter(Filter): Parameters ---------- - bins : str, or sequence of str - The particles to tally represented as strings ('neutron', 'photon', - 'electron', 'positron'). + bins : str, int, openmc.ParticleType, or sequence + The particle types to tally represented as names, PDG numbers, or types. filter_id : int Unique identifier for the filter @@ -763,11 +761,16 @@ class ParticleFilter(Filter): @Filter.bins.setter def bins(self, bins): - cv.check_type('bins', bins, Sequence, str) + if isinstance(bins, (str, Integral, openmc.ParticleType)): + bins = [bins] + else: + cv.check_type('bins', bins, Sequence, + (str, Integral, openmc.ParticleType)) bins = np.atleast_1d(bins) - for edge in bins: - cv.check_value('filter bin', edge, _PARTICLES) - self._bins = bins + normalized = [] + for entry in bins: + normalized.append(str(openmc.ParticleType(entry))) + self._bins = np.array(normalized, dtype=str) @classmethod def from_hdf5(cls, group, **kwargs): diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 02c7784d1c..79310bb37f 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -28,7 +28,7 @@ class _SourceSite(Structure): ('wgt', c_double), ('delayed_group', c_int), ('surf_id', c_int), - ('particle', c_int), + ('particle', c_int32), ('parent_nuclide', c_int), ('parent_id', c_int64), ('progeny_id', c_int64)] diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 5cf496107b..dc81c3487b 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -132,6 +132,9 @@ _dll.openmc_meshsurface_filter_set_translation.errcheck = _error_handler _dll.openmc_new_filter.argtypes = [c_char_p, POINTER(c_int32)] _dll.openmc_new_filter.restype = c_int _dll.openmc_new_filter.errcheck = _error_handler +_dll.openmc_particle_filter_get_bins.argtypes = [c_int32, POINTER(c_int32)] +_dll.openmc_particle_filter_get_bins.restype = c_int +_dll.openmc_particle_filter_get_bins.errcheck = _error_handler _dll.openmc_spatial_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)] _dll.openmc_spatial_legendre_filter_get_order.restype = c_int _dll.openmc_spatial_legendre_filter_get_order.errcheck = _error_handler @@ -402,8 +405,8 @@ class MeshFilter(Filter): translation : Iterable of float 3-D coordinates of the translation vector rotation : Iterable of float - The rotation matrix or angles of the filter mesh. This can either be - a fully specified 3 x 3 rotation matrix or an Iterable of length 3 + The rotation matrix or angles of the filter mesh. This can either be + a fully specified 3 x 3 rotation matrix or an Iterable of length 3 with the angles in degrees about the x, y, and z axes, respectively. """ @@ -454,7 +457,7 @@ class MeshFilter(Filter): else: raise ValueError( f'Invalid size of rotation matrix: {rot_size}') - + @rotation.setter def rotation(self, rotation_data): flat_rotation = np.asarray(rotation_data, dtype=float).flatten() @@ -598,9 +601,9 @@ class ParticleFilter(Filter): @property def bins(self): - particle_i = np.zeros((self.n_bins,), dtype=c_int) + particle_i = np.zeros((self.n_bins,), dtype=c_int32) _dll.openmc_particle_filter_get_bins( - self._index, particle_i.ctypes.data_as(POINTER(c_int))) + self._index, particle_i.ctypes.data_as(POINTER(c_int32))) return [ParticleType(i) for i in particle_i] diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 3a720f98a3..19e6f74d7a 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -1,5 +1,5 @@ from collections.abc import Mapping, Sequence -from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, +from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, c_void_p, create_string_buffer, c_size_t) from math import sqrt import sys @@ -47,7 +47,8 @@ _dll.openmc_mesh_bounding_box.argtypes = [ _dll.openmc_mesh_bounding_box.restype = c_int _dll.openmc_mesh_bounding_box.errcheck = _error_handler _dll.openmc_mesh_material_volumes.argtypes = [ - c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double] + c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double, + c_void_p] _dll.openmc_mesh_material_volumes.restype = c_int _dll.openmc_mesh_material_volumes.errcheck = _error_handler _dll.openmc_mesh_get_plot_bins.argtypes = [ @@ -188,6 +189,7 @@ class Mesh(_FortranObjectWithID): n_samples: int | tuple[int, int, int] = 10_000, max_materials: int = 4, output: bool = True, + bounding_boxes: bool = False, ) -> MeshMaterialVolumes: """Determine volume of materials in each mesh element. @@ -213,6 +215,11 @@ class Mesh(_FortranObjectWithID): Estimated maximum number of materials in any given mesh element. output : bool, optional Whether or not to show output. + bounding_boxes : bool, optional + Whether or not to compute an axis-aligned bounding box for each + (mesh element, material) combination. When enabled, the bounding + box encloses the ray-estimator prisms used for the volume + estimation. Returns ------- @@ -243,23 +250,36 @@ class Mesh(_FortranObjectWithID): table_size = slot_factor*max_materials materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32) volumes = np.zeros((n, table_size), dtype=np.float64) + bboxes = None + if bounding_boxes: + bboxes = np.empty((n, table_size, 6), dtype=np.float64) + bboxes[..., 0:3] = np.inf + bboxes[..., 3:6] = -np.inf # Run material volume calculation while True: try: + bboxes_ptr = None + if bboxes is not None: + bboxes_ptr = bboxes.ctypes.data_as(POINTER(c_double)) with quiet_dll(output): _dll.openmc_mesh_material_volumes( - self._index, nx, ny, nz, table_size, materials, volumes) + self._index, nx, ny, nz, table_size, materials, + volumes, bboxes_ptr) except AllocationError: # Increase size of result array and try again table_size *= 2 materials = np.full((n, table_size), EMPTY_SLOT, dtype=np.int32) volumes = np.zeros((n, table_size), dtype=np.float64) + if bounding_boxes: + bboxes = np.empty((n, table_size, 6), dtype=np.float64) + bboxes[..., 0:3] = np.inf + bboxes[..., 3:6] = -np.inf else: # If no error, break out of loop break - return MeshMaterialVolumes(materials, volumes) + return MeshMaterialVolumes(materials, volumes, bboxes) def get_plot_bins( self, diff --git a/openmc/lib/weight_windows.py b/openmc/lib/weight_windows.py index ed442d33ff..2b26d3b55f 100644 --- a/openmc/lib/weight_windows.py +++ b/openmc/lib/weight_windows.py @@ -53,11 +53,11 @@ _dll.openmc_weight_windows_get_energy_bounds.argtypes = [c_int32, POINTER(POINTE _dll.openmc_weight_windows_get_energy_bounds.restype = c_int _dll.openmc_weight_windows_get_energy_bounds.errcheck = _error_handler -_dll.openmc_weight_windows_set_particle.argtypes = [c_int32, c_int] +_dll.openmc_weight_windows_set_particle.argtypes = [c_int32, c_int32] _dll.openmc_weight_windows_set_particle.restype = c_int _dll.openmc_weight_windows_set_particle.errcheck = _error_handler -_dll.openmc_weight_windows_get_particle.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_weight_windows_get_particle.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_weight_windows_get_particle.restype = c_int _dll.openmc_weight_windows_get_particle.errcheck = _error_handler @@ -201,16 +201,13 @@ class WeightWindows(_FortranObjectWithID): @property def particle(self): - val = c_int() + val = c_int32() _dll.openmc_weight_windows_get_particle(self._index, val) return ParticleType(val.value) @particle.setter def particle(self, p): - if isinstance(p, str): - p = ParticleType.from_string(p) - else: - p = ParticleType(p) + p = ParticleType(p) _dll.openmc_weight_windows_set_particle(self._index, int(p)) @property @@ -304,10 +301,10 @@ class WeightWindows(_FortranObjectWithID): ---------- tally : openmc.lib.Tally The tally used to create the WeightWindows instance. - particle : openmc.ParticleType or str, optional + particle : openmc.ParticleType or str or int, optional The particle type to use for the WeightWindows instance. Should be - specified as an instance of ParticleType or as a string with a value of - 'neutron' or 'photon'. + specified as an instance of ParticleType, a PDG number, or as a + name. Returns ------- @@ -317,7 +314,8 @@ class WeightWindows(_FortranObjectWithID): Raises ------ ValueError - If the particle parameter is not an instance of ParticleType or a string. + If the particle parameter is not an instance of ParticleType, a string, + or an integer PDG number. ValueError If the particle parameter is not a valid particle type (i.e., not 'neutron' or 'photon'). @@ -328,12 +326,13 @@ class WeightWindows(_FortranObjectWithID): If the tally does not have a MeshFilter. """ # do some checks on particle value - if not isinstance(particle, (ParticleType, str)): - raise ValueError(f"Parameter 'particle' must be {ParticleType} or one of ('neutron', 'photon').") + if not isinstance(particle, (ParticleType, str, int)): + raise ValueError( + f"Parameter 'particle' must be {ParticleType} or one of ('neutron', 'photon')." + ) # convert particle type if needed - if isinstance(particle, str): - particle = ParticleType.from_string(particle) + particle = ParticleType(particle) if particle not in (ParticleType.NEUTRON, ParticleType.PHOTON): raise ValueError('Weight windows can only be applied for neutrons or photons') diff --git a/openmc/mesh.py b/openmc/mesh.py index ce5218b5e9..efe1c20a16 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -17,6 +17,7 @@ import openmc import openmc.checkvalue as cv from openmc.checkvalue import PathLike from openmc.utility_funcs import change_directory +from .bounding_box import BoundingBox from ._xml import get_elem_list, get_text from .mixin import IDManagerMixin from .surface import _BOUNDARY_TYPES @@ -40,6 +41,11 @@ class MeshMaterialVolumes(Mapping): Array of shape (elements, max_materials) storing material IDs volumes : numpy.ndarray Array of shape (elements, max_materials) storing material volumes + bboxes : numpy.ndarray, optional + Array of shape (elements, max_materials, 6) storing axis-aligned + bounding boxes for each (element, material) combination with ordering + (xmin, ymin, zmin, xmax, ymax, zmax). Bounding boxes enclose the + ray-estimator prisms used to compute volumes. See Also -------- @@ -64,9 +70,30 @@ class MeshMaterialVolumes(Mapping): [(2, 31.87963824195591), (1, 6.129949130817542)] """ - def __init__(self, materials: np.ndarray, volumes: np.ndarray): + def __init__( + self, + materials: np.ndarray, + volumes: np.ndarray, + bboxes: np.ndarray | None = None + ): self._materials = materials self._volumes = volumes + self._bboxes = bboxes + + if self._bboxes is not None: + if self._bboxes.shape[:2] != self._materials.shape: + raise ValueError( + 'bboxes must have shape (elements, max_materials, 6) ' + 'matching materials/volumes.' + ) + if self._bboxes.shape[2] != 6: + raise ValueError( + 'bboxes must have shape (elements, max_materials, 6).' + ) + + @property + def has_bounding_boxes(self) -> bool: + return self._bboxes is not None @property def num_elements(self) -> int: @@ -92,7 +119,11 @@ class MeshMaterialVolumes(Mapping): volumes[indices] = self._volumes[indices, i] return volumes - def by_element(self, index_elem: int) -> list[tuple[int | None, float]]: + def by_element( + self, + index_elem: int, + include_bboxes: bool = False + ) -> list[tuple[int | None, float] | tuple[int | None, float, BoundingBox | None]]: """Get a list of volumes for each material within a specific element. Parameters @@ -102,15 +133,32 @@ class MeshMaterialVolumes(Mapping): Returns ------- - list of tuple of (material ID, volume) + list of tuple + If ``include_bboxes`` is False (default), returns tuples of + (material ID, volume). If ``include_bboxes`` is True, returns + tuples of (material ID, volume, bounding box). """ table_size = self._volumes.shape[1] - return [ - (m if m > -1 else None, self._volumes[index_elem, i]) - for i in range(table_size) - if (m := self._materials[index_elem, i]) != -2 - ] + if include_bboxes and self._bboxes is None: + raise ValueError('Bounding boxes were not computed for this object.') + + results = [] + for i in range(table_size): + m = self._materials[index_elem, i] + if m == -2: + continue + mat_id = m if m > -1 else None + vol = self._volumes[index_elem, i] + + if include_bboxes: + vals = self._bboxes[index_elem, i] + bbox = BoundingBox(vals[0:3], vals[3:6]) + results.append((mat_id, vol, bbox)) + else: + results.append((mat_id, vol)) + + return results def save(self, filename: PathLike): """Save material volumes to a .npz file. @@ -120,8 +168,10 @@ class MeshMaterialVolumes(Mapping): filename : path-like Filename where data will be saved """ - np.savez_compressed( - filename, materials=self._materials, volumes=self._volumes) + kwargs = {'materials': self._materials, 'volumes': self._volumes} + if self._bboxes is not None: + kwargs['bboxes'] = self._bboxes + np.savez_compressed(filename, **kwargs) @classmethod def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes: @@ -134,7 +184,8 @@ class MeshMaterialVolumes(Mapping): """ filedata = np.load(filename) - return cls(filedata['materials'], filedata['volumes']) + bboxes = filedata['bboxes'] if 'bboxes' in filedata.files else None + return cls(filedata['materials'], filedata['volumes'], bboxes) class MeshBase(IDManagerMixin, ABC): @@ -153,11 +204,17 @@ class MeshBase(IDManagerMixin, ABC): Unique identifier for the mesh name : str Name of the mesh + lower_left : Iterable of float + The lower-left coordinates + upper_right : Iterable of float + The upper-right coordinates bounding_box : openmc.BoundingBox Axis-aligned bounding box of the mesh as defined by the upper-right and lower-left coordinates. indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] + n_elements : int + Number of elements in the mesh """ next_id = 1 @@ -179,6 +236,16 @@ class MeshBase(IDManagerMixin, ABC): self._name = name else: self._name = '' + + @property + @abstractmethod + def lower_left(self): + pass + + @property + @abstractmethod + def upper_right(self): + pass @property def bounding_box(self) -> openmc.BoundingBox: @@ -188,6 +255,11 @@ class MeshBase(IDManagerMixin, ABC): @abstractmethod def indices(self): pass + + @property + @abstractmethod + def n_elements(self): + pass def __repr__(self): string = type(self).__name__ + '\n' @@ -366,6 +438,7 @@ class MeshBase(IDManagerMixin, ABC): model: openmc.Model, n_samples: int | tuple[int, int, int] = 10_000, max_materials: int = 4, + bounding_boxes: bool = False, **kwargs ) -> MeshMaterialVolumes: """Determine volume of materials in each mesh element. @@ -388,6 +461,11 @@ class MeshBase(IDManagerMixin, ABC): the x, y, and z dimensions. max_materials : int, optional Estimated maximum number of materials in any given mesh element. + bounding_boxes : bool, optional + Whether to compute an axis-aligned bounding box for each + (mesh element, material) combination. When enabled, the bounding + box encloses the ray-estimator prisms used for the volume + estimation. **kwargs : dict Keyword arguments passed to :func:`openmc.lib.init` @@ -419,7 +497,8 @@ class MeshBase(IDManagerMixin, ABC): # Compute material volumes volumes = mesh.material_volumes( - n_samples, max_materials, output=kwargs['output']) + n_samples, max_materials, output=kwargs['output'], + bounding_boxes=bounding_boxes) # Restore original tallies model.tallies = original_tallies @@ -557,10 +636,19 @@ class StructuredMesh(MeshBase): s0 = (slice(0, -1),)*ndim + (slice(None),) s1 = (slice(1, None),)*ndim + (slice(None),) return (vertices[s0] + vertices[s1]) / 2 + + @property + def n_elements(self): + return np.prod(self.dimension) @property def num_mesh_cells(self): - return np.prod(self.dimension) + warnings.warn( + "The 'num_mesh_cells' attribute is deprecated and will be removed in a future version. " + "Use 'n_elements' instead.", + FutureWarning, stacklevel=2 + ) + return self.n_elements def write_data_to_vtk(self, filename: PathLike, @@ -822,10 +910,10 @@ class StructuredMesh(MeshBase): """ cv.check_type('data label', label, str) - if dataset.size != self.num_mesh_cells: + if dataset.size != self.n_elements: raise ValueError( f"The size of the dataset '{label}' ({dataset.size}) should be" - f" equal to the number of mesh cells ({self.num_mesh_cells})" + f" equal to the number of mesh cells ({self.n_elements})" ) # accept a flat array as-is, assuming it is in the correct order diff --git a/openmc/model/model.py b/openmc/model/model.py index 6e4c1c5856..67da2f73a3 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1160,6 +1160,16 @@ class Model: y_min = (origin[y] - 0.5*width[1]) * axis_scaling_factor[axis_units] y_max = (origin[y] + 0.5*width[1]) * axis_scaling_factor[axis_units] + # Determine whether any materials contains macroscopic data and if so, + # set energy mode accordingly and check that mg cross sections path is accessible + for mat in self.geometry.get_all_materials().values(): + if mat._macroscopic is not None: + self.settings.energy_mode = 'multi-group' + if 'mg_cross_sections' not in openmc.config: + raise RuntimeError("'mg_cross_sections' path must be set in " + "openmc.config before plotting.") + break + # Get ID map from the C API id_map = self.id_map( origin=origin, @@ -1179,8 +1189,8 @@ class Model: # Convert ID map to RGB image img = id_map_to_rgb( - id_map=id_map, - color_by=color_by, + id_map=id_map, + color_by=color_by, colors=colors, overlap_color=overlap_color ) @@ -1217,7 +1227,7 @@ class Model: extent=(x_min, x_max, y_min, y_max), **contour_kwargs ) - + # If only showing outline, set the axis limits and aspect explicitly if outline == 'only': axes.set_xlim(x_min, x_max) @@ -1685,6 +1695,87 @@ class Model: self.geometry.get_all_materials().values() ) + def _auto_generate_mgxs_lib( + self, + model: openmc.model.model, + groups: openmc.mgxs.EnergyGroups, + correction: str | none, + directory: pathlike, + ) -> openmc.mgxs.Library: + """ + Automatically generate a multi-group cross section libray from a model + with the specified group structure. + + Parameters + ---------- + groups : openmc.mgxs.EnergyGroups + Energy group structure for the MGXS. + nparticles : int + Number of particles to simulate per batch when generating MGXS. + mgxs_path : str + Filename for the MGXS HDF5 file. + correction : str + Transport correction to apply to the MGXS. Options are None and + "P0". + directory : str + Directory to run the simulation in, so as to contain XML files. + + Returns + ------- + mgxs_lib : openmc.mgxs.Library + OpenMC MGXS Library object + """ + + # Initialize MGXS library with a finished OpenMC geometry object + mgxs_lib = openmc.mgxs.Library(model.geometry) + + # Pick energy group structure + mgxs_lib.energy_groups = groups + + # Disable transport correction + mgxs_lib.correction = correction + + # Specify needed cross sections for random ray + if correction == 'P0': + mgxs_lib.mgxs_types = [ + 'nu-transport', 'absorption', 'nu-fission', 'fission', + 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi', + 'kappa-fission' + ] + elif correction is None: + mgxs_lib.mgxs_types = [ + 'total', 'absorption', 'nu-fission', 'fission', + 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi', + 'kappa-fission' + ] + + # Specify a "material" domain type for the cross section tally filters + mgxs_lib.domain_type = "material" + + # Specify the domains over which to compute multi-group cross sections + mgxs_lib.domains = model.geometry.get_all_materials().values() + + # Do not compute cross sections on a nuclide-by-nuclide basis + mgxs_lib.by_nuclide = False + + # Check the library - if no errors are raised, then the library is satisfactory. + mgxs_lib.check_library_for_openmc_mgxs() + + # Construct all tallies needed for the multi-group cross section library + mgxs_lib.build_library() + + # Create a "tallies.xml" file for the MGXS Library + mgxs_lib.add_to_tallies(model.tallies, merge=True) + + # Run + statepoint_filename = model.run(cwd=directory) + + # Load MGXS + with openmc.StatePoint(statepoint_filename) as sp: + mgxs_lib.load_from_statepoint(sp) + + return mgxs_lib + def _create_mgxs_sources( self, groups: openmc.mgxs.EnergyGroups, @@ -1848,52 +1939,8 @@ class Model: model.geometry.root_universe = infinite_universe # Add MGXS Tallies - - # Initialize MGXS library with a finished OpenMC geometry object - mgxs_lib = openmc.mgxs.Library(model.geometry) - - # Pick energy group structure - mgxs_lib.energy_groups = groups - - # Disable transport correction - mgxs_lib.correction = correction - - # Specify needed cross sections for random ray - if correction == 'P0': - mgxs_lib.mgxs_types = [ - 'nu-transport', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - elif correction is None: - mgxs_lib.mgxs_types = [ - 'total', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - - # Specify a "cell" domain type for the cross section tally filters - mgxs_lib.domain_type = "material" - - # Specify the cell domains over which to compute multi-group cross sections - mgxs_lib.domains = model.geometry.get_all_materials().values() - - # Do not compute cross sections on a nuclide-by-nuclide basis - mgxs_lib.by_nuclide = False - - # Check the library - if no errors are raised, then the library is satisfactory. - mgxs_lib.check_library_for_openmc_mgxs() - - # Construct all tallies needed for the multi-group cross section library - mgxs_lib.build_library() - - # Create a "tallies.xml" file for the MGXS Library - mgxs_lib.add_to_tallies(model.tallies, merge=True) - - # Run - statepoint_filename = model.run(cwd=directory) - - # Load MGXS - with openmc.StatePoint(statepoint_filename) as sp: - mgxs_lib.load_from_statepoint(sp) + mgxs_lib = self._auto_generate_mgxs_lib( + model, groups, correction, directory) # Create a MGXS File which can then be written to disk mgxs_set = mgxs_lib.get_xsdata(domain=material, xsdata_name=name) @@ -2055,48 +2102,8 @@ class Model: model.settings.output = {'summary': True, 'tallies': False} # Add MGXS Tallies - - # Initialize MGXS library with a finished OpenMC geometry object - mgxs_lib = openmc.mgxs.Library(model.geometry) - - # Pick energy group structure - mgxs_lib.energy_groups = groups - - # Disable transport correction - mgxs_lib.correction = correction - - # Specify needed cross sections for random ray - if correction == 'P0': - mgxs_lib.mgxs_types = ['nu-transport', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'] - elif correction is None: - mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'] - - # Specify a "cell" domain type for the cross section tally filters - mgxs_lib.domain_type = "material" - - # Specify the cell domains over which to compute multi-group cross sections - mgxs_lib.domains = model.geometry.get_all_materials().values() - - # Do not compute cross sections on a nuclide-by-nuclide basis - mgxs_lib.by_nuclide = False - - # Check the library - if no errors are raised, then the library is satisfactory. - mgxs_lib.check_library_for_openmc_mgxs() - - # Construct all tallies needed for the multi-group cross section library - mgxs_lib.build_library() - - # Create a "tallies.xml" file for the MGXS Library - mgxs_lib.add_to_tallies(model.tallies, merge=True) - - # Run - statepoint_filename = model.run(cwd=directory) - - # Load MGXS - with openmc.StatePoint(statepoint_filename) as sp: - mgxs_lib.load_from_statepoint(sp) + mgxs_lib = self._auto_generate_mgxs_lib( + model, groups, correction, directory) names = [mat.name for mat in mgxs_lib.domains] @@ -2146,52 +2153,8 @@ class Model: model.settings.output = {'summary': True, 'tallies': False} # Add MGXS Tallies - - # Initialize MGXS library with a finished OpenMC geometry object - mgxs_lib = openmc.mgxs.Library(model.geometry) - - # Pick energy group structure - mgxs_lib.energy_groups = groups - - # Disable transport correction - mgxs_lib.correction = correction - - # Specify needed cross sections for random ray - if correction == 'P0': - mgxs_lib.mgxs_types = [ - 'nu-transport', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - elif correction is None: - mgxs_lib.mgxs_types = [ - 'total', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - - # Specify a "cell" domain type for the cross section tally filters - mgxs_lib.domain_type = "material" - - # Specify the cell domains over which to compute multi-group cross sections - mgxs_lib.domains = model.geometry.get_all_materials().values() - - # Do not compute cross sections on a nuclide-by-nuclide basis - mgxs_lib.by_nuclide = False - - # Check the library - if no errors are raised, then the library is satisfactory. - mgxs_lib.check_library_for_openmc_mgxs() - - # Construct all tallies needed for the multi-group cross section library - mgxs_lib.build_library() - - # Create a "tallies.xml" file for the MGXS Library - mgxs_lib.add_to_tallies(model.tallies, merge=True) - - # Run - statepoint_filename = model.run(cwd=directory) - - # Load MGXS - with openmc.StatePoint(statepoint_filename) as sp: - mgxs_lib.load_from_statepoint(sp) + mgxs_lib = self._auto_generate_mgxs_lib( + model, groups, correction, directory) names = [mat.name for mat in mgxs_lib.domains] @@ -2617,5 +2580,3 @@ class SearchResult: def total_batches(self) -> int: """Total number of active batches used across all evaluations.""" return sum(self.batches) - - diff --git a/openmc/particle_restart.py b/openmc/particle_restart.py index 317d8761b6..44864107d9 100644 --- a/openmc/particle_restart.py +++ b/openmc/particle_restart.py @@ -1,6 +1,7 @@ import h5py import openmc.checkvalue as cv +from .particle_type import ParticleType _VERSION_PARTICLE_RESTART = 2 @@ -28,8 +29,8 @@ class Particle: Type of simulation (criticality or fixed source) id : long Identifier of the particle - type : int - Particle type (1 = neutron, 2 = photon, 3 = electron, 4 = positron) + type : openmc.ParticleType + Particle type weight : float Weight of the particle energy : float @@ -52,7 +53,7 @@ class Particle: self.energy = f['energy'][()] self.generations_per_batch = f['generations_per_batch'][()] self.id = f['id'][()] - self.type = f['type'][()] + self.type = ParticleType(f['type'][()]) self.n_particles = f['n_particles'][()] self.run_mode = f['run_mode'][()].decode() self.uvw = f['uvw'][()] diff --git a/openmc/particle_type.py b/openmc/particle_type.py new file mode 100644 index 0000000000..9d77134430 --- /dev/null +++ b/openmc/particle_type.py @@ -0,0 +1,228 @@ +from numbers import Integral + +from openmc.data import gnds_name, zam, ATOMIC_SYMBOL + + +_PDG_NAME = { + 2112: 'neutron', + 22: 'photon', + 11: 'electron', + -11: 'positron', + 2212: 'H1', +} + +_ALIAS_PDG = { + 'neutron': 2112, + 'n': 2112, + 'photon': 22, + 'gamma': 22, + 'electron': 11, + 'positron': -11, + 'proton': 2212, + 'p': 2212, + 'h1': 2212, + 'deuteron': 1000010020, + 'd': 1000010020, + 'h2': 1000010020, + 'triton': 1000010030, + 't': 1000010030, + 'h3': 1000010030, + 'alpha': 1000020040, + 'he4': 1000020040, +} + +_LEGACY_PARTICLE_INDEX = { + 0: 2112, + 1: 22, + 2: 11, + 3: -11, +} + + +class ParticleType: + """Particle type defined by a PDG number. + + ParticleType uses the Particle Data Group (PDG) Monte Carlo numbering scheme + to uniquely identify particle types. This includes elementary particles + (neutrons, photons, etc.) and nuclear codes for isotopes. + + Parameters + ---------- + value : str, int, or ParticleType + The particle identifier. Can be: + + - A string name (e.g., 'neutron', 'photon', 'He4', 'U235') + - An integer PDG number (e.g., 2112 for neutron) + - A string with PDG prefix (e.g., 'pdg:2112') + - An existing ParticleType instance + + Attributes + ---------- + pdg_number : int + The PDG number for this particle type + zam : tuple of int or None + For nuclear particles, the (Z, A, m) tuple where Z is atomic number, + A is mass number, and m is metastable state. None for elementary particles. + is_nucleus : bool + Whether this particle is a nucleus (ion) + + Examples + -------- + >>> neutron = ParticleType('neutron') + >>> neutron.pdg_number + 2112 + >>> he4 = ParticleType('He4') + >>> he4.zam + (2, 4, 0) + >>> ParticleType(2112) == ParticleType('neutron') + True + + """ + + __slots__ = ('_pdg_number',) + + def __init__(self, value: 'str | int | ParticleType'): + if isinstance(value, ParticleType): + pdg = value._pdg_number + elif isinstance(value, str): + pdg = self._pdg_number_from_string(value) + elif isinstance(value, Integral): + pdg = int(value) + # Handle legacy particle indices (0, 1, 2, 3) + if pdg in _LEGACY_PARTICLE_INDEX: + pdg = _LEGACY_PARTICLE_INDEX[pdg] + else: + raise TypeError(f"Cannot create ParticleType from {type(value).__name__}") + + self._pdg_number = pdg + + def __eq__(self, other): + if isinstance(other, ParticleType): + return self._pdg_number == other._pdg_number + if isinstance(other, Integral): + return self._pdg_number == int(other) + if isinstance(other, str): + try: + return self._pdg_number == ParticleType(other)._pdg_number + except (ValueError, TypeError): + return False + return NotImplemented + + def __hash__(self) -> int: + return hash(self._pdg_number) + + def __int__(self) -> int: + return self._pdg_number + + @property + def pdg_number(self) -> int: + return self._pdg_number + + @staticmethod + def _pdg_number_from_string(value: str) -> int: + """Parse a string to get a PDG number. + + Parameters + ---------- + value : str + Particle identifier string + + Returns + ------- + int + PDG number + + Raises + ------ + ValueError + If string cannot be parsed as a valid particle identifier + + """ + s = value.strip() + if not s: + raise ValueError('Particle identifier cannot be empty.') + + lower = s.lower() + if lower.startswith('pdg:'): + code_str = lower[4:] + try: + return int(code_str) + except ValueError: + raise ValueError(f'Invalid PDG number: {code_str}') + + if lower in _ALIAS_PDG: + return _ALIAS_PDG[lower] + + # Assume it is a GNDS nuclide name + Z, A, m = zam(s) + if Z <= 0 or Z > 999 or A <= 0 or A > 999 or m < 0 or m > 9: + raise ValueError('Invalid Z/A/m for nuclear PDG number.') + return 1000000000 + Z * 10000 + A * 10 + m + + def __repr__(self) -> str: + return f'' + + def __str__(self) -> str: + """Return a canonical string representation of the particle type. + + Returns + ------- + str + Canonical name (e.g., 'neutron', 'He4', 'pdg:12345') + + """ + if self._pdg_number in _PDG_NAME: + return _PDG_NAME[self._pdg_number] + + if (zam_tuple := self.zam) is not None: + Z, A, m = zam_tuple + if Z <= 0 or Z > max(ATOMIC_SYMBOL) or A <= 0 or A > 999: + raise ValueError(f"Invalid nuclear PDG number: {self._pdg_number}") + return gnds_name(Z, A, m) + + return f'pdg:{self._pdg_number}' + + @property + def zam(self) -> 'tuple[int, int, int] | None': + """Return the (Z, A, m) tuple for nuclear particles. + + Returns + ------- + tuple of int or None + For nuclear particles, returns (Z, A, m) where Z is atomic number, + A is mass number, and m is metastable state. Returns None for + elementary particles. + + """ + if self._pdg_number < 1000000000: + return None + Z = (self._pdg_number // 10000) % 1000 + A = (self._pdg_number // 10) % 1000 + m = self._pdg_number % 10 + if Z <= 0 or A <= 0: + return None + else: + return (Z, A, m) + + @property + def is_nucleus(self) -> bool: + """Return whether this particle is a nucleus. + + Returns + ------- + bool + True if the particle is a nucleus (ion), False otherwise + + """ + return self.zam is not None + + +# Define common particle constants +ParticleType.NEUTRON = ParticleType(2112) +ParticleType.PHOTON = ParticleType(22) +ParticleType.ELECTRON = ParticleType(11) +ParticleType.POSITRON = ParticleType(-11) +ParticleType.PROTON = ParticleType(2212) +ParticleType.DEUTERON = ParticleType(1000010020) +ParticleType.TRITON = ParticleType(1000010030) +ParticleType.ALPHA = ParticleType(1000020040) diff --git a/openmc/settings.py b/openmc/settings.py index a022b045f2..3cf662e311 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -4,13 +4,14 @@ import itertools from math import ceil from numbers import Integral, Real from pathlib import Path +import traceback import lxml.etree as ET import warnings import openmc import openmc.checkvalue as cv from openmc.checkvalue import PathLike -from openmc.stats.multivariate import MeshSpatial, Box, PolarAzimuthal, Isotropic +from openmc.stats.multivariate import MeshSpatial from ._xml import clean_indentation, get_elem_list, get_text from .mesh import _read_meshes, RegularMesh, MeshBase from .source import SourceBase, MeshSource, IndependentSource @@ -467,6 +468,16 @@ class Settings: for key, value in kwargs.items(): setattr(self, key, value) + def __setattr__(self, name: str, value): + if not name.startswith('_'): + try: + getattr(self, name) + except AttributeError as e: + msg, = traceback.format_exception_only(e) + msg = msg.strip().split(maxsplit=1)[-1] + warnings.warn(msg, stacklevel=2) + super().__setattr__(name, value) + @property def run_mode(self) -> str: return self._run_mode.value diff --git a/openmc/source.py b/openmc/source.py index 84d8a9619d..a11cfd6c8d 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -1,7 +1,6 @@ from __future__ import annotations from abc import ABC, abstractmethod from collections.abc import Iterable, Sequence -from enum import IntEnum from numbers import Real from pathlib import Path import warnings @@ -19,6 +18,8 @@ from openmc.stats.multivariate import UnitSphere, Spatial from openmc.stats.univariate import Univariate from ._xml import get_elem_list, get_text from .mesh import MeshBase, StructuredMesh, UnstructuredMesh +from .particle_type import ParticleType +from .statepoint import _VERSION_STATEPOINT from .utility_funcs import input_path @@ -265,8 +266,8 @@ class IndependentSource(SourceBase): time distribution of source sites strength : float Strength of the source - particle : {'neutron', 'photon', 'electron', 'positron'} - Source particle type + particle : str or int or openmc.ParticleType + Source particle type (name, PDG number, or type) domains : iterable of openmc.Cell, openmc.Material, or openmc.Universe Domains to reject based on, i.e., if a sampled spatial location is not within one of these domains, it will be rejected. @@ -302,10 +303,9 @@ class IndependentSource(SourceBase): type : str Indicator of source type: 'independent' - .. versionadded:: 0.14.0 - - particle : {'neutron', 'photon', 'electron', 'positron'} - Source particle type + .. versionadded:: 0.14.0 + particle : str or int or openmc.ParticleType + Source particle type (alias, PDG number, or GNDS nuclide name) constraints : dict Constraints on sampled source particles. Valid keys include 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', @@ -320,7 +320,7 @@ class IndependentSource(SourceBase): energy: openmc.stats.Univariate | None = None, time: openmc.stats.Univariate | None = None, strength: float = 1.0, - particle: str = 'neutron', + particle: str | int | ParticleType = 'neutron', domains: Sequence[openmc.Cell | openmc.Material | openmc.Universe] | None = None, constraints: dict[str, Any] | None = None @@ -405,14 +405,12 @@ class IndependentSource(SourceBase): self._time = time @property - def particle(self): + def particle(self) -> ParticleType: return self._particle @particle.setter def particle(self, particle): - cv.check_value('source particle', particle, - ['neutron', 'photon', 'electron', 'positron']) - self._particle = particle + self._particle = ParticleType(particle) def populate_xml_element(self, element): """Add necessary source information to an XML element @@ -423,7 +421,7 @@ class IndependentSource(SourceBase): XML element containing source data """ - element.set("particle", self.particle) + element.set("particle", str(self.particle)) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -572,10 +570,10 @@ class MeshSource(SourceBase): s = np.asarray(s) if isinstance(self.mesh, StructuredMesh): - if s.size != self.mesh.num_mesh_cells: + if s.size != self.mesh.n_elements: raise ValueError( f'The length of the source array ({s.size}) does not match ' - f'the number of mesh elements ({self.mesh.num_mesh_cells}).') + f'the number of mesh elements ({self.mesh.n_elements}).') # If user gave a multidimensional array, flatten in the order # of the mesh indices @@ -898,76 +896,6 @@ class FileSource(SourceBase): return cls(**kwargs) -class ParticleType(IntEnum): - """ - IntEnum class representing a particle type. Type - values mirror those found in the C++ class. - """ - NEUTRON = 0 - PHOTON = 1 - ELECTRON = 2 - POSITRON = 3 - - @classmethod - def from_string(cls, value: str): - """ - Constructs a ParticleType instance from a string. - - Parameters - ---------- - value : str - The string representation of the particle type. - - Returns - ------- - The corresponding ParticleType instance. - """ - try: - return cls[value.upper()] - except KeyError: - raise ValueError( - f"Invalid string for creation of {cls.__name__}: {value}") - - @classmethod - def from_pdg_number(cls, pdg_number: int) -> ParticleType: - """Constructs a ParticleType instance from a PDG number. - - The Particle Data Group at LBNL publishes a Monte Carlo particle - numbering scheme as part of the `Review of Particle Physics - <10.1103/PhysRevD.110.030001>`_. This method maps PDG numbers to the - corresponding :class:`ParticleType`. - - Parameters - ---------- - pdg_number : int - The PDG number of the particle type. - - Returns - ------- - The corresponding ParticleType instance. - """ - try: - return { - 2112: ParticleType.NEUTRON, - 22: ParticleType.PHOTON, - 11: ParticleType.ELECTRON, - -11: ParticleType.POSITRON, - }[pdg_number] - except KeyError: - raise ValueError(f"Unrecognized PDG number: {pdg_number}") - - def __repr__(self) -> str: - """ - Returns a string representation of the ParticleType instance. - - Returns: - str: The lowercase name of the ParticleType instance. - """ - return self.name.lower() - - # needed for < Python 3.11 - def __str__(self) -> str: - return self.__repr__() class SourceParticle: @@ -992,8 +920,8 @@ class SourceParticle: Delayed group particle was created in (neutrons only) surf_id : int Surface ID where particle is at, if any. - particle : ParticleType - Type of the particle + particle : ParticleType or str or int + Type of the particle (type, name, or PDG number) """ @@ -1006,7 +934,7 @@ class SourceParticle: wgt: float = 1.0, delayed_group: int = 0, surf_id: int = 0, - particle: ParticleType = ParticleType.NEUTRON + particle: ParticleType | str | int = ParticleType.NEUTRON ): self.r = tuple(r) @@ -1018,9 +946,16 @@ class SourceParticle: self.surf_id = surf_id self.particle = particle + @property + def particle(self) -> ParticleType: + return self._particle + + @particle.setter + def particle(self, particle): + self._particle = ParticleType(particle) + def __repr__(self): - name = self.particle.name.lower() - return f'' + return f'' def to_tuple(self) -> tuple: """Return source particle attributes as a tuple @@ -1032,7 +967,7 @@ class SourceParticle: """ return (self.r, self.u, self.E, self.time, self.wgt, - self.delayed_group, self.surf_id, self.particle.value) + self.delayed_group, self.surf_id, self.particle.pdg_number) def write_source_file( @@ -1116,12 +1051,7 @@ class ParticleList(list): particles = [] with mcpl.MCPLFile(filename) as f: for particle in f.particles: - # Determine particle type based on the PDG number - try: - particle_type = ParticleType.from_pdg_number( - particle.pdgcode) - except ValueError: - particle_type = "UNKNOWN" + particle_type = ParticleType(particle.pdgcode) # Create a source particle instance. Note that MCPL stores # energy in MeV and time in ms. @@ -1179,7 +1109,7 @@ class ParticleList(list): # Extract the attributes of the source particles into a list of tuples data = [(sp.r[0], sp.r[1], sp.r[2], sp.u[0], sp.u[1], sp.u[2], sp.E, sp.time, sp.wgt, sp.delayed_group, sp.surf_id, - sp.particle.name.lower()) for sp in self] + str(sp.particle)) for sp in self] # Define the column names for the DataFrame columns = ['x', 'y', 'z', 'u_x', 'u_y', 'u_z', 'E', 'time', 'wgt', @@ -1226,6 +1156,7 @@ class ParticleList(list): kwargs.setdefault('mode', 'w') with h5py.File(filename, **kwargs) as fh: fh.attrs['filetype'] = np.bytes_("source") + fh.attrs['version'] = np.array([_VERSION_STATEPOINT, 2]) fh.create_dataset('source_bank', data=arr, dtype=source_dtype) @@ -1337,7 +1268,7 @@ def read_collision_track_mcpl(file_path): data['material_id'].append(int(values_dict.get('material_id', 0))) data['universe_id'].append(int(values_dict.get('universe_id', 0))) data['n_collision'].append(int(values_dict.get('n_collision', 0))) - data['particle'].append(ParticleType.from_pdg_number(p.pdgcode)) + data['particle'].append(ParticleType(p.pdgcode)) data['parent_id'].append(int(values_dict.get('parent_id', 0))) data['progeny_id'].append(int(values_dict.get('progeny_id', 0))) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 272367f6f2..94258c833c 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -397,7 +397,7 @@ class Discrete(Univariate): def merge( cls, dists: Sequence[Discrete], - probs: Sequence[int] + probs: Sequence[float] ): """Merge multiple discrete distributions into a single distribution @@ -1897,7 +1897,7 @@ class Mixture(Univariate): def combine_distributions( - dists: Sequence[Univariate], + dists: Sequence[Discrete | Tabular], probs: Sequence[float] ): """Combine distributions with specified probabilities @@ -1912,41 +1912,40 @@ def combine_distributions( Parameters ---------- - dists : iterable of openmc.stats.Univariate + dists : sequence of openmc.stats.Discrete or openmc.stats.Tabular Distributions to combine - probs : iterable of float + probs : sequence of float Probability (or intensity) of each distribution """ - # Get copy of distribution list so as not to modify the argument - dist_list = deepcopy(dists) + for i, dist in enumerate(dists): + cv.check_type(f'dists[{i}]', dist, (Discrete, Tabular)) + cv.check_type(f'probs[{i}]', probs[i], Real) + cv.check_greater_than(f'probs[{i}]', probs[i], 0.0) # Get list of discrete/continuous distribution indices - discrete_index = [i for i, d in enumerate(dist_list) if isinstance(d, Discrete)] - cont_index = [i for i, d in enumerate(dist_list) if isinstance(d, Tabular)] + discrete_index = [i for i, d in enumerate(dists) if isinstance(d, Discrete)] + cont_index = [i for i, d in enumerate(dists) if isinstance(d, Tabular)] - # Apply probabilites to continuous distributions - for i in cont_index: - dist = dist_list[i] - dist._p *= probs[i] + cont_dists = [dists[i] for i in cont_index] + cont_probs = [probs[i] for i in cont_index] if discrete_index: # Create combined discrete distribution - dist_discrete = [dist_list[i] for i in discrete_index] + dist_discrete = [dists[i] for i in discrete_index] discrete_probs = [probs[i] for i in discrete_index] combined_dist = Discrete.merge(dist_discrete, discrete_probs) - - # Replace multiple discrete distributions with merged - for idx in reversed(discrete_index): - dist_list.pop(idx) - dist_list.append(combined_dist) - - # Combine discrete and continuous if present - if len(dist_list) > 1: - probs = [1.0]*len(dist_list) - dist_list[:] = [Mixture(probs, dist_list.copy())] - - return dist_list[0] + if cont_index: + return Mixture(cont_probs + [1.0], cont_dists + [combined_dist]) + else: + return combined_dist + else: + if len(cont_dists) == 1: + dist = cont_dists[0] + return Tabular(dist.x, dist.p * cont_probs[0], + dist.interpolation, bias=dist.bias) + else: + return Mixture(cont_probs, cont_dists) def check_bias_support(parent: Univariate, bias: Univariate | None): diff --git a/openmc/tallies.py b/openmc/tallies.py index 7b1ef02198..74d4722af4 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1982,7 +1982,7 @@ class Tally(IDManagerMixin): # Expand the columns into Pandas MultiIndices for readability if pd.__version__ >= '0.16': - columns = copy.deepcopy(df.columns.values) + columns = copy.deepcopy(list(df.columns.values)) # Convert all elements in columns list to tuples for i, column in enumerate(columns): diff --git a/openmc/tracks.py b/openmc/tracks.py index 81646e7d2e..44e3c0eba6 100644 --- a/openmc/tracks.py +++ b/openmc/tracks.py @@ -4,7 +4,8 @@ from collections.abc import Sequence import h5py from .checkvalue import check_filetype_version -from .source import SourceParticle, ParticleType +from .particle_type import ParticleType +from .source import SourceParticle from pathlib import Path @@ -25,7 +26,7 @@ states : numpy.ndarray """ def _particle_track_repr(self): - return f"" + return f"" ParticleTrack.__repr__ = _particle_track_repr @@ -92,8 +93,8 @@ class Track(Sequence): Parameters ---------- - particle : {'neutron', 'photon', 'electron', 'positron'} - Matching particle type + particle : str or int or openmc.ParticleType + Matching particle type (name, PDG number, or type) state_filter : function Function that takes a state (structured datatype) and returns a bool depending on some criteria. @@ -126,7 +127,7 @@ class Track(Sequence): for t in self: # Check for matching particle if particle is not None: - if t.particle.name.lower() != particle: + if t.particle != ParticleType(particle): continue # Apply arbitrary state filter @@ -184,7 +185,7 @@ class Track(Sequence): def sources(self): sources = [] for particle_track in self: - particle_type = ParticleType(particle_track.particle) + particle_type = particle_track.particle state = particle_track.states[0] sources.append( SourceParticle( diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py index 5d52a579a1..7797986df0 100644 --- a/openmc/weight_windows.py +++ b/openmc/weight_windows.py @@ -10,13 +10,12 @@ import numpy as np import h5py import openmc -from openmc.filter import _PARTICLES from openmc.mesh import MeshBase, RectilinearMesh, CylindricalMesh, SphericalMesh, UnstructuredMesh import openmc.checkvalue as cv from openmc.checkvalue import PathLike from ._xml import get_elem_list, get_text, clean_indentation from .mixin import IDManagerMixin -from .utility_funcs import change_directory +from .particle_type import ParticleType class WeightWindows(IDManagerMixin): @@ -51,7 +50,7 @@ class WeightWindows(IDManagerMixin): A list of values for which each successive pair constitutes a range of energies in [eV] for a single bin. If no energy bins are provided, the maximum and minimum energy for the data available at runtime. - particle_type : {'neutron', 'photon'} + particle_type : str or int or openmc.ParticleType Particle type the weight windows apply to survival_ratio : float Ratio of the survival weight to the lower weight window bound for @@ -116,7 +115,7 @@ class WeightWindows(IDManagerMixin): upper_ww_bounds: Iterable[float] | None = None, upper_bound_ratio: float | None = None, energy_bounds: Iterable[Real] | None = None, - particle_type: str = 'neutron', + particle_type: str | int | openmc.ParticleType = 'neutron', survival_ratio: float = 3.0, max_lower_bound_ratio: float | None = None, max_split: int = 10, @@ -213,13 +212,15 @@ class WeightWindows(IDManagerMixin): self._mesh = mesh @property - def particle_type(self) -> str: + def particle_type(self) -> ParticleType: return self._particle_type @particle_type.setter - def particle_type(self, pt: str): - cv.check_value('Particle type', pt, _PARTICLES) - self._particle_type = pt + def particle_type(self, pt): + ptype = ParticleType(pt) + if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}: + raise ValueError("Weight windows can only be applied for neutrons or photons") + self._particle_type = ptype @property def energy_bounds(self) -> Iterable[Real]: @@ -329,7 +330,7 @@ class WeightWindows(IDManagerMixin): subelement.text = str(self.mesh.id) subelement = ET.SubElement(element, 'particle_type') - subelement.text = self.particle_type + subelement.text = str(self.particle_type) if self.energy_bounds is not None: subelement = ET.SubElement(element, 'energy_bounds') @@ -494,7 +495,7 @@ class WeightWindowGenerator: A list of values for which each successive pair constitutes a range of energies in [eV] for a single bin. If no energy bins are provided, the maximum and minimum energy for the data available at runtime. - particle_type : {'neutron', 'photon'} + particle_type : str or int or openmc.ParticleType Particle type the weight windows apply to method : {'magic', 'fw_cadis'} The weight window generation methodology applied during an update. @@ -513,7 +514,7 @@ class WeightWindowGenerator: energy_bounds : Iterable of Real A list of values for which each successive pair constitutes a range of energies in [eV] for a single bin - particle_type : {'neutron', 'photon'} + particle_type : openmc.ParticleType Particle type the weight windows apply to method : {'magic', 'fw_cadis'} The weight window generation methodology applied during an update. @@ -534,7 +535,7 @@ class WeightWindowGenerator: self, mesh: openmc.MeshBase, energy_bounds: Sequence[float] | None = None, - particle_type: str = 'neutron', + particle_type: str | int | openmc.ParticleType = 'neutron', method: str = 'magic', max_realizations: int = 1, update_interval: int = 1, @@ -555,7 +556,7 @@ class WeightWindowGenerator: def __repr__(self): string = type(self).__name__ + '\n' string += f'\t{"Mesh":<20}=\t{self.mesh.id}\n' - string += f'\t{"Particle:":<20}=\t{self.particle_type}\n' + string += f'\t{"Particle:":<20}=\t{str(self.particle_type)}\n' string += f'\t{"Energy Bounds:":<20}=\t{self.energy_bounds}\n' string += f'\t{"Method":<20}=\t{self.method}\n' string += f'\t{"Max Realizations:":<20}=\t{self.max_realizations}\n' @@ -586,13 +587,15 @@ class WeightWindowGenerator: self._energy_bounds = eb @property - def particle_type(self) -> str: + def particle_type(self) -> ParticleType: return self._particle_type @particle_type.setter - def particle_type(self, pt: str): - cv.check_value('particle type', pt, ('neutron', 'photon')) - self._particle_type = pt + def particle_type(self, pt): + ptype = ParticleType(pt) + if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}: + raise ValueError("Weight windows can only be applied for neutrons or photons") + self._particle_type = ptype @property def method(self) -> str: @@ -695,7 +698,7 @@ class WeightWindowGenerator: subelement = ET.SubElement(element, 'energy_bounds') subelement.text = ' '.join(str(e) for e in self.energy_bounds) particle_elem = ET.SubElement(element, 'particle_type') - particle_elem.text = self.particle_type + particle_elem.text = str(self.particle_type) realizations_elem = ET.SubElement(element, 'max_realizations') realizations_elem.text = str(self.max_realizations) update_interval_elem = ET.SubElement(element, 'update_interval') @@ -730,7 +733,7 @@ class WeightWindowGenerator: mesh_id = int(get_text(elem, 'mesh')) mesh = meshes[mesh_id] - + energy_bounds = get_elem_list(elem, "energy_bounds, float") particle_type = get_text(elem, 'particle_type') diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index a2320e0b46..d77066fb0e 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -31,13 +31,13 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (p.material() == MATERIAL_VOID) return; - int photon = static_cast(ParticleType::photon); + int photon = ParticleType::photon().transport_index(); if (p.E() < settings::energy_cutoff[photon]) return; // Get bremsstrahlung data for this material and particle type BremsstrahlungData* mat; - if (p.type() == ParticleType::positron) { + if (p.type() == ParticleType::positron()) { mat = &model::materials[p.material()]->ttb_->positron; } else { mat = &model::materials[p.material()]->ttb_->electron; @@ -119,7 +119,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) } // Create secondary photon - p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon); + p.create_secondary(p.wgt(), p.u(), w, ParticleType::photon()); *E_lost += w; } } diff --git a/src/cell.cpp b/src/cell.cpp index aceed31cab..ebe28c3d2c 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -660,7 +660,7 @@ Region::Region(std::string region_spec, int32_t cell_id) if (token == OP_UNION) { simple_ = false; // Ensure intersections have precedence over unions - add_precedence(); + enforce_precedence(); break; } } @@ -703,7 +703,7 @@ void Region::apply_demorgan( //! precedence than unions using parentheses. //============================================================================== -int64_t Region::add_parentheses(int64_t start) +void Region::add_parentheses(int64_t start) { int32_t start_token = expression_[start]; // Add left parenthesis and set new position to be after parenthesis @@ -712,14 +712,6 @@ int64_t Region::add_parentheses(int64_t start) } expression_.insert(expression_.begin() + start - 1, OP_LEFT_PAREN); - // Keep track of return iterator distance. If we don't encounter a left - // parenthesis, we return an iterator corresponding to wherever the right - // parenthesis is inserted. If a left parenthesis is encountered, an iterator - // corresponding to the left parenthesis is returned. Also note that we keep - // track of a *distance* instead of an iterator because the underlying memory - // allocation may change. - std::size_t return_it_dist = 0; - // Add right parenthesis // While the start iterator is within the bounds of infix while (start + 1 < expression_.size()) { @@ -733,7 +725,6 @@ int64_t Region::add_parentheses(int64_t start) // in the region, when the operator is an intersection then include the // operator and next surface if (expression_[start] == OP_LEFT_PAREN) { - return_it_dist = start; int depth = 1; do { start++; @@ -750,54 +741,73 @@ int64_t Region::add_parentheses(int64_t start) --start; } expression_.insert(expression_.begin() + start, OP_RIGHT_PAREN); - if (return_it_dist > 0) { - return return_it_dist; - } else { - return start - 1; - } + return; } } } - // If we get here a right parenthesis hasn't been placed, - // return iterator + // If we get here a right parenthesis hasn't been placed expression_.push_back(OP_RIGHT_PAREN); - if (return_it_dist > 0) { - return return_it_dist; - } else { - return start - 1; - } } +//============================================================================== +//! Add parentheses to enforce operator precedence in region expressions +//! +//! This function ensures that intersection operators have higher precedence +//! than union operators by adding parentheses where needed. For example: +//! "1 2 | 3" becomes "(1 2) | 3" +//! "1 | 2 3" becomes "1 | (2 3)" +//! +//! The algorithm uses stacks to track the current operator type and its +//! position at each parenthesis depth level. When it encounters a different +//! operator at the same depth, it adds parentheses to group the +//! higher-precedence operations. //============================================================================== -void Region::add_precedence() +void Region::enforce_precedence() { - int32_t current_op = 0; - std::size_t current_dist = 0; + // Stack tracking the operator type at each depth (0 = no operator seen yet) + vector op_stack = {0}; - for (int64_t i = 0; i < expression_.size(); i++) { + // Stack tracking where the operator sequence started at each depth + vector pos_stack = {0}; + + for (int64_t i = 0; i < expression_.size(); ++i) { int32_t token = expression_[i]; - if (token == OP_UNION || token == OP_INTERSECTION) { - if (current_op == 0) { - // Set the current operator if is hasn't been set - current_op = token; - current_dist = i; - } else if (token != current_op) { - // If the current operator doesn't match the token, add parenthesis to - // assert precedence - if (current_op == OP_INTERSECTION) { - i = add_parentheses(current_dist); - } else { - i = add_parentheses(i); - } - current_op = 0; - current_dist = 0; + if (token == OP_LEFT_PAREN) { + // Entering a new parenthesis level - push new tracking state + op_stack.push_back(0); + pos_stack.push_back(0); + continue; + } else if (token == OP_RIGHT_PAREN) { + // Exiting a parenthesis level - pop tracking state (keep at least one) + if (op_stack.size() > 1) { + op_stack.pop_back(); + pos_stack.pop_back(); + } + continue; + } + + if (token == OP_UNION || token == OP_INTERSECTION) { + if (op_stack.back() == 0) { + // First operator at this depth - record it and its position + op_stack.back() = token; + pos_stack.back() = i; + } else if (token != op_stack.back()) { + // Encountered a different operator at the same depth - need to add + // parentheses to enforce precedence. Intersection has higher + // precedence, so we parenthesize the intersection terms. + if (op_stack.back() == OP_INTERSECTION) { + add_parentheses(pos_stack.back()); + } else { + add_parentheses(i); + } + + // Restart the scan since we modified the expression + i = -1; // Will be incremented to 0 by the for loop + op_stack = {0}; + pos_stack = {0}; } - } else if (token > OP_COMPLEMENT) { - // If the token is a parenthesis reset the current operator - current_op = 0; - current_dist = 0; } } } diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 80f8f7de14..e25e08d748 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -426,6 +426,11 @@ SpatialBox::SpatialBox(pugi::xml_node node, bool fission) upper_right_ = Position {params[3], params[4], params[5]}; } +SpatialBox::SpatialBox(Position lower_left, Position upper_right, bool fission) + : lower_left_(lower_left), upper_right_(upper_right), + only_fissionable_(fission) +{} + std::pair SpatialBox::sample(uint64_t* seed) const { Position xi {prn(seed), prn(seed), prn(seed)}; diff --git a/src/material.cpp b/src/material.cpp index 54caa38409..072e6decad 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -819,9 +819,9 @@ void Material::calculate_xs(Particle& p) const p.macro_xs().fission = 0.0; p.macro_xs().nu_fission = 0.0; - if (p.type() == ParticleType::neutron) { + if (p.type().is_neutron()) { this->calculate_neutron_xs(p); - } else if (p.type() == ParticleType::photon) { + } else if (p.type().is_photon()) { this->calculate_photon_xs(p); } } @@ -829,7 +829,7 @@ void Material::calculate_xs(Particle& p) const void Material::calculate_neutron_xs(Particle& p) const { // Find energy index on energy grid - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); int i_grid = std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing; diff --git a/src/mcpl_interface.cpp b/src/mcpl_interface.cpp index 256f3343fc..30c41ec5af 100644 --- a/src/mcpl_interface.cpp +++ b/src/mcpl_interface.cpp @@ -321,25 +321,7 @@ inline void ensure_mcpl_ready_or_fatal() SourceSite mcpl_particle_to_site(const mcpl_particle_repr_t* particle_repr) { SourceSite site; - switch (particle_repr->pdgcode) { - case 2112: - site.particle = ParticleType::neutron; - break; - case 22: - site.particle = ParticleType::photon; - break; - case 11: - site.particle = ParticleType::electron; - break; - case -11: - site.particle = ParticleType::positron; - break; - default: - fatal_error(fmt::format( - "MCPL: Encountered unexpected PDG code {} when converting to SourceSite.", - particle_repr->pdgcode)); - break; - } + site.particle = ParticleType {particle_repr->pdgcode}; // Copy position and direction site.r.x = particle_repr->position[0]; @@ -368,7 +350,6 @@ vector mcpl_source_sites(std::string path) } size_t n_particles_in_file = g_mcpl_api->hdr_nparticles(mcpl_file); - size_t n_skipped = 0; if (n_particles_in_file > 0) { sites.reserve(n_particles_in_file); } @@ -381,31 +362,16 @@ vector mcpl_source_sites(std::string path) path, sites.size(), n_particles_in_file)); break; } - if (p_repr->pdgcode == 2112 || p_repr->pdgcode == 22 || - p_repr->pdgcode == 11 || p_repr->pdgcode == -11) { - sites.push_back(mcpl_particle_to_site(p_repr)); - } else { - n_skipped++; - } + sites.push_back(mcpl_particle_to_site(p_repr)); } g_mcpl_api->close_file(mcpl_file); - if (n_skipped > 0 && n_particles_in_file > 0) { - double percent_skipped = - 100.0 * static_cast(n_skipped) / n_particles_in_file; - warning(fmt::format( - "MCPL: Skipped {} of {} total particles ({:.1f}%) in file '{}' because " - "their type is not supported by OpenMC.", - n_skipped, n_particles_in_file, percent_skipped, path)); - } - if (sites.empty()) { if (n_particles_in_file > 0) { fatal_error(fmt::format( - "MCPL file '{}' contained {} particles, but none were of the supported " - "types (neutron, photon, electron, positron). OpenMC cannot proceed " - "without source particles.", + "MCPL file '{}' contained {} particles, but no particles could be " + "read.", path, n_particles_in_file)); } else { fatal_error(fmt::format( @@ -461,22 +427,7 @@ void write_mcpl_source_bank_internal(mcpl_outfile_t* file_id, p_repr.ekin = site.E * 1e-6; p_repr.time = site.time * 1e3; p_repr.weight = site.wgt; - switch (site.particle) { - case ParticleType::neutron: - p_repr.pdgcode = 2112; - break; - case ParticleType::photon: - p_repr.pdgcode = 22; - break; - case ParticleType::electron: - p_repr.pdgcode = 11; - break; - case ParticleType::positron: - p_repr.pdgcode = -11; - break; - default: - continue; - } + p_repr.pdgcode = site.particle.pdg_number(); g_mcpl_api->add_particle(file_id, &p_repr); } } @@ -633,22 +584,7 @@ void write_mcpl_collision_track_internal(mcpl_outfile_t* file_id, p_repr.ekin = site.E * 1e-6; p_repr.time = site.time * 1e3; p_repr.weight = site.wgt; - switch (site.particle) { - case ParticleType::neutron: - p_repr.pdgcode = 2112; - break; - case ParticleType::photon: - p_repr.pdgcode = 22; - break; - case ParticleType::electron: - p_repr.pdgcode = 11; - break; - case ParticleType::positron: - p_repr.pdgcode = -11; - break; - default: - continue; - } + p_repr.pdgcode = site.particle.pdg_number(); g_mcpl_api->add_particle(file_id, &p_repr); } } else { diff --git a/src/mesh.cpp b/src/mesh.cpp index 9a7a8e7517..a72b1441f7 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1,6 +1,8 @@ #include "openmc/mesh.h" #include // for copy, equal, min, min_element #include +#include // for uint64_t +#include // for memcpy #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers #include // for ceil #include // for size_t @@ -140,6 +142,63 @@ inline bool atomic_cas_int32(int32_t* ptr, int32_t& expected, int32_t desired) #endif } +// Helper function equivalent to std::bit_cast in C++20 +template +inline To bit_cast_value(const From& value) +{ + To out; + std::memcpy(&out, &value, sizeof(To)); + return out; +} + +inline void atomic_update_double(double* ptr, double value, bool is_min) +{ +#if defined(__GNUC__) || defined(__clang__) + using may_alias_uint64_t [[gnu::may_alias]] = uint64_t; + auto* bits_ptr = reinterpret_cast(ptr); + uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST); + double current = bit_cast_value(current_bits); + while (is_min ? (value < current) : (value > current)) { + uint64_t desired_bits = bit_cast_value(value); + uint64_t expected_bits = current_bits; + if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits, + false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) { + return; + } + current_bits = expected_bits; + current = bit_cast_value(current_bits); + } + +#elif defined(_MSC_VER) + auto* bits_ptr = reinterpret_cast(ptr); + long long current_bits = *bits_ptr; + double current = bit_cast_value(current_bits); + while (is_min ? (value < current) : (value > current)) { + long long desired_bits = bit_cast_value(value); + long long old_bits = + _InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits); + if (old_bits == current_bits) { + return; + } + current_bits = old_bits; + current = bit_cast_value(current_bits); + } + +#else +#error "No compare-and-swap implementation available for this compiler." +#endif +} + +inline void atomic_max_double(double* ptr, double value) +{ + atomic_update_double(ptr, value, false); +} + +inline void atomic_min_double(double* ptr, double value) +{ + atomic_update_double(ptr, value, true); +} + namespace detail { //============================================================================== @@ -147,7 +206,7 @@ namespace detail { //============================================================================== void MaterialVolumes::add_volume( - int index_elem, int index_material, double volume) + int index_elem, int index_material, double volume, const BoundingBox* bbox) { // This method handles adding elements to the materials hash table, // implementing open addressing with linear probing. Consistency across @@ -166,10 +225,18 @@ void MaterialVolumes::add_volume( // Non-atomic read of current material int32_t current_val = *slot_ptr; - // Found the desired material; accumulate volume + // Found the desired material; accumulate volume and bbox if (current_val == index_material) { #pragma omp atomic this->volumes(index_elem, slot) += volume; + if (bbox) { + atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x); + atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y); + atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z); + atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x); + atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y); + atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z); + } return; } @@ -185,6 +252,14 @@ void MaterialVolumes::add_volume( if (claimed_slot || (expected_val == index_material)) { #pragma omp atomic this->volumes(index_elem, slot) += volume; + if (bbox) { + atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x); + atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y); + atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z); + atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x); + atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y); + atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z); + } return; } } @@ -195,7 +270,7 @@ void MaterialVolumes::add_volume( } void MaterialVolumes::add_volume_unsafe( - int index_elem, int index_material, double volume) + int index_elem, int index_material, double volume, const BoundingBox* bbox) { // Linear probe for (int attempt = 0; attempt < table_size_; ++attempt) { @@ -207,9 +282,23 @@ void MaterialVolumes::add_volume_unsafe( // Read current material int32_t current_val = this->materials(index_elem, slot); - // Found the desired material; accumulate volume + // Found the desired material; accumulate volume and bbox if (current_val == index_material) { this->volumes(index_elem, slot) += volume; + if (bbox) { + this->bboxes(index_elem, slot, 0) = + std::min(this->bboxes(index_elem, slot, 0), bbox->min.x); + this->bboxes(index_elem, slot, 1) = + std::min(this->bboxes(index_elem, slot, 1), bbox->min.y); + this->bboxes(index_elem, slot, 2) = + std::min(this->bboxes(index_elem, slot, 2), bbox->min.z); + this->bboxes(index_elem, slot, 3) = + std::max(this->bboxes(index_elem, slot, 3), bbox->max.x); + this->bboxes(index_elem, slot, 4) = + std::max(this->bboxes(index_elem, slot, 4), bbox->max.y); + this->bboxes(index_elem, slot, 5) = + std::max(this->bboxes(index_elem, slot, 5), bbox->max.z); + } return; } @@ -217,6 +306,20 @@ void MaterialVolumes::add_volume_unsafe( if (current_val == EMPTY) { this->materials(index_elem, slot) = index_material; this->volumes(index_elem, slot) += volume; + if (bbox) { + this->bboxes(index_elem, slot, 0) = + std::min(this->bboxes(index_elem, slot, 0), bbox->min.x); + this->bboxes(index_elem, slot, 1) = + std::min(this->bboxes(index_elem, slot, 1), bbox->min.y); + this->bboxes(index_elem, slot, 2) = + std::min(this->bboxes(index_elem, slot, 2), bbox->min.z); + this->bboxes(index_elem, slot, 3) = + std::max(this->bboxes(index_elem, slot, 3), bbox->max.x); + this->bboxes(index_elem, slot, 4) = + std::max(this->bboxes(index_elem, slot, 4), bbox->max.y); + this->bboxes(index_elem, slot, 5) = + std::max(this->bboxes(index_elem, slot, 5), bbox->max.z); + } return; } } @@ -333,6 +436,12 @@ vector Mesh::volumes() const void Mesh::material_volumes(int nx, int ny, int nz, int table_size, int32_t* materials, double* volumes) const +{ + this->material_volumes(nx, ny, nz, table_size, materials, volumes, nullptr); +} + +void Mesh::material_volumes(int nx, int ny, int nz, int table_size, + int32_t* materials, double* volumes, double* bboxes) const { if (mpi::master) { header("MESH MATERIAL VOLUMES CALCULATION", 7); @@ -351,7 +460,8 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, timer.start(); // Create object for keeping track of materials/volumes - detail::MaterialVolumes result(materials, volumes, table_size); + detail::MaterialVolumes result(materials, volumes, bboxes, table_size); + bool compute_bboxes = bboxes != nullptr; // Determine bounding box auto bbox = this->bounding_box(); @@ -370,13 +480,13 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, #pragma omp parallel { // Preallocate vector for mesh indices and length fractions and particle - std::vector bins; - std::vector length_fractions; + vector bins; + vector length_fractions; Particle p; SourceSite site; site.E = 1.0; - site.particle = ParticleType::neutron; + site.particle = ParticleType::neutron(); for (int axis = 0; axis < 3; ++axis) { // Set starting position and direction @@ -453,12 +563,36 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, if (i_material != C_NONE) { i_material = model::materials[i_material]->id(); } + double cumulative_frac = 0.0; for (int i_bin = 0; i_bin < bins.size(); i_bin++) { int mesh_index = bins[i_bin]; double length = distance * length_fractions[i_bin]; + double volume = length * d1 * d2; - // Add volume to result - result.add_volume(mesh_index, i_material, length * d1 * d2); + if (compute_bboxes) { + double axis_start = r0[axis] + distance * cumulative_frac; + double axis_end = axis_start + length; + cumulative_frac += length_fractions[i_bin]; + + Position contrib_min = site.r; + Position contrib_max = site.r; + + contrib_min[ax1] = site.r[ax1] - 0.5 * d1; + contrib_max[ax1] = site.r[ax1] + 0.5 * d1; + contrib_min[ax2] = site.r[ax2] - 0.5 * d2; + contrib_max[ax2] = site.r[ax2] + 0.5 * d2; + contrib_min[axis] = std::min(axis_start, axis_end); + contrib_max[axis] = std::max(axis_start, axis_end); + + BoundingBox contrib_bbox {contrib_min, contrib_max}; + contrib_bbox &= bbox; + + result.add_volume( + mesh_index, i_material, volume, &contrib_bbox); + } else { + // Add volume to result + result.add_volume(mesh_index, i_material, volume); + } } if (distance == max_distance) @@ -505,10 +639,15 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, // Combine results from multiple MPI processes if (mpi::n_procs > 1) { int total = this->n_bins() * table_size; + int total_bbox = total * 6; if (mpi::master) { // Allocate temporary buffer for receiving data - std::vector mats(total); - std::vector vols(total); + vector mats(total); + vector vols(total); + vector recv_bboxes; + if (compute_bboxes) { + recv_bboxes.resize(total_bbox); + } for (int i = 1; i < mpi::n_procs; ++i) { // Receive material indices and volumes from process i @@ -516,6 +655,10 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, MPI_STATUS_IGNORE); MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm, MPI_STATUS_IGNORE); + if (compute_bboxes) { + MPI_Recv(recv_bboxes.data(), total_bbox, MPI_DOUBLE, i, i, + mpi::intracomm, MPI_STATUS_IGNORE); + } // Combine with existing results; we can call thread unsafe version of // add_volume because each thread is operating on a different element @@ -524,7 +667,18 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, for (int k = 0; k < table_size; ++k) { int index = index_elem * table_size + k; if (mats[index] != EMPTY) { - result.add_volume_unsafe(index_elem, mats[index], vols[index]); + if (compute_bboxes) { + int bbox_index = index * 6; + BoundingBox slot_bbox { + {recv_bboxes[bbox_index + 0], recv_bboxes[bbox_index + 1], + recv_bboxes[bbox_index + 2]}, + {recv_bboxes[bbox_index + 3], recv_bboxes[bbox_index + 4], + recv_bboxes[bbox_index + 5]}}; + result.add_volume_unsafe( + index_elem, mats[index], vols[index], &slot_bbox); + } else { + result.add_volume_unsafe(index_elem, mats[index], vols[index]); + } } } } @@ -533,6 +687,9 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, // Send material indices and volumes to process 0 MPI_Send(materials, total, MPI_INT32_T, 0, mpi::rank, mpi::intracomm); MPI_Send(volumes, total, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm); + if (compute_bboxes) { + MPI_Send(bboxes, total_bbox, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm); + } } } @@ -2428,14 +2585,14 @@ extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur) } extern "C" int openmc_mesh_material_volumes(int32_t index, int nx, int ny, - int nz, int table_size, int32_t* materials, double* volumes) + int nz, int table_size, int32_t* materials, double* volumes, double* bboxes) { if (int err = check_mesh(index)) return err; try { model::meshes[index]->material_volumes( - nx, ny, nz, table_size, materials, volumes); + nx, ny, nz, table_size, materials, volumes, bboxes); } catch (const std::exception& e) { set_errmsg(e.what()); if (starts_with(e.what(), "Mesh")) { @@ -3485,9 +3642,6 @@ void LibMesh::initialize() // assuming that unstructured meshes used in OpenMC are 3D n_dimension_ = 3; - if (length_multiplier_ > 0.0) { - libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); - } // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh. // Otherwise assume that it is prepared by its owning application if (unique_m_) { @@ -3537,7 +3691,11 @@ Position LibMesh::centroid(int bin) const { const auto& elem = this->get_element_from_bin(bin); auto centroid = elem.vertex_average(); - return {centroid(0), centroid(1), centroid(2)}; + if (length_multiplier_ > 0.0) { + return length_multiplier_ * Position(centroid(0), centroid(1), centroid(2)); + } else { + return {centroid(0), centroid(1), centroid(2)}; + } } int LibMesh::n_vertices() const @@ -3548,7 +3706,11 @@ int LibMesh::n_vertices() const Position LibMesh::vertex(int vertex_id) const { const auto node_ref = m_->node_ref(vertex_id); - return {node_ref(0), node_ref(1), node_ref(2)}; + if (length_multiplier_ > 0.0) { + return length_multiplier_ * Position(node_ref(0), node_ref(1), node_ref(2)); + } else { + return {node_ref(0), node_ref(1), node_ref(2)}; + } } std::vector LibMesh::connectivity(int elem_id) const @@ -3689,6 +3851,11 @@ int LibMesh::get_bin(Position r) const // look-up a tet using the point locator libMesh::Point p(r.x, r.y, r.z); + if (length_multiplier_ > 0.0) { + // Scale the point down + p /= length_multiplier_; + } + // quick rejection check if (!bbox_.contains_point(p)) { return -1; @@ -3722,22 +3889,32 @@ const libMesh::Elem& LibMesh::get_element_from_bin(int bin) const double LibMesh::volume(int bin) const { - return this->get_element_from_bin(bin).volume(); + return this->get_element_from_bin(bin).volume() * length_multiplier_ * + length_multiplier_ * length_multiplier_; } -AdaptiveLibMesh::AdaptiveLibMesh( - libMesh::MeshBase& input_mesh, double length_multiplier) - : LibMesh(input_mesh, length_multiplier), num_active_(m_->n_active_elem()) +AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, + double length_multiplier, + const std::set& block_ids) + : LibMesh(input_mesh, length_multiplier), block_ids_(block_ids), + block_restrict_(!block_ids_.empty()), + num_active_( + block_restrict_ + ? std::distance(m_->active_subdomain_set_elements_begin(block_ids_), + m_->active_subdomain_set_elements_end(block_ids_)) + : m_->n_active_elem()) { // if the mesh is adaptive elements aren't guaranteed by libMesh to be // contiguous in ID space, so we need to map from bin indices (defined over // active elements) to global dof ids bin_to_elem_map_.reserve(num_active_); elem_to_bin_map_.resize(m_->n_elem(), -1); - for (auto it = m_->active_elements_begin(); it != m_->active_elements_end(); - it++) { - auto elem = *it; - + auto begin = block_restrict_ + ? m_->active_subdomain_set_elements_begin(block_ids_) + : m_->active_elements_begin(); + auto end = block_restrict_ ? m_->active_subdomain_set_elements_end(block_ids_) + : m_->active_elements_end(); + for (const auto& elem : libMesh::as_range(begin, end)) { bin_to_elem_map_.push_back(elem->id()); elem_to_bin_map_[elem->id()] = bin_to_elem_map_.size() - 1; } @@ -3770,6 +3947,27 @@ void AdaptiveLibMesh::write(const std::string& filename) const this->id_)); } +int AdaptiveLibMesh::get_bin(Position r) const +{ + // look-up a tet using the point locator + libMesh::Point p(r.x, r.y, r.z); + + if (length_multiplier_ > 0.0) { + // Scale the point down + p /= length_multiplier_; + } + + // quick rejection check + if (!bbox_.contains_point(p)) { + return -1; + } + + const auto& point_locator = pl_.at(thread_num()); + + const auto elem_ptr = (*point_locator)(p, &block_ids_); + return elem_ptr ? get_bin_from_element(elem_ptr) : -1; +} + int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const { int bin = elem_to_bin_map_[elem->id()]; diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index ed734d401e..34f87d1798 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -244,7 +244,7 @@ void MgxsInterface::read_header(const std::string& path_cross_sections) void put_mgxs_header_data_to_globals() { // Get the minimum and maximum energies - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); data::energy_min[neutron] = data::mg.energy_bins_.back(); data::energy_max[neutron] = data::mg.energy_bins_.front(); diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 5ae6e30ee2..69e603a7c6 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -379,7 +379,7 @@ void Nuclide::create_derived( auto pprod = xt::view(xs_[t], xt::range(j, j + n), XS_PHOTON_PROD); for (const auto& p : rx->products_) { - if (p.particle_ == ParticleType::photon) { + if (p.particle_.is_photon()) { for (int k = 0; k < n; ++k) { double E = grid_[t].energy[k + j]; @@ -501,7 +501,7 @@ void Nuclide::create_derived( void Nuclide::init_grid() { - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); double E_min = data::energy_min[neutron]; double E_max = data::energy_max[neutron]; int M = settings::n_log_bins; @@ -552,7 +552,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const for (int i = 1; i < rx->products_.size(); ++i) { // Skip any non-neutron products const auto& product = rx->products_[i]; - if (product.particle_ != ParticleType::neutron) + if (!product.particle_.is_neutron()) continue; // Evaluate yield diff --git a/src/output.cpp b/src/output.cpp index 80e2b10ab8..ae2daaffc1 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -155,21 +155,21 @@ std::string time_stamp() void print_particle(Particle& p) { // Display particle type and ID. - switch (p.type()) { - case ParticleType::neutron: + switch (p.type().pdg_number()) { + case PDG_NEUTRON: fmt::print("Neutron "); break; - case ParticleType::photon: + case PDG_PHOTON: fmt::print("Photon "); break; - case ParticleType::electron: + case PDG_ELECTRON: fmt::print("Electron "); break; - case ParticleType::positron: + case PDG_POSITRON: fmt::print("Positron "); break; default: - fmt::print("Unknown Particle "); + fmt::print("Particle {} ", p.type().str()); } fmt::print("{}\n", p.id()); diff --git a/src/particle.cpp b/src/particle.cpp index b9f9c8f86b..a1176abc79 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -48,17 +48,19 @@ double Particle::speed() const if (settings::run_CE) { // Determine mass in eV/c^2 double mass; - switch (this->type()) { - case ParticleType::neutron: + switch (this->type().pdg_number()) { + case PDG_NEUTRON: mass = MASS_NEUTRON_EV; break; - case ParticleType::photon: + case PDG_PHOTON: mass = 0.0; break; - case ParticleType::electron: - case ParticleType::positron: + case PDG_ELECTRON: + case PDG_POSITRON: mass = MASS_ELECTRON_EV; break; + default: + fatal_error("Unsupported particle for speed calculation."); } // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<E() * (this->E() + 2 * mass)) / @@ -77,7 +79,11 @@ bool Particle::create_secondary( { // If energy is below cutoff for this particle, don't create secondary // particle - if (E < settings::energy_cutoff[static_cast(type)]) { + int idx = type.transport_index(); + if (idx == C_NONE) { + return false; + } + if (E < settings::energy_cutoff[idx]) { return false; } @@ -235,7 +241,8 @@ void Particle::event_advance() boundary() = distance_to_boundary(*this); // Sample a distance to collision - if (type() == ParticleType::electron || type() == ParticleType::positron) { + if (type() == ParticleType::electron() || + type() == ParticleType::positron()) { collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0; } else if (macro_xs().total == 0.0) { collision_distance() = INFINITY; @@ -244,7 +251,7 @@ void Particle::event_advance() } double speed = this->speed(); - double time_cutoff = settings::time_cutoff[static_cast(type())]; + double time_cutoff = settings::time_cutoff[type().transport_index()]; double distance_cutoff = (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY; @@ -269,8 +276,7 @@ void Particle::event_advance() } // Score track-length estimate of k-eff - if (settings::run_mode == RunMode::EIGENVALUE && - type() == ParticleType::neutron) { + if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) { keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission; } @@ -331,8 +337,7 @@ void Particle::event_cross_surface() void Particle::event_collide() { // Score collision estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE && - type() == ParticleType::neutron) { + if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) { keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total; } @@ -370,8 +375,7 @@ void Particle::event_collide() } } - if (!model::active_pulse_height_tallies.empty() && - type() == ParticleType::photon) { + if (!model::active_pulse_height_tallies.empty() && type().is_photon()) { pht_collision_energy(); } @@ -442,7 +446,7 @@ void Particle::event_revive_from_secondary() // Subtract secondary particle energy from interim pulse-height results if (!model::active_pulse_height_tallies.empty() && - this->type() == ParticleType::photon) { + this->type().is_photon()) { // Since the birth cell of the particle has not been set we // have to determine it before the energy of the secondary particle can be // removed from the pulse-height of this cell. @@ -525,7 +529,7 @@ void Particle::pht_collision_energy() // If the energy of the particle is below the cutoff, it will not be sampled // so its energy is added to the pulse-height in the cell - int photon = static_cast(ParticleType::photon); + int photon = ParticleType::photon().transport_index(); if (E() < settings::energy_cutoff[photon]) { pht_storage()[index] += E(); } @@ -824,7 +828,7 @@ void Particle::write_restart() const break; } write_dataset(file_id, "id", id()); - write_dataset(file_id, "type", static_cast(type())); + write_dataset(file_id, "type", type().pdg_number()); int64_t i = current_work(); if (settings::run_mode == RunMode::EIGENVALUE) { @@ -878,37 +882,6 @@ void Particle::update_neutron_xs( //============================================================================== // Non-method functions //============================================================================== - -std::string particle_type_to_str(ParticleType type) -{ - switch (type) { - case ParticleType::neutron: - return "neutron"; - case ParticleType::photon: - return "photon"; - case ParticleType::electron: - return "electron"; - case ParticleType::positron: - return "positron"; - } - UNREACHABLE(); -} - -ParticleType str_to_particle_type(std::string str) -{ - if (str == "neutron") { - return ParticleType::neutron; - } else if (str == "photon") { - return ParticleType::photon; - } else if (str == "electron") { - return ParticleType::electron; - } else if (str == "positron") { - return ParticleType::positron; - } else { - throw std::invalid_argument {fmt::format("Invalid particle name: {}", str)}; - } -} - void add_surf_source_to_bank(Particle& p, const Surface& surf) { if (simulation::current_batch <= settings::n_inactive || diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 6b7778211a..f02fcb94b5 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -31,6 +31,16 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) hid_t file_id = file_open(settings::path_particle_restart, 'r'); // Read data from file + bool legacy_particle_codes = true; + if (attribute_exists(file_id, "version")) { + array version; + read_attribute(file_id, "version", version); + if (version[0] > VERSION_PARTICLE_RESTART[0] || + (version[0] == VERSION_PARTICLE_RESTART[0] && version[1] >= 1)) { + legacy_particle_codes = false; + } + } + read_dataset(file_id, "current_batch", simulation::current_batch); read_dataset(file_id, "generations_per_batch", settings::gen_per_batch); read_dataset(file_id, "current_generation", simulation::current_gen); @@ -45,7 +55,8 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) read_dataset(file_id, "id", p.id()); int type; read_dataset(file_id, "type", type); - p.type() = static_cast(type); + p.type() = legacy_particle_codes ? legacy_particle_index_to_type(type) + : ParticleType {type}; read_dataset(file_id, "weight", p.wgt()); read_dataset(file_id, "energy", p.E()); read_dataset(file_id, "xyz", p.r()); diff --git a/src/particle_type.cpp b/src/particle_type.cpp new file mode 100644 index 0000000000..fe8f9fe2a0 --- /dev/null +++ b/src/particle_type.cpp @@ -0,0 +1,246 @@ +#include "openmc/particle_type.h" + +#include +#include +#include + +#include "openmc/string_utils.h" + +namespace openmc { +namespace { + +constexpr const char* ATOMIC_SYMBOL[] = {"", "H", "He", "Li", "Be", "B", "C", + "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", + "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", + "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", + "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", + "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", + "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", + "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", + "Es", "Fm", "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", + "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"}; + +constexpr int MAX_Z = + static_cast(sizeof(ATOMIC_SYMBOL) / sizeof(ATOMIC_SYMBOL[0])) - 1; + +bool is_integer_string(const std::string& s) +{ + if (s.empty()) + return false; + size_t i = 0; + if (s[0] == '-' || s[0] == '+') { + if (s.size() == 1) + return false; + i = 1; + } + for (; i < s.size(); ++i) { + if (!std::isdigit(static_cast(s[i]))) + return false; + } + return true; +} + +int atomic_number_from_symbol(std::string_view symbol) +{ + for (int z = 1; z <= MAX_Z; ++z) { + if (symbol == ATOMIC_SYMBOL[z]) { + return z; + } + } + return 0; +} + +bool parse_gnds_nuclide(std::string_view name, int& Z, int& A, int& m) +{ + if (name.empty()) + return false; + + size_t pos = 0; + if (!std::isupper(static_cast(name[pos]))) + return false; + + std::string symbol; + symbol += name[pos++]; + if (pos < name.size() && + std::islower(static_cast(name[pos]))) { + symbol += name[pos++]; + } + + if (pos >= name.size() || + !std::isdigit(static_cast(name[pos]))) { + return false; + } + + size_t a_start = pos; + while ( + pos < name.size() && std::isdigit(static_cast(name[pos]))) { + ++pos; + } + A = std::stoi(std::string {name.substr(a_start, pos - a_start)}); + if (A <= 0 || A > 999) + return false; + + m = 0; + if (pos < name.size()) { + if (name[pos] != '_' || pos + 2 >= name.size() || name[pos + 1] != 'm') { + return false; + } + pos += 2; + size_t m_start = pos; + while (pos < name.size() && + std::isdigit(static_cast(name[pos]))) { + ++pos; + } + if (m_start == pos) + return false; + m = std::stoi(std::string {name.substr(m_start, pos - m_start)}); + if (m < 0 || m > 9) + return false; + } + + if (pos != name.size()) + return false; + + Z = atomic_number_from_symbol(symbol); + return Z != 0; +} + +// Helper to convert nuclear PDG number to nuclide name +std::string nuclide_name_from_pdg(int32_t pdg) +{ + int32_t code = pdg; + int m = code % 10; + int A = (code / 10) % 1000; + int Z = (code / 10000) % 1000; + + if (Z <= 0 || Z > MAX_Z || A <= 0 || A > 999) { + throw std::invalid_argument { + "Invalid nuclear PDG number: " + std::to_string(pdg)}; + } + + std::string name = ATOMIC_SYMBOL[Z] + std::to_string(A); + if (m > 0) { + name += "_m" + std::to_string(m); + } + return name; +} + +} // namespace + +//============================================================================== +// ParticleType member function implementations +//============================================================================== + +ParticleType::ParticleType(std::string_view str) +{ + std::string s {str}; + strtrim(s); + if (s.empty()) { + throw std::invalid_argument {"Particle string is empty."}; + } + + std::string lower = s; + to_lower(lower); + + // Check for pdg: prefix + if (starts_with(lower, "pdg:")) { + std::string value_str = lower.substr(4); + if (!is_integer_string(value_str)) { + throw std::invalid_argument {"Invalid PDG number: " + value_str}; + } + pdg_number_ = std::stoi(value_str); + return; + } + + // Check for known particle names + if (lower == "neutron" || lower == "n") { + pdg_number_ = PDG_NEUTRON; + return; + } + if (lower == "photon" || lower == "gamma") { + pdg_number_ = PDG_PHOTON; + return; + } + if (lower == "electron") { + pdg_number_ = PDG_ELECTRON; + return; + } + if (lower == "positron") { + pdg_number_ = PDG_POSITRON; + return; + } + if (lower == "proton" || lower == "p" || lower == "h1") { + pdg_number_ = PDG_PROTON; + return; + } + if (lower == "deuteron" || lower == "d" || lower == "h2") { + pdg_number_ = PDG_DEUTERON; + return; + } + if (lower == "triton" || lower == "t" || lower == "h3") { + pdg_number_ = PDG_TRITON; + return; + } + if (lower == "alpha" || lower == "he4") { + pdg_number_ = PDG_ALPHA; + return; + } + + // Check for integer string + if (is_integer_string(s)) { + pdg_number_ = std::stoi(s); + return; + } + + // Try to parse as GNDS nuclide name + int Z = 0; + int A = 0; + int m = 0; + if (!parse_gnds_nuclide(s, Z, A, m)) { + throw std::invalid_argument {"Invalid nuclide name: " + s}; + } + pdg_number_ = 1000000000 + Z * 10000 + A * 10 + m; +} + +std::string ParticleType::str() const +{ + if (pdg_number_ == PDG_NEUTRON) + return "neutron"; + if (pdg_number_ == PDG_PHOTON) + return "photon"; + if (pdg_number_ == PDG_ELECTRON) + return "electron"; + if (pdg_number_ == PDG_POSITRON) + return "positron"; + if (pdg_number_ == PDG_PROTON) + return "proton"; + + if (is_nucleus()) { + return nuclide_name_from_pdg(pdg_number_); + } + + return "pdg:" + std::to_string(pdg_number_); +} + +//============================================================================== +// Free function implementations +//============================================================================== + +ParticleType legacy_particle_index_to_type(int index) +{ + switch (index) { + case 0: + return ParticleType {PDG_NEUTRON}; + case 1: + return ParticleType {PDG_PHOTON}; + case 2: + return ParticleType {PDG_ELECTRON}; + case 3: + return ParticleType {PDG_POSITRON}; + default: + throw std::invalid_argument { + "Invalid legacy particle index: " + std::to_string(index)}; + } +} + +} // namespace openmc diff --git a/src/photon.cpp b/src/photon.cpp index 4926e3eae6..951acb9fbd 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -293,7 +293,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_group(rgroup); // Truncate the bremsstrahlung data at the cutoff energy - int photon = static_cast(ParticleType::photon); + int photon = ParticleType::photon().transport_index(); const auto& E {electron_energy}; double cutoff = settings::energy_cutoff[photon]; if (cutoff > E(0)) { @@ -805,7 +805,7 @@ void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const if (shell.transitions.empty()) { Direction u = isotropic_direction(p.current_seed()); double E = shell.binding_energy; - p.create_secondary(p.wgt(), u, E, ParticleType::photon); + p.create_secondary(p.wgt(), u, E, ParticleType::photon()); continue; } @@ -833,12 +833,13 @@ void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const holes[n_holes++] = transition.secondary_subshell; // Create auger electron - p.create_secondary(p.wgt(), u, transition.energy, ParticleType::electron); + p.create_secondary( + p.wgt(), u, transition.energy, ParticleType::electron()); } else { // Radiative transition -- get X-ray energy // Create fluorescent photon - p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon); + p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon()); } } } diff --git a/src/physics.cpp b/src/physics.cpp index d03c85195e..6bdd3f2857 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -46,27 +46,29 @@ void collision(Particle& p) ++(p.n_collision()); // Sample reaction for the material the particle is in - switch (p.type()) { - case ParticleType::neutron: + switch (p.type().pdg_number()) { + case PDG_NEUTRON: sample_neutron_reaction(p); break; - case ParticleType::photon: + case PDG_PHOTON: sample_photon_reaction(p); break; - case ParticleType::electron: + case PDG_ELECTRON: sample_electron_reaction(p); break; - case ParticleType::positron: + case PDG_POSITRON: sample_positron_reaction(p); break; + default: + fatal_error("Unsupported particle PDG for collision sampling."); } if (settings::weight_window_checkpoint_collision) apply_weight_windows(p); // Kill particle if energy falls below cutoff - int type = static_cast(p.type()); - if (p.E() < settings::energy_cutoff[type]) { + int type = p.type().transport_index(); + if (type != C_NONE && p.E() < settings::energy_cutoff[type]) { p.wgt() = 0.0; } @@ -75,11 +77,11 @@ void collision(Particle& p) std::string msg; if (p.event() == TallyEvent::KILL) { msg = fmt::format(" Killed. Energy = {} eV.", p.E()); - } else if (p.type() == ParticleType::neutron) { + } else if (p.type().is_neutron()) { msg = fmt::format(" {} with {}. Energy = {} eV.", reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_, p.E()); - } else if (p.type() == ParticleType::photon) { + } else if (p.type().is_photon()) { msg = fmt::format(" {} with {}. Energy = {} eV.", reaction_name(p.event_mt()), to_element(data::nuclides[p.event_nuclide()]->name_), p.E()); @@ -208,7 +210,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Initialize fission site object with particle data SourceSite site; site.r = p.r(); - site.particle = ParticleType::neutron; + site.particle = ParticleType::neutron(); site.time = p.time(); site.wgt = 1. / weight; site.surf_id = 0; @@ -218,7 +220,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Reject site if it exceeds time cutoff if (site.delayed_group > 0) { - double t_cutoff = settings::time_cutoff[static_cast(site.particle)]; + double t_cutoff = settings::time_cutoff[site.particle.transport_index()]; if (site.time > t_cutoff) { continue; } @@ -287,7 +289,7 @@ void sample_photon_reaction(Particle& p) // Kill photon if below energy cutoff -- an extra check is made here because // photons with energy below the cutoff may have been produced by neutrons // reactions or atomic relaxation - int photon = static_cast(ParticleType::photon); + int photon = ParticleType::photon().transport_index(); if (p.E() < settings::energy_cutoff[photon]) { p.E() = 0.0; p.wgt() = 0.0; @@ -337,13 +339,13 @@ void sample_photon_reaction(Particle& p) // Create Compton electron double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b; - int electron = static_cast(ParticleType::electron); + int electron = ParticleType::electron().transport_index(); if (E_electron >= settings::energy_cutoff[electron]) { double mu_electron = (alpha - alpha_out * p.mu()) / std::sqrt(alpha * alpha + alpha_out * alpha_out - 2.0 * alpha * alpha_out * p.mu()); Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron()); } // Allow electrons to fill orbital and produce Auger electrons and @@ -416,7 +418,7 @@ void sample_photon_reaction(Particle& p) u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi); // Create secondary electron - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron()); // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons @@ -441,12 +443,11 @@ void sample_photon_reaction(Particle& p) // Create secondary electron Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron()); // Create secondary positron u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron); - + p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron()); p.event() = TallyEvent::ABSORB; p.event_mt() = PAIR_PROD; p.wgt() = 0.0; @@ -481,8 +482,8 @@ void sample_positron_reaction(Particle& p) Direction u = isotropic_direction(p.current_seed()); // Create annihilation photon pair traveling in opposite directions - p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon); - p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon); + p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon()); + p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon()); p.E() = 0.0; p.wgt() = 0.0; @@ -607,7 +608,7 @@ void sample_photon_product( continue; for (int j = 0; j < rx->products_.size(); ++j) { - if (rx->products_[j].particle_ == ParticleType::photon) { + if (rx->products_[j].particle_.is_photon()) { // For fission, artificially increase the photon yield to account // for delayed photons double f = 1.0; @@ -1096,7 +1097,7 @@ void sample_fission_neutron( rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed); // resample if energy is greater than maximum neutron energy - constexpr int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); if (site->E < data::energy_max[neutron]) break; @@ -1156,7 +1157,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) if (std::floor(yield) == yield && yield > 0) { // If yield is integral, create exactly that many secondary particles for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { - p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron); + p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron()); } } else { // Otherwise, change weight of particle based on yield @@ -1212,7 +1213,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide) } // Create the secondary photon - bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon); + bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon()); // Tag secondary particle with parent nuclide if (created_photon && settings::use_decay_photons) { diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 4c28cb1795..866d4d728e 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -136,7 +136,7 @@ void create_fission_sites(Particle& p) // Initialize fission site object with particle data SourceSite site; site.r = p.r(); - site.particle = ParticleType::neutron; + site.particle = ParticleType::neutron(); site.time = p.time(); site.wgt = 1. / weight; @@ -171,7 +171,7 @@ void create_fission_sites(Particle& p) site.time -= std::log(prn(p.current_seed())) / decay_rate; // Reject site if it exceeds time cutoff - double t_cutoff = settings::time_cutoff[static_cast(site.particle)]; + double t_cutoff = settings::time_cutoff[site.particle.transport_index()]; if (site.time > t_cutoff) { continue; } diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index ec14795dd2..c2effaa5d4 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -109,18 +109,21 @@ void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) // Add scattering + fission source int material = srh.material(); + double density_mult = srh.density_mult(); if (material != MATERIAL_VOID) { double inverse_k_eff = 1.0 / k_eff_; for (int g_out = 0; g_out < negroups_; g_out++) { - double sigma_t = sigma_t_[material * negroups_ + g_out]; + double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult; double scatter_source = 0.0; double fission_source = 0.0; for (int g_in = 0; g_in < negroups_; g_in++) { double scalar_flux = srh.scalar_flux_old(g_in); - double sigma_s = - sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; - double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; + double sigma_s = sigma_s_[material * negroups_ * negroups_ + + g_out * negroups_ + g_in] * + density_mult; + double nu_sigma_f = + nu_sigma_f_[material * negroups_ + g_in] * density_mult; double chi = chi_[material * negroups_ + g_out]; scatter_source += sigma_s * scalar_flux; @@ -198,7 +201,8 @@ void FlatSourceDomain::set_flux_to_flux_plus_source( source_regions_.volume_sq(sr); } } else { - double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g]; + double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g] * + source_regions_.density_mult(sr); source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume); source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g); } @@ -332,7 +336,8 @@ void FlatSourceDomain::compute_k_eff() double sr_fission_source_new = 0; for (int g = 0; g < negroups_; g++) { - double nu_sigma_f = nu_sigma_f_[material * negroups_ + g]; + double nu_sigma_f = nu_sigma_f_[material * negroups_ + g] * + source_regions_.density_mult(sr); sr_fission_source_old += nu_sigma_f * source_regions_.scalar_flux_old(sr, g); sr_fission_source_new += @@ -562,7 +567,8 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const // to get the total source strength in the expected units. double sigma_t = 1.0; if (material != MATERIAL_VOID) { - sigma_t = sigma_t_[material * negroups_ + g]; + sigma_t = + sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr); } simulation_external_source_strength += source_regions_.external_source(sr, g) * sigma_t * volume; @@ -618,7 +624,9 @@ void FlatSourceDomain::random_ray_tally() // source strength. double volume = source_regions_.volume(sr) * simulation_volume_; - double material = source_regions_.material(sr); + int material = source_regions_.material(sr); + double density_mult = source_regions_.density_mult(sr); + for (int g = 0; g < negroups_; g++) { double flux = source_regions_.scalar_flux_new(sr, g) * source_normalization_factor; @@ -634,19 +642,22 @@ void FlatSourceDomain::random_ray_tally() case SCORE_TOTAL: if (material != MATERIAL_VOID) { - score = flux * volume * sigma_t_[material * negroups_ + g]; + score = + flux * volume * sigma_t_[material * negroups_ + g] * density_mult; } break; case SCORE_FISSION: if (material != MATERIAL_VOID) { - score = flux * volume * sigma_f_[material * negroups_ + g]; + score = + flux * volume * sigma_f_[material * negroups_ + g] * density_mult; } break; case SCORE_NU_FISSION: if (material != MATERIAL_VOID) { - score = flux * volume * nu_sigma_f_[material * negroups_ + g]; + score = flux * volume * nu_sigma_f_[material * negroups_ + g] * + density_mult; } break; @@ -654,10 +665,15 @@ void FlatSourceDomain::random_ray_tally() score = 1.0; break; + case SCORE_KAPPA_FISSION: + score = flux * volume * kappa_fission_[material * negroups_ + g] * + density_mult; + break; + default: fatal_error("Invalid score specified in tallies.xml. Only flux, " - "total, fission, nu-fission, and events are supported in " - "random ray mode."); + "total, fission, nu-fission, kappa-fission, and events " + "are supported in random ray mode."); break; } // Apply score to the appropriate tally bin @@ -913,7 +929,8 @@ void FlatSourceDomain::output_to_vtk() const for (int g = 0; g < negroups_; g++) { int64_t source_element = fsr * negroups_ + g; float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g); - double sigma_f = sigma_f_[mat * negroups_ + g]; + double sigma_f = sigma_f_[mat * negroups_ + g] * + source_regions_.density_mult(fsr); total_fission += sigma_f * flux; } } @@ -934,7 +951,8 @@ void FlatSourceDomain::output_to_vtk() const // multiply it back to get the true external source. double sigma_t = 1.0; if (mat != MATERIAL_VOID) { - sigma_t = sigma_t_[mat * negroups_ + g]; + sigma_t = sigma_t_[mat * negroups_ + g] * + source_regions_.density_mult(fsr); } total_external += source_regions_.external_source(fsr, g) * sigma_t; } @@ -1152,6 +1170,10 @@ void FlatSourceDomain::flatten_xs() } chi_.push_back(chi); + double kappa_fission = + m.get_xs(MgxsType::KAPPA_FISSION, g_out, NULL, NULL, NULL, t, a); + kappa_fission_.push_back(kappa_fission); + for (int g_in = 0; g_in < negroups_; g_in++) { double sigma_s = m.get_xs(MgxsType::NU_SCATTER, g_in, &g_out, NULL, NULL, t, a); @@ -1167,6 +1189,7 @@ void FlatSourceDomain::flatten_xs() nu_sigma_f_.push_back(0); sigma_f_.push_back(0); chi_.push_back(0); + kappa_fission_.push_back(0); for (int g_in = 0; g_in < negroups_; g_in++) { sigma_s_.push_back(0); } @@ -1244,7 +1267,8 @@ void FlatSourceDomain::set_adjoint_sources() continue; } for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; + double sigma_t = + sigma_t_[material * negroups_ + g] * source_regions_.density_mult(sr); source_regions_.external_source(sr, g) /= sigma_t; } } @@ -1495,6 +1519,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( handle.material() = material; + handle.density_mult() = cell.density_mult(gs.cell_instance()); + // Store the mesh index (if any) assigned to this source region handle.mesh() = mesh_idx; @@ -1523,7 +1549,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // Divide external source term by sigma_t if (material != C_NONE) { for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; + double sigma_t = + sigma_t_[material * negroups_ + g] * handle.density_mult(); handle.external_source(g) /= sigma_t; } } @@ -1598,6 +1625,7 @@ void FlatSourceDomain::apply_transport_stabilization() #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions(); sr++) { int material = source_regions_.material(sr); + double density_mult = source_regions_.density_mult(sr); if (material == MATERIAL_VOID) { continue; } @@ -1605,9 +1633,10 @@ void FlatSourceDomain::apply_transport_stabilization() // Only apply stabilization if the diagonal (in-group) scattering XS is // negative double sigma_s = - sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g]; + sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g] * + density_mult; if (sigma_s < 0.0) { - double sigma_t = sigma_t_[material * negroups_ + g]; + double sigma_t = sigma_t_[material * negroups_ + g] * density_mult; double phi_new = source_regions_.scalar_flux_new(sr, g); double phi_old = source_regions_.scalar_flux_old(sr, g); diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index 47ffbb727b..02f4c9e235 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -43,12 +43,13 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) // Add scattering + fission source int material = srh.material(); + double density_mult = srh.density_mult(); if (material != MATERIAL_VOID) { double inverse_k_eff = 1.0 / k_eff_; MomentMatrix invM = srh.mom_matrix().inverse(); for (int g_out = 0; g_out < negroups_; g_out++) { - double sigma_t = sigma_t_[material * negroups_ + g_out]; + double sigma_t = sigma_t_[material * negroups_ + g_out] * density_mult; double scatter_flat = 0.0f; double fission_flat = 0.0f; @@ -61,9 +62,11 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) MomentArray flux_linear = srh.flux_moments_old(g_in); // Handles for cross sections - double sigma_s = - sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; - double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; + double sigma_s = sigma_s_[material * negroups_ * negroups_ + + g_out * negroups_ + g_in] * + density_mult; + double nu_sigma_f = + nu_sigma_f_[material * negroups_ + g_in] * density_mult; double chi = chi_[material * negroups_ + g_out]; // Compute source terms for flat and linear components of the flux diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index c19d136a4a..1b61d8c207 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -435,7 +435,8 @@ void RandomRay::attenuate_flux_flat_source( // MOC incoming flux attenuation + source contribution/attenuation equation for (int g = 0; g < negroups_; g++) { - float sigma_t = domain_->sigma_t_[material * negroups_ + g]; + float sigma_t = + domain_->sigma_t_[material * negroups_ + g] * srh.density_mult(); float tau = sigma_t * distance; float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau) float new_delta_psi = (angular_flux_[g] - srh.source(g)) * exponential; @@ -558,7 +559,8 @@ void RandomRay::attenuate_flux_linear_source( for (int g = 0; g < negroups_; g++) { // Compute tau, the optical thickness of the ray segment - float sigma_t = domain_->sigma_t_[material * negroups_ + g]; + float sigma_t = + domain_->sigma_t_[material * negroups_ + g] * srh.density_mult(); float tau = sigma_t * distance; // If tau is very small, set it to zero to avoid numerical issues. diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index d475b2593e..7bab3a9b1b 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -28,17 +28,13 @@ void openmc_run_random_ray() // Run forward simulation ////////////////////////////////////////////////////////// - // Check if adjoint calculation is needed. If it is, we will run the forward - // calculation first and then the adjoint calculation later. - bool adjoint_needed = FlatSourceDomain::adjoint_; - - // Configure the domain for forward simulation - FlatSourceDomain::adjoint_ = false; - - // If we're going to do an adjoint simulation afterwards, report that this is - // the initial forward flux solve. - if (adjoint_needed && mpi::master) - header("FORWARD FLUX SOLVE", 3); + if (mpi::master) { + if (FlatSourceDomain::adjoint_) { + FlatSourceDomain::adjoint_ = false; + openmc::print_adjoint_header(); + FlatSourceDomain::adjoint_ = true; + } + } // Initialize OpenMC general data structures openmc_simulation_init(); @@ -53,76 +49,20 @@ void openmc_run_random_ray() // Initialize fixed sources, if present sim.apply_fixed_sources_and_mesh_domains(); - // Begin main simulation timer - simulation::time_total.start(); - - // Execute random ray simulation + // Run initial random ray simulation sim.simulate(); - // End main simulation timer - simulation::time_total.stop(); - - // Normalize and save the final forward flux - double source_normalization_factor = - sim.domain()->compute_fixed_source_normalization_factor() / - (settings::n_batches - settings::n_inactive); - -#pragma omp parallel for - for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) { - sim.domain()->source_regions_.scalar_flux_final(se) *= - source_normalization_factor; - } - - // Finalize OpenMC - openmc_simulation_finalize(); - - // Output all simulation results - sim.output_simulation_results(); - ////////////////////////////////////////////////////////// // Run adjoint simulation (if enabled) ////////////////////////////////////////////////////////// - if (!adjoint_needed) { - return; + if (sim.adjoint_needed_) { + // Setup for adjoint simulation + sim.prepare_adjoint_simulation(); + + // Run adjoint simulation + sim.simulate(); } - - reset_timers(); - - // Configure the domain for adjoint simulation - FlatSourceDomain::adjoint_ = true; - - if (mpi::master) - header("ADJOINT FLUX SOLVE", 3); - - // Initialize OpenMC general data structures - openmc_simulation_init(); - - sim.domain()->k_eff_ = 1.0; - - // Initialize adjoint fixed sources, if present - sim.prepare_fixed_sources_adjoint(); - - // Transpose scattering matrix - sim.domain()->transpose_scattering_matrix(); - - // Swap nu_sigma_f and chi - sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_); - - // Begin main simulation timer - simulation::time_total.start(); - - // Execute random ray simulation - sim.simulate(); - - // End main simulation timer - simulation::time_total.stop(); - - // Finalize OpenMC - openmc_simulation_finalize(); - - // Output all simulation results - sim.output_simulation_results(); } // Enforces restrictions on inputs in random ray mode. While there are @@ -143,11 +83,12 @@ void validate_random_ray_inputs() case SCORE_FISSION: case SCORE_NU_FISSION: case SCORE_EVENTS: + case SCORE_KAPPA_FISSION: break; default: fatal_error( - "Invalid score specified. Only flux, total, fission, nu-fission, and " - "event scores are supported in random ray mode."); + "Invalid score specified. Only flux, total, fission, nu-fission, " + "kappa-fission, and event scores are supported in random ray mode."); } } @@ -353,6 +294,17 @@ void openmc_reset_random_ray() RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } +void print_adjoint_header() +{ + if (!FlatSourceDomain::adjoint_) + // If we're going to do an adjoint simulation afterwards, report that this + // is the initial forward flux solve. + header("FORWARD FLUX SOLVE", 3); + else + // Otherwise report that we are doing the adjoint simulation + header("ADJOINT FLUX SOLVE", 3); +} + //============================================================================== // RandomRaySimulation implementation //============================================================================== @@ -383,6 +335,16 @@ RandomRaySimulation::RandomRaySimulation() // Convert OpenMC native MGXS into a more efficient format // internal to the random ray solver domain_->flatten_xs(); + + // Check if adjoint calculation is needed. If it is, we will run the forward + // calculation first and then the adjoint calculation later. + adjoint_needed_ = FlatSourceDomain::adjoint_; + + // Adjoint is always false for the forward calculation + FlatSourceDomain::adjoint_ = false; + + // The first simulation is run after initialization + is_first_simulation_ = true; } void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() @@ -403,8 +365,41 @@ void RandomRaySimulation::prepare_fixed_sources_adjoint() } } +void RandomRaySimulation::prepare_adjoint_simulation() +{ + // Configure the domain for adjoint simulation + FlatSourceDomain::adjoint_ = true; + + // Reset k-eff + domain_->k_eff_ = 1.0; + + // Initialize adjoint fixed sources, if present + prepare_fixed_sources_adjoint(); + + // Transpose scattering matrix + domain_->transpose_scattering_matrix(); + + // Swap nu_sigma_f and chi + domain_->nu_sigma_f_.swap(domain_->chi_); +} + void RandomRaySimulation::simulate() { + if (!is_first_simulation_) { + if (mpi::master && adjoint_needed_) + openmc::print_adjoint_header(); + + // Reset the timers and reinitialize the general OpenMC datastructures if + // this is after the first simulation + reset_timers(); + + // Initialize OpenMC general data structures + openmc_simulation_init(); + } + + // Begin main simulation timer + simulation::time_total.start(); + // Random ray power iteration loop while (simulation::current_batch < settings::n_batches) { // Initialize the current batch @@ -493,6 +488,31 @@ void RandomRaySimulation::simulate() } // End random ray power iteration loop domain_->count_external_source_regions(); + + // End main simulation timer + simulation::time_total.stop(); + + // Normalize and save the final flux + double source_normalization_factor = + domain_->compute_fixed_source_normalization_factor() / + (settings::n_batches - settings::n_inactive); + +#pragma omp parallel for + for (uint64_t se = 0; se < domain_->n_source_elements(); se++) { + domain_->source_regions_.scalar_flux_final(se) *= + source_normalization_factor; + } + + // Finalize OpenMC + openmc_simulation_finalize(); + + // Output all simulation results + output_simulation_results(); + + // Toggle that the simulation object has been initialized after the first + // simulation + if (is_first_simulation_) + is_first_simulation_ = false; } void RandomRaySimulation::output_simulation_results() const diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 3b06f0ed09..15c65221aa 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -11,10 +11,11 @@ namespace openmc { //============================================================================== SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) : negroups_(sr.scalar_flux_old_.size()), material_(&sr.material_), - is_small_(&sr.is_small_), n_hits_(&sr.n_hits_), - is_linear_(sr.source_gradients_.size() > 0), lock_(&sr.lock_), - volume_(&sr.volume_), volume_t_(&sr.volume_t_), volume_sq_(&sr.volume_sq_), - volume_sq_t_(&sr.volume_sq_t_), volume_naive_(&sr.volume_naive_), + density_mult_(&sr.density_mult_), is_small_(&sr.is_small_), + n_hits_(&sr.n_hits_), is_linear_(sr.source_gradients_.size() > 0), + lock_(&sr.lock_), volume_(&sr.volume_), volume_t_(&sr.volume_t_), + volume_sq_(&sr.volume_sq_), volume_sq_t_(&sr.volume_sq_t_), + volume_naive_(&sr.volume_naive_), position_recorded_(&sr.position_recorded_), external_source_present_(&sr.external_source_present_), position_(&sr.position_), centroid_(&sr.centroid_), @@ -70,6 +71,7 @@ void SourceRegionContainer::push_back(const SourceRegion& sr) // Scalar fields material_.push_back(sr.material_); + density_mult_.push_back(sr.density_mult_); is_small_.push_back(sr.is_small_); n_hits_.push_back(sr.n_hits_); lock_.push_back(sr.lock_); @@ -123,6 +125,7 @@ void SourceRegionContainer::assign( // Clear existing data n_source_regions_ = 0; material_.clear(); + density_mult_.clear(); is_small_.clear(); n_hits_.clear(); lock_.clear(); @@ -180,6 +183,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) SourceRegionHandle handle; handle.negroups_ = negroups(); handle.material_ = &material(sr); + handle.density_mult_ = &density_mult(sr); handle.is_small_ = &is_small(sr); handle.n_hits_ = &n_hits(sr); handle.is_linear_ = is_linear(); diff --git a/src/reaction.cpp b/src/reaction.cpp index d96790c6d4..9ac5a1f528 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -70,9 +70,8 @@ Reaction::Reaction( if (settings::use_decay_photons) { // Remove photon products for D1S method - products_.erase( - std::remove_if(products_.begin(), products_.end(), - [](const auto& p) { return p.particle_ == ParticleType::photon; }), + products_.erase(std::remove_if(products_.begin(), products_.end(), + [](const auto& p) { return p.particle_.is_photon(); }), products_.end()); // Determine product for D1S method diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 3ba2c0cfb0..ee560d6077 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -26,7 +26,7 @@ ReactionProduct::ReactionProduct(hid_t group) // Read particle type std::string temp; read_attribute(group, "particle", temp); - particle_ = str_to_particle_type(temp); + particle_ = ParticleType {temp}; // Read emission mode and decay rate read_attribute(group, "emission_mode", temp); @@ -42,7 +42,7 @@ ReactionProduct::ReactionProduct(hid_t group) if (emission_mode_ == EmissionMode::delayed) { if (attribute_exists(group, "decay_rate")) { read_attribute(group, "decay_rate", decay_rate_); - } else if (particle_ == ParticleType::neutron) { + } else if (particle_.is_neutron()) { warning(fmt::format("Decay rate doesn't exist for delayed neutron " "emission ({}).", object_name(group))); @@ -85,7 +85,7 @@ ReactionProduct::ReactionProduct(hid_t group) ReactionProduct::ReactionProduct(const ChainNuclide::Product& product) { - particle_ = ParticleType::photon; + particle_ = ParticleType::photon(); emission_mode_ = EmissionMode::delayed; // Get chain nuclide object for radionuclide diff --git a/src/simulation.cpp b/src/simulation.cpp index e0ff12864f..5eb3a21dce 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -701,7 +701,7 @@ void initialize_data() for (const auto& nuc : data::nuclides) { if (nuc->grid_.size() >= 1) { - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); data::energy_min[neutron] = std::max(data::energy_min[neutron], nuc->grid_[0].energy.front()); data::energy_max[neutron] = @@ -712,7 +712,7 @@ void initialize_data() if (settings::photon_transport) { for (const auto& elem : data::elements) { if (elem->energy_.size() >= 1) { - int photon = static_cast(ParticleType::photon); + int photon = ParticleType::photon().transport_index(); int n = elem->energy_.size(); data::energy_min[photon] = std::max(data::energy_min[photon], std::exp(elem->energy_(1))); @@ -725,9 +725,9 @@ void initialize_data() // Determine if minimum/maximum energy for bremsstrahlung is greater/less // than the current minimum/maximum if (data::ttb_e_grid.size() >= 1) { - int photon = static_cast(ParticleType::photon); - int electron = static_cast(ParticleType::electron); - int positron = static_cast(ParticleType::positron); + int photon = ParticleType::photon().transport_index(); + int electron = ParticleType::electron().transport_index(); + int positron = ParticleType::positron().transport_index(); int n_e = data::ttb_e_grid.size(); const std::vector charged = {electron, positron}; @@ -751,7 +751,7 @@ void initialize_data() // grid has not been allocated if (nuc->grid_.size() > 0) { double max_E = nuc->grid_[0].energy.back(); - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); if (max_E == data::energy_max[neutron]) { write_message(7, "Maximum neutron transport energy: {} eV for {}", data::energy_max[neutron], nuc->name_); @@ -768,7 +768,7 @@ void initialize_data() for (auto& nuc : data::nuclides) { nuc->init_grid(); } - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); simulation::log_spacing = std::log(data::energy_max[neutron] / data::energy_min[neutron]) / settings::n_log_bins; diff --git a/src/source.cpp b/src/source.cpp index b30b964d39..8bd9c78935 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -37,6 +37,20 @@ namespace openmc { +namespace { + +void validate_particle_type(ParticleType type, const std::string& context) +{ + if (type.is_transportable()) + return; + + fatal_error( + fmt::format("Unsupported source particle type '{}' (PDG {}) in {}.", + type.str(), type.pdg_number(), context)); +} + +} // namespace + //============================================================================== // Global variables //============================================================================== @@ -284,22 +298,15 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) { // Check for particle type if (check_for_node(node, "particle")) { - auto temp_str = get_node_value(node, "particle", true, true); - if (temp_str == "neutron") { - particle_ = ParticleType::neutron; - } else if (temp_str == "photon") { - particle_ = ParticleType::photon; + auto temp_str = get_node_value(node, "particle", false, true); + particle_ = ParticleType(temp_str); + if (particle_ == ParticleType::photon() || + particle_ == ParticleType::electron() || + particle_ == ParticleType::positron()) { settings::photon_transport = true; - } else if (temp_str == "electron") { - particle_ = ParticleType::electron; - settings::photon_transport = true; - } else if (temp_str == "positron") { - particle_ = ParticleType::positron; - settings::photon_transport = true; - } else { - fatal_error(std::string("Unknown source particle type: ") + temp_str); } } + validate_particle_type(particle_, "IndependentSource"); // Check for external source file if (check_for_node(node, "file")) { @@ -390,7 +397,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const // Sample energy and time for neutron and photon sources if (settings::solver_type != SolverType::RANDOM_RAY) { // Check for monoenergetic source above maximum particle energy - auto p = static_cast(particle_); + auto p = particle_.transport_index(); auto energy_ptr = dynamic_cast(energy_.get()); if (energy_ptr) { auto energies = xt::adapt(energy_ptr->x()); @@ -472,6 +479,11 @@ void FileSource::load_sites_from_file(const std::string& path) // Close file file_close(file_id); } + + // Make sure particles in source file have valid types + for (const auto& site : this->sites_) { + validate_particle_type(site.particle, "FileSource"); + } } SourceSite FileSource::sample(uint64_t* seed) const @@ -585,6 +597,11 @@ MeshSource::MeshSource(pugi::xml_node node) : Source(node) std::make_unique(mesh_idx, elem_index)); } + // Make sure sources use valid particle types + for (const auto& src : sources_) { + validate_particle_type(src->particle_type(), "MeshSource"); + } + // the number of source distributions should either be one or equal to the // number of mesh elements if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) { diff --git a/src/state_point.cpp b/src/state_point.cpp index 8ccebeb05a..455299529b 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -22,6 +22,7 @@ #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" #include "openmc/output.h" +#include "openmc/particle_type.h" #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/tallies/derivative.h" @@ -632,6 +633,7 @@ void write_h5_source_point(const char* filename, span source_bank, if (mpi::master || parallel) { file_id = file_open(filename_.c_str(), 'w', true); write_attribute(file_id, "filetype", "source"); + write_attribute(file_id, "version", VERSION_STATEPOINT); } // Get pointer to source bank and write to file @@ -677,6 +679,16 @@ std::string dtype_member_names(hid_t dtype_id) void read_source_bank( hid_t group_id, vector& sites, bool distribute) { + bool legacy_particle_codes = true; + if (attribute_exists(group_id, "version")) { + array version; + read_attribute(group_id, "version", version); + if (version[0] > VERSION_STATEPOINT[0] || + (version[0] == VERSION_STATEPOINT[0] && version[1] >= 2)) { + legacy_particle_codes = false; + } + } + hid_t banktype = h5banktype(true); // Open the dataset @@ -738,6 +750,12 @@ void read_source_bank( H5Sclose(memspace); H5Dclose(dset); H5Tclose(banktype); + + if (legacy_particle_codes) { + for (auto& site : sites) { + site.particle = legacy_particle_index_to_type(site.particle.pdg_number()); + } + } } void write_unstructured_mesh_results() diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index eef1d1e63c..031068f3a8 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -13,7 +13,7 @@ void ParticleFilter::from_xml(pugi::xml_node node) // Convert to vector of ParticleType vector types; for (auto& p : particles) { - types.push_back(str_to_particle_type(p)); + types.emplace_back(p); } this->set_particles(types); } @@ -47,7 +47,7 @@ void ParticleFilter::to_statepoint(hid_t filter_group) const Filter::to_statepoint(filter_group); vector particles; for (auto p : particles_) { - particles.push_back(particle_type_to_str(p)); + particles.push_back(p.str()); } write_dataset(filter_group, "bins", particles); } @@ -55,10 +55,10 @@ void ParticleFilter::to_statepoint(hid_t filter_group) const std::string ParticleFilter::text_label(int bin) const { const auto& p = particles_.at(bin); - return fmt::format("Particle: {}", particle_type_to_str(p)); + return fmt::format("Particle: {}", p.str()); } -extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[]) +extern "C" int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]) { if (int err = verify_filter(idx)) return err; @@ -68,7 +68,7 @@ extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[]) if (pf) { const auto& particles = pf->particles(); for (int i = 0; i < particles.size(); i++) { - bins[i] = static_cast(particles[i]); + bins[i] = particles[i].pdg_number(); } } else { set_errmsg("The filter at the specified index is not a ParticleFilter"); diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index a134cd7929..9d37a6ac74 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -231,12 +231,12 @@ Tally::Tally(pugi::xml_node node) const auto& f = model::tally_filters[particle_filter_index].get(); auto pf = dynamic_cast(f); for (auto p : pf->particles()) { - if (p != ParticleType::neutron) { + if (!p.is_neutron()) { warning(fmt::format( "Particle filter other than NEUTRON used with " "photon transport turned off. All tallies for particle type {}" " will have no scores", - static_cast(p))); + p.str())); } } } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 0cefb43e0e..7b1e5a9c4f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -328,10 +328,10 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux, //! Helper function to obtain reaction Q value for photons and charged particles double get_reaction_q_value(const Particle& p) { - if (p.type() == ParticleType::photon && p.event_mt() == PAIR_PROD) { + if (p.type().is_photon() && p.event_mt() == PAIR_PROD) { // pair production return -2 * MASS_ELECTRON_EV; - } else if (p.type() == ParticleType::positron) { + } else if (p.type() == ParticleType::positron()) { // positron annihilation return 2 * MASS_ELECTRON_EV; } else { @@ -344,7 +344,7 @@ double get_reaction_q_value(const Particle& p) double score_particle_heating(const Particle& p, const Tally& tally, double flux, int rxn_bin, int i_nuclide, double atom_density) { - if (p.type() == ParticleType::neutron) + if (p.type().is_neutron()) return score_neutron_heating( p, tally, flux, rxn_bin, i_nuclide, atom_density); if (i_nuclide == -1 || i_nuclide == p.event_nuclide() || @@ -584,8 +584,6 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, // Get the pre-collision energy of the particle. auto E = p.E_last(); - using Type = ParticleType; - for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; auto score_index = start_index + i; @@ -598,9 +596,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, case SCORE_TOTAL: if (i_nuclide >= 0) { - if (p.type() == Type::neutron) { + if (p.type().is_neutron()) { score = p.neutron_xs(i_nuclide).total * atom_density * flux; - } else if (p.type() == Type::photon) { + } else if (p.type().is_photon()) { score = p.photon_xs(i_nuclide).total * atom_density * flux; } } else { @@ -609,7 +607,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, break; case SCORE_INVERSE_VELOCITY: - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; // Score inverse velocity in units of s/cm. @@ -617,11 +615,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, break; case SCORE_SCATTER: - if (p.type() != Type::neutron && p.type() != Type::photon) + if (!p.type().is_neutron() && !p.type().is_photon()) continue; if (i_nuclide >= 0) { - if (p.type() == Type::neutron) { + if (p.type().is_neutron()) { const auto& micro = p.neutron_xs(i_nuclide); score = (micro.total - micro.absorption) * atom_density * flux; } else { @@ -629,7 +627,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, score = (micro.coherent + micro.incoherent) * atom_density * flux; } } else { - if (p.type() == Type::neutron) { + if (p.type().is_neutron()) { score = (p.macro_xs().total - p.macro_xs().absorption) * flux; } else { score = (p.macro_xs().coherent + p.macro_xs().incoherent) * flux; @@ -638,11 +636,11 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, break; case SCORE_ABSORPTION: - if (p.type() != Type::neutron && p.type() != Type::photon) + if (!p.type().is_neutron() && !p.type().is_photon()) continue; if (i_nuclide >= 0) { - if (p.type() == Type::neutron) { + if (p.type().is_neutron()) { score = p.neutron_xs(i_nuclide).absorption * atom_density * flux; } else { const auto& xs = p.photon_xs(i_nuclide); @@ -650,7 +648,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, (xs.total - xs.coherent - xs.incoherent) * atom_density * flux; } } else { - if (p.type() == Type::neutron) { + if (p.type().is_neutron()) { score = p.macro_xs().absorption * flux; } else { score = @@ -806,7 +804,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, // this loop. for (auto d = 1; d < rxn.products_.size(); ++d) { const auto& product = rxn.products_[d]; - if (product.particle_ != Type::neutron) + if (!product.particle_.is_neutron()) continue; auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); @@ -860,7 +858,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, // this loop. for (auto d = 1; d < rxn.products_.size(); ++d) { const auto& product = rxn.products_[d]; - if (product.particle_ != Type::neutron) + if (!product.particle_.is_neutron()) continue; auto yield = @@ -911,7 +909,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, continue; case ELASTIC: - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; if (i_nuclide >= 0) { @@ -942,7 +940,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, break; case SCORE_IFP_TIME_NUM: - if ((p.type() == Type::neutron) && (p.fission())) { + if (p.type().is_neutron() && p.fission()) { const auto& lifetime = simulation::ifp_source_lifetime_bank[p.current_work() - 1][0]; score = lifetime * p.wgt_last(); @@ -950,7 +948,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, break; case SCORE_IFP_BETA_NUM: - if ((p.type() == Type::neutron) && (p.fission())) { + if (p.type().is_neutron() && p.fission()) { const auto& delayed_group = simulation::ifp_source_delayed_group_bank[p.current_work() - 1][0]; if (delayed_group > 0) { @@ -968,7 +966,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, break; case SCORE_IFP_DENOM: - if ((p.type() == Type::neutron) && (p.fission())) + if (p.type().is_neutron() && p.fission()) score = p.wgt_last(); break; @@ -984,7 +982,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (!simulation::need_depletion_rx) goto default_case; - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; int m; @@ -1017,7 +1015,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, case INCOHERENT: case PHOTOELECTRIC: case PAIR_PROD: - if (p.type() != Type::photon) + if (!p.type().is_photon()) continue; if (i_nuclide >= 0) { @@ -1047,7 +1045,7 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, // The default block is really only meant for redundant neutron reactions // (e.g. 444, 901) - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; // Any other cross section has to be calculated on-the-fly @@ -1101,8 +1099,6 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, p.neutron_xs(p.event_nuclide()).total : 0.0; - using Type = ParticleType; - for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; auto score_index = start_index + i; @@ -1113,7 +1109,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, // All events score to a flux bin. We actually use a collision estimator // in place of an analog one since there is no way to count 'events' // exactly for the flux - if (p.type() == Type::neutron || p.type() == Type::photon) { + if (p.type().is_neutron() || p.type().is_photon()) { score = flux * p.wgt_last() / p.macro_xs().total; } else { score = 0.; @@ -1127,7 +1123,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, break; case SCORE_INVERSE_VELOCITY: - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; // All events score to an inverse velocity bin. We actually use a @@ -1137,7 +1133,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, break; case SCORE_SCATTER: - if (p.type() != Type::neutron && p.type() != Type::photon) + if (!p.type().is_neutron() && !p.type().is_photon()) continue; // Skip any event where the particle didn't scatter @@ -1149,7 +1145,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, break; case SCORE_NU_SCATTER: - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; // Only analog estimators are available. @@ -1174,7 +1170,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, break; case SCORE_ABSORPTION: - if (p.type() != Type::neutron && p.type() != Type::photon) + if (!p.type().is_neutron() && !p.type().is_photon()) continue; if (settings::survival_biasing) { @@ -1403,7 +1399,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, // this loop. for (auto d = 1; d < rxn.products_.size(); ++d) { const auto& product = rxn.products_[d]; - if (product.particle_ != Type::neutron) + if (!product.particle_.is_neutron()) continue; auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); @@ -1495,7 +1491,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, continue; case ELASTIC: - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; // Check if event MT matches @@ -1524,7 +1520,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, if (!simulation::need_depletion_rx) goto default_case; - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; // Check if the event MT matches @@ -1537,7 +1533,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, case INCOHERENT: case PHOTOELECTRIC: case PAIR_PROD: - if (p.type() != Type::photon) + if (!p.type().is_photon()) continue; if (score_bin == PHOTOELECTRIC) { @@ -1564,7 +1560,7 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, // The default block is really only meant for redundant neutron reactions // (e.g. 444, 901) - if (p.type() != Type::neutron) + if (!p.type().is_neutron()) continue; // Any other score is assumed to be a MT number. Thus, we just need @@ -2288,10 +2284,7 @@ void score_analog_tally_ce(Particle& p) // Since electrons/positrons are not transported, we assign a flux of zero. // Note that the heating score does NOT use the flux and will be non-zero for // electrons/positrons. - double flux = - (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) - ? 1.0 - : 0.0; + double flux = (p.type().is_neutron() || p.type().is_photon()) ? 1.0 : 0.0; for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2419,7 +2412,7 @@ void score_tracklength_tally_general( if (j == C_NONE) { // Determine log union grid index if (i_log_union == C_NONE) { - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); i_log_union = std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing; } @@ -2516,7 +2509,7 @@ void score_collision_tally(Particle& p) { // Determine the collision estimate of the flux double flux = 0.0; - if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) { + if (p.type().is_neutron() || p.type().is_photon()) { flux = p.wgt_last() / p.macro_xs().total; } @@ -2550,7 +2543,7 @@ void score_collision_tally(Particle& p) if (j == C_NONE) { // Determine log union grid index if (i_log_union == C_NONE) { - int neutron = static_cast(ParticleType::neutron); + int neutron = ParticleType::neutron().transport_index(); i_log_union = std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing; } diff --git a/src/track_output.cpp b/src/track_output.cpp index f4344d50f7..e86f774f04 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -122,7 +122,7 @@ void finalize_particle_track(Particle& p) int offset = 0; for (auto& track_i : p.tracks()) { offsets.push_back(offset); - particles.push_back(static_cast(track_i.particle)); + particles.push_back(track_i.particle.pdg_number()); offset += track_i.states.size(); tracks.insert(tracks.end(), track_i.states.begin(), track_i.states.end()); } diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 4838e45915..5ca4addbbf 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -59,7 +59,7 @@ void apply_weight_windows(Particle& p) return; // WW on photon and neutron only - if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon) + if (!p.type().is_neutron() && !p.type().is_photon()) return; // skip dead or no energy @@ -179,7 +179,7 @@ WeightWindows::WeightWindows(pugi::xml_node node) // get the particle type auto particle_type_str = std::string(get_node_value(node, "particle_type")); - particle_type_ = openmc::str_to_particle_type(particle_type_str); + particle_type_ = ParticleType {particle_type_str}; // Determine associated mesh int32_t mesh_id = std::stoi(get_node_value(node, "mesh")); @@ -252,7 +252,7 @@ WeightWindows* WeightWindows::from_hdf5( std::string particle_type; read_dataset(ww_group, "particle_type", particle_type); - wws->particle_type_ = openmc::str_to_particle_type(particle_type); + wws->particle_type_ = ParticleType {particle_type}; read_dataset(ww_group, "energy_bounds", wws->energy_bounds_); @@ -284,7 +284,10 @@ void WeightWindows::set_defaults() { // set energy bounds to the min/max energy supported by the data if (energy_bounds_.size() == 0) { - int p_type = static_cast(particle_type_); + int p_type = particle_type_.transport_index(); + if (p_type == C_NONE) { + fatal_error("Weight windows particle is not supported for transport."); + } energy_bounds_.push_back(data::energy_min[p_type]); energy_bounds_.push_back(data::energy_max[p_type]); } @@ -345,10 +348,9 @@ void WeightWindows::set_energy_bounds(span bounds) void WeightWindows::set_particle_type(ParticleType p_type) { - if (p_type != ParticleType::neutron && p_type != ParticleType::photon) - fatal_error( - fmt::format("Particle type '{}' cannot be applied to weight windows.", - particle_type_to_str(p_type))); + if (!p_type.is_neutron() && !p_type.is_photon()) + fatal_error(fmt::format( + "Particle type '{}' cannot be applied to weight windows.", p_type.str())); particle_type_ = p_type; } @@ -608,8 +610,7 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value, if (p_it == particles.end()) { auto msg = fmt::format("Particle type '{}' not present on Filter {} for " "Tally {} used to update WeightWindows {}", - particle_type_to_str(this->particle_type_), pf->id(), tally->id(), - this->id()); + this->particle_type_.str(), pf->id(), tally->id(), this->id()); fatal_error(msg); } @@ -818,8 +819,7 @@ void WeightWindows::to_hdf5(hid_t group) const hid_t ww_group = create_group(group, fmt::format("weight_windows_{}", id())); write_dataset(ww_group, "mesh", this->mesh()->id()); - write_dataset( - ww_group, "particle_type", openmc::particle_type_to_str(particle_type_)); + write_dataset(ww_group, "particle_type", particle_type_.str()); write_dataset(ww_group, "energy_bounds", energy_bounds_); write_dataset(ww_group, "lower_ww_bounds", lower_ww_); write_dataset(ww_group, "upper_ww_bounds", upper_ww_); @@ -846,8 +846,8 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) max_realizations_, active_batches); warning(msg); } - auto tmp_str = get_node_value(node, "particle_type", true, true); - auto particle_type = str_to_particle_type(tmp_str); + auto tmp_str = get_node_value(node, "particle_type", false, true); + auto particle_type = ParticleType {tmp_str}; update_interval_ = std::stoi(get_node_value(node, "update_interval")); on_the_fly_ = get_node_value_bool(node, "on_the_fly"); @@ -856,7 +856,10 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) if (check_for_node(node, "energy_bounds")) { e_bounds = get_node_array(node, "energy_bounds"); } else { - int p_type = static_cast(particle_type); + int p_type = particle_type.transport_index(); + if (p_type == C_NONE) { + fatal_error("Weight windows particle is not supported for transport."); + } e_bounds.push_back(data::energy_min[p_type]); e_bounds.push_back(data::energy_max[p_type]); } @@ -1110,23 +1113,25 @@ extern "C" int openmc_weight_windows_get_energy_bounds( return 0; } -extern "C" int openmc_weight_windows_set_particle(int32_t index, int particle) +extern "C" int openmc_weight_windows_set_particle( + int32_t index, int32_t particle) { if (int err = verify_ww_index(index)) return err; const auto& wws = variance_reduction::weight_windows.at(index); - wws->set_particle_type(static_cast(particle)); + wws->set_particle_type(ParticleType {particle}); return 0; } -extern "C" int openmc_weight_windows_get_particle(int32_t index, int* particle) +extern "C" int openmc_weight_windows_get_particle( + int32_t index, int32_t* particle) { if (int err = verify_ww_index(index)) return err; const auto& wws = variance_reduction::weight_windows.at(index); - *particle = static_cast(wws->particle_type()); + *particle = wws->particle_type().pdg_number(); return 0; } diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 5f87db9eac..20341b420f 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -6,6 +6,7 @@ set(TEST_NAMES test_math test_mcpl_stat_sum test_mesh + test_region # Add additional unit test files here ) diff --git a/tests/cpp_unit_tests/test_distribution.cpp b/tests/cpp_unit_tests/test_distribution.cpp index 2024b812e6..f7c480c8d1 100644 --- a/tests/cpp_unit_tests/test_distribution.cpp +++ b/tests/cpp_unit_tests/test_distribution.cpp @@ -1,4 +1,6 @@ #include "openmc/distribution.h" +#include "openmc/distribution_spatial.h" +#include "openmc/position.h" #include "openmc/random_lcg.h" #include #include @@ -79,3 +81,14 @@ TEST_CASE("Test alias sampling method for pugixml constructor") REQUIRE(dist.alias()[i] == correct_alias[i]); } } + +TEST_CASE("Test construction of SpatialBox with parameters") +{ + openmc::Position ll {-1, -2, -3}; + openmc::Position ur {30, 15, 5}; + openmc::SpatialBox box(ll, ur); + + REQUIRE(box.lower_left() == openmc::Position {-1, -2, -3}); + REQUIRE(box.upper_right() == openmc::Position {30, 15, 5}); + REQUIRE_FALSE(box.only_fissionable()); +} diff --git a/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp b/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp index 909830e035..d0c26f26c2 100644 --- a/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp +++ b/tests/cpp_unit_tests/test_mcpl_stat_sum.cpp @@ -25,7 +25,7 @@ TEST_CASE("MCPL stat:sum field") // Initialize test particles for (int i = 0; i < 100; ++i) { - source_bank[i].particle = openmc::ParticleType::neutron; + source_bank[i].particle = openmc::ParticleType::neutron(); source_bank[i].r = {i * 0.1, i * 0.2, i * 0.3}; source_bank[i].u = {0.0, 0.0, 1.0}; source_bank[i].E = 2.0e6; // 2 MeV @@ -63,7 +63,7 @@ TEST_CASE("MCPL stat:sum field") // Initialize particles for (int i = 0; i < count; ++i) { - source_bank[i].particle = openmc::ParticleType::neutron; + source_bank[i].particle = openmc::ParticleType::neutron(); source_bank[i].r = {0.0, 0.0, 0.0}; source_bank[i].u = {0.0, 0.0, 1.0}; source_bank[i].E = 1.0e6; diff --git a/tests/cpp_unit_tests/test_region.cpp b/tests/cpp_unit_tests/test_region.cpp new file mode 100644 index 0000000000..b3d9a14271 --- /dev/null +++ b/tests/cpp_unit_tests/test_region.cpp @@ -0,0 +1,101 @@ +#include + +#include "openmc/cell.h" +#include "openmc/surface.h" + +#include + +namespace { + +// Helper class to set up and tear down test surfaces +class SurfaceFixture { +public: + SurfaceFixture() + { + pugi::xml_document doc; + pugi::xml_node surf_node = doc.append_child("surface"); + surf_node.set_name("surface"); + surf_node.append_attribute("id") = "0"; + surf_node.append_attribute("type") = "x-plane"; + surf_node.append_attribute("coeffs") = "1"; + + for (int i = 1; i < 10; ++i) { + surf_node.attribute("id") = i; + openmc::model::surfaces.push_back( + std::make_unique(surf_node)); + openmc::model::surface_map[i] = i - 1; + } + } + + ~SurfaceFixture() + { + openmc::model::surfaces.clear(); + openmc::model::surface_map.clear(); + } +}; + +} // anonymous namespace + +TEST_CASE("Test region simplification") +{ + SurfaceFixture fixture; + + SECTION("Original bug case from issue #3685") + { + // Input: "-1 2 (-3 4) | (-5 6)" was being incorrectly interpreted + auto region = openmc::Region("(-1 2 (-3 4) | (-5 6))", 0); + REQUIRE(region.str() == " ( ( -1 2 ( -3 4 ) ) | ( -5 6 ) )"); + } + + SECTION("Simple union - no extra parentheses needed") + { + auto region = openmc::Region("1 | 2", 0); + REQUIRE(region.str() == " 1 | 2"); + } + + SECTION("Intersection then union") + { + // Intersection should have higher precedence, so (1 2) grouped + auto region = openmc::Region("1 2 | 3", 0); + REQUIRE(region.str() == " ( 1 2 ) | 3"); + } + + SECTION("Union then intersection") + { + // The (2 3) intersection should be grouped + auto region = openmc::Region("1 | 2 3", 0); + REQUIRE(region.str() == " 1 | ( 2 3 )"); + } + + SECTION("Nested parentheses preserved") + { + // These parentheses are meaningful and should be preserved + auto region = openmc::Region("(1 | 2) (3 | 4)", 0); + REQUIRE(region.str() == " ( 1 | 2 ) ( 3 | 4 )"); + } + + SECTION("Deep nesting") + { + auto region = openmc::Region("((1 2) | (3 4)) 5", 0); + REQUIRE(region.str() == " ( ( 1 2 ) | ( 3 4 ) ) 5"); + } + + SECTION("Multiple unions") + { + auto region = openmc::Region("1 | 2 | 3", 0); + REQUIRE(region.str() == " 1 | 2 | 3"); + } + + SECTION("Multiple intersections") + { + auto region = openmc::Region("1 2 3", 0); + // Simple cell - no operators in output + REQUIRE(region.str() == " 1 2 3"); + } + + SECTION("Complex mixed expression") + { + auto region = openmc::Region("1 2 | 3 4 | 5 6", 0); + REQUIRE(region.str() == " ( 1 2 ) | ( 3 4 ) | ( 5 6 )"); + } +} diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index 1bb00129d0..9595765d73 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -243,4 +243,4 @@ class DummyOperator(TransportOperator): Maps cell name to index in global geometry. """ - return self.volume, self.nuc_list, self.local_mats, self.burnable_mats + return self.volume, self.nuc_list, self.local_mats, self.burnable_mats, {"1": ""} diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/__init__.py b/tests/regression_tests/random_ray_auto_convert_kappa_fission/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat new file mode 100644 index 0000000000..9f3a827f68 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/inputs_true.dat @@ -0,0 +1,73 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + + + -0.63 -0.63 -1 0.63 0.63 1 + + + true + + + multi-group + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 + 150.0 + + + + + + linear + + + 2 2 + -0.63 -0.63 + 0.63 0.63 + + + + + 1 + + + 61 + kappa-fission + + + diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/results_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/results_true.dat new file mode 100644 index 0000000000..27b32d787b --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/infinite_medium/results_true.dat @@ -0,0 +1,5 @@ +k-combined: +7.479770E-01 1.624548E-02 +tally 1: +2.965503E+08 +1.762157E+16 diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat new file mode 100644 index 0000000000..edf68f7e25 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/inputs_true.dat @@ -0,0 +1,73 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + + + -0.63 -0.63 -1 0.63 0.63 1 + + + true + + + multi-group + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 + 150.0 + + + + + + linear + + + 2 2 + -0.63 -0.63 + 0.63 0.63 + + + + + 1 + + + 97 + kappa-fission + + + diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/results_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/results_true.dat new file mode 100644 index 0000000000..5a379db7e3 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/material_wise/results_true.dat @@ -0,0 +1,5 @@ +k-combined: +7.372542E-01 6.967831E-03 +tally 1: +2.909255E+08 +1.693395E+16 diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat new file mode 100644 index 0000000000..edf68f7e25 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/inputs_true.dat @@ -0,0 +1,73 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + + + -0.63 -0.63 -1 0.63 0.63 1 + + + true + + + multi-group + + + + -0.63 -0.63 -1.0 0.63 0.63 1.0 + + + 30.0 + 150.0 + + + + + + linear + + + 2 2 + -0.63 -0.63 + 0.63 0.63 + + + + + 1 + + + 97 + kappa-fission + + + diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/results_true.dat b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/results_true.dat new file mode 100644 index 0000000000..a60e5aa3c6 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/stochastic_slab/results_true.dat @@ -0,0 +1,5 @@ +k-combined: +6.413334E-01 2.083132E-02 +tally 1: +2.541463E+08 +1.297259E+16 diff --git a/tests/regression_tests/random_ray_auto_convert_kappa_fission/test.py b/tests/regression_tests/random_ray_auto_convert_kappa_fission/test.py new file mode 100644 index 0000000000..6decf165a7 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_kappa_fission/test.py @@ -0,0 +1,60 @@ +import os + +import openmc +from openmc.examples import pwr_pin_cell +from openmc import RegularMesh +from openmc.utility_funcs import change_directory +import pytest + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +@pytest.mark.parametrize("method", ["material_wise", "stochastic_slab", "infinite_medium"]) +def test_random_ray_auto_convert(method): + with change_directory(method): + openmc.reset_auto_ids() + + # Start with a normal continuous energy model + model = pwr_pin_cell() + + # Convert to a multi-group model + model.convert_to_multigroup( + method=method, groups='CASMO-2', nparticles=100, + overwrite_mgxs_library=False, mgxs_path="mgxs.h5" + ) + + # Convert to a random ray model + model.convert_to_random_ray() + + # Set the number of particles + model.settings.particles = 100 + + # Overlay a basic 2x2 mesh + n = 2 + mesh = RegularMesh() + mesh.dimension = (n, n) + bbox = model.geometry.bounding_box + mesh.lower_left = (bbox.lower_left[0], bbox.lower_left[1]) + mesh.upper_right = (bbox.upper_right[0], bbox.upper_right[1]) + model.settings.random_ray['source_region_meshes'] = [ + (mesh, [model.geometry.root_universe])] + + # Set the source shape to linear + model.settings.random_ray['source_shape'] = 'linear' + + # Set a material tally + t = openmc.Tally(name = 'KF Tally') + t.filters = [openmc.MaterialFilter(bins=1)] + t.scores = ['kappa-fission'] + model.tallies.append(t) + + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/random_ray_cell_density/__init__.py b/tests/regression_tests/random_ray_cell_density/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat b/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat new file mode 100644 index 0000000000..0dd354cf02 --- /dev/null +++ b/tests/regression_tests/random_ray_cell_density/eigen/inputs_true.dat @@ -0,0 +1,109 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.126 0.126 + 10 10 + -0.63 -0.63 + +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 + + + 1.26 1.26 + 2 2 + -1.26 -1.26 + +2 2 +2 5 + + + + + + + + + + + + + + + + + + + + + eigenvalue + 100 + 10 + 5 + multi-group + + 100.0 + 20.0 + + + -1.26 -1.26 -1 1.26 1.26 1 + + + true + + + + + 2 2 + -1.26 -1.26 + 1.26 1.26 + + + 1 + + + 1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0 + + + 1 2 + flux fission nu-fission + analog + + + diff --git a/tests/regression_tests/random_ray_cell_density/eigen/results_true.dat b/tests/regression_tests/random_ray_cell_density/eigen/results_true.dat new file mode 100644 index 0000000000..3e6ce40392 --- /dev/null +++ b/tests/regression_tests/random_ray_cell_density/eigen/results_true.dat @@ -0,0 +1,171 @@ +k-combined: +7.606488E-01 9.025978E-03 +tally 1: +6.209503E-01 +7.732732E-02 +4.335729E-01 +3.769026E-02 +1.055230E+00 +2.232538E-01 +4.077492E-01 +3.335434E-02 +1.188353E-01 +2.833644E-03 +2.892214E-01 +1.678476E-02 +2.715245E-01 +1.488194E-02 +1.736084E-02 +6.079207E-05 +4.225281E-02 +3.600946E-04 +3.700491E-01 +2.796457E-02 +2.369715E-02 +1.145314E-04 +5.767413E-02 +6.784132E-04 +1.223083E+00 +3.058118E-01 +2.790791E-02 +1.593458E-04 +6.792309E-02 +9.438893E-04 +4.095779E+00 +3.377358E+00 +1.251201E-02 +3.153331E-05 +3.096010E-02 +1.930723E-04 +2.608638E+00 +1.361111E+00 +7.059742E-02 +9.968342E-04 +1.963632E-01 +7.711968E-03 +1.164050E+00 +2.710289E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +5.325486E-01 +5.675375E-02 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.058117E-01 +1.905099E-02 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.645650E-01 +4.428410E-02 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.350036E+00 +3.718611E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.538916E+00 +2.521670E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.201812E+00 +9.698678E-01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +6.776093E-01 +9.211106E-02 +2.527022E-01 +1.280619E-02 +6.150266E-01 +7.585593E-02 +4.194143E-01 +3.528585E-02 +6.303447E-02 +7.972233E-04 +1.534133E-01 +4.722259E-03 +2.743408E-01 +1.520611E-02 +8.961769E-03 +1.621446E-05 +2.181115E-02 +9.604444E-05 +3.871710E-01 +3.063162E-02 +1.288801E-02 +3.392550E-05 +3.136683E-02 +2.009537E-04 +1.259769E+00 +3.241545E-01 +1.478360E-02 +4.466102E-05 +3.598077E-02 +2.645508E-04 +4.032169E+00 +3.273866E+00 +6.215187E-03 +7.776862E-06 +1.537904E-02 +4.761620E-05 +2.583647E+00 +1.335275E+00 +3.533526E-02 +2.498005E-04 +9.828323E-02 +1.932572E-03 +7.851145E-01 +1.234783E-01 +2.966638E-01 +1.763751E-02 +7.220204E-01 +1.044737E-01 +4.505673E-01 +4.067697E-02 +6.823342E-02 +9.332472E-04 +1.660665E-01 +5.527980E-03 +2.832534E-01 +1.626181E-02 +9.297014E-03 +1.749663E-05 +2.262707E-02 +1.036392E-04 +4.056699E-01 +3.370803E-02 +1.358439E-02 +3.774180E-05 +3.306168E-02 +2.235591E-04 +1.279744E+00 +3.345476E-01 +1.512005E-02 +4.668112E-05 +3.679962E-02 +2.765169E-04 +3.917888E+00 +3.090981E+00 +6.082371E-03 +7.449625E-06 +1.505040E-02 +4.561259E-05 +2.520906E+00 +1.271106E+00 +3.477796E-02 +2.419730E-04 +9.673314E-02 +1.872015E-03 diff --git a/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat b/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat new file mode 100644 index 0000000000..e90f25973e --- /dev/null +++ b/tests/regression_tests/random_ray_cell_density/fs/inputs_true.dat @@ -0,0 +1,244 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + + + + + fixed source + 90 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + 500.0 + 100.0 + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + true + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_cell_density/fs/results_true.dat b/tests/regression_tests/random_ray_cell_density/fs/results_true.dat new file mode 100644 index 0000000000..e12f476e09 --- /dev/null +++ b/tests/regression_tests/random_ray_cell_density/fs/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +6.203077E-01 +7.706659E-02 +tally 2: +3.203448E-02 +2.058679E-04 +tally 3: +2.091970E-03 +8.765168E-07 diff --git a/tests/regression_tests/random_ray_cell_density/test.py b/tests/regression_tests/random_ray_cell_density/test.py new file mode 100644 index 0000000000..cb6062cdc1 --- /dev/null +++ b/tests/regression_tests/random_ray_cell_density/test.py @@ -0,0 +1,43 @@ +import os + +import openmc +from openmc.examples import random_ray_lattice, random_ray_three_region_cube +from openmc.utility_funcs import change_directory +import pytest + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +@pytest.mark.parametrize("run_mode", ["eigen", "fs"]) +def test_random_ray_basic(run_mode): + with change_directory(run_mode): + if run_mode == "eigen": + openmc.reset_auto_ids() + model = random_ray_lattice() + # Double the densities of the lower-left fuel pin -> cell instances [0, 9). + for id, cell in model.geometry.get_all_cells().items(): + if cell.fill.name == "UO2 fuel": + cell.density = [((i < 8) + 1.0) for i in range(24)] + + # Gold file was generated with manually scaled fuel cross sections. + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() + else: + openmc.reset_auto_ids() + model = random_ray_three_region_cube() + # Increase the density in the source region. + for id, cell in model.geometry.get_all_cells().items(): + if cell.fill.name == "source": + cell.density = 1e3 + + # Gold file was generated with manually scaled source cross sections. + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/source_dlopen/source_sampling.cpp b/tests/regression_tests/source_dlopen/source_sampling.cpp index fc61ef1fdd..678eea4be0 100644 --- a/tests/regression_tests/source_dlopen/source_sampling.cpp +++ b/tests/regression_tests/source_dlopen/source_sampling.cpp @@ -10,7 +10,7 @@ class CustomSource : public openmc::Source { { openmc::SourceSite particle; // wgt - particle.particle = openmc::ParticleType::neutron; + particle.particle = openmc::ParticleType::neutron(); particle.wgt = 1.0; // position diff --git a/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp b/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp index bf49af4c38..81f3770c62 100644 --- a/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp +++ b/tests/regression_tests/source_parameterized_dlopen/parameterized_source_sampling.cpp @@ -2,36 +2,37 @@ #include "openmc/source.h" class CustomSource : public openmc::Source { - public: - CustomSource(double energy) : energy_(energy) { } +public: + CustomSource(double energy) : energy_(energy) {} - // Samples from an instance of this class. - openmc::SourceSite sample(uint64_t* seed) const - { - openmc::SourceSite particle; - // wgt - particle.particle = openmc::ParticleType::neutron; - particle.wgt = 1.0; - // position - particle.r.x = 0.0; - particle.r.y = 0.0; - particle.r.z = 0.0; - // angle - particle.u = {1.0, 0.0, 0.0}; - particle.E = this->energy_; - particle.delayed_group = 0; + // Samples from an instance of this class. + openmc::SourceSite sample(uint64_t* seed) const + { + openmc::SourceSite particle; + // wgt + particle.particle = openmc::ParticleType::neutron(); + particle.wgt = 1.0; + // position + particle.r.x = 0.0; + particle.r.y = 0.0; + particle.r.z = 0.0; + // angle + particle.u = {1.0, 0.0, 0.0}; + particle.E = this->energy_; + particle.delayed_group = 0; - return particle; - } + return particle; + } - private: - double energy_; +private: + double energy_; }; -// A function to create a unique pointer to an instance of this class when generated -// via a plugin call using dlopen/dlsym. -// You must have external C linkage here otherwise dlopen will not find the file -extern "C" std::unique_ptr openmc_create_source(std::string parameter) +// A function to create a unique pointer to an instance of this class when +// generated via a plugin call using dlopen/dlsym. You must have external C +// linkage here otherwise dlopen will not find the file +extern "C" std::unique_ptr openmc_create_source( + std::string parameter) { double energy = std::stod(parameter); return std::make_unique(energy); diff --git a/tests/regression_tests/surface_source/surface_source_true.h5 b/tests/regression_tests/surface_source/surface_source_true.h5 index d343d12aed..02e94c90b9 100644 Binary files a/tests/regression_tests/surface_source/surface_source_true.h5 and b/tests/regression_tests/surface_source/surface_source_true.h5 differ diff --git a/tests/regression_tests/surface_source_write/test.py b/tests/regression_tests/surface_source_write/test.py index f144eb82a7..094df1f8b8 100644 --- a/tests/regression_tests/surface_source_write/test.py +++ b/tests/regression_tests/surface_source_write/test.py @@ -634,7 +634,7 @@ def return_surface_source_data(filepath): wgt = point.wgt delayed_group = point.delayed_group surf_id = point.surf_id - particle = point.particle + particle = point.particle.pdg_number key = ( f"{r[0]:.10e} {r[1]:.10e} {r[2]:.10e} {u[0]:.10e} {u[1]:.10e} {u[2]:.10e}" f"{e:.10e} {time:.10e} {wgt:.10e} {delayed_group} {surf_id} {particle}" diff --git a/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py b/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py index 8166ba68c8..366ce18efb 100644 --- a/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py +++ b/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py @@ -131,7 +131,7 @@ def test_write_data_to_vtk_size_mismatch(mesh): mesh : openmc.StructuredMesh The mesh to test """ - right_size = mesh.num_mesh_cells + right_size = mesh.n_elements data = np.random.random(right_size + 1) # Error message has \ in to escape characters that are otherwise recognized @@ -139,7 +139,7 @@ def test_write_data_to_vtk_size_mismatch(mesh): # string when using the match argument as that uses regular expression expected_error_msg = ( fr"The size of the dataset 'label' \({len(data)}\) should be equal to " - fr"the number of mesh cells \({mesh.num_mesh_cells}\)" + fr"the number of mesh cells \({mesh.n_elements}\)" ) with pytest.raises(ValueError, match=expected_error_msg): mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data}) diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index eace1976ef..2f17138805 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -113,6 +113,11 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts assert atoms[0] == pytest.approx(n0) assert atoms[1] / atoms[0] == pytest.approx(0.5, rel=tolerance) + # Check that material name is preserved in depletion results + step_result = results[0] + mat_from_results = step_result.get_material(f"{w.id}") + assert mat_from_results.name == 'tungsten' + def test_decay(run_in_tmpdir): """Test decay-only timesteps where no transport solve is performed""" diff --git a/tests/unit_tests/test_deplete_decay.py b/tests/unit_tests/test_deplete_decay.py index 6e7b0b101e..db96fcbe2e 100644 --- a/tests/unit_tests/test_deplete_decay.py +++ b/tests/unit_tests/test_deplete_decay.py @@ -82,3 +82,8 @@ def test_deplete_decay_step_fissionable(run_in_tmpdir): _, u238 = results.get_atoms(f"{mat.id}", "U238") assert u238[1] == pytest.approx(original_atoms) + + # Check that material name is preserved in depletion results + step_result = results[0] + mat_from_results = step_result.get_material(f"{mat.id}") + assert mat_from_results.name == 'I do not decay.' diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 6463eaa208..558cb434be 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -59,8 +59,9 @@ def test_results_save(run_in_tmpdir): burn_list = full_burn_list[2*comm.rank: 2*comm.rank + 2] nuc_list = ["na", "nb"] + name_list = {mat: "" for mat in full_burn_list} op.get_results_info.return_value = ( - vol_dict, nuc_list, burn_list, full_burn_list) + vol_dict, nuc_list, burn_list, full_burn_list, name_list) # Construct end-of-step concentrations x1 = [rng.random(2), rng.random(2)] @@ -130,7 +131,8 @@ def test_results_save_without_rates(run_in_tmpdir): vol_dict = {"0": 1.0} nuc_list = ["na"] burn_list = ["0"] - op.get_results_info.return_value = (vol_dict, nuc_list, burn_list, burn_list) + name_list = {mat: "" for mat in burn_list} + op.get_results_info.return_value = (vol_dict, nuc_list, burn_list, burn_list, name_list) x = [np.array([1.0])] rates = ReactionRates(burn_list, nuc_list, ["ra"]) diff --git a/tests/unit_tests/test_element.py b/tests/unit_tests/test_element.py index d91cfaf6e2..82e526ec38 100644 --- a/tests/unit_tests/test_element.py +++ b/tests/unit_tests/test_element.py @@ -45,11 +45,11 @@ def test_expand_no_isotopes(): def test_expand_ta(): - ref = {'Ta180': 0.01201, 'Ta181': 99.98799} + ref = {'Ta181': 100.0} element = openmc.Element('Ta') for isotope in element.expand(100.0, 'ao'): assert isotope[1] == approx(ref[isotope[0]]) - + def test_expand_exceptions(): """ Test that correct exceptions are raised for invalid input """ diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 9aca8b5965..c5855a7b05 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -1,6 +1,8 @@ from math import pi from tempfile import TemporaryDirectory from pathlib import Path +import itertools +import random import h5py import numpy as np @@ -691,6 +693,49 @@ def test_mesh_material_volumes_serialize(): assert new_volumes.by_element(3) == [(2, 1.0)] +def test_mesh_material_volumes_serialize_with_bboxes(): + materials = np.array([ + [1, -1, -2], + [-1, -2, -2], + [2, 1, -2], + [2, -2, -2] + ]) + volumes = np.array([ + [0.5, 0.5, 0.0], + [1.0, 0.0, 0.0], + [0.5, 0.5, 0.0], + [1.0, 0.0, 0.0] + ]) + + # (xmin, ymin, zmin, xmax, ymax, zmax) + bboxes = np.empty((4, 3, 6)) + bboxes[..., 0:3] = np.inf + bboxes[..., 3:6] = -np.inf + bboxes[0, 0] = [-1.0, -2.0, -3.0, 1.0, 2.0, 3.0] # material 1 + bboxes[0, 1] = [-5.0, -6.0, -7.0, 5.0, 6.0, 7.0] # void + bboxes[1, 0] = [0.0, 0.0, 0.0, 10.0, 1.0, 2.0] # void + bboxes[2, 0] = [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0] # material 2 + bboxes[2, 1] = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0] # material 1 + bboxes[3, 0] = [-2.0, -2.0, -2.0, 2.0, 2.0, 2.0] # material 2 + + mmv = openmc.MeshMaterialVolumes(materials, volumes, bboxes) + with TemporaryDirectory() as tmpdir: + path = f'{tmpdir}/volumes_bboxes.npz' + mmv.save(path) + loaded = openmc.MeshMaterialVolumes.from_npz(path) + + assert loaded.has_bounding_boxes + first = loaded.by_element(0, include_bboxes=True)[0][2] + assert isinstance(first, openmc.BoundingBox) + np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0)) + np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0)) + + second = loaded.by_element(0, include_bboxes=True)[1][2] + assert isinstance(second, openmc.BoundingBox) + np.testing.assert_array_equal(second.lower_left, (-5.0, -6.0, -7.0)) + np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0)) + + def test_mesh_material_volumes_boundary_conditions(sphere_model): """Test the material volumes method using a regular mesh that overlaps with a vacuum boundary condition.""" @@ -718,6 +763,53 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model): assert evaluated[1] == pytest.approx(expected[1], rel=1e-2) +def test_mesh_material_volumes_bounding_boxes(): + # Create a model with 8 spherical cells at known locations with random radii + box = openmc.model.RectangularParallelepiped( + -10, 10, -10, 10, -10, 10, boundary_type='vacuum') + + mat = openmc.Material() + mat.add_nuclide('H1', 1.0) + + sph_cells = [] + for x, y, z in itertools.product((-5., 5.), repeat=3): + mat_i = mat.clone() + sph = openmc.Sphere(x, y, z, r=random.uniform(0.5, 1.5)) + sph_cells.append(openmc.Cell(region=-sph, fill=mat_i)) + background = openmc.Cell(region=-box & openmc.Intersection([~c.region for c in sph_cells])) + + model = openmc.Model() + model.geometry = openmc.Geometry(sph_cells + [background]) + model.settings.particles = 1000 + model.settings.batches = 10 + + # Create a one-element mesh that encompasses the entire geometry + mesh = openmc.RegularMesh() + mesh.lower_left = (-10., -10., -10.) + mesh.upper_right = (10., 10., 10.) + mesh.dimension = (1, 1, 1) + + # Run material volume calculation with bounding boxes + n_samples = (400, 400, 400) + mmv = mesh.material_volumes(model, n_samples, max_materials=10, bounding_boxes=True) + assert mmv.has_bounding_boxes + + # Create a mapping of material ID to bounding box + bbox_by_mat = { + mat_id: bbox + for mat_id, vol, bbox in mmv.by_element(0, include_bboxes=True) + if mat_id is not None and vol > 0.0 + } + + # Match the mesh ray spacing used for the bounding box estimator. + tol = 0.5 * mesh.bounding_box.width[0] / n_samples[0] + for cell in sph_cells: + bbox = bbox_by_mat[cell.fill.id] + cell_bbox = cell.bounding_box + np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol) + np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol) + + def test_raytrace_mesh_infinite_loop(run_in_tmpdir): # Create a model with one large spherical cell sphere = openmc.Sphere(r=100, boundary_type='vacuum') diff --git a/tests/unit_tests/test_particle_type.py b/tests/unit_tests/test_particle_type.py new file mode 100644 index 0000000000..67e385ea63 --- /dev/null +++ b/tests/unit_tests/test_particle_type.py @@ -0,0 +1,298 @@ +"""Unit tests for ParticleType class.""" + +import pytest +from openmc import ParticleType + + +# Tests for creating ParticleType instances + +def test_create_from_int(): + """Test creation from PDG number.""" + p = ParticleType(2112) + assert p.pdg_number == 2112 + assert int(p) == 2112 + + +def test_create_from_string_name(): + """Test creation from particle name.""" + p = ParticleType('neutron') + assert p.pdg_number == 2112 + + p = ParticleType('photon') + assert p.pdg_number == 22 + + +def test_create_from_string_aliases(): + """Test creation from particle aliases.""" + assert ParticleType('n').pdg_number == 2112 + assert ParticleType('gamma').pdg_number == 22 + assert ParticleType('p').pdg_number == 2212 + assert ParticleType('proton').pdg_number == 2212 + assert ParticleType('d').pdg_number == 1000010020 + assert ParticleType('t').pdg_number == 1000010030 + + +def test_create_from_string_nuclide(): + """Test creation from GNDS nuclide name.""" + p = ParticleType('He4') + assert p.pdg_number == 1000020040 + + p = ParticleType('U235') + assert p.pdg_number == 1000922350 + + p = ParticleType('Am242_m1') + assert p.pdg_number == 1000952421 + + +def test_create_from_string_pdg_prefix(): + """Test creation with pdg: prefix.""" + p = ParticleType('pdg:2112') + assert p.pdg_number == 2112 + + p = ParticleType('PDG:22') + assert p.pdg_number == 22 + + +def test_create_from_particle_type(): + """Test creation from existing ParticleType.""" + p1 = ParticleType(2112) + p2 = ParticleType(p1) + assert p1 == p2 + assert p1 is not p2 # Different instances + + +def test_legacy_particle_indices(): + """Test backward compatibility with legacy indices 0-3.""" + assert ParticleType(0) == ParticleType.NEUTRON + assert ParticleType(1) == ParticleType.PHOTON + assert ParticleType(2) == ParticleType.ELECTRON + assert ParticleType(3) == ParticleType.POSITRON + + +def test_create_invalid_type(): + """Test creation with invalid type raises TypeError.""" + with pytest.raises(TypeError): + ParticleType([2112]) + + with pytest.raises(TypeError): + ParticleType({'pdg': 2112}) + + +def test_create_invalid_string(): + """Test creation with invalid string raises ValueError.""" + with pytest.raises(ValueError): + ParticleType('') + + with pytest.raises(ValueError): + ParticleType('pdg:invalid') + + +def test_create_case_insensitive(): + """Test that string parsing is case insensitive for aliases.""" + assert ParticleType('NEUTRON').pdg_number == 2112 + assert ParticleType('Neutron').pdg_number == 2112 + assert ParticleType('PHOTON').pdg_number == 22 + + +# Tests for equality and comparison + +def test_equality_same_pdg(): + """Test that instances with same PDG are equal.""" + p1 = ParticleType(2112) + p2 = ParticleType(2112) + assert p1 == p2 + + +def test_equality_int(): + """Test equality comparison with int.""" + p = ParticleType(2112) + assert p == 2112 + assert 2112 == p + + +def test_equality_string(): + """Test equality comparison with string.""" + p = ParticleType(2112) + assert p == 'neutron' + assert p == 'pdg:2112' + assert p == 'n' + + +def test_inequality(): + """Test inequality comparisons.""" + p1 = ParticleType(2112) + p2 = ParticleType(22) + assert p1 != p2 + assert p1 != 22 + assert p1 != 'photon' + + +def test_equality_with_constants(): + """Test equality with class constants.""" + p = ParticleType(2112) + assert p == ParticleType.NEUTRON + assert ParticleType.NEUTRON == p + + +def test_equality_invalid_string(): + """Test equality with invalid string returns False.""" + p = ParticleType(2112) + assert not (p == 'invalid_particle') + assert p != 'invalid_particle' + + +# Tests for hashing behavior + +def test_hash_consistency(): + """Test that equal instances hash to the same value.""" + p1 = ParticleType(2112) + p2 = ParticleType(2112) + assert hash(p1) == hash(p2) + + +def test_set_deduplication(): + """Test that equal instances deduplicate in sets.""" + p1 = ParticleType(2112) + p2 = ParticleType(2112) + s = {p1, p2} + assert len(s) == 1 + + +def test_dict_key(): + """Test use as dictionary key.""" + p1 = ParticleType(2112) + p2 = ParticleType(2112) + d = {p1: 'neutron'} + assert d[p2] == 'neutron' + + +def test_hash_different_particles(): + """Test that different particles have different hashes (usually).""" + p1 = ParticleType(2112) + p2 = ParticleType(22) + # Different PDG numbers should (almost always) have different hashes + assert hash(p1) != hash(p2) + + +# Tests for properties and computed attributes + +def test_pdg_number_property(): + """Test pdg_number property.""" + p = ParticleType(2112) + assert p.pdg_number == 2112 + assert isinstance(p.pdg_number, int) + + +def test_int_conversion(): + """Test __int__ conversion.""" + p = ParticleType(1000020040) + assert int(p) == 1000020040 + + +def test_zam_elementary(): + """Test zam property for elementary particles.""" + assert ParticleType.NEUTRON.zam is None + assert ParticleType.PHOTON.zam is None + assert ParticleType.ELECTRON.zam is None + assert ParticleType.POSITRON.zam is None + + +def test_zam_nucleus(): + """Test zam property for nuclear particles.""" + he4 = ParticleType('He4') + assert he4.zam == (2, 4, 0) + + u235 = ParticleType('U235') + Z, A, m = u235.zam + assert Z == 92 + assert A == 235 + assert m == 0 + + +def test_zam_metastable(): + """Test zam property for metastable nuclei.""" + # Am242m has m=1 + am242m = ParticleType('Am242_m1') + Z, A, m = am242m.zam + assert Z == 95 + assert A == 242 + assert m == 1 + + +def test_is_nucleus_false(): + """Test is_nucleus for elementary particles.""" + assert not ParticleType.NEUTRON.is_nucleus + assert not ParticleType.PHOTON.is_nucleus + assert not ParticleType.ELECTRON.is_nucleus + assert not ParticleType.POSITRON.is_nucleus + + +def test_is_nucleus_true(): + """Test is_nucleus for nuclear particles.""" + assert ParticleType.ALPHA.is_nucleus + assert ParticleType('He4').is_nucleus + assert ParticleType('U235').is_nucleus + assert ParticleType.DEUTERON.is_nucleus + assert ParticleType.TRITON.is_nucleus + + +# Tests for __str__ and __repr__ + +def test_str_elementary(): + """Test string representation of elementary particles.""" + assert str(ParticleType.NEUTRON) == 'neutron' + assert str(ParticleType.PHOTON) == 'photon' + assert str(ParticleType.ELECTRON) == 'electron' + assert str(ParticleType.POSITRON) == 'positron' + assert str(ParticleType.PROTON) == 'H1' + + +def test_str_nucleus(): + """Test string representation of nuclei.""" + assert str(ParticleType.ALPHA) == 'He4' + assert str(ParticleType('U235')) == 'U235' + assert str(ParticleType.DEUTERON) == 'H2' + assert str(ParticleType.TRITON) == 'H3' + + +def test_str_arbitrary_pdg(): + """Test string representation of arbitrary PDG number.""" + # PDG number that doesn't match any known particle + p = ParticleType(12345) + assert str(p) == 'pdg:12345' + + +def test_repr(): + """Test repr includes PDG number.""" + p = ParticleType.NEUTRON + repr_str = repr(p) + assert 'ParticleType' in repr_str + assert 'PDG=2112' in repr_str + assert 'neutron' in repr_str + + +def test_repr_nucleus(): + """Test repr for nuclear particles.""" + p = ParticleType.ALPHA + repr_str = repr(p) + assert 'ParticleType' in repr_str + assert 'PDG=1000020040' in repr_str + assert 'He4' in repr_str + + +def test_create_from_numpy_int(): + """Test creation from numpy integer types.""" + import numpy as np + p = ParticleType(np.int32(2112)) + assert p.pdg_number == 2112 + + p = ParticleType(np.int64(22)) + assert p.pdg_number == 22 + + +def test_equality_numpy_int(): + """Test equality comparison with numpy integer types.""" + import numpy as np + p = ParticleType(2112) + assert p == np.int32(2112) + assert p == np.int64(2112) diff --git a/tests/unit_tests/test_source_file.py b/tests/unit_tests/test_source_file.py index 41906c80f8..3e84fbe517 100644 --- a/tests/unit_tests/test_source_file.py +++ b/tests/unit_tests/test_source_file.py @@ -42,7 +42,7 @@ def test_source_file(run_in_tmpdir): assert np.all(arr['E'] == n - np.arange(n)) assert np.all(arr['wgt'] == 1.0) assert np.all(arr['delayed_group'] == 0) - assert np.all(arr['particle'] == 0) + assert np.all(arr['particle'] == 2112) # PDG number for neutron # Ensure sites read in are consistent sites = openmc.ParticleList.from_hdf5('test_source.h5') @@ -64,7 +64,7 @@ def test_source_file(run_in_tmpdir): dgs = np.array([s.delayed_group for s in sites]) assert np.all(dgs == 0) p_types = np.array([s.particle for s in sites]) - assert np.all(p_types == 0) + assert np.all(p_types == 2112) # PDG number for neutron # Ensure a ParticleList item is a SourceParticle site = sites[0] diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index e92e9e1352..507e85743d 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -633,12 +633,27 @@ def test_combine_distributions(): assert len(mixed.distribution) == 2 assert len(mixed.probability) == 2 + # Single tabular returns a tabular distribution with scaled probabilities + t_single = openmc.stats.Tabular([0.0, 1.0], [2.0, 0.0]) + scaled = openmc.stats.combine_distributions([t_single], [0.25]) + assert isinstance(scaled, openmc.stats.Tabular) + assert scaled.p == pytest.approx([0.5, 0.0]) + + # Mixture with biased tabular should preserve unbiased mean via weights + bias = openmc.stats.Tabular([0.0, 1.0], [2.0, 0.0]) + t_biased = openmc.stats.Tabular([0.0, 1.0], [1.0, 1.0], bias=bias) + d1 = openmc.stats.delta_function(0.0) + mixed = openmc.stats.combine_distributions([t_biased, d1], [0.5, 0.5]) + assert isinstance(mixed, openmc.stats.Mixture) + samples, weights = mixed.sample(10_000) + assert_sample_mean(samples*weights, 0.25) + # Combine 1 discrete and 2 tabular -- the tabular distributions should # combine to produce a uniform distribution with mean 0.5. The combined # distribution should have a mean of 0.25. t1 = openmc.stats.Tabular([0., 1.], [2.0, 0.0]) t2 = openmc.stats.Tabular([0., 1.], [0.0, 2.0]) - d1 = openmc.stats.Discrete([0.0], [1.0]) + d1 = openmc.stats.delta_function(0.0) combined = openmc.stats.combine_distributions([t1, t2, d1], [0.25, 0.25, 0.5]) assert combined.integral() == pytest.approx(1.0) diff --git a/tests/unit_tests/test_tracks.py b/tests/unit_tests/test_tracks.py index fa86641598..74da83552f 100644 --- a/tests/unit_tests/test_tracks.py +++ b/tests/unit_tests/test_tracks.py @@ -70,7 +70,7 @@ def test_tracks(sphere_model, particle, run_in_tmpdir): particle_track = track.particle_tracks[0] assert isinstance(particle_track, openmc.ParticleTrack) - assert particle_track.particle.name.lower() == particle + assert str(particle_track.particle).lower() == particle assert isinstance(particle_track.states, np.ndarray) # Sanity checks on actual data @@ -140,8 +140,10 @@ def test_filter(sphere_model, run_in_tmpdir): assert matches == tracks matches = tracks.filter(state_filter=lambda s: s['E'] > 0.0) assert matches == tracks - matches = tracks.filter(particle='bunnytron') + matches = tracks.filter(particle='proton') assert matches == [] + with pytest.raises(ValueError): + tracks.filter(particle='bunnytron') def test_write_to_vtk(sphere_model): diff --git a/tests/unit_tests/test_universe.py b/tests/unit_tests/test_universe.py index efe8552a64..6fbf1cc383 100644 --- a/tests/unit_tests/test_universe.py +++ b/tests/unit_tests/test_universe.py @@ -1,3 +1,5 @@ +from pathlib import Path + import lxml.etree as ET import numpy as np import openmc @@ -104,6 +106,43 @@ def test_plot(run_in_tmpdir, sphere_model): plt.close('all') +def test_mg_plot(run_in_tmpdir): + # Create a simple universe with macroscopic data + h2o_data = openmc.Macroscopic('LWTR') + water = openmc.Material(name='Water') + water.set_density('macro', 1.0) + water.add_macroscopic(h2o_data) + sph = openmc.Sphere(r=10, boundary_type="vacuum") + cell = openmc.Cell(region=-sph, fill=water) + univ = openmc.Universe(cells=[cell]) + + # Create MGXS library and export to HDF5 + groups = openmc.mgxs.EnergyGroups([1e-5, 20.0e6]) + h2o_xsdata = openmc.XSdata('LWTR', groups) + h2o_xsdata.order = 0 + h2o_xsdata.set_total([1.0]) + h2o_xsdata.set_absorption([0.5]) + scatter_matrix = np.array([[[0.5]]]) + scatter_matrix = np.rollaxis(scatter_matrix, 0, 3) + h2o_xsdata.set_scatter_matrix(scatter_matrix) + mg_library = openmc.MGXSLibrary(groups) + mg_library.add_xsdatas([h2o_xsdata]) + mgxs_path = Path.cwd() / 'mgxs.h5' + mg_library.export_to_hdf5(mgxs_path) + + # Set MG cross sections in config and plot + with openmc.config.patch('mg_cross_sections', mgxs_path): + univ.plot(width=(200, 200), basis='yz', color_by='cell') + univ.plot(width=(200, 200), basis='yz', color_by='material') + + with pytest.raises(RuntimeError): + univ.plot(width=(200, 200), basis='yz', color_by='cell') + + # Close plots to avoid warning + import matplotlib.pyplot as plt + plt.close('all') + + def test_get_nuclides(uo2): c = openmc.Cell(fill=uo2) univ = openmc.Universe(cells=[c]) diff --git a/tests/unit_tests/weightwindows/test_wwinp_reader.py b/tests/unit_tests/weightwindows/test_wwinp_reader.py index 637463bb2e..601517d47c 100644 --- a/tests/unit_tests/weightwindows/test_wwinp_reader.py +++ b/tests/unit_tests/weightwindows/test_wwinp_reader.py @@ -52,7 +52,7 @@ n_mesh.z_grid = np.array([-100.0, n_e_bounds = (np.array([0.0, 100000.0, 146780.0]),) -n_particles = ('neutron',) +n_particles = [openmc.ParticleType.NEUTRON] # expected results - neutron and photon data np_mesh = openmc.RectilinearMesh() @@ -63,7 +63,7 @@ np_mesh.z_grid = n_mesh.z_grid np_e_bounds = (np.array([0.0, 100000.0, 146780.0, 215440.0]), np.array([0.0, 1.0E8])) -np_particles = ('neutron', 'photon') +np_particles = [openmc.ParticleType.NEUTRON, openmc.ParticleType.PHOTON] # expected results - photon data only p_mesh = openmc.RectilinearMesh() @@ -74,7 +74,7 @@ p_mesh.y_grid = np_mesh.y_grid p_mesh.z_grid = np.array([-50.0, 50.0]) p_e_bounds = (np.array([0.0, 100000.0, 146780.0, 215440.0, 316230.0]),) -p_particles = ('photon',) +p_particles = [openmc.ParticleType.PHOTON] expected_results = [('wwinp_n', n_mesh, n_particles, n_e_bounds), ('wwinp_np', np_mesh, np_particles, np_e_bounds),