diff --git a/.gitignore b/.gitignore index dd8dfb14a..32ef7919a 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 61301bdb4..7594ed5f2 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/CMakeLists.txt b/CMakeLists.txt index d5af01610..61d2cc6c5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -532,7 +532,7 @@ else() endif() if(OPENMC_USE_DAGMC) - target_compile_definitions(libopenmc PRIVATE OPENMC_DAGMC_ENABLED) + target_compile_definitions(libopenmc PUBLIC OPENMC_DAGMC_ENABLED) target_link_libraries(libopenmc dagmc-shared) if(OPENMC_USE_UWUW) diff --git a/docs/Makefile b/docs/Makefile index a93338df3..e320432a2 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 000000000..75702eff7 --- /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 1777e795a..e924aa6e2 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 3a2e7fb93..61740239b 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 cda031389..5199c9438 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 f723db06e..cda0307ca 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 60ff0afae..5cd18bef1 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 fb0215916..7091757f2 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 b3911c8b3..f06a1eba4 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 21981fdaa..19f1e5d5f 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/random_ray.rst b/docs/source/usersguide/random_ray.rst index 57d3abfbd..ff319d669 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -1116,10 +1116,15 @@ external source is present in the problem. Simulation settings (e.g., number of rays, batches, etc.) will be identical for both calculations. At the conclusion of the run, all results (e.g., tallies, plots, etc.) will be derived from the adjoint flux rather than the forward flux but are not labeled -any differently. The initial forward flux solution will not be stored or -available in the final statepoint file. Those wishing to do analysis requiring -both the forward and adjoint solutions will need to run two separate -simulations and load both statepoint files. +any differently. When an initial forward solve is performed (i.e., when no +user-specified adjoint source is present), its output files are also written to +disk with a ``forward`` infix, so they are not overwritten by the subsequent +adjoint solve. This applies to the statepoint, ``tallies.out``, and any voxel +plots, e.g., ``statepoint.forward.N.h5`` and ``tallies.forward.out``; the +adjoint solve keeps the usual file names. This allows analyses requiring both +the forward and adjoint solutions to be performed from a single run. When +generating FW-CADIS weight windows, no weight window file is written for the +forward solve, as only the final adjoint-derived weight windows are meaningful. .. note:: Use of the automated diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index e8514b856..6faeb59e1 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 6b78145a4..9f6987d74 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/constants.h b/include/openmc/constants.h index 7baed25c5..26b4a224c 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -368,6 +368,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY }; enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID }; enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY }; enum class RandomRaySampleMethod { PRNG, HALTON, S2 }; +enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT }; //============================================================================== // Geometry Constants diff --git a/include/openmc/error.h b/include/openmc/error.h index d73795aee..1f6e15f49 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 107cc7d1f..43cf9bb58 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 dcc0e21fe..d3bccca7b 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 f72948f6e..44e82fd23 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 f97d31384..ba8ac84a1 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/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 6f51af34d..09414fd44 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -76,7 +76,10 @@ public: //---------------------------------------------------------------------------- // Static Data members static bool volume_normalized_flux_tallies_; - static bool adjoint_; // If the user wants outputs based on the adjoint flux + // If the user wants outputs based on the adjoint flux + static bool adjoint_requested_; + // The solve currently being executed + static RandomRaySolve solve_; static bool fw_cadis_local_; static double diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index ccd2cbe47..e186c549f 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -22,7 +22,7 @@ public: void apply_fixed_sources_and_mesh_domains(); void prepare_fw_fixed_sources_adjoint(); void prepare_local_fixed_sources_adjoint(); - void prepare_adjoint_simulation(bool fw_adjoint); + void prepare_adjoint_simulation(bool from_forward); void simulate(); void output_simulation_results() const; void instability_check( diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 65f2e6bc4..1d2bbe1e8 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -51,10 +51,10 @@ inline void hash_combine(size_t& seed, const size_t v) // every iteration. struct TallyTask { int tally_idx; - int filter_idx; + int64_t filter_idx; int score_idx; int score_type; - TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type) + TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type) : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx), score_type(score_type) {} @@ -690,7 +690,7 @@ private: // Private Methods // Helper function for indexing - inline int index(int64_t sr, int g) const { return sr * negroups_ + g; } + inline int64_t index(int64_t sr, int g) const { return sr * negroups_ + g; } }; } // namespace openmc diff --git a/include/openmc/shared_array.h b/include/openmc/shared_array.h index b309ca3f1..2829a2eb9 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 e307b1ed2..51b54a1d1 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/tallies/tally.h b/include/openmc/tallies/tally.h index 66e6ff764..ae604b732 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -97,9 +97,9 @@ public: //! Given already-set filters, set the stride lengths void set_strides(); - int32_t strides(int i) const { return strides_[i]; } + int64_t strides(int i) const { return strides_[i]; } - int32_t n_filter_bins() const { return n_filter_bins_; } + int64_t n_filter_bins() const { return n_filter_bins_; } bool multiply_density() const { return multiply_density_; } @@ -184,9 +184,9 @@ private: vector filters_; //!< Filter indices in global filters array //! Index strides assigned to each filter to support 1D indexing. - vector strides_; + vector strides_; - int32_t n_filter_bins_ {0}; + int64_t n_filter_bins_ {0}; //! Whether to multiply by atom density for reaction rates bool multiply_density_ {true}; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index d1aed2831..4303ddb7a 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -41,7 +41,7 @@ public: FilterBinIter& operator++(); - int index_ {1}; + int64_t index_ {1}; double weight_ {1.}; vector& filter_matches_; diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index a5d404133..d0b385d16 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 c22e54e7d..6e4840990 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 d49043b0a..87e0cf9ae 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 000000000..fe0036cbc --- /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 000000000..66031b2f0 --- /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/data/njoy.py b/openmc/data/njoy.py index ddc68cc37..7c99637f0 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -259,14 +259,14 @@ broadr / %%%%%%%%%%%%%%%%%%%%%%% Doppler broaden XS %%%%%%%%%%%%%%%%%%%%%%%%%%%% _TEMPLATE_HEATR = """ heatr / %%%%%%%%%%%%%%%%%%%%%%%%% Add heating kerma %%%%%%%%%%%%%%%%%%%%%%%%%%%% {nendf} {nheatr_in} {nheatr} / -{mat} 4 0 0 0 / +{mat} 4 0 0 0 0 {ed}/ 302 318 402 444 / """ _TEMPLATE_HEATR_LOCAL = """ heatr / %%%%%%%%%%%%%%%%% Add heating kerma (local photons) %%%%%%%%%%%%%%%%%%%% {nendf} {nheatr_in} {nheatr_local} / -{mat} 4 0 0 1 / +{mat} 4 0 0 1 0 {ed}/ 302 318 402 444 / """ @@ -411,7 +411,7 @@ def make_pendf(filename, pendf='pendf', **kwargs): def make_ace(filename, temperatures=None, acer=True, xsdir=None, output_dir=None, pendf=False, error=0.001, broadr=True, heatr=True, gaspr=True, purr=True, evaluation=None, - smoothing=True, **kwargs): + smoothing=True, displacement_energy=None, **kwargs): """Generate incident neutron ACE file from an ENDF file File names can be passed to @@ -463,6 +463,9 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, indicates which evaluation should be used. smoothing : bool, optional If the smoothing option (ACER card 6) is on (True) or off (False). + displacement_energy : float, optional + Threshold displacement energy in [eV]. Used in HEATR to calculate + damage-energy cross section. When None, use NJOY defaults. **kwargs Keyword arguments passed to :func:`openmc.data.njoy.run` @@ -515,6 +518,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, # heatr if heatr: + ed = displacement_energy if displacement_energy is not None else 0 nheatr_in = nlast nheatr_local = nheatr_in + 1 tapeout[nheatr_local] = (output_dir / "heatr_local") if heatr is True \ diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 42bb958ca..7049f7d31 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 6dbb3adb2..6c19531e7 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 adb0d3dbc..0cf0141d2 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 44d6ac273..196682394 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 eded9c3fa..3a92efda2 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 e425e6f83..b48feef1b 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 aeece7acf..531184095 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 a11cfd6c8..5947540d6 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 098487c06..c769fb7d7 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 5bbda4830..d166f1247 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 7f5b498ad..cbdd13127 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 ba8658f10..5a3452422 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 ddb61385f..bf654f4f9 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 e799a340e..9e56c58e8 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 81f08fa04..ddacc2bd9 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/mesh.cpp b/src/mesh.cpp index 181af8466..89d1ff24b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -6,7 +6,8 @@ #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers #include // for ceil #include // for size_t -#include // for accumulate +#include +#include // for accumulate #include #ifdef _MSC_VER @@ -1080,13 +1081,26 @@ int StructuredMesh::get_bin(Position r) const int StructuredMesh::n_bins() const { - return std::accumulate( - shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>()); + // Bin indices are stored as 32-bit ints in the tally system. + int64_t n = 1; + for (int i = 0; i < n_dimension_; ++i) + n *= shape_[i]; + if (n > std::numeric_limits::max()) { + fatal_error(fmt::format( + "Mesh {} has too many bins ({}) for 32-bit tally indexing", id_, n)); + } + return static_cast(n); } int StructuredMesh::n_surface_bins() const { - return 4 * n_dimension_ * n_bins(); + // Surface bin indices are stored as 32-bit ints in the tally system. + int64_t n = static_cast(n_bins()) * 4 * n_dimension_; + if (n > std::numeric_limits::max()) { + fatal_error(fmt::format( + "Mesh {} has too many surface bins ({}) for tally indexing", id_, n)); + } + return static_cast(n); } tensor::Tensor StructuredMesh::count_sites( diff --git a/src/output.cpp b/src/output.cpp index 725e8d069..04a8caa60 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -621,8 +621,15 @@ void write_tallies() if (model::tallies.empty()) return; + // Tag tallies.out written during the forward solve of an adjoint run + const char* forward = + (FlatSourceDomain::solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) + ? "forward." + : ""; + // Set filename for tallies_out - std::string filename = fmt::format("{}tallies.out", settings::path_output); + std::string filename = + fmt::format("{}tallies.{}out", settings::path_output, forward); // Open the tallies.out file. std::ofstream tallies_out; diff --git a/src/particle.cpp b/src/particle.cpp index 779ae18da..10e8f213d 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 dcc771329..4a8b4d368 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 b9a1136af..96df1c5be 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 29457569b..f2a81fc1c 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/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 06c6ef14d..83128fdaa 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -30,7 +30,8 @@ namespace openmc { RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ { RandomRayVolumeEstimator::HYBRID}; bool FlatSourceDomain::volume_normalized_flux_tallies_ {false}; -bool FlatSourceDomain::adjoint_ {false}; +bool FlatSourceDomain::adjoint_requested_ {false}; +RandomRaySolve FlatSourceDomain::solve_ {RandomRaySolve::FORWARD}; bool FlatSourceDomain::fw_cadis_local_ {false}; double FlatSourceDomain::diagonal_stabilization_rho_ {1.0}; std::unordered_map>> @@ -556,7 +557,7 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const // If we are in adjoint mode of a fixed source problem, the external // source is already normalized, such that all resulting fluxes are // also normalized. - if (adjoint_) { + if (solve_ == RandomRaySolve::ADJOINT) { return 1.0; } @@ -724,7 +725,7 @@ void FlatSourceDomain::random_ray_tally() for (int i = 0; i < model::tallies.size(); i++) { Tally& tally {*model::tallies[i]}; #pragma omp parallel for - for (int bin = 0; bin < tally.n_filter_bins(); bin++) { + for (int64_t bin = 0; bin < tally.n_filter_bins(); bin++) { for (int score_idx = 0; score_idx < tally.n_scores(); score_idx++) { auto score_type = tally.scores_[score_idx]; if (score_type == SCORE_FLUX) { @@ -795,6 +796,12 @@ void FlatSourceDomain::output_to_vtk() const double z_delta = width.z / Nz; std::string filename = openmc_plot->path_plot(); + // Tag plots written during the forward solve of an adjoint run + if (solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) { + auto dot = filename.find_last_of('.'); + filename = filename.substr(0, dot) + ".forward" + filename.substr(dot); + } + // Perform sanity checks on file size uint64_t bytes = Nx * Ny * Nz * (negroups_ + 1 + 1 + 1) * sizeof(float); write_message(5, "Processing plot {}: {}... (Estimated size is {} MB)", @@ -1002,9 +1009,10 @@ void FlatSourceDomain::output_to_vtk() const void FlatSourceDomain::apply_external_source_to_source_region( int src_idx, SourceRegionHandle& srh) { - auto s = (adjoint_ && !model::adjoint_sources.empty()) - ? model::adjoint_sources[src_idx].get() - : model::external_sources[src_idx].get(); + auto s = + (solve_ == RandomRaySolve::ADJOINT && !model::adjoint_sources.empty()) + ? model::adjoint_sources[src_idx].get() + : model::external_sources[src_idx].get(); auto is = dynamic_cast(s); auto discrete = dynamic_cast(is->energy()); double strength_factor = is->strength(); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 24650afb8..00fff99a7 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -189,7 +189,7 @@ void validate_random_ray_inputs() // Validate adjoint sources /////////////////////////////////////////////////////////////////// - if (FlatSourceDomain::adjoint_ && !model::adjoint_sources.empty()) { + if (FlatSourceDomain::adjoint_requested_ && !model::adjoint_sources.empty()) { for (int i = 0; i < model::adjoint_sources.size(); i++) { Source* s = model::adjoint_sources[i].get(); @@ -289,7 +289,8 @@ void openmc_finalize_random_ray() { FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID; FlatSourceDomain::volume_normalized_flux_tallies_ = false; - FlatSourceDomain::adjoint_ = false; + FlatSourceDomain::adjoint_requested_ = false; + FlatSourceDomain::solve_ = RandomRaySolve::FORWARD; FlatSourceDomain::fw_cadis_local_ = false; FlatSourceDomain::fw_cadis_local_targets_.clear(); FlatSourceDomain::mesh_domain_map_.clear(); @@ -356,19 +357,17 @@ void RandomRaySimulation::prepare_local_fixed_sources_adjoint() } } -void RandomRaySimulation::prepare_adjoint_simulation(bool fw_adjoint) +void RandomRaySimulation::prepare_adjoint_simulation(bool from_forward) { reset_timers(); if (mpi::master) header("ADJOINT FLUX SOLVE", 3); - if (fw_adjoint) { - // Forward simulation has already been run; - // Configure the domain for adjoint simulation and - // re-initialize OpenMC general data structures - FlatSourceDomain::adjoint_ = true; - + if (from_forward) { + // The forward solve has already run. Re-initialize OpenMC's general data + // structures for the adjoint solve and derive the adjoint source from the + // forward flux. openmc_simulation_init(); prepare_fw_fixed_sources_adjoint(); @@ -603,7 +602,8 @@ void RandomRaySimulation::print_results_random_ray( } fmt::print(" Volume Estimator Type = {}\n", estimator); - std::string adjoint_true = (FlatSourceDomain::adjoint_) ? "ON" : "OFF"; + std::string adjoint_true = + (FlatSourceDomain::solve_ == RandomRaySolve::ADJOINT) ? "ON" : "OFF"; fmt::print(" Adjoint Flux Mode = {}\n", adjoint_true); std::string shape; @@ -675,60 +675,49 @@ void RandomRaySimulation::print_results_random_ray( void openmc_run_random_ray() { - ////////////////////////////////////////////////////////// - // Run forward simulation - ////////////////////////////////////////////////////////// + using namespace openmc; - // Check if adjoint calculation is needed, and if local adjoint source(s) - // are present. If an adjoint calculation is needed and no sources are - // specified, we will run a forward calculation first to calculate adjoint - // sources for global variance reduction, then perform an adjoint - // calculation later. - bool adjoint_needed = openmc::FlatSourceDomain::adjoint_; - bool fw_adjoint = openmc::model::adjoint_sources.empty() && adjoint_needed; + // Determine which solves to run. If adjoint results are requested and no + // user-defined adjoint source is present, an initial forward solve is needed + // to construct the adjoint source from the forward flux (FW-CADIS). If the + // user has defined an adjoint source, the forward solve is skipped and only + // the adjoint solve is run. + const bool run_adjoint = FlatSourceDomain::adjoint_requested_; + const bool have_adjoint_source = !model::adjoint_sources.empty(); + const bool run_forward = !(run_adjoint && have_adjoint_source); - // If we're going to do an adjoint simulation with forward-weighted adjoint - // sources afterwards, report that this is the initial forward flux solve. - if (!adjoint_needed || fw_adjoint) { - // Configure the domain for forward simulation - openmc::FlatSourceDomain::adjoint_ = false; - - if (adjoint_needed && openmc::mpi::master) - openmc::header("FORWARD FLUX SOLVE", 3); + // Set the initial solve type + if (!run_forward) { + FlatSourceDomain::solve_ = RandomRaySolve::ADJOINT; + } else if (run_adjoint) { + FlatSourceDomain::solve_ = RandomRaySolve::FORWARD_FOR_ADJOINT; } else { - // Configure domain for adjoint simulation (later) - openmc::FlatSourceDomain::adjoint_ = true; + FlatSourceDomain::solve_ = RandomRaySolve::FORWARD; } // Initialize OpenMC general data structures openmc_simulation_init(); // Validate that inputs meet requirements for random ray mode - if (openmc::mpi::master) - openmc::validate_random_ray_inputs(); + if (mpi::master) + validate_random_ray_inputs(); // Initialize Random Ray Simulation Object - openmc::RandomRaySimulation sim; + RandomRaySimulation sim; - if (!adjoint_needed || fw_adjoint) { - // Initialize fixed sources, if present + // Run the forward solve + if (run_forward) { + // When an adjoint solve follows, report this as the initial forward solve + if (run_adjoint && mpi::master) + header("FORWARD FLUX SOLVE", 3); sim.apply_fixed_sources_and_mesh_domains(); - - // Execute random ray simulation sim.simulate(); } - ////////////////////////////////////////////////////////// - // Run adjoint simulation (if enabled) - ////////////////////////////////////////////////////////// - - if (!adjoint_needed) { - return; + // Run the adjoint solve + if (run_adjoint) { + FlatSourceDomain::solve_ = RandomRaySolve::ADJOINT; + sim.prepare_adjoint_simulation(run_forward); + sim.simulate(); } - - // Setup for adjoint simulation - sim.prepare_adjoint_simulation(fw_adjoint); - - // Execute random ray simulation - sim.simulate(); } diff --git a/src/settings.cpp b/src/settings.cpp index 8ae252ae1..15a7b9c27 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -322,7 +322,7 @@ void get_run_parameters(pugi::xml_node node_base) get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies"); } if (check_for_node(random_ray_node, "adjoint")) { - FlatSourceDomain::adjoint_ = + FlatSourceDomain::adjoint_requested_ = get_node_value_bool(random_ray_node, "adjoint"); } if (check_for_node(random_ray_node, "sample_method")) { @@ -1247,6 +1247,7 @@ void read_settings_xml(pugi::xml_node root) // read weight windows from file if (check_for_node(root, "weight_windows_file")) { weight_windows_file = get_node_value(root, "weight_windows_file"); + weight_windows_on = true; } // read settings for weight windows value, this will override diff --git a/src/simulation.cpp b/src/simulation.cpp index 89aa9ca0e..03f40a726 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -12,10 +12,12 @@ #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" #include "openmc/random_lcg.h" +#include "openmc/random_ray/flat_source_domain.h" #include "openmc/settings.h" #include "openmc/source.h" #include "openmc/state_point.h" @@ -40,6 +42,7 @@ #include #include +#include #include //============================================================================== @@ -200,9 +203,12 @@ int openmc_simulation_finalize() if (settings::output_tallies && mpi::master) write_tallies(); - // If weight window generators are present in this simulation, - // write a weight windows file - if (variance_reduction::weight_windows_generators.size() > 0) { + // If weight window generators are present in this simulation, write a + // weight windows file. This is skipped during the forward solve of an + // adjoint (FW-CADIS) run, where only the adjoint-derived weight windows + // are meaningful. + if (variance_reduction::weight_windows_generators.size() > 0 && + FlatSourceDomain::solve_ != RandomRaySolve::FORWARD_FOR_ADJOINT) { openmc_weight_windows_export(); } @@ -350,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 //============================================================================== @@ -882,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); + } } } @@ -914,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; @@ -947,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 @@ -979,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); @@ -996,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 @@ -1071,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; @@ -1141,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 bee20d5c3..12bbdf84e 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)); } @@ -523,9 +530,16 @@ void FileSource::load_sites_from_file(const std::string& path) file_close(file_id); } - // Make sure particles in source file have valid types + // Make sure particles in source file have valid types. If any particle is a + // photon, electron, or positron, enable photon transport so that the + // appropriate cross sections are loaded. for (const auto& site : this->sites_) { validate_particle_type(site.particle, "FileSource"); + if (site.particle == ParticleType::photon() || + site.particle == ParticleType::electron() || + site.particle == ParticleType::positron()) { + settings::photon_transport = true; + } } } @@ -666,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/state_point.cpp b/src/state_point.cpp index da1c141a2..02c68bf60 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -22,6 +22,7 @@ #include "openmc/nuclide.h" #include "openmc/output.h" #include "openmc/particle_type.h" +#include "openmc/random_ray/flat_source_domain.h" #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/tallies/derivative.h" @@ -46,9 +47,15 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) // Determine width for zero padding int w = std::to_string(settings::n_max_batches).size(); + // Tag statepoints written during the forward solve of an adjoint run + const char* forward = + (FlatSourceDomain::solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) + ? "forward." + : ""; + // Set filename for state point - filename_ = fmt::format("{0}statepoint.{1:0{2}}.h5", settings::path_output, - simulation::current_batch, w); + filename_ = fmt::format("{0}statepoint.{3}{1:0{2}}.h5", + settings::path_output, simulation::current_batch, w, forward); } // If a file name was specified, ensure it has .h5 file extension diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index e4ea56e59..085633247 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -512,7 +512,7 @@ void Tally::set_strides() // longest stride. auto n = filters_.size(); strides_.resize(n, 0); - int stride = 1; + int64_t stride = 1; for (int i = n - 1; i >= 0; --i) { strides_[i] = stride; stride *= model::tally_filters[filters_[i]]->n_bins(); @@ -887,7 +887,7 @@ void Tally::accumulate() if (higher_moments_) { #pragma omp parallel for // filter bins (specific cell, energy bins) - for (int i = 0; i < results_.shape(0); ++i) { + for (int64_t i = 0; i < results_.shape(0); ++i) { // score bins (flux, total reaction rate, fission reaction rate, etc.) for (int j = 0; j < results_.shape(1); ++j) { double val = results_(i, j, TallyResult::VALUE) * norm; @@ -902,7 +902,7 @@ void Tally::accumulate() } else { #pragma omp parallel for // filter bins (specific cell, energy bins) - for (int i = 0; i < results_.shape(0); ++i) { + for (int64_t i = 0; i < results_.shape(0); ++i) { // score bins (flux, total reaction rate, fission reaction rate, etc.) for (int j = 0; j < results_.shape(1); ++j) { double val = results_(i, j, TallyResult::VALUE) * norm; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index b7947afab..241cbf008 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -158,7 +158,7 @@ void score_fission_delayed_dg(int i_tally, int d_bin, double score, dg_match.bins_[i_bin] = d_bin; // Determine the filter scoring index - auto filter_index = 0; + int64_t filter_index = 0; double filter_weight = 1.; for (auto i = 0; i < tally.filters().size(); ++i) { auto i_filt = tally.filters(i); @@ -449,7 +449,7 @@ void score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) (score_bin == SCORE_PROMPT_NU_FISSION && g == 0)) { // Find the filter scoring index for this filter combination - int filter_index = 0; + int64_t filter_index = 0; double filter_weight = 1.0; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); @@ -497,7 +497,7 @@ void score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) } else { // Find the filter index and weight for this filter combination - int filter_index = 0; + int64_t filter_index = 0; double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); @@ -578,8 +578,8 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) //! collision estimator. void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, - int filter_index, double filter_weight, int i_nuclide, double atom_density, - double flux) + int64_t filter_index, double filter_weight, int i_nuclide, + double atom_density, double flux) { Tally& tally {*model::tallies[i_tally]}; @@ -1112,8 +1112,8 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, //! is not used for analog tallies. void score_general_ce_analog(Particle& p, int i_tally, int start_index, - int filter_index, double filter_weight, int i_nuclide, double atom_density, - double flux) + int64_t filter_index, double filter_weight, int i_nuclide, + double atom_density, double flux) { Tally& tally {*model::tallies[i_tally]}; @@ -1615,8 +1615,8 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, //! argument is really just used for filter weights. void score_general_mg(Particle& p, int i_tally, int start_index, - int filter_index, double filter_weight, int i_nuclide, double atom_density, - double flux) + int64_t filter_index, double filter_weight, int i_nuclide, + double atom_density, double flux) { auto& tally {*model::tallies[i_tally]}; @@ -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/tallies/trigger.cpp b/src/tallies/trigger.cpp index 6c54edd5e..7515f01a9 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -28,7 +28,7 @@ KTrigger keff_trigger; //============================================================================== std::pair get_tally_uncertainty( - int i_tally, int score_index, int filter_index) + int i_tally, int score_index, int64_t filter_index) { const auto& tally {model::tallies[i_tally]}; @@ -71,7 +71,7 @@ void check_tally_triggers(double& ratio, int& tally_id, int& score) continue; const auto& results = t.results_; - for (auto filter_index = 0; filter_index < results.shape(0); + for (int64_t filter_index = 0; filter_index < results.shape(0); ++filter_index) { // Compute the tally uncertainty metrics. auto uncert_pair = diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 0614110cd..0d268ce17 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -791,7 +791,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) if (method_string == "magic") { method_ = WeightWindowUpdateMethod::MAGIC; if (settings::solver_type == SolverType::RANDOM_RAY && - FlatSourceDomain::adjoint_) { + FlatSourceDomain::adjoint_requested_) { fatal_error("Random ray weight window generation with MAGIC cannot be " "done in adjoint mode."); } @@ -800,7 +800,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) if (settings::solver_type != SolverType::RANDOM_RAY) { fatal_error("FW-CADIS can only be run in random ray solver mode."); } - FlatSourceDomain::adjoint_ = true; + FlatSourceDomain::adjoint_requested_ = true; if (check_for_node(node, "targets")) { FlatSourceDomain::fw_cadis_local_ = true; targets_ = get_node_array(node, "targets"); @@ -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/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index ce7e539ea..3eb88f023 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -6,6 +6,7 @@ set(TEST_NAMES test_math test_mcpl_stat_sum test_mesh + test_ray test_region test_tensor # Add additional unit test files here diff --git a/tests/cpp_unit_tests/test_distribution.cpp b/tests/cpp_unit_tests/test_distribution.cpp index e46088a63..479e7c8a7 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 1ad7c4b70..7467aa1fe 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/cpp_unit_tests/test_ray.cpp b/tests/cpp_unit_tests/test_ray.cpp new file mode 100644 index 000000000..c043c4d25 --- /dev/null +++ b/tests/cpp_unit_tests/test_ray.cpp @@ -0,0 +1,117 @@ +#include +#include + +#include +#include +#include + +#include "openmc/geometry.h" +#include "openmc/geometry_aux.h" +#include "openmc/ray.h" +#include "openmc/settings.h" +#include "openmc/surface.h" + +namespace { + +using Catch::Matchers::WithinAbs; + +constexpr double DISTANCE_TOLERANCE = 1.0e-12; + +class SphericalShellFixture { +public: + SphericalShellFixture() + : run_mode_(openmc::settings::run_mode), + root_universe_(openmc::model::root_universe), + n_coord_levels_(openmc::model::n_coord_levels) + { + openmc::settings::run_mode = openmc::RunMode::PLOTTING; + + constexpr auto geometry = R"( + + + + + + + )"; + + auto result = document_.load_string(geometry); + REQUIRE(result); + openmc::read_geometry_xml(document_.document_element()); + openmc::finalize_geometry(); + openmc::finalize_cell_densities(); + } + + ~SphericalShellFixture() + { + openmc::free_memory_geometry(); + openmc::free_memory_surfaces(); + openmc::model::universe_level_counts.clear(); + openmc::model::root_universe = root_universe_; + openmc::model::n_coord_levels = n_coord_levels_; + openmc::settings::run_mode = run_mode_; + } + +private: + pugi::xml_document document_; + openmc::RunMode run_mode_; + int root_universe_; + int n_coord_levels_; +}; + +class RecordingRay : public openmc::Ray { +public: + using Ray::Ray; + + void on_intersection() override + { + intersections.emplace_back( + traversal_distance_, boundary().surface_index() + 1); + } + + openmc::vector> intersections; +}; + +} // namespace + +TEST_CASE_METHOD(SphericalShellFixture, "Trace a single chord through a sphere") +{ + constexpr double impact_parameter = 1.5; + const double half_chord = + std::sqrt(4.0 - impact_parameter * impact_parameter); + + RecordingRay ray({-3.0, impact_parameter, 0.0}, {1.0, 0.0, 0.0}); + ray.trace(); + + REQUIRE(ray.intersections.size() == 2); + CHECK(ray.intersections[0].second == 2); + CHECK_THAT(ray.intersections[0].first, WithinAbs(0.0, DISTANCE_TOLERANCE)); + CHECK(ray.intersections[1].second == 2); + CHECK_THAT(ray.intersections[1].first, + WithinAbs(2.0 * half_chord, DISTANCE_TOLERANCE)); +} + +TEST_CASE_METHOD( + SphericalShellFixture, "Trace both segments of a spherical shell") +{ + constexpr double impact_parameter = 0.5; + const double outer = std::sqrt(4.0 - impact_parameter * impact_parameter); + const double inner = std::sqrt(1.0 - impact_parameter * impact_parameter); + + RecordingRay ray({-3.0, impact_parameter, 0.0}, {1.0, 0.0, 0.0}); + ray.trace(); + + REQUIRE(ray.intersections.size() == 4); + CHECK(ray.intersections[0].second == 2); + CHECK_THAT(ray.intersections[0].first, WithinAbs(0.0, DISTANCE_TOLERANCE)); + CHECK(ray.intersections[1].second == 1); + CHECK_THAT( + ray.intersections[1].first, WithinAbs(outer - inner, DISTANCE_TOLERANCE)); + CHECK(ray.intersections[2].second == 1); + CHECK_THAT( + ray.intersections[2].first, WithinAbs(outer + inner, DISTANCE_TOLERANCE)); + CHECK(ray.intersections[3].second == 2); + CHECK_THAT( + ray.intersections[3].first, WithinAbs(2.0 * outer, DISTANCE_TOLERANCE)); +} diff --git a/tests/cpp_unit_tests/test_tally.cpp b/tests/cpp_unit_tests/test_tally.cpp index 964d30cc4..03322e764 100644 --- a/tests/cpp_unit_tests/test_tally.cpp +++ b/tests/cpp_unit_tests/test_tally.cpp @@ -1,3 +1,4 @@ +#include "openmc/tallies/filter_energy.h" #include "openmc/tallies/tally.h" #include @@ -46,10 +47,34 @@ TEST_CASE("Test add/set_filter") REQUIRE(model::filter_map[cell_filter->id()] == tally->filters(0)); REQUIRE(model::filter_map[particle_filter->id()] == tally->filters(1)); - // set filters with a duplicate filter, should only add the filter to the tally once + // set filters with a duplicate filter, should only add the filter to the + // tally once filters = {cell_filter, cell_filter}; tally->set_filters(filters); REQUIRE(tally->filters().size() == 1); REQUIRE(model::filter_map[cell_filter->id()] == tally->filters(0)); +} +// Regression test for 64-bit tally filter-bin counts (mesh x groups > 2^31). +TEST_CASE("Tally filter-bin count does not overflow 32 bits") +{ + // Two energy filters whose bin counts multiply to 2.5e9, above INT32_MAX. + constexpr int64_t bins_per_filter = 50000; + + // Only the bin count matters here, so the edge values are an arbitrary ramp. + std::vector edges(bins_per_filter + 1); + for (int64_t i = 0; i < bins_per_filter + 1; ++i) + edges[i] = static_cast(i); + + Tally* tally = Tally::create(); + for (int i = 0; i < 2; ++i) { + Filter* filter = Filter::create("energy"); + dynamic_cast(filter)->set_bins(edges); + tally->add_filter(filter); + } + tally->set_strides(); + + // set_strides() previously accumulated this product in a 32-bit int. + REQUIRE(tally->n_filter_bins() == bins_per_filter * bins_per_filter); + REQUIRE(tally->n_filter_bins() > 2147483647); } \ No newline at end of file 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 000000000..e69de29bb 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 000000000..4dd4f0bd7 --- /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 000000000..58c4e306d --- /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 000000000..0f409badd --- /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 1d3b063bd..c5b6bd889 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/testing_harness.py b/tests/testing_harness.py index 8d156bd64..c1977556d 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -143,7 +143,8 @@ class TestHarness: def _cleanup(self): """Delete statepoints, tally, and test files.""" output = glob.glob('statepoint.*.h5') - output += ['tallies.out', 'results_test.dat', 'summary.h5'] + output += ['tallies.out', 'tallies.forward.out'] + output += ['results_test.dat', 'summary.h5'] output += glob.glob('volume_*.h5') for f in output: if os.path.exists(f): diff --git a/tests/unit_tests/test_data_decay.py b/tests/unit_tests/test_data_decay.py index a48c55325..f81f857e0 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 4f1880014..32be5eb89 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 14db68913..320d09eda 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 39c532c54..bc1a078b0 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 f62306e83..4cfae0df2 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 89dfc03dd..12d363376 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 f1192bd8f..796453708 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 000000000..1f27cc6f2 --- /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 266bd4976..5e617919d 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 000000000..e5f08245d --- /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_file.py b/tests/unit_tests/test_source_file.py index 3e84fbe51..f19ea74a6 100644 --- a/tests/unit_tests/test_source_file.py +++ b/tests/unit_tests/test_source_file.py @@ -145,3 +145,27 @@ def test_source_file_transport(run_in_tmpdir): # Try running OpenMC model.run() + + +def test_source_file_photon_transport(run_in_tmpdir): + # Create a source file containing a photon. Note that photon_transport is + # not explicitly enabled in the settings -- it should be turned on + # automatically because the source file contains a photon. + particle = openmc.SourceParticle(E=1.0e6, particle='photon') + openmc.write_source_file([particle], 'photon_source.h5') + + # Create simple model to use the photon source file + model = openmc.Model() + al = openmc.Material() + al.add_element('Al', 1.0) + al.set_density('g/cm3', 2.7) + sph = openmc.Sphere(r=10.0, boundary_type='vacuum') + cell = openmc.Cell(fill=al, region=-sph) + model.geometry = openmc.Geometry([cell]) + model.settings.source = openmc.FileSource(path='photon_source.h5') + model.settings.particles = 10 + model.settings.batches = 3 + model.settings.run_mode = 'fixed source' + + # Running OpenMC should succeed + model.run() diff --git a/tests/unit_tests/test_source_tokamak.py b/tests/unit_tests/test_source_tokamak.py new file mode 100644 index 000000000..f926c2007 --- /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 4e067ffe2..e4419252e 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