OpenMC/openmc/model/model.py

from __future__ import annotations
from collections.abc import Callable, Iterable, Sequence
import copy
from dataclasses import dataclass, field
from functools import cache
from pathlib import Path
import math
from numbers import Integral, Real
import random
import re
from tempfile import NamedTemporaryFile, TemporaryDirectory
from typing import Any, Protocol
import warnings

import h5py
import lxml.etree as ET
import numpy as np
from scipy.optimize import curve_fit

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, 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


# Protocol for a function that is passed to search_keff
class ModelModifier(Protocol):
    def __call__(self, val: float, **kwargs: Any) -> None:
        ...


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.

    This class can be used to store instances of :class:`openmc.Geometry`,
    :class:`openmc.Materials`, :class:`openmc.Settings`,
    :class:`openmc.Tallies`, and :class:`openmc.Plots`, thus making a complete
    model. The :meth:`Model.export_to_xml` method will export XML files for all
    attributes that have been set. If the :attr:`Model.materials` attribute is
    not set, it will attempt to create a ``materials.xml`` file based on all
    materials appearing in the geometry.

    .. versionchanged:: 0.13.0
        The model information can now be loaded in to OpenMC directly via
        openmc.lib

    Parameters
    ----------
    geometry : openmc.Geometry, optional
        Geometry information
    materials : openmc.Materials, optional
        Materials information
    settings : openmc.Settings, optional
        Settings information
    tallies : openmc.Tallies, optional
        Tallies information
    plots : openmc.Plots, optional
        Plot information
    description : str, optional
        A description of the model

    Attributes
    ----------
    geometry : openmc.Geometry
        Geometry information
    materials : openmc.Materials
        Materials information
    settings : openmc.Settings
        Settings information
    tallies : openmc.Tallies
        Tallies information
    plots : openmc.Plots
        Plot information
    description : str
        A description of the model

    """

    def __init__(
        self,
        geometry: openmc.Geometry | None = None,
        materials: openmc.Materials | None = None,
        settings: openmc.Settings | None = None,
        tallies: openmc.Tallies | None = None,
        plots: openmc.Plots | None = None,
        description: str = '',
    ):
        self.geometry = openmc.Geometry() if geometry is None else geometry
        self.materials = openmc.Materials() if materials is None else materials
        self.settings = openmc.Settings() if settings is None else settings
        self.tallies = openmc.Tallies() if tallies is None else tallies
        self.plots = openmc.Plots() if plots is None else plots
        self.description = description

    @property
    def geometry(self) -> openmc.Geometry:
        return self._geometry

    @geometry.setter
    def geometry(self, geometry):
        check_type('geometry', geometry, openmc.Geometry)
        self._geometry = geometry

    @property
    def materials(self) -> openmc.Materials:
        return self._materials

    @materials.setter
    def materials(self, materials):
        check_type('materials', materials, Iterable, openmc.Material)
        if isinstance(materials, openmc.Materials):
            self._materials = materials
        else:
            if not hasattr(self, '_materials'):
                self._materials = openmc.Materials()
            del self._materials[:]
            for mat in materials:
                self._materials.append(mat)

    @property
    def settings(self) -> openmc.Settings:
        return self._settings

    @settings.setter
    def settings(self, settings):
        check_type('settings', settings, openmc.Settings)
        self._settings = settings

    @property
    def tallies(self) -> openmc.Tallies:
        return self._tallies

    @tallies.setter
    def tallies(self, tallies):
        check_type('tallies', tallies, Iterable, openmc.Tally)
        if isinstance(tallies, openmc.Tallies):
            self._tallies = tallies
        else:
            if not hasattr(self, '_tallies'):
                self._tallies = openmc.Tallies()
            del self._tallies[:]
            for tally in tallies:
                self._tallies.append(tally)

    @property
    def plots(self) -> openmc.Plots:
        return self._plots

    @plots.setter
    def plots(self, plots):
        check_type('plots', plots, Iterable, openmc.PlotBase)
        if isinstance(plots, openmc.Plots):
            self._plots = plots
        else:
            if not hasattr(self, '_plots'):
                self._plots = openmc.Plots()
            del self._plots[:]
            for plot in plots:
                self._plots.append(plot)

    @property
    def description(self) -> str:
        return self._description

    @description.setter
    def description(self, description):
        check_type('description', description, str)
        self._description = description

    @property
    def bounding_box(self) -> openmc.BoundingBox:
        return self.geometry.bounding_box

    @property
    def is_initialized(self) -> bool:
        try:
            import openmc.lib
            return openmc.lib.is_initialized
        except ImportError:
            return False

    @property
    @cache
    def _materials_by_id(self) -> dict:
        """Dictionary mapping material ID --> material"""
        if self.materials:
            mats = self.materials
        else:
            mats = self.geometry.get_all_materials().values()
        return {mat.id: mat for mat in mats}

    @property
    @cache
    def _cells_by_id(self) -> dict:
        """Dictionary mapping cell ID --> cell"""
        cells = self.geometry.get_all_cells()
        return {cell.id: cell for cell in cells.values()}

    @property
    @cache
    def _cells_by_name(self) -> dict[int, openmc.Cell]:
        # Get the names maps, but since names are not unique, store a set for
        # each name key. In this way when the user requests a change by a name,
        # the change will be applied to all of the same name.
        result = {}
        for cell in self.geometry.get_all_cells().values():
            if cell.name not in result:
                result[cell.name] = set()
            result[cell.name].add(cell)
        return result

    @property
    @cache
    def _materials_by_name(self) -> dict[int, openmc.Material]:
        if self.materials is None:
            mats = self.geometry.get_all_materials().values()
        else:
            mats = self.materials
        result = {}
        for mat in mats:
            if mat.name not in result:
                result[mat.name] = set()
            result[mat.name].add(mat)
        return result

    # TODO: This should really get incorporated in lower-level calls to
    # get_all_materials, but right now it requires information from the Model object
    def _get_all_materials(self) -> dict[int, openmc.Material]:
        """Get all materials including those in DAGMC universes

        Returns
        -------
        dict
            Dictionary mapping material ID to material instances
        """
        # Get all materials from the Geometry object
        materials = self.geometry.get_all_materials()

        # Account for materials in DAGMC universes
        for cell in self.geometry.get_all_cells().values():
            if isinstance(cell.fill, openmc.DAGMCUniverse):
                names = cell.fill.material_names
                materials.update({
                    mat.id: mat for mat in self.materials if mat.name in names
                })

        return materials

    def add_kinetics_parameters_tallies(self, num_groups: int | None = None):
        """Add tallies for calculating kinetics parameters using the IFP method.

        This method adds tallies to the model for calculating two kinetics
        parameters, the generation time and the effective delayed neutron
        fraction (beta effective). After a model is run, these parameters can be
        determined through the :meth:`openmc.StatePoint.ifp_results` method.

        Parameters
        ----------
        num_groups : int, optional
            Number of precursor groups to filter the delayed neutron fraction.
            If None, only the total effective delayed neutron fraction is
            tallied.

        """
        if not any('ifp-time-numerator' in t.scores for t in self.tallies):
            gen_time_tally = openmc.Tally(name='IFP time numerator')
            gen_time_tally.scores = ['ifp-time-numerator']
            self.tallies.append(gen_time_tally)
        if not any('ifp-beta-numerator' in t.scores for t in self.tallies):
            beta_tally = openmc.Tally(name='IFP beta numerator')
            beta_tally.scores = ['ifp-beta-numerator']
            if num_groups is not None:
                beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))]
            self.tallies.append(beta_tally)
        if not any('ifp-denominator' in t.scores for t in self.tallies):
            denom_tally = openmc.Tally(name='IFP denominator')
            denom_tally.scores = ['ifp-denominator']
            self.tallies.append(denom_tally)

    # TODO: This should also be incorporated into lower-level calls in
    # settings.py, but it requires information about the tallies currently
    # on the active Model
    def _assign_fw_cadis_tally_IDs(self):
        # Verify that all tallies assigned as targets on WeightWindowGenerators
        # exist within model.tallies. If this is the case, convert the .targets
        # attribute of each WeightWindowGenerator to a sequence of tally IDs.
        if len(self.settings.weight_window_generators) == 0:
            return

        # List of valid tally IDs
        reference_tally_ids = np.asarray([tal.id for tal in self.tallies])

        for wwg in self.settings.weight_window_generators:
            # Only proceeds if the "targets" attribute is an openmc.Tallies,
            # which means it hasn't been checked against model.tallies.
            if isinstance(wwg.targets, openmc.Tallies):
                id_vec = []
                for tal in wwg.targets:
                    # check against model tallies for equivalence
                    id_next = None
                    for reference_tal in self.tallies:
                        if tal == reference_tal:
                            id_next = reference_tal.id
                            break

                    if id_next is None:
                        raise RuntimeError(
                            f'Local FW-CADIS target tally {tal.id} not found on model.tallies!')
                    else:
                        id_vec.append(id_next)

                wwg.targets = id_vec

            elif isinstance(wwg.targets, np.ndarray):
                invalid = wwg.targets[~np.isin(wwg.targets, reference_tally_ids)]
                if len(invalid) > 0:
                    raise RuntimeError(
                        f'Local FW-CADIS target tally IDs {invalid} not found on model.tallies!')

    @classmethod
    def from_xml(
        cls,
        geometry: PathLike = "geometry.xml",
        materials: PathLike = "materials.xml",
        settings: PathLike = "settings.xml",
        tallies: PathLike = "tallies.xml",
        plots: PathLike = "plots.xml",
    ) -> Model:
        """Create model from existing XML files

        Parameters
        ----------
        geometry : PathLike
            Path to geometry.xml file
        materials : PathLike
            Path to materials.xml file
        settings : PathLike
            Path to settings.xml file
        tallies : PathLike
            Path to tallies.xml file

            .. versionadded:: 0.13.0
        plots : PathLike
            Path to plots.xml file

            .. versionadded:: 0.13.0

        Returns
        -------
        openmc.model.Model
            Model created from XML files

        """
        materials = openmc.Materials.from_xml(materials)
        geometry = openmc.Geometry.from_xml(geometry, materials)
        settings = openmc.Settings.from_xml(settings)
        tallies = openmc.Tallies.from_xml(
            tallies) if Path(tallies).exists() else None
        plots = openmc.Plots.from_xml(plots) if Path(plots).exists() else None
        return cls(geometry, materials, settings, tallies, plots)

    @classmethod
    def from_model_xml(cls, path: PathLike = "model.xml") -> Model:
        """Create model from single XML file

        .. versionadded:: 0.13.3

        Parameters
        ----------
        path : PathLike
            Path to model.xml file
        """
        parser = ET.XMLParser(huge_tree=True)
        tree = ET.parse(path, parser=parser)
        root = tree.getroot()

        model = cls()

        desc_elem = root.find('description')
        if desc_elem is not None and desc_elem.text:
            model.description = desc_elem.text

        meshes = {}
        model.settings = openmc.Settings.from_xml_element(
            root.find('settings'), meshes)
        model.materials = openmc.Materials.from_xml_element(
            root.find('materials'))
        model.geometry = openmc.Geometry.from_xml_element(
            root.find('geometry'), model.materials)

        if root.find('tallies') is not None:
            model.tallies = openmc.Tallies.from_xml_element(
                root.find('tallies'), meshes)

        if root.find('plots') is not None:
            model.plots = openmc.Plots.from_xml_element(root.find('plots'))

        return model

    def init_lib(
        self,
        threads: int | None = None,
        geometry_debug: bool = False,
        restart_file: PathLike | None = None,
        tracks: bool = False,
        output: bool = True,
        event_based: bool | None = None,
        intracomm=None,
        directory: PathLike | None = None,
    ):
        """Initializes the model in memory via the C API

        .. versionadded:: 0.13.0

        Parameters
        ----------
        threads : int, optional
            Number of OpenMP threads. If OpenMC is compiled with OpenMP
            threading enabled, the default is implementation-dependent but is
            usually equal to the number of hardware threads available
            (or a value set by the :envvar:`OMP_NUM_THREADS` environment
            variable).
        geometry_debug : bool, optional
            Turn on geometry debugging during simulation. Defaults to False.
        restart_file : PathLike, optional
            Path to restart file to use
        tracks : bool, optional
            Enables the writing of particles tracks. The number of particle
            tracks written to tracks.h5 is limited to 1000 unless
            Settings.max_tracks is set. Defaults to False.
        output : bool
            Capture OpenMC output from standard out
        event_based : None or bool, optional
            Turns on event-based parallelism if True. If None, the value in
            the Settings will be used.
        intracomm : mpi4py.MPI.Intracomm or None, optional
            MPI intracommunicator
        directory : PathLike or None, optional
            Directory to write XML files to. Defaults to None.
        """

        import openmc.lib

        # TODO: right now the only way to set most of the above parameters via
        # the C API are at initialization time despite use-cases existing to
        # set them for individual runs. For now this functionality is exposed
        # where it exists (here in init), but in the future the functionality
        # should be exposed so that it can be accessed via model.run(...)

        args = _process_CLI_arguments(
            volume=False, geometry_debug=geometry_debug,
            restart_file=restart_file, threads=threads, tracks=tracks,
            event_based=event_based, path_input=directory)

        # Args adds the openmc_exec command in the first entry; remove it
        args = args[1:]

        self.finalize_lib()

        # The Model object needs to be aware of the communicator so it can
        # use it in certain cases, therefore lets store the communicator
        if intracomm is not None:
            self._intracomm = intracomm
        else:
            self._intracomm = DummyCommunicator()

        if self._intracomm.rank == 0:
            if directory is not None:
                self.export_to_xml(directory=directory)
            else:
                self.export_to_xml()
        self._intracomm.barrier()

        # We cannot pass DummyCommunicator to openmc.lib.init so pass instead
        # the user-provided intracomm which will either be None or an mpi4py
        # communicator
        openmc.lib.init(args=args, intracomm=intracomm, output=output)

    def sync_dagmc_universes(self):
        """Synchronize all DAGMC universes in the current geometry.

        This method iterates over all DAGMC universes in the geometry and
        synchronizes their cells with the current material assignments. Requires
        that the model has been initialized via :meth:`Model.init_lib`.
        Synchronized DAGMC cells can then be edited and exported as nested
        `<cell>` overrides inside each `<dagmc_universe>` element.

        .. versionadded:: 0.15.1

        """
        if self.is_initialized:
            if self.materials:
                materials = self.materials
            else:
                materials = list(self.geometry.get_all_materials().values())
            for univ in self.geometry.get_all_universes().values():
                if isinstance(univ, openmc.DAGMCUniverse):
                    univ.sync_dagmc_cells(materials)
        else:
            raise ValueError("The model must be initialized before calling "
                             "this method")

    def finalize_lib(self):
        """Finalize simulation and free memory allocated for the C API

        .. versionadded:: 0.13.0

        """

        import openmc.lib

        openmc.lib.finalize()

    def deplete(
        self,
        method: str = "cecm",
        final_step: bool = True,
        operator_kwargs: dict | None = None,
        directory: PathLike = ".",
        output: bool = True,
        **integrator_kwargs,
    ):
        """Deplete model using specified timesteps/power

        .. versionchanged:: 0.13.0
            The *final_step*, *operator_kwargs*, *directory*, and *output*
            arguments were added.

        Parameters
        ----------
        timesteps : iterable of float or iterable of tuple
            Array of timesteps. Note that values are not cumulative. The units are
            specified by the `timestep_units` argument when `timesteps` is an
            iterable of float. Alternatively, units can be specified for each step
            by passing an iterable of (value, unit) tuples.
        method : str
             Integration method used for depletion (e.g., 'cecm', 'predictor').
             Defaults to 'cecm'.
        final_step : bool, optional
            Indicate whether or not a transport solve should be run at the end
            of the last timestep. Defaults to running this transport solve.
        operator_kwargs : dict
            Keyword arguments passed to the depletion operator initializer
            (e.g., :func:`openmc.deplete.Operator`)
        directory : PathLike, optional
            Directory to write XML files to. If it doesn't exist already, it
            will be created. Defaults to the current working directory
        output : bool
            Capture OpenMC output from standard out
        integrator_kwargs : dict
            Remaining keyword arguments passed to the depletion integrator
            (e.g., :class:`openmc.deplete.CECMIntegrator`).

        """

        if operator_kwargs is None:
            op_kwargs = {}
        elif isinstance(operator_kwargs, dict):
            op_kwargs = operator_kwargs
        else:
            raise ValueError("operator_kwargs must be a dict or None")

        # Import openmc.deplete here so the Model can be used even if the
        # shared library is unavailable.
        import openmc.deplete as dep

        # Store whether or not the library was initialized when we started
        started_initialized = self.is_initialized

        with change_directory(directory):
            with openmc.lib.quiet_dll(output):
                # TODO: Support use of IndependentOperator too
                depletion_operator = dep.CoupledOperator(self, **op_kwargs)

            # Tell depletion_operator.finalize NOT to clear C API memory when
            # it is done
            depletion_operator.cleanup_when_done = False

            # Set up the integrator
            check_value('method', method,
                        dep.integrators.integrator_by_name.keys())
            integrator_class = dep.integrators.integrator_by_name[method]
            integrator = integrator_class(depletion_operator, **integrator_kwargs)

            # Now perform the depletion
            with openmc.lib.quiet_dll(output):
                integrator.integrate(final_step)

            # Now make the python Materials match the C API material data
            for mat_id, mat in self._materials_by_id.items():
                if mat.depletable:
                    # Get the C data
                    c_mat = openmc.lib.materials[mat_id]
                    nuclides, densities = c_mat._get_densities()
                    # And now we can remove isotopes and add these ones in
                    mat.nuclides.clear()
                    for nuc, density in zip(nuclides, densities):
                        mat.add_nuclide(nuc, density)
                    mat.set_density('atom/b-cm', sum(densities))

            # If we didnt start intialized, we should cleanup after ourselves
            if not started_initialized:
                depletion_operator.cleanup_when_done = True
                depletion_operator.finalize()

    def _link_geometry_to_filters(self):
        """Establishes a link between distribcell filters and the geometry"""
        for tally in self.tallies:
            for f in tally.filters:
                if isinstance(f, openmc.DistribcellFilter):
                    f._geometry = self.geometry

    def export_to_xml(self, directory: PathLike = '.', remove_surfs: bool = False,
                      nuclides_to_ignore: Iterable[str] | None = None):
        """Export model to separate XML files.

        Parameters
        ----------
        directory : PathLike
            Directory to write XML files to. If it doesn't exist already, it
            will be created.
        remove_surfs : bool
            Whether or not to remove redundant surfaces from the geometry when
            exporting.

            .. versionadded:: 0.13.1
        nuclides_to_ignore : list of str
            Nuclides to ignore when exporting to XML.

        """
        # Create directory if required
        d = Path(directory)
        if not d.is_dir():
            d.mkdir(parents=True, exist_ok=True)

        self._assign_fw_cadis_tally_IDs()
        self.settings.export_to_xml(d)
        self.geometry.export_to_xml(d, remove_surfs=remove_surfs)

        # If a materials collection was specified, export it. Otherwise, look
        # for all materials in the geometry and use that to automatically build
        # a collection.
        if self.materials:
            self.materials.export_to_xml(d, nuclides_to_ignore=nuclides_to_ignore)
        else:
            materials = openmc.Materials(self.geometry.get_all_materials()
                                         .values())
            materials.export_to_xml(d, nuclides_to_ignore=nuclides_to_ignore)

        if self.tallies:
            self.tallies.export_to_xml(d)
        if self.plots:
            self.plots.export_to_xml(d)

        self._link_geometry_to_filters()

    def export_to_model_xml(self, path: PathLike = 'model.xml', remove_surfs: bool = False,
                            nuclides_to_ignore: Iterable[str] | None = None):
        """Export model to a single XML file.

        .. versionadded:: 0.13.3

        Parameters
        ----------
        path : str or PathLike
            Location of the XML file to write (default is 'model.xml'). Can be a
            directory or file path.
        remove_surfs : bool
            Whether or not to remove redundant surfaces from the geometry when
            exporting.
        nuclides_to_ignore : list of str
            Nuclides to ignore when exporting to XML.

        """
        xml_path = Path(path)
        # if the provided path doesn't end with the XML extension, assume the
        # input path is meant to be a directory. If the directory does not
        # exist, create it and place a 'model.xml' file there.
        if not str(xml_path).endswith('.xml'):
            if not xml_path.exists():
                xml_path.mkdir(parents=True, exist_ok=True)
            elif not xml_path.is_dir():
                raise FileExistsError(f"File exists and is not a directory: '{xml_path}'")
            xml_path /= 'model.xml'
        # if this is an XML file location and the file's parent directory does
        # not exist, create it before continuing
        elif not xml_path.parent.exists():
            xml_path.parent.mkdir(parents=True, exist_ok=True)

        if remove_surfs:
            warnings.warn("remove_surfs kwarg will be deprecated soon, please "
                          "set the Geometry.merge_surfaces attribute instead.")
            self.geometry.merge_surfaces = True

        # Link FW-CADIS WeightWindowGenerator target tallies, if present
        self._assign_fw_cadis_tally_IDs()

        # provide a memo to track which meshes have been written
        mesh_memo = set()
        settings_element = self.settings.to_xml_element(mesh_memo)
        geometry_element = self.geometry.to_xml_element()

        xml.clean_indentation(geometry_element, level=1)
        xml.clean_indentation(settings_element, level=1)

        # If a materials collection was specified, export it. Otherwise, look
        # for all materials in the geometry and use that to automatically build
        # a collection.
        if self.materials:
            materials = self.materials
        else:
            materials = openmc.Materials(self.geometry.get_all_materials()
                                         .values())

        with open(xml_path, 'w', encoding='utf-8', errors='xmlcharrefreplace') as fh:
            # write the XML header
            fh.write("<?xml version='1.0' encoding='utf-8'?>\n")
            fh.write("<model>\n")
            if self.description:
                description_element = ET.Element('description')
                description_element.text = self.description
                fh.write("  ")
                fh.write(ET.tostring(description_element, encoding="unicode"))
                fh.write("\n")
            # Write the materials collection to the open XML file first.
            # This will write the XML header also
            materials._write_xml(fh, False, level=1,
                                 nuclides_to_ignore=nuclides_to_ignore)
            # Write remaining elements as a tree
            fh.write(ET.tostring(geometry_element, encoding="unicode"))
            fh.write(ET.tostring(settings_element, encoding="unicode"))

            if self.tallies:
                tallies_element = self.tallies.to_xml_element(mesh_memo)
                xml.clean_indentation(
                    tallies_element, level=1, trailing_indent=self.plots)
                fh.write(ET.tostring(tallies_element, encoding="unicode"))
            if self.plots:
                plots_element = self.plots.to_xml_element()
                xml.clean_indentation(
                    plots_element, level=1, trailing_indent=False)
                fh.write(ET.tostring(plots_element, encoding="unicode"))
            fh.write("</model>\n")

        self._link_geometry_to_filters()

    def import_properties(self, filename: PathLike):
        """Import physical properties

        .. versionchanged:: 0.13.0
            This method now updates values as loaded in memory with the C API

        Parameters
        ----------
        filename : PathLike
            Path to properties HDF5 file

        See Also
        --------
        openmc.lib.export_properties

        """
        import openmc.lib

        cells = self.geometry.get_all_cells()
        materials = self.geometry.get_all_materials()

        with h5py.File(filename, 'r') as fh:
            cells_group = fh['geometry/cells']

            # Make sure number of cells matches
            n_cells = fh['geometry'].attrs['n_cells']
            if n_cells != len(cells):
                raise ValueError("Number of cells in properties file doesn't "
                                 "match current model.")

            # Update temperatures and densities for cells filled with materials
            for name, group in cells_group.items():
                cell_id = int(name.split()[1])
                cell = cells[cell_id]
                if cell.fill_type in ('material', 'distribmat'):
                    temperature = group['temperature'][()]
                    cell.temperature = temperature
                    if self.is_initialized:
                        lib_cell = openmc.lib.cells[cell_id]
                        if temperature.size > 1:
                            for i, T in enumerate(temperature):
                                lib_cell.set_temperature(T, i)
                        else:
                            lib_cell.set_temperature(temperature[0])

                    if group['density']:
                      density = group['density'][()]
                      if density.size > 1:
                          cell.density = [rho for rho in density]
                      else:
                          cell.density = density
                      if self.is_initialized:
                          lib_cell = openmc.lib.cells[cell_id]
                          if density.size > 1:
                              for i, rho in enumerate(density):
                                  lib_cell.set_density(rho, i)
                          else:
                              lib_cell.set_density(density[0])

            # Make sure number of materials matches
            mats_group = fh['materials']
            n_cells = mats_group.attrs['n_materials']
            if n_cells != len(materials):
                raise ValueError("Number of materials in properties file "
                                 "doesn't match current model.")

            # Update material densities
            for name, group in mats_group.items():
                mat_id = int(name.split()[1])
                atom_density = group.attrs['atom_density']
                materials[mat_id].set_density('atom/b-cm', atom_density)
                if self.is_initialized:
                    C_mat = openmc.lib.materials[mat_id]
                    C_mat.set_density(atom_density, 'atom/b-cm')

    def run(
        self,
        particles: int | None = None,
        threads: int | None = None,
        geometry_debug: bool = False,
        restart_file: PathLike | None = None,
        tracks: bool = False,
        output: bool = True,
        cwd: PathLike = ".",
        openmc_exec: PathLike = "openmc",
        mpi_args: Iterable[str] = None,
        event_based: bool | None = None,
        export_model_xml: bool = True,
        apply_tally_results: bool = False,
        **export_kwargs,
    ) -> Path:
        """Run OpenMC

        If the C API has been initialized, then the C API is used, otherwise,
        this method creates the XML files and runs OpenMC via a system call. In
        both cases this method returns the path to the last statepoint file
        generated.

        .. versionchanged:: 0.12
            Instead of returning the final k-effective value, this function now
            returns the path to the final statepoint written.

        .. versionchanged:: 0.13.0
            This method can utilize the C API for execution

        Parameters
        ----------
        particles : int, optional
            Number of particles to simulate per generation.
        threads : int, optional
            Number of OpenMP threads. If OpenMC is compiled with OpenMP
            threading enabled, the default is implementation-dependent but is
            usually equal to the number of hardware threads available (or a
            value set by the :envvar:`OMP_NUM_THREADS` environment variable).
        geometry_debug : bool, optional
            Turn on geometry debugging during simulation. Defaults to False.
        restart_file : str or PathLike
            Path to restart file to use
        tracks : bool, optional
            Enables the writing of particles tracks. The number of particle
            tracks written to tracks.h5 is limited to 1000 unless
            Settings.max_tracks is set. Defaults to False.
        output : bool, optional
            Capture OpenMC output from standard out
        cwd : PathLike, optional
            Path to working directory to run in. Defaults to the current working
            directory.
        openmc_exec : str, optional
            Path to OpenMC executable. Defaults to 'openmc'.
        mpi_args : list of str, optional
            MPI execute command and any additional MPI arguments to pass, e.g.
            ['mpiexec', '-n', '8'].
        event_based : None or bool, optional
            Turns on event-based parallelism if True. If None, the value in the
            Settings will be used.
        export_model_xml : bool, optional
            Exports a single model.xml file rather than separate files. Defaults
            to True.

            .. versionadded:: 0.13.3
        apply_tally_results : bool
            Whether to apply results of the final statepoint file to the
            model's tally objects.

            .. versionadded:: 0.15.1
        **export_kwargs
            Keyword arguments passed to either :meth:`Model.export_to_model_xml`
            or :meth:`Model.export_to_xml`.

        Returns
        -------
        Path
            Path to the last statepoint written by this run (None if no
            statepoint was written)

        """

        # Setting tstart here ensures we don't pick up any pre-existing
        # statepoint files in the output directory -- just in case there are
        # differences between the system clock and the filesystem, we get the
        # time of a just-created temporary file
        with NamedTemporaryFile() as fp:
            tstart = Path(fp.name).stat().st_mtime
        last_statepoint = None

        # Operate in the provided working directory
        with change_directory(cwd):
            if self.is_initialized:
                # Handle the run options as applicable
                # First dont allow ones that must be set via init
                for arg_name, arg, default in zip(
                    ['threads', 'geometry_debug', 'restart_file', 'tracks'],
                    [threads, geometry_debug, restart_file, tracks],
                    [None, False, None, False]
                ):
                    if arg != default:
                        msg = f"{arg_name} must be set via Model.is_initialized(...)"
                        raise ValueError(msg)

                init_particles = openmc.lib.settings.particles
                if particles is not None:
                    if isinstance(particles, Integral) and particles > 0:
                        openmc.lib.settings.particles = particles

                init_event_based = openmc.lib.settings.event_based
                if event_based is not None:
                    openmc.lib.settings.event_based = event_based

                # Then run using the C API
                openmc.lib.run(output)

                # Reset changes for the openmc.run kwargs handling
                openmc.lib.settings.particles = init_particles
                openmc.lib.settings.event_based = init_event_based

            else:
                # Then run via the command line
                if export_model_xml:
                    self.export_to_model_xml(**export_kwargs)
                else:
                    self.export_to_xml(**export_kwargs)
                path_input = export_kwargs.get("path", None)
                openmc.run(particles, threads, geometry_debug, restart_file,
                           tracks, output, Path('.'), openmc_exec, mpi_args,
                           event_based, path_input)

            # Get output directory and return the last statepoint written
            if self.settings.output and 'path' in self.settings.output:
                output_dir = Path(self.settings.output['path'])
            else:
                output_dir = Path.cwd()
            for sp in output_dir.glob('statepoint.*.h5'):
                mtime = sp.stat().st_mtime
                if mtime >= tstart:  # >= allows for poor clock resolution
                    tstart = mtime
                    last_statepoint = sp

        if apply_tally_results:
            self.apply_tally_results(last_statepoint)

        return last_statepoint

    def calculate_volumes(
        self,
        threads: int | None = None,
        output: bool = True,
        cwd: PathLike = ".",
        openmc_exec: PathLike = "openmc",
        mpi_args: list[str] | None = None,
        apply_volumes: bool = True,
        export_model_xml: bool = True,
        **export_kwargs,
    ):
        """Runs an OpenMC stochastic volume calculation and, if requested,
        applies volumes to the model

        .. versionadded:: 0.13.0

        Parameters
        ----------
        threads : int, optional
            Number of OpenMP threads. If OpenMC is compiled with OpenMP
            threading enabled, the default is implementation-dependent but is
            usually equal to the number of hardware threads available (or a
            value set by the :envvar:`OMP_NUM_THREADS` environment variable).
            This currenty only applies to the case when not using the C API.
        output : bool, optional
            Capture OpenMC output from standard out
        openmc_exec : str, optional
            Path to OpenMC executable. Defaults to 'openmc'.
            This only applies to the case when not using the C API.
        mpi_args : list of str, optional
            MPI execute command and any additional MPI arguments to pass,
            e.g. ['mpiexec', '-n', '8'].
            This only applies to the case when not using the C API.
        cwd : str, optional
            Path to working directory to run in. Defaults to the current
            working directory.
        apply_volumes : bool, optional
            Whether apply the volume calculation results from this calculation
            to the model. Defaults to applying the volumes.
        export_model_xml : bool, optional
            Exports a single model.xml file rather than separate files. Defaults
            to True.
        **export_kwargs
            Keyword arguments passed to either :meth:`Model.export_to_model_xml`
            or :meth:`Model.export_to_xml`.

        """

        if len(self.settings.volume_calculations) == 0:
            # Then there is no volume calculation specified
            raise ValueError("The Settings.volume_calculations attribute must"
                             " be specified before executing this method!")

        with change_directory(cwd):
            if self.is_initialized:
                if threads is not None:
                    msg = "Threads must be set via Model.is_initialized(...)"
                    raise ValueError(msg)
                if mpi_args is not None:
                    msg = "The MPI environment must be set otherwise such as" \
                        "with the call to mpi4py"
                    raise ValueError(msg)

                # Compute the volumes
                openmc.lib.calculate_volumes(output)

            else:
                if export_model_xml:
                    self.export_to_model_xml(**export_kwargs)
                else:
                    self.export_to_xml(**export_kwargs)
                path_input = export_kwargs.get("path", None)
                openmc.calculate_volumes(
                    threads=threads, output=output, openmc_exec=openmc_exec,
                    mpi_args=mpi_args, path_input=path_input
                )

            # Now we apply the volumes
            if apply_volumes:
                # Load the results and add them to the model
                for i, vol_calc in enumerate(self.settings.volume_calculations):
                    vol_calc.load_results(f"volume_{i + 1}.h5")
                    # First add them to the Python side
                    if vol_calc.domain_type == "material" and self.materials:
                        for material in self.materials:
                            if material.id in vol_calc.volumes:
                                material.add_volume_information(vol_calc)
                    else:
                        self.geometry.add_volume_information(vol_calc)

                    # And now repeat for the C API
                    if self.is_initialized and vol_calc.domain_type == 'material':
                        # Then we can do this in the C API
                        for domain_id in vol_calc.ids:
                            openmc.lib.materials[domain_id].volume = \
                                vol_calc.volumes[domain_id].n


    def _set_plot_defaults(
        self,
        origin: Sequence[float] | None,
        width: Sequence[float] | None,
        pixels: int | Sequence[int],
        basis: str
    ):
        _check_pixels(pixels)

        x, y, _ = _BASIS_INDICES[basis]

        bb = self.bounding_box
        # checks to see if bounding box contains -inf or inf values
        if np.isinf(bb.extent[basis]).any():
            if origin is None:
                origin = (0, 0, 0)
            if width is None:
                width = (10, 10)
        else:
            if origin is None:
                # if nan values in the bb.center they get replaced with 0.0
                # this happens when the bounding_box contains inf values
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore", RuntimeWarning)
                    origin = np.nan_to_num(bb.center)
            if width is None:
                bb_width = bb.width
                width = (bb_width[x], bb_width[y])

        if isinstance(pixels, int):
            aspect_ratio = width[0] / width[1]
            pixels_y = math.sqrt(pixels / aspect_ratio)
            pixels = (int(pixels / pixels_y), int(pixels_y))

        return origin, width, pixels

    def id_map(
        self,
        origin: Sequence[float] | None = None,
        width: Sequence[float] | None = None,
        pixels: int | Sequence[int] = 40000,
        basis: str = 'xy',
        color_overlaps: bool = False,
        **init_kwargs
    ) -> np.ndarray:
        """Generate an ID map for domains based on the plot parameters

        If the model is not yet initialized, it will be initialized with
        openmc.lib. If the model is initialized, the model will remain
        initialized after this method call exits.

        .. versionadded:: 0.15.3

        Parameters
        ----------
        origin : Sequence[float], optional
            Origin of the plot. If unspecified, this argument defaults to the
            center of the bounding box if the bounding box does not contain inf
            values for the provided basis, otherwise (0.0, 0.0, 0.0).
        width : Sequence[float], optional
            Width of the plot. If unspecified, this argument defaults to the
            width of the bounding box if the bounding box does not contain inf
            values for the provided basis, otherwise (10.0, 10.0).
        pixels : int | Sequence[int], optional
            If an iterable of ints is provided then this directly sets the
            number of pixels to use in each basis direction. If a single int is
            provided then this sets the total number of pixels in the plot and
            the number of pixels in each basis direction is calculated from this
            total and the image aspect ratio based on the width argument.
        basis : {'xy', 'yz', 'xz'}, optional
            Basis of the plot.
        color_overlaps : bool, optional
            Whether to assign unique IDs (-3) to overlapping regions. If False,
            overlapping regions will be assigned the ID of the lowest-numbered
            cell that occupies that region. Defaults to False.
        **init_kwargs
            Keyword arguments passed to :meth:`Model.init_lib`.

        Returns
        -------
        id_map : numpy.ndarray
            A NumPy array with shape (vertical pixels, horizontal pixels, 3) of
            OpenMC property IDs with dtype int32. The last dimension of the
            array contains cell IDs, cell instances, and material IDs (in that
            order).
        """
        ids, _ = self.slice_data(
            origin=origin,
            width=width,
            pixels=pixels,
            basis=basis,
            show_overlaps=color_overlaps,
            level=-1,
            include_properties=False,
            **init_kwargs,
        )
        return ids

    def slice_data(
        self,
        origin: Sequence[float] | None = None,
        width: Sequence[float] | None = None,
        pixels: int | Sequence[int] = 40000,
        basis: str = 'xy',
        u_span: Sequence[float] | None = None,
        v_span: Sequence[float] | None = None,
        show_overlaps: bool = False,
        level: int = -1,
        filter: openmc.Filter | None = None,
        include_properties: bool = True,
        **init_kwargs
    ) -> tuple[np.ndarray, np.ndarray | None]:
        """Generate geometry and property data for a 2D plot slice.

        This method combines the functionality of :meth:`id_map` and property
        mapping into a single call, avoiding duplicate geometry lookups. It also
        supports filter bin index lookup for tally visualization.

        .. versionadded:: 0.16.0

        Parameters
        ----------
        origin : Sequence[float], optional
            Origin of the plot. If unspecified, this argument defaults to the
            center of the bounding box if the bounding box does not contain inf
            values for the provided basis, otherwise (0.0, 0.0, 0.0).
        width : Sequence[float], optional
            Width of the plot. If unspecified, this argument defaults to the
            width of the bounding box if the bounding box does not contain inf
            values for the provided basis, otherwise (10.0, 10.0).
        pixels : int | Sequence[int], optional
            If an iterable of ints is provided then this directly sets the
            number of pixels to use in each basis direction. If a single int is
            provided then this sets the total number of pixels in the plot and
            the number of pixels in each basis direction is calculated from this
            total and the image aspect ratio based on the width argument.
        basis : {'xy', 'yz', 'xz'}, optional
            Basis of the plot.
        u_span : Sequence[float], optional
            Full-width span vector for an oriented slice (3 values). Mutually
            exclusive with width.
        v_span : Sequence[float], optional
            Full-height span vector for an oriented slice (3 values). Mutually
            exclusive with width.
        show_overlaps : bool, optional
            Whether to identify and assign unique IDs (-3) to overlapping
            regions. If False, overlapping regions will be assigned the ID of
            the lowest-numbered cell that occupies that region. Defaults to
            False.
        level : int, optional
            Universe level to plot (-1 for deepest). Defaults to -1.
        filter : openmc.Filter, optional
            If provided, the information for each pixel also includes an index
            in the filter corresponding to the pixel position.
        include_properties : bool, optional
            Whether to include temperature/density data. Defaults to True.
        **init_kwargs
            Keyword arguments passed to :meth:`Model.init_lib`.

        Returns
        -------
        geom_data : numpy.ndarray
            Shape (v_res, h_res, 3) or (v_res, h_res, 4) int32 array. Contains
            [cell_id, cell_instance, material_id] when no filter, or [cell_id,
            cell_instance, material_id, filter_bin] with filter.
        property_data : numpy.ndarray or None
            Shape (v_res, h_res, 2) float64 array with [temperature, density],
            or None if include_properties=False.
        """
        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.")

        if u_span is not None or v_span is not None:
            if u_span is None or v_span is None:
                raise ValueError("Both u_span and v_span must be provided.")
            if origin is None:
                origin = (0.0, 0.0, 0.0)
            if isinstance(pixels, int):
                u_norm = np.linalg.norm(u_span)
                v_norm = np.linalg.norm(v_span)
                aspect_ratio = u_norm / v_norm
                pixels_y = math.sqrt(pixels / aspect_ratio)
                pixels = (int(pixels / pixels_y), int(pixels_y))
        else:
            origin, width, pixels = self._set_plot_defaults(
                origin, width, pixels, basis)

        # Silence output by default. Also set arguments to start in volume
        # calculation mode to avoid loading cross sections
        init_kwargs.setdefault('output', False)
        init_kwargs.setdefault('args', ['-c'])

        # If filter does not already appear in the model, temporarily add a
        # tally with the filter
        original_length = len(self.tallies)
        if filter is not None:
            filter_ids = {f.id for t in self.tallies for f in t.filters}
            if filter.id not in filter_ids:
                # Create temporary tally while preserving ID assignment
                next_id = openmc.Tally.next_id
                temp_tally = openmc.Tally()
                temp_tally.filters = [filter]
                temp_tally.scores = ['flux']
                self.tallies.append(temp_tally)
                openmc.Tally.used_ids.remove(temp_tally.id)
                openmc.Tally.next_id = next_id

        with openmc.lib.TemporarySession(self, **init_kwargs):
            geom_data, property_data = openmc.lib.slice_data(
                origin=origin,
                width=width,
                basis=basis,
                u_span=u_span,
                v_span=v_span,
                pixels=pixels,
                show_overlaps=show_overlaps,
                level=level,
                filter=filter,
                include_properties=include_properties,
            )

        # If filter was temporarily added, remove it
        if len(self.tallies) > original_length:
            self.tallies.pop()

        return geom_data, property_data

    @add_plot_params
    def plot(
        self,
        origin: Sequence[float] | None = None,
        width: Sequence[float] | None = None,
        pixels: int | Sequence[int] = 40000,
        basis: str = 'xy',
        color_by: str = 'cell',
        colors: dict | None = None,
        seed: int | None = None,
        axes=None,
        legend: bool = False,
        axis_units: str = 'cm',
        outline: bool | str = False,
        show_overlaps: bool = False,
        overlap_color: Sequence[int] | str = (255, 0, 0),
        n_samples: int | None = None,
        plane_tolerance: float = 1.,
        legend_kwargs: dict | None = None,
        source_kwargs: dict | None = None,
        contour_kwargs: dict | None = None,
        **kwargs,
    ):
        """Display a slice plot of the model.

        .. versionadded:: 0.15.1
        """
        import matplotlib.patches as mpatches
        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))
        legend_kwargs.setdefault('loc', 2)
        legend_kwargs.setdefault('borderaxespad', 0.0)
        if source_kwargs is None:
            source_kwargs = {}
        source_kwargs.setdefault('marker', 'x')

        # Set indices using basis and create axis labels
        x, y, z = _BASIS_INDICES[basis]
        xlabel, ylabel = f'{basis[0]} [{axis_units}]', f'{basis[1]} [{axis_units}]'

        # Determine extents of plot
        origin, width, pixels = self._set_plot_defaults(
            origin, width, pixels, basis)

        axis_scaling_factor = {'km': 0.00001, 'm': 0.01, 'cm': 1, 'mm': 10}

        x_min = (origin[x] - 0.5*width[0]) * axis_scaling_factor[axis_units]
        x_max = (origin[x] + 0.5*width[0]) * axis_scaling_factor[axis_units]
        y_min = (origin[y] - 0.5*width[1]) * axis_scaling_factor[axis_units]
        y_max = (origin[y] + 0.5*width[1]) * axis_scaling_factor[axis_units]

        # Determine whether any materials contains macroscopic data and if so,
        # set energy mode accordingly and check that mg cross sections path is accessible
        for mat in self.geometry.get_all_materials().values():
            if mat._macroscopic is not None:
                self.settings.energy_mode = 'multi-group'
                if 'mg_cross_sections' not in openmc.config:
                    raise RuntimeError("'mg_cross_sections' path must be set in "
                                       "openmc.config before plotting.")
                break

        # Get plot IDs from the C API
        id_map, _ = self.slice_data(
            origin=origin,
            width=width,
            pixels=pixels,
            basis=basis,
            show_overlaps=show_overlaps,
            include_properties=False,
        )

        # Generate colors if not provided
        if colors is None and seed is not None:
            # Use the colorize method to generate random colors
            plot = openmc.SlicePlot()
            plot.color_by = color_by
            plot.colorize(self.geometry, seed=seed)
            colors = plot.colors

        # Convert ID map to RGB image
        img = id_map_to_rgb(
            id_map=id_map,
            color_by=color_by,
            colors=colors,
            overlap_color=overlap_color
        )

        # Create a figure sized such that the size of the axes within
        # exactly matches the number of pixels specified
        if axes is None:
            px = 1/plt.rcParams['figure.dpi']
            fig, axes = plt.subplots()
            axes.set_xlabel(xlabel)
            axes.set_ylabel(ylabel)
            params = fig.subplotpars
            width_px = pixels[0]*px/(params.right - params.left)
            height_px = pixels[1]*px/(params.top - params.bottom)
            fig.set_size_inches(width_px, height_px)

        if outline:
            # Combine R, G, B values into a single int for contour detection
            rgb = (img * 256).astype(int)
            image_value = (rgb[..., 0] << 16) + \
                (rgb[..., 1] << 8) + (rgb[..., 2])

            # Set default arguments for contour()
            if contour_kwargs is None:
                contour_kwargs = {}
            contour_kwargs.setdefault('colors', 'k')
            contour_kwargs.setdefault('linestyles', 'solid')
            contour_kwargs.setdefault('algorithm', 'serial')

            axes.contour(
                image_value,
                origin="upper",
                levels=np.unique(image_value),
                extent=(x_min, x_max, y_min, y_max),
                **contour_kwargs
            )

            # If only showing outline, set the axis limits and aspect explicitly
            if outline == 'only':
                axes.set_xlim(x_min, x_max)
                axes.set_ylim(y_min, y_max)
                axes.set_aspect('equal')

        # Add legend showing which colors represent which material or cell
        if legend:
            if colors is None or len(colors) == 0:
                raise ValueError("Must pass 'colors' dictionary if you "
                                 "are adding a legend via legend=True.")

            if color_by == "cell":
                expected_key_type = openmc.Cell
            else:
                expected_key_type = openmc.Material

            patches = []
            for key, color in colors.items():
                if isinstance(key, int):
                    raise TypeError(
                        "Cannot use IDs in colors dict for auto legend.")
                elif not isinstance(key, expected_key_type):
                    raise TypeError(
                        "Color dict key type does not match color_by")

                # this works whether we're doing cells or materials
                label = key.name if key.name != '' else key.id

                # matplotlib takes RGB on 0-1 scale rather than 0-255
                if len(color) == 3 and not isinstance(color, str):
                    scaled_color = (
                        color[0]/255, color[1]/255, color[2]/255)
                else:
                    scaled_color = color

                key_patch = mpatches.Patch(color=scaled_color, label=label)
                patches.append(key_patch)

            axes.legend(handles=patches, **legend_kwargs)

        # Plot image and return the axes
        if outline != 'only':
            axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs)

        if n_samples:
            # Sample external source particles
            particles = self.sample_external_source(n_samples)

            # Get points within tolerance of the slice plane
            slice_value = origin[z]
            xs = []
            ys = []
            tol = plane_tolerance
            for particle in particles:
                if (slice_value - tol < particle.r[z] < slice_value + tol):
                    xs.append(particle.r[x] * axis_scaling_factor[axis_units])
                    ys.append(particle.r[y] * axis_scaling_factor[axis_units])
            axes.scatter(xs, ys, **source_kwargs)

        return axes

    def sample_external_source(
        self,
        n_samples: int = 1000,
        prn_seed: int | None = None,
        as_array: bool = False,
        **init_kwargs
    ) -> openmc.ParticleList | np.ndarray:
        """Sample external source and return source particles.

        .. versionadded:: 0.15.1

        Parameters
        ----------
        n_samples : int
            Number of samples
        prn_seed : int
            Pseudorandom number generator (PRNG) seed; if None, one will be
            generated randomly.
        as_array : bool
            If True, return a numpy structured array instead of a
            :class:`~openmc.ParticleList`.
        **init_kwargs
            Keyword arguments passed to :func:`openmc.lib.init`

        Returns
        -------
        openmc.ParticleList or numpy.ndarray
            List of sampled source particles, or a structured array when
            *as_array* is True.
        """
        import openmc.lib

        # Silence output by default. Also set arguments to start in volume
        # calculation mode to avoid loading cross sections
        init_kwargs.setdefault('output', False)
        init_kwargs.setdefault('args', ['-c'])

        with openmc.lib.TemporarySession(self, **init_kwargs):
            return openmc.lib.sample_external_source(
                n_samples=n_samples, prn_seed=prn_seed, as_array=as_array
            )

    def apply_tally_results(self, statepoint: PathLike | openmc.StatePoint):
        """Apply results from a statepoint to tally objects on the Model

        Parameters
        ----------
        statepoint : PathLike or openmc.StatePoint
            Statepoint file used to update tally results
        """
        self.tallies.add_results(statepoint)

    def plot_geometry(
        self,
        output: bool = True,
        cwd: PathLike = ".",
        openmc_exec: PathLike = "openmc",
        export_model_xml: bool = True,
        **export_kwargs,
    ):
        """Creates plot images as specified by the Model.plots attribute

        .. versionadded:: 0.13.0

        Parameters
        ----------
        output : bool, optional
            Capture OpenMC output from standard out
        cwd : PathLike, optional
            Path to working directory to run in. Defaults to the current
            working directory.
        openmc_exec : PathLike, optional
            Path to OpenMC executable. Defaults to 'openmc'.
            This only applies to the case when not using the C API.
        export_model_xml : bool, optional
            Exports a single model.xml file rather than separate files. Defaults
            to True.
        **export_kwargs
            Keyword arguments passed to either :meth:`Model.export_to_model_xml`
            or :meth:`Model.export_to_xml`.

        """

        if len(self.plots) == 0:
            # Then there is no volume calculation specified
            raise ValueError("The Model.plots attribute must be specified "
                             "before executing this method!")

        with change_directory(cwd):
            if self.is_initialized:
                # Compute the volumes
                openmc.lib.plot_geometry(output)
            else:
                if export_model_xml:
                    self.export_to_model_xml(**export_kwargs)
                else:
                    self.export_to_xml(**export_kwargs)
                path_input = export_kwargs.get("path", None)
                openmc.plot_geometry(output=output, openmc_exec=openmc_exec,
                                     path_input=path_input)

    def _change_py_lib_attribs(
        self,
        names_or_ids: Iterable[str] | Iterable[int],
        value: float | Iterable[float],
        obj_type: str,
        attrib_name: str,
        density_units: str = "atom/b-cm",
    ):
        # Method to do the same work whether it is a cell or material and
        # a temperature or volume
        check_type('names_or_ids', names_or_ids, Iterable, (Integral, str))
        check_type('obj_type', obj_type, str)
        obj_type = obj_type.lower()
        check_value('obj_type', obj_type, ('material', 'cell'))
        check_value('attrib_name', attrib_name,
                    ('temperature', 'volume', 'density', 'rotation',
                     'translation'))
        # The C API only allows setting density units of atom/b-cm and g/cm3
        check_value('density_units', density_units, ('atom/b-cm', 'g/cm3'))
        # The C API has no way to set cell volume or material temperature
        # so lets raise exceptions as needed
        if obj_type == 'cell' and attrib_name == 'volume':
            raise NotImplementedError(
                'Setting a Cell volume is not supported!')
        if obj_type == 'material' and attrib_name == 'temperature':
            raise NotImplementedError(
                'Setting a material temperature is not supported!')

        # And some items just dont make sense
        if obj_type == 'cell' and attrib_name == 'density':
            raise ValueError('Cannot set a Cell density!')
        if obj_type == 'material' and attrib_name in ('rotation',
                                                      'translation'):
            raise ValueError('Cannot set a material rotation/translation!')

        # Set the
        if obj_type == 'cell':
            by_name = self._cells_by_name
            by_id = self._cells_by_id
            if self.is_initialized:
                obj_by_id = openmc.lib.cells
        else:
            by_name = self._materials_by_name
            by_id = self._materials_by_id
            if self.is_initialized:
                obj_by_id = openmc.lib.materials
        # Get the list of ids to use if converting from names and accepting
        # only values that have actual ids
        ids = []
        for name_or_id in names_or_ids:
            if isinstance(name_or_id, Integral):
                if name_or_id in by_id:
                    ids.append(int(name_or_id))
                else:
                    cap_obj = obj_type.capitalize()
                    msg = f'{cap_obj} ID {name_or_id} " \
                        "is not present in the model!'
                    raise InvalidIDError(msg)
            elif isinstance(name_or_id, str):
                if name_or_id in by_name:
                    # Then by_name[name_or_id] is a list so we need to add all
                    # entries
                    ids.extend([obj.id for obj in by_name[name_or_id]])
                else:
                    cap_obj = obj_type.capitalize()
                    msg = f'{cap_obj} {name_or_id} " \
                        "is not present in the model!'
                    raise InvalidIDError(msg)

        # Now perform the change to both python and C API
        for id_ in ids:
            obj = by_id[id_]
            if attrib_name == 'density':
                obj.set_density(density_units, value)
            else:
                setattr(obj, attrib_name, value)
            # Next lets keep what is in C API memory up to date as well
            if self.is_initialized:
                lib_obj = obj_by_id[id_]
                if attrib_name == 'density':
                    lib_obj.set_density(value, density_units)
                elif attrib_name == 'temperature':
                    lib_obj.set_temperature(value)
                else:
                    setattr(lib_obj, attrib_name, value)

    def rotate_cells(
        self, names_or_ids: Iterable[str] | Iterable[int], vector: Iterable[float]
    ):
        """Rotate the identified cell(s) by the specified rotation vector.
        The rotation is only applied to cells filled with a universe.

        .. note:: If applying this change to a name that is not unique, then
                  the change will be applied to all objects of that name.

        .. versionadded:: 0.13.0

        Parameters
        ----------
        names_or_ids : Iterable of str or int
            The cell names (if str) or id (if int) that are to be translated
            or rotated. This parameter can include a mix of names and ids.
        vector : Iterable of float
            The rotation vector of length 3 to apply. This array specifies the
            angles in degrees about the x, y, and z axes, respectively.

        """

        self._change_py_lib_attribs(names_or_ids, vector, 'cell', 'rotation')

    def translate_cells(
        self, names_or_ids: Iterable[str] | Iterable[int], vector: Iterable[float]
    ):
        """Translate the identified cell(s) by the specified translation vector.
        The translation is only applied to cells filled with a universe.

        .. note:: If applying this change to a name that is not unique, then
                  the change will be applied to all objects of that name.

        .. versionadded:: 0.13.0

        Parameters
        ----------
        names_or_ids : Iterable of str or int
            The cell names (if str) or id (if int) that are to be translated
            or rotated. This parameter can include a mix of names and ids.
        vector : Iterable of float
            The translation vector of length 3 to apply. This array specifies
            the x, y, and z dimensions of the translation.

        """

        self._change_py_lib_attribs(names_or_ids, vector, 'cell',
                                    'translation')

    def update_densities(
        self,
        names_or_ids: Iterable[str] | Iterable[int],
        density: float,
        density_units: str = "atom/b-cm",
    ):
        """Update the density of a given set of materials to a new value

        .. note:: If applying this change to a name that is not unique, then
                  the change will be applied to all objects of that name.

        .. versionadded:: 0.13.0

        Parameters
        ----------
        names_or_ids : Iterable of str or int
            The material names (if str) or id (if int) that are to be updated.
            This parameter can include a mix of names and ids.
        density : float
            The density to apply in the units specified by `density_units`
        density_units : {'atom/b-cm', 'g/cm3'}, optional
            Units for `density`. Defaults to 'atom/b-cm'

        """

        self._change_py_lib_attribs(names_or_ids, density, 'material',
                                    'density', density_units)

    def update_cell_temperatures(
        self, names_or_ids: Iterable[str] | Iterable[int], temperature: float
    ):
        """Update the temperature of a set of cells to the given value

        .. note:: If applying this change to a name that is not unique, then
                  the change will be applied to all objects of that name.

        .. versionadded:: 0.13.0

        Parameters
        ----------
        names_or_ids : Iterable of str or int
            The cell names (if str) or id (if int) that are to be updated.
            This parameter can include a mix of names and ids.
        temperature : float
            The temperature to apply in units of Kelvin

        """

        self._change_py_lib_attribs(names_or_ids, temperature, 'cell',
                                    'temperature')

    def update_material_volumes(
        self, names_or_ids: Iterable[str] | Iterable[int], volume: float
    ):
        """Update the volume of a set of materials to the given value

        .. note:: If applying this change to a name that is not unique, then
                  the change will be applied to all objects of that name.

        .. versionadded:: 0.13.0

        Parameters
        ----------
        names_or_ids : Iterable of str or int
            The material names (if str) or id (if int) that are to be updated.
            This parameter can include a mix of names and ids.
        volume : float
            The volume to apply in units of cm^3

        """

        self._change_py_lib_attribs(names_or_ids, volume, 'material', 'volume')

    def differentiate_depletable_mats(self, diff_volume_method: str = None):
        """Assign distribmats for each depletable material

        .. versionadded:: 0.14.0

        .. versionchanged:: 0.15.1
            diff_volume_method default is None, do not set volumes on the new
            material ovjects. Is now a convenience method for
            differentiate_mats(diff_volume_method, depletable_only=True)

        Parameters
        ----------
        diff_volume_method : str
            Specifies how the volumes of the new materials should be found.
            - None: Do not assign volumes to the new materials (Default)
            - 'divide equally': Divide the original material volume equally between the new materials
            - 'match cell': Set the volume of the material to the volume of the cell they fill
        """
        self.differentiate_mats(diff_volume_method, depletable_only=True)

    def differentiate_mats(self, diff_volume_method: str = None, depletable_only: bool = True):
        """Assign distribmats for each material

        .. versionadded:: 0.15.1

        Parameters
        ----------
        diff_volume_method : str
            Specifies how the volumes of the new materials should be found.
            - None: Do not assign volumes to the new materials (Default)
            - 'divide equally': Divide the original material volume equally between the new materials
            - 'match cell': Set the volume of the material to the volume of the cell they fill
        depletable_only : bool
            Default is True, only depletable materials will be differentiated. If False, all materials will be
            differentiated.
        """
        check_value('volume differentiation method', diff_volume_method, ("divide equally", "match cell", None))

        # Count the number of instances for each cell and material
        self.geometry.determine_paths(instances_only=True)

        # Get list of materials
        if self.materials:
            materials = self.materials
        else:
            materials = list(self.geometry.get_all_materials().values())

        # Find all or depletable_only materials which have multiple instance
        distribmats = set()
        for mat in materials:
            # Differentiate all materials with multiple instances
            diff_mat = mat.num_instances > 1
            # If depletable_only is True, differentiate only depletable materials
            if depletable_only:
                diff_mat = diff_mat and mat.depletable
            if diff_mat:
                # Assign volumes to the materials according to requirements
                if diff_volume_method == "divide equally":
                    if mat.volume is None:
                        raise RuntimeError(
                            "Volume not specified for "
                            f"material with ID={mat.id}.")
                    else:
                        mat.volume /= mat.num_instances
                elif diff_volume_method == "match cell":
                    for cell in self.geometry.get_all_material_cells().values():
                        if cell.fill == mat:
                            if not cell.volume:
                                raise ValueError(
                                    f"Volume of cell ID={cell.id} not specified. "
                                    "Set volumes of cells prior to using "
                                    "diff_volume_method='match cell'.")
                distribmats.add(mat)

        if not distribmats:
            return

        # Assign distribmats to cells
        for cell in self.geometry.get_all_material_cells().values():
            if cell.fill in distribmats:
                mat = cell.fill

                # Clone materials
                if cell.num_instances > 1:
                    cell.fill = [mat.clone() for _ in range(cell.num_instances)]
                else:
                    cell.fill = mat.clone()

                # For 'match cell', assign volumes based on the cells
                if diff_volume_method == 'match cell':
                    if cell.fill_type == 'distribmat':
                        for clone_mat in cell.fill:
                            clone_mat.volume = cell.volume
                    else:
                        cell.fill.volume = cell.volume

        if self.materials is not None:
            self.materials = openmc.Materials(
                self.geometry.get_all_materials().values()
            )

    @staticmethod
    def _auto_generate_mgxs_lib(
        model: openmc.model.Model,
        groups: openmc.mgxs.EnergyGroups,
        correction: str | None,
        directory: PathLike,
    ) -> openmc.mgxs.Library:
        """
        Automatically generate a multi-group cross section libray from a model
        with the specified group structure.

        Parameters
        ----------
        model : openmc.Model
            The model to generate the MGXS library from.
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        correction : str
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        directory : str
            Directory to run the simulation in, so as to contain XML files.

        Returns
        -------
        mgxs_lib : openmc.mgxs.Library
            OpenMC MGXS Library object
        """

        # Initialize MGXS library with a finished OpenMC geometry object
        mgxs_lib = openmc.mgxs.Library(model.geometry)

        # Pick energy group structure
        mgxs_lib.energy_groups = groups

        # Disable transport correction
        mgxs_lib.correction = correction

        # Specify needed cross sections for random ray
        if correction == 'P0':
            mgxs_lib.mgxs_types = [
                'nu-transport', 'absorption', 'nu-fission', 'fission',
                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi',
                'kappa-fission'
            ]
        elif correction is None:
            mgxs_lib.mgxs_types = [
                'total', 'absorption', 'nu-fission', 'fission',
                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi',
                'kappa-fission'
            ]

        # Specify a "material" domain type for the cross section tally filters
        mgxs_lib.domain_type = "material"

        # Specify the domains over which to compute multi-group cross sections
        mgxs_lib.domains = model.geometry.get_all_materials().values()

        # Do not compute cross sections on a nuclide-by-nuclide basis
        mgxs_lib.by_nuclide = False

        # Check the library - if no errors are raised, then the library is satisfactory.
        mgxs_lib.check_library_for_openmc_mgxs()

        # Construct all tallies needed for the multi-group cross section library
        mgxs_lib.build_library()

        # Create a "tallies.xml" file for the MGXS Library
        mgxs_lib.add_to_tallies(model.tallies, merge=True)

        # Run
        statepoint_filename = model.run(cwd=directory)

        # Load MGXS
        with openmc.StatePoint(statepoint_filename) as sp:
            mgxs_lib.load_from_statepoint(sp)

        return mgxs_lib

    def _create_mgxs_sources(
        self,
        groups: openmc.mgxs.EnergyGroups,
        spatial_dist: openmc.stats.Spatial,
        source_energy: openmc.stats.Univariate | None = None,
    ) -> list[openmc.IndependentSource]:
        """Create a list of independent sources to use with MGXS generation.

        Note that in all cases, a discrete source that is uniform over all
        energy groups is created (strength = 0.01) to ensure that total cross
        sections are generated for all energy groups. In the case that the user
        has provided a source_energy distribution as an argument, an additional
        source (strength = 0.99) is created using that energy distribution. If
        the user has not provided a source_energy distribution, but the model
        has sources defined, and all of those sources are of IndependentSource
        type, then additional sources are created based on the model's existing
        sources, keeping their energy distributions but replacing their
        spatial/angular distributions, with their combined strength being 0.99.
        If the user has not provided a source_energy distribution and no sources
        are defined on the model and the run mode is 'eigenvalue', then a
        default Watt spectrum source (strength = 0.99) is added.

        Parameters
        ----------
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        spatial_dist : openmc.stats.Spatial
            Spatial distribution to use for all sources.
        source_energy : openmc.stats.Univariate, optional
            Energy distribution to use when generating MGXS data, replacing any
            existing sources in the model.

        Returns
        -------
        list[openmc.IndependentSource]
            A list of independent sources to use for MGXS generation.
        """
        # Make a discrete source that is uniform over the bins of the group structure
        midpoints = []
        strengths = []
        for i in range(groups.num_groups):
            bounds = groups.get_group_bounds(i+1)
            midpoints.append((bounds[0] + bounds[1]) / 2.0)
            strengths.append(1.0)

        uniform_energy = openmc.stats.Discrete(x=midpoints, p=strengths)
        uniform_distribution = openmc.IndependentSource(spatial_dist, energy=uniform_energy, strength=0.01)
        sources = [uniform_distribution]

        # If the user provided an energy distribution, use that
        if source_energy is not None:
            user_energy = openmc.IndependentSource(
                space=spatial_dist, energy=source_energy, strength=0.99)
            sources.append(user_energy)

        # If the user did not provide an energy distribution, create sources
        # based on what is in their model, keeping the energy spectrum but
        # replacing the spatial/angular distributions. We only do this if ALL
        # sources are of IndependentSource type, as we can't pull the energy
        # distribution from e.g. CompiledSource or FileSource types.
        else:
            if self.settings.source is not None:
                for src in self.settings.source:
                    if not isinstance(src, openmc.IndependentSource):
                        break
                else:
                    n_user_sources = len(self.settings.source)
                    for src in self.settings.source:
                        # Create a new IndependentSource with adjusted strength, space, and angle
                        user_source = openmc.IndependentSource(
                            space=spatial_dist,
                            energy=src.energy,
                            strength=0.99 / n_user_sources
                        )
                        sources.append(user_source)
            else:
                # No user sources defined. If we are in eigenvalue mode, then use the default Watt spectrum.
                if self.settings.run_mode == 'eigenvalue':
                    watt_energy = openmc.stats.Watt()
                    watt_source = openmc.IndependentSource(
                        space=spatial_dist, energy=watt_energy, strength=0.99)
                    sources.append(watt_source)

        return sources

    @staticmethod
    def _isothermal_infinite_media_mgxs(
        material: openmc.Material,
        groups: openmc.mgxs.EnergyGroups,
        settings: openmc.Settings,
        correction: str | None,
        directory: PathLike,
        source: openmc.IndependentSource,
        temperature: float | None = None,
    ) -> openmc.XSdata:
        """Generate a single MGXS set for one material, where the geometry is an
        infinite medium composed of that material at an isothermal temperature value.

        Parameters
        ----------
        material : openmc.Material
            The material to generate MGXS for
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        settings : openmc.Settings
            Settings for the generation run, used verbatim except for the
            fields owned by the infinite medium method.
        correction : str
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        directory : str
            Directory to run the simulation in, so as to contain XML files.
        source : openmc.IndependentSource
            Source to use when generating MGXS.
        temperature : float, optional
            The isothermal temperature value to apply to the material. If not specified,
            defaults to the temperature in the material.

        Returns
        -------
        data : openmc.XSdata
            The material MGXS for the given temperature isotherm.
        """
        model = openmc.Model()

        # Set materials on the model
        model.materials = [material]
        if temperature is not None:
          model.materials[-1].temperature = temperature

        # The provided settings are used verbatim, except for the fields
        # owned by the infinite medium method
        model.settings = copy.deepcopy(settings)
        model.settings.source = source

        # Geometry
        box = openmc.model.RectangularPrism(
            100000.0, 100000.0, boundary_type='reflective')
        name = material.name
        infinite_cell = openmc.Cell(name=name, fill=model.materials[-1], region=-box)
        infinite_universe = openmc.Universe(name=name, cells=[infinite_cell])
        model.geometry.root_universe = infinite_universe

        # Generate MGXS
        mgxs_lib = Model._auto_generate_mgxs_lib(
                model, groups, correction, directory)

        if temperature is not None:
            return mgxs_lib.get_xsdata(domain=material, xsdata_name=name,
                                       temperature=temperature)
        else:
            return mgxs_lib.get_xsdata(domain=material, xsdata_name=name)

    def _generate_infinite_medium_mgxs(
        self,
        groups: openmc.mgxs.EnergyGroups,
        settings: openmc.Settings,
        mgxs_path: PathLike,
        correction: str | None,
        directory: PathLike,
        source_energy: openmc.stats.Univariate | None = None,
        temperatures: Sequence[float] | None = None,
    ) -> None:
        """Generate a MGXS library by running multiple OpenMC simulations, each
        representing an infinite medium simulation of a single isolated
        material. A discrete source is used to sample particles, with an equal
        strength spread across each of the energy groups. This is a highly naive
        method that ignores all spatial self shielding effects and all resonance
        shielding effects between materials. If temperature data points are provided,
        isothermal cross sections are generated at each temperature point for
        each material to build a temperature interpolation table.

        Note that in all cases, a discrete source that is uniform over all
        energy groups is created (strength = 0.01) to ensure that total cross
        sections are generated for all energy groups. In the case that the user
        has provided a source_energy distribution as an argument, an additional
        source (strength = 0.99) is created using that energy distribution. If
        the user has not provided a source_energy distribution, but the model
        has sources defined, and all of those sources are of IndependentSource
        type, then additional sources are created based on the model's existing
        sources, keeping their energy distributions but replacing their
        spatial/angular distributions, with their combined strength being 0.99.
        If the user has not provided a source_energy distribution and no sources
        are defined on the model and the run mode is 'eigenvalue', then a
        default Watt spectrum source (strength = 0.99) is added.

        Parameters
        ----------
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        settings : openmc.Settings
            Settings for the generation run(s), used verbatim except for the
            fields owned by the infinite medium method.
        mgxs_path : str
            Filename for the MGXS HDF5 file.
        correction : str
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        directory : str
            Directory to run the simulation in, so as to contain XML files.
        source_energy : openmc.stats.Univariate, optional
            Energy distribution to use when generating MGXS data, replacing any
            existing sources in the model.
        temperatures : Sequence[float], optional
            A list of temperatures to generate MGXS at. Each infinite material region
            is isothermal at a given temperature data point for cross
            section generation.
        """

        src = self._create_mgxs_sources(
            groups,
            spatial_dist=openmc.stats.Point(),
            source_energy=source_energy
        )

        if temperatures is None:
            mgxs_sets = []
            for material in self.materials:
                xs_data = Model._isothermal_infinite_media_mgxs(
                    material,
                    groups,
                    settings,
                    correction,
                    directory,
                    src
                )
                mgxs_sets.append(xs_data)

            # Write the file to disk.
            mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
            for mgxs_set in mgxs_sets:
                mgxs_file.add_xsdata(mgxs_set)
            mgxs_file.export_to_hdf5(mgxs_path)
        else:
            # Build a series of XSData objects, one for each isothermal temperature value.
            raw_mgxs_sets = {}
            for temperature in temperatures:
                raw_mgxs_sets[temperature] = []
                for material in self.materials:
                    xs_data = Model._isothermal_infinite_media_mgxs(
                        material,
                        groups,
                        settings,
                        correction,
                        directory,
                        src,
                        temperature
                    )
                    raw_mgxs_sets[temperature].append(xs_data)

            # Unpack the isothermal XSData objects and build a single XSData object per material.
            mgxs_sets = []
            for m in range(len(self.materials)):
                mgxs_sets.append(openmc.XSdata(self.materials[m].name, groups,
                                               temperatures=temperatures))
                mgxs_sets[-1].order = 0
                for temperature in temperatures:
                    mgxs_sets[-1].add_temperature_data(raw_mgxs_sets[temperature][m])

            # Write the file to disk.
            mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
            for mgxs_set in mgxs_sets:
                mgxs_file.add_xsdata(mgxs_set)
            mgxs_file.export_to_hdf5(mgxs_path)

    @staticmethod
    def _create_stochastic_slab_geometry(
        materials: Sequence[openmc.Material],
        cell_thickness: float = 1.0,
        num_repeats: int = 100,
    ) -> tuple[openmc.Geometry, openmc.stats.Box]:
        """Create a geometry representing a stochastic "sandwich" of materials in a
        layered slab geometry. To reduce the impact of the order of materials in
        the slab, the materials are applied to 'num_repeats' different randomly
        positioned layers of 'cell_thickness' each.

        Parameters
        ----------
        materials : list of openmc.Material
            List of materials to assign. Each material will appear exactly num_repeats times,
            then the ordering is randomly shuffled.
        cell_thickness : float, optional
            Thickness of each lattice cell in x (default 1.0 cm).
        num_repeats : int, optional
            Number of repeats for each material (default 100).

        Returns
        -------
        geometry : openmc.Geometry
            The constructed geometry.
        box : openmc.stats.Box
            A spatial sampling distribution covering the full slab domain.
        """
        if not materials:
            raise ValueError("At least one material must be provided.")

        num_materials = len(materials)
        total_cells = num_materials * num_repeats
        total_width = total_cells * cell_thickness

        # Generate an infinite cell/universe for each material
        universes = []
        for i in range(num_materials):
            cell = openmc.Cell(fill=materials[i])
            universes.append(openmc.Universe(cells=[cell]))

        # Make a list of randomized material idx assignments for the stochastic slab
        assignments = list(range(num_materials)) * num_repeats
        random.seed(42)
        random.shuffle(assignments)

        # Create a list of the (randomized) universe assignments to be used
        # when defining the problem lattice.
        lattice_entries = [universes[m] for m in assignments]

        # Create the RectLattice for the 1D material variation in x.
        lattice = openmc.RectLattice()
        lattice.pitch = (cell_thickness, total_width, total_width)
        lattice.lower_left = (0.0, 0.0, 0.0)
        lattice.universes = [[lattice_entries]]
        lattice.outer = universes[0]

        # Define the six outer surfaces with reflective boundary conditions
        rpp = openmc.model.RectangularParallelepiped(
            0.0, total_width, 0.0, total_width, 0.0, total_width,
            boundary_type='reflective'
        )

        # Create an outer cell that fills with the lattice.
        outer_cell = openmc.Cell(fill=lattice, region=-rpp)

        # Build the geometry
        geometry = openmc.Geometry([outer_cell])

        # Define the spatial distribution that covers the full cubic domain
        box = openmc.stats.Box(*outer_cell.bounding_box)

        return geometry, box

    @staticmethod
    def _isothermal_stochastic_slab_mgxs(
        stoch_geom: openmc.Geometry,
        groups: openmc.mgxs.EnergyGroups,
        settings: openmc.Settings,
        correction: str | None,
        directory: PathLike,
        source: openmc.IndependentSource,
        temperature: float | None = None,
    ) -> dict[str, openmc.XSdata]:
        """Generate MGXS assuming a stochastic "sandwich" of materials in a layered
        slab geometry. If a temperature is specified, all materials in the slab have
        their temperatures set to be isothermal at this temperature.

        Parameters
        ----------
        stoch_geom : openmc.Geometry
            The stochastic slab geometry.
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        settings : openmc.Settings
            Settings for the generation run, used verbatim except for the
            fields owned by the stochastic slab method.
        correction : str
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        directory : str
            Directory to run the simulation in, so as to contain XML files.
        source : openmc.IndependentSource
            Source to use when generating MGXS.
        temperature : float, optional
            The isothermal temperature value to apply to the materials in the
            slab. If not specified, defaults to the temperature in the materials.

        Returns
        -------
        data : dict[str, openmc.XSdata]
            A dictionary where the key is the name of the material and the value is the isothermal MGXS.
        """

        model = openmc.Model()
        model.geometry = stoch_geom

        if temperature is not None:
            for material in model.geometry.get_all_materials().values():
                material.temperature = temperature

        # The provided settings are used verbatim, except for the fields
        # owned by the stochastic slab method
        model.settings = copy.deepcopy(settings)
        model.settings.source = source

        # Generate MGXS
        mgxs_lib = Model._auto_generate_mgxs_lib(
                model, groups, correction, directory)

        # Fetch all of the isothermal results.
        if temperature is not None:
            return {
                mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
                                               temperature=temperature)
                    for mat in mgxs_lib.domains
            }
        else:
            return {
                mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name)
                    for mat in mgxs_lib.domains
            }

    def _generate_stochastic_slab_mgxs(
        self,
        groups: openmc.mgxs.EnergyGroups,
        settings: openmc.Settings,
        mgxs_path: PathLike,
        correction: str | None,
        directory: PathLike,
        source_energy: openmc.stats.Univariate | None = None,
        temperatures: Sequence[float] | None = None,
    ) -> None:
        """Generate MGXS assuming a stochastic "sandwich" of materials in a layered
        slab geometry. While geometry-specific spatial shielding effects are not
        captured, this method can be useful when the geometry has materials only
        found far from the source region that the "material_wise" method would
        not be capable of generating cross sections for. Conversely, this method
        will generate cross sections for all materials in the problem regardless
        of type. If this is a fixed source problem, a discrete source is used to
        sample particles, with an equal strength spread across each of the
        energy groups. If temperature data points are provided,
        isothermal cross sections are generated at each temperature point for
        the stochastic slab to build a temperature interpolation table.

        Parameters
        ----------
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        settings : openmc.Settings
            Settings for the generation run(s), used verbatim except for the
            fields owned by the stochastic slab method.
        mgxs_path : str
            Filename for the MGXS HDF5 file.
        correction : str
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        directory : str
            Directory to run the simulation in, so as to contain XML files.
        source_energy : openmc.stats.Univariate, optional
            Energy distribution to use when generating MGXS data, replacing any
            existing sources in the model. In all cases, a discrete source that
            is uniform over all energy groups is created (strength = 0.01) to
            ensure that total cross sections are generated for all energy
            groups. In the case that the user has provided a source_energy
            distribution as an argument, an additional source (strength = 0.99)
            is created using that energy distribution. If the user has not
            provided a source_energy distribution, but the model has sources
            defined, and all of those sources are of IndependentSource type,
            then additional sources are created based on the model's existing
            sources, keeping their energy distributions but replacing their
            spatial/angular distributions, with their combined strength being
            0.99. If the user has not provided a source_energy distribution and
            no sources are defined on the model and the run mode is
            'eigenvalue', then a default Watt spectrum source (strength = 0.99)
            is added.
        temperatures : Sequence[float], optional
            A list of temperatures to generate MGXS at. Each infinite material region
            is isothermal at a given temperature data point for cross
            section generation.
        """

        # Stochastic slab geometry
        geo, spatial_distribution = Model._create_stochastic_slab_geometry(
            self.materials)

        src = self._create_mgxs_sources(
            groups,
            spatial_dist=spatial_distribution,
            source_energy=source_energy
        )

        if temperatures is None:
            mgxs_sets = Model._isothermal_stochastic_slab_mgxs(
                geo,
                groups,
                settings,
                correction,
                directory,
                src
            ).values()

            # Write the file to disk.
            mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
            for mgxs_set in mgxs_sets:
                mgxs_file.add_xsdata(mgxs_set)
            mgxs_file.export_to_hdf5(mgxs_path)
        else:
            # Build a series of XSData objects, one for each isothermal temperature value.
            raw_mgxs_sets = {}
            for temperature in temperatures:
                raw_mgxs_sets[temperature] = Model._isothermal_stochastic_slab_mgxs(
                    geo,
                    groups,
                    settings,
                    correction,
                    directory,
                    src,
                    temperature
                )

            # Unpack the isothermal XSData objects and build a single XSData object per material.
            mgxs_sets = []
            for mat in self.materials:
                mgxs_sets.append(openmc.XSdata(mat.name, groups, temperatures=temperatures))
                mgxs_sets[-1].order = 0
                for temperature in temperatures:
                    mgxs_sets[-1].add_temperature_data(raw_mgxs_sets[temperature][mat.name])

            # Write the file to disk.
            mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
            for mgxs_set in mgxs_sets:
                mgxs_file.add_xsdata(mgxs_set)
            mgxs_file.export_to_hdf5(mgxs_path)

    @staticmethod
    def _isothermal_materialwise_mgxs(
        input_model: openmc.Model,
        groups: openmc.mgxs.EnergyGroups,
        settings: openmc.Settings,
        correction: str | None,
        directory: PathLike,
        temperature: float | None = None,
    ) -> dict[str, openmc.XSdata]:
        """Generate a material-wise MGXS library for the model by running the
        original continuous energy OpenMC simulation. If a temperature is
        specified, each material in the input model is set to that temperature.
        Otherwise, the original material temperatures are used. If temperature
        data points are provided, isothermal cross sections are generated at
        each temperature point for the whole model to build a temperature
        interpolation table.

        Parameters
        ----------
        input_model : openmc.Model
            The model to use when computing material-wise MGXS.
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        settings : openmc.Settings
            Settings for the generation run, used verbatim.
        correction : str
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        directory : str
            Directory to run the simulation in, so as to contain XML files.
        temperature : float, optional
            The isothermal temperature value to apply to the materials in the
            input model. If not specified, defaults to the temperatures in the
            materials.

        Returns
        -------
        data : dict[str, openmc.XSdata]
            A dictionary where the key is the name of the material and the value is the isothermal MGXS.
        """
        model = copy.deepcopy(input_model)
        model.tallies = openmc.Tallies()

        if temperature is not None:
            for material in model.geometry.get_all_materials().values():
                material.temperature = temperature

        # The provided settings are used verbatim
        model.settings = copy.deepcopy(settings)

        # Generate MGXS
        mgxs_lib = Model._auto_generate_mgxs_lib(
                model, groups, correction, directory)

        # Fetch all of the isothermal results.
        if temperature is not None:
            return {
                mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name,
                                               temperature=temperature)
                    for mat in mgxs_lib.domains
            }
        else:
            return {
                mat.name : mgxs_lib.get_xsdata(domain=mat, xsdata_name=mat.name)
                    for mat in mgxs_lib.domains
            }

    def _generate_material_wise_mgxs(
        self,
        groups: openmc.mgxs.EnergyGroups,
        settings: openmc.Settings,
        mgxs_path: PathLike,
        correction: str | None,
        directory: PathLike,
        temperatures: Sequence[float] | None = None,
    ) -> None:
        """Generate a material-wise MGXS library for the model by running the
        original continuous energy OpenMC simulation of the full material
        geometry and source, and tally MGXS data for each material. This method
        accurately conserves reaction rates totaled over the entire simulation
        domain. However, when the geometry has materials only found far from the
        source region, it is possible the Monte Carlo solver may not be able to
        score any tallies to these material types, thus resulting in zero cross
        section values for these materials. For such cases, the "stochastic
        slab" method may be more appropriate.

        Parameters
        ----------
        groups : openmc.mgxs.EnergyGroups
            Energy group structure for the MGXS.
        settings : openmc.Settings
            Settings for the generation run(s), used verbatim.
        mgxs_path : PathLike
            Filename for the MGXS HDF5 file.
        correction : str
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        directory : PathLike
            Directory to run the simulation in, so as to contain XML files.
        temperatures : Sequence[float], optional
            A list of temperatures to generate MGXS at. Each infinite material region
            is isothermal at a given temperature data point for cross
            section generation.
        """
        if temperatures is None:
            mgxs_sets = Model._isothermal_materialwise_mgxs(
                self, groups, settings, correction, directory).values()

            # Write the file to disk.
            mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
            for mgxs_set in mgxs_sets:
                mgxs_file.add_xsdata(mgxs_set)
            mgxs_file.export_to_hdf5(mgxs_path)
        else:
            # Build a series of XSData objects, one for each isothermal temperature value.
            raw_mgxs_sets = {}
            for temperature in temperatures:
                raw_mgxs_sets[temperature] = Model._isothermal_materialwise_mgxs(
                    self, groups, settings, correction, directory, temperature)

            # Unpack the isothermal XSData objects and build a single XSData object per material.
            mgxs_sets = []
            for mat in self.materials:
                mgxs_sets.append(openmc.XSdata(mat.name, groups, temperatures=temperatures))
                mgxs_sets[-1].order = 0
                for temperature in temperatures:
                    mgxs_sets[-1].add_temperature_data(raw_mgxs_sets[temperature][mat.name])

            # Write the file to disk.
            mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
            for mgxs_set in mgxs_sets:
                mgxs_file.add_xsdata(mgxs_set)
            mgxs_file.export_to_hdf5(mgxs_path)

    def convert_to_multigroup(
        self,
        method: str = "material_wise",
        groups: str | Sequence[float] | openmc.mgxs.EnergyGroups = "CASMO-2",
        nparticles: int | None = None,
        overwrite_mgxs_library: bool = False,
        mgxs_path: PathLike = "mgxs.h5",
        correction: str | None = None,
        source_energy: openmc.stats.Univariate | None = None,
        temperatures: Sequence[float] | None = None,
        temperature_settings: dict | None = None,
        **kwargs,
    ):
        """Convert all materials from continuous energy to multigroup.

        If no MGXS data library file is found, generate one using one or more
        continuous energy Monte Carlo simulations.

        Parameters
        ----------
        method : {"material_wise", "stochastic_slab", "infinite_medium"}, optional
            Method to generate the MGXS.
        groups : openmc.mgxs.EnergyGroups, str, or sequence of float, optional
            Energy group structure for the MGXS. Can be an
            :class:`openmc.mgxs.EnergyGroups` object, a string name of a
            predefined group structure from :data:`openmc.mgxs.GROUP_STRUCTURES`
            (e.g., ``"CASMO-2"``), or a sequence of floats specifying energy
            bin boundaries in eV (e.g., ``[0.0, 1e6]`` for a single group).
            Defaults to ``"CASMO-2"``.
        nparticles : int, optional
            Number of particles to simulate per batch when generating MGXS.
            Defaults to 2000.

            .. deprecated:: 0.15.4
                Pass ``particles`` as a keyword argument instead.
        overwrite_mgxs_library : bool, optional
            Whether to overwrite an existing MGXS library file.
        mgxs_path : str, optional
            Path to the mgxs.h5 library file.
        correction : str, optional
            Transport correction to apply to the MGXS. Options are None and
            "P0".
        source_energy : openmc.stats.Univariate, optional
            Energy distribution to use when generating MGXS data, replacing any
            existing sources in the model. In all cases, a discrete source that
            is uniform over all energy groups is created (strength = 0.01) to
            ensure that total cross sections are generated for all energy
            groups. In the case that the user has provided a source_energy
            distribution as an argument, an additional source (strength = 0.99)
            is created using that energy distribution. If the user has not
            provided a source_energy distribution, but the model has sources
            defined, and all of those sources are of IndependentSource type,
            then additional sources are created based on the model's existing
            sources, keeping their energy distributions but replacing their
            spatial/angular distributions, with their combined strength being
            0.99. If the user has not provided a source_energy distribution and
            no sources are defined on the model and the run mode is
            'eigenvalue', then a default Watt spectrum source (strength = 0.99)
            is added. Note that this argument is only used when using the
            "stochastic_slab" or "infinite_medium" MGXS generation methods.
        temperatures : Sequence[float], optional
            A list of temperatures to generate MGXS at. Each infinite material region
            is isothermal at a given temperature data point for cross
            section generation.
        temperature_settings : dict, optional
            A dictionary of temperature settings to use when generating MGXS.
            Valid entries for temperature_settings are the same as the valid
            entries in openmc.Settings.temperature_settings.

            .. deprecated:: 0.15.4
                Pass ``temperature`` as a keyword argument instead.
        **kwargs
            :class:`openmc.Settings` attributes used to customize the continuous
            energy simulation(s) that generate the MGXS library. Only the
            attributes given override the generation defaults. For example,
            ``model.convert_to_multigroup(particles=100_000)`` adjusts only the
            particle count. The run mode cannot be overridden, as it is
            determined by the generation method. The surrogate-geometry methods
            also construct their own sources and always disable
            ``create_fission_neutrons`` so that fission is treated as capture.
            A ``weight_windows_file`` is applied during ``"material_wise"``
            generation and ignored with a warning by the other methods; see the
            user guide for the weight window "bootstrapping" workflow this
            enables. Cannot be combined with the deprecated ``nparticles`` or
            ``temperature_settings`` arguments.

            .. versionadded:: 0.15.4
        """
        if not isinstance(groups, openmc.mgxs.EnergyGroups):
            groups = openmc.mgxs.EnergyGroups(groups)

        check_value('method', method,
                    ('material_wise', 'stochastic_slab', 'infinite_medium'))

        # Keyword arguments are Settings attributes applied as overrides on
        # the generation defaults
        settings = openmc.Settings(**kwargs) if kwargs else None

        # The model may reference its materials only through the geometry.
        # The materials are converted in place and library-wide attributes
        # like cross_sections must persist on the model afterwards, so
        # populate the explicit collection if it is empty (sorted by ID for
        # reproducibility, since geometry traversal order is arbitrary)
        if not self.materials:
            self.materials = openmc.Materials(sorted(
                self.geometry.get_all_materials().values(),
                key=lambda mat: mat.id))

        if nparticles is not None or temperature_settings is not None:
            warnings.warn(
                'The "nparticles" and "temperature_settings" arguments are '
                'deprecated. Pass "particles" and "temperature" as keyword '
                'arguments instead.', FutureWarning)
            if settings is not None:
                raise ValueError(
                    'The deprecated "nparticles" and "temperature_settings" '
                    'arguments cannot be combined with Settings keyword '
                    'arguments.')

        # Resolve the settings for the MGXS generation run(s) in three layers,
        # with later layers taking precedence: the model's own settings
        # ("material_wise") or a fresh Settings object (surrogate methods), then
        # the generation defaults, then the user's keyword-argument overrides.
        user_settings = settings
        if method == 'material_wise':
            settings = copy.deepcopy(self.settings)
        else:
            settings = openmc.Settings()
            settings.temperature = copy.deepcopy(self.settings.temperature)

        settings.batches = 100 if method == 'infinite_medium' else 200
        if method != 'infinite_medium':
            settings.inactive = 100
        settings.particles = 2000
        settings.output = {'summary': True, 'tallies': False}
        if nparticles is not None:
            settings.particles = nparticles
        if temperature_settings is not None:
            settings.temperature = temperature_settings

        if user_settings is not None:
            # The surrogate-geometry methods construct their own sources
            if method != "material_wise" and len(user_settings.source) > 0:
                warnings.warn(
                    'The given "source" setting is ignored by the '
                    f'"{method}" MGXS generation method, which constructs '
                    'its own sources.')
            # Apply the settings attributes passed by the caller
            for name in kwargs:
                if hasattr(type(settings), name):
                    setattr(settings, name, copy.deepcopy(
                        getattr(user_settings, name)))

        # The run mode is the one attribute that cannot be detected as
        # user-populated (a fresh Settings object defaults it to 'eigenvalue'),
        # so it is owned by the generation method: "material_wise" always takes
        # it from the model, while the surrogate-geometry methods always run in
        # fixed source mode. The surrogate-geometry methods also treat fission
        # as capture (nu-fission is still tallied)
        if method == 'material_wise':
            settings.run_mode = self.settings.run_mode
        else:
            settings.run_mode = 'fixed source'
            settings.create_fission_neutrons = False

        # A weight windows file on the generation settings is loaded and
        # applied (specifying a file turns weight windows on) during the
        # "material_wise" method's continuous energy simulation of the
        # original geometry, allowing materials far from the source --
        # which an analog simulation may struggle to reach -- to still be
        # tallied, and thus obtain nonzero cross sections. The
        # "stochastic_slab" and "infinite_medium" methods use simplified
        # surrogate geometries for which weight windows defined over the
        # original geometry are neither applicable nor needed.
        if settings.weight_windows_file is not None and \
                method != "material_wise":
            warnings.warn(
                'The "weight_windows_file" setting is only applicable to '
                'the "material_wise" MGXS generation method and will be '
                f'ignored for the "{method}" method.'
            )
            settings.weight_windows_file = None

        # Do all work (including MGXS generation) in a temporary directory
        # to avoid polluting the working directory with residual XML files
        with TemporaryDirectory() as tmpdir:

            # Determine if there are DAGMC universes in the model. If so, we need to synchronize
            # the dagmc materials with cells.
            # TODO: Can this be done without having to init/finalize?
            for univ in self.geometry.get_all_universes().values():
                if isinstance(univ, openmc.DAGMCUniverse):
                    # Initialize in stochastic volume mode (non-transport mode)
                    # This mode doesn't require
                    # valid transport settings like particles/batches
                    original_run_mode = self.settings.run_mode
                    self.settings.run_mode = 'volume'
                    self.init_lib(directory=tmpdir)
                    self.sync_dagmc_universes()
                    self.finalize_lib()
                    # Restore original run mode
                    self.settings.run_mode = original_run_mode
                    break

            # Temporarily replace each material's name with a unique, valid HDF5
            # dataset name (its name plus ID) for use as its MGXS library entry
            # and macroscopic. The ID keeps the name unique even when materials
            # share a name; the original names are restored at the end.
            original_names = [material.name for material in self.materials]
            for material in self.materials:
                base = material.name if material.name and material.name.strip() \
                    else "material"
                material.name = re.sub(r'[^a-zA-Z0-9]', '_', base) + f"_{material.id}"

            # If needed, generate the needed MGXS data library file
            if not Path(mgxs_path).is_file() or overwrite_mgxs_library:
                if method == "infinite_medium":
                    self._generate_infinite_medium_mgxs(
                        groups, settings, mgxs_path, correction, tmpdir,
                        source_energy, temperatures)
                elif method == "material_wise":
                    self._generate_material_wise_mgxs(
                        groups, settings, mgxs_path, correction, tmpdir,
                        temperatures)
                elif method == "stochastic_slab":
                    self._generate_stochastic_slab_mgxs(
                        groups, settings, mgxs_path, correction, tmpdir,
                        source_energy, temperatures)
                else:
                    raise ValueError(
                        f'MGXS generation method "{method}" not recognized')
            else:
                print(f'Existing MGXS library file "{mgxs_path}" will be used')

            # Convert all continuous energy materials to multigroup
            self.materials.cross_sections = mgxs_path
            for material in self.materials:
                material.set_density('macro', 1.0)
                material._nuclides = []
                material._sab = []
                material.add_macroscopic(material.name)

            self.settings.energy_mode = 'multi-group'

            # Restore the user's original material names.
            for material, name in zip(self.materials, original_names):
                material.name = name

    def convert_to_random_ray(self):
        """Convert a multigroup model to use random ray.

        This method determines values for the needed settings and adds them to
        the settings.random_ray dictionary so as to enable random ray mode. The
        settings that are populated are:

        - 'ray_source' (openmc.IndependentSource): Where random ray starting
          points are sampled from.
        - 'distance_inactive' (float): The "dead zone" distance at the beginning
          of the ray.
        - 'distance_active' (float): The "active" distance of the ray
        - 'particles' (int): Number of rays to simulate

        The method will determine reasonable defaults for each of the above
        variables based on analysis of the model's geometry. The function will
        have no effect if the random ray dictionary is already defined in the
        model settings.
        """
        # If the random ray dictionary is already set, don't overwrite it
        if self.settings.random_ray:
            warnings.warn("Random ray conversion skipped as "
                          "settings.random_ray dictionary is already set.")
            return

        if self.settings.energy_mode != 'multi-group':
            raise ValueError(
                "Random ray conversion failed: energy mode must be "
                "'multi-group'. Use convert_to_multigroup() first."
            )

        # Helper function for detecting infinity
        def _replace_infinity(value):
            if np.isinf(value):
                return 1.0 if value > 0 else -1.0
            return value

        # Get a bounding box for sampling rays. We can utilize the geometry's bounding box
        # though for 2D problems we need to detect the infinities and replace them with an
        # arbitrary finite value.
        bounding_box = self.geometry.bounding_box
        lower_left = [_replace_infinity(v) for v in bounding_box.lower_left]
        upper_right = [_replace_infinity(v) for v in bounding_box.upper_right]
        uniform_dist_ray = openmc.stats.Box(lower_left, upper_right)
        rr_source = openmc.IndependentSource(space=uniform_dist_ray)
        self.settings.random_ray['ray_source'] = rr_source

        # For the dead zone and active length, a reasonable guess is the larger of either:
        # 1) The maximum chord length through the geometry (as defined by its bounding box)
        # 2) 30 cm
        # Then, set the active length to be 5x longer than the dead zone length, for the sake of efficiency.
        chord_length = np.array(upper_right) - np.array(lower_left)
        max_length = max(np.linalg.norm(chord_length), 30.0)

        self.settings.random_ray['distance_inactive'] = max_length
        self.settings.random_ray['distance_active'] = 5 * max_length

        # Take a wild guess as to how many rays are needed
        self.settings.particles = 2 * int(max_length)

    def keff_search(
        self,
        func: ModelModifier,
        x0: float,
        x1: float,
        target: float = 1.0,
        k_tol: float = 1e-4,
        sigma_final: float = 3e-4,
        p: float = 0.5,
        q: float = 0.95,
        memory: int = 4,
        x_min: float | None = None,
        x_max: float | None = None,
        b0: int | None = None,
        b_min: int = 20,
        b_max: int | None = None,
        maxiter: int = 50,
        output: bool = False,
        func_kwargs: dict[str, Any] | None = None,
        run_kwargs: dict[str, Any] | None = None,
    ) -> SearchResult:
        r"""Perform a keff search on a model parametrized by a single variable.

        This method uses the GRsecant method described in a paper by `Price and
        Roskoff <https://doi.org/10.1016/j.pnucene.2023.104731>`_. The GRsecant
        method is a modification of the secant method that accounts for
        uncertainties in the function evaluations. The method uses a weighted
        linear fit of the most recent function evaluations to predict the next
        point to evaluate. It also adaptively changes the number of batches to
        meet the target uncertainty value at each iteration.

        The target uncertainty for iteration :math:`n+1` is determined by the
        following equation (following Eq. (8) in the paper):

        .. math::
            \sigma_{i+1} = q \sigma_\text{final} \left ( \frac{ \min \left \{
            \left\lvert k_i - k_\text{target} \right\rvert : k=0,1,\dots,n
            \right \} }{k_\text{tol}} \right )^p

        where :math:`q` is a multiplicative factor less than 1, given as the
        ``sigma_factor`` parameter below.

        Parameters
        ----------
        func : ModelModifier
            Function that takes the parameter to be searched and makes a
            modification to the model.
        x0 : float
            First guess for the parameter passed to `func`
        x1 : float
            Second guess for the parameter passed to `func`
        target : float, optional
            keff value to search for
        k_tol : float, optional
            Stopping criterion on the function value; the absolute value must be
            within ``k_tol`` of zero to be accepted.
        sigma_final : float, optional
            Maximum accepted k-effective uncertainty for the stopping criterion.
        p : float, optional
            Exponent used in the stopping criterion.
        q : float, optional
            Multiplicative factor used in the stopping criterion.
        memory : int, optional
            Number of most-recent points used in the weighted linear fit of
            ``f(x) = a + b x`` to predict the next point.
        x_min : float, optional
            Minimum allowed value for the parameter ``x``.
        x_max : float, optional
            Maximum allowed value for the parameter ``x``.
        b0 : int, optional
            Number of active batches to use for the initial function
            evaluations. If None, uses the model's current setting.
        b_min : int, optional
            Minimum number of active batches to use in a function evaluation.
        b_max : int, optional
            Maximum number of active batches to use in a function evaluation.
        maxiter : int, optional
            Maximum number of iterations to perform.
        output : bool, optional
            Whether or not to display output showing iteration progress.
        func_kwargs : dict, optional
            Keyword-based arguments to pass to the `func` function.
        run_kwargs : dict, optional
            Keyword arguments to pass to :meth:`openmc.Model.run` or
            :meth:`openmc.lib.run`.

        Returns
        -------
        SearchResult
            Result object containing the estimated root (parameter value) and
            evaluation history (parameters, means, standard deviations, and
            batches), plus convergence status and termination reason.

        """
        import openmc.lib

        check_type('model modifier', func, Callable)
        check_type('target', target, Real)
        if memory < 2:
            raise ValueError("memory must be ≥ 2")
        func_kwargs = {} if func_kwargs is None else dict(func_kwargs)
        run_kwargs = {} if run_kwargs is None else dict(run_kwargs)
        run_kwargs.setdefault('output', False)

        # Create lists to store the history of evaluations
        xs: list[float] = []
        fs: list[float] = []
        ss: list[float] = []
        gs: list[int] = []
        count = 0

        # Helper function to evaluate f and store results
        def eval_at(x: float, batches: int) -> tuple[float, float]:
            # Modify the model with the current guess
            func(x, **func_kwargs)

            # Change the number of batches and run the model
            batches += self.settings.inactive
            if openmc.lib.is_initialized:
                openmc.lib.settings.set_batches(batches)
                openmc.lib.reset()
                openmc.lib.run(**run_kwargs)
                sp_filepath = f'statepoint.{batches}.h5'
            else:
                self.settings.batches = batches
                sp_filepath = self.run(**run_kwargs)

            # Extract keff and its uncertainty
            with openmc.StatePoint(sp_filepath) as sp:
                keff = sp.keff

            if output:
                nonlocal count
                count += 1
                print(f'Iteration {count}: {batches=}, {x=:.6g}, {keff=:.5f}')

            xs.append(float(x))
            fs.append(float(keff.n - target))
            ss.append(float(keff.s))
            gs.append(int(batches))
            return fs[-1], ss[-1]

        # Default b0 to current model settings if not explicitly provided
        if b0 is None:
            b0 = self.settings.batches - self.settings.inactive

        # Perform the search (inlined GRsecant) in a temporary directory
        with TemporaryDirectory() as tmpdir:
            if not openmc.lib.is_initialized:
                run_kwargs.setdefault('cwd', tmpdir)

            # ---- Seed with two evaluations
            f0, s0 = eval_at(x0, b0)
            if abs(f0) <= k_tol and s0 <= sigma_final:
                return SearchResult(x0, xs, fs, ss, gs, True, "converged")
            f1, s1 = eval_at(x1, b0)
            if abs(f1) <= k_tol and s1 <= sigma_final:
                return SearchResult(x1, xs, fs, ss, gs, True, "converged")

            for _ in range(maxiter - 2):
                # ------ Step 1: propose next x via GRsecant
                m = min(memory, len(xs))

                # Perform a curve fit on f(x) = a + bx accounting for
                # uncertainties. This is equivalent to minimizing the function
                # in Equation (A.14)
                (a, b), _ = curve_fit(
                    lambda x, a, b: a + b*x,
                    xs[-m:], fs[-m:], sigma=ss[-m:], absolute_sigma=True
                )
                x_new = float(-a / b)

                # Clamp x_new to the bounds if provided
                if x_min is not None:
                    x_new = max(x_new, x_min)
                if x_max is not None:
                    x_new = min(x_new, x_max)

                # ------ Step 2: choose target σ for next run (Eq. 8 + clamp)

                min_abs_f = float(np.min(np.abs(fs)))
                base = q * sigma_final
                ratio = min_abs_f / k_tol if k_tol > 0 else 1.0
                sig = base * (ratio ** p)
                sig_target = max(sig, base)

                # ------ Step 3: choose generations to hit σ_target (Appendix C)

                # Use at least two past points for regression
                if len(gs) >= 2 and np.var(np.log(gs)) > 0.0:
                    # Perform a curve fit based on Eq. (C.3) to solve for ln(k).
                    # Note that unlike in the paper, we do not leave r as an
                    # undetermined parameter and choose r=0.5.
                    (ln_k,), _ = curve_fit(
                        lambda ln_b, ln_k: ln_k - 0.5*ln_b,
                        np.log(gs[-4:]), np.log(ss[-4:]),
                    )
                    k = float(np.exp(ln_k))
                else:
                    k = float(ss[-1] * math.sqrt(gs[-1]))

                b_new = (k / sig_target) ** 2

                # Clamp and round up to integer
                b_new = max(b_min, math.ceil(b_new))
                if b_max is not None:
                    b_new = min(b_new, b_max)

                # Evaluate at proposed x with batches determined above
                f_new, s_new = eval_at(x_new, b_new)

                # Termination based on both criteria (|f| and σ)
                if abs(f_new) <= k_tol and s_new <= sigma_final:
                    return SearchResult(x_new, xs, fs, ss, gs, True, "converged")

            return SearchResult(xs[-1], xs, fs, ss, gs, False, "maxiter")


@dataclass
class SearchResult:
    """Result of a GRsecant keff search.

    Attributes
    ----------
    root : float
        Estimated parameter value where f(x) = 0 at termination.
    parameters : list[float]
        Parameter values (x) evaluated during the search, in order.
    keffs : list[float]
        Estimated keff values for each evaluation.
    stdevs : list[float]
        One-sigma uncertainties of keff for each evaluation.
    batches : list[int]
        Number of active batches used for each evaluation.
    converged : bool
        Whether both |f| <= k_tol and sigma <= sigma_final were met.
    flag : str
        Reason for termination (e.g., "converged", "maxiter").
    """
    root: float
    parameters: list[float] = field(repr=False)
    means: list[float] = field(repr=False)
    stdevs: list[float] = field(repr=False)
    batches: list[int] = field(repr=False)
    converged: bool
    flag: str

    @property
    def function_calls(self) -> int:
        """Number of function evaluations performed."""
        return len(self.parameters)

    @property
    def total_batches(self) -> int:
        """Total number of active batches used across all evaluations."""
        return sum(self.batches)