diff --git a/.gitignore b/.gitignore index dd8dfb14a9..32ef7919ab 100644 --- a/.gitignore +++ b/.gitignore @@ -25,6 +25,7 @@ examples/**/*.xml # Documentation builds docs/build +docs/doxygen/xml docs/source/_images/*.pdf docs/source/_images/*.aux docs/source/pythonapi/generated/ diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 61301bdb49..7594ed5f2b 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -7,6 +7,8 @@ build: jobs: post_checkout: - git fetch --unshallow || true + - cd docs/doxygen && doxygen && cd - + sphinx: configuration: docs/source/conf.py diff --git a/docs/Makefile b/docs/Makefile index a93338df31..e320432a20 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -45,6 +45,7 @@ help: clean: -rm -rf $(BUILDDIR)/* -rm -rf source/pythonapi/generated/ + -rm -rf doxygen/xml html: $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html diff --git a/docs/doxygen/Doxyfile b/docs/doxygen/Doxyfile new file mode 100644 index 0000000000..75702eff7c --- /dev/null +++ b/docs/doxygen/Doxyfile @@ -0,0 +1,13 @@ +# Doxyfile 1.9.1 + +# This file describes the settings to be used by the documentation system +# doxygen (www.doxygen.org) for a project. + +# Difference with default Doxyfile 1.9.1 +PROJECT_NAME = OpenMC +QUIET = YES +WARN_IF_UNDOCUMENTED = NO +INPUT = ../../include/openmc/capi.h +GENERATE_HTML = NO +GENERATE_LATEX = NO +GENERATE_XML = YES diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index 1777e795ad..e924aa6e2a 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -46,43 +46,20 @@ Type Definitions Functions --------- -.. c:function:: int openmc_calculate_volumes() +.. + Once documentation is complete in capi.h, use: + .. doxygenfile:: capi.h + to populate this documentation without using + .. doxygenfunction:: + for every function. - Run a stochastic volume calculation +.. doxygenfunction:: openmc_calculate_volumes - :return: Return status (negative if an error occurred) - :rtype: int +.. doxygenfunction:: openmc_cell_get_fill -.. c:function:: int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n) +.. doxygenfunction:: openmc_cell_get_id - Get the fill for a cell - - :param int32_t index: Index in the cells array - :param int* type: Type of the fill - :param int32_t** indices: Array of material indices for cell - :param int32_t* n: Length of indices array - :return: Return status (negative if an error occurred) - :rtype: int - -.. c:function:: int openmc_cell_get_id(int32_t index, int32_t* id) - - Get the ID of a cell - - :param int32_t index: Index in the cells array - :param int32_t* id: ID of the cell - :return: Return status (negative if an error occurred) - :rtype: int - -.. c:function:: int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T) - - Get the temperature of a cell - - :param int32_t index: Index in the cells array - :param int32_t* instance: Which instance of the cell. If a null pointer is passed, the temperature - of the first instance is returned. - :param double* T: temperature of the cell - :return: Return status (negative if an error occurred) - :rtype: int +.. doxygenfunction:: openmc_cell_get_temperature .. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* density) diff --git a/docs/source/conf.py b/docs/source/conf.py index 3a2e7fb93b..61740239bf 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -11,7 +11,10 @@ # All configuration values have a default; values that are commented out # serve to show the default. -import sys, os +import os +from pathlib import Path +import subprocess +import sys # Determine if we're on Read the Docs server on_rtd = os.environ.get('READTHEDOCS', None) == 'True' @@ -37,6 +40,7 @@ sys.path.insert(0, os.path.abspath('../..')) # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. extensions = [ + 'breathe', 'sphinx.ext.autodoc', 'sphinx.ext.napoleon', 'sphinx.ext.autosummary', @@ -47,6 +51,8 @@ extensions = [ ] if not on_rtd: extensions.append('sphinxcontrib.rsvgconverter') + doxygen_dir = Path(__file__).parents[1] / 'doxygen' + subprocess.run(['doxygen'], cwd=doxygen_dir, check=True) # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -117,6 +123,11 @@ pygments_style = 'tango' # A list of ignored prefixes for module index sorting. #modindex_common_prefix = [] +# -- Options breathe + doxygen ------------------------------------------------- + +breathe_projects = {"OpenMC": "../doxygen/xml"} +breathe_default_project = "OpenMC" +breathe_domain_by_file_pattern = {"*capi.h": "c"} # -- Options for HTML output --------------------------------------------------- diff --git a/docs/source/devguide/contributing.rst b/docs/source/devguide/contributing.rst index cda0313892..5199c9438f 100644 --- a/docs/source/devguide/contributing.rst +++ b/docs/source/devguide/contributing.rst @@ -111,7 +111,8 @@ The TC consists of the following individuals: - `Paul Romano `_ - `Patrick Shriwise `_ - `Adam Nelson `_ -- `Benoit Forget `_ +- `Jonathan Shimwell `_ +- `John Tramm `_ The Project Lead is Paul Romano. diff --git a/docs/source/devguide/docbuild.rst b/docs/source/devguide/docbuild.rst index f723db06ea..cda0307cae 100644 --- a/docs/source/devguide/docbuild.rst +++ b/docs/source/devguide/docbuild.rst @@ -14,6 +14,11 @@ Python API. That is, from the root directory of the OpenMC repository: python -m pip install ".[docs]" +The OpenMC documentation also uses Doxygen to automatically generate its +C/C++ API documentation directly from the docstrings available in the source +code. You will need to have a working installation of Doxygen to generate the +documentation locally. + ----------------------------------- Building Documentation as a Webpage ----------------------------------- diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index 60ff0afaef..5cd18bef16 100644 --- a/docs/source/io_formats/geometry.rst +++ b/docs/source/io_formats/geometry.rst @@ -425,13 +425,13 @@ Each ```` element can have the following attributes or sub-eleme material) assignment. Required. :temperature: - Temperature(s) in [K] to assign to the cell. Must be ≥ 0. Multiple - space-separated values may be given. + Temperature(s) in [K] to assign to the cell. Must be greater than or equal + to 0. Multiple space-separated values may be given. *Default*: None :density: - Density in [g/cm³] to assign to the cell. Must be > 0. Requires a non-void + Density in [g/cm³] to assign to the cell. Must be greater than 0. Requires a non-void material fill. Multiple space-separated values may be given. *Default*: None diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index fb02159169..7091757f23 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -777,9 +777,9 @@ attributes/sub-elements: *Default*: 1.0 :type: - 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. + Indicator of source type. One of ``independent``, ``file``, ``compiled``, + ``mesh``, or ``tokamak``. The type of the source will be determined by this + attribute if it is present. :particle: The source particle type, specified as a PDG number or a string alias (e.g., @@ -1015,6 +1015,80 @@ attributes/sub-elements: mesh element and follows the format for :ref:`source_element`. The number of ```` sub-elements should correspond to the number of mesh elements. + For a source with ``type="tokamak"``, the spatial distribution is described by + a Miller-style flux-surface parameterization and the following sub-elements + are used instead of the ``space`` element: + + :major_radius: + The major radius :math:`R_0` of the plasma in [cm]. + + :minor_radius: + The minor radius :math:`a` of the plasma in [cm]. Must be smaller than + ``major_radius``. + + :elongation: + The plasma elongation :math:`\kappa` (must be > 0). + + :triangularity: + The plasma triangularity :math:`\delta` (must be in [-1, 1]). Negative + values describe negative-triangularity plasmas. + + :shafranov_shift: + The Shafranov shift :math:`\Delta` in [cm] (must be >= 0 and less than + ``minor_radius``/2). + + :r_over_a: + A list of normalized minor-radius grid points :math:`r/a`. Must be strictly + increasing, start at 0, and end at 1. + + :emission_density: + A list of neutron emission densities :math:`S(r)` evaluated at each + ``r_over_a`` grid point (arbitrary units, must be non-negative). Only the + shape matters, since the profile is normalized internally. Values are + interpolated linearly between grid points and the profile is refined on an + internal grid for radial sampling. Must have the same length as + ``r_over_a`` and contain at least one positive value. + + :phi_start: + The starting toroidal angle in [rad]. + + *Default*: 0.0 + + :phi_extent: + The toroidal angle extent in [rad]. The source is sampled uniformly in + :math:`[\phi_\text{start},\ \phi_\text{start} + \phi_\text{extent}]`. + + *Default*: :math:`2\pi` + + :n_alpha: + The number of poloidal-angle grid points used to build the sampling CDFs + (must be > 2). Larger values reduce discretization bias; values below 51 + produce a warning. + + *Default*: 101 + + :vertical_shift: + A vertical shift of the plasma center in [cm]. + + *Default*: 0.0 + + :energy: + For a tokamak source, one or more ``energy`` sub-elements specify the + neutron energy distribution(s). Either a single distribution is given (used + at all radii) or exactly one distribution per ``r_over_a`` grid point is + given, in which case the energy is sampled from one of the two + distributions bracketing the sampled radius, selected stochastically with + probability proportional to the proximity of the radius to each grid point + (stochastic interpolation). Each follows the format of a univariate + probability distribution (see :ref:`univariate`). + + :time: + An optional ``time`` sub-element specifying the time distribution of source + particles, following the format of a univariate probability distribution + (see :ref:`univariate`). + + *Default*: particles are born at :math:`t=0` + .. note:: Biased sampling can be applied to the spatial and energy distributions of a source by using the ```` sub-element (see :ref:`univariate` for details on how to specify bias distributions). diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index b3911c8b3b..f06a1eba49 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -26,6 +26,7 @@ Simulation Settings openmc.FileSource openmc.CompiledSource openmc.MeshSource + openmc.TokamakSource openmc.SourceParticle openmc.VolumeCalculation openmc.Settings diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst index 21981fdaa8..19f1e5d5f9 100644 --- a/docs/source/usersguide/decay_sources.rst +++ b/docs/source/usersguide/decay_sources.rst @@ -132,8 +132,9 @@ can be run:: r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes) If not specified otherwise, a photon transport calculation is run at each time -in the depletion schedule. That means in the case above, we would see three -photon transport calculations. To specify specific times at which photon +in the depletion schedule for which a decay photon source exists. Times without +a decay photon source, such as the initial state of a model containing only +stable nuclides, are omitted. To specify particular times at which photon transport calculations should be run, pass the ``photon_time_indices`` argument. For example, if we wanted to run a photon transport calculation only on the last time (after the 5 hour decay), we would run:: @@ -141,6 +142,19 @@ time (after the 5 hour decay), we would run:: r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes, photon_time_indices=[2]) +To attribute photon tally results to their parent radionuclides, set +``by_parent_nuclide=True``. This automatically adds a +:class:`openmc.ParentNuclideFilter` to every photon tally that does not already +have one. The filter bins are the union of radionuclides contributing to the +prepared decay photon sources. The resulting bins can be used directly when +inspecting the tally results:: + + r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes, + photon_time_indices=[2], by_parent_nuclide=True) + + photon_tally = r2s.results['photon_tallies'][2][0] + tally_by_parent = photon_tally.get_pandas_dataframe() + After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager` instance will have a ``results`` dictionary that allows you to directly access results from each of the steps. It will also write out all the output files into @@ -148,12 +162,13 @@ a directory that is named "r2s_/". The ``output_dir`` argument to the :meth:`~openmc.deplete.R2SManager.run` method enables you to override the default output directory name if desired. -The :meth:`~openmc.deplete.R2SManager.run` method actually runs three +The :meth:`~openmc.deplete.R2SManager.run` method actually runs four lower-level methods under the hood:: r2s.step1_neutron_transport(...) r2s.step2_activation(...) - r2s.step3_photon_transport(...) + r2s.step3_photon_source(...) + r2s.step4_photon_transport(...) For users looking for more control over the calculation, these lower-level methods can be used in lieu of the :meth:`openmc.deplete.R2SManager.run` method. @@ -255,4 +270,3 @@ relevant tallies. This can be done with the aid of the # Apply time correction factors tally = d1s.apply_time_correction(dose_tally, factors, time_index) - diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index e8514b8561..6faeb59e19 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -291,6 +291,53 @@ example, the following would generate a photon source:: For a full list of all classes related to statistical distributions, see :ref:`pythonapi_stats`. +Tokamak Plasma Sources +---------------------- + +For fusion applications, the :class:`openmc.TokamakSource` class provides a +native parametric neutron source for tokamak plasmas. Rather than specifying +spatial, angular, and energy distributions separately, the source is defined by +the plasma geometry (using a `Miller-style flux-surface parameterization +`_) and a radial emission profile. Source +sites are sampled directly from the plasma volume without rejection. + +The plasma shape is described by the major radius :math:`R_0`, minor radius +:math:`a`, elongation :math:`\kappa`, triangularity :math:`\delta`, and +Shafranov shift :math:`\Delta`. The neutron birth profile is given as an +emission density :math:`S(r/a)` tabulated on a normalized minor-radius grid that +runs from 0 (magnetic axis) to 1 (last closed flux surface); only the shape of +the profile matters, since it is normalized internally. The emission density is +linearly interpolated between the supplied points and refined internally for +radial sampling. For example:: + + import numpy as np + + r_over_a = np.linspace(0.0, 1.0, 50) + emission = (1.0 - r_over_a**2)**2 # peaked on-axis profile + + source = openmc.TokamakSource( + major_radius=620.0, # cm + minor_radius=200.0, # cm + elongation=1.8, + triangularity=0.45, + shafranov_shift=10.0, # cm + r_over_a=r_over_a, + emission_density=emission, + energy=openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=20000.0), + ) + + settings.source = source + +The ``energy`` argument accepts either a single +:class:`~openmc.stats.Univariate` distribution applied at all radii, or a +sequence with one distribution per ``r_over_a`` grid point to model a +radially-varying neutron spectrum (energies are then sampled by stochastic +interpolation between the two distributions bracketing the sampled radius). A +time distribution can be given with the ``time`` argument; by default, particles +are born at :math:`t=0`. The toroidal extent can be restricted with +``phi_start`` and ``phi_extent`` to model a sector of the plasma, and +``vertical_shift`` translates the plasma center along the z-axis. + File-based Sources ------------------ diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 6b78145a4c..9f6987d74e 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -9,14 +9,41 @@ extern "C" { #endif +//! Run a stochastic volume calculation +// +//! \return Status (negative if an error occurred) int openmc_calculate_volumes(); + int openmc_cell_filter_get_bins( int32_t index, const int32_t** cells, int32_t* n); + +//! Get the fill for a cell +// +//! \param index Index in the cells array +//! \param type Type of the fill +//! \param indices Array of material indices for cell +//! \param n Length of indices array +//! \return Status (negative if an error occurred) int openmc_cell_get_fill( int32_t index, int* type, int32_t** indices, int32_t* n); + +//! Get the ID of a cell +// +//! \param index Index in the cells array +//! \param id ID of the cell +//! \return Status (negative if an error occurred) int openmc_cell_get_id(int32_t index, int32_t* id); + +//! Get the temperature of a cell +// +//! \param index Index in the cells array +//! \param instance Which instance of the cell. If a null pointer is +//! passed, the temperature of the first instance is returned. +//! \param T temperature of the cell +//!\return Status (negative if an error occurred) int openmc_cell_get_temperature( int32_t index, const int32_t* instance, double* T); + int openmc_cell_get_density( int32_t index, const int32_t* instance, double* rho); int openmc_cell_get_translation(int32_t index, double xyz[]); @@ -280,7 +307,7 @@ int openmc_zernike_filter_set_params( int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]); //! Sets the mesh and energy grid for CMFD reweight -//! \param[in] meshtyally_id id of CMFD Mesh Tally +//! \param[in] meshtally_id id of CMFD Mesh Tally //! \param[in] cmfd_indices indices storing spatial and energy dimensions of //! CMFD problem \param[in] norm CMFD normalization factor void openmc_initialize_mesh_egrid( diff --git a/include/openmc/error.h b/include/openmc/error.h index d73795aee2..1f6e15f49c 100644 --- a/include/openmc/error.h +++ b/include/openmc/error.h @@ -64,14 +64,14 @@ void write_message( int level, const std::string& message, const Params&... fmt_args) { if (settings::verbosity >= level) { - write_message(fmt::format(message, fmt_args...)); + write_message(fmt::format(fmt::runtime(message), fmt_args...)); } } template void write_message(const std::string& message, const Params&... fmt_args) { - write_message(fmt::format(message, fmt_args...)); + write_message(fmt::format(fmt::runtime(message), fmt_args...)); } #ifdef OPENMC_MPI diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 107cc7d1f3..43cf9bb585 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -3,9 +3,12 @@ #include #include +#include +#include #include "openmc/array.h" #include "openmc/constants.h" +#include "openmc/random_ray/source_region.h" // For hash_combine #include "openmc/vector.h" namespace openmc { @@ -13,6 +16,34 @@ namespace openmc { class BoundaryInfo; class GeometryState; +//============================================================================== +//! OverlapKey to store cell and universe data of a single overlap, along with +//! a functor for hashing an OverlapKey into an unordered_map. +//============================================================================== + +struct OverlapKey { + int universe_id; + int cell1_id; + int cell2_id; + + bool operator==(const OverlapKey& other) const + { + return universe_id == other.universe_id && cell1_id == other.cell1_id && + cell2_id == other.cell2_id; + } +}; + +struct OverlapKeyHash { + std::size_t operator()(const OverlapKey& k) const + { + size_t seed = 0; + hash_combine(seed, k.universe_id); + hash_combine(seed, k.cell1_id); + hash_combine(seed, k.cell2_id); + return seed; + } +}; + //============================================================================== // Global variables //============================================================================== @@ -24,6 +55,10 @@ extern "C" int n_coord_levels; //!< Number of CSG coordinate levels extern vector overlap_check_count; +// Overlap data structures get cleared every slice_data run +extern vector overlap_keys; +extern std::unordered_map overlap_key_index; + } // namespace model //============================================================================== @@ -38,8 +73,7 @@ inline bool coincident(double d1, double d2) //============================================================================== //! Check for overlapping cells at a particle's position. //============================================================================== - -bool check_cell_overlap(GeometryState& p, bool error = true); +int check_cell_overlap(GeometryState& p, bool error = true); //============================================================================== //! Get the cell instance for a particle at the specified universe level diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index dcc0e21fe2..d3bccca7b5 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -213,6 +213,20 @@ double exprel(double x); //! \return log(1+x)/x without loss of precision near 0 double log1prel(double x); +//! Evaluate the cylindrical Bessel function of the first kind J_n(x) +//! +//! Uses std::cyl_bessel_j where available (e.g., libstdc++). On standard +//! library implementations lacking the C++17 special math functions (e.g., +//! libc++ on Apple Clang/LLVM), falls back to the ascending power series, +//! which converges to machine precision for the small arguments (|x| <= 2) +//! used in OpenMC. Unlike std::cyl_bessel_j, negative arguments are handled +//! via the parity relation J_n(-x) = (-1)^n J_n(x). +//! +//! \param n Non-negative integer order of the Bessel function +//! \param x Real argument +//! \return J_n(x) +double cyl_bessel_j(int n, double x); + //! Helper function to get index and interpolation function on an incident //! energy grid //! @@ -233,5 +247,17 @@ void get_energy_index( double standard_normal_cdf(double z); +//============================================================================== +//! Return true if two floating-point values are approximately equal within a +//! combined relative and absolute tolerance. +//! +//! \param a first floating point value +//! \param b second floating point value +//! \param rel_tol relative tolerance +//! \param abs_tol absolute tolerance +//! \return true if a and b are approximately equal, false otherwise +//============================================================================== +bool isclose(double a, double b, double rel_tol, double abs_tol); + } // namespace openmc #endif // OPENMC_MATH_FUNCTIONS_H diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index f72948f6eb..44e82fd235 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -582,10 +582,8 @@ public: // Methods and accessors // Cross section caches - NuclideMicroXS& neutron_xs(int i) - { - return neutron_xs_[i]; - } // Microscopic neutron cross sections + // Microscopic neutron cross sections + NuclideMicroXS& neutron_xs(int i) { return neutron_xs_[i]; } const NuclideMicroXS& neutron_xs(int i) const { return neutron_xs_[i]; } // Microscopic photon cross sections diff --git a/include/openmc/plot.h b/include/openmc/plot.h index f97d313847..ba8ac84a1c 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -5,6 +5,7 @@ #include #include #include +#include #include "openmc/tensor.h" #include "pugixml.hpp" @@ -155,7 +156,7 @@ struct IdData { // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); - void set_overlap(size_t y, size_t x); + void set_overlap(size_t y, size_t x, int overlap_idx); // Members tensor::Tensor data_; //!< 2D array of cell & material ids @@ -168,7 +169,7 @@ struct PropertyData { // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); - void set_overlap(size_t y, size_t x); + void set_overlap(size_t y, size_t x, int overlap_idx); // Members tensor::Tensor data_; //!< 2D array of temperature & density data @@ -181,7 +182,7 @@ struct RasterData { // Methods void set_value(size_t y, size_t x, const Particle& p, int level, Filter* filter = nullptr, FilterMatch* match = nullptr); - void set_overlap(size_t y, size_t x); + void set_overlap(size_t y, size_t x, int overlap_idx); // Members tensor::Tensor @@ -278,8 +279,11 @@ T SlicePlotBase::get_map(int32_t filter_index) const if (found_cell) { data.set_value(y, x, p, j, filter, &match); } - if (show_overlaps_ && check_cell_overlap(p, false)) { - data.set_overlap(y, x); + if (show_overlaps_) { + int overlap_idx = check_cell_overlap(p, false); + if (overlap_idx >= 0) { + data.set_overlap(y, x, overlap_idx); + } } } // inner for } diff --git a/include/openmc/shared_array.h b/include/openmc/shared_array.h index b309ca3f1c..2829a2eb93 100644 --- a/include/openmc/shared_array.h +++ b/include/openmc/shared_array.h @@ -114,6 +114,32 @@ public: data_[size_++] = value; } + //! Increase the size of the container by count elements without assigning + //! values to the new elements. Existing elements are preserved if the + //! container needs to grow. This does not perform any thread safety checks. + // + //! \param count The number of elements to append + //! \return The starting index of the appended range + int64_t extend_uninitialized(int64_t count) + { + int64_t offset = size_; + int64_t new_size = size_ + count; + if (new_size > capacity_) { + int64_t new_capacity = capacity_ == 0 ? 8 : capacity_; + while (new_capacity < new_size) { + new_capacity *= 2; + } + unique_ptr new_data = make_unique(new_capacity); + if (size_ > 0) { + std::copy_n(data_.get(), size_, new_data.get()); + } + data_ = std::move(new_data); + capacity_ = new_capacity; + } + size_ = new_size; + return offset; + } + //! Return the number of elements in the container int64_t size() { return size_; } int64_t size() const { return size_; } diff --git a/include/openmc/source.h b/include/openmc/source.h index e307b1ed21..51b54a1d10 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -7,9 +7,12 @@ #include #include #include +#include // for pair #include "pugixml.hpp" +#include "openmc/array.h" +#include "openmc/distribution.h" #include "openmc/distribution_multi.h" #include "openmc/distribution_spatial.h" #include "openmc/memory.h" @@ -260,6 +263,139 @@ private: vector> sources_; //!< Source distributions }; +//============================================================================== +//! Parametric tokamak plasma neutron source +//! +//! This source samples neutron positions from a tokamak plasma geometry using +//! Miller-style flux surface parameterization with user-specified emission +//! profiles and energy distributions. +//! +//! Flux surface parameterization: +//! R = R0 + r*cos(alpha + delta*sin(alpha)) + Delta*(1 - (r/a)^2) +//! Z = kappa * r * sin(alpha) +//! +//! The sampling algorithm: +//! 1. Sample minor radius r from precomputed CDF of S(r) * Jacobian +//! 2. Sample poloidal angle alpha from conditional P(alpha|r) using mixture +//! of precomputed CDFs weighted by functions of r +//! 3. Sample energy and time from user-provided distribution(s) +//! 4. Sample isotropic direction +//! 5. Sample toroidal angle phi uniformly in [phi_start, phi_start + +//! phi_extent] +//! 6. Transform (r, alpha, phi) to Cartesian (x, y, z), applying the optional +//! vertical shift of the plasma center +//! +//! The user provides the emission density S(r) directly (e.g., from transport +//! codes like TRANSP, ASTRA, etc.), allowing full flexibility in reaction +//! physics calculations. S(r) is a profile in arbitrary units sampled on the +//! r_over_a grid; only its shape matters, since it is normalized internally. +//! Energy distributions can be specified as either a single distribution for +//! all r, or one distribution per radial point. +//============================================================================== + +class TokamakSource : public Source { +public: + // Constructors + explicit TokamakSource(pugi::xml_node node); + + //! Sample from the tokamak source distribution + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + SourceSite sample(uint64_t* seed) const override; + +private: + //========================================================================== + // Private methods + + //! Precompute data structures for efficient sampling + void precompute_sampling_distributions(); + + //! Sample minor radius from marginal CDF + //! \param seed Pseudorandom seed pointer + //! \return Sampled r/a value + double sample_r_over_a(uint64_t* seed) const; + + //! Sample poloidal angle given r using mixture of precomputed CDFs + //! \param r_norm Normalized minor radius r/a + //! \param seed Pseudorandom seed pointer + //! \return Sampled poloidal angle alpha [rad] + double sample_poloidal_angle(double r_norm, uint64_t* seed) const; + + //! Compute the k-th mixture weight w_k(r) * I_hat_k for poloidal sampling + //! \param k Basis function index (0-5) + //! \param r Normalized minor radius r/a + //! \return Mixture weight for component k + double mixture_weight(int k, double r) const; + + //! Sample energy from the distribution(s) + //! \param r_norm Normalized minor radius r/a (for distribution selection) + //! \param seed Pseudorandom seed pointer + //! \return (Sampled energy [eV], importance weight) + std::pair sample_energy(double r_norm, uint64_t* seed) const; + + //! Transform from flux coordinates (r, alpha, phi) to Cartesian (x, y, z) + //! \param r Minor radius [cm] + //! \param alpha Poloidal angle [rad] + //! \param phi Toroidal angle [rad] + //! \return Position in Cartesian coordinates [cm] + Position flux_to_cartesian(double r, double alpha, double phi) const; + + //========================================================================== + // Data members + + // Emission profile (input) + vector r_over_a_; //!< Normalized minor radius grid points + vector emission_density_; //!< Emission density S(r) at grid points + + // Energy distribution(s): either 1 for all r, or one per r point + vector> energy_dists_; + + // Time distribution (defaults to a delta distribution at t=0) + UPtrDist time_; + + // Angular distribution (isotropic) + UPtrAngle angle_; + + // Tokamak geometry parameters + double major_radius_; //!< Major radius R0 [cm] + double minor_radius_; //!< Minor radius a [cm] + double elongation_; //!< Elongation kappa + double triangularity_; //!< Triangularity delta + double shafranov_shift_; //!< Shafranov shift Delta [cm] + double vertical_shift_; //!< Vertical shift of plasma center [cm] + + // Normalized geometry parameters (precomputed for efficiency) + double epsilon_; //!< Inverse aspect ratio a/R0 + double delta_tilde_; //!< Normalized Shafranov shift Delta/a + + // Toroidal angle bounds + double phi_start_; //!< Starting toroidal angle [rad] + double phi_extent_; //!< Toroidal angle extent [rad] + + // Precomputed distribution for radial sampling + unique_ptr radial_dist_; + + // Coefficients of the radial geometric polynomial: A*r - B*r^2 - C*r^3 + // Also used as the analytical normalization for poloidal mixture weights + double radial_poly_a_; //!< 1 + ε·Δ̃ + double radial_poly_b_; //!< (3/8)·c₁·ε + double radial_poly_c_; //!< 2·ε·Δ̃ + + // Precomputed Bernstein basis functions for poloidal angle sampling. + // Using the factorization f(r_tilde, alpha) = R_tilde x J_tilde where: + // R_tilde = b0*(1-r)^2 + 2*b1*r*(1-r) + b2*r^2 (Bernstein quadratic) + // J_tilde = b3*(1-r) + b4*r (Bernstein linear) + // The product gives 6 non-negative basis functions g_k(alpha) with + // weights w_k(r_tilde) that are products of Bernstein polynomials. + // Distributions are tabulated on [0, pi] exploiting up-down symmetry. + static constexpr int N_POLOIDAL_BASIS = 6; //!< Number of basis functions + int n_alpha_; //!< Number of poloidal angle grid points + array, N_POLOIDAL_BASIS> + poloidal_dists_; //!< Distributions for each basis function g_k(alpha) + array + poloidal_integrals_; //!< Integrals of g_k(alpha) over [0, pi] +}; + //============================================================================== // Functions //============================================================================== diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index a5d404133c..d0b385d169 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -52,8 +52,12 @@ struct WeightWindow { double weight_cutoff {DEFAULT_WEIGHT_CUTOFF}; int max_split {10}; - //! Whether the weight window is in a valid state - bool is_valid() const { return lower_weight >= 0.0; } + //! Whether the weight window is in a valid state. A non-positive lower + //! bound indicates that no weight window information exists at this + //! location (generators mark such cells with -1, and a lower bound of zero + //! conventionally turns the weight window game off in a cell, as in MCNP + //! wwinp files), in which case no weight window game is played. + bool is_valid() const { return lower_weight > 0.0; } //! Adjust the weight window by a constant factor void scale(double factor) diff --git a/openmc/data/data.py b/openmc/data/data.py index c22e54e7dc..6e48409907 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -1,14 +1,23 @@ +from __future__ import annotations + import itertools import json import os import re from pathlib import Path from math import sqrt, log +from typing import TYPE_CHECKING, Literal from warnings import warn from endf.data import (ATOMIC_NUMBER, ATOMIC_SYMBOL, ELEMENT_SYMBOL, EV_PER_MEV, K_BOLTZMANN, gnds_name, zam) +import openmc +from openmc.checkvalue import PathLike + +if TYPE_CHECKING: + from openmc.deplete import Chain + gnds_name.__module__ = __name__ zam.__module__ = __name__ @@ -296,25 +305,47 @@ def atomic_weight(element): raise ValueError(f"No naturally-occurring isotopes for element '{element}'.") -def half_life(isotope): +def half_life( + isotope: str, + chain_file: Literal[False] | None | PathLike | Chain = False +) -> float | None: """Return half-life of isotope in seconds or None if isotope is stable - Half-life values are from the `ENDF/B-VIII.0 decay sublibrary - `_. + By default, half-life values are from the `ENDF/B-VIII.0 decay sublibrary + `_. A depletion chain can + also be used as the source of half-life values. .. versionadded:: 0.13.1 + .. versionchanged:: 0.15.4 + Added the ``chain_file`` argument. + Parameters ---------- isotope : str Name of isotope, e.g., 'Pu239' + chain_file : False, None, PathLike, or openmc.deplete.Chain, optional + Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is + used. If ``None``, the chain specified by + ``openmc.config['chain_file']`` is used when available. If a path or + :class:`openmc.deplete.Chain` is given, that chain is used. For ``None`` + or an explicit chain, nuclides absent from the chain fall back to + ENDF/B-VIII.0 data. Returns ------- - float - Half-life of isotope in [s] + float or None + Half-life of isotope in [s], or None if the isotope is stable """ + if chain_file is not False: + if chain_file is not None or openmc.config.get('chain_file') is not None: + # Local import avoids a circular dependency + from openmc.deplete.chain import _get_chain + chain = _get_chain(chain_file) + if isotope in chain: + return chain[isotope].half_life + global _HALF_LIFE if not _HALF_LIFE: # Load ENDF/B-VIII.0 data from JSON file @@ -324,7 +355,10 @@ def half_life(isotope): return _HALF_LIFE.get(isotope.lower()) -def decay_constant(isotope): +def decay_constant( + isotope: str, + chain_file: Literal[False] | None | PathLike | Chain = False +) -> float: """Return decay constant of isotope in [s^-1] Decay constants are based on half-life values from the @@ -333,10 +367,20 @@ def decay_constant(isotope): .. versionadded:: 0.13.1 + .. versionchanged:: 0.15.4 + Added the ``chain_file`` argument. + Parameters ---------- isotope : str Name of isotope, e.g., 'Pu239' + chain_file : False, None, PathLike, or openmc.deplete.Chain, optional + Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is + used. If ``None``, the chain specified by + ``openmc.config['chain_file']`` is used when available. If a path or + :class:`openmc.deplete.Chain` is given, that chain is used. For ``None`` + or an explicit chain, nuclides absent from the chain fall back to + ENDF/B-VIII.0 data. Returns ------- @@ -348,7 +392,7 @@ def decay_constant(isotope): openmc.data.half_life """ - t = half_life(isotope) + t = half_life(isotope, chain_file) return _LOG_TWO / t if t else 0.0 @@ -496,5 +540,3 @@ def isotopes(element: str) -> list[tuple[str, float]]: result.append(kv) return result - - diff --git a/openmc/data/dose/dose.py b/openmc/data/dose/dose.py index d49043b0a6..87e0cf9aef 100644 --- a/openmc/data/dose/dose.py +++ b/openmc/data/dose/dose.py @@ -4,50 +4,76 @@ import numpy as np import openmc.checkvalue as cv -_FILES = { - ('icrp74', 'neutron'): Path('icrp74') / 'neutrons.txt', - ('icrp74', 'photon'): Path('icrp74') / 'photons.txt', - ('icrp116', 'electron'): Path('icrp116') / 'electrons.txt', - ('icrp116', 'helium'): Path('icrp116') / 'helium_ions.txt', - ('icrp116', 'mu-'): Path('icrp116') / 'negative_muons.txt', - ('icrp116', 'pi-'): Path('icrp116') / 'negative_pions.txt', - ('icrp116', 'neutron'): Path('icrp116') / 'neutrons.txt', - ('icrp116', 'photon'): Path('icrp116') / 'photons.txt', - ('icrp116', 'photon kerma'): Path('icrp116') / 'photons_kerma.txt', - ('icrp116', 'mu+'): Path('icrp116') / 'positive_muons.txt', - ('icrp116', 'pi+'): Path('icrp116') / 'positive_pions.txt', - ('icrp116', 'positron'): Path('icrp116') / 'positrons.txt', - ('icrp116', 'proton'): Path('icrp116') / 'protons.txt', +_FULL_GEOMETRIES = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO') +_LIMITED_GEOMETRIES = ('AP', 'PA', 'ISO') + +_TABLES = { + ('icrp74', 'effective', 'neutron'): ( + Path('icrp74') / 'neutrons.txt', _FULL_GEOMETRIES), + ('icrp74', 'effective', 'photon'): ( + Path('icrp74') / 'photons.txt', _FULL_GEOMETRIES), + ('icrp116', 'effective', 'electron'): ( + Path('icrp116') / 'electrons.txt', _LIMITED_GEOMETRIES), + ('icrp116', 'effective', 'helium'): ( + Path('icrp116') / 'helium_ions.txt', _LIMITED_GEOMETRIES), + ('icrp116', 'effective', 'mu-'): ( + Path('icrp116') / 'negative_muons.txt', _LIMITED_GEOMETRIES), + ('icrp116', 'effective', 'pi-'): ( + Path('icrp116') / 'negative_pions.txt', _LIMITED_GEOMETRIES), + ('icrp116', 'effective', 'neutron'): ( + Path('icrp116') / 'neutrons.txt', _FULL_GEOMETRIES), + ('icrp116', 'effective', 'photon'): ( + Path('icrp116') / 'photons.txt', _FULL_GEOMETRIES), + ('icrp116', 'effective', 'photon kerma'): ( + Path('icrp116') / 'photons_kerma.txt', _FULL_GEOMETRIES), + ('icrp116', 'effective', 'mu+'): ( + Path('icrp116') / 'positive_muons.txt', _LIMITED_GEOMETRIES), + ('icrp116', 'effective', 'pi+'): ( + Path('icrp116') / 'positive_pions.txt', _LIMITED_GEOMETRIES), + ('icrp116', 'effective', 'positron'): ( + Path('icrp116') / 'positrons.txt', _LIMITED_GEOMETRIES), + ('icrp116', 'effective', 'proton'): ( + Path('icrp116') / 'protons.txt', _FULL_GEOMETRIES), + ('icrp74', 'ambient', 'neutron'): ( + Path('icrp74') / 'neutrons_H10.txt', None), + ('icrp74', 'ambient', 'photon'): ( + Path('icrp74') / 'photons_H10.txt', None), } _DOSE_TABLES = {} -def _load_dose_icrp(data_source: str, particle: str): - """Load effective dose tables from text files. +def _load_dose_table(data_source: str, dose_quantity: str, particle: str): + """Load dose tables from text files. Parameters ---------- data_source : {'icrp74', 'icrp116'} The dose conversion data source to use + dose_quantity : {'effective', 'ambient'} + Dose quantity to load. 'ambient' corresponds to ambient dose + equivalent H*(10). particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'} Incident particle """ - path = Path(__file__).parent / _FILES[data_source, particle] + key = (data_source, dose_quantity, particle) + path = Path(__file__).parent / _TABLES[key][0] data = np.loadtxt(path, skiprows=3, encoding='utf-8') data[:, 0] *= 1e6 # Change energies to eV - _DOSE_TABLES[data_source, particle] = data + _DOSE_TABLES[key] = data -def dose_coefficients(particle, geometry='AP', data_source='icrp116'): - """Return effective dose conversion coefficients. +def dose_coefficients( + particle, geometry='AP', data_source='icrp116', dose_quantity='effective' +): + """Return dose conversion coefficients. - This function provides fluence (and air kerma) to effective or ambient dose - (H*(10)) conversion coefficients for various types of external exposures - based on values in ICRP publications. Corrected values found in a - corrigendum are used rather than the values in the original report. - Available libraries include `ICRP Publication 74 + This function provides fluence (and air kerma) to effective dose or ambient + dose equivalent (H*(10)) conversion coefficients for various types of + external exposures based on values in ICRP publications. Corrected values + found in a corrigendum are used rather than the values in the original + report. Available libraries include `ICRP Publication 74 ` and `ICRP Publication 116 `. @@ -63,45 +89,58 @@ def dose_coefficients(particle, geometry='AP', data_source='icrp116'): particle : {'neutron', 'photon', 'photon kerma', 'electron', 'positron'} Incident particle geometry : {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'} - Irradiation geometry assumed. Refer to ICRP-116 (Section 3.2) for the - meaning of the options here. + Irradiation geometry assumed for effective dose coefficients. Refer to + ICRP-116 (Section 3.2) for the meaning of the options here. This + argument does not apply when ``dose_quantity`` is 'ambient'. data_source : {'icrp74', 'icrp116'} - The data source for the effective dose conversion coefficients. + The data source for the dose conversion coefficients. + dose_quantity : {'effective', 'ambient'} + Dose quantity to return. 'effective' returns effective dose + coefficients; 'ambient' returns ambient dose equivalent (H*(10)) + coefficients. Returns ------- energy : numpy.ndarray Energies at which dose conversion coefficients are given dose_coeffs : numpy.ndarray - Effective dose coefficients in [pSv cm^2] at provided energies. For - 'photon kerma', the coefficients are given in [Sv/Gy]. + Dose coefficients in [pSv cm^2] at provided energies. For 'photon + kerma', the coefficients are given in [Sv/Gy]. """ - cv.check_value('geometry', geometry, {'AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO'}) + cv.check_value('geometry', geometry, _FULL_GEOMETRIES) cv.check_value('data_source', data_source, {'icrp74', 'icrp116'}) + cv.check_value('dose_quantity', dose_quantity, {'effective', 'ambient'}) - if (data_source, particle) not in _FILES: - available_particles = sorted({p for (ds, p) in _FILES if ds == data_source}) + key = (data_source, dose_quantity, particle) + if key not in _TABLES: + available_particles = sorted( + p for ds, dq, p in _TABLES + if ds == data_source and dq == dose_quantity + ) msg = ( - f"'{particle}' has no dose data in data source {data_source}. " - f"Available particles for {data_source} are: {available_particles}" + f"'{particle}' has no {dose_quantity} dose data in data source " + f"{data_source}. Available particles for {data_source} " + f"with dose quantity {dose_quantity} are: {available_particles}" ) raise ValueError(msg) - elif (data_source, particle) not in _DOSE_TABLES: - _load_dose_icrp(data_source, particle) + elif key not in _DOSE_TABLES: + _load_dose_table(data_source, dose_quantity, particle) - # Get all data for selected particle - data = _DOSE_TABLES[data_source, particle] + data = _DOSE_TABLES[key] + columns = _TABLES[key][1] - # Determine index for selected geometry - if particle in ('neutron', 'photon', 'proton', 'photon kerma'): - columns = ('AP', 'PA', 'LLAT', 'RLAT', 'ROT', 'ISO') + if columns is None: + if geometry != 'AP': + raise ValueError( + "Irradiation geometry is not defined for ambient dose " + "equivalent coefficients. Use the default geometry='AP'." + ) + index = 0 else: - columns = ('AP', 'PA', 'ISO') - index = columns.index(geometry) + index = columns.index(geometry) - # Pull out energy and dose from table energy = data[:, 0].copy() dose_coeffs = data[:, index + 1].copy() return energy, dose_coeffs diff --git a/openmc/data/dose/icrp74/neutrons_H10.txt b/openmc/data/dose/icrp74/neutrons_H10.txt new file mode 100644 index 0000000000..fe0036cbc9 --- /dev/null +++ b/openmc/data/dose/icrp74/neutrons_H10.txt @@ -0,0 +1,50 @@ +Neutrons: Ambient per fluence, in units of pSv cm², for monoenergetic particles incident. + +Energy (MeV) Dose + 1.00E-09 6.60 + 1.00E-08 9.00 + 2.53E-08 10.6 + 1.00E-07 12.9 + 2.00E-07 13.5 + 5.00E-07 13.6 + 1.00E-06 13.3 + 2.00E-06 12.9 + 5.00E-06 12.0 + 1.00E-05 11.3 + 2.00E-05 10.6 + 5.00E-05 9.90 + 1.00E-04 9.40 + 2.00E-04 8.90 + 5.00E-04 8.30 + 1.00E-03 7.90 + 2.00E-03 7.70 + 5.00E-03 8.00 + 1.00E-02 10.5 + 2.00E-02 16.6 + 3.00E-02 23.7 + 5.00E-02 41.1 + 7.00E-02 60.0 + 1.00E-01 88.0 + 1.50E-01 132 + 2.00E-01 170 + 3.00E-01 233 + 5.00E-01 322 + 7.00E-01 375 + 9.00E-01 400 + 1 416 + 1.2 425 + 2 420 + 3 412 + 4 408 + 5 405 + 6 400 + 7 405 + 8 409 + 9 420 + 10 440 + 12 480 + 14 520 + 15 540 + 16 555 + 18 570 + 20 600 \ No newline at end of file diff --git a/openmc/data/dose/icrp74/photons_H10.txt b/openmc/data/dose/icrp74/photons_H10.txt new file mode 100644 index 0000000000..66031b2f0e --- /dev/null +++ b/openmc/data/dose/icrp74/photons_H10.txt @@ -0,0 +1,28 @@ +Photons: Ambient dose (H*10) per fluence, in units of pSv cm² + +Energy (MeV) Dose +0.010 0.061 +0.015 0.83 +0.020 1.05 +0.030 0.81 +0.040 0.64 +0.050 0.55 +0.060 0.51 +0.080 0.53 +0.100 0.61 +0.150 0.89 +0.200 1.20 +0.300 1.80 +0.400 2.38 +0.500 2.93 +0.600 3.44 +0.800 4.38 +1 5.20 +1.5 6.90 +2 8.60 +3 11.1 +4 13.4 +5 15.5 +6 17.6 +8 21.6 +10 25.6 \ No newline at end of file diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 42bb958caf..7049f7d314 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -233,7 +233,8 @@ def get_microxs_and_flux( model.export_to_model_xml() comm.barrier() # Reinitialize with tallies - openmc.lib.init(intracomm=comm) + output = run_kwargs.get('output', True) if run_kwargs else True + openmc.lib.init(intracomm=comm, output=output) with TemporaryDirectory() as temp_dir: # Indicate to run in temporary directory unless being executed through diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 6dbb3adb2c..6c19531e7e 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -4,6 +4,7 @@ from contextlib import nullcontext import copy from datetime import datetime import json +from numbers import Integral from pathlib import Path import numpy as np @@ -156,6 +157,7 @@ class R2SManager: mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, operator_kwargs: dict | None = None, + by_parent_nuclide: bool = False, ): """Run the R2S calculation. @@ -179,7 +181,7 @@ class R2SManager: timesteps. For example, if two timesteps are specified, the array of times would contain three entries, and [2] would indicate computing photon results at the last time. A value of None indicates to run - photon transport for each time. + photon transport at each time that has a decay photon source. output_dir : PathLike, optional Path to directory where R2S calculation outputs will be saved. If not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is @@ -207,6 +209,11 @@ class R2SManager: operator_kwargs : dict, optional Additional keyword arguments passed to :class:`openmc.deplete.IndependentOperator`. + by_parent_nuclide : bool, optional + Whether to score photon tallies separately for each parent + radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to + tallies that do not already contain one, with bins determined from + the prepared decay photon sources. Returns ------- @@ -256,9 +263,13 @@ class R2SManager: timesteps, source_rates, timestep_units, output_dir / 'activation', operator_kwargs=operator_kwargs ) - self.step3_photon_transport( + self.step3_photon_source( photon_time_indices, bounding_boxes, output_dir / 'photon_transport', - mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs + mat_vol_kwargs=mat_vol_kwargs, + ) + self.step4_photon_transport( + output_dir / 'photon_transport', run_kwargs=run_kwargs, + by_parent_nuclide=by_parent_nuclide, ) return output_dir @@ -347,7 +358,7 @@ class R2SManager: # Run neutron transport and get fluxes and micros. Run via openmc.lib to # maintain a consistent parallelism strategy with the activation step. - with TemporarySession(): + with TemporarySession(output=False): self.results['fluxes'], self.results['micros'] = get_microxs_and_flux( self.neutron_model, domains, **micro_kwargs) @@ -438,23 +449,20 @@ class R2SManager: # Get depletion results self.results['depletion_results'] = Results(output_path) - def step3_photon_transport( + def step3_photon_source( self, time_indices: Sequence[int] | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, output_dir: PathLike = 'photon_transport', mat_vol_kwargs: dict | None = None, - run_kwargs: dict | None = None, ): - """Run the photon transport step. + """Create decay photon sources. - This step performs photon transport calculations using decay photon - sources created from the activated materials. For each specified time, - it creates appropriate photon sources and runs a transport calculation. - In mesh-based mode, the sources are created using the mesh material - volumes, while in cell-based mode, they are created using bounding boxes - for each cell. This step will populate the 'photon_tallies' key in the - results dictionary. + This step creates decay photon sources from the activated materials for + each specified time. In mesh-based mode, the sources are created using + mesh material volumes, while in cell-based mode, they are created using + bounding boxes for each cell. This step will populate the + 'photon_sources' key in the results dictionary. Parameters ---------- @@ -464,40 +472,54 @@ class R2SManager: timesteps. For example, if two timesteps are specified, the array of times would contain three entries, and [2] would indicate computing photon results at the last time. A value of None indicates to run - photon transport for each time. + photon transport at each time that has a decay photon source. bounding_boxes : dict[int, openmc.BoundingBox], optional Dictionary mapping cell IDs to bounding boxes used for spatial source sampling in cell-based R2S calculations. Required if method is 'cell-based'. output_dir : PathLike, optional - Path to directory where photon transport outputs will be saved. + Path to directory where photon source outputs will be saved. mat_vol_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.MeshBase.material_volumes`. - run_kwargs : dict, optional - Additional keyword arguments passed to :meth:`openmc.Model.run` - during the photon transport step. By default, output is disabled. """ + # Do not retain sources from an earlier successful call if this source + # preparation attempt fails. + self.results.pop('photon_sources', None) + # TODO: Automatically determine bounding box for each cell if bounding_boxes is None and self.method == 'cell-based': raise ValueError("bounding_boxes must be provided for cell-based " "R2S calculations.") - # Set default run arguments if not provided - if run_kwargs is None: - run_kwargs = {} - run_kwargs.setdefault('output', False) - - # Write out JSON file with tally IDs that can be used for loading - # results output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) - # Get default time indices if not provided + # Determine and validate time indices before preparing source data. + n_steps = len(self.results['depletion_results']) + implicit_time_indices = time_indices is None if time_indices is None: - n_steps = len(self.results['depletion_results']) time_indices = list(range(n_steps)) + else: + time_indices = list(time_indices) + if not time_indices: + raise ValueError('time_indices must contain at least one index') + + normalized_indices = [] + for index in time_indices: + if isinstance(index, bool) or not isinstance(index, Integral): + raise TypeError('time_indices must contain only integers') + index = int(index) + if index < -n_steps or index >= n_steps: + raise IndexError( + f'Photon time index {index} is out of range for ' + f'{n_steps} depletion results') + normalized_index = index % n_steps + normalized_indices.append(normalized_index) + + # Remove duplicates while preserving order + time_indices = list(dict.fromkeys(normalized_indices)) # Check whether the photon model is different neutron_univ = self.neutron_model.geometry.root_universe @@ -523,13 +545,6 @@ class R2SManager: self.results['mesh_material_volumes_photon'] = photon_mmv_list - if comm.rank == 0: - tally_ids = [tally.id for tally in self.photon_model.tallies] - with open(output_dir / 'tally_ids.json', 'w') as f: - json.dump(tally_ids, f) - - self.results['photon_tallies'] = {} - # Get dictionary of cells in the photon model if different_photon_model: photon_cells = self.photon_model.geometry.get_all_cells() @@ -547,16 +562,100 @@ class R2SManager: continue work_items.append((cell, original_mat, bounding_boxes[cell.id])) - # Ensure photon transport is enabled in settings + # Create decay photon sources for each time index + photon_sources = { + time_index: self._create_photon_sources(time_index, work_items) + for time_index in time_indices + } + + # Determine if any times have no decay photon sources. If the user + # didn't specify any specific time indices, remove those times from the + # photon_sources dictionary. If the user did specify time indices, raise + # an error if any of those times have no decay photon sources. + empty_indices = [ + time_index for time_index, sources in photon_sources.items() + if not sources + ] + if implicit_time_indices: + for time_index in empty_indices: + del photon_sources[time_index] + if not photon_sources: + raise RuntimeError( + 'No decay photon sources were found at any depletion time') + elif empty_indices: + indices = ', '.join(str(index) for index in empty_indices) + raise RuntimeError( + f'No decay photon source was found for requested time ' + f'indices: {indices}') + + self.results['photon_sources'] = photon_sources + + def step4_photon_transport( + self, + output_dir: PathLike = 'photon_transport', + run_kwargs: dict | None = None, + by_parent_nuclide: bool = False, + ): + """Run photon transport using prepared decay photon sources. + + This step runs a photon transport calculation for each source list + created by :meth:`step3_photon_source`. It will populate the + 'photon_tallies' key in the results dictionary. + + Parameters + ---------- + output_dir : PathLike, optional + Path to directory where photon transport outputs will be saved. + run_kwargs : dict, optional + Additional keyword arguments passed to :meth:`openmc.Model.run`. + By default, output is disabled. + by_parent_nuclide : bool, optional + Whether to score photon tallies separately for each parent + radionuclide. A :class:`~openmc.ParentNuclideFilter` is added to + tallies that do not already contain one, with bins determined from + the prepared decay photon sources. + """ + if 'photon_sources' not in self.results: + raise RuntimeError( + 'Photon sources must be created with step3_photon_source ' + 'before running photon transport.') + photon_sources = self.results['photon_sources'] + if not photon_sources: + raise RuntimeError( + 'No decay photon sources are available for transport') + + if by_parent_nuclide: + radionuclides = sorted({ + nuclide + for sources in photon_sources.values() + for source in sources + for nuclide in source.energy.nuclides + }) + + if radionuclides: + parent_filter = openmc.ParentNuclideFilter(radionuclides) + for tally in self.photon_model.tallies: + if not tally.contains_filter(openmc.ParentNuclideFilter): + tally.filters.append(parent_filter) + + if run_kwargs is None: + run_kwargs = {} + run_kwargs.setdefault('output', False) + + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + if comm.rank == 0: + tally_ids = [tally.id for tally in self.photon_model.tallies] + with open(output_dir / 'tally_ids.json', 'w') as f: + json.dump(tally_ids, f) + + self.results['photon_tallies'] = {} + + # Ensure photon transport is enabled in settings. self.photon_model.settings.photon_transport = True - for time_index in time_indices: - # Convert time_index (which may be negative) to a normal index - if time_index < 0: - time_index += len(self.results['depletion_results']) - - # Build decay photon sources and assign to the photon model - sources = self._create_photon_sources(time_index, work_items) + for time_index, sources in photon_sources.items(): self.photon_model.settings.source = sources # Run photon transport calculation diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index adb0d3dbc0..0cf0141d24 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -1,12 +1,17 @@ +from __future__ import annotations + import numbers import bisect import math from collections.abc import Iterable +from typing import Literal from warnings import warn import h5py import numpy as np +import openmc +from .chain import Chain, _get_chain from .stepresult import StepResult, VERSION_RESULTS import openmc.checkvalue as cv from openmc.data import atomic_mass, AVOGADRO @@ -103,7 +108,8 @@ class Results(list): mat: Material | str, units: str = "Bq/cm3", by_nuclide: bool = False, - volume: float | None = None + volume: float | None = None, + chain_file: Literal[False] | None | PathLike | Chain = None ) -> tuple[np.ndarray, np.ndarray | list[dict]]: """Get activity of material over time. @@ -115,22 +121,31 @@ class Results(list): Material object or material id to evaluate units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'} Specifies the type of activity to return, options include total - activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity [Bq/cm3]. + activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity + [Bq/cm3]. by_nuclide : bool Specifies if the activity should be returned for the material as a whole or per nuclide. Default is False. volume : float, optional Volume of the material. If not passed, defaults to using the :attr:`Material.volume` attribute. + chain_file : False, None, PathLike, or openmc.deplete.Chain, optional + Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is + used. If ``None``, the chain specified by + ``openmc.config['chain_file']`` is used when available. If a path or + :class:`openmc.deplete.Chain` is given, that chain is used. For + ``None`` or an explicit chain, nuclides absent from the chain fall + back to ENDF/B-VIII.0 data. + + .. versionadded:: 0.15.4 Returns ------- times : numpy.ndarray Array of times in [s] activities : numpy.ndarray or List[dict] - Array of total activities if by_nuclide = False (default) - or list of dictionaries of activities by nuclide if - by_nuclide = True. + Array of total activities if by_nuclide = False (default) or list of + dictionaries of activities by nuclide if by_nuclide = True. """ if isinstance(mat, Material): @@ -140,6 +155,13 @@ class Results(list): else: raise TypeError('mat should be of type openmc.Material or str') + if chain_file is not False: + if chain_file is None: + if openmc.config.get('chain_file') is not None: + chain_file = _get_chain(None) + else: + chain_file = _get_chain(chain_file) + times = np.empty_like(self, dtype=float) if by_nuclide: activities = [None] * len(self) @@ -149,7 +171,8 @@ class Results(list): # Evaluate activity for each depletion time for i, result in enumerate(self): times[i] = result.time[0] - activities[i] = result.get_material(mat_id).get_activity(units, by_nuclide, volume) + activities[i] = result.get_material(mat_id).get_activity( + units, by_nuclide, volume, chain_file=chain_file) return times, activities diff --git a/openmc/lib/plot.py b/openmc/lib/plot.py index 44d6ac273d..1966823947 100644 --- a/openmc/lib/plot.py +++ b/openmc/lib/plot.py @@ -87,7 +87,7 @@ _dll.openmc_slice_data.errcheck = _error_handler def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None, - pixels=None, show_overlaps=False, level=-1, filter=None, + pixels=None, show_overlaps=False, level=None, filter=None, include_properties=True): """Generate a 2D raster of geometry and property data for plotting. @@ -111,7 +111,7 @@ def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None, show_overlaps : bool, optional Whether to detect overlapping cells level : int, optional - Universe level (-1 for deepest) + Universe level (None for deepest) filter : openmc.lib.Filter, optional Filter for bin index lookup include_properties : bool, optional @@ -127,6 +127,12 @@ def slice_data(origin, width=None, basis='xy', u_span=None, v_span=None, Array of shape (v_res, h_res, 2) with float64 dtype containing [temperature, density], or None if include_properties=False """ + # Set deepest level as default + if level is None: + level = -1 + if not isinstance(level, int): + raise TypeError("level must be an integer.") + if pixels is None: raise ValueError("pixels must be specified.") if len(pixels) != 2: @@ -255,6 +261,52 @@ def property_map(plot): return prop_data +_dll.openmc_slice_data_overlap_count.argtypes = [POINTER(c_size_t)] +_dll.openmc_slice_data_overlap_count.restype = c_int +_dll.openmc_slice_data_overlap_count.errcheck = _error_handler + +_dll.openmc_slice_data_overlap_info.argtypes = [c_size_t, POINTER(c_int32)] +_dll.openmc_slice_data_overlap_info.restype = c_int +_dll.openmc_slice_data_overlap_info.errcheck = _error_handler + + +# Python wrappings for overlap functions +def slice_data_overlap_count() -> int: + """Return the number of unique overlaps from the last slice plot. + + Returns + ------- + int + Number of unique overlapping cell pairs detected by the most recent + :func:`slice_data` call with overlap checking enabled. + """ + count = c_size_t() + _dll.openmc_slice_data_overlap_count(count) + return count.value + + +def slice_data_overlap_info() -> np.ndarray: + """Return identifying information for overlaps from the last slice plot. + + Returns + ------- + numpy.ndarray + Array of shape ``(n, 3)`` with int32 dtype, where ``n`` is the number + of unique overlaps detected by the most recent :func:`slice_data` call + with overlap checking enabled. Each row contains ``[universe_id, + cell1_id, cell2_id]``. + """ + n = slice_data_overlap_count() + overlap_info = np.empty((n, 3), dtype=np.int32) + + if n > 0: + _dll.openmc_slice_data_overlap_info( + n, + overlap_info.ctypes.data_as(POINTER(c_int32)), + ) + return overlap_info + + _dll.openmc_get_plot_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_plot_index.restype = c_int _dll.openmc_get_plot_index.errcheck = _error_handler diff --git a/openmc/material.py b/openmc/material.py index eded9c3faa..3a92efda22 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -8,7 +8,7 @@ from pathlib import Path import re import sys import tempfile -from typing import Sequence, Dict +from typing import TYPE_CHECKING, Literal, Sequence, Dict import warnings import lxml.etree as ET @@ -28,6 +28,9 @@ from openmc.data.data import _get_element_symbol, JOULE_PER_EV from openmc.data.function import Tabulated1D from openmc.data import mass_energy_absorption_coefficient, dose_coefficients +if TYPE_CHECKING: + from openmc.deplete import Chain + # Units for density supported by OpenMC DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum', @@ -1385,8 +1388,13 @@ class Material(IDManagerMixin): return densities - def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False, - volume: float | None = None) -> dict[str, float] | float: + def get_activity( + self, + units: str = 'Bq/cm3', + by_nuclide: bool = False, + volume: float | None = None, + chain_file: Literal[False] | None | PathLike | Chain = None + ) -> dict[str, float] | float: """Return the activity of the material or each nuclide within. .. versionadded:: 0.13.1 @@ -1405,13 +1413,22 @@ class Material(IDManagerMixin): :attr:`Material.volume` attribute. .. versionadded:: 0.13.3 + chain_file : False, None, PathLike, or openmc.deplete.Chain, optional + Source of half-life values. If ``False``, only ENDF/B-VIII.0 data is + used. If ``None``, the chain specified by + ``openmc.config['chain_file']`` is used when available. If a path or + :class:`openmc.deplete.Chain` is given, that chain is used. For + ``None`` or an explicit chain, nuclides absent from the chain fall + back to ENDF/B-VIII.0 data. + + .. versionadded:: 0.15.4 Returns ------- Union[dict, float] - If by_nuclide is True then a dictionary whose keys are nuclide - names and values are activity is returned. Otherwise the activity - of the material is returned as a float. + If by_nuclide is True then a dictionary whose keys are nuclide names + and values are activity is returned. Otherwise the activity of the + material is returned as a float. """ cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3', 'Ci', 'Ci/m3'}) @@ -1440,7 +1457,8 @@ class Material(IDManagerMixin): activity = {} for nuclide, atoms_per_bcm in self.get_nuclide_atom_densities().items(): - inv_seconds = openmc.data.decay_constant(nuclide) + inv_seconds = openmc.data.decay_constant( + nuclide, chain_file=chain_file) activity[nuclide] = inv_seconds * 1e24 * atoms_per_bcm * multiplier return activity if by_nuclide else sum(activity.values()) diff --git a/openmc/model/model.py b/openmc/model/model.py index 733503fb39..0541138660 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -21,7 +21,8 @@ import openmc import openmc._xml as xml from openmc.dummy_comm import DummyCommunicator from openmc.executor import _process_CLI_arguments -from openmc.checkvalue import check_type, check_value, PathLike +from openmc.checkvalue import (check_type, check_value, check_greater_than, + check_length, PathLike) from openmc.exceptions import InvalidIDError from openmc.plots import add_plot_params, _BASIS_INDICES, id_map_to_rgb from openmc.utility_funcs import change_directory @@ -33,6 +34,15 @@ class ModelModifier(Protocol): ... +def _check_pixels(pixels: int | Sequence[int]) -> None: + if isinstance(pixels, Integral): + check_greater_than('pixels', pixels, 0) + else: + check_length('pixels', pixels, 2) + for p in pixels: + check_greater_than('pixels', p, 0) + + class Model: """Model container. @@ -1050,6 +1060,8 @@ class Model: pixels: int | Sequence[int], basis: str ): + _check_pixels(pixels) + x, y, _ = _BASIS_INDICES[basis] bb = self.bounding_box @@ -1212,6 +1224,8 @@ class Model: """ import openmc.lib + _check_pixels(pixels) + if width is not None and (u_span is not None or v_span is not None): raise ValueError("width is mutually exclusive with u_span/v_span.") @@ -1301,7 +1315,11 @@ class Model: import matplotlib.pyplot as plt check_type('n_samples', n_samples, int | None) + if n_samples is not None: + check_greater_than('n_samples', n_samples, 0, equality=True) check_type('plane_tolerance', plane_tolerance, Real) + check_greater_than('plane_tolerance', plane_tolerance, 0.0) + if legend_kwargs is None: legend_kwargs = {} legend_kwargs.setdefault('bbox_to_anchor', (1.05, 1)) diff --git a/openmc/plots.py b/openmc/plots.py index aeece7acf6..531184095d 100644 --- a/openmc/plots.py +++ b/openmc/plots.py @@ -383,7 +383,7 @@ def id_map_to_rgb( """ # Initialize RGB array with white background (values between 0 and 1 for matplotlib) img = np.ones(id_map.shape, dtype=float) - + # Get the appropriate index based on color_by if color_by == 'cell': id_index = 0 # Cell IDs are in the first channel @@ -391,14 +391,14 @@ def id_map_to_rgb( id_index = 2 # Material IDs are in the third channel else: raise ValueError("color_by must be either 'cell' or 'material'") - + # Get all unique IDs in the plot unique_ids = np.unique(id_map[:, :, id_index]) - + # Generate default colors if not provided if colors is None: colors = {} - + # Convert colors dict to use IDs as keys color_map = {} for key, color in colors.items(): @@ -406,13 +406,13 @@ def id_map_to_rgb( color_map[key.id] = color else: color_map[key] = color - + # Generate random colors for IDs not in color_map rng = np.random.RandomState(1) for uid in unique_ids: if uid > 0 and uid not in color_map: color_map[uid] = rng.randint(0, 256, (3,)) - + # Apply colors to each pixel for uid in unique_ids: if uid == -1: # Background/void @@ -432,7 +432,7 @@ def id_map_to_rgb( rgb = color mask = id_map[:, :, id_index] == uid img[mask] = np.array(rgb) / 255.0 - + return img class PlotBase(IDManagerMixin): diff --git a/openmc/source.py b/openmc/source.py index a11cfd6c8d..5947540d6b 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -1,7 +1,7 @@ from __future__ import annotations from abc import ABC, abstractmethod from collections.abc import Iterable, Sequence -from numbers import Real +from numbers import Integral, Real from pathlib import Path import warnings from typing import Any @@ -46,7 +46,7 @@ class SourceBase(ABC): Attributes ---------- - type : {'independent', 'file', 'compiled', 'mesh'} + type : {'independent', 'file', 'compiled', 'mesh', 'tokamak'} Indicator of source type. strength : float Strength of the source @@ -203,6 +203,8 @@ class SourceBase(ABC): return FileSource.from_xml_element(elem) elif source_type == 'mesh': return MeshSource.from_xml_element(elem, meshes) + elif source_type == 'tokamak': + return TokamakSource.from_xml_element(elem) else: raise ValueError( f'Source type {source_type} is not recognized') @@ -896,6 +898,430 @@ class FileSource(SourceBase): return cls(**kwargs) +class TokamakSource(SourceBase): + r"""A source representing neutron emission from a tokamak plasma. + + This source samples neutron positions from a tokamak plasma geometry using + Miller-style flux surface parameterization. The user provides an emission + profile S(r/a) as a function of normalized minor radius, along with one or + more energy distributions. + + The flux surface parameterization is + + .. math:: + + \begin{aligned} + R &= R_0 + r \cos\left(\alpha + \delta \sin\alpha\right) + + \Delta \left[1 - \left(\frac{r}{a}\right)^2\right] \\ + Z &= Z_\mathrm{shift} + \kappa r \sin\alpha + \end{aligned} + + where :math:`R_0` is major radius, :math:`a` is minor radius, + :math:`\kappa` is elongation, :math:`\delta` is triangularity, + :math:`\Delta` is the Shafranov shift, and :math:`Z_\mathrm{shift}` is + the vertical shift. + + .. versionadded:: 0.15.4 + + Parameters + ---------- + major_radius : float + Major radius R0 in [cm] + minor_radius : float + Minor radius a in [cm] + elongation : float + Plasma elongation κ (must be > 0) + triangularity : float + Plasma triangularity δ (must be in [-1, 1]) + shafranov_shift : float + Shafranov shift Δ in [cm] (must be >= 0 and < a/2) + r_over_a : numpy.ndarray + Normalized minor radius grid points, must start at 0 and end at 1 + emission_density : numpy.ndarray + Emission density S(r) at each r/a point (arbitrary units, must be >= 0). + Values are linearly interpolated between grid points and refined on an + internal grid for radial sampling. Must have the same length as + ``r_over_a`` and contain at least one positive value. + energy : openmc.stats.Univariate or Sequence[openmc.stats.Univariate] + Energy distribution(s). Either a single distribution used at all radii, + or one distribution per ``r_over_a`` grid point. When one distribution + per grid point is given, the energy of a sampled particle is drawn from + one of the two distributions bracketing its sampled radius, selected + stochastically with probability proportional to the proximity of the + radius to each grid point (stochastic interpolation). + time : openmc.stats.Univariate, optional + Time distribution of the source. If None, particles are born at + :math:`t=0`, matching the default behavior of + :class:`openmc.IndependentSource`. + phi_start : float + Starting toroidal angle in [rad] (default: 0) + phi_extent : float + Toroidal angle extent in [rad] (default: 2π) + n_alpha : int + Number of poloidal angle grid points for CDF sampling (default: 101) + vertical_shift : float + Vertical shift of the plasma center in [cm] (default: 0) + strength : float + Strength of the source (default: 1.0) + constraints : dict + Constraints on sampled source particles. See :class:`SourceBase` for + valid keys and values. + + Attributes + ---------- + major_radius : float + Major radius R0 in [cm] + minor_radius : float + Minor radius a in [cm] + elongation : float + Plasma elongation κ + triangularity : float + Plasma triangularity δ + shafranov_shift : float + Shafranov shift Δ in [cm] + r_over_a : numpy.ndarray + Normalized minor radius grid points + emission_density : numpy.ndarray + Emission density S(r) at each r/a point + energy : list of openmc.stats.Univariate + Energy distribution(s) + time : openmc.stats.Univariate or None + Time distribution of the source + phi_start : float + Starting toroidal angle in [rad] + phi_extent : float + Toroidal angle extent in [rad] + n_alpha : int + Number of poloidal angle grid points + vertical_shift : float + Vertical shift of the plasma center in [cm] + strength : float + Strength of the source + type : str + Indicator of source type: 'tokamak' + constraints : dict + Constraints on sampled source particles + + """ + + def __init__( + self, + major_radius: float, + minor_radius: float, + elongation: float, + triangularity: float, + shafranov_shift: float, + r_over_a: Sequence[float], + emission_density: Sequence[float], + energy: Univariate | Sequence[Univariate], + time: Univariate | None = None, + phi_start: float = 0.0, + phi_extent: float = 2.0 * np.pi, + n_alpha: int = 101, + vertical_shift: float = 0.0, + strength: float = 1.0, + constraints: dict[str, Any] | None = None + ): + super().__init__(strength=strength, constraints=constraints) + self.major_radius = major_radius + self.minor_radius = minor_radius + self.elongation = elongation + self.triangularity = triangularity + self.shafranov_shift = shafranov_shift + self.r_over_a = r_over_a + self.emission_density = emission_density + self.phi_start = phi_start + self.phi_extent = phi_extent + self.n_alpha = n_alpha + self.vertical_shift = vertical_shift + self.energy = energy + self.time = time + + self._validate() + + def _validate(self): + """Validate relationships between tokamak source parameters.""" + if self.minor_radius >= self.major_radius: + raise ValueError( + f"minor_radius ({self.minor_radius}) must be smaller than " + f"major_radius ({self.major_radius})") + if self.shafranov_shift >= 0.5 * self.minor_radius: + raise ValueError( + f"shafranov_shift ({self.shafranov_shift}) must be smaller " + f"than half the minor_radius ({0.5 * self.minor_radius})") + if len(self.emission_density) != len(self.r_over_a): + raise ValueError( + f"emission_density (length {len(self.emission_density)}) must " + f"have the same length as r_over_a (length {len(self.r_over_a)})") + if not np.any(self.emission_density > 0.0): + raise ValueError("emission_density must contain a positive value") + if len(self.energy) not in (1, len(self.r_over_a)): + raise ValueError( + f"Number of energy distributions ({len(self.energy)}) must be " + f"either 1 or equal to the number of r_over_a grid points " + f"({len(self.r_over_a)})") + + @property + def type(self) -> str: + return "tokamak" + + @property + def major_radius(self) -> float: + return self._major_radius + + @major_radius.setter + def major_radius(self, value: float): + cv.check_type('major radius', value, Real) + cv.check_greater_than('major radius', value, 0.0) + self._major_radius = value + + @property + def minor_radius(self) -> float: + return self._minor_radius + + @minor_radius.setter + def minor_radius(self, value: float): + cv.check_type('minor radius', value, Real) + cv.check_greater_than('minor radius', value, 0.0) + self._minor_radius = value + + @property + def elongation(self) -> float: + return self._elongation + + @elongation.setter + def elongation(self, value: float): + cv.check_type('elongation', value, Real) + cv.check_greater_than('elongation', value, 0.0) + self._elongation = value + + @property + def triangularity(self) -> float: + return self._triangularity + + @triangularity.setter + def triangularity(self, value: float): + cv.check_type('triangularity', value, Real) + cv.check_greater_than('triangularity', value, -1.0, equality=True) + cv.check_less_than('triangularity', value, 1.0, equality=True) + self._triangularity = value + + @property + def shafranov_shift(self) -> float: + return self._shafranov_shift + + @shafranov_shift.setter + def shafranov_shift(self, value: float): + cv.check_type('Shafranov shift', value, Real) + cv.check_greater_than('Shafranov shift', value, 0.0, equality=True) + self._shafranov_shift = value + + @property + def r_over_a(self) -> np.ndarray: + return self._r_over_a + + @r_over_a.setter + def r_over_a(self, value: Sequence[float]): + value = np.asarray(value, dtype=float) + if value.ndim != 1 or len(value) < 2: + raise ValueError("r_over_a must be a 1-D array with at least 2 points") + if value[0] != 0.0: + raise ValueError("r_over_a must start at 0") + if value[-1] != 1.0: + raise ValueError("r_over_a must end at 1") + if not np.all(np.diff(value) > 0): + raise ValueError("r_over_a must be strictly increasing") + self._r_over_a = value + + @property + def emission_density(self) -> np.ndarray: + return self._emission_density + + @emission_density.setter + def emission_density(self, value: Sequence[float]): + value = np.asarray(value, dtype=float) + if value.ndim != 1: + raise ValueError("emission_density must be a 1-D array") + if np.any(value < 0): + raise ValueError("emission_density values cannot be negative") + self._emission_density = value + + @property + def energy(self) -> list[Univariate]: + return self._energy + + @energy.setter + def energy(self, value: Univariate | Sequence[Univariate]): + if isinstance(value, Univariate): + self._energy = [value] + else: + cv.check_iterable_type('energy distributions', value, Univariate) + self._energy = list(value) + + @property + def time(self) -> Univariate | None: + return self._time + + @time.setter + def time(self, value: Univariate | None): + if value is not None: + cv.check_type('time distribution', value, Univariate) + self._time = value + + @property + def phi_start(self) -> float: + return self._phi_start + + @phi_start.setter + def phi_start(self, value: float): + cv.check_type('phi_start', value, Real) + self._phi_start = value + + @property + def phi_extent(self) -> float: + return self._phi_extent + + @phi_extent.setter + def phi_extent(self, value: float): + cv.check_type('phi_extent', value, Real) + cv.check_greater_than('phi_extent', value, 0.0) + cv.check_less_than('phi_extent', value, 2.0 * np.pi, equality=True) + self._phi_extent = value + + @property + def n_alpha(self) -> int: + return self._n_alpha + + @n_alpha.setter + def n_alpha(self, value: int): + cv.check_type('n_alpha', value, Integral) + cv.check_greater_than('n_alpha', value, 2) + if value < 51: + warnings.warn( + "n_alpha values below 51 may introduce noticeable " + "discretization bias in tokamak source sampling", stacklevel=2) + self._n_alpha = value + + @property + def vertical_shift(self) -> float: + return self._vertical_shift + + @vertical_shift.setter + def vertical_shift(self, value: float): + cv.check_type('vertical shift', value, Real) + self._vertical_shift = value + + def populate_xml_element(self, element): + """Add necessary tokamak source information to an XML element + + Returns + ------- + element : lxml.etree._Element + XML element containing source data + + """ + self._validate() + + # Geometry parameters + ET.SubElement(element, "major_radius").text = str(self.major_radius) + ET.SubElement(element, "minor_radius").text = str(self.minor_radius) + ET.SubElement(element, "elongation").text = str(self.elongation) + ET.SubElement(element, "triangularity").text = str(self.triangularity) + ET.SubElement(element, "shafranov_shift").text = str(self.shafranov_shift) + + # Toroidal angle bounds + ET.SubElement(element, "phi_start").text = str(self.phi_start) + ET.SubElement(element, "phi_extent").text = str(self.phi_extent) + + # Poloidal sampling resolution + ET.SubElement(element, "n_alpha").text = str(self.n_alpha) + + # Vertical shift + if self.vertical_shift != 0.0: + ET.SubElement(element, "vertical_shift").text = str(self.vertical_shift) + + # Emission profile + ET.SubElement(element, "r_over_a").text = ' '.join(str(r) for r in self.r_over_a) + ET.SubElement(element, "emission_density").text = ' '.join(str(s) for s in self.emission_density) + + # Energy distribution(s) + for dist in self.energy: + element.append(dist.to_xml_element('energy')) + + # Time distribution + if self.time is not None: + element.append(self.time.to_xml_element('time')) + + @classmethod + def from_xml_element(cls, elem: ET.Element) -> TokamakSource: + """Generate tokamak source from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + + Returns + ------- + openmc.TokamakSource + Source generated from XML element + + """ + # Read geometry parameters + major_radius = float(get_text(elem, 'major_radius')) + minor_radius = float(get_text(elem, 'minor_radius')) + elongation = float(get_text(elem, 'elongation')) + triangularity = float(get_text(elem, 'triangularity')) + shafranov_shift = float(get_text(elem, 'shafranov_shift')) + + # Read optional parameters + phi_start_text = get_text(elem, 'phi_start') + phi_start = float(phi_start_text) if phi_start_text else 0.0 + + phi_extent_text = get_text(elem, 'phi_extent') + phi_extent = float(phi_extent_text) if phi_extent_text else 2.0 * np.pi + + n_alpha_text = get_text(elem, 'n_alpha') + n_alpha = int(n_alpha_text) if n_alpha_text else 101 + + vertical_shift_text = get_text(elem, 'vertical_shift') + vertical_shift = float(vertical_shift_text) if vertical_shift_text else 0.0 + + # Read emission profile + r_over_a = np.array([float(x) for x in get_text(elem, 'r_over_a').split()]) + emission_density = np.array([float(x) for x in get_text(elem, 'emission_density').split()]) + + # Read energy distributions + energy = [Univariate.from_xml_element(e) for e in elem.findall('energy')] + if len(energy) == 1: + energy = energy[0] + + # Read time distribution + time_elem = elem.find('time') + time = Univariate.from_xml_element(time_elem) if time_elem is not None else None + + # Read constraints and strength + constraints = cls._get_constraints(elem) + strength_text = get_text(elem, 'strength') + strength = float(strength_text) if strength_text else 1.0 + + return cls( + major_radius=major_radius, + minor_radius=minor_radius, + elongation=elongation, + triangularity=triangularity, + shafranov_shift=shafranov_shift, + r_over_a=r_over_a, + emission_density=emission_density, + energy=energy, + time=time, + phi_start=phi_start, + phi_extent=phi_extent, + n_alpha=n_alpha, + vertical_shift=vertical_shift, + strength=strength, + constraints=constraints + ) class SourceParticle: diff --git a/pyproject.toml b/pyproject.toml index 098487c06e..c769fb7d76 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,6 +40,7 @@ dependencies = [ [project.optional-dependencies] depletion-mpi = ["mpi4py"] docs = [ + "breathe", "sphinx", "sphinxcontrib-katex", "sphinx-numfig", diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index 5bbda48305..d166f12479 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -38,6 +38,9 @@ void VacuumBC::handle_particle(Particle& p, const Surface& surf) const void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const { Direction u = surf.reflect(p.r(), p.u(), &p); + + // normalize reflected u to ensure no floating point error leads to + // unnormalized directions u /= u.norm(); // Handle the effects of the surface albedo on the particle's weight. @@ -53,6 +56,9 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const void WhiteBC::handle_particle(Particle& p, const Surface& surf) const { Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed()); + + // normalize outgoing u to ensure no floating point error leads to + // unnormalized directions u /= u.norm(); // Handle the effects of the surface albedo on the particle's weight. @@ -241,6 +247,10 @@ void RotationalPeriodicBC::handle_particle( new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_]; new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_]; + // normalize new_u to ensure no floating point error leads to unnormalized + // directions + new_u /= new_u.norm(); + // Handle the effects of the surface albedo on the particle's weight. BoundaryCondition::handle_albedo(p, surf); diff --git a/src/distribution.cpp b/src/distribution.cpp index 7f5b498add..cbdd131273 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -552,14 +552,9 @@ double Tabular::sample_unbiased(uint64_t* seed) const double c = prn(seed); // Find first CDF bin which is above the sampled value - double c_i = c_[0]; - int i; - std::size_t n = c_.size(); - for (i = 0; i < n - 1; ++i) { - if (c <= c_[i + 1]) - break; - c_i = c_[i + 1]; - } + auto c_iter = std::lower_bound(c_.begin() + 1, c_.end(), c); + int i = std::distance(c_.begin(), c_iter) - 1; + double c_i = c_[i]; // Determine bounding PDF values double x_i = x_[i]; diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index ba8658f10e..5a34524222 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -262,9 +262,11 @@ std::pair SphericalIndependent::sample(uint64_t* seed) const MeshSpatial::MeshSpatial(pugi::xml_node node) { - if (get_node_value(node, "type", true, true) != "mesh") { - fatal_error(fmt::format( - "Incorrect spatial type '{}' for a MeshSpatial distribution")); + auto spatial_type = get_node_value(node, "type", true, true); + if (spatial_type != "mesh") { + fatal_error( + fmt::format("Incorrect spatial type '{}' for a MeshSpatial distribution", + spatial_type)); } // No in-tet distributions implemented, could include distributions for the diff --git a/src/geometry.cpp b/src/geometry.cpp index ddb61385f1..bf654f4f9b 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -26,16 +26,22 @@ int n_coord_levels; vector overlap_check_count; +vector overlap_keys; +std::unordered_map overlap_key_index; + } // namespace model //============================================================================== // Non-member functions //============================================================================== -bool check_cell_overlap(GeometryState& p, bool error) +int check_cell_overlap(GeometryState& p, bool error) { int n_coord = p.n_coord(); + // If no overlap found, return a nonphysical index + int overlap_index = -1; + // Loop through each coordinate level for (int j = 0; j < n_coord; j++) { Universe& univ = *model::universes[p.coord(j).universe()]; @@ -44,21 +50,40 @@ bool check_cell_overlap(GeometryState& p, bool error) for (auto index_cell : univ.cells_) { Cell& c = *model::cells[index_cell]; if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) { +#pragma omp atomic + ++model::overlap_check_count[index_cell]; if (index_cell != p.coord(j).cell()) { if (error) { fatal_error( fmt::format("Overlapping cells detected: {}, {} on universe {}", c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_)); } - return true; + + // With no fatal error (plotter is calling), now adds overlaps and + // ensures order does not matter when making overlap key + int cell_a = model::cells[index_cell]->id_; + int cell_b = model::cells[p.coord(j).cell()]->id_; + int a = std::min(cell_a, cell_b); + int b = std::max(cell_a, cell_b); + OverlapKey key {univ.id_, a, b}; +#pragma omp critical(overlap_key_update) + { + auto it = model::overlap_key_index.find(key); + if (it != model::overlap_key_index.end()) { + overlap_index = it->second; // already exists, reuse index + } else { + int idx = int(model::overlap_keys.size()); + model::overlap_keys.push_back(key); + model::overlap_key_index[key] = idx; + overlap_index = idx; + } + } + break; } -#pragma omp atomic - ++model::overlap_check_count[index_cell]; } } } - - return false; + return overlap_index; } //============================================================================== diff --git a/src/lattice.cpp b/src/lattice.cpp index e799a340e6..9e56c58e8e 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -10,6 +10,7 @@ #include "openmc/geometry.h" #include "openmc/geometry_aux.h" #include "openmc/hdf5_interface.h" +#include "openmc/math_functions.h" #include "openmc/string_utils.h" #include "openmc/vector.h" #include "openmc/xml_interface.h" @@ -260,25 +261,33 @@ std::pair> RectLattice::distance( // Determine the oncoming edge. double x0 {copysign(0.5 * pitch_[0], u.x)}; double y0 {copysign(0.5 * pitch_[1], u.y)}; - double z0; - double d = std::min( - u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY); + // Evaluate the distance to each oncoming edge independently. Comparing these + // distances directly (rather than reconstructing the crossing position) + // avoids the floating-point cancellation that occurs for large pitches. + double dx = u.x != 0.0 ? (x0 - x) / u.x : INFTY; + double dy = u.y != 0.0 ? (y0 - y) / u.y : INFTY; + double dz = INFTY; if (is_3d_) { - z0 = copysign(0.5 * pitch_[2], u.z); - d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY); + double z0 {copysign(0.5 * pitch_[2], u.z)}; + dz = u.z != 0.0 ? (z0 - z) / u.z : INFTY; } - // Determine which lattice boundaries are being crossed + // The distance to the nearest lattice boundary is the smallest axial + // distance. + double d = std::min({dx, dy, dz}); + + // Determine which lattice boundaries are being crossed. The axis attaining + // the minimum is exactly equal to d, so at least one translation is always + // set for a finite crossing; a near-equal second axis indicates a corner + // crossing. array lattice_trans = {0, 0, 0}; - if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION) + if (isclose(d, dx, FP_COINCIDENT, FP_PRECISION)) lattice_trans[0] = copysign(1, u.x); - if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION) + if (isclose(d, dy, FP_COINCIDENT, FP_PRECISION)) lattice_trans[1] = copysign(1, u.y); - if (is_3d_) { - if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION) - lattice_trans[2] = copysign(1, u.z); - } + if (is_3d_ && isclose(d, dz, FP_COINCIDENT, FP_PRECISION)) + lattice_trans[2] = copysign(1, u.z); return {d, lattice_trans}; } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 81f08fa04e..ddacc2bd9b 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -1,5 +1,7 @@ #include "openmc/math_functions.h" +#include // for numeric_limits + #include "openmc/external/Faddeeva.hh" #include "openmc/constants.h" @@ -944,6 +946,46 @@ double log1prel(double x) } } +double cyl_bessel_j(int n, double x) +{ + // Handle negative arguments via the parity relation + // J_n(-x) = (-1)^n J_n(x); std::cyl_bessel_j has a domain error for x < 0. + double sign = 1.0; + if (x < 0.0) { + x = -x; + if (n % 2 == 1) + sign = -1.0; + } + +#if defined(__cpp_lib_math_special_functions) && \ + __cpp_lib_math_special_functions >= 201603L + return sign * std::cyl_bessel_j(static_cast(n), x); +#else + // Ascending power series (e.g., Abramowitz & Stegun eq. 9.1.10): + // J_n(x) = sum_{m=0}^inf (-1)^m / (m! (m+n)!) * (x/2)^(2m+n) + // The term ratio is -(x/2)^2 / (m*(m+n)), so for |x| <= 2 the series + // converges to machine precision within ~20 terms. + double half_x = 0.5 * x; + + // First term: (x/2)^n / n! + double term = 1.0; + for (int k = 1; k <= n; ++k) { + term *= half_x / k; + } + + double sum = term; + double neg_half_x_sq = -half_x * half_x; + for (int m = 1; m <= 50; ++m) { + term *= neg_half_x_sq / (m * (m + n)); + sum += term; + if (std::abs(term) <= + std::numeric_limits::epsilon() * std::abs(sum)) + break; + } + return sign * sum; +#endif +} + // Helper function to get index and interpolation function on an incident energy // grid void get_energy_index( @@ -959,4 +1001,12 @@ void get_energy_index( } } +// Return true if two floating-point values are approximately equal within a +// combined relative and absolute tolerance. +bool isclose(double a, double b, double rel_tol, double abs_tol) +{ + return std::abs(a - b) <= + std::max(rel_tol * std::max(std::abs(a), std::abs(b)), abs_tol); +} + } // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index 779ae18da9..10e8f213dd 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -338,13 +338,14 @@ void Particle::event_cross_surface() boundary().lattice_translation()[2] != 0) { // Particle crosses lattice boundary + int i_lattice = coord(boundary().coord_level() - 1).lattice(); bool verbose = settings::verbosity >= 10 || trace(); cross_lattice(*this, boundary(), verbose); event() = TallyEvent::LATTICE; // Score cell to cell partial currents if (!model::active_surface_tallies.empty()) { - auto& lat {*model::lattices[lowest_coord().lattice()]}; + auto& lat {*model::lattices[i_lattice]}; bool is_valid; Direction normal = lat.get_normal(boundary().lattice_translation(), is_valid); @@ -553,15 +554,30 @@ void Particle::event_death() finalize_particle_track(*this); } -// Contribute tally reduction variables to global accumulator + // Contribute tally reduction variables to global accumulator + const auto k_absorption = keff_tally_absorption(); + const auto k_collision = keff_tally_collision(); + const auto k_tracklength = keff_tally_tracklength(); + const auto leakage = keff_tally_leakage(); + + if (settings::run_mode == RunMode::EIGENVALUE) { + if (k_absorption != 0.0) { #pragma omp atomic - global_tally_absorption += keff_tally_absorption(); + global_tally_absorption += k_absorption; + } + if (k_collision != 0.0) { #pragma omp atomic - global_tally_collision += keff_tally_collision(); + global_tally_collision += k_collision; + } + if (k_tracklength != 0.0) { #pragma omp atomic - global_tally_tracklength += keff_tally_tracklength(); + global_tally_tracklength += k_tracklength; + } + } + if (leakage != 0.0) { #pragma omp atomic - global_tally_leakage += keff_tally_leakage(); + global_tally_leakage += leakage; + } // Reset particle tallies once accumulated keff_tally_absorption() = 0.0; diff --git a/src/physics.cpp b/src/physics.cpp index dcc7713292..4a8b4d368e 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -359,8 +359,8 @@ void sample_photon_reaction(Particle& p) // Allow electrons to fill orbital and produce Auger electrons and // fluorescent photons. Since Compton subshell data does not match atomic // relaxation data, use the mapping between the data to find the subshell - if (settings::atomic_relaxation && i_shell >= 0 && - element.subshell_map_[i_shell] >= 0) { + if (settings::atomic_relaxation && element.has_atomic_relaxation_ && + i_shell >= 0 && element.subshell_map_[i_shell] >= 0) { element.atomic_relaxation(element.subshell_map_[i_shell], p); } diff --git a/src/plot.cpp b/src/plot.cpp index b9a1136afd..96df1c5beb 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -73,7 +73,7 @@ void IdData::set_value(size_t y, size_t x, const Particle& p, int level, } } -void IdData::set_overlap(size_t y, size_t x) +void IdData::set_overlap(size_t y, size_t x, int /*overlap_idx*/) { for (size_t k = 0; k < data_.shape(2); ++k) data_(y, x, k) = OVERLAP; @@ -88,13 +88,10 @@ void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level, { Cell* c = model::cells.at(p.lowest_coord().cell()).get(); data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; - if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) { - Material* m = model::materials.at(p.material()).get(); - data_(y, x, 1) = m->density_gpcc_; - } + data_(y, x, 1) = c->density(p.cell_instance()); } -void PropertyData::set_overlap(size_t y, size_t x) +void PropertyData::set_overlap(size_t y, size_t x, int /*overlap_idx*/) { data_(y, x) = OVERLAP; } @@ -150,18 +147,18 @@ void RasterData::set_value(size_t y, size_t x, const Particle& p, int level, // set density (g/cm³) if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) { Material* m = model::materials.at(p.material()).get(); - property_data_(y, x, 1) = m->density_gpcc_; + property_data_(y, x, 1) = c->density(p.cell_instance()); } } -void RasterData::set_overlap(size_t y, size_t x) +void RasterData::set_overlap(size_t y, size_t x, int overlap_idx) { - // Set cell, instance, and material to OVERLAP, but preserve filter bin - id_data_(y, x, 0) = OVERLAP; + // Set cell, instance, and material to OVERLAP, but preserve filter bin for + // tally plotting. Cell encodes the overlap index as a negative number so that + // it can be used to look up overlap information in the plotter. + id_data_(y, x, 0) = OVERLAP - overlap_idx - 1; id_data_(y, x, 1) = OVERLAP; id_data_(y, x, 2) = OVERLAP; - // Note: id_data_(y, x, 3) is NOT overwritten - preserves filter bin for tally - // plotting property_data_(y, x, 0) = OVERLAP; property_data_(y, x, 1) = OVERLAP; @@ -1991,10 +1988,12 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3], plot_params.show_overlaps_ = color_overlaps; plot_params.slice_level_ = level; + // Clear overlap data structures on new slice call + model::overlap_keys.clear(); + model::overlap_key_index.clear(); + // Use get_map to generate data auto data = plot_params.get_map(filter_index); - - // Copy geometry data std::copy(data.id_data_.begin(), data.id_data_.end(), geom_data); // Copy property data if requested @@ -2010,6 +2009,32 @@ extern "C" int openmc_slice_data(const double origin[3], const double u_span[3], return 0; } +// Gets the number of overlaps that we need data for +extern "C" int openmc_slice_data_overlap_count(size_t* count) +{ + if (!count) { + set_errmsg("Null pointer passed for overlap count."); + return OPENMC_E_INVALID_ARGUMENT; + } + *count = model::overlap_keys.size(); + + return 0; +} + +// Plotter pre-allocates array size based on what is returned with +// overlap_count; populates an array of size 3*count +extern "C" int openmc_slice_data_overlap_info( + size_t count, int32_t* overlap_info) +{ + for (size_t i = 0; i < count; ++i) { + overlap_info[i * 3] = model::overlap_keys[i].universe_id; + overlap_info[i * 3 + 1] = model::overlap_keys[i].cell1_id; + overlap_info[i * 3 + 2] = model::overlap_keys[i].cell2_id; + } + + return 0; +} + extern "C" int openmc_get_plot_index(int32_t id, int32_t* index) { auto it = model::plot_map.find(id); diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index 29457569b9..f2a81fc1c7 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -12,6 +12,45 @@ constexpr uint64_t prn_mult {6364136223846793005ULL}; // multiplication constexpr uint64_t prn_add {1442695040888963407ULL}; // additive factor, c uint64_t prn_stride {DEFAULT_STRIDE}; // stride between particles +namespace { + +struct SkipAheadCoefficients { + uint64_t multiplier; + uint64_t increment; +}; + +SkipAheadCoefficients future_seed_coefficients(uint64_t n) +{ + // The algorithm here to determine the parameters used to skip ahead is + // described in F. Brown, "Random Number Generation with Arbitrary Stride," + // Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in + // O(log2(N)) operations instead of O(N). Basically, it computes parameters G + // and C which can then be used to find x_N = G*x_0 + C mod 2^M. + + // Initialize constants + uint64_t g {prn_mult}; + uint64_t c {prn_add}; + uint64_t g_new {1}; + uint64_t c_new {0}; + + while (n > 0) { + // Check if the least significant bit is 1. + if (n & 1) { + g_new *= g; + c_new = c_new * g + c; + } + c *= (g + 1); + g *= g; + + // Move bits right, dropping least significant bit. + n >>= 1; + } + + return {g_new, c_new}; +} + +} // namespace + //============================================================================== // PRN //============================================================================== @@ -69,9 +108,10 @@ uint64_t init_seed(int64_t id, int offset) void init_particle_seeds(int64_t id, uint64_t* seeds) { + auto [multiplier, increment] = + future_seed_coefficients(static_cast(id) * prn_stride); for (int i = 0; i < N_STREAMS; i++) { - seeds[i] = - future_seed(static_cast(id) * prn_stride, master_seed + i); + seeds[i] = multiplier * (master_seed + i) + increment; } } @@ -90,33 +130,9 @@ void advance_prn_seed(int64_t n, uint64_t* seed) uint64_t future_seed(uint64_t n, uint64_t seed) { - // The algorithm here to determine the parameters used to skip ahead is - // described in F. Brown, "Random Number Generation with Arbitrary Stride," - // Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in - // O(log2(N)) operations instead of O(N). Basically, it computes parameters G - // and C which can then be used to find x_N = G*x_0 + C mod 2^M. - - // Initialize constants - uint64_t g {prn_mult}; - uint64_t c {prn_add}; - uint64_t g_new {1}; - uint64_t c_new {0}; - - while (n > 0) { - // Check if the least significant bit is 1. - if (n & 1) { - g_new *= g; - c_new = c_new * g + c; - } - c *= (g + 1); - g *= g; - - // Move bits right, dropping least significant bit. - n >>= 1; - } - // With G and C, we can now find the new seed. - return g_new * seed + c_new; + auto [multiplier, increment] = future_seed_coefficients(n); + return multiplier * seed + increment; } //============================================================================== diff --git a/src/simulation.cpp b/src/simulation.cpp index 900a383ab3..03f40a726e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -12,6 +12,7 @@ #include "openmc/material.h" #include "openmc/message_passing.h" #include "openmc/nuclide.h" +#include "openmc/openmp_interface.h" #include "openmc/output.h" #include "openmc/particle.h" #include "openmc/photon.h" @@ -41,6 +42,7 @@ #include #include +#include #include //============================================================================== @@ -354,6 +356,94 @@ int64_t simulation_tracks_completed {0}; } // namespace simulation +namespace { + +//! Collect thread-local secondary banks into the shared secondary bank in +//! sorted order. +//! +//! \param thread_banks Secondary banks produced by each OpenMP thread +void collect_sorted_history_secondary_banks( + vector>& thread_banks) +{ + // Count the total number of all secondary sites produced + int64_t n_collected = 0; + for (const auto& bank : thread_banks) { + n_collected += bank.size(); + } + + // Count the expected number of progeny from per-parent progeny counts + int64_t n_progeny = 0; + for (int64_t count : simulation::progeny_per_particle) { + n_progeny += count; + } + + if (n_collected != n_progeny) { + fatal_error("Mismatch detected between sum of all particle progeny and " + "secondary bank size during collection."); + } + + // Convert per-parent progeny counts to offsets into the sorted bank + std::exclusive_scan(simulation::progeny_per_particle.begin(), + simulation::progeny_per_particle.end(), + simulation::progeny_per_particle.begin(), 0); + + // Allocate the shared bank once for the complete generation + simulation::shared_secondary_bank_write.resize(0); + simulation::shared_secondary_bank_write.extend_uninitialized(n_progeny); + + // Place each secondary according to its parent and progeny identifiers + for (const auto& bank : thread_banks) { + for (const auto& site : bank) { + if (site.parent_id < 0 || + site.parent_id >= + static_cast(simulation::progeny_per_particle.size())) { + fatal_error(fmt::format("Invalid parent_id {} for banked site " + "(expected range [0, {})).", + site.parent_id, simulation::progeny_per_particle.size())); + } + int64_t idx = + simulation::progeny_per_particle[site.parent_id] + site.progeny_id; + if (idx < 0 || idx >= n_progeny) { + fatal_error("Mismatch detected between sum of all particle progeny and " + "secondary bank size during collection."); + } + simulation::shared_secondary_bank_write[idx] = site; + } + } +} + +//! Collect particle-local secondary banks into the shared secondary bank. +//! +//! \param n_particles Number of particles in the active event-based buffer +void collect_event_secondary_banks(int64_t n_particles) +{ + // Compute offsets for each particle's local secondary bank. + vector offsets(n_particles); + int64_t total = 0; + for (int64_t i = 0; i < n_particles; ++i) { + offsets[i] = total; + total += simulation::particles[i].local_secondary_bank().size(); + } + + // Extend the shared bank once for all collected secondaries + int64_t bank_offset = + simulation::shared_secondary_bank_write.extend_uninitialized(total); + + // Copy each local bank into its assigned range and clear the local storage +#pragma omp parallel for schedule(static) + for (int64_t i = 0; i < n_particles; ++i) { + auto& local_bank = simulation::particles[i].local_secondary_bank(); + if (!local_bank.empty()) { + std::copy(local_bank.cbegin(), local_bank.cend(), + simulation::shared_secondary_bank_write.data() + bank_offset + + offsets[i]); + local_bank.clear(); + } + } +} + +} // namespace + //============================================================================== // Non-member functions //============================================================================== @@ -886,11 +976,14 @@ void transport_history_based_single_particle(Particle& p) void transport_history_based() { -#pragma omp parallel for schedule(runtime) - for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) { +#pragma omp parallel + { Particle p; - initialize_particle_track(p, i_work, false); - transport_history_based_single_particle(p); +#pragma omp for schedule(runtime) + for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) { + initialize_particle_track(p, i_work, false); + transport_history_based_single_particle(p); + } } } @@ -918,30 +1011,27 @@ void transport_history_based_shared_secondary() std::fill(simulation::progeny_per_particle.begin(), simulation::progeny_per_particle.end(), 0); + vector> thread_banks(num_threads()); + // Phase 1: Transport primary particles and deposit first generation of // secondaries in the shared secondary bank #pragma omp parallel { - vector thread_bank; + auto& thread_bank = thread_banks[thread_num()]; + Particle p; #pragma omp for schedule(runtime) for (int64_t i = 1; i <= simulation::work_per_rank; i++) { - Particle p; initialize_particle_track(p, i, false); transport_history_based_single_particle(p); for (auto& site : p.local_secondary_bank()) { thread_bank.push_back(site); } - } - - // Drain thread-local bank into the shared secondary bank (once per thread) -#pragma omp critical(SharedSecondaryBank) - { - for (auto& site : thread_bank) { - simulation::shared_secondary_bank_write.thread_unsafe_append(site); - } + p.local_secondary_bank().clear(); } } + collect_sorted_history_secondary_banks(thread_banks); + thread_banks.clear(); simulation::simulation_tracks_completed += settings::n_particles; @@ -951,10 +1041,6 @@ void transport_history_based_shared_secondary() int64_t alive_secondary = 1; while (alive_secondary) { - // Sort the shared secondary bank by parent ID then progeny ID to - // ensure reproducibility. - sort_bank(simulation::shared_secondary_bank_write, false); - // Synchronize the shared secondary bank amongst all MPI ranks, such // that each MPI rank has an approximately equal number of secondary // tracks. Also reports the total number of secondaries alive across @@ -983,16 +1069,17 @@ void transport_history_based_shared_secondary() simulation::shared_secondary_bank_read.size()); std::fill(simulation::progeny_per_particle.begin(), simulation::progeny_per_particle.end(), 0); + thread_banks.resize(num_threads()); // Transport all secondary tracks from the shared secondary bank #pragma omp parallel { - vector thread_bank; + auto& thread_bank = thread_banks[thread_num()]; + Particle p; #pragma omp for schedule(runtime) for (int64_t i = 1; i <= simulation::shared_secondary_bank_read.size(); i++) { - Particle p; initialize_particle_track(p, i, true); SourceSite& site = simulation::shared_secondary_bank_read[i - 1]; p.event_revive_from_secondary(site); @@ -1000,18 +1087,14 @@ void transport_history_based_shared_secondary() for (auto& secondary_site : p.local_secondary_bank()) { thread_bank.push_back(secondary_site); } - } - - // Drain thread-local bank into the shared secondary bank (once per - // thread) -#pragma omp critical(SharedSecondaryBank) - { - for (auto& secondary_site : thread_bank) { - simulation::shared_secondary_bank_write.thread_unsafe_append( - secondary_site); - } + p.local_secondary_bank().clear(); } } // End of transport loop over tracks in shared secondary bank + simulation::shared_secondary_bank_write = + std::move(simulation::shared_secondary_bank_read); + simulation::shared_secondary_bank_read = SharedArray(); + collect_sorted_history_secondary_banks(thread_banks); + thread_banks.clear(); n_generation_depth++; simulation::simulation_tracks_completed += alive_secondary; } // End of loop over secondary generations @@ -1075,13 +1158,7 @@ void transport_event_based_shared_secondary() process_transport_events(); process_death_events(n_particles); - // Collect secondaries from all particle buffers into shared bank - for (int64_t i = 0; i < n_particles; i++) { - for (auto& site : simulation::particles[i].local_secondary_bank()) { - simulation::shared_secondary_bank_write.thread_unsafe_append(site); - } - simulation::particles[i].local_secondary_bank().clear(); - } + collect_event_secondary_banks(n_particles); remaining_work -= n_particles; source_offset += n_particles; @@ -1145,13 +1222,7 @@ void transport_event_based_shared_secondary() process_transport_events(); process_death_events(n_particles); - // Collect secondaries from all particle buffers into shared bank - for (int64_t i = 0; i < n_particles; i++) { - for (auto& site : simulation::particles[i].local_secondary_bank()) { - simulation::shared_secondary_bank_write.thread_unsafe_append(site); - } - simulation::particles[i].local_secondary_bank().clear(); - } + collect_event_secondary_banks(n_particles); sec_remaining -= n_particles; sec_offset += n_particles; diff --git a/src/source.cpp b/src/source.cpp index 1e783b623a..12bbdf84e4 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -4,7 +4,9 @@ #define HAS_DYNAMIC_LINKING #endif -#include // for move +#include // for max +#include // for sin, cos, abs +#include // for move #ifdef HAS_DYNAMIC_LINKING #include // for dlopen, dlsym, dlclose, dlerror @@ -16,17 +18,20 @@ #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/cell.h" +#include "openmc/constants.h" #include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/geometry.h" #include "openmc/hdf5_interface.h" #include "openmc/material.h" +#include "openmc/math_functions.h" #include "openmc/mcpl_interface.h" #include "openmc/memory.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" +#include "openmc/random_dist.h" #include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" @@ -100,6 +105,8 @@ unique_ptr Source::create(pugi::xml_node node) return make_unique(node); } else if (source_type == "mesh") { return make_unique(node); + } else if (source_type == "tokamak") { + return make_unique(node); } else { fatal_error(fmt::format("Invalid source type '{}' found.", source_type)); } @@ -673,6 +680,471 @@ SourceSite MeshSource::sample(uint64_t* seed) const return source(element)->sample_with_constraints(seed); } +//============================================================================== +// TokamakSource implementation +//============================================================================== + +TokamakSource::TokamakSource(pugi::xml_node node) : Source(node) +{ + // Read geometry parameters + major_radius_ = std::stod(get_node_value(node, "major_radius")); + minor_radius_ = std::stod(get_node_value(node, "minor_radius")); + elongation_ = std::stod(get_node_value(node, "elongation")); + triangularity_ = std::stod(get_node_value(node, "triangularity")); + shafranov_shift_ = std::stod(get_node_value(node, "shafranov_shift")); + + // Read optional vertical shift + if (check_for_node(node, "vertical_shift")) { + vertical_shift_ = std::stod(get_node_value(node, "vertical_shift")); + } else { + vertical_shift_ = 0.0; + } + + // Read optional toroidal angle bounds + if (check_for_node(node, "phi_start")) { + phi_start_ = std::stod(get_node_value(node, "phi_start")); + } else { + phi_start_ = 0.0; + } + if (check_for_node(node, "phi_extent")) { + phi_extent_ = std::stod(get_node_value(node, "phi_extent")); + } else { + phi_extent_ = 2.0 * PI; + } + if (check_for_node(node, "n_alpha")) { + n_alpha_ = std::stoi(get_node_value(node, "n_alpha")); + } else { + n_alpha_ = 101; // Default + } + + // Read emission profile + r_over_a_ = get_node_array(node, "r_over_a"); + emission_density_ = get_node_array(node, "emission_density"); + + // Read energy distribution(s) + for (auto energy_node : node.children("energy")) { + energy_dists_.push_back(distribution_from_xml(energy_node)); + } + + // Read optional time distribution; default to a delta distribution at t=0 + // for the same behavior as IndependentSource + if (check_for_node(node, "time")) { + time_ = distribution_from_xml(node.child("time")); + } else { + double T[] {0.0}; + double p[] {1.0}; + time_ = UPtrDist {new Discrete {T, p, 1}}; + } + + // Validate inputs + if (emission_density_.size() != r_over_a_.size()) { + fatal_error("TokamakSource: emission_density and r_over_a must have the " + "same length."); + } + if (r_over_a_.size() < 2) { + fatal_error( + "TokamakSource: At least 2 radial points are required for profiles."); + } + if (r_over_a_.front() != 0.0) { + fatal_error("TokamakSource: r_over_a must start at 0."); + } + if (r_over_a_.back() != 1.0) { + fatal_error("TokamakSource: r_over_a must end at 1."); + } + for (size_t i = 1; i < r_over_a_.size(); ++i) { + if (r_over_a_[i] <= r_over_a_[i - 1]) { + fatal_error("TokamakSource: r_over_a must be strictly increasing."); + } + } + for (size_t i = 0; i < emission_density_.size(); ++i) { + if (emission_density_[i] < 0.0) { + fatal_error("TokamakSource: emission_density values cannot be negative."); + } + } + if (major_radius_ <= 0.0) { + fatal_error("TokamakSource: major_radius must be > 0."); + } + if (minor_radius_ <= 0.0) { + fatal_error("TokamakSource: minor_radius must be > 0."); + } + if (minor_radius_ >= major_radius_) { + fatal_error("TokamakSource: minor_radius must be less than major_radius."); + } + if (elongation_ <= 0.0) { + fatal_error("TokamakSource: elongation must be > 0."); + } + if (triangularity_ < -1.0 || triangularity_ > 1.0) { + fatal_error("TokamakSource: triangularity must be in the range [-1, 1]."); + } + if (shafranov_shift_ < 0.0) { + fatal_error("TokamakSource: shafranov_shift must be >= 0."); + } + if (shafranov_shift_ >= 0.5 * minor_radius_) { + fatal_error("TokamakSource: shafranov_shift must be less than half the " + "minor radius."); + } + if (phi_extent_ <= 0.0 || phi_extent_ > 2.0 * PI) { + fatal_error("TokamakSource: phi_extent must be > 0 and <= 2*pi."); + } + if (n_alpha_ <= 2) { + fatal_error("TokamakSource: n_alpha must be > 2."); + } + if (n_alpha_ < 51) { + warning("TokamakSource: n_alpha values below 51 may introduce noticeable " + "discretization bias in source sampling."); + } + if (energy_dists_.empty()) { + fatal_error("TokamakSource: At least one energy distribution is required."); + } + if (energy_dists_.size() != 1 && energy_dists_.size() != r_over_a_.size()) { + fatal_error("TokamakSource: energy distributions must be either 1 (for all " + "r) or match the number of r_over_a points."); + } + + // Compute normalized geometry parameters + epsilon_ = minor_radius_ / major_radius_; + delta_tilde_ = shafranov_shift_ / minor_radius_; + + // Initialize isotropic angular distribution + angle_ = UPtrAngle {new Isotropic()}; + + precompute_sampling_distributions(); +} + +void TokamakSource::precompute_sampling_distributions() +{ + // Use precomputed normalized geometry parameters + double eps = epsilon_; // Inverse aspect ratio (a/R0) + double Dt = delta_tilde_; // Normalized Shafranov shift (Delta/a) + double delta = triangularity_; + + //========================================================================== + // RADIAL CDF (computed first since it's simpler and sampled first) + //========================================================================== + // The marginal radial PDF is obtained by analytically integrating the joint + // distribution f(r_tilde, alpha) over alpha. The result is: + // + // p(r_tilde) ~ S(r_tilde) * [(1 + eps*Dt)*r_tilde + // - (3/8)*c1*eps*r_tilde^2 + // - 2*eps*Dt*r_tilde^3] + // + // where the Bessel function coefficients are: + // c0 = J_0(delta) + J_2(delta) + // c1 = (J_1(2*delta) + J_3(2*delta)) / c0 + // + // For delta -> 0, c0 -> 1 and c1 -> 0, giving the circular cross-section + // limit. + + // Compute Bessel function coefficients. openmc::cyl_bessel_j handles + // negative arguments (negative triangularity) via the parity relation + // J_n(-x) = (-1)^n * J_n(x). + double J0_d = cyl_bessel_j(0, delta); + double J2_d = cyl_bessel_j(2, delta); + double J1_2d = cyl_bessel_j(1, 2.0 * delta); + double J3_2d = cyl_bessel_j(3, 2.0 * delta); + double c0 = J0_d + J2_d; + double c1 = (J1_2d + J3_2d) / c0; + + // Coefficients for the radial polynomial: A*r - B*r^2 - C*r^3 + radial_poly_a_ = 1.0 + eps * Dt; + radial_poly_b_ = 0.375 * c1 * eps; // 3/8 * c1 * eps + radial_poly_c_ = 2.0 * eps * Dt; + + // Build a refined radial grid that retains the user-specified grid points. + // The emission density is interpreted as linear-linear between those points. + constexpr int MIN_SUBINTERVALS = 8; + constexpr double MAX_GRID_SPACING = 1.0e-3; + vector radial_grid {r_over_a_.front()}; + vector radial_emission {emission_density_.front()}; + for (size_t i = 1; i < r_over_a_.size(); ++i) { + double r_lo = r_over_a_[i - 1]; + double r_hi = r_over_a_[i]; + double s_lo = emission_density_[i - 1]; + double s_hi = emission_density_[i]; + int n_subintervals = std::max(MIN_SUBINTERVALS, + static_cast(std::ceil((r_hi - r_lo) / MAX_GRID_SPACING))); + for (int j = 1; j <= n_subintervals; ++j) { + double t = static_cast(j) / n_subintervals; + radial_grid.push_back(r_lo + t * (r_hi - r_lo)); + radial_emission.push_back(s_lo + t * (s_hi - s_lo)); + } + } + + vector radial_pdf(radial_grid.size()); + for (size_t i = 0; i < radial_grid.size(); ++i) { + double r = radial_grid[i]; + // p(r) ~ S(r) * [A*r - B*r^2 - C*r^3] + double geometric_factor = + radial_poly_a_ * r - radial_poly_b_ * r * r - radial_poly_c_ * r * r * r; + radial_pdf[i] = radial_emission[i] * std::max(0.0, geometric_factor); + } + + // Check that the refined profile contains positive probability mass before + // constructing the normalized tabular distribution. + double total = 0.0; + for (size_t i = 1; i < radial_grid.size(); ++i) { + total += 0.5 * (radial_pdf[i - 1] + radial_pdf[i]) * + (radial_grid[i] - radial_grid[i - 1]); + } + if (total <= 0.0) { + fatal_error( + "TokamakSource: Integrated emission density is zero or negative. " + "Check emission_density profile."); + } + radial_dist_ = make_unique(radial_grid.data(), radial_pdf.data(), + radial_grid.size(), Interpolation::lin_lin); + + //========================================================================== + // POLOIDAL CDFs (for conditional sampling of alpha given r) + //========================================================================== + // The conditional distribution P(alpha | r) is a mixture: + // P(alpha | r) ~ sum_k w_k(r) * I_hat_k * p_k(alpha) + // where: + // - w_k(r) are the "dynamic" Bernstein weight functions (depend on r) + // - I_hat_k are the "static" normalized integrals (precomputed constants) + // - p_k(alpha) are the normalized basis distributions (precomputed CDFs) + // + // The static weights I_hat_k = I_k / (2*pi*c0) are: + // I_hat_0 = 1 + eps*Dt + // I_hat_1 = 1 + eps*Dt - (3/16)*c1*eps + // I_hat_2 = 1 - (3/8)*c1*eps + // I_hat_3 = 1 + eps*Dt + // I_hat_4 = 1 + (1/2)*eps*Dt - (3/16)*c1*eps + // I_hat_5 = 1 - eps*Dt - (3/8)*c1*eps + + // Compute static weights analytically + poloidal_integrals_[0] = 1.0 + eps * Dt; + poloidal_integrals_[1] = 1.0 + eps * Dt - 0.1875 * c1 * eps; // 3/16 = 0.1875 + poloidal_integrals_[2] = 1.0 - 0.375 * c1 * eps; // 3/8 = 0.375 + poloidal_integrals_[3] = 1.0 + eps * Dt; + poloidal_integrals_[4] = 1.0 + 0.5 * eps * Dt - 0.1875 * c1 * eps; + poloidal_integrals_[5] = 1.0 - eps * Dt - 0.375 * c1 * eps; + + // Build the alpha grid on [0, pi] (half domain due to up-down symmetry) + int n_alpha = n_alpha_; + vector alpha_grid(n_alpha); + double dalpha = PI / (n_alpha - 1); + for (int i = 0; i < n_alpha; ++i) { + alpha_grid[i] = i * dalpha; + } + + // Compute basis function values g_k(alpha) for tabular distributions + // Using Bernstein form: + // R_tilde = b0*(1-r)^2 + 2*b1*r*(1-r) + b2*r^2 + // J_tilde = b3*(1-r) + b4*r + // with: + // b0(alpha) = 1 + eps*Dt + // b1(alpha) = b0 + (eps/2)*cos(psi), psi = alpha + delta*sin(alpha) + // b2(alpha) = 1 + eps*cos(psi) + // b3(alpha) = cos(delta*sin(alpha)) + // + (delta/4)*(cos(alpha - delta*sin(alpha)) + // - cos(3*alpha + delta*sin(alpha))) + // b4(alpha) = b3(alpha) - 2*Dt*cos(alpha) + + array, N_POLOIDAL_BASIS> basis; + for (int k = 0; k < N_POLOIDAL_BASIS; ++k) { + basis[k].resize(n_alpha); + } + + for (int i = 0; i < n_alpha; ++i) { + double alpha = alpha_grid[i]; + double sin_alpha = std::sin(alpha); + double cos_alpha = std::cos(alpha); + double delta_sin_alpha = delta * sin_alpha; + double psi = alpha + delta_sin_alpha; + double cos_psi = std::cos(psi); + + // Bernstein coefficients b0-b4 + double b0 = 1.0 + eps * Dt; + double b1 = b0 + 0.5 * eps * cos_psi; + double b2 = 1.0 + eps * cos_psi; + double b3 = + std::cos(delta_sin_alpha) + 0.25 * delta * + (std::cos(alpha - delta_sin_alpha) - + std::cos(3.0 * alpha + delta_sin_alpha)); + double b4 = b3 - 2.0 * Dt * cos_alpha; + + // 6 basis functions g_k(alpha) = b_i * b_j + basis[0][i] = b0 * b3; // w0 = (1-r)^3 + basis[1][i] = b1 * b3; // w1 = 2*r*(1-r)^2 + basis[2][i] = b2 * b3; // w2 = r^2*(1-r) + basis[3][i] = b0 * b4; // w3 = r*(1-r)^2 + basis[4][i] = b1 * b4; // w4 = 2*r^2*(1-r) + basis[5][i] = b2 * b4; // w5 = r^3 + } + + // Build a linear-linear distribution for each basis function p_k(alpha) + for (int k = 0; k < N_POLOIDAL_BASIS; ++k) { + poloidal_dists_[k] = make_unique( + alpha_grid.data(), basis[k].data(), n_alpha, Interpolation::lin_lin); + } +} + +double TokamakSource::sample_r_over_a(uint64_t* seed) const +{ + return radial_dist_->sample(seed).first; +} + +double TokamakSource::mixture_weight(int k, double r) const +{ + double s = 1.0 - r; + switch (k) { + case 0: + return s * s * s * poloidal_integrals_[0]; + case 1: + return 2.0 * r * s * s * poloidal_integrals_[1]; + case 2: + return r * r * s * poloidal_integrals_[2]; + case 3: + return r * s * s * poloidal_integrals_[3]; + case 4: + return 2.0 * r * r * s * poloidal_integrals_[4]; + case 5: + return r * r * r * poloidal_integrals_[5]; + default: + UNREACHABLE(); + } +} + +double TokamakSource::sample_poloidal_angle(double r_norm, uint64_t* seed) const +{ + // Sample from the conditional distribution P(alpha | r_tilde) using + // mixture sampling with 6 precomputed basis distributions. + // + // The conditional is: P(alpha | r) ~ sum_k w_k(r) * I_hat_k * p_k(alpha) + // where: + // - w_k(r) are the "dynamic" Bernstein weight functions + // - I_hat_k are the "static" normalized integrals (precomputed in + // poloidal_integrals_) + // - p_k(alpha) are the normalized, precomputed basis distributions + // + // The normalization sum_k w_k(r) * I_hat_k equals the radial geometric + // polynomial evaluated at r, which is known analytically. + // + // Algorithm: + // 1. Compute total from analytical normalization + // 2. Lazily evaluate mixture weights with early exit to select component k + // 3. Sample alpha from the selected basis distribution + + // Analytical normalization: sum_k w_k(r) * I_hat_k + double total = + radial_poly_a_ - radial_poly_b_ * r_norm - radial_poly_c_ * r_norm * r_norm; + double xi = prn(seed) * total; + + // Sample component via lazy evaluation with early exit + // Order optimized for peaked emission profiles: 0, 1, 4, 5, 3, 2 + constexpr int order[] = {0, 1, 4, 5, 3, 2}; + double cumsum = 0.0; + int component = order[N_POLOIDAL_BASIS - 1]; + for (int i = 0; i < N_POLOIDAL_BASIS; ++i) { + cumsum += mixture_weight(order[i], r_norm); + if (xi < cumsum) { + component = order[i]; + break; + } + } + + // Sample alpha from [0, pi] + double alpha = poloidal_dists_[component]->sample(seed).first; + + // Exploit up-down symmetry: randomly flip to [pi, 2*pi] with 50% probability + // This is equivalent to flipping the sign of Z in the final position + if (prn(seed) >= 0.5) { + alpha = 2.0 * PI - alpha; + } + return alpha; +} + +std::pair TokamakSource::sample_energy( + double r_norm, uint64_t* seed) const +{ + if (energy_dists_.size() == 1) { + // Single distribution for all r + return energy_dists_[0]->sample(seed); + } + + // Multiple distributions: stochastic selection between bracketing r points + // Find the interval containing r_norm + size_t i = lower_bound_index(r_over_a_.begin(), r_over_a_.end(), r_norm); + + // Handle boundary cases + if (i >= energy_dists_.size() - 1) { + return energy_dists_.back()->sample(seed); + } + + // Stochastic interpolation: randomly select one of the two bracketing + // distributions based on distance to each + double t = (r_norm - r_over_a_[i]) / (r_over_a_[i + 1] - r_over_a_[i]); + size_t idx = (prn(seed) < t) ? i + 1 : i; + return energy_dists_[idx]->sample(seed); +} + +Position TokamakSource::flux_to_cartesian( + double r, double alpha, double phi) const +{ + // Flux surface parameterization: + // R = R0 + r*cos(alpha + delta*sin(alpha)) + Delta*(1 - (r/a)^2) + // Z = kappa * r * sin(alpha) + // x = R * cos(phi) + // y = R * sin(phi) + // z = Z + + double psi = alpha + triangularity_ * std::sin(alpha); + double r_over_a_sq = (r * r) / (minor_radius_ * minor_radius_); + + double R = + major_radius_ + r * std::cos(psi) + shafranov_shift_ * (1.0 - r_over_a_sq); + double Z = elongation_ * r * std::sin(alpha); + + double x = R * std::cos(phi); + double y = R * std::sin(phi); + double z = Z; + + return {x, y, z}; +} + +SourceSite TokamakSource::sample(uint64_t* seed) const +{ + SourceSite site; + site.particle = ParticleType::neutron(); + site.wgt = 1.0; + site.delayed_group = 0; + + // 1. Sample r/a from radial CDF + double r_norm = sample_r_over_a(seed); + double r = r_norm * minor_radius_; + + // 2. Sample poloidal angle from conditional distribution P(alpha|r) + double alpha = sample_poloidal_angle(r_norm, seed); + + // 3. Sample toroidal angle uniformly in [phi_start, phi_start + phi_extent] + double phi = phi_start_ + phi_extent_ * prn(seed); + + // 4. Convert to Cartesian coordinates + site.r = flux_to_cartesian(r, alpha, phi); + + // 4a. Apply vertical shift if non-zero + if (vertical_shift_ != 0.0) { + site.r.z += vertical_shift_; + } + + // 5. Sample isotropic direction + site.u = angle_->sample(seed).first; + + // 6. Sample energy from distribution(s), applying the importance weight so + // that biased distributions are handled correctly + auto [E, E_wgt] = sample_energy(r_norm, seed); + site.E = E; + + // 7. Sample particle creation time + auto [time, time_wgt] = time_->sample(seed); + site.time = time; + + site.wgt *= E_wgt * time_wgt; + + return site; +} + //============================================================================== // Non-member functions //============================================================================== diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 7bce6ec137..241cbf0084 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2730,6 +2730,9 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) int orig_cell = p.coord(0).cell(); double orig_E_last = p.E_last(); + // Set particle in top level + p.n_coord() = 1; + for (auto i_tally : tallies) { auto& tally {*model::tallies[i_tally]}; @@ -2740,7 +2743,6 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) for (auto cell_id : cells) { // Temporarily change cell of particle - p.n_coord() = 1; p.coord(0).cell() = cell_id; // Determine index of cell in model::pulse_height_cells @@ -2756,20 +2758,20 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) // we skip the assume_separate break below. auto filter_iter = FilterBinIter(tally, p); auto end = FilterBinIter(tally, true, &p.filter_matches()); - if (filter_iter == end) - continue; + if (filter_iter != end) { - // Loop over filter bins. - for (; filter_iter != end; ++filter_iter) { - auto filter_index = filter_iter.index_; - auto filter_weight = filter_iter.weight_; + // Loop over filter bins. + for (; filter_iter != end; ++filter_iter) { + auto filter_index = filter_iter.index_; + auto filter_weight = filter_iter.weight_; - // Loop over scores. - for (auto score_index = 0; score_index < tally.scores_.size(); - ++score_index) { + // Loop over scores. + for (auto score_index = 0; score_index < tally.scores_.size(); + ++score_index) { #pragma omp atomic - tally.results_(filter_index, score_index, TallyResult::VALUE) += - filter_weight; + tally.results_(filter_index, score_index, TallyResult::VALUE) += + filter_weight; + } } } @@ -2777,10 +2779,10 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) for (auto& match : p.filter_matches()) match.bins_present_ = false; } - // Restore cell/energy - p.n_coord() = orig_n_coord; - p.coord(0).cell() = orig_cell; - p.E_last() = orig_E_last; } + // Restore cell/energy + p.n_coord() = orig_n_coord; + p.coord(0).cell() = orig_cell; + p.E_last() = orig_E_last; } } // namespace openmc diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 416d1421c3..704d4de66e 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -837,7 +837,8 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) ratio_)); if (ratio_ <= 1.0) fatal_error(fmt::format("Invalid weight window ratio '{}' (<= 1.0) " - "specified for weight window generation")); + "specified for weight window generation", + ratio_)); // create a matching weight windows object auto wws = WeightWindows::create(); diff --git a/tests/cpp_unit_tests/test_distribution.cpp b/tests/cpp_unit_tests/test_distribution.cpp index e46088a630..479e7c8a7a 100644 --- a/tests/cpp_unit_tests/test_distribution.cpp +++ b/tests/cpp_unit_tests/test_distribution.cpp @@ -82,6 +82,31 @@ TEST_CASE("Test alias sampling method for pugixml constructor") } } +TEST_CASE("Test sampling a large linear-linear tabular distribution") +{ + constexpr int n_points = 10001; + constexpr int n_samples = 200000; + openmc::vector x(n_points); + openmc::vector p(n_points); + for (int i = 0; i < n_points; ++i) { + x[i] = static_cast(i) / (n_points - 1); + p[i] = 2.0 * x[i]; + } + + openmc::Tabular dist( + x.data(), p.data(), n_points, openmc::Interpolation::lin_lin); + uint64_t seed = openmc::init_seed(0, 0); + + double mean = 0.0; + for (int i = 0; i < n_samples; ++i) { + mean += dist.sample(&seed).first; + } + mean /= n_samples; + + // The normalized PDF is 2x on [0, 1], which has a mean of 2/3. + REQUIRE_THAT(mean, Catch::Matchers::WithinAbs(2.0 / 3.0, 0.003)); +} + TEST_CASE("Test construction of SpatialBox with parameters") { openmc::Position ll {-1, -2, -3}; diff --git a/tests/cpp_unit_tests/test_math.cpp b/tests/cpp_unit_tests/test_math.cpp index 1ad7c4b709..7467aa1fe7 100644 --- a/tests/cpp_unit_tests/test_math.cpp +++ b/tests/cpp_unit_tests/test_math.cpp @@ -46,6 +46,21 @@ TEST_CASE("Test t_percentile") } } +TEST_CASE("Test cylindrical Bessel functions") +{ + constexpr double x = 2.0; + constexpr double expected[] {0.22389077914123567, 0.57672480775687339, + 0.35283402861563773, 0.12894324947440205}; + + for (int n = 0; n < 4; ++n) { + REQUIRE_THAT(openmc::cyl_bessel_j(n, x), + Catch::Matchers::WithinRel(expected[n], 1.0e-14)); + REQUIRE_THAT(openmc::cyl_bessel_j(n, -x), + Catch::Matchers::WithinRel( + (n % 2 == 0 ? 1.0 : -1.0) * expected[n], 1.0e-14)); + } +} + TEST_CASE("Test calc_pn") { // The reference solutions come from scipy.special.eval_legendre @@ -355,3 +370,24 @@ TEST_CASE("Test broaden_wmp_polynomials") REQUIRE_THAT(ref_val, Catch::Matchers::Approx(test_val)); } } + +TEST_CASE("Test isclose") +{ + using openmc::isclose; + + // Identical values are always close, regardless of tolerances. + REQUIRE(isclose(1.0, 1.0, 0.0, 0.0)); + REQUIRE(isclose(0.0, 0.0, 0.0, 0.0)); + + // Absolute tolerance governs comparisons near zero. + REQUIRE(isclose(0.0, 1e-15, 0.0, 1e-14)); + REQUIRE_FALSE(isclose(0.0, 1e-13, 0.0, 1e-14)); + + // Relative tolerance scales with the magnitude of the operands. + REQUIRE(isclose(1.0e12, 1.0e12 + 1.0, 1e-12, 0.0)); + REQUIRE_FALSE(isclose(1.0e12, 1.0e12 + 10.0, 1e-12, 0.0)); + + // The looser of the two tolerances wins. + REQUIRE(isclose(1.0, 1.0 + 1e-13, 0.0, 1e-12)); + REQUIRE(isclose(1.0e6, 1.0e6 + 1e-4, 1e-9, 0.0)); +} diff --git a/tests/regression_tests/lattice_large_pitch/__init__.py b/tests/regression_tests/lattice_large_pitch/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/lattice_large_pitch/inputs_true.dat b/tests/regression_tests/lattice_large_pitch/inputs_true.dat new file mode 100644 index 0000000000..4dd4f0bd7b --- /dev/null +++ b/tests/regression_tests/lattice_large_pitch/inputs_true.dat @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + 100.0 100.0 + 2 + 3 3 + -150.0 -150.0 + +1 2 1 +2 1 2 +1 2 1 + + + + + + + + fixed source + 1000 + 10 + + + + 6 6 + -150.0 -150.0 + 150.0 150.0 + + + 1 + + + 1 + flux + + + diff --git a/tests/regression_tests/lattice_large_pitch/results_true.dat b/tests/regression_tests/lattice_large_pitch/results_true.dat new file mode 100644 index 0000000000..58c4e306dc --- /dev/null +++ b/tests/regression_tests/lattice_large_pitch/results_true.dat @@ -0,0 +1,73 @@ +tally 1: +1.749265E+01 +3.275925E+01 +4.159430E+01 +1.799568E+02 +7.027454E+01 +4.967568E+02 +6.789716E+01 +4.652258E+02 +4.272935E+01 +1.861129E+02 +2.073342E+01 +4.374279E+01 +4.028297E+01 +1.641642E+02 +6.734263E+01 +4.627171E+02 +1.037657E+02 +1.078682E+03 +1.108394E+02 +1.240927E+03 +7.236770E+01 +5.302034E+02 +4.356361E+01 +1.921990E+02 +6.731536E+01 +4.615977E+02 +9.850684E+01 +9.810272E+02 +5.164863E+02 +2.670336E+04 +5.097770E+02 +2.604770E+04 +1.064847E+02 +1.137536E+03 +7.052986E+01 +5.003220E+02 +7.065046E+01 +5.036723E+02 +1.044890E+02 +1.103250E+03 +5.121544E+02 +2.632389E+04 +5.065331E+02 +2.570752E+04 +1.044794E+02 +1.101249E+03 +6.569142E+01 +4.361120E+02 +4.739812E+01 +2.313289E+02 +7.402284E+01 +5.591689E+02 +1.066905E+02 +1.149521E+03 +1.038665E+02 +1.086813E+03 +7.160008E+01 +5.217744E+02 +4.219705E+01 +1.816100E+02 +2.089838E+01 +4.464899E+01 +4.154271E+01 +1.745866E+02 +7.081253E+01 +5.049997E+02 +6.673273E+01 +4.509038E+02 +4.184624E+01 +1.771546E+02 +1.853898E+01 +3.538608E+01 diff --git a/tests/regression_tests/lattice_large_pitch/test.py b/tests/regression_tests/lattice_large_pitch/test.py new file mode 100644 index 0000000000..0f409badde --- /dev/null +++ b/tests/regression_tests/lattice_large_pitch/test.py @@ -0,0 +1,70 @@ +"""Regression test for rectangular lattices with large pitch values. + +Large pitches used to trigger a segmentation fault in ``RectLattice::distance`` +because the boundary-crossing check compared a reconstructed crossing position +against an absolute tolerance that did not scale with the geometry (see #3852). +This test transports particles across a lattice with a large pitch to ensure the +crossing logic remains robust. + +""" + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + model = openmc.Model() + + # Large pitch that previously caused floating-point cancellation + pitch = 100.0 + n = 3 + + air = openmc.Material() + air.set_density('g/cm3', 0.001) + air.add_nuclide('N14', 1.0) + metal = openmc.Material() + metal.set_density('g/cm3', 7.0) + metal.add_nuclide('Fe56', 1.0) + + metal_cell = openmc.Cell(fill=metal) + metal_uni = openmc.Universe(cells=[metal_cell]) + air_cell = openmc.Cell(fill=air) + air_uni = openmc.Universe(cells=[air_cell]) + + lattice = openmc.RectLattice() + lattice.lower_left = (-pitch*n/2, -pitch*n/2) + lattice.pitch = (pitch, pitch) + lattice.outer = air_uni + lattice.universes = [ + [metal_uni, air_uni, metal_uni], + [air_uni, metal_uni, air_uni], + [metal_uni, air_uni, metal_uni], + ] + + box = openmc.model.RectangularPrism(pitch*n, pitch*n, boundary_type='vacuum') + root_cell = openmc.Cell(region=-box, fill=lattice) + model.geometry = openmc.Geometry([root_cell]) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 10 + model.settings.particles = 1000 + + mesh = openmc.RegularMesh() + mesh.dimension = (6, 6) + mesh.lower_left = (-pitch*n/2, -pitch*n/2) + mesh.upper_right = (pitch*n/2, pitch*n/2) + mesh_filter = openmc.MeshFilter(mesh) + tally = openmc.Tally(tally_id=1) + tally.filters = [mesh_filter] + tally.scores = ['flux'] + model.tallies = [tally] + + return model + + +def test_lattice_large_pitch(model): + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/weightwindows/test.py b/tests/regression_tests/weightwindows/test.py index 1d3b063bd4..c5b6bd889d 100644 --- a/tests/regression_tests/weightwindows/test.py +++ b/tests/regression_tests/weightwindows/test.py @@ -119,6 +119,61 @@ def test_weightwindows(shared_secondary, subdir): test.main() +def test_zero_bound_windows_play_no_game(tmp_path): + # A weight window lower bound of zero means no weight window information + # exists there (MCNP wwinp files use zero to turn the game off in a cell), + # so transport must proceed as if weight windows were disabled. Previously, + # zero-bound windows demanded a split at every checkpoint (weight/0 -> + # max_split), multiplying the particle population until terminated by the + # split or weight cutoff limits. + model = build_model(False) + for ww in model.settings.weight_windows: + ww.lower_ww_bounds = np.zeros_like(ww.lower_ww_bounds) + ww.upper_ww_bounds = np.zeros_like(ww.upper_ww_bounds) + sp_zero = model.run(cwd=tmp_path / 'zero_windows') + + model.settings.weight_windows_on = False + sp_off = model.run(cwd=tmp_path / 'windows_off') + + with openmc.StatePoint(sp_zero) as sp: + flux_zero = list(sp.tallies.values())[0].mean + with openmc.StatePoint(sp_off) as sp: + flux_off = list(sp.tallies.values())[0].mean + + np.testing.assert_allclose(flux_zero, flux_off, rtol=1e-12) + + +def test_zero_and_negative_bounds_equivalent(tmp_path): + # Zero and negative lower bounds both mean that no weight window + # information exists in a cell (generators mark such cells with -1, and + # MCNP wwinp files use zero), so they must produce identical transport. + # Unlike the all-zero case above, here particles are born under valid + # windows and encounter the no-information region in flight; previously a + # zero lower bound in that situation demanded a split at every checkpoint + # in the cell (weight/0 -> max_split), multiplying the particle population, + # while -1 played no game. + def run_with(bound_value, subdir): + model = build_model(False) + for ww in model.settings.weight_windows: + lb = np.array(ww.lower_ww_bounds, copy=True) + ub = np.array(ww.upper_ww_bounds, copy=True) + lb[3:, :, :, :] = bound_value + ub[3:, :, :, :] = bound_value + ww.lower_ww_bounds = lb + ww.upper_ww_bounds = ub + return model.run(cwd=tmp_path / subdir) + + sp_zero = run_with(0.0, 'zero_region') + sp_negative = run_with(-1.0, 'negative_region') + + with openmc.StatePoint(sp_zero) as sp: + flux_zero = list(sp.tallies.values())[0].mean + with openmc.StatePoint(sp_negative) as sp: + flux_negative = list(sp.tallies.values())[0].mean + + np.testing.assert_allclose(flux_zero, flux_negative, rtol=1e-12) + + def test_wwinp_cylindrical(): ww = openmc.WeightWindowsList.from_wwinp('ww_n_cyl.txt')[0] diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py index a48c553255..f81f857e0c 100644 --- a/tests/unit_tests/test_data_decay.py +++ b/tests/unit_tests/test_data_decay.py @@ -155,16 +155,16 @@ def test_decay_photon_energy(): with pytest.raises(DataError): openmc.data.decay_photon_energy('I135') - # Set chain file to simple chain - openmc.config['chain_file'] = Path(__file__).parents[1] / "chain_simple.xml" + # Temporarily Set chain file to simple chain + with openmc.config.patch('chain_file', Path(__file__).parents[1] / 'chain_simple.xml'): - # Check strength of I135 source and presence of specific spectral line - src = openmc.data.decay_photon_energy('I135') - assert isinstance(src, openmc.stats.Discrete) - assert src.integral() == pytest.approx(3.920996223799345e-05) - assert 1260409. in src.x + # Check strength of I135 source and presence of specific spectral line + src = openmc.data.decay_photon_energy('I135') + assert isinstance(src, openmc.stats.Discrete) + assert src.integral() == pytest.approx(3.920996223799345e-05) + assert 1260409. in src.x - # Check Xe135 source, which should be tabular - src = openmc.data.decay_photon_energy('Xe135') - assert isinstance(src, openmc.stats.Tabular) - assert src.integral() == pytest.approx(2.076506258964966e-05) + # Check Xe135 source, which should be tabular + src = openmc.data.decay_photon_energy('Xe135') + assert isinstance(src, openmc.stats.Tabular) + assert src.integral() == pytest.approx(2.076506258964966e-05) diff --git a/tests/unit_tests/test_data_dose.py b/tests/unit_tests/test_data_dose.py index 4f18800143..32be5eb89c 100644 --- a/tests/unit_tests/test_data_dose.py +++ b/tests/unit_tests/test_data_dose.py @@ -34,6 +34,20 @@ def test_dose_coefficients(): assert energy[-1] == approx(20.0e6) assert dose[-1] == approx(338.0) + energy, dose = dose_coefficients( + 'neutron', data_source='icrp74', dose_quantity='ambient') + assert energy[0] == approx(1e-3) + assert dose[0] == approx(6.60) + assert energy[-1] == approx(20.0e6) + assert dose[-1] == approx(600) + + energy, dose = dose_coefficients( + 'photon', data_source='icrp74', dose_quantity='ambient') + assert energy[0] == approx(0.01e6) + assert dose[0] == approx(0.061) + assert energy[-1] == approx(10e6) + assert dose[-1] == approx(25.6) + # Invalid particle/geometry should raise an exception with raises(ValueError): dose_coefficients('slime', 'LAT') @@ -41,6 +55,14 @@ def test_dose_coefficients(): dose_coefficients('neutron', 'ZZ') with raises(ValueError): dose_coefficients('neutron', data_source='icrp7000') + with raises(ValueError): + dose_coefficients('neutron', dose_quantity='banana') + with raises(ValueError): + dose_coefficients( + 'neutron', data_source='icrp116', dose_quantity='ambient') + with raises(ValueError): + dose_coefficients( + 'neutron', 'ISO', data_source='icrp74', dose_quantity='ambient') with raises(ValueError) as excinfo: dose_coefficients("photons", data_source="icrp116") expected_particles = [ @@ -57,7 +79,8 @@ def test_dose_coefficients(): "proton", ] expected_msg = ( - "'photons' has no dose data in data source icrp116. " - f"Available particles for icrp116 are: {expected_particles}" + "'photons' has no effective dose data in data source icrp116. " + "Available particles for icrp116 with dose quantity effective are: " + f"{expected_particles}" ) assert str(excinfo.value) == expected_msg diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py index 14db689130..320d09eda5 100644 --- a/tests/unit_tests/test_data_misc.py +++ b/tests/unit_tests/test_data_misc.py @@ -7,6 +7,7 @@ from pathlib import Path import numpy as np import pytest import openmc.data +from openmc.deplete import Chain, Nuclide def test_data_library(tmpdir): @@ -134,7 +135,8 @@ def test_zam(): with pytest.raises(ValueError): openmc.data.zam('Am242-m1') -def test_half_life(): + +def test_half_life(tmp_path): assert openmc.data.half_life('H2') is None assert openmc.data.half_life('U235') == pytest.approx(2.22102e16) assert openmc.data.half_life('Am242') == pytest.approx(57672.0) @@ -143,3 +145,32 @@ def test_half_life(): assert openmc.data.decay_constant('U235') == pytest.approx(log(2.0)/2.22102e16) assert openmc.data.decay_constant('Am242') == pytest.approx(log(2.0)/57672.0) assert openmc.data.decay_constant('Am242_m1') == pytest.approx(log(2.0)/4449622000.0) + + # Create minimal chain with H3 and Am242 to test half-life and decay + # constant retrieval from chain file + chain = Chain() + h3 = Nuclide("H3") + h3.half_life = 1.0 + chain.add_nuclide(h3) + am242 = Nuclide("Am242") + chain.add_nuclide(am242) + + assert openmc.data.half_life('H3', chain_file=chain) == 1.0 + assert openmc.data.decay_constant('H3', chain_file=chain) == pytest.approx(log(2.0)) + + # Nuclides that are present but stable in the chain should not fall back to + # ENDF/B-VIII.0 data. + assert openmc.data.half_life('Am242', chain_file=chain) is None + assert openmc.data.decay_constant('Am242', chain_file=chain) == 0.0 + + # Nuclides missing from the chain fall back to ENDF/B-VIII.0 data. + assert openmc.data.half_life('U235', chain_file=chain) == pytest.approx(2.22102e16) + + chain_path = tmp_path / "chain.xml" + chain.export_to_xml(chain_path) + assert openmc.data.half_life('H3', chain_file=chain_path) == 1.0 + + endf_h3 = openmc.data.half_life('H3') + with openmc.config.patch('chain_file', chain_path): + assert openmc.data.half_life('H3', chain_file=None) == 1.0 + assert openmc.data.half_life('H3', chain_file=False) == endf_h3 diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 39c532c549..bc1a078b05 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -6,6 +6,7 @@ from pathlib import Path import numpy as np import pytest +import openmc import openmc.deplete @@ -38,6 +39,36 @@ def test_get_activity(res): np.testing.assert_allclose(a_xe135, a_xe135_ref) +def test_get_activity_chain_file(res, tmp_path): + """Tests evaluating activity with chain half-life data""" + _, a_endf = res.get_activity("1", by_nuclide=True, chain_file=False) + xe135_endf = np.array([a["Xe135"] for a in a_endf]) + + chain = openmc.deplete.Chain() + xe135 = openmc.deplete.Nuclide("Xe135") + xe135.half_life = openmc.data.half_life("Xe135") / 2.0 + chain.add_nuclide(xe135) + + t_chain, a_chain = res.get_activity("1", by_nuclide=True, chain_file=chain) + xe135_chain = np.array([a["Xe135"] for a in a_chain]) + + t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0]) + np.testing.assert_allclose(t_chain, t_ref) + np.testing.assert_allclose(xe135_chain, 2.0 * xe135_endf) + + chain_path = tmp_path / "chain.xml" + chain.export_to_xml(chain_path) + with openmc.config.patch('chain_file', chain_path): + _, a_config = res.get_activity("1", by_nuclide=True) + xe135_config = np.array([a["Xe135"] for a in a_config]) + np.testing.assert_allclose(xe135_config, xe135_chain) + + stable_chain = openmc.deplete.Chain() + stable_chain.add_nuclide(openmc.deplete.Nuclide("Xe135")) + _, a_stable = res.get_activity("1", by_nuclide=True, chain_file=stable_chain) + assert all(a["Xe135"] == 0.0 for a in a_stable) + + def test_get_atoms(res): """Tests evaluating single nuclide concentration.""" t, n = res.get_atoms("1", "Xe135") diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index f62306e834..4cfae0df28 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -926,6 +926,22 @@ def test_property_map(lib_init): assert np.allclose(expected_properties, properties, atol=1e-04) +def test_slice_data(lib_init): + expected_properties = np.array( + [[(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)], + [ (293.6, 6.55), (293.6, 10.29769), (293.6, 6.55)], + [(293.6, 0.740582), (293.6, 6.55), (293.6, 0.740582)]], dtype='float') + origin = (0.0, 0.0, 0.0) + _, properties = openmc.lib.slice_data( + origin, + width=(1.26, 1.26), + basis='xy', + pixels=(3, 3), + include_properties=True + ) + assert np.allclose(expected_properties, properties, atol=1e-04) + + def test_solid_raytrace_plot(lib_init, pincell_model): # Ensure plot mapping can be accessed and grows after allocation n0 = len(openmc.lib.plots) diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 89dfc03ddd..12d3633767 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -7,7 +7,7 @@ import numpy as np import openmc from openmc.data import decay_photon_energy -from openmc.deplete import Chain +from openmc.deplete import Chain, Nuclide import openmc.examples import openmc.model import openmc.stats @@ -614,6 +614,35 @@ def test_get_activity(): assert m4.get_activity(units='Ci/m3') == pytest.approx(ci/m3) +def test_get_activity_chain_file(tmp_path): + m = openmc.Material() + m.add_nuclide("H3", 1.0) + m.set_density('g/cm3', 1.0) + + chain = Chain() + h3 = Nuclide("H3") + h3.half_life = 1.0 + chain.add_nuclide(h3) + + atoms_per_bcm = m.get_nuclide_atom_densities()["H3"] + expected = np.log(2.0) * 1e24 * atoms_per_bcm + + assert m.get_activity(chain_file=chain) == pytest.approx(expected) + + chain_path = tmp_path / "chain.xml" + chain.export_to_xml(chain_path) + assert m.get_activity(chain_file=chain_path) == pytest.approx(expected) + + endf_activity = m.get_activity(chain_file=False) + with openmc.config.patch('chain_file', chain_path): + assert m.get_activity() == pytest.approx(expected) + assert m.get_activity(chain_file=False) == pytest.approx(endf_activity) + + stable_chain = Chain() + stable_chain.add_nuclide(Nuclide("H3")) + assert m.get_activity(chain_file=stable_chain) == 0.0 + + def test_get_decay_heat(): # Set chain file for testing openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml' diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 47a5214cc3..303adcd31e 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -660,6 +660,29 @@ def test_model_plot(): plt.close('all') +def test_model_plot_invalid_inputs(): + surface = openmc.Sphere(r=10.0, boundary_type="vacuum") + cell = openmc.Cell(region=-surface) + model = openmc.Model(openmc.Geometry([cell])) + + with pytest.raises(ValueError): + model.plot(n_samples=-1) + with pytest.raises(TypeError): + model.plot(n_samples=1.5) + with pytest.raises(ValueError): + model.plot(plane_tolerance=0.0) + with pytest.raises(TypeError): + model.plot(plane_tolerance='1') + with pytest.raises(ValueError): + model.plot(pixels=-1) + with pytest.raises(ValueError): + model.plot(pixels=(0, 100)) + with pytest.raises(ValueError): + model.plot(pixels=(100,)) + with pytest.raises(ValueError): + model.slice_data(u_span=(2, 0, 0), v_span=(0, 2, 0), pixels=-1) + + def test_model_id_map_initialization(run_in_tmpdir): model = openmc.examples.pwr_assembly() model.init_lib(output=False) diff --git a/tests/unit_tests/test_pulse_height.py b/tests/unit_tests/test_pulse_height.py new file mode 100644 index 0000000000..1f27cc6f26 --- /dev/null +++ b/tests/unit_tests/test_pulse_height.py @@ -0,0 +1,59 @@ +import numpy as np +import pytest +import openmc + + +@pytest.fixture +def model(): + openmc.reset_auto_ids() + model = openmc.Model() + + # Define materials + NaI = openmc.Material() + NaI.set_density('g/cm3', 3.7) + NaI.add_element('Na', 1.0) + NaI.add_element('I', 1.0) + + # Define geometry: two spheres in each other + s1 = openmc.Sphere(r=1) + s2 = openmc.Sphere(r=2, boundary_type='vacuum') + inner_sphere = openmc.Cell(name='inner sphere', fill=NaI, region=-s1) + outer_sphere = openmc.Cell(name='outer sphere', fill=NaI, region=+s1 & -s2) + model.geometry = openmc.Geometry([inner_sphere, outer_sphere]) + + # Define settings + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 10000 + model.settings.photon_transport = True + model.settings.source = openmc.IndependentSource( + energy=openmc.stats.delta_function(1e6), + particle='photon' + ) + + # Define tallies + energy_filter = openmc.EnergyFilter([1e3, 1e7]) + + tally1 = openmc.Tally() + tally1.scores = ['pulse-height'] + cell_filter1 = openmc.CellFilter([inner_sphere, outer_sphere]) + tally1.filters = [cell_filter1, energy_filter] + + tally2 = openmc.Tally() + tally2.scores = ['pulse-height'] + cell_filter2 = openmc.CellFilter([outer_sphere, inner_sphere]) + tally2.filters = [cell_filter2, energy_filter] + + model.tallies = [tally1, tally2] + return model + + +def test_pulse_height(model, run_in_tmpdir): + sp_path = model.run() + sp = openmc.StatePoint(sp_path) + t1 = sp.tallies[1].mean.squeeze() + t2 = sp.tallies[2].mean.squeeze() + + np.testing.assert_array_equal(t1, t2[::-1]) + + diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py index 266bd49763..5e617919d9 100644 --- a/tests/unit_tests/test_r2s.py +++ b/tests/unit_tests/test_r2s.py @@ -33,8 +33,8 @@ def simple_model_and_mesh(): # Simple settings with a point source settings = openmc.Settings() - settings.batches = 10 - settings.particles = 1000 + settings.batches = 2 + settings.particles = 250 settings.run_mode = 'fixed source' settings.source = openmc.IndependentSource() model = openmc.Model(geometry, settings=settings) @@ -46,6 +46,16 @@ def simple_model_and_mesh(): return model, (c1, c2), mesh +@pytest.fixture +def source_stage_manager(simple_model_and_mesh): + model, (c1, c2), _ = simple_model_and_mesh + r2s = R2SManager(model, [c1, c2]) + r2s.results['depletion_results'] = [None, None] + r2s.results['activation_materials'] = [c1.fill, c2.fill] + bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box} + return r2s, bounding_boxes + + def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): model, (c1, c2), mesh = simple_model_and_mesh @@ -59,9 +69,9 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): outdir = r2s.run( timesteps=[(1.0, 'd')], source_rates=[1.0], - photon_time_indices=[1], output_dir=tmp_path, chain_file=chain, + micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']}, ) # Check directories and files exist @@ -73,7 +83,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): assert (act / 'depletion_results.h5').exists() pt = Path(outdir) / 'photon_transport' assert (pt / 'tally_ids.json').exists() - assert (pt / 'time_1' / 'statepoint.10.h5').exists() + assert not (pt / 'time_0').exists() + assert (pt / 'time_1' / 'statepoint.2.h5').exists() # Basic results structure checks assert len(r2s.results['fluxes']) == 2 @@ -82,6 +93,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): assert len(r2s.results['mesh_material_volumes'][0]) == 2 assert len(r2s.results['activation_materials']) == 2 assert len(r2s.results['depletion_results']) == 2 + assert list(r2s.results['photon_sources']) == [1] + assert r2s.results['photon_sources'][1] # Check activation materials amats = r2s.results['activation_materials'] @@ -125,6 +138,7 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path): photon_time_indices=[1], output_dir=tmp_path, chain_file=chain, + micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']}, ) # Check that per-mesh MMV files were written @@ -137,7 +151,7 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path): assert (act / 'depletion_results.h5').exists() pt = Path(outdir) / 'photon_transport' assert (pt / 'tally_ids.json').exists() - assert (pt / 'time_1' / 'statepoint.10.h5').exists() + assert (pt / 'time_1' / 'statepoint.2.h5').exists() # Two meshes, each with 1 element containing both materials → # 2 element-material combinations per mesh, 4 total @@ -167,6 +181,9 @@ def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path): def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): model, (c1, c2), _ = simple_model_and_mesh + tally = openmc.Tally() + tally.scores = ['flux'] + model.tallies = [tally] # Use cell-based domains r2s = R2SManager(model, [c1, c2]) @@ -180,9 +197,11 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): timesteps=[(1.0, 'd')], source_rates=[1.0], photon_time_indices=[1], + by_parent_nuclide=True, output_dir=tmp_path, bounding_boxes=bounding_boxes, - chain_file=chain + chain_file=chain, + micro_kwargs={'nuclides': ['Ni58'], 'reactions': ['(n,p)']}, ) # Check directories and files exist @@ -193,13 +212,15 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): assert (act / 'depletion_results.h5').exists() pt = Path(outdir) / 'photon_transport' assert (pt / 'tally_ids.json').exists() - assert (pt / 'time_1' / 'statepoint.10.h5').exists() + assert (pt / 'time_1' / 'statepoint.2.h5').exists() # Basic results structure checks assert len(r2s.results['fluxes']) == 2 assert len(r2s.results['micros']) == 2 assert len(r2s.results['activation_materials']) == 2 assert len(r2s.results['depletion_results']) == 2 + assert r2s.photon_model.tallies[0].contains_filter( + openmc.ParentNuclideFilter) # Check activation materials amats = r2s.results['activation_materials'] @@ -217,3 +238,88 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): assert len(r2s_loaded.results['micros']) == 2 assert len(r2s_loaded.results['activation_materials']) == 2 assert len(r2s_loaded.results['depletion_results']) == 2 + + +def test_step4_requires_photon_sources(simple_model_and_mesh, tmp_path): + model, (c1, c2), _ = simple_model_and_mesh + r2s = R2SManager(model, [c1, c2]) + output_dir = tmp_path / 'photon' + + with pytest.raises(RuntimeError, match='step3_photon_source'): + r2s.step4_photon_transport(output_dir) + + r2s.results['photon_sources'] = {} + with pytest.raises(RuntimeError, match='No decay photon sources'): + r2s.step4_photon_transport(output_dir) + + assert not output_dir.exists() + + +def test_default_photon_times_skip_empty_sources( + source_stage_manager, tmp_path, monkeypatch +): + r2s, bounding_boxes = source_stage_manager + source = object() + sources_by_time = {0: [], 1: [source]} + monkeypatch.setattr( + r2s, '_create_photon_sources', + lambda time_index, work_items: sources_by_time[time_index]) + + r2s.step3_photon_source( + bounding_boxes=bounding_boxes, output_dir=tmp_path) + + assert r2s.results['photon_sources'] == {1: [source]} + + +def test_explicit_empty_photon_source_fails( + source_stage_manager, tmp_path, monkeypatch +): + r2s, bounding_boxes = source_stage_manager + source = object() + sources_by_time = {0: [], 1: [source]} + monkeypatch.setattr( + r2s, '_create_photon_sources', + lambda time_index, work_items: sources_by_time[time_index]) + r2s.results['photon_sources'] = {99: [source]} + + with pytest.raises(RuntimeError, match='requested time indices: 0'): + r2s.step3_photon_source( + [0, 1], bounding_boxes, output_dir=tmp_path) + + assert 'photon_sources' not in r2s.results + + +def test_default_photon_times_require_a_source( + source_stage_manager, tmp_path, monkeypatch +): + r2s, bounding_boxes = source_stage_manager + monkeypatch.setattr( + r2s, '_create_photon_sources', + lambda time_index, work_items: []) + + with pytest.raises(RuntimeError, match='at any depletion time'): + r2s.step3_photon_source( + bounding_boxes=bounding_boxes, output_dir=tmp_path) + + assert 'photon_sources' not in r2s.results + + +@pytest.mark.parametrize( + ('time_indices', 'exception'), + [ + ([], ValueError), + ([2], IndexError), + ([-3], IndexError), + ([1.0], TypeError), + ], +) +def test_photon_time_index_validation( + source_stage_manager, tmp_path, time_indices, exception +): + r2s, bounding_boxes = source_stage_manager + + with pytest.raises(exception): + r2s.step3_photon_source( + time_indices, bounding_boxes, output_dir=tmp_path) + + assert 'photon_sources' not in r2s.results diff --git a/tests/unit_tests/test_slice_data_overlap_info.py b/tests/unit_tests/test_slice_data_overlap_info.py new file mode 100644 index 0000000000..e5f08245d7 --- /dev/null +++ b/tests/unit_tests/test_slice_data_overlap_info.py @@ -0,0 +1,89 @@ +import pytest +import numpy as np +import openmc +import openmc.lib + +# Sentinel value matching _OVERLAP in plotmodel.py and OVERLAP in plot.cpp +_OVERLAP = -3 + + +@pytest.fixture(scope='module') +def overlap_model(): + openmc.reset_auto_ids() + + # Three cylinders: cyl1 and cyl2 overlap near x=0, cyl2 and cyl3 overlap + # near x=4. This gives us two spatially distinct overlap regions in one model. + mat1 = openmc.Material(components={'H1': 1.0}) + mat2 = openmc.Material(components={'H1': 1.0}) + mat3 = openmc.Material(components={'H1': 1.0}) + + # cyl1 and cyl2 overlap on the left, cyl2 and cyl3 overlap on the right + cyl1 = openmc.ZCylinder(x0=-2.0, r=2.5) + cyl2 = openmc.ZCylinder(x0=0.0, r=2.5) + cyl3 = openmc.ZCylinder(x0=2.0, r=2.5) + boundary = openmc.Sphere(r=20.0, boundary_type='vacuum') + cell1 = openmc.Cell(region=-cyl1, fill=mat1) + cell2 = openmc.Cell(region=-cyl2, fill=mat2) + cell3 = openmc.Cell(region=-cyl3, fill=mat3) + cell_outside = openmc.Cell(region=+cyl1 & +cyl2 & +cyl3 & -boundary) + geometry = openmc.Geometry([cell1, cell2, cell3, cell_outside]) + + settings = openmc.Settings() + settings.run_mode = 'fixed source' + settings.particles = 100 + settings.batches = 1 + model = openmc.Model(geometry=geometry, settings=settings) + + with openmc.lib.TemporarySession(model, args=['-s', '1']): + yield + + +def run_slice(origin=(0.0, 0.0, 0.0), width=(10.0, 6.0), show_overlaps=True): + # Helper that runs a slice over a region covering both overlap zones + geom_data, _ = openmc.lib.slice_data( + origin=origin, + width=width, + basis='xy', + pixels=(100, 60), + show_overlaps=show_overlaps, + include_properties=False, + ) + return geom_data + + +def test_overlaps_enabled(overlap_model): + # Run a single slice with overlap detection enabled and check all + # expected properties in one pass. + geom_data = run_slice() + overlap_info = openmc.lib.slice_data_overlap_info() + n = overlap_info.shape[0] + cell_ids = geom_data[:, :, 0] + + # cell_ids should contain values more negative than _OVERLAP; RasterData + # encodes each unique overlap as OVERLAP - overlap_idx - 1 into slot 2. + assert np.any(cell_ids < _OVERLAP) + + # overlap_keys should have 2 entries for the two distinct overlapping + # cylinder pairs in this model. + assert n == 2, f"Expected exactly 2 overlap entries, got {n}" + + # Each entry is a (universe_id, cell1_id, cell2_id) triple; verify values. + for i in range(n): + universe_id = int(overlap_info[i, 0]) + cell1_id = int(overlap_info[i, 1]) + cell2_id = int(overlap_info[i, 2]) + assert universe_id == 1 + assert cell1_id in {1, 2, 3} + assert cell2_id in {1, 2, 3} + assert cell1_id != cell2_id + + +def test_overlaps_disabled(overlap_model): + # With show_overlaps=False, set_overlap is never called and overlap_keys + # is never written to, so the image and map should both be clean. + geom_data = run_slice(show_overlaps=False) + overlap_info = openmc.lib.slice_data_overlap_info() + n = overlap_info.shape[0] + + assert not np.any(geom_data[:, :, 2] < _OVERLAP) + assert n == 0 diff --git a/tests/unit_tests/test_source_tokamak.py b/tests/unit_tests/test_source_tokamak.py new file mode 100644 index 0000000000..f926c2007b --- /dev/null +++ b/tests/unit_tests/test_source_tokamak.py @@ -0,0 +1,263 @@ +import numpy as np +import pytest + +import openmc +import openmc.stats + +from tests.unit_tests import assert_sample_mean + + +def make_source(**kwargs): + """Build a valid TokamakSource, overriding defaults via kwargs.""" + r_over_a = np.linspace(0.0, 1.0, 10) + params = dict( + major_radius=620.0, + minor_radius=200.0, + elongation=1.8, + triangularity=0.45, + shafranov_shift=10.0, + r_over_a=r_over_a, + emission_density=(1.0 - r_over_a**2), + energy=openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=2.0e4), + ) + params.update(kwargs) + return openmc.TokamakSource(**params) + + +def test_tokamak_source_roundtrip(): + src = make_source( + phi_start=0.1, phi_extent=np.pi, n_alpha=51, vertical_shift=5.0, + strength=2.0, time=openmc.stats.Uniform(0.0, 1e-6)) + + elem = src.to_xml_element() + assert elem.get('type') == 'tokamak' + + new = openmc.SourceBase.from_xml_element(elem) + assert isinstance(new, openmc.TokamakSource) + assert new.major_radius == src.major_radius + assert new.minor_radius == src.minor_radius + assert new.elongation == src.elongation + assert new.triangularity == src.triangularity + assert new.shafranov_shift == src.shafranov_shift + assert new.phi_start == src.phi_start + assert new.phi_extent == src.phi_extent + assert new.n_alpha == src.n_alpha + assert new.vertical_shift == src.vertical_shift + assert new.strength == src.strength + np.testing.assert_allclose(new.r_over_a, src.r_over_a) + np.testing.assert_allclose(new.emission_density, src.emission_density) + assert len(new.energy) == 1 + assert isinstance(new.time, openmc.stats.Uniform) + assert new.time.a == src.time.a + assert new.time.b == src.time.b + + +def test_tokamak_source_default_time(): + src = make_source() + assert src.time is None + + new = openmc.SourceBase.from_xml_element(src.to_xml_element()) + assert new.time is None + + with pytest.raises(TypeError): + make_source(time=1.0) + + +def test_tokamak_source_multiple_energies(): + r_over_a = np.linspace(0.0, 1.0, 5) + energies = [openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=kt) + for kt in (1.0e4, 1.5e4, 2.0e4, 2.5e4, 3.0e4)] + src = make_source(r_over_a=r_over_a, + emission_density=np.ones_like(r_over_a), + energy=energies) + assert len(src.energy) == len(r_over_a) + + new = openmc.SourceBase.from_xml_element(src.to_xml_element()) + assert len(new.energy) == len(r_over_a) + + +@pytest.mark.parametrize("kwargs, match", [ + (dict(minor_radius=700.0), "smaller than major_radius"), + (dict(shafranov_shift=150.0), "half the minor_radius"), + (dict(emission_density=np.ones(5)), "same length as r_over_a"), + (dict(energy=[openmc.stats.muir(14.08e6, 5.0, 2.0e4)] * 2), + "Number of energy distributions"), + (dict(r_over_a=np.linspace(0.1, 1.0, 10)), "must start at 0"), + (dict(r_over_a=np.linspace(0.0, 0.9, 10)), "must end at 1"), + (dict(emission_density=-np.linspace(0.0, 1.0, 10)), "cannot be negative"), + (dict(emission_density=np.zeros(10)), "must contain a positive value"), +]) +def test_tokamak_source_invalid(kwargs, match): + with pytest.raises(ValueError, match=match): + make_source(**kwargs) + + +@pytest.mark.parametrize("value", [-1.5, 1.5]) +def test_tokamak_source_invalid_triangularity(value): + with pytest.raises(ValueError): + make_source(triangularity=value) + + +def test_tokamak_source_invalid_n_alpha(): + with pytest.raises(ValueError): + make_source(n_alpha=2) + + +@pytest.mark.parametrize(("attribute", "value", "match"), [ + ("minor_radius", 700.0, "smaller than major_radius"), + ("shafranov_shift", 150.0, "half the minor_radius"), + ("emission_density", np.ones(5), "same length as r_over_a"), + ("energy", [openmc.stats.delta_function(1.0)] * 2, + "Number of energy distributions"), +]) +def test_tokamak_source_mutation_validation(attribute, value, match): + src = make_source() + setattr(src, attribute, value) + with pytest.raises(ValueError, match=match): + src.to_xml_element() + + +@pytest.mark.parametrize(("emission_density", "expected_mean"), [ + ([1.0, 1.0], 2.0 / 3.0), + ([1.0, 0.0], 0.5), +]) +def test_tokamak_source_radial_sampling( + run_in_tmpdir, emission_density, expected_mean +): + """Check radial sampling for profiles on the coarsest valid grid.""" + major_radius = 620.0 + minor_radius = 200.0 + src = make_source( + major_radius=major_radius, + minor_radius=minor_radius, + elongation=1.0, + triangularity=0.0, + shafranov_shift=0.0, + r_over_a=[0.0, 1.0], + emission_density=emission_density, + energy=openmc.stats.delta_function(14.07e6), + ) + + sphere = openmc.Sphere(r=2000.0, boundary_type='vacuum') + model = openmc.Model( + geometry=openmc.Geometry([openmc.Cell(region=-sphere)]), + settings=openmc.Settings( + particles=100, batches=1, run_mode='fixed source', source=src), + ) + + sites = model.sample_external_source(20_000) + xyz = np.array([site.r for site in sites]) + major_r = np.hypot(xyz[:, 0], xyz[:, 1]) + r_over_a = np.hypot(major_r - major_radius, xyz[:, 2]) / minor_radius + assert_sample_mean(r_over_a, expected_mean) + + +def test_tokamak_source_poloidal_sampling(run_in_tmpdir): + """Check linear-linear poloidal sampling on a coarse internal grid.""" + major_radius = 620.0 + minor_radius = 200.0 + with pytest.warns(UserWarning, match="below 51"): + src = make_source( + major_radius=major_radius, + minor_radius=minor_radius, + elongation=1.0, + triangularity=0.0, + shafranov_shift=0.0, + emission_density=np.ones(10), + n_alpha=3, + energy=openmc.stats.delta_function(14.07e6), + ) + + sphere = openmc.Sphere(r=2000.0, boundary_type='vacuum') + model = openmc.Model( + geometry=openmc.Geometry([openmc.Cell(region=-sphere)]), + settings=openmc.Settings( + particles=100, batches=1, run_mode='fixed source', source=src), + ) + + sites = model.sample_external_source(100_000) + xyz = np.array([site.r for site in sites]) + major_r = np.hypot(xyz[:, 0], xyz[:, 1]) + + # With three alpha points, linear interpolation of cos(alpha) on [0, pi] + # gives 1 - 2*alpha/pi. Integrating the resulting density gives this mean. + expected_R = ( + major_radius + + 2.0 * minor_radius**2 / (np.pi**2 * major_radius) + ) + assert_sample_mean(major_r, expected_R) + + +@pytest.mark.flaky(reruns=1) +def test_tokamak_source_sampling(run_in_tmpdir): + """Exercise the compiled C++ sampling path and check invariants. + + Sampled moments are compared against direct numerical quadrature of the + exact source density S(r)*R*|J|, where the Jacobian J of the flux-surface + map is computed from analytic partial derivatives. This check is + independent of the Bernstein-mixture factorization used by the + implementation. + """ + R0, a, kappa, delta = 620.0, 200.0, 1.8, -0.5 + shafranov, zshift = 40.0, 25.0 + phi_start, phi_extent = 0.5, np.pi / 2 + + # Fine grids to make discretization error negligible relative to + # statistical uncertainty + r_over_a = np.linspace(0.0, 1.0, 200) + src = make_source( + major_radius=R0, minor_radius=a, elongation=kappa, + triangularity=delta, shafranov_shift=shafranov, vertical_shift=zshift, + r_over_a=r_over_a, emission_density=1.0 - r_over_a**2, + phi_start=phi_start, phi_extent=phi_extent, n_alpha=201, + energy=openmc.stats.delta_function(14.07e6)) + + sphere = openmc.Sphere(r=2000.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere) + settings = openmc.Settings( + particles=100, batches=1, run_mode='fixed source', source=src) + model = openmc.Model(geometry=openmc.Geometry([cell]), settings=settings) + + n_samples = 20_000 + sites = model.sample_external_source(n_samples) + + xyz = np.array([s.r for s in sites]) + R = np.hypot(xyz[:, 0], xyz[:, 1]) + z = xyz[:, 2] + + # Energy, weight, and time invariants + assert np.all([s.E == 14.07e6 for s in sites]) + assert np.all([s.wgt == 1.0 for s in sites]) + assert np.all([s.time == 0.0 for s in sites]) + + # Toroidal angle within the requested sector + phi = np.arctan2(xyz[:, 1], xyz[:, 0]) + assert phi.min() >= phi_start + assert phi.max() <= phi_start + phi_extent + + # Positions bounded by the last closed flux surface + assert R.min() >= R0 - a + assert R.max() <= R0 + a + shafranov + assert np.abs(z - zshift).max() <= kappa * a + + # Up-down symmetry about the vertical shift + assert_sample_mean(z, zshift) + + # Reference moments by 2D quadrature of the exact density + r = np.linspace(0.0, 1.0, 1001)[:, np.newaxis] # r/a + alpha = np.linspace(0.0, 2 * np.pi, 2001)[np.newaxis, :] + psi = alpha + delta * np.sin(alpha) + R_map = R0 + a * r * np.cos(psi) + shafranov * (1.0 - r**2) + Z_map = kappa * a * r * np.sin(alpha) + dR_dr = a * np.cos(psi) - 2.0 * shafranov * r + dR_da = -a * r * np.sin(psi) * (1.0 + delta * np.cos(alpha)) + dZ_dr = kappa * a * np.sin(alpha) * np.ones_like(psi) + dZ_da = kappa * a * r * np.cos(alpha) + jac = np.abs(dR_dr * dZ_da - dR_da * dZ_dr) + dens = (1.0 - r**2) * R_map * jac + norm = dens.sum() + expected_R = (R_map * dens).sum() / norm + expected_z2 = (Z_map**2 * dens).sum() / norm + + assert_sample_mean(R, expected_R) + assert_sample_mean((z - zshift)**2, expected_z2) diff --git a/tests/unit_tests/test_surface_flux.py b/tests/unit_tests/test_surface_flux.py index 4e067ffe23..e4419252ed 100644 --- a/tests/unit_tests/test_surface_flux.py +++ b/tests/unit_tests/test_surface_flux.py @@ -98,6 +98,50 @@ def test_surface_filter_flux_angled(two_cell_model, run_in_tmpdir): assert flux_mean == pytest.approx(1.0 / mu) +def test_surface_tally_during_lattice_crossing(run_in_tmpdir): + openmc.reset_auto_ids() + model = openmc.Model() + + xmin = openmc.XPlane(-1.0, boundary_type="vacuum") + xmax = openmc.XPlane(1.0, boundary_type="vacuum") + ymin = openmc.YPlane(-1.0, boundary_type="vacuum") + ymax = openmc.YPlane(1.0, boundary_type="vacuum") + zmin = openmc.ZPlane(-1.0, boundary_type="vacuum") + zmax = openmc.ZPlane(1.0, boundary_type="vacuum") + + inner_cell = openmc.Cell() + inner_univ = openmc.Universe(cells=[inner_cell]) + + tile_cell = openmc.Cell(fill=inner_univ) + tile_univ = openmc.Universe(cells=[tile_cell]) + + lattice = openmc.RectLattice() + lattice.lower_left = (-1.0, -1.0) + lattice.pitch = (1.0, 2.0) + lattice.universes = [[tile_univ, tile_univ]] + + root_cell = openmc.Cell( + fill=lattice, region=+xmin & -xmax & +ymin & -ymax & +zmin & -zmax) + model.geometry = openmc.Geometry([root_cell]) + + src = openmc.IndependentSource() + src.space = openmc.stats.Point((-0.5, 0.0, 0.0)) + src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 5 + model.settings.source = src + + current_tally = openmc.Tally() + current_tally.filters = [openmc.SurfaceFilter(xmax)] + current_tally.scores = ['current'] + model.tallies = [current_tally] + + model.run(apply_tally_results=True) + assert current_tally.mean.flat[0] == pytest.approx(1.0) + + def test_cellfrom_filter_flux_directional(two_cell_model, run_in_tmpdir): """SurfaceFilter + CellFromFilter + flux scores only the correct direction.""" model, xmid, cell1, cell2 = two_cell_model