crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

fix conflicts

Browse source
This commit is contained in:
church89 2024-02-07 12:58:50 +01:00
commit ad8c7ed819
226 changed files with 5252 additions and 1641 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.github/pull_request_template.md vendored
View file

@ -13,7 +13,7 @@ Fixes # (issue)
# Checklist
- [ ] I have performed a self-review of my own code
- [ ] I have run [clang-format](https://docs.openmc.org/en/latest/devguide/styleguide.html#automatic-formatting) on any C++ source files (if applicable)
- [ ] I have run [clang-format](https://docs.openmc.org/en/latest/devguide/styleguide.html#automatic-formatting) (version 15) on any C++ source files (if applicable)
- [ ] I have followed the [style guidelines](https://docs.openmc.org/en/latest/devguide/styleguide.html#python) for Python source files (if applicable)
- [ ] I have made corresponding changes to the documentation (if applicable)
- [ ] I have added tests that prove my fix is effective or that my feature works (if applicable)

8
.github/workflows/ci.yml vendored
View file

@ -127,6 +127,14 @@ jobs:
echo "$HOME/NJOY2016/build" >> $GITHUB_PATH
$GITHUB_WORKSPACE/tools/ci/gha-install.sh
- name: cache-xs
uses: actions/cache@v3
with:
path: |
~/nndc_hdf5
~/endf-b-vii.1
key: ${{ runner.os }}-build-xs-cache
- name: before
shell: bash
run: $GITHUB_WORKSPACE/tools/ci/gha-before-script.sh

9
.github/workflows/format-check.yml vendored
View file

@ -1,6 +1,13 @@
name: C++ Format Check
on: pull_request
on:
# allow workflow to be run manually
workflow_dispatch:
pull_request:
branches:
- develop
- master
jobs:
cpp-linter:

13
CMakeLists.txt
View file

@ -3,8 +3,8 @@ project(openmc C CXX)
# Set version numbers
set(OPENMC_VERSION_MAJOR 0)
set(OPENMC_VERSION_MINOR 13)
set(OPENMC_VERSION_RELEASE 4)
set(OPENMC_VERSION_MINOR 14)
set(OPENMC_VERSION_RELEASE 1)
set(OPENMC_VERSION ${OPENMC_VERSION_MAJOR}.${OPENMC_VERSION_MINOR}.${OPENMC_VERSION_RELEASE})
configure_file(include/openmc/version.h.in "${CMAKE_BINARY_DIR}/include/openmc/version.h" @ONLY)
@ -120,6 +120,7 @@ if(OPENMC_USE_DAGMC)
message(FATAL_ERROR "Discovered DAGMC Version: ${DAGMC_VERSION}. \
Please update DAGMC to version 3.2.0 or greater.")
endif()
message(STATUS "Found DAGMC: ${DAGMC_DIR} (version ${DAGMC_VERSION})")
endif()
#===============================================================================
@ -154,6 +155,11 @@ if(NOT DEFINED HDF5_PREFER_PARALLEL)
endif()
find_package(HDF5 REQUIRED COMPONENTS C HL)
# Remove HDF5 transitive dependencies that are system libraries
list(FILTER HDF5_LIBRARIES EXCLUDE REGEX ".*lib(pthread|dl|m).*")
message(STATUS "HDF5 Libraries: ${HDF5_LIBRARIES}")
if(HDF5_IS_PARALLEL)
if(NOT OPENMC_USE_MPI)
message(FATAL_ERROR "Parallel HDF5 was detected, but MPI was not enabled.\
@ -405,6 +411,7 @@ list(APPEND libopenmc_SOURCES
src/tallies/filter_energyfunc.cpp
src/tallies/filter_legendre.cpp
src/tallies/filter_material.cpp
src/tallies/filter_materialfrom.cpp
src/tallies/filter_mesh.cpp
src/tallies/filter_meshsurface.cpp
src/tallies/filter_mu.cpp
@ -490,7 +497,7 @@ endif()
# target_link_libraries treats any arguments starting with - but not -l as
# linker flags. Thus, we can pass both linker flags and libraries together.
target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}
xtensor gsl::gsl-lite-v1 fmt::fmt)
xtensor gsl::gsl-lite-v1 fmt::fmt ${CMAKE_DL_LIBS})
if(TARGET pugixml::pugixml)
target_link_libraries(libopenmc pugixml::pugixml)

2
LICENSE
View file

@ -1,4 +1,4 @@
Copyright (c) 2011-2023 Massachusetts Institute of Technology, UChicago Argonne
Copyright (c) 2011-2024 Massachusetts Institute of Technology, UChicago Argonne
LLC, and OpenMC contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of

2
README.md
View file

@ -1,7 +1,7 @@
# OpenMC Monte Carlo Particle Transport Code
[![License](https://img.shields.io/badge/license-MIT-green)](https://docs.openmc.org/en/latest/license.html)
[![GitHub Actions build status (Linux)](https://github.com/openmc-dev/openmc/workflows/CI/badge.svg?branch=develop)](https://github.com/openmc-dev/openmc/actions?query=workflow%3ACI)
[![GitHub Actions build status (Linux)](https://github.com/openmc-dev/openmc/actions/workflows/ci.yml/badge.svg?branch=develop)](https://github.com/openmc-dev/openmc/actions/workflows/ci.yml)
[![Code Coverage](https://coveralls.io/repos/github/openmc-dev/openmc/badge.svg?branch=develop)](https://coveralls.io/github/openmc-dev/openmc?branch=develop)
[![dockerhub-publish-develop-dagmc](https://github.com/openmc-dev/openmc/workflows/dockerhub-publish-develop-dagmc/badge.svg)](https://github.com/openmc-dev/openmc/actions?query=workflow%3Adockerhub-publish-develop-dagmc)
[![dockerhub-publish-develop](https://github.com/openmc-dev/openmc/workflows/dockerhub-publish-develop/badge.svg)](https://github.com/openmc-dev/openmc/actions?query=workflow%3Adockerhub-publish-develop)

6
docs/source/conf.py
View file

@ -62,16 +62,16 @@ master_doc = 'index'
# General information about the project.
project = 'OpenMC'
copyright = '2011-2023, Massachusetts Institute of Technology, UChicago Argonne LLC, and OpenMC contributors'
copyright = '2011-2024, Massachusetts Institute of Technology, UChicago Argonne LLC, and OpenMC contributors'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = "0.13"
version = "0.14"
# The full version, including alpha/beta/rc tags.
release = "0.13.4"
release = "0.14.1-dev"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.

4
docs/source/devguide/styleguide.rst
View file

@ -29,6 +29,10 @@ whenever a file is saved. For example, `Visual Studio Code
<https://code.visualstudio.com/docs/cpp/cpp-ide#_code-formatting>`_ includes
support for running clang-format.
.. note::
OpenMC's CI uses `clang-format` version 15. A different version of `clang-format`
may produce different line changes and as a result fail the CI test.
Miscellaneous
-------------

31
docs/source/devguide/tests.rst
View file

@ -59,6 +59,31 @@ If you want to view testing output on failure run::
ctest --output-on-failure
Possible Reasons for Test Failures
----------------------------------
You may find that when you run the test suite, not everything passes. First,
make sure you have satisfied all the prerequisites above. After you have done
that, consider the following:
- When building OpenMC, make sure you run CMake with
``-DCMAKE_BUILD_TYPE=Debug``. Building with a release build will result in
some test failures due to differences in which compiler optimizations are
used.
- Because tallies involve the sum of many floating point numbers, the
non-associativity of floating point numbers can result in different answers
especially when the number of threads is high (different order of operations).
Thus, if you are running on a CPU with many cores, you may need to limit the
number of OpenMP threads used. It is recommended to set the
:envvar:`OMP_NUM_THREADS` environment variable to 2.
- Recent versions of NumPy use instruction dispatch that may generate different
results depending the particular ISA that you are running on. To avoid issues,
you may need to disable AVX512 instructions. This can be done by setting the
:envvar:`NPY_DISABLE_CPU_FEATURES` environment variable to "AVX512F
AVX512_SKX". When NumPy/SciPy are built against OpenBLAS, you may also need to
limit the number of threads that OpenBLAS uses internally; this can be done by
setting the :envvar:`OPENBLAS_NUM_THREADS` environment variable to 1.
Generating XML Inputs
---------------------
@ -109,6 +134,12 @@ following files to your new test directory:
compiler options during openmc configuration and build (e.g., no MPI, no
debug/optimization).
For tests using the Python API, both the **inputs_true.dat** and
**results_true.dat** files can be generated automatically in the correct format
via::
pytest --update <name-of-test>
In addition to this description, please see the various types of tests that are
already included in the test suite to see how to create them. If all is
implemented correctly, the new test will automatically be discovered by pytest.

89
docs/source/io_formats/settings.rst
View file

@ -276,6 +276,16 @@ then, OpenMC will only use up to the :math:`P_1` data.
.. note:: This element is not used in the continuous-energy
:ref:`energy_mode`.
--------------------------------------
``<max_write_lost_particles>`` Element
--------------------------------------
This ``<max_write_lost_particles>`` element indicates the maximum number of
particle restart files (per MPI process) to write for lost particles.
*Default*: None
.. _mesh_element:
------------------
@ -463,6 +473,8 @@ pseudo-random number generator.
*Default*: 1
.. _source_element:
--------------------
``<source>`` Element
--------------------
@ -481,7 +493,8 @@ attributes/sub-elements:
*Default*: 1.0
:type:
Indicator of source type. One of ``independent``, ``file``, or ``compiled``.
Indicator of source type. One of ``independent``, ``file``, ``compiled``, or ``mesh``.
The type of the source will be determined by this attribute if it is present.
:particle:
The source particle type, either ``neutron`` or ``photon``.
@ -654,6 +667,14 @@ attributes/sub-elements:
*Default*: false
:mesh:
For mesh sources, this indicates the ID of the corresponding mesh.
:source:
For mesh sources, this sub-element specifies the source for an individual
mesh element and follows the format for :ref:`source_element`. The number of
``<source>`` sub-elements should correspond to the number of mesh elements.
.. _univariate:
Univariate Probability Distributions
@ -1058,8 +1079,14 @@ The ``<volume_calc>`` element indicates that a stochastic volume calculation
should be run at the beginning of the simulation. This element has the following
sub-elements/attributes:
:cells:
The unique IDs of cells for which the volume should be estimated.
:domain_type:
The type of each domain for the volume calculation ("cell", "material", or
"universe").
*Default*: None
:domain_ids:
The unique IDs of domains for which the volume should be estimated.
*Default*: None
@ -1069,16 +1096,41 @@ sub-elements/attributes:
*Default*: None
:lower_left:
The lower-left Cartesian coordinates of a bounding box that is used to
sample points within.
The lower-left Cartesian coordinates of a bounding box that is used to
sample points within.
*Default*: None
*Default*: None
:upper_right:
The upper-right Cartesian coordinates of a bounding box that is used to
sample points within.
The upper-right Cartesian coordinates of a bounding box that is used to
sample points within.
*Default*: None
*Default*: None
:threshold:
Presence of a ``<threshold>`` sub-element indicates that the volume
calculation will be halted based on a threshold on the error. It has the
following sub-elements/attributes:
:type:
The type of the trigger. Accepted options are "variance", "std_dev",
and "rel_err".
:variance:
Variance of the mean, :math:`\sigma^2`
:std_dev:
Standard deviation of the mean, :math:`\sigma`
:rel_err:
Relative error of the mean, :math:`\frac{\sigma}{\mu}`
*Default*: None
:threshold:
The trigger's convergence criterion for the given type.
*Default*: None
----------------------------
``<weight_windows>`` Element
@ -1201,6 +1253,25 @@ mesh-based weight windows.
*Default*: 5.0
---------------------------------------
``<weight_window_checkpoints>`` Element
---------------------------------------
The ``<weight_window_checkpoints>`` element indicates the checkpoints for weight
window split/roulette (surface, collision or both). This element has the
following sub-elements/attributes:
:surface:
If set to "true", weight window checks will be performed at surface
crossings.
*Default*: False
:collision:
If set to "true", weight window checks will be performed at collisions.
*Default*: True
--------------------------------------
``<weight_windows_file>`` Element
--------------------------------------

8
docs/source/io_formats/summary.rst
View file

@ -60,9 +60,11 @@ The current version of the summary file format is 6.0.
- **coefficients** (*double[]*) -- Array of coefficients that define
the surface. See :ref:`surface_element` for what coefficients are
defined for each surface type.
- **boundary_condition** (*char[]*) -- Boundary condition applied to
the surface. Can be 'transmission', 'vacuum', 'reflective', or
'periodic'.
- **boundary_type** (*char[]*) -- Boundary condition applied to
the surface. Can be 'transmission', 'vacuum', 'reflective',
'periodic', or 'white'.
- **albedo** (*double*) -- Boundary albedo as a positive multiplier
of particle weight. If absent, it is assumed to be 1.0.
- **geom_type** (*char[]*) -- Type of geometry used to create the cell.
Either 'csg' or 'dagmc'.

2
docs/source/license.rst
View file

@ -4,7 +4,7 @@
License Agreement
=================
Copyright © 2011-2023 Massachusetts Institute of Technology, UChicago Argonne
Copyright © 2011-2024 Massachusetts Institute of Technology, UChicago Argonne
LLC, and OpenMC contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of

3
docs/source/pythonapi/base.rst
View file

@ -25,6 +25,7 @@ Simulation Settings
openmc.IndependentSource
openmc.FileSource
openmc.CompiledSource
openmc.MeshSource
openmc.SourceParticle
openmc.VolumeCalculation
openmc.Settings
@ -118,6 +119,7 @@ Constructing Tallies
openmc.Filter
openmc.UniverseFilter
openmc.MaterialFilter
openmc.MaterialFromFilter
openmc.CellFilter
openmc.CellFromFilter
openmc.CellBornFilter
@ -160,6 +162,7 @@ Geometry Plotting
:template: myclass.rst
openmc.Plot
openmc.ProjectionPlot
openmc.Plots
Running OpenMC

3
docs/source/pythonapi/capi.rst
View file

@ -54,11 +54,14 @@ Classes
:template: myclass.rst
Cell
CylindricalMesh
EnergyFilter
MaterialFilter
Material
MeshFilter
MeshSurfaceFilter
Nuclide
RectilinearMesh
RegularMesh
SphericalMesh
Tally

4
docs/source/pythonapi/model.rst
View file

@ -11,8 +11,6 @@ Convenience Functions
:template: myfunction.rst
openmc.model.borated_water
openmc.model.hexagonal_prism
openmc.model.rectangular_prism
openmc.model.subdivide
openmc.model.pin
@ -26,9 +24,11 @@ Composite Surfaces
openmc.model.CruciformPrism
openmc.model.CylinderSector
openmc.model.HexagonalPrism
openmc.model.IsogonalOctagon
openmc.model.Polygon
openmc.model.RectangularParallelepiped
openmc.model.RectangularPrism
openmc.model.RightCircularCylinder
openmc.model.XConeOneSided
openmc.model.YConeOneSided

284
docs/source/releasenotes/0.14.0.rst Normal file
View file

@ -0,0 +1,284 @@
====================
What's New in 0.14.0
====================
.. currentmodule:: openmc
-------
Summary
-------
This release of OpenMC includes many bug fixes, performance improvements, and
several notable new features. Some of the highlights include projection plots,
pulse height tallies for photons, weight window generation, and an ability to
specify continuous removal or feed of nuclides/elements during depletion.
Additionally, one of the longstanding annoyances of depletion calculations,
namely the need to include initial "dilute" nuclides, has been eliminated. There
are also a wide array of general improvements in the Python API.
------------------------------------
Compatibility Notes and Deprecations
------------------------------------
- The :class:`openmc.deplete.MicroXS` has been completely redesigned and
improved. See further comments below under "New Features". (`#2572
<https://github.com/openmc-dev/openmc/pull/2572>`_, `#2579
<https://github.com/openmc-dev/openmc/pull/2579>`_, `#2595
<https://github.com/openmc-dev/openmc/pull/2595>`_, `#2700
<https://github.com/openmc-dev/openmc/pull/2700>`_)
- The ``rectangular_prism`` function has been replaced by the
:class:`openmc.model.RectangularPrism` class and the ``hexagonal_prism``
function has been replaced by the :class:`openmc.model.HexagonalPrism` class.
Note that whereas the ``rectangular_prism`` and ``hexagonal_prism`` functions
returned a region representing the interior of the prism, the new
:class:`~openmc.model.RectangularPrism` and
:class:`~openmc.model.HexagonalPrism` classes return composite surfaces, so
you need to use the unary ``-`` or ``+`` operators to obtain a region that can
be assigned to a cell. (`#2739
<https://github.com/openmc-dev/openmc/pull/2739>`_)
- The ``Source`` class has been refactored and split up into three separate
classes: :class:`~openmc.IndependentSource`, :class:`~openmc.FileSource`, and
:class:`~openmc.CompiledSource`. (`#2524
<https://github.com/openmc-dev/openmc/pull/2524>`_)
- The ``vertices`` and ``centroids`` attributes on mesh classes now always
return Cartesian coordinates and the shape of the returned arrays has changed
to allow `ijk` indexing using a tuple (i.e., `xyz = vertices[i, j, k]`).
(`#2711 <https://github.com/openmc-dev/openmc/pull/2711>`_)
- The :attr:`openmc.Material.decay_photon_energy` attribute has been replaced by
the :meth:`openmc.Material.get_decay_photon_energy` method. (`#2715
<https://github.com/openmc-dev/openmc/pull/2715>`_)
------------
New Features
------------
- A new :class:`openmc.ProjectionPlot` class enables the generation of orthographic or
perspective projection plots. (`#1926
<https://github.com/openmc-dev/openmc/pull/1926>`_)
- The :class:`openmc.model.RightCircularCylinder` class now supports optional
filleted edges. (`#2309 <https://github.com/openmc-dev/openmc/pull/2309>`_)
- Continuous removal or feed of nuclides/elements between materials can now be
modeled during depletion via the
:meth:`openmc.deplete.abc.Integrator.add_transfer_rate` method. (`#2358
<https://github.com/openmc-dev/openmc/pull/2358>`_, `#2564
<https://github.com/openmc-dev/openmc/pull/2564>`_, `#2626
<https://github.com/openmc-dev/openmc/pull/2626>`_)
- The MAGIC method for global weight window generation has been implemented as
part of the C++ API. (`#2359
<https://github.com/openmc-dev/openmc/pull/2359>`_)
- A new capability for pulse height tallies (currently limited to photons) has
been added and can be used via the "pulse-height" tally score. (`#2452
<https://github.com/openmc-dev/openmc/pull/2452>`_)
- A :class:`openmc.model.CruciformPrism` class has been added that provides a
generalized cruciform prism composite surface. (`#2457
<https://github.com/openmc-dev/openmc/pull/2457>`_)
- Type hints have been added in various places throughout the Python API.
(`#2462 <https://github.com/openmc-dev/openmc/pull/2462>`_, `#2467
<https://github.com/openmc-dev/openmc/pull/2467>`_, `#2468
<https://github.com/openmc-dev/openmc/pull/2468>`_, `#2470
<https://github.com/openmc-dev/openmc/pull/2470>`_, `#2471
<https://github.com/openmc-dev/openmc/pull/2471>`_, `#2601
<https://github.com/openmc-dev/openmc/pull/2601>`_)
- Voxel plots can now be generated through the :meth:`openmc.Plot.to_vtk`
method. (`#2464 <https://github.com/openmc-dev/openmc/pull/2464>`_)
- The :class:`openmc.mgxs.EnergyGroups` class now allows you to alternatively
pass a string of the group structure name (e.g., "CCFE-709") instead of the
energy group boundaries. (`#2466
<https://github.com/openmc-dev/openmc/pull/2466>`_)
- Several enhancements have been made to the :meth:`openmc.Universe.plot` method
(addition of axis labels with units, ability to show legend and/or outlines, automatic
determination of origin/width, ability to pass total number of pixels).
(`#2472 <https://github.com/openmc-dev/openmc/pull/2472>`_, `#2482
<https://github.com/openmc-dev/openmc/pull/2482>`_, `#2483
<https://github.com/openmc-dev/openmc/pull/2483>`_, `#2492
<https://github.com/openmc-dev/openmc/pull/2492>`_, `#2513
<https://github.com/openmc-dev/openmc/pull/2513>`_, `#2575
<https://github.com/openmc-dev/openmc/pull/2575>`_)
- Functionality in the Python dealing with bounding boxes now relies on a new
:class:`openmc.BoundingBox` class. (`#2475
<https://github.com/openmc-dev/openmc/pull/2475>`_)
- Users now have more flexibility in specifying nuclides and reactions in the
:func:`openmc.plot_xs` function. (`#2478
<https://github.com/openmc-dev/openmc/pull/2478>`_)
- The import time of the :mod:`openmc` Python module has been improved by
deferring the import of matplotlib. (`#2488
<https://github.com/openmc-dev/openmc/pull/2488>`_)
- Mesh clases in the Python API now support a ``bounding_box`` property. (`#2507
<https://github.com/openmc-dev/openmc/pull/2507>`_, `#2620
<https://github.com/openmc-dev/openmc/pull/2620>`_, `#2621
<https://github.com/openmc-dev/openmc/pull/2621>`_)
- The ``Source`` class has been refactored and split up into three separate
classes: :class:`~openmc.IndependentSource`, :class:`~openmc.FileSource`, and
:class:`~openmc.CompiledSource`. (`#2524
<https://github.com/openmc-dev/openmc/pull/2524>`_)
- Support was added for curvilinear elements when exporting cylindrical and
spherical meshes to VTK. (`#2533
<https://github.com/openmc-dev/openmc/pull/2533>`_)
- The :class:`openmc.Tally` class now has a
:attr:`~openmc.Tally.multiply_density` attribute that indicates whether
reaction rate tallies should include the number density of the nuclide of
interest. (`#2539 <https://github.com/openmc-dev/openmc/pull/2539>`_)
- The :func:`~openmc.wwinp_to_wws` function now supports ``wwinp`` files with
cylindrical or spherical meshes. (`#2556
<https://github.com/openmc-dev/openmc/pull/2556>`_)
- Depletion no longer relies on adding initial "dilute" nuclides to each
depletable material in order to compute reaction rates. (`#2559
<https://github.com/openmc-dev/openmc/pull/2559>`_, `#2568
<https://github.com/openmc-dev/openmc/pull/2568>`_)
- The :class:`openmc.deplete.Results` class now has
:meth:`~openmc.deplete.Results.get_mass` (`#2565
<https://github.com/openmc-dev/openmc/pull/2565>`_),
:meth:`~openmc.deplete.Results.get_activity` (`#2617
<https://github.com/openmc-dev/openmc/pull/2617>`_), and
:meth:`~openmc.deplete.Results.get_decay_heat` (`#2625
<https://github.com/openmc-dev/openmc/pull/2625>`_) methods.
- The :meth:`openmc.deplete.StepResult.save` method now supports a ``path``
argument. (`#2567 <https://github.com/openmc-dev/openmc/pull/2567>`_)
- The :class:`openmc.deplete.MicroXS` has been completely redesigned and
improved. First, it no longer relies on the :mod:`openmc.mgxs` module, no
longer subclasses :class:`pandas.DataFrame`, and doesn't require adding
initial "dilute" nuclides into material compositions. It now enables users to
specify an energy group structure to collect multigroup cross sections,
specify nuclides/reactions, and works with mesh domains in addition to the
existing domains. A new :func:`openmc.deplete.get_microxs_and_flux` function
was added that improves the workflow for calculating microscopic cross
sections along with fluxes. Altogether, these changes make it straightforward
to switch between coupled and independent operators for depletion/activation
calculations. (`#2572 <https://github.com/openmc-dev/openmc/pull/2572>`_,
`#2579 <https://github.com/openmc-dev/openmc/pull/2579>`_, `#2595
<https://github.com/openmc-dev/openmc/pull/2595>`_, `#2700
<https://github.com/openmc-dev/openmc/pull/2700>`_)
- The :class:`openmc.Geometry` class now has ``merge_surfaces`` and
``surface_precision`` arguments. (`#2602
<https://github.com/openmc-dev/openmc/pull/2602>`_)
- Several predefined energy group structures have been added ("MPACT-51",
"MPACT-60", "MPACT-69", "SCALE-252"). (`#2614
<https://github.com/openmc-dev/openmc/pull/2614>`_)
- When running a depletion calculation, you are now allowed to include nuclides
in the initial material compositions that do not have neutron cross sections
(decay-only nuclides). (`#2616
<https://github.com/openmc-dev/openmc/pull/2616>`_)
- The :class:`~openmc.CylindricalMesh` and :class:`~openmc.SphericalMesh`
classes can now be fully formed using the constructor. (`#2619
<https://github.com/openmc-dev/openmc/pull/2619>`_)
- A time cutoff can now be specified in the :attr:`openmc.Settings.cutoff`
attribute. (`#2631 <https://github.com/openmc-dev/openmc/pull/2631>`_)
- The :meth:`openmc.Material.add_element` method now supports a
``cross_sections`` argument that allows a cross section data source to be
specified. (`#2633 <https://github.com/openmc-dev/openmc/pull/2633>`_)
- The :class:`~openmc.Cell` class now has a :meth:`~openmc.Cell.plot` method.
(`#2648 <https://github.com/openmc-dev/openmc/pull/2648>`_)
- The :class:`~openmc.Geometry` class now has a :meth:`~openmc.Geometry.plot`
method. (`#2661 <https://github.com/openmc-dev/openmc/pull/2661>`_)
- When weight window checks are performed can now be explicitly specified with
the :attr:`openmc.Settings.weight_window_checkpoints` attribute. (`#2670
<https://github.com/openmc-dev/openmc/pull/2670>`_)
- The :class:`~openmc.Settings` class now has a
:attr:`~openmc.Settings.max_write_lost_particles` attribute that can limit the
number of lost particle files written. (`#2688
<https://github.com/openmc-dev/openmc/pull/2688>`_)
- The :class:`~openmc.deplete.CoupledOperator` class now has a
``diff_volume_method`` argument that specifies how the volume of new materials
should be determined. (`#2691
<https://github.com/openmc-dev/openmc/pull/2691>`_)
- The :meth:`openmc.DAGMCUniverse.bounding_region` method now has a
``padding_distance`` argument. (`#2701
<https://github.com/openmc-dev/openmc/pull/2701>`_)
- A new :meth:`openmc.Material.get_decay_photon_energy` method replaces the
:attr:`decay_photon_energy` attribute and includes an ability to eliminate
low-importance points. This is facilitated by a new
:meth:`openmc.stats.Discrete.clip` method. (`#2715
<https://github.com/openmc-dev/openmc/pull/2715>`_)
- The :meth:`openmc.model.Model.differentiate_depletable_mats` method allows
depletable materials to be differentiated independent of the depletion
calculation itself. (`#2718
<https://github.com/openmc-dev/openmc/pull/2718>`_)
- Albedos can now be specified on surface boundary conditions. (`#2724
<https://github.com/openmc-dev/openmc/pull/2724>`_)
---------
Bug Fixes
---------
- Enable use of NCrystal materials in plot_xs (`#2435 <https://github.com/openmc-dev/openmc/pull/2435>`_)
- Avoid segfault from extern "C" std::string (`#2455 <https://github.com/openmc-dev/openmc/pull/2455>`_)
- Fix several issues with the Model class (`#2465 <https://github.com/openmc-dev/openmc/pull/2465>`_)
- Provide alternative batch estimation message (`#2479 <https://github.com/openmc-dev/openmc/pull/2479>`_)
- Correct index check for remove_tally (`#2494 <https://github.com/openmc-dev/openmc/pull/2494>`_)
- Support for NCrystal material in from_xml_element (`#2496 <https://github.com/openmc-dev/openmc/pull/2496>`_)
- Fix compilation with gcc 5 (`#2498 <https://github.com/openmc-dev/openmc/pull/2498>`_)
- Fixed in the Tally::add_filter method (`#2501 <https://github.com/openmc-dev/openmc/pull/2501>`_)
- Fix meaning of "masking" for plots (`#2510 <https://github.com/openmc-dev/openmc/pull/2510>`_)
- Fix description of statepoint.batches in Settings class (`#2514 <https://github.com/openmc-dev/openmc/pull/2514>`_)
- Reorder list initialization of Plot constructor (`#2519 <https://github.com/openmc-dev/openmc/pull/2519>`_)
- Added mkdir to cwd argument in Model.run (`#2523 <https://github.com/openmc-dev/openmc/pull/2523>`_)
- Fix export of spherical coordinates in SphericalMesh (`#2538 <https://github.com/openmc-dev/openmc/pull/2538>`_)
- Add virtual destructor on PlottableInterface (`#2541 <https://github.com/openmc-dev/openmc/pull/2541>`_)
- Ensure parent directory is created during depletion (`#2543 <https://github.com/openmc-dev/openmc/pull/2543>`_)
- Fix potential out-of-bounds access in TimeFilter (`#2532 <https://github.com/openmc-dev/openmc/pull/2532>`_)
- Remove use of sscanf for reading surface coefficients (`#2574 <https://github.com/openmc-dev/openmc/pull/2574>`_)
- Fix torus intersection bug (`#2589 <https://github.com/openmc-dev/openmc/pull/2589>`_)
- Multigroup per-thread cache fixes (`#2591 <https://github.com/openmc-dev/openmc/pull/2591>`_)
- Bank surface source particles in all active cycles (`#2592 <https://github.com/openmc-dev/openmc/pull/2592>`_)
- Fix for muir standard deviation (`#2598 <https://github.com/openmc-dev/openmc/pull/2598>`_)
- Check for zero fission cross section (`#2600 <https://github.com/openmc-dev/openmc/pull/2600>`_)
- XML read fixes in Plot classes (`#2623 <https://github.com/openmc-dev/openmc/pull/2623>`_)
- Added infinity check in VolumeCalculation (`#2634 <https://github.com/openmc-dev/openmc/pull/2634>`_)
- Fix sampling issue in Mixture distributions (`#2658 <https://github.com/openmc-dev/openmc/pull/2658>`_)
- Prevent segfault in distance to boundary calculation (`#2659 <https://github.com/openmc-dev/openmc/pull/2659>`_)
- Several CylindricalMesh fixes (`#2676
<https://github.com/openmc-dev/openmc/pull/2676>`_, `#2680
<https://github.com/openmc-dev/openmc/pull/2680>`_, `#2684
<https://github.com/openmc-dev/openmc/pull/2684>`_, `#2710
<https://github.com/openmc-dev/openmc/pull/2710>`_)
- Add type checks on Intersection, Union, Complement (`#2685 <https://github.com/openmc-dev/openmc/pull/2685>`_)
- Fixed typo in CF4Integrator docstring (`#2704 <https://github.com/openmc-dev/openmc/pull/2704>`_)
- Ensure property setters are used in CylindricalMesh and SphericalMesh (`#2709 <https://github.com/openmc-dev/openmc/pull/2709>`_)
- Fix sample_external_source bug (`#2713 <https://github.com/openmc-dev/openmc/pull/2713>`_)
- Fix localization issue affecting openmc-plotter (`#2723 <https://github.com/openmc-dev/openmc/pull/2723>`_)
- Correct openmc.lib wrapper for evaluate_legendre (`#2729 <https://github.com/openmc-dev/openmc/pull/2729>`_)
- Bug fix in Region.from_expression during tokenization (`#2733 <https://github.com/openmc-dev/openmc/pull/2733>`_)
- Fix bug in temperature interpolation (`#2734 <https://github.com/openmc-dev/openmc/pull/2734>`_)
- Check for invalid domain IDs in volume calculations (`#2742 <https://github.com/openmc-dev/openmc/pull/2742>`_)
- Skip boundary condition check for volume calculations (`#2743 <https://github.com/openmc-dev/openmc/pull/2743>`_)
- Fix loop over coordinates for source domain rejection (`#2751 <https://github.com/openmc-dev/openmc/pull/2751>`_)
------------
Contributors
------------
- `April Novak <https://github.com/aprilnovak>`_
- `Baptiste Mouginot <https://github.com/bam241>`_
- `Ben Collins <https://github.com/bscollin>`_
- `Chritopher Billingham <https://github.com/cadarache2014>`_
- `Christopher Fichtlscherer <https://github.com/cfichtlscherer>`_
- `Christina Cai <https://github.com/christinacai123>`_
- `Lorenzo Chierici <https://github.com/church89>`_
- `Huw Rhys Jones <https://github.com/dubway420>`_
- `Emilio Castro <https://github.com/ecasglez>`_
- `Erik Knudsen <https://github.com/ebknudsen>`_
- `Ethan Peterson <https://github.com/eepeterson>`_
- `Egor Afanasenko <https://github.com/egor1abs>`_
- `Paul Wilson <https://github.com/gonuke>`_
- `Gavin Ridley <https://github.com/gridley>`_
- `Hunter Belanger <https://github.com/HunterBelanger>`_
- `Jack Fletcher <https://github.com/j-fletcher>`_
- `John Vincent Cauilan <https://github.com/johvincau>`_
- `Josh May <https://github.com/joshmay1>`_
- `John Tramm <https://github.com/jtramm>`_
- `Kevin McLaughlin <https://github.com/kevinm387>`_
- `Yue Jin <https://github.com/kingyue737>`_
- `Lewis Gross <https://github.com/lewisgross1296>`_
- `Luke Labrie-Cleary <https://github.com/LukeLabie>`_
- `Patrick Myers <https://github.com/myerspat>`_
- `Nicola Rizzi <https://github.com/nicriz>`_
- `Yuvraj Jain <https://github.com/nutcasev15>`_
- `Paul Romano <https://github.com/paulromano>`_
- `Patrick Shriwise <https://github.com/pshriwise>`_
- `Rosie Barker <https://github.com/rlbarker>`_
- `Jonathan Shimwell <https://github.com/Shimwell>`_
- `John Tchakerian <https://github.com/stchaker>`_
- `Travis Labossiere-Hickman <https://github.com/tjlaboss>`_
- `Xinyan Wang <https://github.com/XinyanBradley>`_
- `Olek Yardas <https://github.com/yardasol>`_
- `Zoe Prieto <https://github.com/zoeprieto>`_

1
docs/source/releasenotes/index.rst
View file

@ -7,6 +7,7 @@ Release Notes
.. toctree::
:maxdepth: 1
0.14.0
0.13.3
0.13.2
0.13.1

45
docs/source/usersguide/geometry.rst
View file

@ -147,12 +147,13 @@ For many regions, a bounding-box can be determined automatically::
While a bounding box can be determined for regions involving half-spaces of
spheres, cylinders, and axis-aligned planes, it generally cannot be determined
if the region involves cones, non-axis-aligned planes, or other exotic
second-order surfaces. For example, the :func:`openmc.model.hexagonal_prism`
function returns the interior region of a hexagonal prism; because it is bounded
by a :class:`openmc.Plane`, trying to get its bounding box won't work::
second-order surfaces. For example, the :class:`openmc.model.HexagonalPrism`
class returns a hexagonal prism surface; because it utilizes a
:class:`openmc.Plane`, trying to get the bounding box of its interior won't
work::
>>> hex = openmc.model.hexagonal_prism()
>>> hex.bounding_box
>>> hex = openmc.model.HexagonalPrism()
>>> (-hex).bounding_box
(array([-0.8660254, -inf, -inf]),
array([ 0.8660254, inf, inf]))
@ -172,13 +173,17 @@ surface. To specify a vacuum boundary condition, simply change the
outer_surface = openmc.Sphere(r=100.0)
outer_surface.boundary_type = 'vacuum'
Reflective and periodic boundary conditions can be set with the strings
'reflective' and 'periodic'. Vacuum and reflective boundary conditions can be
applied to any type of surface. Periodic boundary conditions can be applied to
pairs of planar surfaces. If there are only two periodic surfaces they will be
matched automatically. Otherwise it is necessary to specify pairs explicitly
using the :attr:`Surface.periodic_surface` attribute as in the following
example::
Reflective, periodic, and white boundary conditions can be set with the
strings 'reflective', 'periodic', and 'white' respectively.
Vacuum, reflective and white boundary conditions can be applied to any
type of surface. The 'white' boundary condition supports diffuse particle
reflection in contrast to specular reflection provided by the 'reflective'
boundary condition.
Periodic boundary conditions can be applied to pairs of planar surfaces.
If there are only two periodic surfaces they will be matched automatically.
Otherwise it is necessary to specify pairs explicitly using the
:attr:`Surface.periodic_surface` attribute as in the following example::
p1 = openmc.Plane(a=0.3, b=5.0, d=1.0, boundary_type='periodic')
p2 = openmc.Plane(a=0.3, b=5.0, d=-1.0, boundary_type='periodic')
@ -196,6 +201,20 @@ lies in the first quadrant of the Cartesian grid. If the geometry instead lies
in the fourth quadrant, the :class:`YPlane` must be replaced by a
:class:`Plane` with the normal vector pointing in the :math:`-y` direction.
Additionally, 'reflective', 'periodic', and 'white' boundary conditions have
an albedo parameter that can be used to modify the importance of particles
that encounter the boundary. The albedo value specifies the ratio between
the particle's importance after interaction with the boundary to its initial
importance. The following example creates a reflective planar surface which
reduces the reflected particles' importance by 33.3%::
x1 = openmc.XPlane(1.0, boundary_type='reflective', albedo=0.667)
# This is equivalent
x1 = openmc.XPlane(1.0)
x1.boundary_type = 'reflective'
x1.albedo = 0.667
.. _usersguide_cells:
-----
@ -410,7 +429,7 @@ code would work::
hexlat.universes = [outer_ring, middle_ring, inner_ring]
If you need to create a hexagonal boundary (composed of six planar surfaces) for
a hexagonal lattice, :func:`openmc.model.hexagonal_prism` can be used.
a hexagonal lattice, :class:`openmc.model.HexagonalPrism` can be used.
.. _usersguide_geom_export:

104
docs/source/usersguide/plots.rst
View file

@ -126,27 +126,29 @@ will depend on the 3D viewer, but should be straightforward.
Projection Plots
----------------
.. image:: ../_images/hexlat_anim.gif
:width: 200px
.. only:: html
The :class:`openmc.ProjectionPlot` class presents an alternative method
of producing 3D visualizations of OpenMC geometries. It was developed to
overcome the primary shortcoming of voxel plots, that an enormous number
of voxels must be employed to capture detailed geometric features.
Projection plots perform volume rendering on material or
cell volumes, with colors specified in the same manner as slice plots.
This is done using the native ray tracing capabilities within OpenMC,
so any geometry in which particles successfully run without overlaps
or leaks will work with projection plots.
.. image:: ../_images/hexlat_anim.gif
:width: 200px
One drawback of projection plots is that particle tracks cannot be overlaid
on them at present. Moreover, checking for overlap regions is not currently possible with projection plots. The image heading this section can
be created by adding the following code to the hexagonal lattice example packaged
with OpenMC, before exporting to plots.xml.
The :class:`openmc.ProjectionPlot` class presents an alternative method of
producing 3D visualizations of OpenMC geometries. It was developed to overcome
the primary shortcoming of voxel plots, that an enormous number of voxels must
be employed to capture detailed geometric features. Projection plots perform
volume rendering on material or cell volumes, with colors specified in the same
manner as slice plots. This is done using the native ray tracing capabilities
within OpenMC, so any geometry in which particles successfully run without
overlaps or leaks will work with projection plots.
One drawback of projection plots is that particle tracks cannot be overlaid on
them at present. Moreover, checking for overlap regions is not currently
possible with projection plots. The image heading this section can be created by
adding the following code to the hexagonal lattice example packaged with OpenMC,
before exporting to plots.xml.
::
r = 5
r = 5
import numpy as np
for i in range(100):
phi = 2 * np.pi * i/100
@ -162,48 +164,46 @@ with OpenMC, before exporting to plots.xml.
thisp.xs[iron] = 1.0
thisp.wireframe_domains = [fuel]
thisp.wireframe_thickness = 2
plot_file.append(thisp)
This generates a sequence of png files which can be joined to form a gif.
Each image specifies a different camera position using some simple periodic
functions to create a perfectly looped gif. :attr:`ProjectionPlot.look_at`
defines where the camera's centerline should point at.
:attr:`ProjectionPlot.camera_position` similarly defines where the camera
is situated in the universe level we seek to plot. The other settings
resemble those employed by :class:`openmc.Plot`, with the exception of
the :class:`ProjectionPlot.set_transparent` method and :attr:`ProjectionPlot.xs`
dictionary. These are used to control volume rendering of material
volumes. "xs" here stands for cross section, and it defines material
opacities in units of inverse centimeters. Setting this value to a
large number would make a material or cell opaque, and setting it to
zero makes a material transparent. Thus, the :class:`ProjectionPlot.set_transparent`
can be used to make all materials in the geometry transparent. From there,
individual material or cell opacities can be tuned to produce the
desired result.
This generates a sequence of png files which can be joined to form a gif. Each
image specifies a different camera position using some simple periodic functions
to create a perfectly looped gif. :attr:`ProjectionPlot.look_at` defines where
the camera's centerline should point at. :attr:`ProjectionPlot.camera_position`
similarly defines where the camera is situated in the universe level we seek to
plot. The other settings resemble those employed by :class:`openmc.Plot`, with
the exception of the :class:`ProjectionPlot.set_transparent` method and
:attr:`ProjectionPlot.xs` dictionary. These are used to control volume rendering
of material volumes. "xs" here stands for cross section, and it defines material
opacities in units of inverse centimeters. Setting this value to a large number
would make a material or cell opaque, and setting it to zero makes a material
transparent. Thus, the :class:`ProjectionPlot.set_transparent` can be used to
make all materials in the geometry transparent. From there, individual material
or cell opacities can be tuned to produce the desired result.
Two camera projections are available when using these plots, perspective
and orthographic. The default, perspective projection,
is a cone of rays passing through each pixel which radiate from the camera
position and span the field of view in the x and y positions. The horizontal
field of view can be set with the :attr: `ProjectionPlot.horizontal_field_of_view` attribute,
which is to be specified in units of degrees. The field of view only influences
behavior in perspective projection mode.
Two camera projections are available when using these plots, perspective and
orthographic. The default, perspective projection, is a cone of rays passing
through each pixel which radiate from the camera position and span the field of
view in the x and y positions. The horizontal field of view can be set with the
:attr: `ProjectionPlot.horizontal_field_of_view` attribute, which is to be
specified in units of degrees. The field of view only influences behavior in
perspective projection mode.
In the orthographic projection, rays follow the same angle but originate from
different points. The horizontal width of this plane of ray starting points
may be set with the :attr: `ProjectionPlot.orthographic_width` element. If this element
is nonzero, the orthographic projection is employed. Left to its default value
of zero, the perspective projection is employed.
different points. The horizontal width of this plane of ray starting points may
be set with the :attr: `ProjectionPlot.orthographic_width` element. If this
element is nonzero, the orthographic projection is employed. Left to its default
value of zero, the perspective projection is employed.
Lastly, projection plots come packaged with wireframe generation that
can target either all surface/cell/material boundaries in the geometry,
or only wireframing around specific regions. In the above example, we
have set only the fuel region from the hexagonal lattice example to have
a wireframe drawn around it. This is accomplished by setting the
:attr: `ProjectionPlot.wireframe_domains`, which may be set to either material
IDs or cell IDs. The :attr:`ProjectionPlot.wireframe_thickness`
attribute sets the wireframe thickness in units of pixels.
Lastly, projection plots come packaged with wireframe generation that can target
either all surface/cell/material boundaries in the geometry, or only wireframing
around specific regions. In the above example, we have set only the fuel region
from the hexagonal lattice example to have a wireframe drawn around it. This is
accomplished by setting the :attr: `ProjectionPlot.wireframe_domains`, which may
be set to either material IDs or cell IDs. The
:attr:`ProjectionPlot.wireframe_thickness` attribute sets the wireframe
thickness in units of pixels.
.. note:: When setting specific material or cell regions to have wireframes
drawn around them, the plot must be colored by materials if wireframing

4
examples/assembly/assembly.py
View file

@ -99,11 +99,11 @@ def assembly_model():
assembly.universes[gt_pos[:, 0], gt_pos[:, 1]] = guide_tube_pin()
# Create outer boundary of the geometry to surround the lattice
outer_boundary = openmc.model.rectangular_prism(
outer_boundary = openmc.model.RectangularPrism(
pitch, pitch, boundary_type='reflective')
# Create a cell filled with the lattice
main_cell = openmc.Cell(fill=assembly, region=outer_boundary)
main_cell = openmc.Cell(fill=assembly, region=-outer_boundary)
# Finally, create geometry by providing a list of cells that fill the root
# universe

4
examples/custom_source/build_xml.py
View file

@ -8,8 +8,8 @@ mats = openmc.Materials([iron])
mats.export_to_xml()
# Create a 5 cm x 5 cm box filled with iron
box = openmc.model.rectangular_prism(10.0, 10.0, boundary_type='vacuum')
cell = openmc.Cell(fill=iron, region=box)
box = openmc.model.RectangularPrism(10.0, 10.0, boundary_type='vacuum')
cell = openmc.Cell(fill=iron, region=-box)
geometry = openmc.Geometry([cell])
geometry.export_to_xml()

4
examples/parameterized_custom_source/build_xml.py
View file

@ -8,8 +8,8 @@ mats = openmc.Materials([iron])
mats.export_to_xml()
# Create a 5 cm x 5 cm box filled with iron
box = openmc.model.rectangular_prism(10.0, 10.0, boundary_type='vacuum')
cell = openmc.Cell(fill=iron, region=box)
box = openmc.model.RectangularPrism(10.0, 10.0, boundary_type='vacuum')
cell = openmc.Cell(fill=iron, region=-box)
geometry = openmc.Geometry([cell])
geometry.export_to_xml()

4
examples/pincell/build_xml.py
View file

@ -43,13 +43,13 @@ clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
# Create a region represented as the inside of a rectangular prism
pitch = 1.25984
box = openmc.rectangular_prism(pitch, pitch, boundary_type='reflective')
box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')
# Create cells, mapping materials to regions
fuel = openmc.Cell(fill=uo2, region=-fuel_or)
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
water = openmc.Cell(fill=borated_water, region=+clad_or & box)
water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
# Create a geometry and export to XML
geometry = openmc.Geometry([fuel, gap, clad, water])

4
examples/pincell_depletion/run_depletion.py
View file

@ -41,13 +41,13 @@ pitch = 1.25984
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
box = openmc.model.rectangular_prism(pitch, pitch, boundary_type='reflective')
box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')
# Define cells
fuel = openmc.Cell(fill=uo2, region=-fuel_or)
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
water = openmc.Cell(fill=borated_water, region=+clad_or & box)
water = openmc.Cell(fill=borated_water, region=+clad_or & -box)
# Define overall geometry
geometry = openmc.Geometry([fuel, gap, clad, water])

4
examples/pincell_multigroup/build_xml.py
View file

@ -90,11 +90,11 @@ fuel_or = openmc.ZCylinder(r=0.54, name='Fuel OR')
# Create a region represented as the inside of a rectangular prism
pitch = 1.26
box = openmc.rectangular_prism(pitch, pitch, boundary_type='reflective')
box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')
# Instantiate Cells
fuel = openmc.Cell(fill=uo2, region=-fuel_or, name='fuel')
moderator = openmc.Cell(fill=water, region=+fuel_or & box, name='moderator')
moderator = openmc.Cell(fill=water, region=+fuel_or & -box, name='moderator')
# Create a geometry with the two cells and export to XML
geometry = openmc.Geometry([fuel, moderator])

41
include/openmc/boundary_condition.h
View file

@ -1,7 +1,10 @@
#ifndef OPENMC_BOUNDARY_CONDITION_H
#define OPENMC_BOUNDARY_CONDITION_H
#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"
#include "openmc/position.h"
#include <fmt/core.h>
namespace openmc {
@ -15,6 +18,8 @@ class Surface;
class BoundaryCondition {
public:
virtual ~BoundaryCondition() = default;
//! Perform tracking operations for a particle that strikes the boundary.
//! \param p The particle that struck the boundary. This class is not meant
//! to directly modify anything about the particle, but it will do so
@ -22,8 +27,44 @@ public:
//! \param surf The specific surface on the boundary the particle struck.
virtual void handle_particle(Particle& p, const Surface& surf) const = 0;
//! Modify the incident particle's weight according to the boundary's albedo.
//! \param p The particle that struck the boundary. This function calculates
//! the reduction in the incident particle's weight as it interacts
//! with a boundary. The lost weight is tallied before the remaining weight
//! is reassigned to the incident particle. Implementations of the
//! handle_particle function typically call this method in its body.
//! \param surf The specific surface on the boundary the particle struck.
void handle_albedo(Particle& p, const Surface& surf) const
{
if (!has_albedo())
return;
double initial_wgt = p.wgt();
// Treat the lost weight fraction as leakage, similar to VacuumBC.
// This ensures the lost weight is tallied properly.
p.wgt() *= (1.0 - albedo_);
p.cross_vacuum_bc(surf);
p.wgt() = initial_wgt * albedo_;
};
//! Return a string classification of this BC.
virtual std::string type() const = 0;
//! Write albedo data of this BC to hdf5.
void to_hdf5(hid_t surf_group) const
{
if (has_albedo()) {
write_string(surf_group, "albedo", fmt::format("{}", albedo_), false);
}
};
//! Set albedo of this BC.
void set_albedo(double albedo) { albedo_ = albedo; }
//! Return if this BC has an albedo.
bool has_albedo() const { return (albedo_ > 0.0); }
private:
double albedo_ = -1.0;
};
//==============================================================================

5
include/openmc/capi.h
View file

@ -93,6 +93,8 @@ int openmc_material_set_id(int32_t index, int32_t id);
int openmc_material_get_name(int32_t index, const char** name);
int openmc_material_set_name(int32_t index, const char* name);
int openmc_material_set_volume(int32_t index, double volume);
int openmc_material_get_depletable(int32_t index, bool* depletable);
int openmc_material_set_depletable(int32_t index, bool depletable);
int openmc_material_filter_get_bins(
int32_t index, const int32_t** bins, size_t* n);
int openmc_material_filter_set_bins(
@ -103,6 +105,9 @@ int openmc_mesh_filter_get_translation(int32_t index, double translation[3]);
int openmc_mesh_filter_set_translation(int32_t index, double translation[3]);
int openmc_mesh_get_id(int32_t index, int32_t* id);
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_material_volumes(int32_t index, int n_sample, int bin,
int result_size, void* result, int* hits, uint64_t* seed);
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);

14
include/openmc/cell.h
View file

@ -40,6 +40,7 @@ constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
//==============================================================================
class Cell;
class GeometryState;
class ParentCell;
class CellInstance;
class Universe;
@ -82,7 +83,7 @@ public:
//! Find the oncoming boundary of this cell.
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const;
Position r, Direction u, int32_t on_surface) const;
//! Get the BoundingBox for this cell.
BoundingBox bounding_box(int32_t cell_id) const;
@ -183,7 +184,7 @@ public:
//! Find the oncoming boundary of this cell.
virtual std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const = 0;
Position r, Direction u, int32_t on_surface, GeometryState* p) const = 0;
//! Write all information needed to reconstruct the cell to an HDF5 group.
//! \param group_id An HDF5 group id.
@ -260,7 +261,8 @@ protected:
//! \param[in] instance of the cell to find parent cells for
//! \param[in] p particle used to do a fast search for parent cells
//! \return parent cells
vector<ParentCell> find_parent_cells(int32_t instance, Particle& p) const;
vector<ParentCell> find_parent_cells(
int32_t instance, GeometryState& p) const;
//! Determine the path to this cell instance in the geometry hierarchy
//! \param[in] instance of the cell to find parent cells for
@ -332,10 +334,10 @@ public:
// Methods
vector<int32_t> surfaces() const override { return region_.surfaces(); }
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const override
std::pair<double, int32_t> distance(Position r, Direction u,
int32_t on_surface, GeometryState* p) const override
{
return region_.distance(r, u, on_surface, p);
return region_.distance(r, u, on_surface);
}
bool contains(Position r, Direction u, int32_t on_surface) const override

1
include/openmc/constants.h
View file

@ -317,7 +317,6 @@ enum class GlobalTally { K_COLLISION, K_ABSORPTION, K_TRACKLENGTH, LEAKAGE };
// Miscellaneous
constexpr int C_NONE {-1};
constexpr int F90_NONE {0}; // TODO: replace usage of this with C_NONE
// Interpolation rules
enum class Interpolation {

8
include/openmc/dagmc.h
View file

@ -42,7 +42,7 @@ public:
double evaluate(Position r) const override;
double distance(Position r, Direction u, bool coincident) const override;
Direction normal(Position r) const override;
Direction reflect(Position r, Direction u, Particle* p) const override;
Direction reflect(Position r, Direction u, GeometryState* p) const override;
inline void to_hdf5_inner(hid_t group_id) const override {};
@ -63,8 +63,8 @@ public:
bool contains(Position r, Direction u, int32_t on_surface) const override;
std::pair<double, int32_t> distance(
Position r, Direction u, int32_t on_surface, Particle* p) const override;
std::pair<double, int32_t> distance(Position r, Direction u,
int32_t on_surface, GeometryState* p) const override;
BoundingBox bounding_box() const override;
@ -143,7 +143,7 @@ public:
//! string of the ID ranges for entities of dimension \p dim
std::string dagmc_ids_for_dim(int dim) const;
bool find_cell(Particle& p) const override;
bool find_cell(GeometryState& p) const override;
void to_hdf5(hid_t universes_group) const override;

7
include/openmc/distribution.h
View file

@ -7,6 +7,7 @@
#include <cstddef> // for size_t
#include "pugixml.hpp"
#include <gsl/gsl-lite.hpp>
#include "openmc/constants.h"
#include "openmc/memory.h" // for unique_ptr
@ -43,9 +44,9 @@ class DiscreteIndex {
public:
DiscreteIndex() {};
DiscreteIndex(pugi::xml_node node);
DiscreteIndex(const double* p, int n);
DiscreteIndex(gsl::span<const double> p);
void assign(const double* p, int n);
void assign(gsl::span<const double> p);
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
@ -77,7 +78,7 @@ private:
class Discrete : public Distribution {
public:
explicit Discrete(pugi::xml_node node);
Discrete(const double* x, const double* p, int n);
Discrete(const double* x, const double* p, size_t n);
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer

2
include/openmc/distribution_multi.h
View file

@ -22,6 +22,8 @@ public:
explicit UnitSphereDistribution(pugi::xml_node node);
virtual ~UnitSphereDistribution() = default;
static unique_ptr<UnitSphereDistribution> create(pugi::xml_node node);
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Direction sampled

15
include/openmc/distribution_spatial.h
View file

@ -19,6 +19,8 @@ public:
//! Sample a position from the distribution
virtual Position sample(uint64_t* seed) const = 0;
static unique_ptr<SpatialDistribution> create(pugi::xml_node node);
};
//==============================================================================
@ -102,16 +104,27 @@ private:
class MeshSpatial : public SpatialDistribution {
public:
explicit MeshSpatial(pugi::xml_node node);
explicit MeshSpatial(int32_t mesh_id, gsl::span<const double> strengths);
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
const Mesh* mesh() const { return model::meshes.at(mesh_idx_).get(); }
//! Sample the mesh for an element and position within that element
//! \param seed Pseudorandom number seed pointer
//! \return Sampled element index and position within that element
std::pair<int32_t, Position> sample_mesh(uint64_t* seed) const;
//! For unstructured meshes, ensure that elements are all linear tetrahedra
void check_element_types() const;
// Accessors
const Mesh* mesh() const { return model::meshes.at(mesh_idx_).get(); }
int32_t n_sources() const { return this->mesh()->n_bins(); }
double total_strength() { return this->elem_idx_dist_.integral(); }
private:
int32_t mesh_idx_ {C_NONE};
DiscreteIndex elem_idx_dist_; //!< Distribution of

5
include/openmc/file_utils.h
View file

@ -16,6 +16,11 @@ bool dir_exists(const std::string& path);
//! \return Whether file exists
bool file_exists(const std::string& filename);
//! Determine directory containing given file
//! \param[in] filename Path to file
//! \return Name of directory containing file
std::string dir_name(const std::string& filename);
// Gets the file extension of whatever string is passed in. This is defined as
// a sequence of strictly alphanumeric characters which follow the last period,
// i.e. at least one alphabet character is present, and zero or more numbers.

16
include/openmc/geometry.h
View file

@ -11,7 +11,7 @@
namespace openmc {
class BoundaryInfo;
class Particle;
class GeometryState;
//==============================================================================
// Global variables
@ -39,7 +39,7 @@ inline bool coincident(double d1, double d2)
//! Check for overlapping cells at a particle's position.
//==============================================================================
bool check_cell_overlap(Particle& p, bool error = true);
bool check_cell_overlap(GeometryState& p, bool error = true);
//==============================================================================
//! Get the cell instance for a particle at the specified universe level
@ -50,7 +50,7 @@ bool check_cell_overlap(Particle& p, bool error = true);
//! should be computed. \return The instance of the cell at the specified level.
//==============================================================================
int cell_instance_at_level(const Particle& p, int level);
int cell_instance_at_level(const GeometryState& p, int level);
//==============================================================================
//! Locate a particle in the geometry tree and set its geometry data fields.
@ -60,20 +60,22 @@ int cell_instance_at_level(const Particle& p, int level);
//! \return True if the particle's location could be found and ascribed to a
//! valid geometry coordinate stack.
//==============================================================================
bool exhaustive_find_cell(Particle& p);
bool neighbor_list_find_cell(Particle& p); // Only usable on surface crossings
bool exhaustive_find_cell(GeometryState& p, bool verbose = false);
bool neighbor_list_find_cell(
GeometryState& p, bool verbose = false); // Only usable on surface crossings
//==============================================================================
//! Move a particle into a new lattice tile.
//==============================================================================
void cross_lattice(Particle& p, const BoundaryInfo& boundary);
void cross_lattice(
GeometryState& p, const BoundaryInfo& boundary, bool verbose = false);
//==============================================================================
//! Find the next boundary a particle will intersect.
//==============================================================================
BoundaryInfo distance_to_boundary(Particle& p);
BoundaryInfo distance_to_boundary(GeometryState& p);
} // namespace openmc

12
include/openmc/interpolate.h
View file

@ -4,6 +4,9 @@
#include <cmath>
#include <vector>
#include <gsl/gsl-lite.hpp>
#include "openmc/error.h"
#include "openmc/search.h"
namespace openmc {
@ -47,14 +50,15 @@ inline double interpolate_lagrangian(gsl::span<const double> xs,
numerator *= (x - xs[idx + j]);
denominator *= (xs[idx + i] - xs[idx + j]);
}
output += (numerator / denominator) * ys[i];
output += (numerator / denominator) * ys[idx + i];
}
return output;
}
double interpolate(gsl::span<const double> xs, gsl::span<const double> ys,
double x, Interpolation i = Interpolation::lin_lin)
inline double interpolate(gsl::span<const double> xs,
gsl::span<const double> ys, double x,
Interpolation i = Interpolation::lin_lin)
{
int idx = lower_bound_index(xs.begin(), xs.end(), x);
@ -92,4 +96,4 @@ double interpolate(gsl::span<const double> xs, gsl::span<const double> ys,
} // namespace openmc
#endif
#endif

17
include/openmc/material.h
View file

@ -94,7 +94,7 @@ public:
// which will get auto-assigned to the next available ID. After creating
// the new material, it is added to openmc::model::materials.
//! \return reference to the cloned material
Material & clone();
Material& clone();
//----------------------------------------------------------------------------
// Accessors
@ -149,6 +149,7 @@ public:
//! Get whether material is fissionable
//! \return Whether material is fissionable
bool fissionable() const { return fissionable_; }
bool& fissionable() { return fissionable_; }
//! Get volume of material
//! \return Volume in [cm^3]
@ -158,6 +159,10 @@ public:
//! \return Temperature in [K]
double temperature() const;
//! Whether or not the material is depletable
bool depletable() const { return depletable_; }
bool& depletable() { return depletable_; }
//! Get pointer to NCrystal material object
//! \return Pointer to NCrystal material object
const NCrystalMat& ncrystal_mat() const { return ncrystal_mat_; };
@ -173,11 +178,8 @@ public:
double density_; //!< Total atom density in [atom/b-cm]
double density_gpcc_; //!< Total atom density in [g/cm^3]
double volume_ {-1.0}; //!< Volume in [cm^3]
bool fissionable_ {
false}; //!< Does this material contain fissionable nuclides
bool depletable_ {false}; //!< Is the material depletable?
vector<bool> p0_; //!< Indicate which nuclides are to be treated with
//!< iso-in-lab scattering
vector<bool> p0_; //!< Indicate which nuclides are to be treated with
//!< iso-in-lab scattering
// To improve performance of tallying, we store an array (direct address
// table) that indicates for each nuclide in data::nuclides the index of the
@ -210,6 +212,9 @@ private:
// Private data members
gsl::index index_;
bool depletable_ {false}; //!< Is the material depletable?
bool fissionable_ {
false}; //!< Does this material contain fissionable nuclides
//! \brief Default temperature for cells containing this material.
//!
//! A negative value indicates no default temperature was specified.

110
include/openmc/mesh.h
View file

@ -9,7 +9,9 @@
#include "hdf5.h"
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
#include <gsl/gsl-lite.hpp>
#include "openmc/error.h"
#include "openmc/memory.h" // for unique_ptr
#include "openmc/particle.h"
#include "openmc/position.h"
@ -68,12 +70,20 @@ extern const libMesh::Parallel::Communicator* libmesh_comm;
class Mesh {
public:
// Types, aliases
struct MaterialVolume {
int32_t material; //!< material index
double volume; //!< volume in [cm^3]
};
// Constructors and destructor
Mesh() = default;
Mesh(pugi::xml_node node);
virtual ~Mesh() = default;
// Methods
//! Perform any preparation needed to support use in mesh filters
virtual void prepare_for_tallies() {};
//! Update a position to the local coordinates of the mesh
virtual void local_coords(Position& r) const {};
@ -81,12 +91,12 @@ public:
//! Return a position in the local coordinates of the mesh
virtual Position local_coords(const Position& r) const { return r; };
//! Sample a mesh volume using a certain seed
//! Sample a position within a mesh element
//
//! \param[in] seed Seed to use for random sampling
//! \param[in] bin Bin value of the tet sampled
//! \return sampled position within tet
virtual Position sample(uint64_t* seed, int32_t bin) const = 0;
//! \param[in] bin Bin value of the mesh element sampled
//! \param[inout] seed Seed to use for random sampling
//! \return sampled position within mesh element
virtual Position sample_element(int32_t bin, uint64_t* seed) const = 0;
//! Determine which bins were crossed by a particle
//
@ -156,9 +166,28 @@ public:
virtual std::string get_mesh_type() const = 0;
//! Determine volume of materials within a single mesh elemenet
//
//! \param[in] n_sample Number of samples within each element
//! \param[in] bin Index of mesh element
//! \param[out] Array of (material index, volume) for desired element
//! \param[inout] seed Pseudorandom number seed
//! \return Number of materials within element
int material_volumes(int n_sample, int bin, gsl::span<MaterialVolume> volumes,
uint64_t* seed) const;
//! Determine volume of materials within a single mesh elemenet
//
//! \param[in] n_sample Number of samples within each element
//! \param[in] bin Index of mesh element
//! \param[inout] seed Pseudorandom number seed
//! \return Vector of (material index, volume) for desired element
vector<MaterialVolume> material_volumes(
int n_sample, int bin, uint64_t* seed) const;
// Data members
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
int id_ {-1}; //!< User-specified ID
int n_dimension_ {-1}; //!< Number of dimensions
};
class StructuredMesh : public Mesh {
@ -183,7 +212,12 @@ public:
}
};
Position sample(uint64_t* seed, int32_t bin) const override;
Position sample_element(int32_t bin, uint64_t* seed) const override
{
return sample_element(get_indices_from_bin(bin), seed);
};
virtual Position sample_element(const MeshIndex& ijk, uint64_t* seed) const;
int get_bin(Position r) const override;
@ -240,6 +274,30 @@ public:
//! \param[in] i Direction index
virtual int get_index_in_direction(double r, int i) const = 0;
//! Get the coordinate for the mesh grid boundary in the positive direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
virtual double positive_grid_boundary(const MeshIndex& ijk, int i) const
{
auto msg =
fmt::format("Attempting to call positive_grid_boundary on a {} mesh.",
get_mesh_type());
fatal_error(msg);
};
//! Get the coordinate for the mesh grid boundary in the negative direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
virtual double negative_grid_boundary(const MeshIndex& ijk, int i) const
{
auto msg =
fmt::format("Attempting to call negative_grid_boundary on a {} mesh.",
get_mesh_type());
fatal_error(msg);
};
//! Get the closest distance from the coordinate r to the grid surface
//! in i direction that bounds mesh cell ijk and that is larger than l
//! The coordinate r does not have to be inside the mesh cell ijk. In
@ -322,18 +380,17 @@ public:
void to_hdf5(hid_t group) const override;
// New methods
//! Get the coordinate for the mesh grid boundary in the positive direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double positive_grid_boundary(const MeshIndex& ijk, int i) const;
double positive_grid_boundary(const MeshIndex& ijk, int i) const override;
//! Get the coordinate for the mesh grid boundary in the negative direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double negative_grid_boundary(const MeshIndex& ijk, int i) const;
double negative_grid_boundary(const MeshIndex& ijk, int i) const override;
//! Count number of bank sites in each mesh bin / energy bin
//
@ -373,18 +430,17 @@ public:
void to_hdf5(hid_t group) const override;
// New methods
//! Get the coordinate for the mesh grid boundary in the positive direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double positive_grid_boundary(const MeshIndex& ijk, int i) const;
double positive_grid_boundary(const MeshIndex& ijk, int i) const override;
//! Get the coordinate for the mesh grid boundary in the negative direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double negative_grid_boundary(const MeshIndex& ijk, int i) const;
double negative_grid_boundary(const MeshIndex& ijk, int i) const override;
//! Return the volume for a given mesh index
double volume(const MeshIndex& ijk) const override;
@ -410,6 +466,8 @@ public:
static const std::string mesh_type;
Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const override;
@ -420,10 +478,16 @@ public:
double volume(const MeshIndex& ijk) const override;
array<vector<double>, 3> grid_;
// grid accessors
double r(int i) const { return grid_[0][i]; }
double phi(int i) const { return grid_[1][i]; }
double z(int i) const { return grid_[2][i]; }
int set_grid();
// Data members
array<vector<double>, 3> grid_;
private:
double find_r_crossing(
const Position& r, const Direction& u, double l, int shell) const;
@ -466,6 +530,8 @@ public:
static const std::string mesh_type;
Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const override;
@ -474,10 +540,15 @@ public:
void to_hdf5(hid_t group) const override;
array<vector<double>, 3> grid_;
double r(int i) const { return grid_[0][i]; }
double theta(int i) const { return grid_[1][i]; }
double phi(int i) const { return grid_[2][i]; }
int set_grid();
// Data members
array<vector<double>, 3> grid_;
private:
double find_r_crossing(
const Position& r, const Direction& u, double l, int shell) const;
@ -621,7 +692,10 @@ public:
// Overridden Methods
Position sample(uint64_t* seed, int32_t bin) const override;
//! Perform any preparation needed to support use in mesh filters
void prepare_for_tallies() override;
Position sample_element(int32_t bin, uint64_t* seed) const override;
void bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins, vector<double>& lengths) const override;
@ -789,7 +863,7 @@ public:
void bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins, vector<double>& lengths) const override;
Position sample(uint64_t* seed, int32_t bin) const override;
Position sample_element(int32_t bin, uint64_t* seed) const override;
int get_bin(Position r) const override;

46
include/openmc/openmp_interface.h
View file

@ -7,6 +7,27 @@
namespace openmc {
//==============================================================================
//! Accessor functions related to number of threads and thread number
//==============================================================================
inline int num_threads()
{
#ifdef _OPENMP
return omp_get_max_threads();
#else
return 1;
#endif
}
inline int thread_num()
{
#ifdef _OPENMP
return omp_get_thread_num();
#else
return 0;
#endif
}
//==============================================================================
//! An object used to prevent concurrent access to a piece of data.
//
@ -29,11 +50,26 @@ public:
#endif
}
// Mutexes cannot be copied. We need to explicitly delete the copy
// constructor and copy assignment operator to ensure the compiler doesn't
// "help" us by implicitly trying to copy the underlying mutexes.
OpenMPMutex(const OpenMPMutex&) = delete;
OpenMPMutex& operator=(const OpenMPMutex&) = delete;
// omp_lock_t objects cannot be deep copied, they can only be shallow
// copied. Thus, while shallow copying of an omp_lock_t object is
// completely valid (provided no race conditions exist), true copying
// of an OpenMPMutex object is not valid due to the action of the
// destructor. However, since locks are fungible, we can simply replace
// copying operations with default construction. This allows storage of
// OpenMPMutex objects within containers that may need to move/copy them
// (e.g., std::vector). It is left to the caller to understand that
// copying of OpenMPMutex does not produce two handles to the same mutex,
// rather, it produces two different mutexes.
// Copy constructor
OpenMPMutex(const OpenMPMutex& other) { OpenMPMutex(); }
// Copy assignment operator
OpenMPMutex& operator=(const OpenMPMutex& other)
{
OpenMPMutex();
return *this;
}
//! Lock the mutex.
//

14
include/openmc/particle.h
View file

@ -5,7 +5,6 @@
//! \brief Particle type
#include <cstdint>
#include <sstream>
#include <string>
#include "openmc/constants.h"
@ -103,17 +102,8 @@ public:
//! mark a particle as lost and create a particle restart file
//! \param message A warning message to display
void mark_as_lost(const char* message);
void mark_as_lost(const std::string& message)
{
mark_as_lost(message.c_str());
}
void mark_as_lost(const std::stringstream& message)
{
mark_as_lost(message.str());
}
virtual void mark_as_lost(const char* message) override;
using GeometryState::mark_as_lost;
//! create a particle restart HDF5 file
void write_restart() const;

439
include/openmc/particle_data.h
View file

@ -194,6 +194,157 @@ struct BoundaryInfo {
lattice_translation {}; //!< which way lattice indices will change
};
/*
* Contains all geometry state information for a particle.
*/
class GeometryState {
public:
GeometryState();
/*
* GeometryState does not store any ID info, so give some reasonable behavior
* here. The Particle class redefines this. This is only here for the error
* reporting behavior that occurs in geometry.cpp. The explanation for
* mark_as_lost is the same.
*/
virtual void mark_as_lost(const char* message);
void mark_as_lost(const std::string& message);
void mark_as_lost(const std::stringstream& message);
// resets all coordinate levels for the particle
void clear()
{
for (auto& level : coord_)
level.reset();
n_coord_ = 1;
}
// Initialize all internal state from position and direction
void init_from_r_u(Position r_a, Direction u_a)
{
clear();
surface() = 0;
material() = C_NONE;
r() = r_a;
u() = u_a;
r_last_current() = r_a;
r_last() = r_a;
u_last() = u_a;
}
// Unique ID. This is not geometric info, but the
// error reporting in geometry.cpp requires this.
// We could save this to implement it in Particle,
// but that would require virtuals.
int64_t& id() { return id_; }
const int64_t& id() const { return id_; }
// Number of current coordinate levels
int& n_coord() { return n_coord_; }
const int& n_coord() const { return n_coord_; }
// Offset for distributed properties
int& cell_instance() { return cell_instance_; }
const int& cell_instance() const { return cell_instance_; }
// Coordinates for all nesting levels
LocalCoord& coord(int i) { return coord_[i]; }
const LocalCoord& coord(int i) const { return coord_[i]; }
const vector<LocalCoord>& coord() const { return coord_; }
// Innermost universe nesting coordinates
LocalCoord& lowest_coord() { return coord_[n_coord_ - 1]; }
const LocalCoord& lowest_coord() const { return coord_[n_coord_ - 1]; }
// Last coordinates on all nesting levels, before crossing a surface
int& n_coord_last() { return n_coord_last_; }
const int& n_coord_last() const { return n_coord_last_; }
int& cell_last(int i) { return cell_last_[i]; }
const int& cell_last(int i) const { return cell_last_[i]; }
// Coordinates of last collision or reflective/periodic surface
// crossing for current tallies
Position& r_last_current() { return r_last_current_; }
const Position& r_last_current() const { return r_last_current_; }
// Previous direction and spatial coordinates before a collision
Position& r_last() { return r_last_; }
const Position& r_last() const { return r_last_; }
Position& u_last() { return u_last_; }
const Position& u_last() const { return u_last_; }
// Accessors for position in global coordinates
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
// Accessors for position in local coordinates
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
// Accessors for direction in global coordinates
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
// Accessors for direction in local coordinates
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
// Surface that the particle is on
int& surface() { return surface_; }
const int& surface() const { return surface_; }
// Boundary information
BoundaryInfo& boundary() { return boundary_; }
#ifdef DAGMC
// DagMC state variables
moab::DagMC::RayHistory& history() { return history_; }
Direction& last_dir() { return last_dir_; }
#endif
// material of current and last cell
int& material() { return material_; }
const int& material() const { return material_; }
int& material_last() { return material_last_; }
const int& material_last() const { return material_last_; }
// temperature of current and last cell
double& sqrtkT() { return sqrtkT_; }
const double& sqrtkT() const { return sqrtkT_; }
double& sqrtkT_last() { return sqrtkT_last_; }
private:
int64_t id_ {-1}; //!< Unique ID
int n_coord_ {1}; //!< number of current coordinate levels
int cell_instance_; //!< offset for distributed properties
vector<LocalCoord> coord_; //!< coordinates for all levels
int n_coord_last_ {1}; //!< number of current coordinates
vector<int> cell_last_; //!< coordinates for all levels
Position r_last_current_; //!< coordinates of the last collision or
//!< reflective/periodic surface crossing for
//!< current tallies
Position r_last_; //!< previous coordinates
Direction u_last_; //!< previous direction coordinates
int surface_ {0}; //!< index for surface particle is on
BoundaryInfo boundary_; //!< Info about the next intersection
int material_ {-1}; //!< index for current material
int material_last_ {-1}; //!< index for last material
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
#ifdef DAGMC
moab::DagMC::RayHistory history_;
Direction last_dir_;
#endif
};
//============================================================================
//! Defines how particle data is laid out in memory
//============================================================================
@ -229,163 +380,112 @@ struct BoundaryInfo {
* Algorithms.” Annals of Nuclear Energy 113 (March 2018): 506–18.
* https://doi.org/10.1016/j.anucene.2017.11.032.
*/
class ParticleData {
public:
//----------------------------------------------------------------------------
// Constructors
ParticleData();
class ParticleData : public GeometryState {
private:
//==========================================================================
// Data members (accessor methods are below)
// Data members -- see public: below for descriptions
// Cross section caches
vector<NuclideMicroXS> neutron_xs_; //!< Microscopic neutron cross sections
vector<ElementMicroXS> photon_xs_; //!< Microscopic photon cross sections
MacroXS macro_xs_; //!< Macroscopic cross sections
CacheDataMG mg_xs_cache_; //!< Multigroup XS cache
vector<NuclideMicroXS> neutron_xs_;
vector<ElementMicroXS> photon_xs_;
MacroXS macro_xs_;
CacheDataMG mg_xs_cache_;
int64_t id_; //!< Unique ID
ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.)
ParticleType type_ {ParticleType::neutron};
int n_coord_ {1}; //!< number of current coordinate levels
int cell_instance_; //!< offset for distributed properties
vector<LocalCoord> coord_; //!< coordinates for all levels
double E_;
double E_last_;
int g_ {0};
int g_last_;
// Particle coordinates before crossing a surface
int n_coord_last_ {1}; //!< number of current coordinates
vector<int> cell_last_; //!< coordinates for all levels
double wgt_ {1.0};
double mu_;
double time_ {0.0};
double time_last_ {0.0};
double wgt_last_ {1.0};
// Energy data
double E_; //!< post-collision energy in eV
double E_last_; //!< pre-collision energy in eV
int g_ {0}; //!< post-collision energy group (MG only)
int g_last_; //!< pre-collision energy group (MG only)
bool fission_ {false};
TallyEvent event_;
int event_nuclide_;
int event_mt_;
int delayed_group_ {0};
// Other physical data
double wgt_ {1.0}; //!< particle weight
double mu_; //!< angle of scatter
double time_ {0.0}; //!< time in [s]
double time_last_ {0.0}; //!< previous time in [s]
int n_bank_ {0};
int n_bank_second_ {0};
double wgt_bank_ {0.0};
int n_delayed_bank_[MAX_DELAYED_GROUPS];
// Other physical data
Position r_last_current_; //!< coordinates of the last collision or
//!< reflective/periodic surface crossing for
//!< current tallies
Position r_last_; //!< previous coordinates
Direction u_last_; //!< previous direction coordinates
double wgt_last_ {1.0}; //!< pre-collision particle weight
int cell_born_ {-1};
// What event took place
bool fission_ {false}; //!< did particle cause implicit fission
TallyEvent event_; //!< scatter, absorption
int event_nuclide_; //!< index in nuclides array
int event_mt_; //!< reaction MT
int delayed_group_ {0}; //!< delayed group
int n_collision_ {0};
// Post-collision physical data
int n_bank_ {0}; //!< number of fission sites banked
int n_bank_second_ {0}; //!< number of secondary particles banked
double wgt_bank_ {0.0}; //!< weight of fission sites banked
int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission
//!< sites banked
// Indices for various arrays
int surface_ {0}; //!< index for surface particle is on
int cell_born_ {-1}; //!< index for cell particle was born in
int material_ {-1}; //!< index for current material
int material_last_ {-1}; //!< index for last material
// Boundary information
BoundaryInfo boundary_;
// Temperature of current cell
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double sqrtkT_last_ {0.0}; //!< last temperature
// Statistical data
int n_collision_ {0}; //!< number of collisions
// Track output
bool write_track_ {false};
// Current PRNG state
uint64_t seeds_[N_STREAMS]; // current seeds
int stream_; // current RNG stream
uint64_t seeds_[N_STREAMS];
int stream_;
// Secondary particle bank
vector<SourceSite> secondary_bank_;
int64_t current_work_; // current work index
int64_t current_work_;
vector<double> flux_derivs_; // for derivatives for this particle
vector<double> flux_derivs_;
vector<FilterMatch> filter_matches_; // tally filter matches
vector<FilterMatch> filter_matches_;
vector<TrackStateHistory> tracks_; // tracks for outputting to file
vector<TrackStateHistory> tracks_;
vector<NuBank> nu_bank_; // bank of most recently fissioned particles
vector<NuBank> nu_bank_;
vector<double> pht_storage_; // interim pulse-height results
vector<double> pht_storage_;
// Global tally accumulators
double keff_tally_absorption_ {0.0};
double keff_tally_collision_ {0.0};
double keff_tally_tracklength_ {0.0};
double keff_tally_leakage_ {0.0};
bool trace_ {false}; //!< flag to show debug information
bool trace_ {false};
double collision_distance_; // distance to particle's next closest collision
double collision_distance_;
int n_event_ {0}; // number of events executed in this particle's history
int n_event_ {0};
// Weight window information
int n_split_ {0}; // Number of times this particle has been split
double ww_factor_ {
0.0}; // Particle-specific factor for on-the-fly weight window adjustment
int n_split_ {0};
double ww_factor_ {0.0};
// DagMC state variables
#ifdef DAGMC
moab::DagMC::RayHistory history_;
Direction last_dir_;
#endif
int64_t n_progeny_ {0}; // Number of progeny produced by this particle
int64_t n_progeny_ {0};
public:
//----------------------------------------------------------------------------
// Constructors
ParticleData();
//==========================================================================
// Methods and accessors
NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; }
// Cross section caches
NuclideMicroXS& neutron_xs(int i)
{
return neutron_xs_[i];
} // Microscopic neutron cross sections
const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; }
// Microscopic photon cross sections
ElementMicroXS& photon_xs(int i) { return photon_xs_[i]; }
// Macroscopic cross sections
MacroXS& macro_xs() { return macro_xs_; }
const MacroXS& macro_xs() const { return macro_xs_; }
// Multigroup macroscopic cross sections
CacheDataMG& mg_xs_cache() { return mg_xs_cache_; }
const CacheDataMG& mg_xs_cache() const { return mg_xs_cache_; }
int64_t& id() { return id_; }
const int64_t& id() const { return id_; }
// Particle type (n, p, e, gamma, etc)
ParticleType& type() { return type_; }
const ParticleType& type() const { return type_; }
int& n_coord() { return n_coord_; }
const int& n_coord() const { return n_coord_; }
int& cell_instance() { return cell_instance_; }
const int& cell_instance() const { return cell_instance_; }
LocalCoord& coord(int i) { return coord_[i]; }
const LocalCoord& coord(int i) const { return coord_[i]; }
const vector<LocalCoord>& coord() const { return coord_; }
LocalCoord& lowest_coord() { return coord_[n_coord_ - 1]; }
const LocalCoord& lowest_coord() const { return coord_[n_coord_ - 1]; }
int& n_coord_last() { return n_coord_last_; }
const int& n_coord_last() const { return n_coord_last_; }
int& cell_last(int i) { return cell_last_[i]; }
const int& cell_last(int i) const { return cell_last_[i]; }
// Current particle energy, energy before collision,
// and corresponding multigroup group indices. Energy
// units are eV.
double& E() { return E_; }
const double& E() const { return E_; }
double& E_last() { return E_last_; }
@ -395,112 +495,119 @@ public:
int& g_last() { return g_last_; }
const int& g_last() const { return g_last_; }
// Statistic weight of particle. Setting to zero
// indicates that the particle is dead.
double& wgt() { return wgt_; }
double wgt() const { return wgt_; }
double& wgt_last() { return wgt_last_; }
const double& wgt_last() const { return wgt_last_; }
bool alive() const { return wgt_ != 0.0; }
// Polar scattering angle after a collision
double& mu() { return mu_; }
const double& mu() const { return mu_; }
// Tracks the time of a particle as it traverses the problem.
// Units are seconds.
double& time() { return time_; }
const double& time() const { return time_; }
double& time_last() { return time_last_; }
const double& time_last() const { return time_last_; }
bool alive() const { return wgt_ != 0.0; }
Position& r_last_current() { return r_last_current_; }
const Position& r_last_current() const { return r_last_current_; }
Position& r_last() { return r_last_; }
const Position& r_last() const { return r_last_; }
Position& u_last() { return u_last_; }
const Position& u_last() const { return u_last_; }
double& wgt_last() { return wgt_last_; }
const double& wgt_last() const { return wgt_last_; }
bool& fission() { return fission_; }
// What event took place, described in greater detail below
TallyEvent& event() { return event_; }
const TallyEvent& event() const { return event_; }
int& event_nuclide() { return event_nuclide_; }
bool& fission() { return fission_; } // true if implicit fission
int& event_nuclide() { return event_nuclide_; } // index of collision nuclide
const int& event_nuclide() const { return event_nuclide_; }
int& event_mt() { return event_mt_; }
int& delayed_group() { return delayed_group_; }
int& event_mt() { return event_mt_; } // MT number of collision
int& delayed_group() { return delayed_group_; } // delayed group
int& n_bank() { return n_bank_; }
int& n_bank_second() { return n_bank_second_; }
double& wgt_bank() { return wgt_bank_; }
int* n_delayed_bank() { return n_delayed_bank_; }
int& n_delayed_bank(int i) { return n_delayed_bank_[i]; }
// Post-collision data
int& n_bank() { return n_bank_; } // number of banked fission sites
int& n_bank_second()
{
return n_bank_second_;
} // number of secondaries banked
double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites
int* n_delayed_bank()
{
return n_delayed_bank_;
} // number of delayed fission sites
int& n_delayed_bank(int i)
{
return n_delayed_bank_[i];
} // number of delayed fission sites
int& surface() { return surface_; }
const int& surface() const { return surface_; }
// Index of cell particle is born in
int& cell_born() { return cell_born_; }
const int& cell_born() const { return cell_born_; }
int& material() { return material_; }
const int& material() const { return material_; }
int& material_last() { return material_last_; }
BoundaryInfo& boundary() { return boundary_; }
double& sqrtkT() { return sqrtkT_; }
const double& sqrtkT() const { return sqrtkT_; }
double& sqrtkT_last() { return sqrtkT_last_; }
// index of the current and last material
// Total number of collisions suffered by particle
int& n_collision() { return n_collision_; }
const int& n_collision() const { return n_collision_; }
// whether this track is to be written
bool& write_track() { return write_track_; }
// RNG state
uint64_t& seeds(int i) { return seeds_[i]; }
uint64_t* seeds() { return seeds_; }
int& stream() { return stream_; }
// secondary particle bank
SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
// Current simulation work index
int64_t& current_work() { return current_work_; }
const int64_t& current_work() const { return current_work_; }
// Used in tally derivatives
double& flux_derivs(int i) { return flux_derivs_[i]; }
const double& flux_derivs(int i) const { return flux_derivs_[i]; }
// Matches of tallies
decltype(filter_matches_)& filter_matches() { return filter_matches_; }
FilterMatch& filter_matches(int i) { return filter_matches_[i]; }
// Tracks to output to file
decltype(tracks_)& tracks() { return tracks_; }
// Bank of recently fissioned particles
decltype(nu_bank_)& nu_bank() { return nu_bank_; }
NuBank& nu_bank(int i) { return nu_bank_[i]; }
// Interim pulse height tally storage
vector<double>& pht_storage() { return pht_storage_; }
// Global tally accumulators
double& keff_tally_absorption() { return keff_tally_absorption_; }
double& keff_tally_collision() { return keff_tally_collision_; }
double& keff_tally_tracklength() { return keff_tally_tracklength_; }
double& keff_tally_leakage() { return keff_tally_leakage_; }
// Shows debug info
bool& trace() { return trace_; }
// Distance to the next collision
double& collision_distance() { return collision_distance_; }
// Number of events particle has undergone
int& n_event() { return n_event_; }
// Number of times variance reduction has caused a particle split
int n_split() const { return n_split_; }
int& n_split() { return n_split_; }
// Particle-specific factor for on-the-fly weight window adjustment
double ww_factor() const { return ww_factor_; }
double& ww_factor() { return ww_factor_; }
#ifdef DAGMC
moab::DagMC::RayHistory& history() { return history_; }
Direction& last_dir() { return last_dir_; }
#endif
// Number of progeny produced by this particle
int64_t& n_progeny() { return n_progeny_; }
// Accessors for position in global coordinates
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
// Accessors for position in local coordinates
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
// Accessors for direction in global coordinates
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
// Accessors for direction in local coordinates
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
//! Gets the pointer to the particle's current PRN seed
uint64_t* current_seed() { return seeds_ + stream_; }
const uint64_t* current_seed() const { return seeds_ + stream_; }
@ -512,14 +619,6 @@ public:
micro.last_E = 0.0;
}
//! resets all coordinate levels for the particle
void clear()
{
for (auto& level : coord_)
level.reset();
n_coord_ = 1;
}
//! Get track information based on particle's current state
TrackState get_track_state() const;

15
include/openmc/plot.h
View file

@ -1,6 +1,7 @@
#ifndef OPENMC_PLOT_H
#define OPENMC_PLOT_H
#include <cmath>
#include <sstream>
#include <unordered_map>
@ -99,8 +100,8 @@ public:
virtual void print_info() const = 0;
const std::string& path_plot() const { return path_plot_; }
const int id() const { return id_; }
const int level() const { return level_; }
int id() const { return id_; }
int level() const { return level_; }
// Public color-related data
PlottableInterface(pugi::xml_node plot_node);
@ -120,7 +121,7 @@ struct IdData {
IdData(size_t h_res, size_t v_res);
// Methods
void set_value(size_t y, size_t x, const Particle& p, int level);
void set_value(size_t y, size_t x, const GeometryState& p, int level);
void set_overlap(size_t y, size_t x);
// Members
@ -132,7 +133,7 @@ struct PropertyData {
PropertyData(size_t h_res, size_t v_res);
// Methods
void set_value(size_t y, size_t x, const Particle& p, int level);
void set_value(size_t y, size_t x, const GeometryState& p, int level);
void set_overlap(size_t y, size_t x);
// Members
@ -200,11 +201,11 @@ T SlicePlotBase::get_map() const
xyz[out_i] = origin_[out_i] + width_[1] / 2. - out_pixel / 2.;
// arbitrary direction
Direction dir = {0.7071, 0.7071, 0.0};
Direction dir = {1. / std::sqrt(2.), 1. / std::sqrt(2.), 0.0};
#pragma omp parallel
{
Particle p;
GeometryState p;
p.r() = xyz;
p.u() = dir;
p.coord(0).universe = model::root_universe;
@ -290,7 +291,7 @@ private:
* find a distance to the boundary in a non-standard surface intersection
* check. It's an exhaustive search over surfaces in the top-level universe.
*/
static int advance_to_boundary_from_void(Particle& p);
static int advance_to_boundary_from_void(GeometryState& p);
/* Checks if a vector of two TrackSegments is equivalent. We define this
* to mean not having matching intersection lengths, but rather having

11
include/openmc/settings.h
View file

@ -59,8 +59,12 @@ extern bool trigger_predict; //!< predict batches for triggers?
extern bool ufs_on; //!< uniform fission site method on?
extern bool urr_ptables_on; //!< use unresolved resonance prob. tables?
extern "C" bool weight_windows_on; //!< are weight windows are enabled?
extern bool write_all_tracks; //!< write track files for every particle?
extern bool write_initial_source; //!< write out initial source file?
extern bool weight_window_checkpoint_surface; //!< enable weight window check
//!< upon surface crossing?
extern bool weight_window_checkpoint_collision; //!< enable weight window check
//!< upon collision?
extern bool write_all_tracks; //!< write track files for every particle?
extern bool write_initial_source; //!< write out initial source file?
// Paths to various files
extern std::string path_cross_sections; //!< path to cross_sections.xml
@ -82,6 +86,9 @@ extern "C" int32_t max_lost_particles; //!< maximum number of lost particles
extern double
rel_max_lost_particles; //!< maximum number of lost particles, relative to the
//!< total number of particles
extern "C" int32_t
max_write_lost_particles; //!< maximum number of lost particles
//!< to be written to files
extern "C" int32_t gen_per_batch; //!< number of generations per batch
extern "C" int64_t n_particles; //!< number of particles per generation

42
include/openmc/source.h
View file

@ -50,6 +50,8 @@ public:
// Methods that can be overridden
virtual double strength() const { return 1.0; }
static unique_ptr<Source> create(pugi::xml_node node);
};
//==============================================================================
@ -101,12 +103,13 @@ private:
class FileSource : public Source {
public:
// Constructors
explicit FileSource(std::string path);
explicit FileSource(const vector<SourceSite>& sites) : sites_ {sites} {}
explicit FileSource(pugi::xml_node node);
explicit FileSource(const std::string& path);
// Methods
SourceSite sample(uint64_t* seed) const override;
void load_sites_from_file(
const std::string& path); //!< Load source sites from file
private:
vector<SourceSite> sites_; //!< Source sites from a file
};
@ -118,7 +121,7 @@ private:
class CompiledSourceWrapper : public Source {
public:
// Constructors, destructors
CompiledSourceWrapper(std::string path, std::string parameters);
CompiledSourceWrapper(pugi::xml_node node);
~CompiledSourceWrapper();
// Defer implementation to custom source library
@ -129,6 +132,8 @@ public:
double strength() const override { return compiled_source_->strength(); }
void setup(const std::string& path, const std::string& parameters);
private:
void* shared_library_; //!< library from dlopen
unique_ptr<Source> compiled_source_;
@ -136,6 +141,35 @@ private:
typedef unique_ptr<Source> create_compiled_source_t(std::string parameters);
//==============================================================================
//! Mesh-based source with different distributions for each element
//==============================================================================
class MeshSource : public Source {
public:
// Constructors
explicit MeshSource(pugi::xml_node node);
//! Sample from the external source distribution
//! \param[inout] seed Pseudorandom seed pointer
//! \return Sampled site
SourceSite sample(uint64_t* seed) const override;
// Properties
double strength() const override { return space_->total_strength(); }
// Accessors
const std::unique_ptr<Source>& source(int32_t i) const
{
return sources_.size() == 1 ? sources_[0] : sources_[i];
}
private:
// Data members
unique_ptr<MeshSpatial> space_; //!< Mesh spatial
vector<std::unique_ptr<Source>> sources_; //!< Source distributions
};
//==============================================================================
// Functions
//==============================================================================

15
include/openmc/surface.h
View file

@ -83,10 +83,10 @@ struct BoundingBox {
class Surface {
public:
int id_; //!< Unique ID
std::string name_; //!< User-defined name
std::shared_ptr<BoundaryCondition> bc_ {nullptr}; //!< Boundary condition
GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC)
int id_; //!< Unique ID
std::string name_; //!< User-defined name
unique_ptr<BoundaryCondition> bc_; //!< Boundary condition
GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC)
bool surf_source_ {false}; //!< Activate source banking for the surface?
explicit Surface(pugi::xml_node surf_node);
@ -105,12 +105,13 @@ public:
//! Determine the direction of a ray reflected from the surface.
//! \param[in] r The point at which the ray is incident.
//! \param[in] u Incident direction of the ray
//! \param[inout] p Pointer to the particle
//! \param[inout] p Pointer to the particle. Only DAGMC uses this.
//! \return Outgoing direction of the ray
virtual Direction reflect(Position r, Direction u, Particle* p) const;
virtual Direction reflect(
Position r, Direction u, GeometryState* p = nullptr) const;
virtual Direction diffuse_reflect(
Position r, Direction u, uint64_t* seed) const;
Position r, Direction u, uint64_t* seed, GeometryState* p = nullptr) const;
//! Evaluate the equation describing the surface.
//!

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

@ -31,6 +31,7 @@ enum class FilterType {
ENERGY_OUT,
LEGENDRE,
MATERIAL,
MATERIALFROM,
MESH,
MESH_SURFACE,
MU,

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

@ -46,7 +46,7 @@ public:
void set_materials(gsl::span<const int32_t> materials);
private:
protected:
//----------------------------------------------------------------------------
// Data members

29
include/openmc/tallies/filter_materialfrom.h Normal file
View file

@ -0,0 +1,29 @@
#ifndef OPENMC_TALLIES_FILTER_MATERIALFROM_H
#define OPENMC_TALLIES_FILTER_MATERIALFROM_H
#include <string>
#include "openmc/tallies/filter_material.h"
namespace openmc {
//==============================================================================
//! Specifies which material particles exit when crossing a surface.
//==============================================================================
class MaterialFromFilter : public MaterialFilter {
public:
//----------------------------------------------------------------------------
// Methods
std::string type_str() const override { return "materialfrom"; }
FilterType type() const override { return FilterType::MATERIALFROM; }
void get_all_bins(const Particle& p, TallyEstimator estimator,
FilterMatch& match) const override;
std::string text_label(int bin) const override;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_MATERIALFROM_H

3
include/openmc/universe.h
View file

@ -9,6 +9,7 @@ namespace openmc {
class DAGUniverse;
#endif
class GeometryState;
class Universe;
class UniversePartitioner;
@ -32,7 +33,7 @@ public:
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const;
virtual bool find_cell(Particle& p) const;
virtual bool find_cell(GeometryState& p) const;
BoundingBox bounding_box() const;

36
include/openmc/volume_calc.h
View file

@ -1,18 +1,23 @@
#ifndef OPENMC_VOLUME_CALC_H
#define OPENMC_VOLUME_CALC_H
#include <algorithm> // for find
#include <cstdint>
#include <string>
#include <vector>
#include "openmc/array.h"
#include "openmc/openmp_interface.h"
#include "openmc/position.h"
#include "openmc/tallies/trigger.h"
#include "openmc/vector.h"
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
#include <gsl/gsl-lite.hpp>
#ifdef _OPENMP
#include <omp.h>
#endif
namespace openmc {
@ -89,6 +94,35 @@ extern vector<VolumeCalculation> volume_calcs;
// Non-member functions
//==============================================================================
//! Reduce vector of indices and hits from each thread to a single copy
//
//! \param[in] local_indices Indices specific to each thread
//! \param[in] local_hits Hit count specific to each thread
//! \param[out] indices Reduced vector of indices
//! \param[out] hits Reduced vector of hits
template<typename T, typename T2>
void reduce_indices_hits(const vector<T>& local_indices,
const vector<T2>& local_hits, vector<T>& indices, vector<T2>& hits)
{
const int n_threads = num_threads();
#pragma omp for ordered schedule(static, 1)
for (int i = 0; i < n_threads; ++i) {
#pragma omp ordered
for (int j = 0; j < local_indices.size(); ++j) {
// Check if this material has been added to the master list and if
// so, accumulate the number of hits
auto it = std::find(indices.begin(), indices.end(), local_indices[j]);
if (it == indices.end()) {
indices.push_back(local_indices[j]);
hits.push_back(local_hits[j]);
} else {
hits[it - indices.begin()] += local_hits[j];
}
}
}
}
void free_memory_volume();
} // namespace openmc

2
man/man1/openmc.1
View file

@ -63,7 +63,7 @@ Indicates the default path to an HDF5 file that contains multi-group cross
section libraries if the user has not specified the <cross_sections> tag in
.I materials.xml\fP.
.SH LICENSE
Copyright \(co 2011-2023 Massachusetts Institute of Technology, UChicago
Copyright \(co 2011-2024 Massachusetts Institute of Technology, UChicago
Argonne LLC, and OpenMC contributors.
.PP
Permission is hereby granted, free of charge, to any person obtaining a copy of

4
openmc/__init__.py
View file

@ -36,7 +36,7 @@ from . import examples
from .config import *
# Import a few names from the model module
from openmc.model import rectangular_prism, hexagonal_prism, Model
from openmc.model import Model
__version__ = '0.13.4-dev'
__version__ = '0.14.1-dev'

110
openmc/bounding_box.py
View file

@ -1,3 +1,4 @@
from __future__ import annotations
from typing import Iterable
import numpy as np
@ -5,10 +6,10 @@ import numpy as np
from .checkvalue import check_length
class BoundingBox(tuple):
class BoundingBox:
"""Axis-aligned bounding box.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -34,17 +35,70 @@ class BoundingBox(tuple):
The width of the x, y and z axis in [cm]
"""
def __new__(cls, lower_left: Iterable[float], upper_right: Iterable[float]):
def __init__(self, lower_left: Iterable[float], upper_right: Iterable[float]):
check_length("lower_left", lower_left, 3, 3)
check_length("upper_right", upper_right, 3, 3)
lower_left = np.array(lower_left, dtype=float)
upper_right = np.array(upper_right, dtype=float)
return tuple.__new__(cls, (lower_left, upper_right))
self._bounds = np.asarray([lower_left, upper_right], dtype=float)
def __repr__(self) -> str:
return "BoundingBox(lower_left={}, upper_right={})".format(
tuple(self.lower_left), tuple(self.upper_right))
def __getitem__(self, key) -> np.ndarray:
return self._bounds[key]
def __len__(self):
return 2
def __setitem__(self, key, val):
self._bounds[key] = val
def __iand__(self, other: BoundingBox) -> BoundingBox:
"""Updates the box be the intersection of itself and another box
Parameters
----------
other : BoundingBox
The box used to resize this box
Returns
-------
An updated bounding box
"""
self.lower_left = np.maximum(self.lower_left, other.lower_left)
self.upper_right = np.minimum(self.upper_right, other.upper_right)
return self
def __and__(self, other: BoundingBox) -> BoundingBox:
new = BoundingBox(*self)
new &= other
return new
def __ior__(self, other: BoundingBox) -> BoundingBox:
"""Updates the box be the union of itself and another box
Parameters
----------
other : BoundingBox
The box used to resize this box
Returns
-------
An updated bounding box
"""
self.lower_left = np.minimum(self.lower_left, other.lower_left)
self.upper_right = np.maximum(self.upper_right, other.upper_right)
return self
def __or__(self, other: BoundingBox) -> BoundingBox:
new = BoundingBox(*self)
new |= other
return new
def __contains__(self, point):
"""Check whether or not a point is in the bounding box"""
return all(point > self.lower_left) and all(point < self.upper_right)
@property
def center(self) -> np.ndarray:
return (self[0] + self[1]) / 2
@ -53,10 +107,20 @@ class BoundingBox(tuple):
def lower_left(self) -> np.ndarray:
return self[0]
@lower_left.setter
def lower_left(self, llc):
check_length('lower_left', llc, 3, 3)
self[0] = llc
@property
def upper_right(self) -> np.ndarray:
return self[1]
@upper_right.setter
def upper_right(self, urc):
check_length('upper_right', urc, 3, 3)
self[1] = urc
@property
def volume(self) -> float:
return np.abs(np.prod(self[1] - self[0]))
@ -87,3 +151,37 @@ class BoundingBox(tuple):
@property
def width(self):
return self.upper_right - self.lower_left
def expand(self, padding_distance: float, inplace: bool = False) -> BoundingBox:
"""Returns an expanded bounding box
Parameters
----------
padding_distance : float
The distance to enlarge the bounding box by
inplace : bool
Whether or not to return a new BoundingBox instance or to modify the
current BoundingBox object.
Returns
-------
An expanded bounding box
"""
if inplace:
self[0] -= padding_distance
self[1] += padding_distance
return self
else:
return BoundingBox(self[0] - padding_distance, self[1] + padding_distance)
@classmethod
def infinite(cls) -> BoundingBox:
"""Create an infinite box. Useful as a starting point for determining
geometry bounds.
Returns
-------
An infinitely large bounding box.
"""
infs = np.full((3,), np.inf)
return cls(-infs, infs)

4
openmc/cell.py
View file

@ -1,8 +1,8 @@
from collections.abc import Iterable
from math import cos, sin, pi
from numbers import Real
import lxml.etree as ET
import lxml.etree as ET
import numpy as np
from uncertainties import UFloat
@ -562,7 +562,7 @@ class Cell(IDManagerMixin):
def plot(self, *args, **kwargs):
"""Display a slice plot of the cell.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------

9
openmc/cmfd.py
View file

@ -1308,9 +1308,6 @@ class CMFDRun:
Whether or not to run an adjoint calculation
"""
# Check for physical adjoint
physical_adjoint = adjoint and self._adjoint_type == 'physical'
# Start timer for build
time_start_buildcmfd = time.time()
@ -1432,9 +1429,6 @@ class CMFDRun:
# indices of the actual problem so that cmfd_flux can be multiplied by
# nfissxs
# Calculate volume
vol = np.prod(self._hxyz, axis=3)
# Reshape phi by number of groups
phi = self._phi.reshape((n, ng))
@ -2148,7 +2142,6 @@ class CMFDRun:
is_accel = self._coremap != _CMFD_NOACCEL
# Logical for determining whether a zero flux "albedo" b.c. should be
# applied
is_zero_flux_alb = abs(self._albedo - _ZERO_FLUX) < _TINY_BIT
x_inds, y_inds, z_inds = np.indices((nx, ny, nz))
# Define slice equivalent to is_accel[0,:,:]
@ -2725,8 +2718,6 @@ class CMFDRun:
# Define flux in each cell
cell_flux = self._flux / dxdydz
# Extract indices of coremap that are accelerated
is_accel = self._coremap != _CMFD_NOACCEL
# Define dhat at left surface for all mesh cells on left boundary
boundary = self._first_x_accel

5
openmc/config.py
View file

@ -4,7 +4,7 @@ from pathlib import Path
import warnings
from openmc.data import DataLibrary
from openmc.data.decay import _DECAY_PHOTON_ENERGY
from openmc.data.decay import _DECAY_ENERGY, _DECAY_PHOTON_ENERGY
__all__ = ["config"]
@ -41,6 +41,7 @@ class _Config(MutableMapping):
os.environ['OPENMC_CHAIN_FILE'] = str(value)
# Reset photon source data since it relies on chain file
_DECAY_PHOTON_ENERGY.clear()
_DECAY_ENERGY.clear()
else:
raise KeyError(f'Unrecognized config key: {key}. Acceptable keys '
'are "cross_sections", "mg_cross_sections" and '
@ -76,7 +77,7 @@ def _default_config():
if (chain_file is None and
config.get('cross_sections') is not None and
config['cross_sections'].exists()
):
):
# Check for depletion chain in cross_sections.xml
data = DataLibrary.from_xml(config['cross_sections'])
for lib in reversed(data.libraries):

5
openmc/data/angle_distribution.py
View file

@ -235,7 +235,6 @@ class AngleDistribution(EqualityMixin):
items = get_cont_record(file_obj)
li = items[2]
nk = items[4]
center_of_mass = (items[3] == 2)
# Check for obsolete energy transformation matrix. If present, just skip
# it and keep reading
@ -259,7 +258,6 @@ class AngleDistribution(EqualityMixin):
mu = []
for i in range(n_energy):
items, al = get_list_record(file_obj)
temperature = items[0]
energy[i] = items[1]
coefficients = np.asarray([1.0] + al)
mu.append(Legendre(coefficients))
@ -273,7 +271,6 @@ class AngleDistribution(EqualityMixin):
mu = []
for i in range(n_energy):
params, f = get_tab1_record(file_obj)
temperature = params[0]
energy[i] = params[1]
if f.n_regions > 1:
raise NotImplementedError('Angular distribution with multiple '
@ -289,7 +286,6 @@ class AngleDistribution(EqualityMixin):
mu = []
for i in range(n_energy_legendre):
items, al = get_list_record(file_obj)
temperature = items[0]
energy_legendre[i] = items[1]
coefficients = np.asarray([1.0] + al)
mu.append(Legendre(coefficients))
@ -300,7 +296,6 @@ class AngleDistribution(EqualityMixin):
energy_tabulated = np.zeros(n_energy_tabulated)
for i in range(n_energy_tabulated):
params, f = get_tab1_record(file_obj)
temperature = params[0]
energy_tabulated[i] = params[1]
if f.n_regions > 1:
raise NotImplementedError('Angular distribution with multiple '

2
openmc/data/decay.py
View file

@ -599,7 +599,7 @@ def decay_photon_energy(nuclide: str) -> Optional[Univariate]:
openmc.stats.Univariate or None
Distribution of energies in [eV] of photons emitted from decay, or None
if no photon source exists. Note that the probabilities represent
intensities, given as [decay/sec].
intensities, given as [Bq].
"""
if not _DECAY_PHOTON_ENERGY:
chain_file = openmc.config.get('chain_file')

4
openmc/data/library.py
View file

@ -1,8 +1,8 @@
import os
import lxml.etree as ET
import pathlib
import h5py
import lxml.etree as ET
import openmc
from openmc._xml import clean_indentation, reorder_attributes
@ -15,7 +15,7 @@ class DataLibrary(list):
cross section data from a single file. The dictionary has keys 'path',
'type', and 'materials'.
.. versionchanged:: 0.13.4
.. versionchanged:: 0.14.0
This class now behaves like a list rather than requiring you to access
the list of libraries through a special attribute.

1
openmc/data/reaction.py
View file

@ -109,7 +109,6 @@ def _get_products(ev, mt):
za = int(params[0])
awr = params[1]
lip = params[2]
law = params[3]
if za == 0:

6
openmc/data/resonance.py
View file

@ -93,9 +93,8 @@ class Resonances:
n_isotope = items[4] # Number of isotopes
ranges = []
for iso in range(n_isotope):
for _ in range(n_isotope):
items = get_cont_record(file_obj)
abundance = items[1]
fission_widths = (items[3] == 1) # fission widths are given?
n_ranges = items[4] # number of resonance energy ranges
@ -424,14 +423,12 @@ class MultiLevelBreitWigner(ResonanceRange):
# Determine penetration and shift corresponding to resonance energy
k = wave_number(A, E)
rho = k*self.channel_radius[l](E)
rhohat = k*self.scattering_radius[l](E)
p[i], s[i] = penetration_shift(l, rho)
# Determine penetration at modified energy for competitive reaction
if gx > 0:
Ex = E + self.q_value[l]*(A + 1)/A
rho = k*self.channel_radius[l](Ex)
rhohat = k*self.scattering_radius[l](Ex)
px[i], sx[i] = penetration_shift(l, rho)
else:
px[i] = sx[i] = 0.0
@ -680,7 +677,6 @@ class ReichMoore(ResonanceRange):
# Determine penetration and shift corresponding to resonance energy
k = wave_number(A, E)
rho = k*self.channel_radius[l](E)
rhohat = k*self.scattering_radius[l](E)
p[i], s[i] = penetration_shift(l, rho)
df['p'] = p

14
openmc/data/resonance_covariance.py
View file

@ -93,10 +93,8 @@ class ResonanceCovariances(Resonances):
n_isotope = items[4] # Number of isotopes
ranges = []
for iso in range(n_isotope):
for _ in range(n_isotope):
items = endf.get_cont_record(file_obj)
abundance = items[1]
fission_widths = (items[3] == 1) # Flag for fission widths
n_ranges = items[4] # Number of resonance energy ranges
for j in range(n_ranges):
@ -378,7 +376,6 @@ class MultiLevelBreitWignerCovariance(ResonanceCovarianceRange):
# Other scatter radius parameters
items = endf.get_cont_record(file_obj)
target_spin = items[0]
lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form
nls = items[4] # number of l-values
@ -387,8 +384,6 @@ class MultiLevelBreitWignerCovariance(ResonanceCovarianceRange):
items = endf.get_cont_record(file_obj)
# Number of short range type resonance covariances
num_short_range = items[4]
# Number of long range type resonance covariances
num_long_range = items[5]
# Read resonance widths, J values, etc
records = []
@ -421,7 +416,6 @@ class MultiLevelBreitWignerCovariance(ResonanceCovarianceRange):
# compact correlations
elif lcomp == 2:
items, values = endf.get_list_record(file_obj)
mean = items
num_res = items[5]
energy = values[0::12]
spin = values[1::12]
@ -613,20 +607,14 @@ class ReichMooreCovariance(ResonanceCovarianceRange):
# Other scatter radius parameters
items = endf.get_cont_record(file_obj)
target_spin = items[0]
lcomp = items[3] # Flag for compatibility 0, 1, 2 - 2 is compact form
nls = items[4] # Number of l-values
# Build covariance matrix for General Resolved Resonance Formats
if lcomp == 1:
items = endf.get_cont_record(file_obj)
# Number of short range type resonance covariances
num_short_range = items[4]
# Number of long range type resonance covariances
num_long_range = items[5]
# Read resonance widths, J values, etc
channel_radius = {}
scattering_radius = {}
records = []
for i in range(num_short_range):
items, values = endf.get_list_record(file_obj)

2
openmc/deplete/_matrix_funcs.py
View file

@ -77,4 +77,4 @@ def leqi_f4(chain, inputs, fission_yields=None):
return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1
+ (dt ** 2 + 2 * dt * dt_l + dt_l ** 2)
/ (12 * dt_l * (dt + dt_l)) * f2
+ (4 * dt * dt_l + 5 * dt_l ** 2) / (12 * dt_l * (dt + dt_l)) * f3)
+ (4 * dt + 5 * dt_l) / (12 * (dt + dt_l)) * f3)

40
openmc/deplete/abc.py
View file

@ -18,7 +18,7 @@ from warnings import warn
from numpy import nonzero, empty, asarray
from uncertainties import ufloat
from openmc.checkvalue import check_type, check_greater_than, check_value
from openmc.checkvalue import checkvalue, check_type, check_greater_than, PathLike
from openmc.mpi import comm
from .stepresult import StepResult
from .chain import Chain
@ -70,8 +70,9 @@ def change_directory(output_dir):
Directory to switch to.
"""
orig_dir = os.getcwd()
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
try:
output_dir.mkdir(parents=True, exist_ok=True)
os.chdir(output_dir)
yield
finally:
@ -545,7 +546,7 @@ class Integrator(ABC):
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
* ``A`` is a :class:`scipy.sparse.csc_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
@ -559,7 +560,7 @@ class Integrator(ABC):
Instance of Batchwise class to perform batch-wise scheme during
transport-depletion simulation.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
"""
@ -776,7 +777,12 @@ class Integrator(ABC):
x, root = self._batchwise.search_for_keff(x, step_index)
return x, root
def integrate(self, final_step=True, output=True):
def integrate(
self,
final_step: bool = True,
output: bool = True,
path: PathLike = 'depletion_results.h5'
):
"""Perform the entire depletion process across all steps
Parameters
@ -790,6 +796,10 @@ class Integrator(ABC):
Indicate whether to display information about progress
.. versionadded:: 0.13.1
path : PathLike
Path to file to write. Defaults to 'depletion_results.h5'.
.. versionadded:: 0.14.1
"""
with change_directory(self.operator.output_dir):
n = self.operator.initial_condition()
@ -820,7 +830,7 @@ class Integrator(ABC):
# Remove actual EOS concentration for next step
n = n_list.pop()
StepResult.save(self.operator, n_list, res_list, [t, t + dt],
source_rate, self._i_res + i, proc_time, root)
source_rate, self._i_res + i, proc_time, root, path)
t += dt
@ -969,7 +979,7 @@ class SIIntegrator(Integrator):
User-supplied functions are expected to have the following signature:
``solver(A, n0, t) -> n1`` where
* ``A`` is a :class:`scipy.sparse.csr_matrix` making up the
* ``A`` is a :class:`scipy.sparse.csc_matrix` making up the
depletion matrix
* ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
for a given material in atoms/cm3
@ -1002,15 +1012,21 @@ class SIIntegrator(Integrator):
self.operator.settings.particles //= self.n_steps
return inherited
def integrate(self, output=True):
def integrate(
self,
output: bool = True,
path: PathLike = "depletion_results.h5"
):
"""Perform the entire depletion process across all steps
Parameters
----------
output : bool, optional
Indicate whether to display information about progress
path : PathLike
Path to file to write. Defaults to 'depletion_results.h5'.
.. versionadded:: 0.13.1
.. versionadded:: 0.14.1
"""
with change_directory(self.operator.output_dir):
n = self.operator.initial_condition()
@ -1040,13 +1056,13 @@ class SIIntegrator(Integrator):
n = n_list.pop()
StepResult.save(self.operator, n_list, res_list, [t, t + dt],
p, self._i_res + i, proc_time)
p, self._i_res + i, proc_time, path)
t += dt
# No final simulation for SIE, use last iteration results
StepResult.save(self.operator, [n], [res_list[-1]], [t, t],
p, self._i_res + len(self), proc_time)
p, self._i_res + len(self), proc_time, path)
self.operator.write_bos_data(self._i_res + len(self))
self.operator.finalize()
@ -1071,7 +1087,7 @@ class DepSystemSolver(ABC):
Parameters
----------
A : scipy.sparse.csr_matrix
A : scipy.sparse.csc_matrix
Sparse transmutation matrix ``A[j, i]`` describing rates at
which isotope ``i`` transmutes to isotope ``j``
n0 : numpy.ndarray

12
openmc/deplete/chain.py
View file

@ -596,7 +596,7 @@ class Chain:
Returns
-------
scipy.sparse.csr_matrix
scipy.sparse.csc_matrix
Sparse matrix representing depletion.
See Also
@ -680,16 +680,16 @@ class Chain:
# Clear set of reactions
reactions.clear()
# Use DOK matrix as intermediate representation, then convert to CSR and return
# Use DOK matrix as intermediate representation, then convert to CSC and return
n = len(self)
matrix_dok = sp.dok_matrix((n, n))
dict.update(matrix_dok, matrix)
return matrix_dok.tocsr()
return matrix_dok.tocsc()
def form_rr_term(self, tr_rates, mats):
"""Function to form the transfer rate term matrices.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -713,7 +713,7 @@ class Chain:
Returns
-------
scipy.sparse.csr_matrix
scipy.sparse.csc_matrix
Sparse matrix representing transfer term.
"""
@ -746,7 +746,7 @@ class Chain:
n = len(self)
matrix_dok = sp.dok_matrix((n, n))
dict.update(matrix_dok, matrix)
return matrix_dok.tocsr()
return matrix_dok.tocsc()
def get_branch_ratios(self, reaction="(n,gamma)"):
"""Return a dictionary with reaction branching ratios

133
openmc/deplete/coupled_operator.py
View file

@ -10,6 +10,7 @@ filesystem.
import copy
from warnings import warn
from typing import Optional
import numpy as np
from uncertainties import ufloat
@ -33,18 +34,19 @@ from .helpers import (
__all__ = ["CoupledOperator", "Operator", "OperatorResult"]
def _find_cross_sections(model):
def _find_cross_sections(model: Optional[str] = None):
"""Determine cross sections to use for depletion
Parameters
----------
model : openmc.model.Model
model : openmc.model.Model, optional
Reactor model
"""
if model.materials and model.materials.cross_sections is not None:
# Prefer info from Model class if available
return model.materials.cross_sections
if model:
if model.materials and model.materials.cross_sections is not None:
# Prefer info from Model class if available
return model.materials.cross_sections
# otherwise fallback to environment variable
cross_sections = openmc.config.get("cross_sections")
@ -67,7 +69,7 @@ def _get_nuclides_with_data(cross_sections):
Returns
-------
nuclides : set of str
Set of nuclide names that have cross secton data
Set of nuclide names that have cross section data
"""
nuclides = set()
@ -164,15 +166,18 @@ class CoupledOperator(OpenMCOperator):
``None`` implies no limit on the depth.
.. versionadded:: 0.12
diff_volume_method : str
Specifies how the volumes of the new materials should be found. Default
is to 'divide equally' which divides the original material volume
equally between the new materials, 'match cell' sets the volume of the
material to volume of the cell they fill.
.. versionadded:: 0.14.0
Attributes
----------
model : openmc.model.Model
OpenMC model object
geometry : openmc.Geometry
OpenMC geometry object
settings : openmc.Settings
OpenMC settings object
output_dir : pathlib.Path
Path to output directory to save results.
round_number : bool
@ -206,8 +211,8 @@ class CoupledOperator(OpenMCOperator):
}
def __init__(self, model, chain_file=None, prev_results=None,
diff_burnable_mats=False, normalization_mode="fission-q",
fission_q=None,
diff_burnable_mats=False, diff_volume_method="divide equally",
normalization_mode="fission-q", fission_q=None,
fission_yield_mode="constant", fission_yield_opts=None,
reaction_rate_mode="direct", reaction_rate_opts=None,
reduce_chain=False, reduce_chain_level=None):
@ -233,8 +238,6 @@ class CoupledOperator(OpenMCOperator):
warn("Fission Q dictionary will not be used")
fission_q = None
self.model = model
self.settings = model.settings
self.geometry = model.geometry
# determine set of materials in the model
if not model.materials:
@ -256,44 +259,26 @@ class CoupledOperator(OpenMCOperator):
'fission_yield_opts': fission_yield_opts
}
# Records how many times the operator has been called
self._n_calls = 0
super().__init__(
model.materials,
cross_sections,
chain_file,
prev_results,
diff_burnable_mats,
fission_q,
helper_kwargs,
reduce_chain,
reduce_chain_level)
materials=model.materials,
cross_sections=cross_sections,
chain_file=chain_file,
prev_results=prev_results,
diff_burnable_mats=diff_burnable_mats,
diff_volume_method=diff_volume_method,
fission_q=fission_q,
helper_kwargs=helper_kwargs,
reduce_chain=reduce_chain,
reduce_chain_level=reduce_chain_level)
def _differentiate_burnable_mats(self):
"""Assign distribmats for each burnable material"""
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Extract all burnable materials which have multiple instances
distribmats = set(
[mat for mat in self.materials
if mat.depletable and mat.num_instances > 1])
for mat in distribmats:
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
"material with ID={}.".format(mat.id))
mat.volume /= mat.num_instances
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
cell.fill = [mat.clone()
for i in range(cell.num_instances)]
self.materials = openmc.Materials(
self.model.geometry.get_all_materials().values()
self.model.differentiate_depletable_mats(
diff_volume_method=self.diff_volume_method
)
def _load_previous_results(self):
@ -367,7 +352,7 @@ class CoupledOperator(OpenMCOperator):
if normalization_mode == "fission-q":
self._normalization_helper = ChainFissionHelper()
elif normalization_mode == "energy-deposition":
score = "heating" if self.settings.photon_transport else "heating-local"
score = "heating" if self.model.settings.photon_transport else "heating-local"
self._normalization_helper = EnergyScoreHelper(score)
else:
self._normalization_helper = SourceRateHelper()
@ -389,8 +374,12 @@ class CoupledOperator(OpenMCOperator):
# Create XML files
if comm.rank == 0:
self.geometry.export_to_xml()
self.settings.export_to_xml()
self.model.geometry.export_to_xml()
self.model.settings.export_to_xml()
if self.model.plots:
self.model.plots.export_to_xml()
if self.model.tallies:
self.model.tallies.export_to_xml()
self._generate_materials_xml()
# Initialize OpenMC library
@ -441,6 +430,16 @@ class CoupledOperator(OpenMCOperator):
# Reset results in OpenMC
openmc.lib.reset()
# The timers are reset only if the operator has been called before.
# This is because we call this method after loading cross sections, and
# no transport has taken place yet. As a result, we only reset the
# timers after the first step so as to correctly report the time spent
# reading cross sections in the first depletion step, and from there
# correctly report all particle tracking rates in multistep depletion
# solvers.
if self._n_calls > 0:
openmc.lib.reset_timers()
self._update_materials_and_nuclides(vec)
# If the source rate is zero, return zero reaction rates without running
@ -461,6 +460,8 @@ class CoupledOperator(OpenMCOperator):
op_result = OperatorResult(keff, rates)
self._n_calls += 1
return copy.deepcopy(op_result)
def _update_materials(self):
@ -520,13 +521,41 @@ class CoupledOperator(OpenMCOperator):
"openmc_simulation_n{}.h5".format(step),
write_source=False)
openmc.lib.reset_timers()
def finalize(self):
"""Finalize a depletion simulation and release resources."""
if self.cleanup_when_done:
openmc.lib.finalize()
# The next few class variables and methods should be removed after one
# release cycle or so. For now, we will provide compatibility to
# accessing CoupledOperator.settings and CoupledOperator.geometry. In
# the future these should stay on the Model class.
var_warning_msg = "The CoupledOperator.{0} variable should be \
accessed through CoupledOperator.model.{0}."
geometry_warning_msg = var_warning_msg.format("geometry")
settings_warning_msg = var_warning_msg.format("settings")
@property
def settings(self):
warn(self.settings_warning_msg, FutureWarning)
return self.model.settings
@settings.setter
def settings(self, new_settings):
warn(self.settings_warning_msg, FutureWarning)
self.model.settings = new_settings
@property
def geometry(self):
warn(self.geometry_warning_msg, FutureWarning)
return self.model.geometry
@geometry.setter
def geometry(self, new_geometry):
warn(self.geometry_warning_msg, FutureWarning)
self.model.geometry = new_geometry
# Retain deprecated name for the time being
def Operator(*args, **kwargs):

4
openmc/deplete/cram.py
View file

@ -75,9 +75,9 @@ class IPFCramSolver(DepSystemSolver):
Final compositions after ``dt``
"""
A = sp.csr_matrix(A * dt, dtype=np.float64)
A = dt * sp.csc_matrix(A, dtype=np.float64)
y = n0.copy()
ident = sp.eye(A.shape[0])
ident = sp.eye(A.shape[0], format='csc')
for alpha, theta in zip(self.alpha, self.theta):
y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
return y * self.alpha0

6
openmc/deplete/helpers.py
View file

@ -268,6 +268,11 @@ class FluxCollapseHelper(ReactionRateHelper):
if self._reactions_direct and self._nuclides_direct is None:
self._rate_tally.nuclides = nuclides
# Make sure nuclide data is loaded
for nuclide in self.nuclides:
if nuclide not in openmc.lib.nuclides:
openmc.lib.load_nuclide(nuclide)
def generate_tallies(self, materials, scores):
"""Produce multigroup flux spectrum tally
@ -589,7 +594,6 @@ class ConstantFissionYieldHelper(FissionYieldHelper):
self._constant_yields[name] = yield_data
continue
# Specific energy not found, use closest energy
distances = [abs(energy - ene) for ene in nuc.yield_energies]
min_E = min(nuc.yield_energies, key=lambda e: abs(e - energy))
self._constant_yields[name] = nuc.yield_data[min_E]

18
openmc/deplete/independent_operator.py
View file

@ -39,7 +39,7 @@ class IndependentOperator(OpenMCOperator):
.. versionadded:: 0.13.1
.. versionchanged:: 0.13.4
.. versionchanged:: 0.14.0
Arguments updated to include list of fluxes and microscopic cross
sections.
@ -130,6 +130,14 @@ class IndependentOperator(OpenMCOperator):
# Validate micro-xs parameters
check_type('materials', materials, openmc.Materials)
check_type('micros', micros, Iterable, MicroXS)
check_type('fluxes', fluxes, Iterable, float)
if not (len(fluxes) == len(micros) == len(materials)):
msg = (f'The length of fluxes ({len(fluxes)}) should be equal to '
f'the length of micros ({len(micros)}) and the length of '
f'materials ({len(materials)}).')
raise ValueError(msg)
if keff is not None:
check_type('keff', keff, tuple, float)
keff = ufloat(*keff)
@ -148,10 +156,10 @@ class IndependentOperator(OpenMCOperator):
self.fluxes = fluxes
super().__init__(
materials,
micros,
chain_file,
prev_results,
materials=materials,
cross_sections=micros,
chain_file=chain_file,
prev_results=prev_results,
fission_q=fission_q,
helper_kwargs=helper_kwargs,
reduce_chain=reduce_chain,

5
openmc/deplete/integrators.py
View file

@ -137,7 +137,7 @@ class CF4Integrator(Integrator):
\mathbf{n}_{i+1} &= \exp \left ( \frac{\mathbf{A}_1}{4} + \frac{\mathbf{A}_2}{6}
+ \frac{\mathbf{A}_3}{6} - \frac{\mathbf{A}_4}{12} \right )
\exp \left ( -\frac{\mathbf{A}_1}{12} + \frac{\mathbf{A}_2}{6} +
\frac{\mathbf{A}_3}{6} - \frac{\mathbf{A}_4}{4} \right ) \mathbf{n}_i.
\frac{\mathbf{A}_3}{6} + \frac{\mathbf{A}_4}{4} \right ) \mathbf{n}_i.
\end{aligned}
"""
_num_stages = 4
@ -360,8 +360,7 @@ class LEQIIntegrator(Integrator):
h_i)} \mathbf{A}_0 + \frac{h_{i-1}}{12 (h_{i-1} + h_i)} \mathbf{A}_1 \\
\mathbf{F}_4 &= \frac{-h_i^2}{12 h_{i-1} (h_{i-1} + h_i)} \mathbf{A}_{-1} +
\frac{h_{i-1}^2 + 2 h_i h_{i-1} + h_i^2}{12 h_{i-1} (h_{i-1} + h_i)}
\mathbf{A}_0 + \frac{5 h_{i-1}^2 + 4 h_i h_{i-1}}{12 h_{i-1}
(h_{i-1} + h_i)} \mathbf{A}_1 \\
\mathbf{A}_0 + \frac{5 h_{i-1} + 4 h_i}{12 (h_{i-1} + h_i)} \mathbf{A}_1 \\
\mathbf{n}_{i+1} &= \exp(h_i \mathbf{F}_4) \exp(h_i \mathbf{F}_3) \mathbf{n}_i
\end{aligned}

164
openmc/deplete/microxs.py
View file

@ -5,8 +5,8 @@ IndependentOperator class for depletion.
"""
from __future__ import annotations
import tempfile
from typing import List, Tuple, Iterable, Optional, Union
from tempfile import TemporaryDirectory
from typing import List, Tuple, Iterable, Optional, Union, Sequence
import pandas as pd
import numpy as np
@ -14,14 +14,30 @@ import numpy as np
from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike
from openmc.exceptions import DataError
from openmc import StatePoint
from openmc.mgxs import GROUP_STRUCTURES
from openmc.data import REACTION_MT
import openmc
from .abc import change_directory
from .chain import Chain, REACTIONS
from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
import openmc.lib
_valid_rxns = list(REACTIONS)
_valid_rxns.append('fission')
def _resolve_chain_file_path(chain_file: str):
# Determine what reactions and nuclides are available in chain
if chain_file is None:
chain_file = openmc.config.get('chain_file')
if 'chain_file' in openmc.config:
raise DataError(
"No depletion chain specified and could not find depletion "
"chain in openmc.config['chain_file']"
)
return chain_file
def get_microxs_and_flux(
model: openmc.Model,
domains,
@ -33,13 +49,13 @@ def get_microxs_and_flux(
) -> Tuple[List[np.ndarray], List[MicroXS]]:
"""Generate a microscopic cross sections and flux from a Model
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
model : openmc.Model
OpenMC model object. Must contain geometry, materials, and settings.
domains : list of openmc.Material or openmc.Cell or openmc.Universe
domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase
Domains in which to tally reaction rates.
nuclides : list of str
Nuclides to get cross sections for. If not specified, all burnable
@ -69,13 +85,7 @@ def get_microxs_and_flux(
original_tallies = model.tallies
# Determine what reactions and nuclides are available in chain
if chain_file is None:
chain_file = openmc.config.get('chain_file')
if chain_file is None:
raise DataError(
"No depletion chain specified and could not find depletion "
"chain in openmc.config['chain_file']"
)
chain_file = _resolve_chain_file_path(chain_file)
chain = Chain.from_xml(chain_file)
if reactions is None:
reactions = chain.reactions
@ -92,7 +102,9 @@ def get_microxs_and_flux(
energy_filter = openmc.EnergyFilter.from_group_structure(energies)
else:
energy_filter = openmc.EnergyFilter(energies)
if isinstance(domains[0], openmc.Material):
if isinstance(domains, openmc.MeshBase):
domain_filter = openmc.MeshFilter(domains)
elif isinstance(domains[0], openmc.Material):
domain_filter = openmc.MaterialFilter(domains)
elif isinstance(domains[0], openmc.Cell):
domain_filter = openmc.CellFilter(domains)
@ -113,7 +125,7 @@ def get_microxs_and_flux(
model.tallies = openmc.Tallies([rr_tally, flux_tally])
# create temporary run
with tempfile.TemporaryDirectory() as temp_dir:
with TemporaryDirectory() as temp_dir:
if run_kwargs is None:
run_kwargs = {}
else:
@ -154,7 +166,7 @@ class MicroXS:
.. versionadded:: 0.13.1
.. versionchanged:: 0.13.4
.. versionchanged:: 0.14.0
Class was heavily refactored and no longer subclasses pandas.DataFrame.
Parameters
@ -187,9 +199,125 @@ class MicroXS:
self.data = data
self.nuclides = nuclides
self.reactions = reactions
self._index_nuc = {nuc: i for i, nuc in enumerate(nuclides)}
self._index_rx = {rx: i for i, rx in enumerate(reactions)}
# TODO: Add a classmethod for generating MicroXS directly from cross section
# data using a known flux spectrum
@classmethod
def from_multigroup_flux(
cls,
energies: Union[Sequence[float], str],
multigroup_flux: Sequence[float],
chain_file: Optional[PathLike] = None,
temperature: float = 293.6,
nuclides: Optional[Sequence[str]] = None,
reactions: Optional[Sequence[str]] = None,
**init_kwargs: dict,
) -> MicroXS:
"""Generated microscopic cross sections from a known flux.
The size of the MicroXS matrix depends on the chain file and cross
sections available. MicroXS entry will be 0 if the nuclide cross section
is not found.
.. versionadded:: 0.14.1
Parameters
----------
energies : iterable of float or str
Energy group boundaries in [eV] or the name of the group structure
multi_group_flux : iterable of float
Energy-dependent multigroup flux values
chain_file : str, optional
Path to the depletion chain XML file that will be used in depletion
simulation. Defaults to ``openmc.config['chain_file']``.
temperature : int, optional
Temperature for cross section evaluation in [K].
nuclides : list of str, optional
Nuclides to get cross sections for. If not specified, all burnable
nuclides from the depletion chain file are used.
reactions : list of str, optional
Reactions to get cross sections for. If not specified, all neutron
reactions listed in the depletion chain file are used.
**init_kwargs : dict
Keyword arguments passed to :func:`openmc.lib.init`
Returns
-------
MicroXS
"""
check_type("temperature", temperature, (int, float))
# if energy is string then use group structure of that name
if isinstance(energies, str):
energies = GROUP_STRUCTURES[energies]
else:
# if user inputs energies check they are ascending (low to high) as
# some depletion codes use high energy to low energy.
if not np.all(np.diff(energies) > 0):
raise ValueError('Energy group boundaries must be in ascending order')
# check dimension consistency
if len(multigroup_flux) != len(energies) - 1:
raise ValueError('Length of flux array should be len(energies)-1')
chain_file_path = _resolve_chain_file_path(chain_file)
chain = Chain.from_xml(chain_file_path)
cross_sections = _find_cross_sections(model=None)
nuclides_with_data = _get_nuclides_with_data(cross_sections)
# If no nuclides were specified, default to all nuclides from the chain
if not nuclides:
nuclides = chain.nuclides
nuclides = [nuc.name for nuc in nuclides]
# Get reaction MT values. If no reactions specified, default to the
# reactions available in the chain file
if reactions is None:
reactions = chain.reactions
mts = [REACTION_MT[name] for name in reactions]
# Normalize multigroup flux
multigroup_flux = np.asarray(multigroup_flux)
multigroup_flux /= multigroup_flux.sum()
# Create 2D array for microscopic cross sections
microxs_arr = np.zeros((len(nuclides), len(mts)))
with TemporaryDirectory() as tmpdir:
# Create a material with all nuclides
mat_all_nucs = openmc.Material()
for nuc in nuclides:
if nuc in nuclides_with_data:
mat_all_nucs.add_nuclide(nuc, 1.0)
mat_all_nucs.set_density("atom/b-cm", 1.0)
# Create simple model containing the above material
surf1 = openmc.Sphere(boundary_type="vacuum")
surf1_cell = openmc.Cell(fill=mat_all_nucs, region=-surf1)
model = openmc.Model()
model.geometry = openmc.Geometry([surf1_cell])
model.settings = openmc.Settings(
particles=1, batches=1, output={'summary': False})
with change_directory(tmpdir):
# Export model within temporary directory
model.export_to_model_xml()
with openmc.lib.run_in_memory(**init_kwargs):
# For each nuclide and reaction, compute the flux-averaged
# cross section
for nuc_index, nuc in enumerate(nuclides):
if nuc not in nuclides_with_data:
continue
lib_nuc = openmc.lib.nuclides[nuc]
for mt_index, mt in enumerate(mts):
xs = lib_nuc.collapse_rate(
mt, temperature, energies, multigroup_flux
)
microxs_arr[nuc_index, mt_index] = xs
return cls(microxs_arr, nuclides, reactions)
@classmethod
def from_csv(cls, csv_file, **kwargs):
@ -222,8 +350,8 @@ class MicroXS:
def __getitem__(self, index):
nuc, rx = index
i_nuc = self.nuclides.index(nuc)
i_rx = self.reactions.index(rx)
i_nuc = self._index_nuc[nuc]
i_rx = self._index_rx[rx]
return self.data[i_nuc, i_rx]
def to_csv(self, *args, **kwargs):

2
openmc/deplete/nuclide.py
View file

@ -8,8 +8,8 @@ from collections.abc import Mapping
from collections import namedtuple, defaultdict
from warnings import warn
from numbers import Real
import lxml.etree as ET
import lxml.etree as ET
import numpy as np
from openmc.checkvalue import check_type

17
openmc/deplete/openmc_operator.py
View file

@ -12,6 +12,7 @@ from typing import List, Tuple, Dict
import numpy as np
import openmc
from openmc.checkvalue import check_value
from openmc.exceptions import DataError
from openmc.mpi import comm
from .abc import TransportOperator, OperatorResult
@ -45,7 +46,6 @@ class OpenMCOperator(TransportOperator):
in the previous results.
diff_burnable_mats : bool, optional
Whether to differentiate burnable materials with multiple instances.
Volumes are divided equally from the original material volume.
fission_q : dict, optional
Dictionary of nuclides and their fission Q values [eV].
helper_kwargs : dict
@ -58,6 +58,14 @@ class OpenMCOperator(TransportOperator):
if ``reduce_chain`` evaluates to true. The default value of
``None`` implies no limit on the depth.
diff_volume_method : str
Specifies how the volumes of the new materials should be found. Default
is to 'divide equally' which divides the original material volume
equally between the new materials, 'match cell' sets the volume of the
material to volume of the cell they fill.
.. versionadded:: 0.14.0
Attributes
----------
materials : openmc.Materials
@ -97,6 +105,7 @@ class OpenMCOperator(TransportOperator):
chain_file=None,
prev_results=None,
diff_burnable_mats=False,
diff_volume_method='divide equally',
fission_q=None,
helper_kwargs=None,
reduce_chain=False,
@ -116,6 +125,10 @@ class OpenMCOperator(TransportOperator):
self.materials = materials
self.cross_sections = cross_sections
check_value('diff volume method', diff_volume_method,
{'divide equally', 'match cell'})
self.diff_volume_method = diff_volume_method
# Reduce the chain to only those nuclides present
if reduce_chain:
init_nuclides = set()
@ -200,7 +213,7 @@ class OpenMCOperator(TransportOperator):
msg = (f"Nuclilde {nuclide} in material {mat.id} is not "
"present in the depletion chain and has no cross "
"section data.")
raise warn(msg)
warn(msg)
if mat.depletable:
burnable_mats.add(str(mat.id))
if mat.volume is None:

2
openmc/deplete/pool.py
View file

@ -66,7 +66,7 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
transfer_rates : openmc.deplete.TransferRates, Optional
Object to perform continuous reprocessing.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
matrix_args: Any, optional
Additional arguments passed to matrix_func

8
openmc/deplete/results.py
View file

@ -104,6 +104,8 @@ class Results(list):
) -> Tuple[np.ndarray, typing.Union[np.ndarray, List[dict]]]:
"""Get activity of material over time.
.. versionadded:: 0.14.0
Parameters
----------
mat : openmc.Material, str
@ -220,6 +222,8 @@ class Results(list):
) -> Tuple[np.ndarray, typing.Union[np.ndarray, List[dict]]]:
"""Get decay heat of material over time.
.. versionadded:: 0.14.0
Parameters
----------
mat : openmc.Material, str
@ -273,7 +277,7 @@ class Results(list):
) -> Tuple[np.ndarray, np.ndarray]:
"""Get mass of nuclides over time from a single material
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -314,7 +318,7 @@ class Results(list):
# Divide by volume to get density
mass /= self[0].volume[mat_id]
elif mass_units == "kg":
mass *= 1e3
mass /= 1e3
return times, mass

2
openmc/deplete/stepresult.py
View file

@ -548,7 +548,7 @@ class StepResult:
path : PathLike
Path to file to write. Defaults to 'depletion_results.h5'.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
"""
# Get indexing terms
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()

2
openmc/deplete/transfer_rates.py
View file

@ -16,7 +16,7 @@ class TransferRates:
An instance of this class can be passed directly to an instance of one of
the :class:`openmc.deplete.Integrator` classes.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------

1
openmc/element.py
View file

@ -1,4 +1,5 @@
import re
import lxml.etree as ET
import openmc.checkvalue as cv

39
openmc/filter.py
View file

@ -3,9 +3,9 @@ from collections.abc import Iterable
import hashlib
from itertools import product
from numbers import Real, Integral
import lxml.etree as ET
import warnings
import lxml.etree as ET
import numpy as np
import pandas as pd
@ -22,7 +22,7 @@ from ._xml import get_text
_FILTER_TYPES = (
'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre',
'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
'collision', 'time'
)
@ -451,7 +451,7 @@ class UniverseFilter(WithIDFilter):
Parameters
----------
bins : openmc.UniverseBase, int, or iterable thereof
The Universes to tally. Either openmc.UniverseBase objects or their
The Universes to tally. Either :class:`openmc.UniverseBase` objects or their
Integral ID numbers can be used.
filter_id : int
Unique identifier for the filter
@ -475,7 +475,31 @@ class MaterialFilter(WithIDFilter):
Parameters
----------
bins : openmc.Material, Integral, or iterable thereof
The Materials to tally. Either openmc.Material objects or their
The material(s) to tally. Either :class:`openmc.Material` objects or their
Integral ID numbers can be used.
filter_id : int
Unique identifier for the filter
Attributes
----------
bins : Iterable of Integral
openmc.Material IDs.
id : int
Unique identifier for the filter
num_bins : Integral
The number of filter bins
"""
expected_type = Material
class MaterialFromFilter(WithIDFilter):
"""Bins tally event locations based on the Material they occurred in.
Parameters
----------
bins : openmc.Material, Integral, or iterable thereof
The material(s) to tally. Either :class:`openmc.Material` objects or their
Integral ID numbers can be used.
filter_id : int
Unique identifier for the filter
@ -499,7 +523,7 @@ class CellFilter(WithIDFilter):
Parameters
----------
bins : openmc.Cell, int, or iterable thereof
The cells to tally. Either openmc.Cell objects or their ID numbers can
The cells to tally. Either :class:`openmc.Cell` objects or their ID numbers can
be used.
filter_id : int
Unique identifier for the filter
@ -1222,8 +1246,8 @@ class RealFilter(Filter):
def get_bin_index(self, filter_bin):
i = np.where(self.bins[:, 1] == filter_bin[1])[0]
if len(i) == 0:
msg = (f'Unable to get the bin index for Filter since '
'"{filter_bin}" is not one of the bins')
msg = ('Unable to get the bin index for Filter since '
f'"{filter_bin}" is not one of the bins')
raise ValueError(msg)
else:
return i[0]
@ -1400,6 +1424,7 @@ class EnergyFilter(RealFilter):
"""
cv.check_value('group_structure', group_structure, openmc.mgxs.GROUP_STRUCTURES.keys())
return cls(openmc.mgxs.GROUP_STRUCTURES[group_structure.upper()])

1
openmc/filter_expansion.py
View file

@ -1,4 +1,5 @@
from numbers import Integral, Real
import lxml.etree as ET
import openmc.checkvalue as cv

24
openmc/geometry.py
View file

@ -293,7 +293,8 @@ class Geometry:
if isinstance(materials, (str, os.PathLike)):
materials = openmc.Materials.from_xml(materials)
tree = ET.parse(path)
parser = ET.XMLParser(huge_tree=True)
tree = ET.parse(path, parser=parser)
root = tree.getroot()
return cls.from_xml_element(root, materials)
@ -391,6 +392,20 @@ class Geometry:
universes.update(self.root_universe.get_all_universes())
return universes
def get_all_nuclides(self) -> typing.List[str]:
"""Return all nuclides within the geometry.
Returns
-------
list
Sorted list of all nuclides in materials appearing in the geometry
"""
all_nuclides = set()
for material in self.get_all_materials().values():
all_nuclides |= set(material.get_nuclides())
return sorted(all_nuclides)
def get_all_materials(self) -> typing.Dict[int, openmc.Material]:
"""Return all materials within the geometry.
@ -738,7 +753,7 @@ class Geometry:
def plot(self, *args, **kwargs):
"""Display a slice plot of the geometry.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -764,8 +779,9 @@ class Geometry:
Assigns colors to specific materials or cells. Keys are instances of
:class:`Cell` or :class:`Material` and values are RGB 3-tuples, RGBA
4-tuples, or strings indicating SVG color names. Red, green, blue,
and alpha should all be floats in the range [0.0, 1.0], for example:
.. code-block:: python
and alpha should all be floats in the range [0.0, 1.0], for
example::
# Make water blue
water = openmc.Cell(fill=h2o)
universe.plot(..., colors={water: (0., 0., 1.))

2
openmc/lattice.py
View file

@ -4,8 +4,8 @@ from copy import deepcopy
from math import sqrt, floor
from numbers import Real
import types
import lxml.etree as ET
import lxml.etree as ET
import numpy as np
import openmc

9
openmc/lib/core.py
View file

@ -174,7 +174,7 @@ def export_properties(filename=None, output=True):
def export_weight_windows(filename="weight_windows.h5", output=True):
"""Export weight windows.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -198,7 +198,7 @@ def export_weight_windows(filename="weight_windows.h5", output=True):
def import_weight_windows(filename='weight_windows.h5', output=True):
"""Import weight windows.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -611,7 +611,7 @@ class _DLLGlobal:
class _FortranObject:
def __repr__(self):
return "{}[{}]".format(type(self).__name__, self._index)
return "<{}(index={})>".format(type(self).__name__, self._index)
class _FortranObjectWithID(_FortranObject):
@ -622,6 +622,9 @@ class _FortranObjectWithID(_FortranObject):
# OutOfBoundsError will be raised here by virtue of referencing self.id
self.id
def __repr__(self):
return "<{}(id={})>".format(type(self).__name__, self.id)
@contextmanager
def quiet_dll(output=True):

7
openmc/lib/filter.py
View file

@ -20,7 +20,7 @@ __all__ = [
'Filter', 'AzimuthalFilter', 'CellFilter', 'CellbornFilter', 'CellfromFilter',
'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter',
'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
'MaterialFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter',
'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter',
'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter',
'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
]
@ -342,6 +342,10 @@ class MaterialFilter(Filter):
_dll.openmc_material_filter_set_bins(self._index, n, bins)
class MaterialFromFilter(Filter):
filter_type = 'materialfrom'
class MeshFilter(Filter):
filter_type = 'mesh'
@ -501,6 +505,7 @@ _FILTER_TYPE_MAP = {
'energyfunction': EnergyFunctionFilter,
'legendre': LegendreFilter,
'material': MaterialFilter,
'materialfrom': MaterialFromFilter,
'mesh': MeshFilter,
'meshsurface': MeshSurfaceFilter,
'mu': MuFilter,

20
openmc/lib/material.py
View file

@ -1,5 +1,5 @@
from collections.abc import Mapping
from ctypes import c_int, c_int32, c_double, c_char_p, POINTER, c_size_t
from ctypes import c_bool, c_int, c_int32, c_double, c_char_p, POINTER, c_size_t
from weakref import WeakValueDictionary
import numpy as np
@ -60,6 +60,12 @@ _dll.openmc_material_set_name.errcheck = _error_handler
_dll.openmc_material_set_volume.argtypes = [c_int32, c_double]
_dll.openmc_material_set_volume.restype = c_int
_dll.openmc_material_set_volume.errcheck = _error_handler
_dll.openmc_material_get_depletable.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_material_get_depletable.restype = c_int
_dll.openmc_material_get_depletable.errcheck = _error_handler
_dll.openmc_material_set_depletable.argtypes = [c_int32, c_bool]
_dll.openmc_material_set_depletable.restype = c_int
_dll.openmc_material_set_depletable.errcheck = _error_handler
_dll.n_materials.argtypes = []
_dll.n_materials.restype = c_size_t
@ -89,6 +95,8 @@ class Material(_FortranObjectWithID):
List of nuclides in the material
densities : numpy.ndarray
Array of densities in atom/b-cm
depletable : bool
Whether this material is marked as depletable
name : str
Name of the material
temperature : float
@ -169,6 +177,16 @@ class Material(_FortranObjectWithID):
def volume(self, volume):
_dll.openmc_material_set_volume(self._index, volume)
@property
def depletable(self):
depletable = c_bool()
_dll.openmc_material_get_depletable(self._index, depletable)
return depletable.value
@depletable.setter
def depletable(self, depletable):
_dll.openmc_material_set_depletable(self._index, depletable)
@property
def nuclides(self):
return self._get_densities()[0]

2
openmc/lib/math.py
View file

@ -103,7 +103,7 @@ def evaluate_legendre(data, x):
"""
data_arr = np.array(data, dtype=np.float64)
return _dll.evaluate_legendre(len(data),
return _dll.evaluate_legendre(len(data)-1,
data_arr.ctypes.data_as(POINTER(c_double)), x)

130
openmc/lib/mesh.py
View file

@ -1,6 +1,8 @@
from collections.abc import Mapping
from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
create_string_buffer)
from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, Structure,
create_string_buffer, c_uint64, c_size_t)
from random import getrandbits
from typing import Optional, List, Tuple, Sequence
from weakref import WeakValueDictionary
import numpy as np
@ -10,9 +12,19 @@ from ..exceptions import AllocationError, InvalidIDError
from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler
from .material import Material
from .plot import _Position
__all__ = ['RegularMesh', 'RectilinearMesh', 'CylindricalMesh', 'SphericalMesh', 'UnstructuredMesh', 'meshes']
class _MaterialVolume(Structure):
_fields_ = [
("material", c_int32),
("volume", c_double)
]
# Mesh functions
_dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32),
POINTER(c_int32)]
@ -24,6 +36,19 @@ _dll.openmc_mesh_get_id.errcheck = _error_handler
_dll.openmc_mesh_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_mesh_set_id.restype = c_int
_dll.openmc_mesh_set_id.errcheck = _error_handler
_dll.openmc_mesh_get_n_elements.argtypes = [c_int32, POINTER(c_size_t)]
_dll.openmc_mesh_get_n_elements.restype = c_int
_dll.openmc_mesh_get_n_elements.errcheck = _error_handler
_dll.openmc_mesh_material_volumes.argtypes = [
c_int32, c_int, c_int, c_int, POINTER(_MaterialVolume),
POINTER(c_int), POINTER(c_uint64)]
_dll.openmc_mesh_material_volumes.restype = c_int
_dll.openmc_mesh_material_volumes.errcheck = _error_handler
_dll.openmc_mesh_get_plot_bins.argtypes = [
c_int32, _Position, _Position, c_int, POINTER(c_int), POINTER(c_int32)
]
_dll.openmc_mesh_get_plot_bins.restype = c_int
_dll.openmc_mesh_get_plot_bins.errcheck = _error_handler
_dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_mesh_index.restype = c_int
_dll.openmc_get_mesh_index.errcheck = _error_handler
@ -123,6 +148,107 @@ class Mesh(_FortranObjectWithID):
def id(self, mesh_id):
_dll.openmc_mesh_set_id(self._index, mesh_id)
@property
def n_elements(self):
n = c_size_t()
_dll.openmc_mesh_get_n_elements(self._index, n)
return n.value
def material_volumes(
self,
n_samples: int = 10_000,
prn_seed: Optional[int] = None
) -> List[List[Tuple[Material, float]]]:
"""Determine volume of materials in each mesh element
.. versionadded:: 0.14.1
Parameters
----------
n_samples : int
Number of samples in each mesh element
prn_seed : int
Pseudorandom number generator (PRNG) seed; if None, one will be
generated randomly.
Returns
-------
List of tuple of (material, volume) for each mesh element. Void volume
is represented by having a value of None in the first element of a
tuple.
"""
if n_samples <= 0:
raise ValueError("Number of samples must be positive")
if prn_seed is None:
prn_seed = getrandbits(63)
prn_seed = c_uint64(prn_seed)
# Preallocate space for MaterialVolume results
size = 16
result = (_MaterialVolume * size)()
hits = c_int() # Number of materials hit in a given element
volumes = []
for i_element in range(self.n_elements):
while True:
try:
_dll.openmc_mesh_material_volumes(
self._index, n_samples, i_element, size, result, hits, prn_seed)
except AllocationError:
# Increase size of result array and try again
size *= 2
result = (_MaterialVolume * size)()
else:
# If no error, break out of loop
break
volumes.append([
(Material(index=r.material), r.volume)
for r in result[:hits.value]
])
return volumes
def get_plot_bins(
self,
origin: Sequence[float],
width: Sequence[float],
basis: str,
pixels: Sequence[int]
) -> np.ndarray:
"""Get mesh bin indices for a rasterized plot.
.. versionadded:: 0.14.1
Parameters
----------
origin : iterable of float
Origin of the plotting view. Should have length 3.
width : iterable of float
Width of the plotting view. Should have length 2.
basis : {'xy', 'xz', 'yz'}
Plotting basis.
pixels : iterable of int
Number of pixels in each direction. Should have length 2.
Returns
-------
2D numpy array with mesh bin indices corresponding to each pixel within
the plotting view.
"""
origin = _Position(*origin)
width = _Position(*width)
basis = {'xy': 1, 'xz': 2, 'yz': 3}[basis]
pixel_array = (c_int*2)(*pixels)
img_data = np.zeros((pixels[1], pixels[0]), dtype=np.dtype('int32'))
_dll.openmc_mesh_get_plot_bins(
self._index, origin, width, basis, pixel_array,
img_data.ctypes.data_as(POINTER(c_int32))
)
return img_data
class RegularMesh(Mesh):
"""RegularMesh stored internally.

2
openmc/lib/nuclide.py
View file

@ -91,7 +91,7 @@ class Nuclide(_FortranObject):
energy : iterable of float
Energy group boundaries in [eV]
flux : iterable of float
Flux in each energt group (not normalized per eV)
Flux in each energy group (not normalized per eV)
Returns
-------

2
openmc/lib/tally.py
View file

@ -183,7 +183,7 @@ class Tally(_FortranObjectWithID):
multiply_density : bool
Whether reaction rates should be multiplied by atom density
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
nuclides : list of str
List of nuclides to score results for
num_realizations : int

4
openmc/lib/weight_windows.py
View file

@ -109,7 +109,7 @@ class WeightWindows(_FortranObjectWithID):
OpenMC library. To obtain a view of a weight windows object with a given ID,
use the :data:`openmc.lib.weight_windows` mapping.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -392,4 +392,4 @@ class _WeightWindowsMapping(Mapping):
def __delitem__(self):
raise NotImplementedError("WeightWindows object remove not implemented")
weight_windows = _WeightWindowsMapping()
weight_windows = _WeightWindowsMapping()

78
openmc/material.py
View file

@ -8,8 +8,8 @@ import re
import typing # imported separately as py3.8 requires typing.Iterable
import warnings
from typing import Optional, List, Union, Dict
import lxml.etree as ET
import lxml.etree as ET
import numpy as np
import h5py
@ -19,7 +19,7 @@ import openmc.checkvalue as cv
from ._xml import clean_indentation, reorder_attributes
from .mixin import IDManagerMixin
from openmc.checkvalue import PathLike
from openmc.stats import Univariate
from openmc.stats import Univariate, Discrete, Mixture
# Units for density supported by OpenMC
@ -93,13 +93,6 @@ class Material(IDManagerMixin):
fissionable_mass : float
Mass of fissionable nuclides in the material in [g]. Requires that the
:attr:`volume` attribute is set.
decay_photon_energy : openmc.stats.Univariate or None
Energy distribution of photons emitted from decay of unstable nuclides
within the material, or None if no photon source exists. The integral of
this distribution is the total intensity of the photon source in
[decay/sec].
.. versionadded:: 0.13.2
ncrystal_cfg : str
NCrystal configuration string
@ -282,15 +275,67 @@ class Material(IDManagerMixin):
@property
def decay_photon_energy(self) -> Optional[Univariate]:
atoms = self.get_nuclide_atoms()
warnings.warn(
"The 'decay_photon_energy' property has been replaced by the "
"get_decay_photon_energy() method and will be removed in a future "
"version.", FutureWarning)
return self.get_decay_photon_energy(0.0)
def get_decay_photon_energy(
self,
clip_tolerance: float = 1e-6,
units: str = 'Bq',
volume: Optional[float] = None
) -> Optional[Univariate]:
r"""Return energy distribution of decay photons from unstable nuclides.
.. versionadded:: 0.14.0
Parameters
----------
clip_tolerance : float
Maximum fraction of :math:`\sum_i x_i p_i` for discrete
distributions that will be discarded.
units : {'Bq', 'Bq/g', 'Bq/cm3'}
Specifies the units on the integral of the distribution.
volume : float, optional
Volume of the material. If not passed, defaults to using the
:attr:`Material.volume` attribute.
Returns
-------
Decay photon energy distribution. The integral of this distribution is
the total intensity of the photon source in the requested units.
"""
cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/cm3'})
if units == 'Bq':
multiplier = volume if volume is not None else self.volume
if multiplier is None:
raise ValueError("volume must be specified if units='Bq'")
elif units == 'Bq/cm3':
multiplier = 1
elif units == 'Bq/g':
multiplier = 1.0 / self.get_mass_density()
dists = []
probs = []
for nuc, num_atoms in atoms.items():
for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items():
source_per_atom = openmc.data.decay_photon_energy(nuc)
if source_per_atom is not None:
dists.append(source_per_atom)
probs.append(num_atoms)
return openmc.data.combine_distributions(dists, probs) if dists else None
probs.append(1e24 * atoms_per_bcm * multiplier)
# If no photon sources, exit early
if not dists:
return None
# Get combined distribution, clip low-intensity values in discrete spectra
combined = openmc.data.combine_distributions(dists, probs)
if isinstance(combined, (Discrete, Mixture)):
combined.clip(clip_tolerance, inplace=True)
return combined
@classmethod
def from_hdf5(cls, group: h5py.Group) -> Material:
@ -1668,7 +1713,7 @@ class Materials(cv.CheckedList):
self._write_xml(fh, nuclides_to_ignore=nuclides_to_ignore)
@classmethod
def from_xml_element(cls, elem) -> Material:
def from_xml_element(cls, elem) -> Materials:
"""Generate materials collection from XML file
Parameters
@ -1695,7 +1740,7 @@ class Materials(cv.CheckedList):
return materials
@classmethod
def from_xml(cls, path: PathLike = 'materials.xml') -> Material:
def from_xml(cls, path: PathLike = 'materials.xml') -> Materials:
"""Generate materials collection from XML file
Parameters
@ -1709,7 +1754,8 @@ class Materials(cv.CheckedList):
Materials collection
"""
tree = ET.parse(path)
parser = ET.XMLParser(huge_tree=True)
tree = ET.parse(path, parser=parser)
root = tree.getroot()
return cls.from_xml_element(root)

297
openmc/mesh.py
View file

@ -1,8 +1,9 @@
from __future__ import annotations
import typing
import warnings
from abc import ABC, abstractmethod
from collections.abc import Iterable
from math import pi
from math import pi, sqrt, atan2
from numbers import Integral, Real
from pathlib import Path
from typing import Optional, Sequence, Tuple
@ -35,6 +36,9 @@ class MeshBase(IDManagerMixin, ABC):
Unique identifier for the mesh
name : str
Name of the mesh
bounding_box : openmc.BoundingBox
Axis-aligned bounding box of the mesh as defined by the upper-right and
lower-left coordinates.
"""
@ -58,6 +62,10 @@ class MeshBase(IDManagerMixin, ABC):
else:
self._name = ''
@property
def bounding_box(self) -> openmc.BoundingBox:
return openmc.BoundingBox(self.lower_left, self.upper_right)
def __repr__(self):
string = type(self).__name__ + '\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@ -170,24 +178,27 @@ class StructuredMesh(MeshBase):
@property
def vertices(self):
"""Return coordinates of mesh vertices.
"""Return coordinates of mesh vertices in Cartesian coordinates. Also
see :meth:`CylindricalMesh.vertices_cylindrical` and
:meth:`SphericalMesh.vertices_spherical` for coordinates in other coordinate
systems.
Returns
-------
vertices : numpy.ndarray
Returns a numpy.ndarray representing the coordinates of the mesh
vertices with a shape equal to (ndim, dim1 + 1, ..., dimn + 1). Can be
unpacked along the first dimension with xx, yy, zz = mesh.vertices.
vertices with a shape equal to (dim1 + 1, ..., dimn + 1, ndim). X, Y, Z values
can be unpacked with xx, yy, zz = np.rollaxis(mesh.vertices, -1).
"""
return self._generate_vertices(*self._grids)
@staticmethod
def _generate_vertices(i_grid, j_grid, k_grid):
"""Returns an array with shape (3, i_grid.size+1, j_grid.size+1, k_grid.size+1)
"""Returns an array with shape (i_grid.size, j_grid.size, k_grid.size, 3)
containing the corner vertices of mesh elements.
"""
return np.stack(np.meshgrid(i_grid, j_grid, k_grid, indexing='ij'), axis=0)
return np.stack(np.meshgrid(i_grid, j_grid, k_grid, indexing='ij'), axis=-1)
@staticmethod
def _generate_edge_midpoints(grids):
@ -203,7 +214,7 @@ class StructuredMesh(MeshBase):
midpoint_grids : list of numpy.ndarray
The edge midpoints for the i, j, and k midpoints of each element in
i, j, k ordering. The shapes of the resulting grids are
[(3, ni-1, nj, nk), (3, ni, nj-1, nk), (3, ni, nj, nk-1)]
[(ni-1, nj, nk, 3), (ni, nj-1, nk, 3), (ni, nj, nk-1, 3)]
"""
# generate a set of edge midpoints for each dimension
midpoint_grids = []
@ -213,7 +224,7 @@ class StructuredMesh(MeshBase):
# each grid is comprised of the mid points for one dimension and the
# corner vertices of the other two
for dims in ((0, 1, 2), (1, 0, 2), (2, 0, 1)):
# compute the midpoints along the first dimension
# compute the midpoints along the last dimension
midpoints = grids[dims[0]][:-1] + 0.5 * np.diff(grids[dims[0]])
coords = (midpoints, grids[dims[1]], grids[dims[2]])
@ -248,15 +259,16 @@ class StructuredMesh(MeshBase):
-------
centroids : numpy.ndarray
Returns a numpy.ndarray representing the mesh element centroid
coordinates with a shape equal to (ndim, dim1, ..., dimn). Can be
unpacked along the first dimension with xx, yy, zz = mesh.centroids.
coordinates with a shape equal to (dim1, ..., dimn, ndim). X,
Y, Z values can be unpacked with xx, yy, zz =
np.rollaxis(mesh.centroids, -1).
"""
ndim = self.n_dimension
vertices = self.vertices
s0 = (slice(None),) + (slice(0, -1),)*ndim
s1 = (slice(None),) + (slice(1, None),)*ndim
# this line ensures that the vertices aren't adjusted by the origin or
# converted to the Cartesian system for cylindrical and spherical meshes
vertices = StructuredMesh.vertices.fget(self)
s0 = (slice(0, -1),)*ndim + (slice(None),)
s1 = (slice(1, None),)*ndim + (slice(None),)
return (vertices[s0] + vertices[s1]) / 2
@property
@ -346,10 +358,8 @@ class StructuredMesh(MeshBase):
import vtk
from vtk.util import numpy_support as nps
vertices = self.cartesian_vertices.T.reshape(-1, 3)
vtkPts = vtk.vtkPoints()
vtkPts.SetData(nps.numpy_to_vtk(vertices, deep=True))
vtkPts.SetData(nps.numpy_to_vtk(np.swapaxes(self.vertices, 0, 2).reshape(-1, 3), deep=True))
vtk_grid = vtk.vtkStructuredGrid()
vtk_grid.SetPoints(vtkPts)
vtk_grid.SetDimensions(*[dim + 1 for dim in self.dimension])
@ -368,7 +378,7 @@ class StructuredMesh(MeshBase):
import vtk
from vtk.util import numpy_support as nps
corner_vertices = self.cartesian_vertices.T.reshape(-1, 3)
corner_vertices = np.swapaxes(self.vertices, 0, 2).reshape(-1, 3)
vtkPts = vtk.vtkPoints()
vtk_grid = vtk.vtkUnstructuredGrid()
@ -410,7 +420,7 @@ class StructuredMesh(MeshBase):
# list of point IDs
midpoint_vertices = self.midpoint_vertices
for edge_grid in midpoint_vertices:
for pnt in edge_grid.T.reshape(-1, 3):
for pnt in np.swapaxes(edge_grid, 0, 2).reshape(-1, 3):
point_ids.append(_insert_point(pnt))
# determine how many elements in each dimension
@ -443,7 +453,7 @@ class StructuredMesh(MeshBase):
# initial offset for corner vertices and midpoint dimension
flat_idx = corner_vertices.shape[0] + sum(n_midpoint_vertices[:dim])
# generate a flat index into the table of point IDs
midpoint_shape = midpoint_vertices[dim].shape[1:]
midpoint_shape = midpoint_vertices[dim].shape[:-1]
flat_idx += np.ravel_multi_index((i+di, j+dj, k+dk),
midpoint_shape,
order='F')
@ -505,9 +515,9 @@ class RegularMesh(StructuredMesh):
upper_right : Iterable of float
The upper-right corner of the structured mesh. If only two coordinate
are given, it is assumed that the mesh is an x-y mesh.
bounding_box: openmc.BoundingBox
Axis-aligned bounding box of the cell defined by the upper-right and lower-
left coordinates
bounding_box : openmc.BoundingBox
Axis-aligned bounding box of the mesh as defined by the upper-right and
lower-left coordinates.
width : Iterable of float
The width of mesh cells in each direction.
indices : Iterable of tuple
@ -599,12 +609,6 @@ class RegularMesh(StructuredMesh):
self._upper_right = None
warnings.warn("Unsetting upper_right attribute.")
@property
def cartesian_vertices(self):
"""Returns vertices in cartesian coordinates. Identical to ``vertices`` for RegularMesh and RectilinearMesh
"""
return self.vertices
@property
def volumes(self):
"""Return Volumes for every mesh cell
@ -664,12 +668,6 @@ class RegularMesh(StructuredMesh):
x1, = self.upper_right
return (np.linspace(x0, x1, nx + 1),)
@property
def bounding_box(self):
return openmc.BoundingBox(
np.array(self.lower_left), np.array(self.upper_right)
)
def __repr__(self):
string = super().__repr__()
string += '{0: <16}{1}{2}\n'.format('\tDimensions', '=\t', self.n_dimension)
@ -1009,6 +1007,9 @@ class RectilinearMesh(StructuredMesh):
indices : Iterable of tuple
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
(2, 1, 1), ...]
bounding_box : openmc.BoundingBox
Axis-aligned bounding box of the mesh as defined by the upper-right and
lower-left coordinates.
"""
@ -1036,7 +1037,7 @@ class RectilinearMesh(StructuredMesh):
@x_grid.setter
def x_grid(self, grid):
cv.check_type('mesh x_grid', grid, Iterable, Real)
self._x_grid = np.asarray(grid)
self._x_grid = np.asarray(grid, dtype=float)
@property
def y_grid(self):
@ -1045,7 +1046,7 @@ class RectilinearMesh(StructuredMesh):
@y_grid.setter
def y_grid(self, grid):
cv.check_type('mesh y_grid', grid, Iterable, Real)
self._y_grid = np.asarray(grid)
self._y_grid = np.asarray(grid, dtype=float)
@property
def z_grid(self):
@ -1054,17 +1055,19 @@ class RectilinearMesh(StructuredMesh):
@z_grid.setter
def z_grid(self, grid):
cv.check_type('mesh z_grid', grid, Iterable, Real)
self._z_grid = np.asarray(grid)
self._z_grid = np.asarray(grid, dtype=float)
@property
def _grids(self):
return (self.x_grid, self.y_grid, self.z_grid)
@property
def cartesian_vertices(self):
"""Returns vertices in cartesian coordiantes. Identical to ``vertices`` for RegularMesh and RectilinearMesh
"""
return self.vertices
def lower_left(self):
return np.array([self.x_grid[0], self.y_grid[0], self.z_grid[0]])
@property
def upper_right(self):
return np.array([self.x_grid[-1], self.y_grid[-1], self.z_grid[-1]])
@property
def volumes(self):
@ -1185,7 +1188,7 @@ class CylindricalMesh(StructuredMesh):
Parameters
----------
r_grid : numpy.ndarray
1-D array of mesh boundary points along the r-axis.
1-D array of mesh boundary points along the r-axis
Requirement is r >= 0.
z_grid : numpy.ndarray
1-D array of mesh boundary points along the z-axis relative to the
@ -1232,9 +1235,9 @@ class CylindricalMesh(StructuredMesh):
upper_right : Iterable of float
The upper-right corner of the structured mesh. If only two coordinate
are given, it is assumed that the mesh is an x-y mesh.
bounding_box: openmc.OpenMC
Axis-aligned cartesian bounding box of cell defined by upper-right and lower-
left coordinates
bounding_box : openmc.BoundingBox
Axis-aligned bounding box of the mesh as defined by the upper-right and
lower-left coordinates.
"""
@ -1249,9 +1252,9 @@ class CylindricalMesh(StructuredMesh):
):
super().__init__(mesh_id, name)
self._r_grid = r_grid
self._phi_grid = phi_grid
self._z_grid = z_grid
self.r_grid = r_grid
self.phi_grid = phi_grid
self.z_grid = z_grid
self.origin = origin
@property
@ -1281,7 +1284,7 @@ class CylindricalMesh(StructuredMesh):
@r_grid.setter
def r_grid(self, grid):
cv.check_type('mesh r_grid', grid, Iterable, Real)
self._r_grid = np.asarray(grid)
self._r_grid = np.asarray(grid, dtype=float)
@property
def phi_grid(self):
@ -1290,7 +1293,7 @@ class CylindricalMesh(StructuredMesh):
@phi_grid.setter
def phi_grid(self, grid):
cv.check_type('mesh phi_grid', grid, Iterable, Real)
self._phi_grid = np.asarray(grid)
self._phi_grid = np.asarray(grid, dtype=float)
@property
def z_grid(self):
@ -1299,7 +1302,7 @@ class CylindricalMesh(StructuredMesh):
@z_grid.setter
def z_grid(self, grid):
cv.check_type('mesh z_grid', grid, Iterable, Real)
self._z_grid = np.asarray(grid)
self._z_grid = np.asarray(grid, dtype=float)
@property
def _grids(self):
@ -1331,10 +1334,6 @@ class CylindricalMesh(StructuredMesh):
self.origin[2] + self.z_grid[-1]
))
@property
def bounding_box(self):
return openmc.BoundingBox(self.lower_left, self.upper_right)
def __repr__(self):
fmt = '{0: <16}{1}{2}\n'
string = super().__repr__()
@ -1357,6 +1356,66 @@ class CylindricalMesh(StructuredMesh):
string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
return string
def get_indices_at_coords(
self,
coords: Sequence[float]
) -> Tuple[int, int, int]:
"""Finds the index of the mesh voxel at the specified x,y,z coordinates.
Parameters
----------
coords : Sequence[float]
The x, y, z axis coordinates
Returns
-------
Tuple[int, int, int]
The r, phi, z indices
"""
r_value_from_origin = sqrt((coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)
if r_value_from_origin < self.r_grid[0] or r_value_from_origin > self.r_grid[-1]:
raise ValueError(
f'The specified x, y ({coords[0]}, {coords[1]}) combine to give an r value of '
f'{r_value_from_origin} from the origin of {self.origin}.which '
f'is outside the origin absolute r grid values {self.r_grid}.'
)
r_index = np.searchsorted(self.r_grid, r_value_from_origin) - 1
z_grid_values = np.array(self.z_grid) + self.origin[2]
if coords[2] < z_grid_values[0] or coords[2] > z_grid_values[-1]:
raise ValueError(
f'The specified z value ({coords[2]}) from the z origin of '
f'{self.origin[-1]} is outside of the absolute z grid range {z_grid_values}.'
)
z_index = np.argmax(z_grid_values > coords[2]) - 1
delta_x = coords[0] - self.origin[0]
delta_y = coords[1] - self.origin[1]
# atan2 returns values in -pi to +pi range
phi_value = atan2(delta_y, delta_x)
if delta_x < 0 and delta_y < 0:
# returned phi_value anticlockwise and negative
phi_value += 2 * pi
if delta_x > 0 and delta_y < 0:
# returned phi_value anticlockwise and negative
phi_value += 2 * pi
phi_grid_values = np.array(self.phi_grid)
if phi_value < phi_grid_values[0] or phi_value > phi_grid_values[-1]:
raise ValueError(
f'The phi value ({phi_value}) resulting from the specified x, y '
f'values is outside of the absolute phi grid range {phi_grid_values}.'
)
phi_index = np.argmax(phi_grid_values > phi_value) - 1
return (r_index, phi_index, z_index)
@classmethod
def from_hdf5(cls, group: h5py.Group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
@ -1439,6 +1498,10 @@ class CylindricalMesh(StructuredMesh):
num=dimension[2]+1
)
origin = (cached_bb.center[0], cached_bb.center[1], z_grid[0])
# make z-grid relative to the origin
z_grid -= origin[2]
mesh = cls(
r_grid=r_grid,
z_grid=z_grid,
@ -1523,19 +1586,37 @@ class CylindricalMesh(StructuredMesh):
return np.multiply.outer(np.outer(V_r, V_p), V_z)
@property
def cartesian_vertices(self):
return self._convert_to_cartesian(self.vertices, self.origin)
def vertices(self):
warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
return self._convert_to_cartesian(self.vertices_cylindrical, self.origin)
@property
def vertices_cylindrical(self):
"""Returns vertices of the mesh in cylindrical coordinates.
"""
return super().vertices
@property
def centroids(self):
warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
return self._convert_to_cartesian(self.centroids_cylindrical, self.origin)
@property
def centroids_cylindrical(self):
"""Returns centroids of the mesh in cylindrical coordinates.
"""
return super().centroids
@staticmethod
def _convert_to_cartesian(arr, origin: Sequence[float]):
"""Converts an array with xyz values in the first dimension (shape (3, ...))
"""Converts an array with r, phi, z values in the last dimension (shape (..., 3))
to Cartesian coordinates.
"""
x = arr[0, ...] * np.cos(arr[1, ...]) + origin[0]
y = arr[0, ...] * np.sin(arr[1, ...]) + origin[1]
arr[0, ...] = x
arr[1, ...] = y
arr[2, ...] += origin[2]
x = arr[..., 0] * np.cos(arr[..., 1]) + origin[0]
y = arr[..., 0] * np.sin(arr[..., 1]) + origin[1]
arr[..., 0] = x
arr[..., 1] = y
arr[..., 2] += origin[2]
return arr
@ -1594,8 +1675,8 @@ class SphericalMesh(StructuredMesh):
The upper-right corner of the structured mesh. If only two coordinate
are given, it is assumed that the mesh is an x-y mesh.
bounding_box : openmc.BoundingBox
Axis-aligned bounding box of the cell defined by the upper-right and lower-
left coordinates
Axis-aligned bounding box of the mesh as defined by the upper-right and
lower-left coordinates.
"""
@ -1610,9 +1691,9 @@ class SphericalMesh(StructuredMesh):
):
super().__init__(mesh_id, name)
self._r_grid = r_grid
self._theta_grid = theta_grid
self._phi_grid = phi_grid
self.r_grid = r_grid
self.theta_grid = theta_grid
self.phi_grid = phi_grid
self.origin = origin
@property
@ -1633,7 +1714,7 @@ class SphericalMesh(StructuredMesh):
def origin(self, coords):
cv.check_type('mesh origin', coords, Iterable, Real)
cv.check_length("mesh origin", coords, 3)
self._origin = np.asarray(coords)
self._origin = np.asarray(coords, dtype=float)
@property
def r_grid(self):
@ -1642,7 +1723,7 @@ class SphericalMesh(StructuredMesh):
@r_grid.setter
def r_grid(self, grid):
cv.check_type('mesh r_grid', grid, Iterable, Real)
self._r_grid = np.asarray(grid)
self._r_grid = np.asarray(grid, dtype=float)
@property
def theta_grid(self):
@ -1651,7 +1732,7 @@ class SphericalMesh(StructuredMesh):
@theta_grid.setter
def theta_grid(self, grid):
cv.check_type('mesh theta_grid', grid, Iterable, Real)
self._theta_grid = np.asarray(grid)
self._theta_grid = np.asarray(grid, dtype=float)
@property
def phi_grid(self):
@ -1660,7 +1741,7 @@ class SphericalMesh(StructuredMesh):
@phi_grid.setter
def phi_grid(self, grid):
cv.check_type('mesh phi_grid', grid, Iterable, Real)
self._phi_grid = np.asarray(grid)
self._phi_grid = np.asarray(grid, dtype=float)
@property
def _grids(self):
@ -1686,10 +1767,6 @@ class SphericalMesh(StructuredMesh):
r = self.r_grid[-1]
return np.array((self.origin[0] + r, self.origin[1] + r, self.origin[2] + r))
@property
def bounding_box(self):
return openmc.BoundingBox(self.lower_left, self.upper_right)
def __repr__(self):
fmt = '{0: <16}{1}{2}\n'
string = super().__repr__()
@ -1800,21 +1877,40 @@ class SphericalMesh(StructuredMesh):
return np.multiply.outer(np.outer(V_r, V_t), V_p)
@property
def cartesian_vertices(self):
return self._convert_to_cartesian(self.vertices, self.origin)
def vertices(self):
warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
return self._convert_to_cartesian(self.vertices_spherical, self.origin)
@property
def vertices_spherical(self):
"""Returns vertices of the mesh in cylindrical coordinates.
"""
return super().vertices
@property
def centroids(self):
warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
return self._convert_to_cartesian(self.centroids_spherical, self.origin)
@property
def centroids_spherical(self):
"""Returns centroids of the mesh in cylindrical coordinates.
"""
return super().centroids
@staticmethod
def _convert_to_cartesian(arr, origin: Sequence[float]):
"""Converts an array with xyz values in the first dimension (shape (3, ...))
"""Converts an array with r, theta, phi values in the last dimension (shape (..., 3))
to Cartesian coordinates.
"""
r_xy = arr[0, ...] * np.sin(arr[1, ...])
x = r_xy * np.cos(arr[2, ...])
y = r_xy * np.sin(arr[2, ...])
z = arr[0, ...] * np.cos(arr[1, ...])
arr[0, ...] = x + origin[0]
arr[1, ...] = y + origin[1]
arr[2, ...] = z + origin[2]
r_xy = arr[..., 0] * np.sin(arr[..., 1])
x = r_xy * np.cos(arr[..., 2])
y = r_xy * np.sin(arr[..., 2])
z = arr[..., 0] * np.cos(arr[..., 1])
arr[..., 0] = x + origin[0]
arr[..., 1] = y + origin[1]
arr[..., 2] = z + origin[2]
return arr
@ -1852,8 +1948,8 @@ class UnstructuredMesh(MeshBase):
library : {'moab', 'libmesh'}
Mesh library used for the unstructured mesh tally
output : bool
Indicates whether or not automatic tally output should
be generated for this mesh
Indicates whether or not automatic tally output should be generated for
this mesh
volumes : Iterable of float
Volumes of the unstructured mesh elements
centroids : numpy.ndarray
@ -1873,6 +1969,10 @@ class UnstructuredMesh(MeshBase):
.. versionadded:: 0.13.1
total_volume : float
Volume of the unstructured mesh in total
bounding_box : openmc.BoundingBox
Axis-aligned bounding box of the mesh as defined by the upper-right and
lower-left coordinates.
"""
_UNSUPPORTED_ELEM = -1
@ -2005,6 +2105,14 @@ class UnstructuredMesh(MeshBase):
self.length_multiplier)
return string
@property
def lower_left(self):
return self.vertices.min(axis=0)
@property
def upper_right(self):
return self.vertices.max(axis=0)
def centroid(self, bin: int):
"""Return the vertex averaged centroid of an element
@ -2090,7 +2198,6 @@ class UnstructuredMesh(MeshBase):
grid.SetPoints(vtk_pnts)
n_skipped = 0
elems = []
for elem_type, conn in zip(self.element_types, self.connectivity):
if elem_type == self._LINEAR_TET:
elem = vtk.vtkTetra()
@ -2104,15 +2211,13 @@ class UnstructuredMesh(MeshBase):
if c == -1:
break
elem.GetPointIds().SetId(i, c)
elems.append(elem)
grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())
if n_skipped > 0:
warnings.warn(f'{n_skipped} elements were not written because '
'they are not of type linear tet/hex')
for elem in elems:
grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds())
# check that datasets are the correct size
datasets_out = []
if datasets is not None:
@ -2225,7 +2330,7 @@ def _read_meshes(elem):
A dictionary with mesh IDs as keys and openmc.MeshBase
instanaces as values
"""
out = dict()
out = {}
for mesh_elem in elem.findall('mesh'):
mesh = MeshBase.from_xml_element(mesh_elem)
out[mesh.id] = mesh

2
openmc/mgxs/groups.py
View file

@ -17,7 +17,7 @@ class EnergyGroups:
The energy group boundaries in [eV] or the name of the group structure
(Must be a valid key in the openmc.mgxs.GROUP_STRUCTURES dictionary).
.. versionchanged:: 0.13.4
.. versionchanged:: 0.14.0
Changed to allow a string specifying the group structure name.
Attributes

9
openmc/mgxs/mdgxs.py
View file

@ -216,7 +216,7 @@ class MDGXS(MGXS):
group_edges = self.energy_groups.group_edges
energy_filter = openmc.EnergyFilter(group_edges)
if self.delayed_groups != None:
if self.delayed_groups is not None:
delayed_filter = openmc.DelayedGroupFilter(self.delayed_groups)
filters = [[energy_filter], [delayed_filter, energy_filter]]
else:
@ -743,9 +743,9 @@ class MDGXS(MGXS):
df.to_csv(filename + '.csv', index=False)
elif format == 'excel':
if self.domain_type == 'mesh':
df.to_excel(filename + '.xls')
df.to_excel(filename + '.xlsx')
else:
df.to_excel(filename + '.xls', index=False)
df.to_excel(filename + '.xlsx', index=False)
elif format == 'pickle':
df.to_pickle(filename + '.pkl')
elif format == 'latex':
@ -1884,9 +1884,6 @@ class DecayRate(MDGXS):
@property
def filters(self):
# Create the non-domain specific Filters for the Tallies
group_edges = self.energy_groups.group_edges
if self.delayed_groups is not None:
delayed_filter = openmc.DelayedGroupFilter(self.delayed_groups)
filters = [[delayed_filter], [delayed_filter]]

4
openmc/mgxs/mgxs.py
View file

@ -2001,9 +2001,9 @@ class MGXS:
df.to_csv(filename + '.csv', index=False)
elif format == 'excel':
if self.domain_type == 'mesh':
df.to_excel(filename + '.xls')
df.to_excel(filename + '.xlsx')
else:
df.to_excel(filename + '.xls', index=False)
df.to_excel(filename + '.xlsx', index=False)
elif format == 'pickle':
df.to_pickle(filename + '.pkl')
elif format == 'latex':

283
openmc/model/funcs.py
View file

@ -1,11 +1,10 @@
from collections.abc import Iterable
from functools import partial
from math import sqrt
from numbers import Real
from operator import attrgetter
from warnings import warn
from openmc import Plane, Cylinder, Universe, Cell
from openmc import Cylinder, Universe, Cell
from .surface_composite import RectangularPrism, HexagonalPrism
from ..checkvalue import (check_type, check_value, check_length,
check_less_than, check_iterable_type)
import openmc.data
@ -108,280 +107,26 @@ def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K',
return out
# Define function to create a plane on given axis
def _plane(axis, name, value, boundary_type='transmission'):
cls = getattr(openmc, f'{axis.upper()}Plane')
return cls(value, name=f'{name} {axis}',
boundary_type=boundary_type)
def rectangular_prism(width, height, axis='z', origin=(0., 0.),
boundary_type='transmission', corner_radius=0.):
"""Get an infinite rectangular prism from four planar surfaces.
.. versionchanged:: 0.11
This function was renamed from `get_rectangular_prism` to
`rectangular_prism`.
Parameters
----------
width: float
Prism width in units of cm. The width is aligned with the y, x,
or x axes for prisms parallel to the x, y, or z axis, respectively.
height: float
Prism height in units of cm. The height is aligned with the z, z,
or y axes for prisms parallel to the x, y, or z axis, respectively.
axis : {'x', 'y', 'z'}
Axis with which the infinite length of the prism should be aligned.
Defaults to 'z'.
origin: Iterable of two floats
Origin of the prism. The two floats correspond to (y,z), (x,z) or
(x,y) for prisms parallel to the x, y or z axis, respectively.
Defaults to (0., 0.).
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surfaces comprising the rectangular prism (default is 'transmission').
corner_radius: float
Prism corner radius in units of cm. Defaults to 0.
Returns
-------
openmc.Region
The inside of a rectangular prism
"""
check_type('width', width, Real)
check_type('height', height, Real)
check_type('corner_radius', corner_radius, Real)
check_value('axis', axis, ['x', 'y', 'z'])
check_type('origin', origin, Iterable, Real)
if axis == 'x':
x1, x2 = 'y', 'z'
elif axis == 'y':
x1, x2 = 'x', 'z'
else:
x1, x2 = 'x', 'y'
# Get cylinder class corresponding to given axis
cyl = getattr(openmc, f'{axis.upper()}Cylinder')
# Create rectangular region
min_x1 = _plane(x1, 'minimum', -width/2 + origin[0],
boundary_type=boundary_type)
max_x1 = _plane(x1, 'maximum', width/2 + origin[0],
boundary_type=boundary_type)
min_x2 = _plane(x2, 'minimum', -height/2 + origin[1],
boundary_type=boundary_type)
max_x2 = _plane(x2, 'maximum', height/2 + origin[1],
boundary_type=boundary_type)
if boundary_type == 'periodic':
min_x1.periodic_surface = max_x1
min_x2.periodic_surface = max_x2
prism = +min_x1 & -max_x1 & +min_x2 & -max_x2
# Handle rounded corners if given
if corner_radius > 0.:
if boundary_type == 'periodic':
raise ValueError('Periodic boundary conditions not permitted when '
'rounded corners are used.')
args = {'r': corner_radius, 'boundary_type': boundary_type}
args[x1 + '0'] = origin[0] - width/2 + corner_radius
args[x2 + '0'] = origin[1] - height/2 + corner_radius
x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
args[x1 + '0'] = origin[0] - width/2 + corner_radius
args[x2 + '0'] = origin[1] - height/2 + corner_radius
x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
args[x1 + '0'] = origin[0] - width/2 + corner_radius
args[x2 + '0'] = origin[1] + height/2 - corner_radius
x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args)
args[x1 + '0'] = origin[0] + width/2 - corner_radius
args[x2 + '0'] = origin[1] - height/2 + corner_radius
x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args)
args[x1 + '0'] = origin[0] + width/2 - corner_radius
args[x2 + '0'] = origin[1] + height/2 - corner_radius
x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args)
x1_min = _plane(x1, 'min', -width/2 + origin[0] + corner_radius,
boundary_type=boundary_type)
x1_max = _plane(x1, 'max', width/2 + origin[0] - corner_radius,
boundary_type=boundary_type)
x2_min = _plane(x2, 'min', -height/2 + origin[1] + corner_radius,
boundary_type=boundary_type)
x2_max = _plane(x2, 'max', height/2 + origin[1] - corner_radius,
boundary_type=boundary_type)
corners = (+x1_min_x2_min & -x1_min & -x2_min) | \
(+x1_min_x2_max & -x1_min & +x2_max) | \
(+x1_max_x2_min & +x1_max & -x2_min) | \
(+x1_max_x2_max & +x1_max & +x2_max)
prism = prism & ~corners
return prism
def get_rectangular_prism(*args, **kwargs):
warn("get_rectangular_prism(...) has been renamed rectangular_prism(...). "
"Future versions of OpenMC will not accept get_rectangular_prism.",
FutureWarning)
return rectangular_prism(*args, **kwargs)
warn("The rectangular_prism(...) function has been replaced by the "
"RectangularPrism(...) class. Future versions of OpenMC will not "
"accept rectangular_prism.", FutureWarning)
return -RectangularPrism(
width=width, height=height, axis=axis, origin=origin,
boundary_type=boundary_type, corner_radius=corner_radius)
def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
boundary_type='transmission', corner_radius=0.):
"""Create a hexagon region from six surface planes.
.. versionchanged:: 0.11
This function was renamed from `get_hexagonal_prism` to
`hexagonal_prism`.
Parameters
----------
edge_length : float
Length of a side of the hexagon in cm
orientation : {'x', 'y'}
An 'x' orientation means that two sides of the hexagon are parallel to
the x-axis and a 'y' orientation means that two sides of the hexagon are
parallel to the y-axis.
origin: Iterable of two floats
Origin of the prism. Defaults to (0., 0.).
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
Boundary condition that defines the behavior for particles hitting the
surfaces comprising the hexagonal prism (default is 'transmission').
corner_radius: float
Prism corner radius in units of cm. Defaults to 0.
Returns
-------
openmc.Region
The inside of a hexagonal prism
"""
l = edge_length
x, y = origin
if orientation == 'y':
right = openmc.XPlane(x + sqrt(3.)/2*l, boundary_type=boundary_type)
left = openmc.XPlane(x - sqrt(3.)/2*l, boundary_type=boundary_type)
c = sqrt(3.)/3.
# y = -x/sqrt(3) + a
upper_right = Plane(a=c, b=1., d=l+x*c+y, boundary_type=boundary_type)
# y = x/sqrt(3) + a
upper_left = Plane(a=-c, b=1., d=l-x*c+y, boundary_type=boundary_type)
# y = x/sqrt(3) - a
lower_right = Plane(a=-c, b=1., d=-l-x*c+y, boundary_type=boundary_type)
# y = -x/sqrt(3) - a
lower_left = Plane(a=c, b=1., d=-l+x*c+y, boundary_type=boundary_type)
prism = -right & +left & -upper_right & -upper_left & \
+lower_right & +lower_left
if boundary_type == 'periodic':
right.periodic_surface = left
upper_right.periodic_surface = lower_left
lower_right.periodic_surface = upper_left
elif orientation == 'x':
top = openmc.YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type)
bottom = openmc.YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type)
c = sqrt(3.)
# y = -sqrt(3)*(x - a)
upper_right = Plane(a=c, b=1., d=c*l+x*c+y, boundary_type=boundary_type)
# y = sqrt(3)*(x + a)
lower_right = Plane(a=-c, b=1., d=-c*l-x*c+y,
boundary_type=boundary_type)
# y = -sqrt(3)*(x + a)
lower_left = Plane(a=c, b=1., d=-c*l+x*c+y, boundary_type=boundary_type)
# y = sqrt(3)*(x + a)
upper_left = Plane(a=-c, b=1., d=c*l-x*c+y, boundary_type=boundary_type)
prism = -top & +bottom & -upper_right & +lower_right & \
+lower_left & -upper_left
if boundary_type == 'periodic':
top.periodic_surface = bottom
upper_right.periodic_surface = lower_left
lower_right.periodic_surface = upper_left
# Handle rounded corners if given
if corner_radius > 0.:
if boundary_type == 'periodic':
raise ValueError('Periodic boundary conditions not permitted when '
'rounded corners are used.')
c = sqrt(3.)/2
t = l - corner_radius/c
# Cylinder with corner radius and boundary type pre-applied
cyl1 = partial(openmc.ZCylinder, r=corner_radius,
boundary_type=boundary_type)
cyl2 = partial(openmc.ZCylinder, r=corner_radius/(2*c),
boundary_type=boundary_type)
if orientation == 'x':
x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t)
x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t)
x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t)
x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t)
x_min_in = cyl1(name='x min in', x0=x-t, y0=y)
x_max_in = cyl1(name='x max in', x0=x+t, y0=y)
x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l)
x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l)
x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l)
x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l)
x_min_out = cyl2(name='x min out', x0=x-l, y0=y)
x_max_out = cyl2(name='x max out', x0=x+l, y0=y)
corners = (+x_min_y_min_in & -x_min_y_min_out |
+x_min_y_max_in & -x_min_y_max_out |
+x_max_y_min_in & -x_max_y_min_out |
+x_max_y_max_in & -x_max_y_max_out |
+x_min_in & -x_min_out |
+x_max_in & -x_max_out)
elif orientation == 'y':
x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2)
x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2)
x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2)
x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2)
y_min_in = cyl1(name='y min in', x0=x, y0=y-t)
y_max_in = cyl1(name='y max in', x0=x, y0=y+t)
x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2)
x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2)
x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2)
x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2)
y_min_out = cyl2(name='y min out', x0=x, y0=y-l)
y_max_out = cyl2(name='y max out', x0=x, y0=y+l)
corners = (+x_min_y_min_in & -x_min_y_min_out |
+x_min_y_max_in & -x_min_y_max_out |
+x_max_y_min_in & -x_max_y_min_out |
+x_max_y_max_in & -x_max_y_max_out |
+y_min_in & -y_min_out |
+y_max_in & -y_max_out)
prism = prism & ~corners
return prism
warn("The hexagonal_prism(...) function has been replaced by the "
"HexagonalPrism(...) class. Future versions of OpenMC will not "
"accept hexagonal_prism.", FutureWarning)
return -HexagonalPrism(
edge_length=edge_length, orientation=orientation, origin=origin,
boundary_type=boundary_type, corner_radius=corner_radius)
def get_hexagonal_prism(*args, **kwargs):

66
openmc/model/model.py
View file

@ -7,10 +7,10 @@ from pathlib import Path
from numbers import Integral
from tempfile import NamedTemporaryFile
import warnings
import lxml.etree as ET
from typing import Optional, Dict
import h5py
import lxml.etree as ET
import openmc
import openmc._xml as xml
@ -253,7 +253,8 @@ class Model:
path : str or PathLike
Path to model.xml file
"""
tree = ET.parse(path)
parser = ET.XMLParser(huge_tree=True)
tree = ET.parse(path, parser=parser)
root = tree.getroot()
model = cls()
@ -499,8 +500,6 @@ class Model:
warnings.warn("remove_surfs kwarg will be deprecated soon, please "
"set the Geometry.merge_surfaces attribute instead.")
self.geometry.merge_surfaces = True
# Can be used to modify tallies in case any surfaces are redundant
redundant_surfaces = self.geometry.remove_redundant_surfaces()
# provide a memo to track which meshes have been written
mesh_memo = set()
@ -611,6 +610,7 @@ class Model:
this method creates the XML files and runs OpenMC via a system call. In
both cases this method returns the path to the last statepoint file
generated.
.. versionchanged:: 0.12
Instead of returning the final k-effective value, this function now
returns the path to the final statepoint written.
@ -788,7 +788,12 @@ class Model:
for i, vol_calc in enumerate(self.settings.volume_calculations):
vol_calc.load_results(f"volume_{i + 1}.h5")
# First add them to the Python side
self.geometry.add_volume_information(vol_calc)
if vol_calc.domain_type == "material" and self.materials:
for material in self.materials:
if material.id in vol_calc.volumes:
material.add_volume_information(vol_calc)
else:
self.geometry.add_volume_information(vol_calc)
# And now repeat for the C API
if self.is_initialized and vol_calc.domain_type == 'material':
@ -1016,3 +1021,54 @@ class Model:
"""
self._change_py_lib_attribs(names_or_ids, volume, 'material', 'volume')
def differentiate_depletable_mats(self, diff_volume_method: str):
"""Assign distribmats for each depletable material
.. versionadded:: 0.14.0
Parameters
----------
diff_volume_method : str
Specifies how the volumes of the new materials should be found.
Default is to 'divide equally' which divides the original material
volume equally between the new materials, 'match cell' sets the
volume of the material to volume of the cell they fill.
"""
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Extract all depletable materials which have multiple instances
distribmats = set(
[mat for mat in self.materials
if mat.depletable and mat.num_instances > 1])
if diff_volume_method == 'divide equally':
for mat in distribmats:
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
f"material with ID={mat.id}.")
mat.volume /= mat.num_instances
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
if diff_volume_method == 'divide equally':
cell.fill = [mat.clone() for _ in range(cell.num_instances)]
elif diff_volume_method == 'match cell':
for _ in range(cell.num_instances):
cell.fill = mat.clone()
if not cell.volume:
raise ValueError(
f"Volume of cell ID={cell.id} not specified. "
"Set volumes of cells prior to using "
"diff_volume_method='match cell'."
)
cell.fill.volume = cell.volume
if self.materials is not None:
self.materials = openmc.Materials(
self.geometry.get_all_materials().values()
)

325
openmc/model/surface_composite.py
View file

@ -1,8 +1,12 @@
from abc import ABC, abstractmethod
from collections.abc import Iterable
from copy import copy
from functools import partial
from math import sqrt, pi, sin, cos, isclose
from numbers import Real
import warnings
import operator
from typing import Sequence
import numpy as np
from scipy.spatial import ConvexHull, Delaunay
@ -50,13 +54,13 @@ class CompositeSurface(ABC):
"""Iterable of attribute names corresponding to underlying surfaces."""
@abstractmethod
def __pos__(self):
"""Return the positive half-space of the composite surface."""
@abstractmethod
def __neg__(self):
def __neg__(self) -> openmc.Region:
"""Return the negative half-space of the composite surface."""
def __pos__(self) -> openmc.Region:
"""Return the positive half-space of the composite surface."""
return ~(-self)
class CylinderSector(CompositeSurface):
"""Infinite cylindrical sector composite surface.
@ -120,10 +124,10 @@ class CylinderSector(CompositeSurface):
**kwargs):
if r2 <= r1:
raise ValueError(f'r2 must be greater than r1.')
raise ValueError('r2 must be greater than r1.')
if theta2 <= theta1:
raise ValueError(f'theta2 must be greater than theta1.')
raise ValueError('theta2 must be greater than theta1.')
phi1 = pi / 180 * theta1
phi2 = pi / 180 * theta2
@ -204,7 +208,7 @@ class CylinderSector(CompositeSurface):
offset.
"""
if theta >= 360. or theta <= 0:
raise ValueError(f'theta must be less than 360 and greater than 0.')
raise ValueError('theta must be less than 360 and greater than 0.')
theta1 = alpha
theta2 = alpha + theta
@ -288,10 +292,10 @@ class IsogonalOctagon(CompositeSurface):
# Side lengths
if r2 > r1 * sqrt(2):
raise ValueError(f'r2 is greater than sqrt(2) * r1. Octagon' +
raise ValueError('r2 is greater than sqrt(2) * r1. Octagon' +
' may be erroneous.')
if r1 > r2 * sqrt(2):
raise ValueError(f'r1 is greater than sqrt(2) * r2. Octagon' +
raise ValueError('r1 is greater than sqrt(2) * r2. Octagon' +
' may be erroneous.')
L_basis_ax = (r2 * sqrt(2) - r1)
@ -836,9 +840,6 @@ class Polygon(CompositeSurface):
def __neg__(self):
return self._region
def __pos__(self):
return ~self._region
@property
def _surface_names(self):
return self._surfnames
@ -1232,7 +1233,7 @@ class CruciformPrism(CompositeSurface):
multiple distances from the center. Equivalent to the 'gcross' derived
surface in Serpent.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -1309,5 +1310,297 @@ class CruciformPrism(CompositeSurface):
)
return openmc.Union(regions)
def __pos__(self):
return ~(-self)
# Define function to create a plane on given axis
def _plane(axis, name, value, boundary_type='transmission', albedo=1.0):
cls = getattr(openmc, f'{axis.upper()}Plane')
return cls(value, name=f'{name} {axis}',
boundary_type=boundary_type, albedo=albedo)
class RectangularPrism(CompositeSurface):
"""Infinite rectangular prism bounded by four planar surfaces.
.. versionadded:: 0.14.0
Parameters
----------
width : float
Prism width in units of [cm]. The width is aligned with the x, x, or z
axes for prisms parallel to the x, y, or z axis, respectively.
height : float
Prism height in units of [cm]. The height is aligned with the x, y, or z
axes for prisms parallel to the x, y, or z axis, respectively.
axis : {'x', 'y', 'z'}
Axis with which the infinite length of the prism should be aligned.
origin : Iterable of two floats
Origin of the prism. The two floats correspond to (y,z), (x,z) or (x,y)
for prisms parallel to the x, y or z axis, respectively.
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
Boundary condition that defines the behavior for particles hitting the
surfaces comprising the rectangular prism.
albedo : float, optional
Albedo of the prism's surfaces as a ratio of particle weight after
interaction with the surface to the initial weight. Values must be
positive. Only applicable if the boundary type is 'reflective',
'periodic', or 'white'.
corner_radius : float
Prism corner radius in units of [cm].
"""
_surface_names = ('min_x1', 'max_x1', 'min_x2', 'max_x2')
def __init__(
self,
width: float,
height: float,
axis: str = 'z',
origin: Sequence[float] = (0., 0.),
boundary_type: str = 'transmission',
albedo: float = 1.,
corner_radius: float = 0.
):
check_type('width', width, Real)
check_type('height', height, Real)
check_type('albedo', albedo, Real)
check_type('corner_radius', corner_radius, Real)
check_value('axis', axis, ('x', 'y', 'z'))
check_type('origin', origin, Iterable, Real)
if axis == 'x':
x1, x2 = 'y', 'z'
elif axis == 'y':
x1, x2 = 'x', 'z'
else:
x1, x2 = 'x', 'y'
# Get cylinder class corresponding to given axis
cyl = getattr(openmc, f'{axis.upper()}Cylinder')
# Create container for boundary arguments
bc_args = {'boundary_type': boundary_type, 'albedo': albedo}
# Create rectangular region
self.min_x1 = _plane(x1, 'minimum', -width/2 + origin[0], **bc_args)
self.max_x1 = _plane(x1, 'maximum', width/2 + origin[0], **bc_args)
self.min_x2 = _plane(x2, 'minimum', -height/2 + origin[1], **bc_args)
self.max_x2 = _plane(x2, 'maximum', height/2 + origin[1], **bc_args)
if boundary_type == 'periodic':
self.min_x1.periodic_surface = self.max_x1
self.min_x2.periodic_surface = self.max_x2
# Handle rounded corners if given
if corner_radius > 0.:
if boundary_type == 'periodic':
raise ValueError('Periodic boundary conditions not permitted when '
'rounded corners are used.')
args = {'r': corner_radius, 'boundary_type': boundary_type, 'albedo': albedo}
args[x1 + '0'] = origin[0] - width/2 + corner_radius
args[x2 + '0'] = origin[1] - height/2 + corner_radius
self.x1_min_x2_min = cyl(name=f'{x1} min {x2} min', **args)
args[x1 + '0'] = origin[0] - width/2 + corner_radius
args[x2 + '0'] = origin[1] + height/2 - corner_radius
self.x1_min_x2_max = cyl(name=f'{x1} min {x2} max', **args)
args[x1 + '0'] = origin[0] + width/2 - corner_radius
args[x2 + '0'] = origin[1] - height/2 + corner_radius
self.x1_max_x2_min = cyl(name=f'{x1} max {x2} min', **args)
args[x1 + '0'] = origin[0] + width/2 - corner_radius
args[x2 + '0'] = origin[1] + height/2 - corner_radius
self.x1_max_x2_max = cyl(name=f'{x1} max {x2} max', **args)
self.x1_min = _plane(x1, 'min', -width/2 + origin[0] + corner_radius,
**bc_args)
self.x1_max = _plane(x1, 'max', width/2 + origin[0] - corner_radius,
**bc_args)
self.x2_min = _plane(x2, 'min', -height/2 + origin[1] + corner_radius,
**bc_args)
self.x2_max = _plane(x2, 'max', height/2 + origin[1] - corner_radius,
**bc_args)
self._surface_names += (
'x1_min_x2_min', 'x1_min_x2_max', 'x1_max_x2_min',
'x1_max_x2_max', 'x1_min', 'x1_max', 'x2_min', 'x2_max'
)
def __neg__(self):
prism = +self.min_x1 & -self.max_x1 & +self.min_x2 & -self.max_x2
# Cut out corners if a corner radius was given
if hasattr(self, 'x1_min'):
corners = (
(+self.x1_min_x2_min & -self.x1_min & -self.x2_min) |
(+self.x1_min_x2_max & -self.x1_min & +self.x2_max) |
(+self.x1_max_x2_min & +self.x1_max & -self.x2_min) |
(+self.x1_max_x2_max & +self.x1_max & +self.x2_max)
)
prism &= ~corners
return prism
class HexagonalPrism(CompositeSurface):
"""Hexagonal prism comoposed of six planar surfaces
.. versionadded:: 0.14.0
Parameters
----------
edge_length : float
Length of a side of the hexagon in [cm]
orientation : {'x', 'y'}
An 'x' orientation means that two sides of the hexagon are parallel to
the x-axis and a 'y' orientation means that two sides of the hexagon are
parallel to the y-axis.
origin : Iterable of two floats
Origin of the prism.
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
Boundary condition that defines the behavior for particles hitting the
surfaces comprising the hexagonal prism.
albedo : float, optional
Albedo of the prism's surfaces as a ratio of particle weight after
interaction with the surface to the initial weight. Values must be
positive. Only applicable if the boundary type is 'reflective',
'periodic', or 'white'.
corner_radius : float
Prism corner radius in units of [cm].
"""
_surface_names = ('plane_max', 'plane_min', 'upper_right', 'upper_left',
'lower_right', 'lower_left')
def __init__(
self,
edge_length: float = 1.,
orientation: str = 'y',
origin: Sequence[float] = (0., 0.),
boundary_type: str = 'transmission',
albedo: float = 1.,
corner_radius: float = 0.
):
check_type('edge_length', edge_length, Real)
check_type('albedo', albedo, Real)
check_type('corner_radius', corner_radius, Real)
check_value('orientation', orientation, ('x', 'y'))
check_type('origin', origin, Iterable, Real)
l = edge_length
x, y = origin
# Create container for boundary arguments
bc_args = {'boundary_type': boundary_type, 'albedo': albedo}
if orientation == 'y':
# Left and right planes
self.plane_max = openmc.XPlane(x + sqrt(3.)/2*l, **bc_args)
self.plane_min = openmc.XPlane(x - sqrt(3.)/2*l, **bc_args)
c = sqrt(3.)/3.
# y = -x/sqrt(3) + a
self.upper_right = openmc.Plane(a=c, b=1., d=l+x*c+y, **bc_args)
# y = x/sqrt(3) + a
self.upper_left = openmc.Plane(a=-c, b=1., d=l-x*c+y, **bc_args)
# y = x/sqrt(3) - a
self.lower_right = openmc.Plane(a=-c, b=1., d=-l-x*c+y, **bc_args)
# y = -x/sqrt(3) - a
self.lower_left = openmc.Plane(a=c, b=1., d=-l+x*c+y, **bc_args)
elif orientation == 'x':
self.plane_max = openmc.YPlane(y + sqrt(3.)/2*l, **bc_args)
self.plane_min = openmc.YPlane(y - sqrt(3.)/2*l, **bc_args)
c = sqrt(3.)
# Upper-right surface: y = -sqrt(3)*(x - a)
self.upper_right = openmc.Plane(a=c, b=1., d=c*l+x*c+y, **bc_args)
# Lower-right surface: y = sqrt(3)*(x + a)
self.lower_right = openmc.Plane(a=-c, b=1., d=-c*l-x*c+y, **bc_args)
# Lower-left surface: y = -sqrt(3)*(x + a)
self.lower_left = openmc.Plane(a=c, b=1., d=-c*l+x*c+y, **bc_args)
# Upper-left surface: y = sqrt(3)*(x + a)
self.upper_left = openmc.Plane(a=-c, b=1., d=c*l-x*c+y, **bc_args)
# Handle periodic boundary conditions
if boundary_type == 'periodic':
self.plane_min.periodic_surface = self.plane_max
self.upper_right.periodic_surface = self.lower_left
self.lower_right.periodic_surface = self.upper_left
# Handle rounded corners if given
if corner_radius > 0.:
if boundary_type == 'periodic':
raise ValueError('Periodic boundary conditions not permitted '
'when rounded corners are used.')
c = sqrt(3.)/2
t = l - corner_radius/c
# Cylinder with corner radius and boundary type pre-applied
cyl1 = partial(openmc.ZCylinder, r=corner_radius, **bc_args)
cyl2 = partial(openmc.ZCylinder, r=corner_radius/(2*c), **bc_args)
if orientation == 'x':
self.x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t)
self.x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t)
self.x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t)
self.x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t)
self.min_in = cyl1(name='x min in', x0=x-t, y0=y)
self.max_in = cyl1(name='x max in', x0=x+t, y0=y)
self.x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l)
self.x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l)
self.x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l)
self.x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l)
self.min_out = cyl2(name='x min out', x0=x-l, y0=y)
self.max_out = cyl2(name='x max out', x0=x+l, y0=y)
elif orientation == 'y':
self.x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2)
self.x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2)
self.x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2)
self.x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2)
self.min_in = cyl1(name='y min in', x0=x, y0=y-t)
self.max_in = cyl1(name='y max in', x0=x, y0=y+t)
self.x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2)
self.x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2)
self.x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2)
self.x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2)
self.min_out = cyl2(name='y min out', x0=x, y0=y-l)
self.max_out = cyl2(name='y max out', x0=x, y0=y+l)
# Add to tuple of surface names
for s in ('in', 'out'):
self._surface_names += (
f'x_min_y_min_{s}', f'x_min_y_max_{s}',
f'x_max_y_min_{s}', f'x_max_y_max_{s}',
f'min_{s}', f'max_{s}')
def __neg__(self) -> openmc.Region:
prism = (
-self.plane_max & +self.plane_min &
-self.upper_right & -self.upper_left &
+self.lower_right & +self.lower_left
)
# Cut out corners if a corner radius was given
if hasattr(self, 'min_in'):
corners = (
+self.x_min_y_min_in & -self.x_min_y_min_out |
+self.x_min_y_max_in & -self.x_min_y_max_out |
+self.x_max_y_min_in & -self.x_max_y_min_out |
+self.x_max_y_max_in & -self.x_max_y_max_out |
+self.min_in & -self.min_out |
+self.max_in & -self.max_out
)
prism &= ~corners
return prism

11
openmc/plots.py
View file

@ -1,10 +1,10 @@
from collections.abc import Iterable, Mapping
from numbers import Integral, Real
from pathlib import Path
import lxml.etree as ET
from typing import Optional
import h5py
import lxml.etree as ET
import numpy as np
import openmc
@ -184,7 +184,7 @@ def _get_plot_image(plot, cwd):
def voxel_to_vtk(voxel_file: PathLike, output: PathLike = 'plot.vti'):
"""Converts a voxel HDF5 file to a VTK file
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -943,7 +943,7 @@ class Plot(PlotBase):
This method runs OpenMC in plotting mode to produce a .vti file.
.. versionadded:: 0.13.4
.. versionadded:: 0.14.0
Parameters
----------
@ -990,6 +990,8 @@ class ProjectionPlot(PlotBase):
projections are more similar to a pinhole camera, and orthographic projections
preserve parallel lines and distances.
.. versionadded:: 0.14.0
Parameters
----------
plot_id : int
@ -1506,6 +1508,7 @@ class Plots(cv.CheckedList):
Plots collection
"""
tree = ET.parse(path)
parser = ET.XMLParser(huge_tree=True)
tree = ET.parse(path, parser=parser)
root = tree.getroot()
return cls.from_xml_element(root)

40
openmc/region.py
View file

@ -30,6 +30,11 @@ class Region(ABC):
def __contains__(self, point):
pass
@property
@abstractmethod
def bounding_box(self) -> BoundingBox:
pass
@abstractmethod
def __str__(self):
pass
@ -105,17 +110,17 @@ class Region(ABC):
# If special character appears immediately after a non-operator,
# create a token with the appropriate half-space
if i_start >= 0:
# When an opening parenthesis appears after a non-operator,
# there's an implicit intersection operator between them
if expression[i] == '(':
tokens.append(' ')
j = int(expression[i_start:i])
if j < 0:
tokens.append(-surfaces[abs(j)])
else:
tokens.append(+surfaces[abs(j)])
# When an opening parenthesis appears after a non-operator,
# there's an implicit intersection operator between them
if expression[i] == '(':
tokens.append(' ')
if expression[i] in '()|~':
# For everything other than intersection, add the operator
# to the list of tokens
@ -413,14 +418,11 @@ class Intersection(Region, MutableSequence):
return '(' + ' '.join(map(str, self)) + ')'
@property
def bounding_box(self):
lower_left = np.array([-np.inf, -np.inf, -np.inf])
upper_right = np.array([np.inf, np.inf, np.inf])
def bounding_box(self) -> BoundingBox:
box = BoundingBox.infinite()
for n in self:
lower_left_n, upper_right_n = n.bounding_box
lower_left[:] = np.maximum(lower_left, lower_left_n)
upper_right[:] = np.minimum(upper_right, upper_right_n)
return BoundingBox(lower_left, upper_right)
box &= n.bounding_box
return box
class Union(Region, MutableSequence):
@ -504,14 +506,12 @@ class Union(Region, MutableSequence):
return '(' + ' | '.join(map(str, self)) + ')'
@property
def bounding_box(self):
lower_left = np.array([np.inf, np.inf, np.inf])
upper_right = np.array([-np.inf, -np.inf, -np.inf])
def bounding_box(self) -> BoundingBox:
bbox = BoundingBox(np.array([np.inf]*3),
np.array([-np.inf]*3))
for n in self:
lower_left_n, upper_right_n = n.bounding_box
lower_left[:] = np.minimum(lower_left, lower_left_n)
upper_right[:] = np.maximum(upper_right, upper_right_n)
return BoundingBox(lower_left, upper_right)
bbox |= n.bounding_box
return bbox
class Complement(Region):
@ -576,7 +576,7 @@ class Complement(Region):
self._node = node
@property
def bounding_box(self):
def bounding_box(self) -> BoundingBox:
# Use De Morgan's laws to distribute the complement operator so that it
# only applies to surface half-spaces, thus allowing us to calculate the
# bounding box in the usual recursive manner.

Some files were not shown because too many files have changed in this diff Show more