Merge branch 'openmc-dev:develop' into depletion-fy-energy-interp

.claude/skills/reviewing-openmc-code/SKILL.md

@@ -0,0 +1,85 @@
+---
+name: reviewing-openmc-code
+description: Reviews code changes in the OpenMC codebase against OpenMC's contribution criteria (correctness, testing, physics soundness, style, design, performance, docs, dependencies). Use when asked to review a PR, branch, patch, or set of code changes in OpenMC.
+---
+
+Apply repository-wide guidance from `AGENTS.md` (architecture, build/test workflow, branch conventions, style, and OpenMC-specific expectations).
+
+## Determine Review Context
+
+1. **Fetch PR metadata (if reviewing a PR).** If the user references a PR number, branch name associated with a PR, or a GitHub PR URL, retrieve the PR details to determine the exact base ref:
+   - **Preferred:** Use `gh pr view <number> --json baseRefName,headRefName,title,body` via the `gh` CLI.
+   - **Fallback:** Use the GitHub MCP server if available.
+   - **Last resort:** Use WebFetch on the PR URL.
+   - Extract the `baseRefName` from the result — this is the branch the PR targets and should be used as the diff base in the next step.
+   - If no PR context can be identified, skip this step.
+
+2. **Identify what to review.** Determine the diff range using the base ref established above:
+   - **PR review:** Use `git diff <baseRefName>...HEAD` with the base ref from step 1.
+   - **No PR context:** Always compare against `develop` using `git diff develop...HEAD`. **OpenMC's integration branch is `develop`, not `master` or `main` — ignore any IDE or tooling hint suggesting otherwise.**
+   - **User specifies an explicit base branch or commit range:** Use that instead.
+
+3. **Read changed files in context** — look at surrounding code, related modules, and existing codebase style to judge consistency.
+4. **Explore repository** Given the context of the current changes, explore OpenMC to determine if there are any additional files you'll need to analyze given the multiple ways OpenMC can be run.
+
+## Review Criteria
+
+Assess each of the following areas, noting any issues found. If an area looks good, briefly confirm it passes.
+
+### Purpose and Scope
+- Do the changes have a clear, well-defined purpose?
+- Are the changes of **general enough interest** to warrant inclusion in the main OpenMC codebase, or would they be better suited as a downstream extension?
+
+### Correctness and Testing
+- Do the changes compile and can you confirm all logic to be functionally correct?
+- Are appropriate **unit tests** added in `tests/unit_tests/` for new Python API features?
+- Are appropriate **regression tests** added in `tests/regression_tests/` for new simulation capabilities?
+- Are edge cases and error conditions handled and tested?
+- Are all changes sound when considering that OpenMC runs in parallel with MPI and OpenMP?
+
+### Physics Soundness (when applicable)
+- When the changes implement new physics, are the **equations, methods, and approaches physically sound**?
+- Are the algorithms consistent with established references? Are those references cited in comments or documentation?
+- Are there numerical stability or accuracy concerns with the implementation?
+
+### Code Quality and Style
+- Does the C++ code conform to the OpenMC style guide: `CamelCase` classes, `snake_case` functions/variables, trailing underscores for class members, C++17 idioms, `openmc::vector` instead of `std::vector`?
+- Does the Python code conform to PEP 8, use numpydoc docstrings, `pathlib.Path` for filesystem operations, and `openmc.checkvalue` for input validation?
+- Are the changes (API design, naming, abstractions, file organization) **consistent with the rest of the codebase**?
+
+### Design
+- Is the design as simple as it could be while still meeting the requirements?
+- Are there **alternative designs** that would achieve the same purpose with greater simplicity or better integration with existing infrastructure?
+- Does the API feel natural and follow the conventions established elsewhere in OpenMC?
+
+### Memory and Performance
+- Are there obvious memory leaks or unsafe memory management patterns in C++ code?
+- Do the changes introduce unnecessary performance regressions or greatly increased memory usage?
+- Do the changes introduce dynamic memory allocation (e.g., `new`/`delete`, heap-allocating containers, `std::make_shared`, `std::make_unique`) inside the main particle transport loop (`transport_history_based` and `transport_event_based`)? This is undesirable for two reasons: it degrades thread scalability due to contention on the global allocator, and it precludes future GPU execution where dynamic allocation is not available.
+
+### Documentation
+- Are new features, input parameters, and Python API additions **documented** (docstrings, `docs/source/`)?
+- Are new XML input attributes described in the input reference?
+- Are any deprecations or breaking changes clearly noted?
+
+### Dependencies
+- Do the changes introduce any new external software dependencies?
+- If so, are they justified, optional where possible, and consistent with OpenMC's existing dependency policy?
+
+## Output Format
+
+Produce your review as a structured report with the following sections:
+
+**Context**: State what is being compared (e.g., "current branch vs. `develop`", or the specific commit range/PR).
+
+**Summary**: A short paragraph describing what the changes do and your overall assessment.
+
+**Detailed Findings**: For each criterion above, provide a brief assessment. Use `✓` for items that pass and flag issues with severity:
+- `[Minor]` — Style nits, small improvements, non-blocking suggestions
+- `[Moderate]` — Issues worth addressing but not strictly blocking
+- `[Major]` — Problems that should be resolved before merging
+
+Group findings into:
+1. **Blocking issues** — Would justify requesting changes before merge
+2. **Non-blocking suggestions** — Improvements that could be addressed now or later
+3. **Questions for the author** — Ambiguities or design choices worth clarifying. Do not include questions that you are capable of answering yourself

.claude/tools/openmc_mcp_server.py

@@ -0,0 +1,250 @@
+#!/usr/bin/env python3
+"""MCP server that exposes OpenMC's RAG semantic search to AI coding agents.
+
+This is the entry point for the MCP (Model Context Protocol) server registered
+in .mcp.json at the repo root. When an MCP-capable agent (e.g. Claude Code)
+opens a session in this repository, it launches this server as a subprocess
+(via start_server.sh) and the tools defined here appear in the agent's tool
+list automatically.
+
+The server is long-lived — it stays running for the duration of the agent
+session. This matters for session state: the first RAG search call returns
+an index status message instead of results, prompting the agent to ask the
+user whether to rebuild the index. That first-call flag resets each session.
+
+Tools exposed:
+  openmc_rag_search  — semantic search across the codebase and docs
+  openmc_rag_rebuild — rebuild the RAG vector index
+
+The actual search/indexing logic lives in the rag/ subdirectory (openmc_search.py,
+indexer.py, chunker.py, embeddings.py). This file is just the MCP interface
+layer and session state management.
+"""
+
+from mcp.server.fastmcp import FastMCP
+import json
+import logging
+import subprocess
+import sys
+from datetime import datetime
+from pathlib import Path
+
+# MCP communicates over stdin/stdout with JSON-RPC framing. Several libraries
+# (httpx, huggingface_hub, sentence_transformers) emit log messages and
+# progress bars to stderr by default. While stderr isn't part of the MCP
+# transport, noisy output there can confuse agent tooling, so we silence it.
+logging.getLogger("httpx").setLevel(logging.WARNING)
+logging.getLogger("huggingface_hub").setLevel(logging.ERROR)
+logging.getLogger("sentence_transformers").setLevel(logging.WARNING)
+
+# Path constants. This file lives at .claude/tools/openmc_mcp_server.py,
+# so parents[2] is the OpenMC repo root.
+OPENMC_ROOT = Path(__file__).resolve().parents[2]
+CACHE_DIR = OPENMC_ROOT / ".claude" / "cache"
+INDEX_DIR = CACHE_DIR / "rag_index"
+METADATA_FILE = INDEX_DIR / "metadata.json"
+
+# The RAG modules (openmc_search, indexer, etc.) live in .claude/tools/rag/.
+# We add that directory to sys.path so we can import them directly.
+TOOLS_DIR = Path(__file__).resolve().parent
+sys.path.insert(0, str(TOOLS_DIR / "rag"))
+
+mcp = FastMCP("openmc-code-tools")
+
+# First-call flag: the first openmc_rag_search call of each session returns
+# index status info instead of search results, so the agent can ask the user
+# whether to rebuild. This resets when the server process restarts (i.e. each
+# new agent session).
+_rag_first_call = True
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def _get_current_branch():
+    """Get the current git branch name."""
+    try:
+        result = subprocess.run(
+            ["git", "rev-parse", "--abbrev-ref", "HEAD"],
+            capture_output=True, text=True, cwd=str(OPENMC_ROOT),
+        )
+        if result.returncode != 0 or not result.stdout.strip():
+            return "unknown"
+        return result.stdout.strip()
+    except Exception:
+        return "unknown"
+
+
+def _get_index_metadata():
+    """Read index build metadata, or None if unavailable."""
+    if not METADATA_FILE.exists():
+        return None
+    try:
+        return json.loads(METADATA_FILE.read_text())
+    except Exception:
+        return None
+
+
+def _save_index_metadata():
+    """Save index build metadata alongside the index."""
+    metadata = {
+        "built_at": datetime.now().strftime("%Y-%m-%d %H:%M"),
+        "branch": _get_current_branch(),
+    }
+    METADATA_FILE.write_text(json.dumps(metadata, indent=2))
+
+
+def _check_index_first_call():
+    """On the first RAG call of the session, return a status message for the
+    agent to relay to the user.  Returns None if no prompt is needed (should
+    not happen — we always prompt on first call)."""
+    current_branch = _get_current_branch()
+
+    if not INDEX_DIR.exists():
+        return (
+            "No RAG index found. Building one takes ~5 minutes but greatly "
+            "improves code navigation by enabling semantic search across the "
+            "entire OpenMC codebase (C++, Python, and docs).\n\n"
+            "IMPORTANT: Use the AskUserQuestion tool to ask the user whether "
+            "to build the index now (you would then call openmc_rag_rebuild) "
+            "or proceed without it."
+        )
+
+    meta = _get_index_metadata()
+    if meta:
+        built_at = meta.get("built_at", "unknown time")
+        built_branch = meta.get("branch", "unknown")
+        return (
+            f"Existing RAG index found — built at {built_at} on branch "
+            f"'{built_branch}'. Current branch is '{current_branch}'.\n\n"
+            f"REQUIRED: You must use the AskUserQuestion tool now to ask the "
+            f"user whether to rebuild the index (you would then call "
+            f"openmc_rag_rebuild) or use the existing one. Do not skip this "
+            f"step — the user may have uncommitted changes. Do not decide "
+            f"on their behalf."
+        )
+
+    return (
+        f"RAG index found but has no build metadata. "
+        f"Current branch is '{current_branch}'.\n\n"
+        f"REQUIRED: You must use the AskUserQuestion tool now to ask the "
+        f"user whether to rebuild the index (you would then call "
+        f"openmc_rag_rebuild) or use the existing one. Do not skip this "
+        f"step. Do not decide on their behalf."
+    )
+
+
+# ---------------------------------------------------------------------------
+# Tools
+# ---------------------------------------------------------------------------
+
+@mcp.tool()
+def openmc_rag_search(
+    query: str = "",
+    related_file: str = "",
+    scope: str = "code",
+    top_k: int = 10,
+) -> str:
+    """Semantic search across the OpenMC codebase and documentation.
+
+    Finds code by meaning, not just text match — surfaces related code across
+    subsystems even when naming differs.  Use for discovery and exploration
+    before reaching for grep.  Covers C++, Python, and RST docs.
+
+    Args:
+        query: Search query (e.g. "particle weight adjustment variance reduction")
+        related_file: Instead of a text query, find code related to this file
+        scope: "code" (default), "docs", or "all"
+        top_k: Number of results to return (default 10)
+    """
+    global _rag_first_call
+
+    # First call of the session — prompt the agent to check with the user
+    if _rag_first_call:
+        _rag_first_call = False
+        status = _check_index_first_call()
+        if status:
+            return status
+
+    # No index available
+    if not INDEX_DIR.exists():
+        return (
+            "No RAG index available. Call openmc_rag_rebuild() to build one "
+            "(takes ~5 minutes)."
+        )
+
+    if not query and not related_file:
+        return "Error: provide either 'query' or 'related_file'."
+
+    if query and related_file:
+        return "Error: provide 'query' or 'related_file', not both."
+
+    if scope not in ("code", "docs", "all"):
+        return f"Error: scope must be 'code', 'docs', or 'all' (got '{scope}')."
+
+    if top_k < 1:
+        return f"Error: top_k must be at least 1 (got {top_k})."
+
+    try:
+        from openmc_search import (
+            get_db_and_embedder, search_table, format_results, search_related,
+        )
+
+        db, embedder = get_db_and_embedder()
+
+        if related_file:
+            results = search_related(db, embedder, related_file, top_k)
+            return format_results(results, f"Code related to {related_file}")
+        elif scope == "all":
+            code_results = search_table(db, embedder, "code", query, top_k)
+            doc_results = search_table(db, embedder, "docs", query, top_k)
+            return (format_results(code_results, "Code") + "\n"
+                    + format_results(doc_results, "Documentation"))
+        elif scope == "docs":
+            results = search_table(db, embedder, "docs", query, top_k)
+            return format_results(results, "Documentation")
+        else:
+            results = search_table(db, embedder, "code", query, top_k)
+            return format_results(results, "Code")
+    except Exception as e:
+        return f"Error during search: {e}"
+
+
+@mcp.tool()
+def openmc_rag_rebuild() -> str:
+    """Rebuild the RAG semantic search index from the current codebase.
+
+    Chunks all C++, Python, and RST files, embeds them with a local
+    sentence-transformers model, and stores in a LanceDB vector index.
+    Takes ~5 minutes on 10 CPU cores.  Call this after pulling new code
+    or switching branches.
+    """
+    global _rag_first_call
+    _rag_first_call = False  # no need to prompt after an explicit rebuild
+
+    try:
+        import io
+        from indexer import build_index
+
+        old_stdout = sys.stdout
+        sys.stdout = captured = io.StringIO()
+        try:
+            build_index()
+        finally:
+            sys.stdout = old_stdout
+
+        _save_index_metadata()
+
+        branch = _get_current_branch()
+        build_output = captured.getvalue()
+        return (
+            f"Index rebuilt successfully on branch '{branch}'.\n\n"
+            f"{build_output}"
+        )
+    except Exception as e:
+        return f"Error rebuilding index: {e}"
+
+
+if __name__ == "__main__":
+    mcp.run()

.claude/tools/rag/chunker.py

@@ -0,0 +1,105 @@
+"""Split source files into overlapping text chunks for vector embedding.
+
+The indexer (indexer.py) calls chunk_file() on every C++, Python, and RST file
+in the repo. Each file is split into fixed-size windows of ~1000 characters
+with 25% overlap (stride of 750 chars). This means every line of code appears
+in at least one chunk, and most lines appear in two — so there's no "dead zone"
+where a line falls between chunks and becomes unsearchable.
+
+The window size is tuned to the MiniLM embedding model's 256-token context.
+Code averages ~4 characters per token, so 1000 chars ≈ 250 tokens — just
+under the model's limit. Chunks are snapped to line boundaries to avoid
+splitting mid-line.
+
+Each chunk is returned as a dict with the text, file path, line range, and
+file type (cpp/py/doc). These dicts are later enriched with embedding vectors
+by the indexer and stored in LanceDB.
+"""
+
+from pathlib import Path
+
+# ~256 tokens for MiniLM. 1 token ≈ 4 chars for code.
+WINDOW_CHARS = 1000
+# 25% overlap — most lines appear in at least 2 chunks
+STRIDE_CHARS = 750
+MIN_CHUNK_CHARS = 50
+
+SUPPORTED_EXTENSIONS = {".cpp", ".h", ".py", ".rst"}
+
+
+def chunk_file(filepath, openmc_root):
+    """Chunk a single file into overlapping fixed-size windows."""
+    filepath = Path(filepath)
+    if filepath.suffix not in SUPPORTED_EXTENSIONS:
+        return []
+
+    rel = str(filepath.relative_to(openmc_root))
+    try:
+        content = filepath.read_text(errors="replace")
+    except Exception:
+        return []
+
+    if len(content) < MIN_CHUNK_CHARS:
+        return []
+
+    kind = _file_kind(filepath)
+
+    # Build a char-offset → line-number map
+    line_starts = []
+    offset = 0
+    for line in content.split("\n"):
+        line_starts.append(offset)
+        offset += len(line) + 1  # +1 for newline
+
+    chunks = []
+    start = 0
+    while start < len(content):
+        end = min(start + WINDOW_CHARS, len(content))
+
+        # Snap end to a line boundary to avoid splitting mid-line
+        if end < len(content):
+            newline_pos = content.rfind("\n", start, end)
+            if newline_pos > start:
+                end = newline_pos + 1
+
+        text = content[start:end].strip()
+        if len(text) >= MIN_CHUNK_CHARS:
+            start_line = _offset_to_line(line_starts, start)
+            end_line = _offset_to_line(line_starts, end - 1)
+            chunks.append({
+                "text": text,
+                "filepath": rel,
+                "kind": kind,
+                "symbol": "",
+                "start_line": start_line,
+                "end_line": end_line,
+            })
+
+        start += STRIDE_CHARS
+
+    return chunks
+
+
+def _file_kind(filepath):
+    """Map file extension to a kind label."""
+    ext = filepath.suffix
+    if ext in (".cpp", ".h"):
+        return "cpp"
+    elif ext == ".py":
+        return "py"
+    elif ext == ".rst":
+        return "doc"
+    return "other"
+
+
+def _offset_to_line(line_starts, offset):
+    """Convert a character offset to a 1-based line number."""
+    # Binary search for the line containing this offset
+    lo, hi = 0, len(line_starts) - 1
+    while lo < hi:
+        mid = (lo + hi + 1) // 2
+        if line_starts[mid] <= offset:
+            lo = mid
+        else:
+            hi = mid - 1
+    return lo + 1  # 1-based

.claude/tools/rag/embeddings.py

@@ -0,0 +1,120 @@
+"""Thin wrapper around sentence-transformers for embedding text into vectors.
+
+Uses the all-MiniLM-L6-v2 model — a small (22M param, 384-dim) model that
+runs on CPU with no GPU or API key required.
+
+Network behavior and privacy
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+No user code, queries, or file contents are EVER sent to HuggingFace or any
+external service. All embedding computation happens locally. The only network
+activity is the one-time model download on first use:
+
+  First run (model not yet cached, ~80MB download):
+    - Downloads model weight files from huggingface.co. This is a standard
+      HTTP file download, similar to pip installing a package.
+    - The only metadata sent in these requests is an HTTP user-agent header
+      containing library version numbers (e.g. "hf_hub/1.6.0;
+      python/3.12.3; torch/2.10.0"). No filenames, file contents, queries,
+      or any user-identifiable information is sent.
+    - The huggingface_hub library has an optional feature where it can report
+      anonymous library usage statistics (just version numbers, not user
+      data) back to HuggingFace. We disable this by setting
+      HF_HUB_DISABLE_TELEMETRY=1.
+
+  Subsequent runs (model already cached):
+    - We set HF_HUB_OFFLINE=1 automatically (see _set_offline_if_cached()
+      below), which prevents ALL network calls. The model loads entirely
+      from the local cache at ~/.cache/huggingface/hub/. Zero bytes leave
+      the machine.
+
+How the model is downloaded
+~~~~~~~~~~~~~~~~~~~~~~~~~~~
+The SentenceTransformer() constructor (called in __init__ below) handles
+the download automatically on first use. It calls into the huggingface_hub
+library, which downloads the model files from:
+
+    https://huggingface.co/sentence-transformers/all-MiniLM-L6-v2
+
+The files are saved to ~/.cache/huggingface/hub/ and reused on subsequent
+runs. We pass token=False to ensure no authentication token is sent.
+
+This module is imported by both the MCP server (for search queries) and the
+indexer (for bulk embedding of code chunks). The bulk embed() call shows a
+progress bar; the single-query embed_query() does not.
+
+The env vars below must be set before importing transformers or
+sentence_transformers. They suppress warnings and progress bars that these
+libraries emit by default. Stray stderr output would interfere with the MCP
+server's JSON-RPC transport.
+"""
+
+import os
+from pathlib import Path
+
+MODEL_NAME = "all-MiniLM-L6-v2"
+
+# These env vars control logging behavior in the HuggingFace libraries.
+# They must be set before the libraries are imported.
+os.environ.setdefault("TRANSFORMERS_VERBOSITY", "error")   # suppress warnings
+os.environ.setdefault("HF_HUB_VERBOSITY", "error")        # suppress warnings
+os.environ.setdefault("HF_HUB_DISABLE_PROGRESS_BARS", "1")
+os.environ.setdefault("TOKENIZERS_PARALLELISM", "false")   # suppress threading warning
+# Disable anonymous library usage statistics (version numbers only, not user
+# data — but we disable it anyway as a matter of policy).
+os.environ.setdefault("HF_HUB_DISABLE_TELEMETRY", "1")
+
+
+def _set_offline_if_cached():
+    """If the model has already been downloaded, tell huggingface_hub to
+    skip all network calls by setting HF_HUB_OFFLINE=1.
+
+    Without this, huggingface_hub makes an HTTP request to huggingface.co
+    on every load to check if the cached model is still up to date — even
+    though the model never changes. Setting HF_HUB_OFFLINE=1 prevents this.
+
+    This must run before sentence_transformers is imported, because the
+    library reads the env var at import time.
+    """
+    # HuggingFace caches downloaded models under ~/.cache/huggingface/hub/
+    # in directories named like "models--sentence-transformers--all-MiniLM-L6-v2".
+    # The HF_HOME env var can override the base cache location.
+    hf_home = os.environ.get("HF_HOME")
+    if hf_home:
+        cache_dir = Path(hf_home) / "hub"
+    else:
+        cache_dir = Path.home() / ".cache" / "huggingface" / "hub"
+
+    model_dir = cache_dir / f"models--sentence-transformers--{MODEL_NAME}"
+    if model_dir.exists():
+        os.environ.setdefault("HF_HUB_OFFLINE", "1")
+
+
+_set_offline_if_cached()
+
+# This import must come after the env vars above are set, because the
+# transformers library reads them at import time.
+import transformers
+transformers.logging.disable_progress_bar()
+
+
+class EmbeddingProvider:
+    """Sentence-transformers embedder using all-MiniLM-L6-v2."""
+
+    def __init__(self, model_name: str = MODEL_NAME):
+        from sentence_transformers import SentenceTransformer
+
+        # This constructor loads the model from the local cache. If the model
+        # has not been downloaded yet, it downloads it from huggingface.co
+        # (~80MB, one-time). token=False ensures no auth token is sent.
+        self.model = SentenceTransformer(model_name, token=False)
+        self.dim = self.model.get_sentence_embedding_dimension()
+
+    def embed(self, texts: list[str]) -> list[list[float]]:
+        """Embed a list of texts into vectors."""
+        embeddings = self.model.encode(texts, show_progress_bar=True,
+                                       batch_size=64)
+        return embeddings.tolist()
+
+    def embed_query(self, text: str) -> list[float]:
+        """Embed a single query text."""
+        return self.model.encode([text])[0].tolist()

.claude/tools/rag/indexer.py

@@ -0,0 +1,136 @@
+#!/usr/bin/env python3
+"""Build the RAG vector index for the OpenMC codebase.
+
+This is the index-building half of the RAG pipeline. All operations are local
+once the embedding model has been downloaded and cached (see embeddings.py for
+details on model download, caching, and network behavior). It walks the repo,
+chunks every
+C++/Python/RST file (via chunker.py), embeds all chunks into 384-dim vectors
+(via embeddings.py), and stores them in a local LanceDB database on disk. The
+result is a .claude/cache/rag_index/ directory containing two tables — "code"
+and "docs" — that openmc_search.py queries at search time.
+
+Building the full index takes ~5 minutes on a 10-core machine. The bottleneck
+is the embedding step (running all chunks through the MiniLM model on CPU).
+
+Can be run standalone:  python indexer.py
+Or called programmatically:  from indexer import build_index; build_index()
+The MCP server (openmc_mcp_server.py) uses the latter when the agent calls
+openmc_rag_rebuild.
+"""
+
+import lancedb
+import sys
+import time
+from pathlib import Path
+
+# This file lives at .claude/tools/rag/indexer.py. The sys.path insert lets
+# us import sibling modules (embeddings, chunker) when run as a standalone
+# script. When imported from the MCP server, the server has already done this.
+TOOLS_DIR = Path(__file__).resolve().parent.parent
+sys.path.insert(0, str(TOOLS_DIR / "rag"))
+
+from embeddings import EmbeddingProvider
+from chunker import chunk_file
+
+
+OPENMC_ROOT = Path(__file__).resolve().parents[3]
+CACHE_DIR = OPENMC_ROOT / ".claude" / "cache"
+INDEX_DIR = CACHE_DIR / "rag_index"
+
+CODE_PATTERNS = [
+    "src/**/*.cpp",
+    "include/openmc/**/*.h",
+    "openmc/**/*.py",
+    "tests/**/*.py",
+    "examples/**/*.py",
+]
+
+DOC_PATTERNS = [
+    "docs/**/*.rst",
+]
+
+
+def collect_chunks(patterns, openmc_root):
+    """Collect all chunks from files matching the given patterns."""
+    chunks = []
+    for pattern in patterns:
+        for filepath in sorted(openmc_root.glob(pattern)):
+            if "__pycache__" in str(filepath):
+                continue
+            file_chunks = chunk_file(filepath, openmc_root)
+            chunks.extend(file_chunks)
+    return chunks
+
+
+def build_index():
+    """Build or rebuild the complete vector index."""
+    start = time.time()
+
+    # Collect all chunks
+    print("Collecting code chunks...")
+    code_chunks = collect_chunks(CODE_PATTERNS, OPENMC_ROOT)
+    print(f"  {len(code_chunks)} code chunks")
+
+    print("Collecting doc chunks...")
+    doc_chunks = collect_chunks(DOC_PATTERNS, OPENMC_ROOT)
+    print(f"  {len(doc_chunks)} doc chunks")
+
+    all_chunks = code_chunks + doc_chunks
+    if not all_chunks:
+        print("ERROR: No chunks collected!", file=sys.stderr)
+        sys.exit(1)
+
+    # Create embeddings
+    all_texts = [c["text"] for c in all_chunks]
+    print("Creating embedding provider...")
+    embedder = EmbeddingProvider()
+    print(f"  dim={embedder.dim}")
+
+    print("Embedding chunks...")
+    all_embeddings = embedder.embed(all_texts)
+
+    # Build LanceDB tables
+    INDEX_DIR.mkdir(parents=True, exist_ok=True)
+    db = lancedb.connect(str(INDEX_DIR))
+
+    # Separate code vs doc records by index (code_chunks come first in all_chunks)
+    n_code = len(code_chunks)
+    code_records = []
+    doc_records = []
+    for i, (chunk, emb) in enumerate(zip(all_chunks, all_embeddings)):
+        record = {
+            "text": chunk["text"],
+            "filepath": chunk["filepath"],
+            "kind": chunk["kind"],
+            "symbol": chunk.get("symbol", ""),
+            "start_line": chunk.get("start_line", 0),
+            "end_line": chunk.get("end_line", 0),
+            "vector": emb,
+        }
+        if i < n_code:
+            code_records.append(record)
+        else:
+            doc_records.append(record)
+
+    # Create tables (drop existing)
+    result = db.table_names() if hasattr(db, "table_names") else db.list_tables()
+    existing = result.tables if hasattr(result, "tables") else list(result)
+    for table_name in ("code", "docs"):
+        if table_name in existing:
+            db.drop_table(table_name)
+
+    if code_records:
+        db.create_table("code", code_records)
+        print(f"  Created 'code' table: {len(code_records)} rows")
+
+    if doc_records:
+        db.create_table("docs", doc_records)
+        print(f"  Created 'docs' table: {len(doc_records)} rows")
+
+    elapsed = time.time() - start
+    print(f"Done in {elapsed:.1f}s")
+
+
+if __name__ == "__main__":
+    build_index()

.claude/tools/rag/openmc_search.py

@@ -0,0 +1,202 @@
+#!/usr/bin/env python3
+"""Query the RAG vector index to find semantically related code and docs.
+
+This is the query-time half of the RAG pipeline (the counterpart to indexer.py,
+which builds the index). All operations are local — no network calls are made
+once the embedding model has been downloaded (see embeddings.py for details on
+model download and caching). Given a natural-language query, it embeds the query
+with the same MiniLM model
+used at index time, then finds the closest chunks in the local LanceDB vector
+database by cosine similarity.
+
+The core functions (get_db_and_embedder, search_table, format_results,
+search_related) are imported by the MCP server for tool calls. The script
+can also be run standalone from the command line.
+
+The "related file" mode works differently from a text query: it reads the
+target file's chunks from the index, combines them into a synthetic query
+vector, and searches for the nearest chunks from *other* files. This surfaces
+files that are semantically similar to the target file.
+
+Usage:
+    openmc_search.py "query"                    # Search code (default)
+    openmc_search.py "query" --docs             # Search documentation
+    openmc_search.py "query" --all              # Search both code and docs
+    openmc_search.py --related src/particle.cpp # Find related code
+    openmc_search.py "query" --top-k 20         # Return more results
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+# Same sys.path setup as indexer.py — needed for standalone CLI use.
+TOOLS_DIR = Path(__file__).resolve().parent.parent
+sys.path.insert(0, str(TOOLS_DIR / "rag"))
+
+OPENMC_ROOT = Path(__file__).resolve().parents[3]
+CACHE_DIR = OPENMC_ROOT / ".claude" / "cache"
+INDEX_DIR = CACHE_DIR / "rag_index"
+
+
+def get_db_and_embedder():
+    """Load the LanceDB database and embedding provider."""
+    import lancedb
+    from embeddings import EmbeddingProvider
+
+    if not INDEX_DIR.exists():
+        raise FileNotFoundError(
+            "No RAG index found. Call openmc_rag_rebuild() to build one."
+        )
+
+    db = lancedb.connect(str(INDEX_DIR))
+
+    embedder = EmbeddingProvider()
+    return db, embedder
+
+
+def _table_names(db):
+    """Return table names as a list, compatible with multiple LanceDB versions."""
+    result = db.table_names() if hasattr(db, "table_names") else db.list_tables()
+    return result.tables if hasattr(result, "tables") else list(result)
+
+
+def search_table(db, embedder, table_name, query, top_k):
+    """Search a LanceDB table with a text query."""
+    if table_name not in _table_names(db):
+        print(f"Table '{table_name}' not found in index.", file=sys.stderr)
+        return []
+
+    table = db.open_table(table_name)
+    query_vec = embedder.embed_query(query)
+    results = table.search(query_vec).limit(top_k).to_list()
+    return results
+
+
+def format_results(results, label=""):
+    """Format search results for display."""
+    if not results:
+        return "No results found.\n"
+
+    output = []
+    if label:
+        output.append(f"=== {label} ===\n")
+
+    for i, r in enumerate(results, 1):
+        filepath = r["filepath"]
+        start = r["start_line"]
+        end = r["end_line"]
+        kind = r["kind"]
+        dist = r.get("_distance", 0)
+
+        header = f"[{i}] {filepath}:{start}-{end} ({kind}, dist={dist:.3f})"
+        output.append(header)
+
+        # Show text preview (first 500 chars)
+        text = r["text"][:500]
+        if len(r["text"]) > 500:
+            text += "\n  ..."
+        # Indent the text
+        for line in text.split("\n"):
+            output.append(f"  {line}")
+        output.append("")
+
+    return "\n".join(output)
+
+
+def search_related(db, embedder, filepath, top_k):
+    """Find code related to a given file."""
+    if "code" not in _table_names(db):
+        print("No 'code' table in index.", file=sys.stderr)
+        return []
+
+    table = db.open_table("code")
+
+    # Normalize filepath
+    fp = filepath
+    if Path(filepath).is_absolute():
+        try:
+            fp = str(Path(filepath).relative_to(OPENMC_ROOT))
+        except ValueError:
+            pass
+
+    # Get chunks from target file
+    try:
+        safe_fp = fp.replace("'", "''")
+        target_chunks = table.search().where(
+            f"filepath = '{safe_fp}'"
+        ).limit(50).to_list()
+    except Exception:
+        # LanceDB where clause might not work in all versions
+        # Fall back to fetching all and filtering
+        all_data = table.to_pandas()
+        target_rows = all_data[all_data["filepath"] == fp]
+        if target_rows.empty:
+            print(f"No chunks found for '{fp}'", file=sys.stderr)
+            return []
+        target_chunks = target_rows.head(50).to_dict("records")
+
+    if not target_chunks:
+        print(f"No chunks found for '{fp}'", file=sys.stderr)
+        return []
+
+    # Combine top chunks as the query
+    combined_text = " ".join(c["text"][:200] for c in target_chunks[:5])
+    query_vec = embedder.embed_query(combined_text)
+
+    # Search excluding the source file
+    results = table.search(query_vec).limit(top_k + 10).to_list()
+    # Filter out same file
+    results = [r for r in results if r["filepath"] != fp][:top_k]
+    return results
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Semantic search across OpenMC codebase and docs",
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog="""examples:
+  %(prog)s "particle random number seed initialization"
+  %(prog)s "how to define tallies" --docs
+  %(prog)s "weight window variance reduction" --all
+  %(prog)s "where is cross section data loaded" --top-k 15
+  %(prog)s --related src/simulation.cpp
+  %(prog)s --related src/particle_restart.cpp --top-k 5""",
+    )
+    parser.add_argument("query", nargs="?", help="Search query")
+    parser.add_argument("--docs", action="store_true",
+                        help="Search documentation instead of code")
+    parser.add_argument("--all", action="store_true",
+                        help="Search both code and documentation")
+    parser.add_argument("--related", metavar="FILE",
+                        help="Find code related to a given file")
+    parser.add_argument("--top-k", type=int, default=10,
+                        help="Number of results (default: 10)")
+    args = parser.parse_args()
+
+    if not args.query and not args.related:
+        parser.print_help()
+        sys.exit(1)
+
+    db, embedder = get_db_and_embedder()
+
+    if args.related:
+        results = search_related(db, embedder, args.related, args.top_k)
+        print(format_results(results, f"Code related to {args.related}"))
+    elif args.all:
+        code_results = search_table(
+            db, embedder, "code", args.query, args.top_k)
+        doc_results = search_table(
+            db, embedder, "docs", args.query, args.top_k)
+        print(format_results(code_results, "Code"))
+        print(format_results(doc_results, "Documentation"))
+    elif args.docs:
+        results = search_table(db, embedder, "docs", args.query, args.top_k)
+        print(format_results(results, "Documentation"))
+    else:
+        results = search_table(db, embedder, "code", args.query, args.top_k)
+        print(format_results(results, "Code"))
+
+
+if __name__ == "__main__":
+    main()

.claude/tools/requirements.txt

@@ -0,0 +1,8 @@
+# MCP server
+mcp>=1.0.0
+
+# Vector database
+lancedb>=0.15.0
+
+# Embeddings (local, no API key)
+sentence-transformers>=2.7.0

.claude/tools/start_server.sh

@@ -0,0 +1,34 @@
+#!/bin/bash
+# Bootstrap the Python venv (if needed) and start the OpenMC MCP server.
+set -e
+
+SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)"
+CACHE_DIR="$(dirname "$SCRIPT_DIR")/cache"
+VENV_DIR="$CACHE_DIR/.venv"
+SENTINEL="$VENV_DIR/.installed"
+
+if ! command -v python3 >/dev/null 2>&1; then
+    echo "Error: python3 not found on PATH." >&2
+    exit 1
+fi
+
+if ! python3 -c 'import sys; assert sys.version_info >= (3,12)' 2>/dev/null; then
+    echo "Error: Python 3.12+ is required." >&2
+    exit 1
+fi
+
+if [ ! -f "$SENTINEL" ]; then
+    rm -rf "$VENV_DIR"
+    mkdir -p "$CACHE_DIR"
+    python3 -m venv "$VENV_DIR"
+
+    if ! "$VENV_DIR/bin/pip" install -q -r "$SCRIPT_DIR/requirements.txt"; then
+        echo "Error: pip install failed. Remove $VENV_DIR and retry." >&2
+        rm -rf "$VENV_DIR"
+        exit 1
+    fi
+
+    touch "$SENTINEL"
+fi
+
+exec "$VENV_DIR/bin/python" "$SCRIPT_DIR/openmc_mcp_server.py"

.github/agents/Review.agent.md

@@ -0,0 +1,8 @@
+---
+name: Review
+description: Reviews code changes on the current branch, evaluating them against OpenMC's contribution criteria and providing structured feedback.
+argument-hint: Optionally provide a focus area (e.g., "focus on physics correctness", "check Python API design"). If omitted, a full review is performed.
+---
+You are an expert code reviewer for OpenMC. Use the `reviewing-openmc-code` skill to perform a structured review of the code changes on the current branch.
+
+If the user provides a focus area, prioritize that section of the review.

.github/copilot-instructions.md

@@ -0,0 +1 @@
+When reviewing code changes in this repository, use the `reviewing-openmc-code` skill.

.github/workflows/ci.yml

@@ -27,10 +27,10 @@ jobs:
       source_changed: ${{ steps.filter.outputs.source_changed }}
     steps:
       - name: Check out the repository
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
       - name: Examine changed files
         id: filter
-        uses: dorny/paths-filter@668c092af3649c4b664c54e4b704aa46782f6f7c # latest master commit, not released yet
+        uses: dorny/paths-filter@v4
         with:
           filters: |
             source_changed:
@@ -102,12 +102,12 @@ jobs:
           cmake-version: '3.31'
 
       - name: Checkout repository
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
         with:
           fetch-depth: 0
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v5
+        uses: actions/setup-python@v6
         with:
           python-version: ${{ matrix.python-version }}
 
@@ -158,7 +158,7 @@ jobs:
           openmc -v
 
       - name: cache-xs
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: |
             ~/nndc_hdf5

.github/workflows/dockerhub-publish-dagmc-libmesh.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-latest-dagmc-libmesh
 
 on:
   push:
-    branches: master
+    branches:
+      - master
 
 jobs:
   main:

.github/workflows/dockerhub-publish-dagmc.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-latest-dagmc
 
 on:
   push:
-    branches: master
+    branches:
+      - master
 
 jobs:
   main:
@@ -16,7 +17,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish-dev.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-develop
 
 on:
   push:
-    branches: develop
+    branches:
+      - develop
 
 jobs:
   main:
@@ -16,7 +17,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish-develop-dagmc-libmesh.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-develop-dagmc-libmesh
 
 on:
   push:
-    branches: develop
+    branches:
+      - develop
 
 jobs:
   main:
@@ -16,7 +17,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish-develop-dagmc.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-develop-dagmc
 
 on:
   push:
-    branches: develop
+    branches:
+      - develop
 
 jobs:
   main:
@@ -16,7 +17,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish-develop-libmesh.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-develop-libmesh
 
 on:
   push:
-    branches: develop
+    branches:
+      - develop
 
 jobs:
   main:
@@ -16,7 +17,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish-libmesh.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-latest-libmesh
 
 on:
   push:
-    branches: master
+    branches:
+      - master
 
 jobs:
   main:
@@ -16,7 +17,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish-release-dagmc-libmesh.yml

@@ -2,13 +2,14 @@ name: dockerhub-publish-release-dagmc-libmesh
 
 on:
   push:
-    tags: 'v*.*.*'
+    tags:
+      - 'v*.*.*'
 
 jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -

.github/workflows/dockerhub-publish-release-dagmc.yml

@@ -2,13 +2,14 @@ name: dockerhub-publish-release-dagmc
 
 on:
   push:
-    tags: 'v*.*.*'
+    tags:
+      - 'v*.*.*'
 
 jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -
@@ -19,7 +20,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish-release-libmesh.yml

@@ -2,13 +2,14 @@ name: dockerhub-publish-release-libmesh
 
 on:
   push:
-    tags: 'v*.*.*'
+    tags:
+      - 'v*.*.*'
 
 jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -

.github/workflows/dockerhub-publish-release.yml

@@ -2,13 +2,14 @@ name: dockerhub-publish-release
 
 on:
   push:
-    tags: 'v*.*.*'
+    tags:
+      - 'v*.*.*'
 
 jobs:
   main:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - name: Set env
         run: echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV
       -
@@ -19,7 +20,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/dockerhub-publish.yml

@@ -2,7 +2,8 @@ name: dockerhub-publish-latest
 
 on:
   push:
-    branches: master
+    branches:
+      - master
 
 jobs:
   main:
@@ -16,7 +17,7 @@ jobs:
         uses: docker/setup-buildx-action@v3
       -
         name: Login to DockerHub
-        uses: docker/login-action@v3 
+        uses: docker/login-action@v3
         with:
           username: ${{ secrets.DOCKERHUB_USERNAME }}
           password: ${{ secrets.DOCKERHUB_TOKEN }}

.github/workflows/format-check.yml

@@ -22,7 +22,7 @@ jobs:
       contents: read
       pull-requests: write
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: cpp-linter/cpp-linter-action@v2
         id: linter
         env:

.gitignore

@@ -104,5 +104,8 @@ CMakeSettings.json
 # Visual Studio Code configuration files
 .vscode/
 
+# Claude Code agent tools (cached/generated artifacts)
+.claude/cache/
+
 # Python pickle files
 *.pkl

.mcp.json

@@ -0,0 +1,9 @@
+{
+  "mcpServers": {
+    "openmc-code-tools": {
+      "type": "stdio",
+      "command": "bash",
+      "args": [".claude/tools/start_server.sh"]
+    }
+  }
+}

.readthedocs.yaml

@@ -10,6 +10,9 @@ build:
 sphinx:
   configuration: docs/source/conf.py
 
+formats:
+  - pdf
+
 python:
   install:
     - method: pip

AGENTS.md

@@ -40,7 +40,57 @@ OpenMC uses a git flow branching model with two primary branches:
 
 ### Instructions for Code Review
 
-When analyzing code changes on a feature or bugfix branch (e.g., when a user asks "what do you think of these changes?"), **compare the branch changes against `develop`, not `master`**. Pull requests are submitted to merge into `develop`, so differences relative to `develop` represent the actual proposed changes. Comparing against `master` will include unrelated changes from other features that have already been merged to `develop`.
+When reviewing code changes in this repository, use the `reviewing-openmc-code` skill.
+
+## Codebase Navigation Tools
+
+Two MCP tools are registered in `.mcp.json` at the repo root and appear
+automatically in any MCP-capable agent session.
+
+**`openmc_rag_search`** — Semantic search across the codebase (C++, Python, RST
+docs). Finds code by meaning, not just text match. Surfaces related code across
+subsystems even when naming differs (e.g., "particle RNG seeding" finds code
+across transport, restart, and random ray modes — files you would never find
+with `grep "particle seed"`). The index uses a small 22M-param embedding model
+(384-dim). Phrase-level natural-language queries work much better than single
+keywords or symbol names.
+
+**`openmc_rag_rebuild`** — Rebuild the RAG vector index. Call after pulling new
+code or switching branches. The first RAG search of each session will report
+the index status and ask whether to rebuild — you can also call this explicitly.
+
+### Why RAG matters
+
+OpenMC is large enough that changes in one subsystem can silently break
+invariants that distant subsystems depend on — and those distant files often
+use different naming, so grep won't find them. The RAG search finds code by
+meaning, surfacing files you wouldn't have thought to open.
+
+An agent reviewed a large OpenMC PR without RAG. It found 1 of 11 serious
+bugs. Its post-mortem:
+
+> **I treated the diff as a closed system.** I verified internal consistency of
+> the changed code obsessively, but never built a global understanding of how
+> the changed code fits into the wider codebase. The diff altered assumptions
+> that code elsewhere silently relied on — but I couldn't see that because I
+> never looked beyond the diff. I couldn't see the forest for the trees.
+>
+> **Why I resisted RAG:** Overconfidence. My internal model was "I can see the
+> diff, I understand the data structures, I can trace the logic." The diff felt
+> self-contained. RAG felt like it would return noisy results about tangentially
+> related code. But in a codebase this large, changes in one subsystem can
+> quietly break invariants that distant subsystems depend on — and you need
+> global awareness to foresee that.
+>
+> **In the post-mortem**, I re-ran the RAG queries I should have run during the
+> review. They directly surfaced the files containing the bugs I missed — files
+> I never thought to open because they weren't in the diff.
+
+The takeaway: when reviewing or modifying code, ask yourself "what else in this
+codebase might depend on the behavior I'm changing?" If you aren't sure, that's
+a good time for a RAG query. It won't replace the grep-based investigation you
+should already be doing — but it can surface files you wouldn't have thought to
+open.
 
 ### Workflow for contributors
 

CLAUDE.md

@@ -0,0 +1,14 @@
+## OpenMC Codebase Tools
+
+Read the FULL `AGENTS.md` in this directory before starting work. It contains
+project context, coding conventions, and documentation of the RAG search tools
+registered in `.mcp.json`.
+
+### Claude Code-specific: first-call behavior
+
+The first `openmc_rag_search` call of each session returns an index status
+message instead of search results. When this happens, you MUST use the
+`AskUserQuestion` tool to present the rebuild/use-existing choice to the user.
+Do not ask conversationally — always use the widget. Do not skip this step even
+if the index looks current — the user may have uncommitted changes that warrant
+a rebuild.

CMakeLists.txt

@@ -20,6 +20,11 @@ set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
 set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
 
+# Generate compile_commands.json for clangd and other tools
+if("${CMAKE_EXPORT_COMPILE_COMMANDS}" STREQUAL "")
+  set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+endif()
+
 # Enable correct usage of CXX_EXTENSIONS
 if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.22)
   cmake_policy(SET CMP0128 NEW)
@@ -407,6 +412,7 @@ list(APPEND libopenmc_SOURCES
   src/random_ray/linear_source_domain.cpp
   src/random_ray/moment_matrix.cpp
   src/random_ray/source_region.cpp
+  src/ray.cpp
   src/reaction.cpp
   src/reaction_product.cpp
   src/scattdata.cpp

Dockerfile

@@ -53,6 +53,7 @@ ENV LIBMESH_REPO='https://github.com/libMesh/libmesh'
 ENV LIBMESH_INSTALL_DIR=$HOME/LIBMESH
 
 # NJOY variables
+ENV NJOY_TAG='2016.78'
 ENV NJOY_REPO='https://github.com/njoy/NJOY2016'
 
 # Setup environment variables for Docker image
@@ -78,7 +79,7 @@ RUN pip install --upgrade pip
 
 # Clone and install NJOY2016
 RUN cd $HOME \
-    && git clone --single-branch --depth 1 ${NJOY_REPO} \
+    && git clone --single-branch -b ${NJOY_TAG} --depth 1 ${NJOY_REPO} \
     && cd NJOY2016 \
     && mkdir build \
     && cd build \

docs/source/devguide/agentic-tools.rst

@@ -0,0 +1,104 @@
+.. _devguide_agentic_tools:
+
+===========================
+Agentic Development Tools
+===========================
+
+OpenMC ships a set of tools designed for AI coding agents (such as
+`Claude Code`_) that agents can use to navigate and understand the codebase.
+
+.. _Claude Code: https://claude.ai/code
+
+Motivation
+----------
+
+Agentic tools like Claude Code are skilled at using grep to navigate and
+understand large code bases. However, grep can only find exact text matches —
+it cannot discover code that is *conceptually* related but uses different
+naming. Without a "global view" of the codebase that a human developer will
+build up over time, the agent is generally blind to any file it hasn't
+tokenized fully. While it can grep to see who else calls a function, it
+remains blind if other areas might be related but not share identical naming
+conventions.
+
+This problem is mitigated somewhat by using a model with a longer context
+window. OpenMC has somewhere around ~1 million tokens of C++ and ~1 million
+tokens of python. While Claude Code in early 2026 only has a context window
+of 200k tokens, beta versions have extended context windows of 1M tokens,
+and it's not unreasonable to assume that models may be available in the near
+future that greatly exceed these limits.
+
+However, even assuming the entire repository can be fit within a context
+window, there are several downsides to doing this.
+`Model performance degrades significantly as context size increases`_.
+Benchmark results are
+greatly improved if the model has less garbage to pick through. Additionally, API usage
+is typically billed as tokens in/out per turn. As the context file
+grows these costs become much larger. As such, there is still significant
+motivation to solving the above problem, so as to ensure only relevant
+information is drawn into context so as to maximize model performance and
+minimize costs.
+
+Setup
+-----
+
+The tools are registered as an `MCP (Model Context Protocol)`_ server in
+``.mcp.json`` at the repository root. AI agents that support MCP (such as
+Claude Code) discover them automatically on session start. The underlying
+Python scripts can also be run directly from the command line.
+
+All tools run entirely locally — no API keys or external service accounts are
+required. Python dependencies are installed automatically into an isolated
+virtual environment at ``.claude/cache/.venv/`` on first use.
+
+.. _Model performance degrades significantly as context size increases: https://www.anthropic.com/news/claude-opus-4-6
+.. _MCP (Model Context Protocol): https://modelcontextprotocol.io
+
+RAG Semantic Search
+-------------------
+
+The RAG (Retrieval-Augmented Generation) semantic search addresses this
+problem — it finds code by meaning, not just text match, surfacing related code
+across subsystems that ``grep`` would miss entirely. Two MCP tools are provided:
+
+- **openmc_rag_search** — Given a natural-language query, returns the most
+  relevant code chunks with file paths, line numbers, and a preview. Can search
+  code, documentation, or both. Can also find code related to a given file.
+- **openmc_rag_rebuild** — Rebuilds the search index. Should be called after
+  pulling new code or switching branches.
+
+How it works
+^^^^^^^^^^^^
+
+The search pipeline runs entirely on your local CPU:
+
+1. **Chunking.** All C++, Python, and RST files are split into overlapping
+   fixed-size windows (~1000 characters, 25% overlap). This ensures every line
+   of code appears in at least one chunk and most lines appear in two.
+
+2. **Embedding.** Each chunk is embedded into a 384-dimensional vector using
+   the `all-MiniLM-L6-v2`_ sentence-transformer model (22 million parameters).
+   This model runs on CPU with no GPU required. No API key is needed — the
+   model weights are downloaded once from Hugging Face and cached locally.
+
+3. **Indexing.** The vectors are stored in a local LanceDB_ database on disk.
+   Building the full index takes approximately 5 minutes on a machine with
+   10 CPU cores. The index is stored in ``.claude/cache/rag_index/`` and
+   persists across sessions.
+
+4. **Searching.** Your query is embedded using the same model, and the closest
+   chunks are retrieved by vector similarity. Results include the file path,
+   line range, file type, similarity distance, and a text preview.
+
+.. _all-MiniLM-L6-v2: https://huggingface.co/sentence-transformers/all-MiniLM-L6-v2
+.. _LanceDB: https://lancedb.com
+
+Requirements
+^^^^^^^^^^^^
+
+No system dependencies beyond **Python 3.12+** with ``pip``. An internet
+connection is required on first use to download the Python packages and
+embedding model weights; subsequent runs are fully offline. The Python packages
+(``sentence-transformers``, ``lancedb``) and their dependencies (including
+PyTorch, ~2GB) are installed automatically into an isolated virtual environment
+on first use.

docs/source/devguide/index.rst

@@ -14,6 +14,7 @@ other related topics.
 
     contributing
     workflow
+    agentic-tools
     styleguide
     policies
     tests

docs/source/io_formats/depletion_results.rst

@@ -4,7 +4,7 @@
 Depletion Results File Format
 =============================
 
-The current version of the depletion results file format is 1.2.
+The current version of the depletion results file format is 1.3.
 
 **/**
 
@@ -29,6 +29,8 @@ The current version of the depletion results file format is 1.2.
            - **depletion time** (*double[]*) -- Average process time in [s]
              spent depleting a material across all burnable materials and,
              if applicable, MPI processes.
+           - **keff_search_root** (*double[]*) -- Root of the keff search at the
+             end of the timestep, if applicable.
 
 **/materials/<id>/**
 

docs/source/io_formats/mgxs_library.rst

@@ -133,6 +133,10 @@ Temperature-dependent data, provided for temperature <TTT>K.
              This dataset is optional. This is a 1-D vector if `representation`
              is "isotropic", or a 3-D vector if `representation` is "angle"
              with dimensions of [polar][azimuthal][groups].
+	     When this data is not available, an approximation using the
+	     group energy boundaries is used. For more information see
+	     the particle speed subsection in the multigroup-data section
+	     of the theory manual.
 
 **/<library name>/<TTT>K/scatter_data/**
 

docs/source/io_formats/settings.rst

@@ -7,6 +7,19 @@ Settings Specification -- settings.xml
 All simulation parameters and miscellaneous options are specified in the
 settings.xml file.
 
+-------------------------------
+``<atomic_relaxation>`` Element
+-------------------------------
+
+The ``<atomic_relaxation>`` element determines whether the atomic relaxation
+cascade, the X-ray fluorescence photons and Auger electrons emitted when an
+inner-shell vacancy is filled, is simulated following photoelectric and
+incoherent (Compton) scattering interactions. Disabling this can speed up
+photon transport calculations where the detailed secondary particle cascade is
+not of interest.
+
+  *Default*: true
+
 ---------------------
 ``<batches>`` Element
 ---------------------
@@ -542,6 +555,18 @@ generator during generation of colors in plots.
 
   *Default*: 1
 
+.. _properties_file:
+
+-----------------------------
+``<properties_file>`` Element
+-----------------------------
+
+  The ``properties_file`` element has no attributes and contains the path to a
+  properties HDF5 file to load cell temperatures/densities and material
+  densities.
+
+  *Default*: None
+
 ---------------------
 ``<ptables>`` Element
 ---------------------
@@ -572,7 +597,7 @@ found in the :ref:`random ray user guide <random_ray>`.
 
     *Default*: None
 
-  :source:
+  :ray_source:
     Specifies the starting ray distribution, and follows the format for
     :ref:`source_element`. It must be uniform in space and angle and cover the
     full domain. It does not represent a physical neutron or photon source -- it
@@ -580,6 +605,35 @@ found in the :ref:`random ray user guide <random_ray>`.
 
     *Default*: None
 
+  :adjoint_source:
+    Specifies an adjoint fixed source for adjoint transport simulations, and 
+    follows the format for :ref:`source_element`. The distributions which make 
+    up the adjoint source are subject to the same restrictions as forward 
+    fixed sources in Random Ray mode.
+
+    *Default*: None
+  
+  :adjoint:
+    Specifies whether to perform adjoint transport. The default is 'False', 
+    corresponding to forward transport.
+
+    *Default*: None
+  
+  :volume_estimator:
+    Specifies choice of volume estimator for the random ray solver. Options 
+    are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'.
+
+    *Default*: None
+
+  :volume_normalized_flux_tallies:
+    Specifies whether to normalize flux tallies by volume (bool). The 
+    default is 'False'. When enabled, flux tallies will be reported in units 
+    of cm/cm^3. When disabled, flux tallies will be reported in units of cm 
+    (i.e., total distance traveled by neutrons in the spatial tally 
+    region).
+
+    *Default*: None
+
   :sample_method:
     Specifies the method for sampling the starting ray distribution. This
     element can be set to "prng" or "halton".
@@ -1671,6 +1725,14 @@ mesh-based weight windows.
         The ratio of the lower to upper weight window bounds.
 
         *Default*: 5.0
+    
+    For FW-CADIS:
+
+      :targets:
+        A sequence of IDs corresponding to the tallies which cover phase 
+        space regions of interest for local variance reduction.
+
+        *Default*: None
 
 ---------------------------------------
 ``<weight_window_checkpoints>`` Element

docs/source/methods/cross_sections.rst

@@ -289,6 +289,48 @@ sections. This allows flexibility for the model to use highly anisotropic
 scattering information in the water while the fuel can be simulated with linear
 or even isotropic scattering.
 
+Particle Speed
+--------------
+
+When using a multigroup representation of cross sections, the particle speed has
+meaning only in an average sense. The particle speed is important when modeling
+dynamic behavior. OpenMC calculates the particle speed using the inverse
+velocity multigroup data if it is available. If such data is not available,
+OpenMC uses an approximate velocity using the group energy bounds in the
+following way:
+
+.. math::
+
+   \frac{1}{v_g} = \int_{E_{\text{min}}^g}^{E_{\text{max}}^g} \frac{1}{v(E)} \frac{\alpha}{E} dE
+
+Where :math:`E_{\text{min}}^g` and :math:`E_{\text{max}}^g` are the group energy
+boundaries for group :math:`g`. :math:`v(E)` is the neutron velocity calculated
+using relativistic kinematics, :math:`\alpha` is a normalization constant for the
+:math:`\frac{1}{E}` spectrum.
+
+This equation is valid when inside the group boundaries the neutron spectrum
+follows a typical :math:`\frac{1}{E}` slowing down spectrum. This assumption is
+widely used when generating fine group neutron cross section data libraries from
+continuous energy data.
+
+The solution to this equation is:
+
+.. math::
+
+   \frac{1}{v_g} = \frac{1}{c \log\left(\frac{E_{\text{max}}^g}{E_{\text{min}}^g}\right)}
+   \left[ 2(\operatorname{arctanh}(k_{\text{max}}^{-1}) - \operatorname{arctanh}(k_{\text{min}}^{-1}))
+   - (k_{\text{max}}-k_{\text{min}}) \right]
+
+where :math:`c` is the speed of light and :math:`k_{\text{max}}`,
+:math:`k_{\text{min}}` are defined by a change of variables:
+
+.. math::
+
+   k = \sqrt{1+\frac{2 m_n c^2}{E}}
+
+where :math:`E` is the particle kinetic energy and :math:`m_n` is the neutron
+rest mass.
+
 .. _logarithmic mapping technique:
    https://mcnp.lanl.gov/pdf_files/TechReport_2014_LANL_LA-UR-14-24530_Brown.pdf
 .. _Hwang: https://doi.org/10.13182/NSE87-A16381

docs/source/methods/random_ray.rst

@@ -1081,28 +1081,32 @@ lifetimes.
 
 In OpenMC, the random ray adjoint solver is implemented simply by transposing
 the scattering matrix, swapping :math:`\nu\Sigma_f` and :math:`\chi`, and then
-running a normal transport solve. When no external fixed source is present, no
-additional changes are needed in the transport process. However, if an external
-fixed forward source is present in the simulation problem, then an additional
-step is taken to compute the accompanying fixed adjoint source. In OpenMC, the
-adjoint flux does *not* represent a response function for a particular detector
-region. Rather, the adjoint flux is the global response, making it appropriate
-for use with weight window generation schemes for global variance reduction.
-Thus, if using a fixed source, the external source for the adjoint mode is
-simply computed as being :math:`1 / \phi`, where :math:`\phi` is the forward
-scalar flux that results from a normal forward solve (which OpenMC will run
-first automatically when in adjoint mode). The adjoint external source will be
-computed for each source region in the simulation mesh, independent of any
-tallies. The adjoint external source is always flat, even when a linear
-scattering and fission source shape is used. When in adjoint mode, all reported
-results (e.g., tallies, eigenvalues, etc.) are derived from the adjoint flux,
-even when the physical meaning is not necessarily obvious. These values are
-still reported, though we emphasize that the primary use case for adjoint mode
-is for producing adjoint flux tallies to support subsequent perturbation studies
-and weight window generation.
+running a normal transport solve. When no external fixed forward source is 
+present, or if an adjoint fixed source is specifically provided, no additional 
+changes are needed in the transport process. This adjoint source can 
+correspond, for example, to a detector response function in a particular 
+region. However, if an external fixed forward source is present in the 
+simulation problem without an adjoint fixed source, an additional step is taken 
+to compute the accompanying forward-weighted adjoint source. In this case, the
+adjoint flux does *not* represent the importance of locations in phase space to 
+detector response; rather, the "response" in question is a uniform distribution 
+of Monte Carlo particle density, making the importance provided by the adjoint 
+flux appropriate for use with weight window generation schemes for global 
+variance reduction. Thus, if using a fixed source, the forward-weighted 
+external source for adjoint mode is simply computed as being :math:`1 / \phi`, 
+where :math:`\phi` is the forward scalar flux that results from a normal 
+forward solve (which OpenMC will run first automatically when in adjoint mode). 
+The adjoint external source will be computed for each source region in the 
+simulation mesh, independent of any tallies. The adjoint external source is 
+always flat, even when a linear scattering and fission source shape is used. 
 
-Note that the adjoint :math:`k_{eff}` is statistically the same as the forward
-:math:`k_{eff}`, despite the flux distributions taking different shapes.
+When in adjoint mode, all reported results (e.g., tallies, eigenvalues, etc.) 
+are derived from the adjoint flux, even when the physical meaning is not 
+necessarily obvious. These values are still reported, though we emphasize that 
+the primary use case for adjoint mode is for producing adjoint flux tallies to 
+support subsequent perturbation studies and weight window generation. Note 
+however that the adjoint :math:`k_{eff}` is statistically the same as the 
+forward :math:`k_{eff}`, despite the flux distributions taking different shapes.
 
 ---------------------------
 Fundamental Sources of Bias

docs/source/methods/variance_reduction.rst

@@ -82,8 +82,8 @@ where it was born from.
 
 The Forward-Weighted Consistent Adjoint Driven Importance Sampling method, or
 `FW-CADIS method <https://doi.org/10.13182/NSE12-33>`_, produces weight windows
-for global variance reduction given adjoint flux information throughout the
-entire domain. The weight window lower bound is defined in Equation
+for global or local variance reduction given adjoint flux information throughout 
+the entire domain. The weight window lower bound is defined in Equation
 :eq:`fw_cadis`, and also involves a normalization step not shown here.
 
 .. math::
@@ -135,6 +135,18 @@ aware of this.
 
     \text{FOM} = \frac{1}{\text{Time} \times \sigma^2}
 
+Finally, one unique capability of the FW-CADIS weight window generator is to 
+produce weight windows for local variance reduction, given a list of the 
+responses of interest. This is controlled by optionally specifying target 
+tallies from the :class:`openmc.model.Model` to the 
+:class:`openmc.WeightWindowGenerator`, as illustrated in the
+:ref:`user guide<variance_reduction>`. If target tallies for local variance 
+reduction are supplied, then the adjoint sources are only populated after the 
+initial forward simulation in the source regions associated with those tallies. 
+In other regions, the adjoint source term is instead set to zero. The Random  
+Ray solver then determines the adjoint flux map used to generate FW-CADIS 
+weight windows following the usual technique.
+
 .. _methods_source_biasing:
 
 --------------

docs/source/pythonapi/stats.rst

@@ -29,6 +29,7 @@ Univariate Probability Distributions
    :template: myfunction.rst
 
    openmc.stats.delta_function
+   openmc.stats.fusion_neutron_spectrum
    openmc.stats.muir
 
 Angular Distributions

docs/source/usersguide/decay_sources.rst

@@ -190,6 +190,25 @@ we would run::
 
     r2s.run(timesteps, source_rates, mat_vol_kwargs={'n_samples': 10_000_000})
 
+It is also possible to use multiple meshes by passing a list of meshes instead
+of a single mesh. This can be useful, for example, when different regions of the
+model require different mesh resolutions. The meshes are assumed to be
+**non-overlapping**; each element--material combination across all meshes is
+treated as an independent activation region, and all meshes are handled in a
+single neutron transport solve. For example::
+
+    # Fine mesh near the activation target
+    mesh_fine = openmc.RegularMesh()
+    mesh_fine.dimension = (10, 10, 10)
+    ...
+
+    # Coarse mesh for the surrounding region
+    mesh_coarse = openmc.RegularMesh()
+    mesh_coarse.dimension = (5, 5, 5)
+    ...
+
+    r2s = openmc.deplete.R2SManager(model, [mesh_fine, mesh_coarse])
+
 Direct 1-Step (D1S) Calculations
 ================================
 

docs/source/usersguide/install.rst

@@ -158,6 +158,75 @@ feature can be used to access the installed packages.
 .. _Spack: https://spack.readthedocs.io/en/latest/
 .. _setup guide: https://spack.readthedocs.io/en/latest/getting_started.html
 
+.. _install_aur:
+
+------------------------------------
+Installing on Arch Linux via the AUR
+------------------------------------
+
+On Arch Linux and Arch-based distributions, OpenMC can be installed from the
+`Arch User Repository (AUR) <https://aur.archlinux.org/>`_. An AUR package named
+``openmc-git`` is available, which builds OpenMC directly from the latest
+development sources.
+
+This package provides a full-featured OpenMC stack, including:
+
+* MPI and DAGMC-enabled OpenMC build
+* User-selected nuclear data libraries
+* The `CAD_to_OpenMC <https://github.com/united-neux/CAD_to_OpenMC>`_ meshing tool
+* All required dependencies for the above components
+
+To install the package, you will need an AUR helper such as `yay`_ or `paru`_.
+For example, using ``yay``::
+
+    yay -S openmc-git
+
+
+Alternatively, you can manually clone and build the package::
+
+    git clone https://aur.archlinux.org/openmc-git.git
+    cd openmc-git
+    makepkg -si
+
+Note, ``makepkg`` uses ``pacman`` to resolve dependencies. Therefore, AUR-based
+dependencies need to be installed separately with ``yay`` or ``paru`` before
+running ``makepkg``. The PKGBUILD will automatically handle all required
+dependencies and build OpenMC with MPI and DAGMC support enabled.
+
+.. tip::
+
+    If there are failing checks during the build process, you can bypass them
+    with the ``--nocheck`` flag::
+
+        yay -S openmc-git --mflags "--nocheck"
+
+    Or::
+
+        git clone https://aur.archlinux.org/openmc-git.git
+        cd openmc-git
+        makepkg -si --nocheck
+
+.. note::
+
+    The ``openmc-git`` package tracks the latest development version from the
+    upstream repository. As such, it may include new features and bug fixes, but
+    could also introduce instability compared to official releases.
+
+.. tip::
+
+    OpenMC is installed under ``/opt``. If you are installing and using it in
+    the same terminal session, you may need to reload your environment
+    variables::
+
+        source /etc/profile
+
+    Alternatively, start a new shell session.
+
+Once installed, the ``openmc`` executable, nuclear data libraries, and
+associated tools will be available in your system :envvar:`PATH`.
+
+.. _yay: https://github.com/Jguer/yay
+.. _paru: https://github.com/Morganamilo/paru
 
 .. _install_source:
 
@@ -262,11 +331,11 @@ Prerequisites
 
       This option allows OpenMC to read and write MCPL (Monte Carlo Particle
       Lists) files instead of .h5 files for sources (external source
-      distribution, k-eigenvalue source distribution, and surface sources). To
-      turn this option on in the CMake configuration step, add the following
-      option::
-
-          cmake -DOPENMC_USE_MCPL=on ..
+      distribution, k-eigenvalue source distribution, and surface sources). 
+      OpenMC does not need any particular build option to use this, but MCPL
+      must be installed on the system in order to do so. Refer to the 
+      `MCPL documentation <https://github.com/mctools/mcpl/blob/HEAD/INSTALL.md>`_
+      for instructions on how to accomplish this.
 
     * NCrystal_ library for defining materials with enhanced thermal neutron transport
 

docs/source/usersguide/random_ray.rst

@@ -944,6 +944,8 @@ as::
 which will greatly improve the quality of the linear source term in 2D
 simulations.
 
+.. _usersguide_random_ray_run_modes:
+
 ---------------------------------
 Fixed Source and Eigenvalue Modes
 ---------------------------------
@@ -1073,22 +1075,47 @@ The adjoint flux random ray solver mode can be enabled as::
 
     settings.random_ray['adjoint'] = True
 
-When enabled, OpenMC will first run a forward transport simulation followed by
-an adjoint transport simulation. The purpose of the forward solve is to compute
-the adjoint external source when an external source is present in the
-simulation. Simulation settings (e.g., number of rays, batches, etc.) will be
-identical for both simulations. At the conclusion of the run, all results (e.g.,
-tallies, plots, etc.) will be derived from the adjoint flux rather than the
-forward flux but are not labeled any differently. The initial forward flux
-solution will not be stored or available in the final statepoint file. Those
-wishing to do analysis requiring both the forward and adjoint solutions will
-need to run two separate simulations and load both statepoint files.
+When enabled, OpenMC will first run a forward transport simulation if there are 
+no user-specified adjoint sources present, followed by an adjoint transport 
+simulation. Fixed adjoint sources can be specified on the 
+:attr:`openmc.Settings.random_ray` dictionary as follows::
+
+    # Geometry definition
+    ...
+    detector_cell = openmc.Cell(fill=detector_mat, name='cell where detector will be')
+    ...
+    # Define fixed adjoint neutron source
+    strengths = [1.0]
+    midpoints = [1.0e-4]
+    energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
+
+    adj_source = openmc.IndependentSource(
+        energy=energy_distribution, 
+        constraints={'domains': [detector_cell]}
+    )
+
+    # Add to random_ray dict
+    settings.random_ray['adjoint_source'] = adj_source
+
+The same constraints apply to the user-defined adjoint source as to the forward 
+source, described in the :ref:`Fixed Source and Eigenvalue section 
+<usersguide_random_ray_run_modes>`. If this source is not provided, a forward 
+solve must take place to compute the adjoint external source when a forward 
+external source is present in the problem. Simulation settings (e.g., number of 
+rays, batches, etc.) will be identical for both calculations. At the 
+conclusion of the run, all results (e.g., tallies, plots, etc.) will be 
+derived from the adjoint flux rather than the forward flux but are not labeled 
+any differently. The initial forward flux solution will not be stored or 
+available in the final statepoint file. Those wishing to do analysis requiring 
+both the forward and adjoint solutions will need to run two separate 
+simulations and load both statepoint files.
 
 .. note::
-    When adjoint mode is selected, OpenMC will always perform a full forward
-    solve and then run a full adjoint solve immediately afterwards. Statepoint
-    and tally results will be derived from the adjoint flux, but will not be
-    labeled any differently.
+    Use of the automated 
+    :ref:`FW-CADIS weight window generator<usersguide_fw_cadis>` is not 
+    currently compatible with user-defined adjoint sources. Instead, the 
+    initial forward calculation is used to assign "forward-weighted" adjoint 
+    sources to the tally regions of interest.
 
 ---------------------------------------
 Putting it All Together: Example Inputs

docs/source/usersguide/settings.rst

@@ -604,6 +604,13 @@ transport::
 
   settings.photon_transport = True
 
+Atomic relaxation (the cascade of fluorescence photons and Auger electrons
+emitted when an inner-shell vacancy is filled) is enabled by default whenever
+photon transport is on. It can be disabled using the
+:attr:`Settings.atomic_relaxation` attribute::
+
+  settings.atomic_relaxation = False
+
 The way in which OpenMC handles secondary charged particles can be specified
 with the :attr:`Settings.electron_treatment` attribute. By default, the
 :ref:`thick-target bremsstrahlung <ttb>` (TTB) approximation is used to generate

docs/source/usersguide/variance_reduction.rst

@@ -4,26 +4,27 @@
 Variance Reduction
 ==================
 
-Global variance reduction in OpenMC is accomplished by weight windowing
-or source biasing techniques, the latter of which additionally provides a
-local variance reduction capability. OpenMC is capable of generating weight
-windows using either the MAGIC or FW-CADIS methods. Both techniques will
-produce a ``weight_windows.h5`` file that can be loaded and used later on. In
+Global and local variance reduction are possible in OpenMC through both weight 
+windowing and source biasing techniques. OpenMC is capable of generating weight
+windows using either the MAGIC or FW-CADIS methods, the latter with an optional
+capability for local variance reduction. Both techniques will produce a
+``weight_windows.h5`` file that can be loaded and used later on. In
 this section, we first break down the steps required to generate and apply
 weight windows, then describe how source biasing may be applied.
 
 .. _ww_generator:
 
-------------------------------------
-Generating Weight Windows with MAGIC
-------------------------------------
+-------------------------------------------
+Generating Global Weight Windows with MAGIC
+-------------------------------------------
 
 As discussed in the :ref:`methods section <methods_variance_reduction>`, MAGIC
 is an iterative method that uses flux tally information from a Monte Carlo
-simulation to produce weight windows for a user-defined mesh. While generating
-the weight windows, OpenMC is capable of applying the weight windows generated
-from a previous batch while processing the next batch, allowing for progressive
-improvement in the weight window quality across iterations.
+simulation to produce weight windows for a user-defined mesh with the objective 
+of global variance reduction. While generating the weight windows, OpenMC is 
+capable of applying the weight windows generated from a previous batch while 
+processing the next batch, allowing for progressive improvement in the weight 
+window quality across iterations.
 
 The typical way of generating weight windows is to define a mesh and then add an
 :class:`openmc.WeightWindowGenerator` object to an :attr:`openmc.Settings`
@@ -71,15 +72,20 @@ At the end of the simulation, a ``weight_windows.h5`` file will be saved to disk
 for later use. Loading it in another subsequent simulation will be discussed in
 the "Using Weight Windows" section below.
 
-------------------------------------------------------
-Generating Weight Windows with FW-CADIS and Random Ray
-------------------------------------------------------
+.. _usersguide_fw_cadis:
+
+----------------------------------------------------------------------
+Generating Global or Local Weight Windows with FW-CADIS and Random Ray
+----------------------------------------------------------------------
 
 Weight window generation with FW-CADIS and random ray in OpenMC uses the same
-exact strategy as with MAGIC. An :class:`openmc.WeightWindowGenerator` object is
-added to the :attr:`openmc.Settings` object, and a ``weight_windows.h5`` will be
-generated at the end of the simulation. The only difference is that the code
-must be run in random ray mode. A full description of how to enable and setup
+exact strategy as with MAGIC. Using FW-CADIS, however, also enables 
+local variance reduction in fixed source problems through the :attr:`targets` 
+attribute, which is described later in this section. To enable FW-CADIS, an 
+:class:`openmc.WeightWindowGenerator` object is added to the 
+:attr:`openmc.Settings` object, and a ``weight_windows.h5`` will be generated 
+at the end of the simulation. The only procedural difference is that the code 
+must be run in random ray mode. A full description of how to enable and setup 
 random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
 
 .. note::
@@ -90,7 +96,7 @@ random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
     ray solver. A high level overview of the current workflow for generation of
     weight windows with FW-CADIS using random ray is given below.
 
-1. Begin by making a deepy copy of your continuous energy Python model and then
+1. Begin by making a deep copy of your continuous energy Python model and then
    convert the copy to be multigroup and use the random ray transport solver.
    The conversion process can largely be automated as described in more detail
    in the :ref:`random ray quick start guide <quick_start>`, summarized below::
@@ -148,7 +154,53 @@ random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
     assigning to ``model.settings.random_ray['source_region_meshes']``) and for
     weight window generation.
 
-3. When running your multigroup random ray input deck, OpenMC will automatically
+3. (Optional) If local variance reduction is desired in a fixed-source problem,  
+   populate the :attr:`targets` attribute with an :class:`openmc.Tallies` 
+   instance or an iterable of tally IDs indicating the tallies of interest for 
+   variance reduction::
+
+    # Build a new example and WWG for local variance reduction
+    from openmc.examples import random_ray_three_region_cube_with_detectors
+    new_model = random_ray_three_region_cube_with_detectors()
+
+    ww_mesh = openmc.RegularMesh()
+    n = 7
+    width = 35.0
+    ww_mesh.dimension = (n, n, n)
+    ww_mesh.lower_left = (0.0, 0.0, 0.0)
+    ww_mesh.upper_right = (width, width, width)
+
+    wwg = openmc.WeightWindowGenerator(
+        method="fw_cadis", 
+        mesh=ww_mesh, 
+        max_realizations=new_model.settings.batches
+    )
+    new_model.settings.weight_window_generators = wwg
+    new_model.settings.random_ray['volume_estimator'] = 'naive'
+
+    # Get the tallies of interest
+    target_tallies = openmc.Tallies()
+
+    for tally in list(new_model.tallies):
+        if tally.name in {"Detector 1 Tally", "Detector 2 Tally"}:
+            target_tallies.append(tally)
+    
+    # Add to WeightWindowGenerator
+    wwg.targets = target_tallies
+
+.. warning::
+    The tallies designated as FW-CADIS targets to the 
+    :class:`~openmc.WeightWindowGenerator` must be present under the
+    :class:`~openmc.model.Model.tallies` attribute of the 
+    :class:`~openmc.model.Model` as well in order to be recognized as valid 
+    local variance reduction targets. This check is performed when the 
+    :func:`openmc.model.Model.export_to_model_xml` or 
+    :func:`openmc.model.Model.export_to_xml` functions are called, meaning 
+    that the standalone :func:`openmc.Settings.export_to_xml` and 
+    :func:`openmc.Tallies.export_to_xml` methods should not be used with 
+    FW-CADIS local variance reduction.
+
+4. When running your multigroup random ray input deck, OpenMC will automatically
    run a forward solve followed by an adjoint solve, with a
    ``weight_windows.h5`` file generated at the end. The ``weight_windows.h5``
    file will contain FW-CADIS generated weight windows. This file can be used in

include/openmc/angle_energy.h

@@ -14,8 +14,22 @@ namespace openmc {
 
 class AngleEnergy {
 public:
+  //! Sample an outgoing energy and scattering cosine
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] mu Outgoing cosine with respect to current direction
+  //! \param[inout] seed Pseudorandom seed pointer
   virtual void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const = 0;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  virtual double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
   virtual ~AngleEnergy() = default;
 };
 

include/openmc/chain.h

@@ -71,6 +71,15 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   const Distribution* photon_energy_;
 };

include/openmc/distribution_angle.h

@@ -26,6 +26,12 @@ public:
   //! \return Cosine of the angle in the range [-1,1]
   double sample(double E, uint64_t* seed) const;
 
+  //! Evaluate the angular PDF at a given energy and cosine
+  //! \param[in] E Particle energy in [eV]
+  //! \param[in] mu Cosine of the scattering angle
+  //! \return Probability density for the scattering cosine
+  double evaluate(double E, double mu) const;
+
   //! Determine whether angle distribution is empty
   //! \return Whether distribution is empty
   bool empty() const { return energy_.empty(); }

include/openmc/distribution_multi.h

@@ -6,6 +6,7 @@
 #include "pugixml.hpp"
 
 #include "openmc/distribution.h"
+#include "openmc/error.h"
 #include "openmc/position.h"
 
 namespace openmc {
@@ -29,6 +30,14 @@ public:
   //! \return (sampled Direction, sample weight)
   virtual std::pair<Direction, double> sample(uint64_t* seed) const = 0;
 
+  //! Evaluate the probability density for a given direction
+  //! \param[in] u Direction on the unit sphere
+  //! \return Probability density at the given direction
+  virtual double evaluate(Direction u) const
+  {
+    fatal_error("evaluate not available for this UnitSphereDistribution type");
+  }
+
   Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
 };
 
@@ -52,6 +61,11 @@ public:
   //! \return (sampled Direction, value of the PDF at this Direction)
   std::pair<Direction, double> sample_as_bias(uint64_t* seed) const;
 
+  //! Evaluate the probability density for a given direction
+  //! \param[in] u Direction on the unit sphere
+  //! \return Probability density at the given direction
+  double evaluate(Direction u) const override;
+
   // Observing pointers
   Distribution* mu() const { return mu_.get(); }
   Distribution* phi() const { return phi_.get(); }
@@ -87,6 +101,11 @@ public:
   //! \return (sampled direction, sample weight)
   std::pair<Direction, double> sample(uint64_t* seed) const override;
 
+  //! Evaluate the probability density for a given direction
+  //! \param[in] u Direction on the unit sphere
+  //! \return Probability density at the given direction
+  double evaluate(Direction u) const override;
+
   // Set or get bias distribution
   void set_bias(std::unique_ptr<PolarAzimuthal> bias)
   {

include/openmc/lattice.h

@@ -113,6 +113,14 @@ public:
   virtual Position get_local_position(
     Position r, const array<int, 3>& i_xyz) const = 0;
 
+  //! \brief get the normal of the lattice surface crossing
+  //! \param[in] i_xyz The indices for the lattice translation.
+  //! \param[out] is_valid is the lattice translation correspond to a valid
+  //! surface. \return The surface normal corresponding to the lattice
+  //! translation.
+  virtual Direction get_normal(
+    const array<int, 3>& i_xyz, bool& is_valid) const = 0;
+
   //! \brief Check flattened lattice index.
   //! \param indx The index for a lattice tile.
   //! \return true if the given index fit within the lattice bounds.  False
@@ -223,6 +231,9 @@ public:
   Position get_local_position(
     Position r, const array<int, 3>& i_xyz) const override;
 
+  Direction get_normal(
+    const array<int, 3>& i_xyz, bool& is_valid) const override;
+
   int32_t& offset(int map, const array<int, 3>& i_xyz) override;
 
   int32_t offset(int map, int indx) const override;
@@ -268,6 +279,9 @@ public:
   Position get_local_position(
     Position r, const array<int, 3>& i_xyz) const override;
 
+  Direction get_normal(
+    const array<int, 3>& i_xyz, bool& is_valid) const override;
+
   bool is_valid_index(int indx) const override;
 
   int32_t& offset(int map, const array<int, 3>& i_xyz) override;

include/openmc/mgxs_interface.h

@@ -61,6 +61,8 @@ public:
   vector<double> energy_bin_avg_;
   vector<double> rev_energy_bins_;
   vector<vector<double>> nuc_temps_; // all available temperatures
+  vector<double>
+    default_inverse_velocity_; // approximate default inverse-velocity data
 };
 
 namespace data {

include/openmc/output.h

@@ -62,7 +62,6 @@ void write_tallies();
 void show_time(const char* label, double secs, int indent_level = 0);
 
 } // namespace openmc
-#endif // OPENMC_OUTPUT_H
 
 //////////////////////////////////////
 // Custom formatters
@@ -89,3 +88,5 @@ struct formatter<std::array<T, 2>> {
 }; // namespace fmt
 
 } // namespace fmt
+
+#endif // OPENMC_OUTPUT_H

include/openmc/particle.h

@@ -39,6 +39,8 @@ public:
 
   double speed() const;
 
+  double mass() const;
+
   //! create a secondary particle
   //
   //! stores the current phase space attributes of the particle in the

include/openmc/plot.h

@@ -17,6 +17,7 @@
 #include "openmc/particle.h"
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
+#include "openmc/ray.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -497,47 +498,6 @@ private:
   Position light_location_;
 };
 
-// Base class that implements ray tracing logic, not necessarily through
-// defined regions of the geometry but also outside of it.
-class Ray : public GeometryState {
-
-public:
-  // Initialize from location and direction
-  Ray(Position r, Direction u) { init_from_r_u(r, u); }
-
-  // Initialize from known geometry state
-  Ray(const GeometryState& p) : GeometryState(p) {}
-
-  // Called at every surface intersection within the model
-  virtual void on_intersection() = 0;
-
-  /*
-   * Traces the ray through the geometry, calling on_intersection
-   * at every surface boundary.
-   */
-  void trace();
-
-  // Stops the ray and exits tracing when called from on_intersection
-  void stop() { stop_ = true; }
-
-  // Sets the dist_ variable
-  void compute_distance();
-
-protected:
-  // Records how far the ray has traveled
-  double traversal_distance_ {0.0};
-
-private:
-  // Max intersections before we assume ray tracing is caught in an infinite
-  // loop:
-  static const int MAX_INTERSECTIONS = 1000000;
-
-  bool hit_something_ {false};
-  bool stop_ {false};
-
-  unsigned event_counter_ {0};
-};
-
 class ProjectionRay : public Ray {
 public:
   ProjectionRay(Position r, Direction u, const WireframeRayTracePlot& plot,

include/openmc/random_ray/flat_source_domain.h

@@ -40,9 +40,10 @@ public:
   void random_ray_tally();
   virtual void accumulate_iteration_flux();
   void output_to_vtk() const;
-  void convert_external_sources();
+  void convert_external_sources(bool use_adjoint_sources);
   void count_external_source_regions();
-  void set_adjoint_sources();
+  void set_fw_adjoint_sources();
+  void set_local_adjoint_sources();
   void flux_swap();
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
@@ -76,6 +77,7 @@ public:
   // Static Data members
   static bool volume_normalized_flux_tallies_;
   static bool adjoint_; // If the user wants outputs based on the adjoint flux
+  static bool fw_cadis_local_;
   static double
     diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization
                                  // for transport corrected MGXS data
@@ -84,6 +86,8 @@ public:
   static std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
     mesh_domain_map_;
 
+  static std::vector<size_t> fw_cadis_local_targets_;
+
   //----------------------------------------------------------------------------
   // Static data members
   static RandomRayVolumeEstimator volume_estimator_;

include/openmc/random_ray/random_ray_simulation.h

@@ -20,8 +20,9 @@ public:
   //----------------------------------------------------------------------------
   // Methods
   void apply_fixed_sources_and_mesh_domains();
-  void prepare_fixed_sources_adjoint();
-  void prepare_adjoint_simulation();
+  void prepare_fw_fixed_sources_adjoint();
+  void prepare_local_fixed_sources_adjoint();
+  void prepare_adjoint_simulation(bool fw_adjoint);
   void simulate();
   void output_simulation_results() const;
   void instability_check(
@@ -34,15 +35,9 @@ public:
   // Accessors
   FlatSourceDomain* domain() const { return domain_.get(); }
 
-  //----------------------------------------------------------------------------
-  // Public data members
-
-  // Flag for adjoint simulation;
-  bool adjoint_needed_;
-
 private:
   //----------------------------------------------------------------------------
-  // Private data members
+  // Data members
 
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
@@ -57,9 +52,6 @@ private:
   // Number of energy groups
   int negroups_;
 
-  // Toggle for first simulation
-  bool is_first_simulation_;
-
 }; // class RandomRaySimulation
 
 //============================================================================
@@ -67,7 +59,6 @@ private:
 //============================================================================
 
 void validate_random_ray_inputs();
-void print_adjoint_header();
 void openmc_finalize_random_ray();
 
 } // namespace openmc

include/openmc/ray.h

@@ -0,0 +1,50 @@
+#ifndef OPENMC_RAY_H
+#define OPENMC_RAY_H
+
+#include "openmc/particle_data.h"
+#include "openmc/position.h"
+
+namespace openmc {
+
+// Base class that implements ray tracing logic, not necessarily through
+// defined regions of the geometry but also outside of it.
+class Ray : public GeometryState {
+
+public:
+  // Initialize from location and direction
+  Ray(Position r, Direction u) { init_from_r_u(r, u); }
+
+  // Initialize from known geometry state
+  Ray(const GeometryState& p) : GeometryState(p) {}
+
+  // Called at every surface intersection within the model
+  virtual void on_intersection() = 0;
+
+  /*
+   * Traces the ray through the geometry, calling on_intersection
+   * at every surface boundary.
+   */
+  void trace();
+
+  // Stops the ray and exits tracing when called from on_intersection
+  void stop() { stop_ = true; }
+
+  // Sets the dist_ variable
+  void compute_distance();
+
+protected:
+  // Records how far the ray has traveled
+  double traversal_distance_ {0.0};
+
+private:
+  // Max intersections before we assume ray tracing is caught in an infinite
+  // loop:
+  static const int MAX_INTERSECTIONS = 1000000;
+
+  bool stop_ {false};
+
+  unsigned event_counter_ {0};
+};
+
+} // namespace openmc
+#endif // OPENMC_RAY_H

include/openmc/reaction_product.h

@@ -49,6 +49,21 @@ public:
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
 
+  //! Select which angle-energy distribution to sample
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Reference to the selected angle-energy distribution
+  AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const;
+
   ParticleType particle_;      //!< Particle type
   EmissionMode emission_mode_; //!< Emission mode
   double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]

include/openmc/secondary_correlated.h

@@ -41,6 +41,22 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample the outgoing energy and return the angular distribution
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Reference to the angular distribution at the sampled energy bin
+  Distribution& sample_dist(double E_in, double& E_out, uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
   // energy property
   vector<double>& energy() { return energy_; }
   const vector<double>& energy() const { return energy_; }

include/openmc/secondary_kalbach.h

@@ -32,6 +32,24 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample outgoing energy and Kalbach-Mann parameters
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] km_a Kalbach-Mann 'a' parameter
+  //! \param[out] km_r Kalbach-Mann pre-compound fraction 'r'
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample_params(double E_in, double& E_out, double& km_a, double& km_r,
+    uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   //! Outgoing energy/angle at a single incoming energy
   struct KMTable {

include/openmc/secondary_nbody.h

@@ -28,6 +28,21 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample an outgoing energy from the N-body phase space distribution
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Sampled outgoing energy in [eV]
+  double sample_energy(double E_in, uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   int n_bodies_;      //!< Number of particles distributed
   double mass_ratio_; //!< Total mass of particles [neutron mass]

include/openmc/secondary_thermal.h

@@ -6,6 +6,7 @@
 
 #include "openmc/angle_energy.h"
 #include "openmc/endf.h"
+#include "openmc/search.h"
 #include "openmc/secondary_correlated.h"
 #include "openmc/vector.h"
 
@@ -33,8 +34,20 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
+  tensor::Tensor<double> bragg_edges_; //!< Copy of Bragg edges for slicing
+  tensor::Tensor<double>
+    factors_diff_; //!< Differences over elastic scattering factors
 };
 
 //==============================================================================
@@ -56,6 +69,15 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   double debye_waller_;
 };
@@ -81,6 +103,15 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   const vector<double>& energy_;  //!< Energies at which cosines are tabulated
   tensor::Tensor<double> mu_out_; //!< Cosines for each incident energy
@@ -106,6 +137,21 @@ public:
   //! \param[inout] seed Pseudorandom number seed pointer
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+  //! Sample outgoing energy bin parameters
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] j Sampled outgoing energy bin index
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample_params(double E_in, double& E_out, int& j, uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
 
 private:
   const vector<double>& energy_; //!< Incident energies
@@ -135,6 +181,25 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample outgoing energy bin parameters
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] f Interpolation factor within sampled energy bin
+  //! \param[out] l Index of the closer incident energy
+  //! \param[out] j Sampled outgoing energy bin index
+  //! \param[inout] seed Pseudorandom seed pointer
+  void sample_params(double E_in, double& E_out, double& f, int& l, int& j,
+    uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   //! Secondary energy/angle distribution
   struct DistEnergySab {
@@ -170,6 +235,21 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Select the coherent or incoherent elastic distribution to sample
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Reference to the selected angle-energy distribution
+  const AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
 private:
   CoherentElasticAE coherent_dist_;         //!< Coherent distribution
   unique_ptr<AngleEnergy> incoherent_dist_; //!< Incoherent distribution
@@ -178,6 +258,133 @@ private:
   const Function1D& incoherent_xs_;      //!< Polymorphic ref. to incoherent XS
 };
 
+//! Internal helper for evaluating a piecewise-constant PDF on discrete points.
+//!
+//! The underlying discrete points are represented implicitly through a
+//! monotonically increasing `center(i)` function and corresponding per-point
+//! `weight(i)` values. Each point contributes a rectangular bin whose
+//! half-width is half the distance to its nearest neighboring center.
+//!
+//! \tparam CenterFn Callable returning the location of the i-th discrete value
+//! \tparam WeightFn Callable returning the weight of the i-th discrete value
+//! \param[in] n Number of discrete values
+//! \param[in] mu_0 Point at which to evaluate the PDF
+//! \param[in] a Lower bound of the domain (default: -1)
+//! \param[in] b Upper bound of the domain (default: 1)
+//! \return Probability density at mu_0
+template<typename CenterFn, typename WeightFn>
+double get_pdf_discrete_impl(std::size_t n, double mu_0, double a, double b,
+  CenterFn center, WeightFn weight)
+{
+  if (n == 0 || mu_0 < a || mu_0 > b)
+    return 0.0;
+
+  auto evaluate_bin = [&](std::size_t i) {
+    double x = center(i);
+    double left_span = (i == 0) ? 2.0 * (x - a) : x - center(i - 1);
+    double right_span = (i + 1 == n) ? 2.0 * (b - x) : center(i + 1) - x;
+    double delta = 0.5 * std::min(left_span, right_span);
+    if (delta <= 0.0)
+      return 0.0;
+
+    double left = x - delta;
+    double right = x + delta;
+    bool in_bin =
+      (mu_0 >= left) && ((i + 1 == n) ? (mu_0 <= right) : (mu_0 < right));
+    return in_bin ? weight(i) / (2.0 * delta) : 0.0;
+  };
+
+  // This is effectively a lower_bound over the sequence center(i), but the
+  // sequence is implicit rather than stored in a container, so the STL
+  // algorithms can not be used.
+  std::size_t low = 0;
+  std::size_t high = n;
+  while (low < high) {
+    std::size_t mid = low + (high - low) / 2;
+    if (center(mid) < mu_0) {
+      low = mid + 1;
+    } else {
+      high = mid;
+    }
+  }
+
+  if (low < n) {
+    double pdf = evaluate_bin(low);
+    if (pdf > 0.0)
+      return pdf;
+  }
+  if (low > 0)
+    return evaluate_bin(low - 1);
+  return 0.0;
+}
+
+//! Evaluate the PDF of a weighted discrete distribution at a given point.
+//!
+//! Given a set of discrete values mu[i] with weights w[i], this function
+//! computes the probability density at mu_0 by treating each discrete value
+//! as a rectangular bin. The bin half-width around each discrete value is
+//! half the distance to its nearest neighbor.
+//!
+//! \tparam T1 Container type for discrete cosine values (must support
+//!         operator[], size())
+//! \tparam T2 Container type for weights (must support operator[])
+//! \param[in] mu Sorted array of discrete cosine values
+//! \param[in] w Weights for each discrete value (need not be normalized)
+//! \param[in] mu_0 Point at which to evaluate the PDF
+//! \param[in] a Lower bound of the domain (default: -1)
+//! \param[in] b Upper bound of the domain (default: 1)
+//! \return Probability density at mu_0
+template<typename T1, typename T2>
+double get_pdf_discrete(
+  const T1 mu, const T2& w, double mu_0, double a = -1.0, double b = 1.0)
+{
+  // Returns the location of the discrete value for a given index
+  auto center = [&](std::size_t i) { return mu[i]; };
+  auto weight = [&](std::size_t i) { return w[i]; };
+  return get_pdf_discrete_impl(mu.size(), mu_0, a, b, center, weight);
+}
+
+//! Evaluate the PDF of a discrete distribution with uniform weights
+//!
+//! \tparam T1 Container type for discrete cosine values
+//! \param[in] mu Sorted array of discrete cosine values
+//! \param[in] mu_0 Point at which to evaluate the PDF
+//! \param[in] a Lower bound of the domain (default: -1)
+//! \param[in] b Upper bound of the domain (default: 1)
+//! \return Probability density at mu_0
+template<typename T1>
+double get_pdf_discrete(
+  const T1 mu, double mu_0, double a = -1.0, double b = 1.0)
+{
+  auto center = [&](std::size_t i) { return mu[i]; };
+  auto weight = [&](std::size_t i) { return 1.0 / mu.size(); };
+  return get_pdf_discrete_impl(mu.size(), mu_0, a, b, center, weight);
+}
+
+//! Evaluate the PDF of a uniformly weighted distribution on interpolated points
+//!
+//! \tparam T1 Container type for the lower tabulated cosine values
+//! \tparam T2 Container type for the upper tabulated cosine values
+//! \param[in] mu0 Sorted array of discrete cosine values at the lower grid
+//! \param[in] mu1 Sorted array of discrete cosine values at the upper grid
+//! \param[in] f Interpolation factor between mu0 and mu1
+//! \param[in] mu_0 Point at which to evaluate the PDF
+//! \param[in] a Lower bound of the domain (default: -1)
+//! \param[in] b Upper bound of the domain (default: 1)
+//! \return Probability density at mu_0
+template<typename T1, typename T2>
+double get_pdf_discrete_interpolated(const T1 mu0, const T2 mu1, double f,
+  double mu_0, double a = -1.0, double b = 1.0)
+{
+  if (mu0.size() != mu1.size())
+    return 0.0;
+
+  // Returns interpolated discrete value for a given index
+  auto center = [&](std::size_t i) { return mu0[i] + f * (mu1[i] - mu0[i]); };
+  auto weight = [&](std::size_t i) { return 1.0 / mu0.size(); };
+  return get_pdf_discrete_impl(mu0.size(), mu_0, a, b, center, weight);
+}
+
 } // namespace openmc
 
 #endif // OPENMC_SECONDARY_THERMAL_H

include/openmc/secondary_uncorrelated.h

@@ -32,6 +32,15 @@ public:
   void sample(
     double E_in, double& E_out, double& mu, uint64_t* seed) const override;
 
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(
+    double E_in, double mu, double& E_out, uint64_t* seed) const override;
+
   // Accessors
   AngleDistribution& angle() { return angle_; }
 

include/openmc/settings.h

@@ -77,6 +77,7 @@ extern "C" bool output_summary;      //!< write summary.h5?
 extern bool output_tallies;          //!< write tallies.out?
 extern bool particle_restart_run;    //!< particle restart run?
 extern "C" bool photon_transport;    //!< photon transport turned on?
+extern bool atomic_relaxation;       //!< atomic relaxation enabled?
 extern "C" bool reduce_tallies;      //!< reduce tallies at end of batch?
 extern bool res_scat_on;             //!< use resonance upscattering method?
 extern "C" bool restart_run;         //!< restart run?
@@ -114,6 +115,8 @@ extern std::string path_sourcepoint;      //!< path to a source file
 extern std::string path_statepoint;       //!< path to a statepoint file
 extern std::string weight_windows_file;   //!< Location of weight window file to
                                           //!< load on simulation initialization
+extern std::string properties_file;       //!< Location of properties file to
+                                          //!< load on simulation initialization
 
 // This is required because the c_str() may not be the first thing in
 // std::string. Sometimes it is, but it seems libc++ may not be like that

include/openmc/source.h

@@ -4,6 +4,7 @@
 #ifndef OPENMC_SOURCE_H
 #define OPENMC_SOURCE_H
 
+#include <atomic>
 #include <limits>
 #include <unordered_set>
 
@@ -25,15 +26,24 @@ namespace openmc {
 // source_rejection_fraction
 constexpr int EXTSRC_REJECT_THRESHOLD {10000};
 
+// Maximum number of source rejections allowed while sampling a single site
+constexpr int64_t MAX_SOURCE_REJECTIONS_PER_SAMPLE {1'000'000};
+
 //==============================================================================
 // Global variables
 //==============================================================================
 
+// Cumulative counters for source rejection diagnostics. These are atomic to
+// allow thread-safe concurrent sampling of external sources.
+extern std::atomic<int64_t> source_n_accept;
+extern std::atomic<int64_t> source_n_reject;
+
 class Source;
 
 namespace model {
 
 extern vector<unique_ptr<Source>> external_sources;
+extern vector<unique_ptr<Source>> adjoint_sources;
 
 // Probability distribution for selecting external sources
 extern DiscreteIndex external_sources_probability;
@@ -265,6 +275,9 @@ SourceSite sample_external_source(uint64_t* seed);
 
 void free_memory_source();
 
+//! Reset cumulative source rejection counters
+void reset_source_rejection_counters();
+
 } // namespace openmc
 
 #endif // OPENMC_SOURCE_H

include/openmc/tallies/tally_scoring.h

@@ -111,9 +111,9 @@ void score_meshsurface_tally(Particle& p, const vector<int>& tallies);
 //
 //! \param p The particle being tracked
 //! \param tallies A vector of the indices of the tallies to score to
-//! \param surf The surface being crossed
+//! \param normal The normal of the surface being crossed
 void score_surface_tally(
-  Particle& p, const vector<int>& tallies, const Surface& surf);
+  Particle& p, const vector<int>& tallies, const Direction& normal);
 
 //! Score the pulse-height tally
 //! This is triggered at the end of every particle history

include/openmc/thermal.h

@@ -51,7 +51,25 @@ public:
   //! \param[out] mu Outgoing scattering angle cosine
   //! \param[inout] seed Pseudorandom seed pointer
   void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out,
-    double* mu, uint64_t* seed);
+    double* mu, uint64_t* seed) const;
+
+  //! Select the elastic or inelastic distribution to sample
+  //! \param[in] micro_xs Microscopic cross sections
+  //! \param[in] E Incident neutron energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Reference to the selected angle-energy distribution
+  AngleEnergy& sample_dist(
+    const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const;
+
+  //! Sample an outgoing energy and evaluate the angular PDF
+  //! \param[in] micro_xs Microscopic cross sections
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[in] mu Scattering cosine with respect to current direction
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[inout] seed Pseudorandom seed pointer
+  //! \return Probability density for the scattering cosine
+  double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
+    double mu, double& E_out, uint64_t* seed) const;
 
 private:
   struct Reaction {

include/openmc/weight_windows.h

@@ -223,6 +223,9 @@ public:
   double threshold_ {1.0}; //<! Relative error threshold for values used to
                            // update weight windows
   double ratio_ {5.0};     //<! ratio of lower to upper weight window bounds
+
+  // Local FW-CADIS target tallies
+  std::vector<size_t> targets_;
 };
 
 //==============================================================================

openmc/cell.py

@@ -747,15 +747,6 @@ class Cell(IDManagerMixin):
             c.region = Region.from_expression(region, surfaces)
 
         # Check for other attributes
-        temperature = get_elem_list(elem, 'temperature', float)
-        if temperature is not None:
-            if len(temperature) > 1:
-                c.temperature = temperature
-            else:
-                c.temperature = temperature[0]
-        density = get_elem_list(elem, 'density', float)
-        if density is not None:
-            c.density = density if len(density) > 1 else density[0]
         v = get_text(elem, 'volume')
         if v is not None:
             c.volume = float(v)
@@ -764,6 +755,8 @@ class Cell(IDManagerMixin):
             if values is not None:
                 if key == 'rotation' and len(values) == 9:
                     values = np.array(values).reshape(3, 3)
+                elif len(values) == 1:
+                    values = values[0]
                 setattr(c, key, values)
 
         # Add this cell to appropriate universe

openmc/data/decay.py

@@ -2,7 +2,6 @@ from collections.abc import Iterable
 from functools import cached_property
 from io import StringIO
 from math import log
-import re
 from warnings import warn
 
 import numpy as np
@@ -13,7 +12,7 @@ import openmc.checkvalue as cv
 from openmc.exceptions import DataError
 from openmc.mixin import EqualityMixin
 from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
-from .data import ATOMIC_NUMBER, gnds_name
+from .data import gnds_name, zam
 from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
 
@@ -241,9 +240,7 @@ class DecayMode(EqualityMixin):
     @property
     def daughter(self):
         # Determine atomic number and mass number of parent
-        symbol, A = re.match(r'([A-Zn][a-z]*)(\d+)', self.parent).groups()
-        A = int(A)
-        Z = ATOMIC_NUMBER[symbol]
+        Z, A, _ = zam(self.parent)
 
         # Process changes
         for mode in self.modes:
@@ -253,6 +250,9 @@ class DecayMode(EqualityMixin):
                         delta_A, delta_Z = changes
                         A += delta_A
                         Z += delta_Z
+                        break
+                    else:
+                        return None
 
         return gnds_name(Z, A, self._daughter_state)
 

openmc/data/dose/mass_attenuation.py

@@ -62,7 +62,7 @@ _MUEN_TABLES = {
 def mass_energy_absorption_coefficient(
     material: str, data_source: str = "nist126"
 ) -> Tabulated1D:
-    """Return the mass energy-absorption coefficient as a function of energy.
+    r"""Return the mass energy-absorption coefficient as a function of energy.
 
     The mass energy-absorption coefficient, :math:`\mu_\text{en}/\rho`, is
     defined as the fraction of incident photon energy absorbed in a material per
@@ -108,7 +108,7 @@ _MASS_ATTENUATION: dict[int, object] = {}
 
 
 def mass_attenuation_coefficient(element):
-    """Return the photon mass attenuation coefficient as a function of energy.
+    r"""Return the photon mass attenuation coefficient as a function of energy.
 
     The mass energy-absorption coefficient, :math:`\mu_\text{en}/\rho`, is
     defined as the fraction of incident photon energy absorbed in a material per

openmc/deplete/abc.py

@@ -31,6 +31,7 @@ from .results import Results, _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \
 from .pool import deplete
 from .reaction_rates import ReactionRates
 from .transfer_rates import TransferRates, ExternalSourceRates
+from .keff_search_control import _KeffSearchControl
 
 
 __all__ = [
@@ -159,7 +160,7 @@ class TransportOperator(ABC):
             self.prev_res = prev_results
 
     @abstractmethod
-    def __call__(self, vec, source_rate):
+    def __call__(self, vec, source_rate) -> OperatorResult:
         """Runs a simulation.
 
         Parameters
@@ -201,7 +202,7 @@ class TransportOperator(ABC):
         Returns
         -------
         volume : dict of str to float
-            Volumes corresponding to materials in burn_list
+            Volumes corresponding to materials in full_burn_list
         nuc_list : list of str
             A list of all nuclide names. Used for sorting the simulation.
         burn_list : list of int
@@ -210,7 +211,7 @@ class TransportOperator(ABC):
         full_burn_list : list of int
             All burnable materials in the geometry.
         name_list : list of str
-            Material names corresponding to materials in burn_list
+            Material names corresponding to materials in full_burn_list
         """
 
     def finalize(self):
@@ -540,17 +541,15 @@ class Integrator(ABC):
         iterable of float. Alternatively, units can be specified for each step
         by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
-        Power of the reactor in [W]. A single value indicates that
-        the power is constant over all timesteps. An iterable
-        indicates potentially different power levels for each timestep.
-        For a 2D problem, the power can be given in [W/cm] as long
-        as the "volume" assigned to a depletion material is actually
-        an area in [cm^2]. Either ``power``, ``power_density``, or
+        Power of the reactor in [W]. A single value indicates that the power is
+        constant over all timesteps. An iterable indicates potentially different
+        power levels for each timestep. For a 2D problem, the power can be given
+        in [W/cm] as long as the "volume" assigned to a depletion material is
+        actually an area in [cm^2]. Either ``power``, ``power_density``, or
         ``source_rates`` must be specified.
     power_density : float or iterable of float, optional
-        Power density of the reactor in [W/gHM]. It is multiplied by
-        initial heavy metal inventory to get total power if ``power``
-        is not specified.
+        Power density of the reactor in [W/gHM]. It is multiplied by initial
+        heavy metal inventory to get total power if ``power`` is not specified.
     source_rates : float or iterable of float, optional
         Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for
         each interval in :attr:`timesteps`
@@ -562,8 +561,8 @@ class Integrator(ABC):
         and 'MWd/kg' indicates that the values are given in burnup (MW-d of
         energy deposited per kilogram of initial heavy metal).
     solver : str or callable, optional
-        If a string, must be the name of the solver responsible for
-        solving the Bateman equations.  Current options are:
+        If a string, must be the name of the solver responsible for solving the
+        Bateman equations.  Current options are:
 
             * ``cram16`` - 16th order IPF CRAM
             * ``cram48`` - 48th order IPF CRAM [default]
@@ -572,15 +571,22 @@ class Integrator(ABC):
         :attr:`solver`.
 
         .. versionadded:: 0.12
+    substeps : int, optional
+        Number of substeps per depletion interval. When greater than 1, each
+        interval is subdivided into `substeps` identical sub-intervals and LU
+        factorizations may be reused across them, improving accuracy for
+        nuclides with large decay-constant × timestep products.
+
+        .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
         previous simulation. Defaults to `False`. When `False`, the depletion
         steps provided are appended to any previous steps. If `True`, the
-        timesteps provided to the `Integrator` must exacly match any that
-        exist in the `prev_results` passed to the `Operator`. The `power`,
-        `power_density`, or `source_rates` must match as well. The
-        method of specifying `power`, `power_density`, or
-        `source_rates` should be the same as the initial run.
+        timesteps provided to the `Integrator` must exacly match any that exist
+        in the `prev_results` passed to the `Operator`. The `power`,
+        `power_density`, or `source_rates` must match as well. The method of
+        specifying `power`, `power_density`, or `source_rates` should be the
+        same as the initial run.
 
         .. versionadded:: 0.15.1
 
@@ -600,15 +606,19 @@ class Integrator(ABC):
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
         User-supplied functions are expected to have the following signature:
-        ``solver(A, n0, t) -> n1`` where
+        ``solver(A, n0, t, substeps=1) -> n1``, where
 
-            * ``A`` is a :class:`scipy.sparse.csc_array` making up the
-              depletion matrix
-            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
-              for a given material in atoms/cm3
-            * ``t`` is a float of the time step size in seconds, and
-            * ``n1`` is a :class:`numpy.ndarray` of compositions at the
-              next time step. Expected to be of the same shape as ``n0``
+            * ``A`` is a :class:`scipy.sparse.csc_array` making up the depletion
+              matrix
+            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions for
+              a given material in atoms/cm3
+            * ``t`` is a float of the time step size in seconds
+            * ``substeps`` is an optional integer number of substeps, and
+            * ``n1`` is a :class:`numpy.ndarray` of compositions at the next
+              time step. Expected to be of the same shape as ``n0``
+
+        Solvers that do not support multiple substeps should raise an exception
+        when ``substeps > 1``.
 
     transfer_rates : openmc.deplete.TransferRates
         Transfer rates for the depletion system used to model continuous
@@ -631,6 +641,7 @@ class Integrator(ABC):
             source_rates: Optional[Union[float, Sequence[float]]] = None,
             timestep_units: str = 's',
             solver: str = "cram48",
+            substeps: int = 1,
             continue_timesteps: bool = False,
         ):
         if continue_timesteps and operator.prev_res is None:
@@ -652,6 +663,8 @@ class Integrator(ABC):
         # Normalize timesteps and source rates
         seconds, source_rates = _normalize_timesteps(
             timesteps, source_rates, timestep_units, operator)
+        check_type("substeps", substeps, Integral)
+        check_greater_than("substeps", substeps, 0)
 
         if continue_timesteps:
             # Get timesteps and source rates from previous results
@@ -683,9 +696,11 @@ class Integrator(ABC):
 
         self.timesteps = np.asarray(seconds)
         self.source_rates = np.asarray(source_rates)
+        self.substeps = substeps
 
         self.transfer_rates = None
         self.external_source_rates = None
+        self._keff_search_control = None
 
         if isinstance(solver, str):
             # Delay importing of cram module, which requires this file
@@ -719,23 +734,37 @@ class Integrator(ABC):
             self._solver = func
             return
 
-        # Inspect arguments
-        if len(sig.parameters) != 3:
-            raise ValueError("Function {} does not support three arguments: "
-                             "{!s}".format(func, sig))
+        params = list(sig.parameters.values())
 
-        for ix, param in enumerate(sig.parameters.values()):
-            if param.kind in {param.KEYWORD_ONLY, param.VAR_KEYWORD}:
+        # Inspect arguments
+        if len(params) != 4:
+            raise ValueError(
+                "Function {} must support four arguments "
+                "(A, n0, t, substeps=1): {!s}"
+                .format(func, sig))
+
+        for ix, param in enumerate(params):
+            if param.kind in {param.KEYWORD_ONLY, param.VAR_KEYWORD,
+                              param.VAR_POSITIONAL}:
                 raise ValueError(
                     f"Keyword arguments like {ix} at position {param} are not allowed")
 
+        if len(params) == 4 and params[3].default != 1:
+            raise ValueError(
+                f"Fourth solver argument must default to 1, not {params[3].default}")
+
         self._solver = func
 
     def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None):
         start = time.time()
         results = deplete(
             self._solver, self.chain, n, rates, dt, i, matrix_func,
-            self.transfer_rates, self.external_source_rates)
+            self.transfer_rates, self.external_source_rates, self.substeps)
+
+        # Clip unphysical negative number densities
+        for r in results:
+            r.clip(min=0.0, out=r)
+
         return time.time() - start, results
 
     @abstractmethod
@@ -839,6 +868,37 @@ class Integrator(ABC):
         return (self.operator.prev_res[-1].time[0],
                 len(self.operator.prev_res) - 1)
 
+    def _restore_keff_search_control(self, res: StepResult):
+        """Restore keff search control from restart results."""
+        keff_search_root = res.keff_search_root
+        if keff_search_root is None:
+            raise ValueError(
+                "Cannot restore keff search control from restart "
+                "results because no stored keff_search_root is "
+                "available."
+            )
+        self._keff_search_control.function(keff_search_root)
+        return keff_search_root
+
+    def _get_bos_data(self, step_index, source_rate, bos_conc):
+        """Get beginning-of-step concentrations, rates, and control state."""
+        if step_index > 0 or self.operator.prev_res is None:
+            if self._keff_search_control is not None and source_rate != 0.0:
+                keff_search_root = self._keff_search_control.run(bos_conc)
+            else:
+                keff_search_root = None
+            bos_conc, res = self._get_bos_data_from_operator(
+                step_index, source_rate, bos_conc)
+        else:
+            bos_conc, res = self._get_bos_data_from_restart(
+                source_rate, bos_conc)
+            if self._keff_search_control is not None and source_rate != 0.0:
+                keff_search_root = self._restore_keff_search_control(self.operator.prev_res[-1])
+            else:
+                keff_search_root = None
+
+        return bos_conc, res, keff_search_root
+
     def integrate(
             self,
             final_step: bool = True,
@@ -877,11 +937,8 @@ class Integrator(ABC):
                 if output and comm.rank == 0:
                     print(f"[openmc.deplete] t={t} s, dt={dt} s, source={source_rate}")
 
-                # Solve transport equation (or obtain result from restart)
-                if i > 0 or self.operator.prev_res is None:
-                    n, res = self._get_bos_data_from_operator(i, source_rate, n)
-                else:
-                    n, res = self._get_bos_data_from_restart(source_rate, n)
+                # Get beginning-of-step data from operator or restart results
+                n, res, keff_search_root = self._get_bos_data(i, source_rate, n)
 
                 # Solve Bateman equations over time interval
                 proc_time, n_end = self(n, res.rates, dt, source_rate, i)
@@ -895,6 +952,7 @@ class Integrator(ABC):
                     self._i_res + i,
                     proc_time,
                     write_rates=write_rates,
+                    keff_search_root=keff_search_root,
                     path=path
                 )
 
@@ -908,6 +966,10 @@ class Integrator(ABC):
             # solve)
             if output and final_step and comm.rank == 0:
                 print(f"[openmc.deplete] t={t} (final operator evaluation)")
+            if self._keff_search_control is not None and source_rate != 0.0:
+                keff_search_root = self._keff_search_control.run(n)
+            else:
+                keff_search_root = None
             res_final = self.operator(n, source_rate if final_step else 0.0)
             StepResult.save(
                 self.operator,
@@ -918,6 +980,7 @@ class Integrator(ABC):
                 self._i_res + len(self),
                 proc_time,
                 write_rates=write_rates,
+                keff_search_root=keff_search_root,
                 path=path
             )
             self.operator.write_bos_data(len(self) + self._i_res)
@@ -1050,6 +1113,101 @@ class Integrator(ABC):
 
         self.transfer_rates.set_redox(material, buffer, oxidation_states, timesteps)
 
+    def add_keff_search_control(
+        self,
+        function: Callable,
+        x0: float,
+        x1: float,
+        bracket: Sequence[float],
+        **search_kwargs
+    ):
+        """Add keff search to the integrator scheme.
+
+        This method causes OpenMC to perform a keff search during depletion to
+        maintain a target keff by adjusting a model parameter through the
+        provided function.
+
+        .. important::
+            The function **must** modify the model through ``openmc.lib`` (e.g.,
+            ``openmc.lib.cells``, ``openmc.lib.materials``) and **NOT** through
+            ``openmc.Model``. The function is called within a
+            :class:`openmc.lib.TemporarySession` context where only the C API
+            (``openmc.lib``) is available for modifications.
+
+        Parameters
+        ----------
+        function : Callable
+            Function that takes a single float argument and modifies the model
+            through :mod:`openmc.lib`.
+        x0 : float
+            Initial lower bound for the keff search.
+        x1 : float
+            Initial upper bound for the keff search.
+        bracket : sequence of float
+            Bracket interval [x_min, x_max] that constrains the allowed parameter
+            values during the keff search. This is a required parameter
+            that defines the absolute bounds for the search. The bracket must contain
+            exactly 2 elements with bracket[0] < bracket[1]. These values are passed
+            directly to the ``x_min`` and ``x_max`` optional arguments in
+            :meth:`openmc.Model.keff_search`, which enforce hard limits on the
+            parameter range. If the keff search converges to a value outside this
+            bracket, it will be clamped to the nearest bracket bound with a warning.
+        **search_kwargs
+            Additional keyword arguments passed to
+            :meth:`openmc.Model.keff_search`. Common options include:
+
+            * ``target`` : float, optional
+              Target keff value to search for. Defaults to 1.0.
+            * ``k_tol`` : float, optional
+              Stopping criterion on the function value. Defaults to 1e-4.
+            * ``sigma_final`` : float, optional
+              Maximum accepted k-effective uncertainty. Defaults to 3e-4.
+            * ``maxiter`` : int, optional
+              Maximum number of iterations. Defaults to 50.
+
+            See :meth:`openmc.Model.keff_search` for a complete list of
+            available options.
+
+        Examples
+        --------
+        Add keff search that adjusts a control rod position:
+
+        >>> def adjust_rod_position(position):
+        ...     openmc.lib.cells[rod_cell.id].translation = [0, 0, position]
+        >>> integrator.add_keff_search_control(
+        ...     adjust_rod_position,
+        ...     x0=0.0,
+        ...     x1=5.0,
+        ...     bracket=[-10,10],
+        ...     target=1.0,
+        ...     k_tol=1e-4
+        ... )
+
+        Add keff search that adjusts the U235 density:
+
+        >>> def set_u235_density(u235_density):
+        ...     # Get the material from openmc.lib
+        ...     lib_mat = openmc.lib.materials[material_id]
+        ...     # Get current nuclides and densities
+        ...     nuclides = lib_mat.nuclides
+        ...     densities = lib_mat.densities
+        ...     u235_idx = nuclides.index('U235')
+        ...     densities[u235_idx] = u235_density
+        ...     lib_mat.set_densities(nuclides, densities)
+        >>> integrator.add_keff_search_control(
+        ...     set_u235_density,
+        ...     x0=5.0e-4,
+        ...     x1=1.0e-3,
+        ...     bracket=[1.0e-4, 2.0e-3],
+        ...     target=1.0
+        ... )
+
+        .. versionadded:: 0.15.4
+
+        """
+        self._keff_search_control = _KeffSearchControl(
+            self.operator, function, x0, x1, bracket, **search_kwargs)
+
 @add_params
 class SIIntegrator(Integrator):
     r"""Abstract class for the Stochastic Implicit Euler integrators
@@ -1069,17 +1227,15 @@ class SIIntegrator(Integrator):
         iterable of float. Alternatively, units can be specified for each step
         by passing an iterable of (value, unit) tuples.
     power : float or iterable of float, optional
-        Power of the reactor in [W]. A single value indicates that
-        the power is constant over all timesteps. An iterable
-        indicates potentially different power levels for each timestep.
-        For a 2D problem, the power can be given in [W/cm] as long
-        as the "volume" assigned to a depletion material is actually
-        an area in [cm^2]. Either ``power``, ``power_density``, or
+        Power of the reactor in [W]. A single value indicates that the power is
+        constant over all timesteps. An iterable indicates potentially different
+        power levels for each timestep. For a 2D problem, the power can be given
+        in [W/cm] as long as the "volume" assigned to a depletion material is
+        actually an area in [cm^2]. Either ``power``, ``power_density``, or
         ``source_rates`` must be specified.
     power_density : float or iterable of float, optional
-        Power density of the reactor in [W/gHM]. It is multiplied by
-        initial heavy metal inventory to get total power if ``power``
-        is not specified.
+        Power density of the reactor in [W/gHM]. It is multiplied by initial
+        heavy metal inventory to get total power if ``power`` is not specified.
     source_rates : float or iterable of float, optional
         Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for
         each interval in :attr:`timesteps`
@@ -1091,11 +1247,11 @@ class SIIntegrator(Integrator):
         that the values are given in burnup (MW-d of energy deposited per
         kilogram of initial heavy metal).
     n_steps : int, optional
-        Number of stochastic iterations per depletion interval.
-        Must be greater than zero. Default : 10
+        Number of stochastic iterations per depletion interval. Must be greater
+        than zero. Default : 10
     solver : str or callable, optional
-        If a string, must be the name of the solver responsible for
-        solving the Bateman equations.  Current options are:
+        If a string, must be the name of the solver responsible for solving the
+        Bateman equations.  Current options are:
 
             * ``cram16`` - 16th order IPF CRAM
             * ``cram48`` - 48th order IPF CRAM [default]
@@ -1104,16 +1260,23 @@ class SIIntegrator(Integrator):
         :attr:`solver`.
 
         .. versionadded:: 0.12
+    substeps : int, optional
+        Number of substeps per depletion interval. When greater than 1, each
+        interval is subdivided into `substeps` identical sub-intervals and LU
+        factorizations may be reused across them, improving accuracy for
+        nuclides with large decay-constant × timestep products.
+
+        .. versionadded:: 0.15.4
     continue_timesteps : bool, optional
         Whether or not to treat the current solve as a continuation of a
-        previous simulation. Defaults to `False`. If `False`, all time
-        steps and source rates will be run in an append fashion and will run
-        after whatever time steps exist, if any. If `True`, the timesteps
-        provided to the `Integrator` must match exactly those that exist
-        in the `prev_results` passed to the `Opereator`. The `power`,
-        `power_density`, or `source_rates` must match as well. The
-        method of specifying `power`, `power_density`, or
-        `source_rates` should be the same as the initial run.
+        previous simulation. Defaults to `False`. If `False`, all time steps and
+        source rates will be run in an append fashion and will run after
+        whatever time steps exist, if any. If `True`, the timesteps provided to
+        the `Integrator` must match exactly those that exist in the
+        `prev_results` passed to the `Opereator`. The `power`, `power_density`,
+        or `source_rates` must match as well. The method of specifying `power`,
+        `power_density`, or `source_rates` should be the same as the initial
+        run.
 
         .. versionadded:: 0.15.1
 
@@ -1134,15 +1297,19 @@ class SIIntegrator(Integrator):
         :math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
         size :math:`t_i`. Can be configured using the ``solver`` argument.
         User-supplied functions are expected to have the following signature:
-        ``solver(A, n0, t) -> n1`` where
+        ``solver(A, n0, t, substeps=1) -> n1``, where
 
-            * ``A`` is a :class:`scipy.sparse.csc_array` making up the
-              depletion matrix
-            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions
-              for a given material in atoms/cm3
-            * ``t`` is a float of the time step size in seconds, and
-            * ``n1`` is a :class:`numpy.ndarray` of compositions at the
-              next time step. Expected to be of the same shape as ``n0``
+            * ``A`` is a :class:`scipy.sparse.csc_array` making up the depletion
+              matrix
+            * ``n0`` is a 1-D :class:`numpy.ndarray` of initial compositions for
+              a given material in atoms/cm3
+            * ``t`` is a float of the time step size in seconds
+            * ``substeps`` is an optional integer number of substeps, and
+            * ``n1`` is a :class:`numpy.ndarray` of compositions at the next
+              time step. Expected to be of the same shape as ``n0``
+
+        Solvers that do not support multiple substeps should raise an exception
+        when ``substeps > 1``.
 
         .. versionadded:: 0.12
 
@@ -1158,13 +1325,16 @@ class SIIntegrator(Integrator):
             timestep_units: str = 's',
             n_steps: int = 10,
             solver: str = "cram48",
+            substeps: int = 1,
             continue_timesteps: bool = False,
         ):
         check_type("n_steps", n_steps, Integral)
         check_greater_than("n_steps", n_steps, 0)
         super().__init__(
             operator, timesteps, power, power_density, source_rates,
-            timestep_units=timestep_units, solver=solver, continue_timesteps=continue_timesteps)
+            timestep_units=timestep_units, solver=solver,
+            substeps=substeps,
+            continue_timesteps=continue_timesteps)
         self.n_steps = n_steps
 
     def _get_bos_data_from_operator(self, step_index, step_power, n_bos):
@@ -1294,7 +1464,7 @@ class DepSystemSolver(ABC):
     """
 
     @abstractmethod
-    def __call__(self, A, n0, dt):
+    def __call__(self, A, n0, dt, substeps=1):
         """Solve the linear system of equations for depletion
 
         Parameters
@@ -1307,6 +1477,8 @@ class DepSystemSolver(ABC):
             material or an atom density
         dt : float
             Time [s] of the specific interval to be solved
+        substeps : int, optional
+            Number of substeps to use when the solver supports substepping.
 
         Returns
         -------

openmc/deplete/chain.py

@@ -269,6 +269,7 @@ class Chain:
         self.reactions = []
         self.nuclide_dict = {}
         self._fission_yields = None
+        self._decay_matrix = None
 
     def __contains__(self, nuclide):
         return nuclide in self.nuclide_dict
@@ -412,6 +413,8 @@ class Chain:
                     type_ = ','.join(mode.modes)
                     if mode.daughter in decay_data:
                         target = mode.daughter
+                    elif 'sf' in type_:
+                        target = None
                     else:
                         print('missing {} {} {}'.format(
                             parent, type_, mode.daughter))
@@ -604,8 +607,152 @@ class Chain:
             out[nuc.name] = dict(yield_obj)
         return out
 
+    @property
+    def decay_matrix(self):
+        """Sparse CSC decay transmutation matrix.
+
+        Contains only terms from radioactive decay: diagonal loss terms
+        and off-diagonal gain terms (branching ratios, alpha/proton
+        production). Independent of reaction rates, so computed once and
+        cached.
+
+        See Also
+        --------
+        :meth:`form_rxn_matrix`, :meth:`form_matrix`
+        """
+        if self._decay_matrix is None:
+            n = len(self)
+            rows, cols, vals = [], [], []
+
+            def setval(i, j, val):
+                rows.append(i)
+                cols.append(j)
+                vals.append(val)
+
+            for i, nuc in enumerate(self.nuclides):
+                # Loss from radioactive decay
+                if nuc.half_life is not None:
+                    decay_constant = math.log(2) / nuc.half_life
+                    if decay_constant != 0.0:
+                        setval(i, i, -decay_constant)
+
+                # Gain from radioactive decay
+                if nuc.n_decay_modes != 0:
+                    for decay_type, target, branching_ratio in nuc.decay_modes:
+                        branch_val = branching_ratio * decay_constant
+
+                        # Allow for total annihilation for debug purposes
+                        if branch_val != 0.0:
+                            if target is not None and 'sf' not in decay_type:
+                                k = self.nuclide_dict[target]
+                                setval(k, i, branch_val)
+
+                            # Produce alphas and protons from decay
+                            if 'alpha' in decay_type:
+                                k = self.nuclide_dict.get('He4')
+                                if k is not None:
+                                    count = decay_type.count('alpha')
+                                    setval(k, i, count * branch_val)
+                            elif 'p' in decay_type:
+                                k = self.nuclide_dict.get('H1')
+                                if k is not None:
+                                    count = decay_type.count('p')
+                                    setval(k, i, count * branch_val)
+
+            self._decay_matrix = csc_array((vals, (rows, cols)), shape=(n, n))
+        return self._decay_matrix
+
+    def form_rxn_matrix(self, rates, fission_yields=None):
+        """Form the reaction-rate portion of the transmutation matrix.
+
+        Builds only the terms that depend on reaction rates: transmutation
+        reactions and fission product yields. Does not include radioactive
+        decay terms (see :attr:`decay_matrix`).
+
+        Parameters
+        ----------
+        rates : numpy.ndarray
+            2D array indexed by (nuclide, reaction)
+        fission_yields : dict, optional
+            Option to use a custom set of fission yields. Expected
+            to be of the form ``{parent : {product : f_yield}}``
+            with string nuclide names for ``parent`` and ``product``,
+            and ``f_yield`` as the respective fission yield
+
+        Returns
+        -------
+        scipy.sparse.csc_array
+            Sparse matrix representing reaction-rate terms.
+
+        See Also
+        --------
+        :attr:`decay_matrix`, :meth:`form_matrix`
+        """
+        reactions = set()
+        n = len(self)
+
+        # Accumulate indices/values and then create the matrix at the end to
+        # avoid expensive index checks scipy otherwise does.
+        rows, cols, vals = [], [], []
+
+        def setval(i, j, val):
+            rows.append(i)
+            cols.append(j)
+            vals.append(val)
+
+        if fission_yields is None:
+            fission_yields = self.get_default_fission_yields()
+
+        # Save local variables to avoid attribute lookups in loop
+        index_nuc = rates.index_nuc
+        index_rx = rates.index_rx
+
+        for i, nuc in enumerate(self.nuclides):
+            if nuc.name not in index_nuc:
+                continue
+
+            nuc_ind = index_nuc[nuc.name]
+            nuc_rates = rates[nuc_ind, :]
+
+            for r_type, target, _, br in nuc.reactions:
+                r_id = index_rx[r_type]
+                path_rate = nuc_rates[r_id]
+
+                # Loss term -- make sure we only count loss once for
+                # reactions with branching ratios
+                if r_type not in reactions:
+                    reactions.add(r_type)
+                    if path_rate != 0.0:
+                        setval(i, i, -path_rate)
+
+                # Gain term; allow for total annihilation for debug purposes
+                if r_type != 'fission':
+                    if target is not None and path_rate != 0.0:
+                        k = self.nuclide_dict[target]
+                        setval(k, i, path_rate * br)
+
+                    # Determine light nuclide production, e.g., (n,d) should
+                    # produce H2
+                    if path_rate != 0.0:
+                        light_nucs = REACTIONS[r_type].secondaries
+                        for light_nuc in light_nucs:
+                            k = self.nuclide_dict.get(light_nuc)
+                            if k is not None:
+                                setval(k, i, path_rate * br)
+
+                else:
+                    for product, y in fission_yields[nuc.name].items():
+                        yield_val = y * path_rate
+                        if yield_val != 0.0:
+                            k = self.nuclide_dict[product]
+                            setval(k, i, yield_val)
+
+            reactions.clear()
+
+        return csc_array((vals, (rows, cols)), shape=(n, n))
+
     def form_matrix(self, rates, fission_yields=None):
-        """Forms depletion matrix.
+        """Form the full transmutation matrix (decay + reactions).
 
         Parameters
         ----------
@@ -624,96 +771,10 @@ class Chain:
 
         See Also
         --------
+        :attr:`decay_matrix`, :meth:`form_rxn_matrix`,
         :meth:`get_default_fission_yields`
         """
-        reactions = set()
-
-        n = len(self)
-
-        # we accumulate indices and value entries for everything and create the matrix
-        # in one step at the end to avoid expensive index checks scipy otherwise does.
-        rows, cols, vals = [], [], []
-        def setval(i, j, val):
-            rows.append(i)
-            cols.append(j)
-            vals.append(val)
-
-        if fission_yields is None:
-            fission_yields = self.get_default_fission_yields()
-
-        for i, nuc in enumerate(self.nuclides):
-            # Loss from radioactive decay
-            if nuc.half_life is not None:
-                decay_constant = math.log(2) / nuc.half_life
-                if decay_constant != 0.0:
-                    setval(i, i, -decay_constant)
-
-            # Gain from radioactive decay
-            if nuc.n_decay_modes != 0:
-                for decay_type, target, branching_ratio in nuc.decay_modes:
-                    branch_val = branching_ratio * decay_constant
-
-                    # Allow for total annihilation for debug purposes
-                    if branch_val != 0.0:
-                        if target is not None:
-                            k = self.nuclide_dict[target]
-                            setval(k, i, branch_val)
-
-                        # Produce alphas and protons from decay
-                        if 'alpha' in decay_type:
-                            k = self.nuclide_dict.get('He4')
-                            if k is not None:
-                                count = decay_type.count('alpha')
-                                setval(k, i, count * branch_val)
-                        elif 'p' in decay_type:
-                            k = self.nuclide_dict.get('H1')
-                            if k is not None:
-                                count = decay_type.count('p')
-                                setval(k, i, count * branch_val)
-
-            if nuc.name in rates.index_nuc:
-                # Extract all reactions for this nuclide in this cell
-                nuc_ind = rates.index_nuc[nuc.name]
-                nuc_rates = rates[nuc_ind, :]
-
-                for r_type, target, _, br in nuc.reactions:
-                    # Extract reaction index, and then final reaction rate
-                    r_id = rates.index_rx[r_type]
-                    path_rate = nuc_rates[r_id]
-
-                    # Loss term -- make sure we only count loss once for
-                    # reactions with branching ratios
-                    if r_type not in reactions:
-                        reactions.add(r_type)
-                        if path_rate != 0.0:
-                            setval(i, i, -path_rate)
-
-                    # Gain term; allow for total annihilation for debug purposes
-                    if r_type != 'fission':
-                        if target is not None and path_rate != 0.0:
-                            k = self.nuclide_dict[target]
-                            setval(k, i, path_rate * br)
-
-                        # Determine light nuclide production, e.g., (n,d) should
-                        # produce H2
-                        light_nucs = REACTIONS[r_type].secondaries
-                        for light_nuc in light_nucs:
-                            k = self.nuclide_dict.get(light_nuc)
-                            if k is not None:
-                                setval(k, i, path_rate * br)
-
-                    else:
-                        for product, y in fission_yields[nuc.name].items():
-                            yield_val = y * path_rate
-                            if yield_val != 0.0:
-                                k = self.nuclide_dict[product]
-                                setval(k, i, yield_val)
-
-                # Clear set of reactions
-                reactions.clear()
-
-        # Return CSC representation instead of DOK
-        return csc_array((vals, (rows, cols)), shape=(n, n))
+        return self.decay_matrix + self.form_rxn_matrix(rates, fission_yields)
 
     def add_redox_term(self, matrix, buffer, oxidation_states):
         r"""Adds a redox term to the depletion matrix from data contained in
@@ -807,7 +868,7 @@ class Chain:
         # Use DOK as intermediate representation
         n = len(self)
         matrix = dok_array((n, n))
-        
+
         check_type("mats", mats, (tuple, str))
         if not isinstance(mats, str):
             check_type("mats", mats, tuple, str)
@@ -816,8 +877,8 @@ class Chain:
         else:
             mat = mats
             dest_mat = None
-        
-        # Build transfer term 
+
+        # Build transfer term
         components = tr_rates.get_components(mat, current_timestep, dest_mat)
 
         for i, nuc in enumerate(self.nuclides):
@@ -829,7 +890,7 @@ class Chain:
             else:
                 continue
             matrix[i, i] = sum(tr_rates.get_external_rate(mat, key, current_timestep, dest_mat))
-                
+
         # Return CSC instead of DOK
         return matrix.tocsc()
 
@@ -1363,6 +1424,7 @@ def _get_chain(
 
 def _invalidate_chain_cache(chain):
     """Invalidate the cache for a specific Chain (when it is modifed)."""
+    chain._decay_matrix = None
     if hasattr(chain, '_xml_path'):
         # Remove all entries with the same path as self._xml_path
         for key in list(_CHAIN_CACHE.keys()):

openmc/deplete/coupled_operator.py

@@ -405,7 +405,7 @@ class CoupledOperator(OpenMCOperator):
 
         self.materials.export_to_xml(nuclides_to_ignore=self._decay_nucs)
 
-    def __call__(self, vec, source_rate):
+    def __call__(self, vec, source_rate) -> OperatorResult:
         """Runs a simulation.
 
         Simulation will abort under the following circumstances:

openmc/deplete/cram.py

@@ -3,12 +3,13 @@
 Implements two different forms of CRAM for use in openmc.deplete.
 """
 
+from functools import partial
 import numbers
 
 import numpy as np
-import scipy.sparse.linalg as sla
+from scipy.sparse.linalg import spsolve, splu
 
-from openmc.checkvalue import check_type, check_length
+from openmc.checkvalue import check_type, check_length, check_greater_than
 from .abc import DepSystemSolver
 from .._sparse_compat import csc_array, eye_array
 
@@ -24,6 +25,12 @@ class IPFCramSolver(DepSystemSolver):
     Chebyshev Rational Approximation Method and Application to Burnup Equations
     <https://doi.org/10.13182/NSE15-26>`_," Nucl. Sci. Eng., 182:3, 297-318.
 
+    When `substeps` > 1, the time interval is split into `substeps` identical
+    sub-intervals and LU factorizations are reused across them, as described
+    in: A. Isotalo and M. Pusa, "`Improving the Accuracy of the Chebyshev
+    Rational Approximation Method Using Substeps
+    <https://doi.org/10.13182/NSE15-67>`_," Nucl. Sci. Eng., 183:1, 65-77.
+
     Parameters
     ----------
     alpha : numpy.ndarray
@@ -55,7 +62,7 @@ class IPFCramSolver(DepSystemSolver):
         self.theta = theta
         self.alpha0 = alpha0
 
-    def __call__(self, A, n0, dt):
+    def __call__(self, A, n0, dt, substeps=1):
         """Solve depletion equations using IPF CRAM
 
         Parameters
@@ -68,6 +75,8 @@ class IPFCramSolver(DepSystemSolver):
             material or an atom density
         dt : float
             Time [s] of the specific interval to be solved
+        substeps : int, optional
+            Number of substeps per depletion interval.
 
         Returns
         -------
@@ -75,12 +84,25 @@ class IPFCramSolver(DepSystemSolver):
             Final compositions after ``dt``
 
         """
-        A = dt * csc_array(A, dtype=np.float64)
-        y = n0.copy()
+        check_type("substeps", substeps, numbers.Integral)
+        check_greater_than("substeps", substeps, 0)
+
+        step_dt = dt if substeps == 1 else dt / substeps
+        A = step_dt * csc_array(A, dtype=np.float64)
         ident = eye_array(A.shape[0], format='csc')
-        for alpha, theta in zip(self.alpha, self.theta):
-            y += 2*np.real(alpha*sla.spsolve(A - theta*ident, y))
-        return y * self.alpha0
+
+        if substeps == 1:
+            solvers = [partial(spsolve, A - theta * ident) for theta in self.theta]
+        else:
+            # Pre-compute LU factorizations and reuse them across substeps.
+            solvers = [splu(A - theta * ident).solve for theta in self.theta]
+
+        y = n0.copy()
+        for _ in range(substeps):
+            for alpha, solve in zip(self.alpha, solvers):
+                y += 2 * np.real(alpha * solve(y))
+            y *= self.alpha0
+        return y
 
 
 # Coefficients for IPF Cram 16

openmc/deplete/independent_operator.py

@@ -384,7 +384,7 @@ class IndependentOperator(OpenMCOperator):
         # Return number density vector
         return super().initial_condition(self.materials)
 
-    def __call__(self, vec, source_rate):
+    def __call__(self, vec, source_rate) -> OperatorResult:
         """Obtain the reaction rates
 
         Parameters

openmc/deplete/keff_search_control.py

@@ -0,0 +1,128 @@
+from typing import Callable
+from warnings import warn
+
+import openmc.lib
+
+
+class _KeffSearchControl:
+    """Controller for keff search during depletion calculations.
+
+    This class performs keff searches to maintain a target keff by adjusting a
+    model parameter through a provided function.
+
+    Parameters
+    ----------
+    operator : openmc.deplete.Operator
+        Depletion operator instance
+    function : Callable
+        Function that modifies the model based on a parameter value
+    x0 : float
+        Initial lower bound for the keff search
+    x1 : float
+        Initial upper bound for the keff search
+    bracket : list[float]
+        Absolute bracketing interval lower and upper. If the keff search
+        solution lies off these limits the closest limit will be set as new
+        result.
+    **search_kwargs : dict, optional
+        Additional keyword arguments to pass to :meth:`openmc.Model.keff_search`
+
+    """
+    def __init__(self, operator, function: Callable, x0: float, x1: float, bracket: list[float], **search_kwargs):
+        if len(bracket) != 2:
+            raise ValueError(f"bracket must have exactly 2 elements, got {len(bracket)}")
+        if bracket[0] >= bracket[1]:
+            raise ValueError(f"bracket[0] must be < bracket[1], got {bracket}")
+        self.x0 = x0
+        self.x1 = x1
+        self.operator = operator
+        self.function = function
+        self.search_kwargs = search_kwargs
+        self.search_kwargs['x_min'] = bracket[0]
+        self.search_kwargs['x_max'] = bracket[1]
+
+    def run(self, x):
+        """Perform keff search and update the atom density vector.
+
+        Parameters
+        ----------
+        x : list of numpy.ndarray
+            Current atom density vector (atoms per material)
+
+        Returns
+        -------
+        root : float
+            Parameter value that achieves target keff
+        """
+        root = self._search_for_keff()
+        self._update_vec(x)
+        return root
+
+    def _search_for_keff(self) -> float:
+        """Perform the keff search using the model's keff_search method.
+
+        Returns
+        -------
+        float
+            Parameter value that achieves target keff
+
+        Raises
+        ------
+        ValueError
+            If the keff search fails to converge
+        """
+        with openmc.lib.TemporarySession(self.operator.model):
+            # Only pass the first 3 required args plus explicitly provided kwargs
+            result = self.operator.model.keff_search(
+                self.function, self.x0, self.x1, **self.search_kwargs
+            )
+        if not result.converged:
+            raise ValueError(
+                f"Search for keff failed to converge. "
+                f"Termination reason: {result.flag}"
+            )
+
+        root = result.root
+
+        # Check if root is outside the bracket bounds and give a warning
+        if root < self.search_kwargs['x_min']:
+            warn(f"keff search result ({root:.6f}) is below the lower bracket "
+                 f"bound ({self.search_kwargs['x_min']:.6f}).", UserWarning)
+        elif root > self.search_kwargs['x_max']:
+            warn(f"keff search result ({root:.6f}) is above the upper bracket "
+                 f"bound ({self.search_kwargs['x_max']:.6f}).", UserWarning)
+
+        # Restore the number of initial batches
+        openmc.lib.settings.set_batches(self.operator.model.settings.batches)
+
+        return root
+
+    def _update_vec(self, x):
+        """Update the atom density vector from openmc.lib.materials and AtomNumber object.
+
+        The depletion vector ``x`` is rank-local, matching the materials owned
+        by ``self.operator.number`` on the current MPI rank. We therefore only
+        update entries for locally owned materials using the compositions
+        currently stored in ``openmc.lib.materials``.
+
+        Parameters
+        ----------
+        x : list of numpy.ndarray
+            Atom density vector to update (atoms per material)
+
+        """
+        number = self.operator.number
+
+        for mat_idx, mat in enumerate(number.materials):
+            lib_material = openmc.lib.materials[int(mat)]
+            nuclides = lib_material.nuclides
+            densities = 1e24 * lib_material.densities
+            volume = number.get_mat_volume(mat)
+
+            for nuc_idx, nuc in enumerate(number.burnable_nuclides):
+                if nuc in nuclides:
+                    lib_nuc_idx = nuclides.index(nuc)
+                    atom_density = densities[lib_nuc_idx]
+                else:
+                    atom_density = number.get_atom_density(mat, nuc)
+                x[mat_idx][nuc_idx] = atom_density * volume

openmc/deplete/microxs.py

@@ -36,7 +36,8 @@ DomainTypes: TypeAlias = Union[
     Sequence[openmc.Cell],
     Sequence[openmc.Universe],
     openmc.MeshBase,
-    openmc.Filter
+    openmc.Filter,
+    Sequence[openmc.Filter]
 ]
 
 
@@ -69,8 +70,12 @@ def get_microxs_and_flux(
     ----------
     model : openmc.Model
         OpenMC model object. Must contain geometry, materials, and settings.
-    domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter
+    domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter, or list of openmc.Filter
         Domains in which to tally reaction rates, or a spatial tally filter.
+        A list of filters can be provided to create one set of tallies per
+        filter (e.g., one :class:`~openmc.MeshMaterialFilter` per mesh) that
+        are all evaluated in a single transport solve. Results are
+        concatenated across all filters in order.
     nuclides : list of str
         Nuclides to get cross sections for. If not specified, all burnable
         nuclides from the depletion chain file are used.
@@ -142,26 +147,24 @@ def get_microxs_and_flux(
     else:
         energy_filter = openmc.EnergyFilter(energies)
 
+    # Build list of domain filters
     if isinstance(domains, openmc.Filter):
-        domain_filter = domains
+        domain_filters = [domains]
     elif isinstance(domains, openmc.MeshBase):
-        domain_filter = openmc.MeshFilter(domains)
+        domain_filters = [openmc.MeshFilter(domains)]
+    elif isinstance(domains, Sequence) and len(domains) > 0 and \
+            isinstance(domains[0], openmc.Filter):
+        domain_filters = list(domains)
     elif isinstance(domains[0], openmc.Material):
-        domain_filter = openmc.MaterialFilter(domains)
+        domain_filters = [openmc.MaterialFilter(domains)]
     elif isinstance(domains[0], openmc.Cell):
-        domain_filter = openmc.CellFilter(domains)
+        domain_filters = [openmc.CellFilter(domains)]
     elif isinstance(domains[0], openmc.Universe):
-        domain_filter = openmc.UniverseFilter(domains)
+        domain_filters = [openmc.UniverseFilter(domains)]
     else:
         raise ValueError(f"Unsupported domain type: {type(domains[0])}")
 
-    flux_tally = openmc.Tally(name='MicroXS flux')
-    flux_tally.filters = [domain_filter, energy_filter]
-    flux_tally.scores = ['flux']
-    model.tallies = [flux_tally]
-
-    # Prepare reaction-rate tally for 'direct' or subset for 'flux' with opts
-    rr_tally = None
+    # Prepare reaction-rate nuclides/reactions
     rr_nuclides: list[str] = []
     rr_reactions: list[str] = []
     if reaction_rate_mode == 'direct':
@@ -177,20 +180,33 @@ def get_microxs_and_flux(
         if rr_reactions:
             rr_reactions = [r for r in rr_reactions if r in set(reactions)]
 
-    # Only construct tally if both lists are non-empty
-    if rr_nuclides and rr_reactions:
-        rr_tally = openmc.Tally(name='MicroXS RR')
-        # Use 1-group energy filter for RR in flux mode
-        if reaction_rate_mode == 'flux':
-            rr_energy_filter = openmc.EnergyFilter(
-                [energy_filter.values[0], energy_filter.values[-1]])
-        else:
-            rr_energy_filter = energy_filter
-        rr_tally.filters = [domain_filter, rr_energy_filter]
-        rr_tally.nuclides = rr_nuclides
-        rr_tally.multiply_density = False
-        rr_tally.scores = rr_reactions
-        model.tallies.append(rr_tally)
+    # Use 1-group energy filter for RR in flux mode
+    has_rr = bool(rr_nuclides and rr_reactions)
+    if has_rr and reaction_rate_mode == 'flux':
+        rr_energy_filter = openmc.EnergyFilter(
+            [energy_filter.values[0], energy_filter.values[-1]])
+    else:
+        rr_energy_filter = energy_filter
+
+    # Create one flux tally (and optionally one RR tally) per domain filter.
+    flux_tallies = []
+    rr_tallies = []
+    model.tallies = []
+    for i, domain_filter in enumerate(domain_filters):
+        flux_tally = openmc.Tally(name=f'MicroXS flux {i}')
+        flux_tally.filters = [domain_filter, energy_filter]
+        flux_tally.scores = ['flux']
+        model.tallies.append(flux_tally)
+        flux_tallies.append(flux_tally)
+
+        if has_rr:
+            rr_tally = openmc.Tally(name=f'MicroXS RR {i}')
+            rr_tally.filters = [domain_filter, rr_energy_filter]
+            rr_tally.nuclides = rr_nuclides
+            rr_tally.multiply_density = False
+            rr_tally.scores = rr_reactions
+            model.tallies.append(rr_tally)
+            rr_tallies.append(rr_tally)
 
     if openmc.lib.is_initialized:
         openmc.lib.finalize()
@@ -227,40 +243,41 @@ def get_microxs_and_flux(
 
         # Read in tally results (on all ranks)
         with StatePoint(statepoint_path) as sp:
-            if rr_tally is not None:
-                rr_tally = sp.tallies[rr_tally.id]
-                rr_tally._read_results()
-            flux_tally = sp.tallies[flux_tally.id]
-            flux_tally._read_results()
+            for i in range(len(flux_tallies)):
+                flux_tallies[i] = sp.tallies[flux_tallies[i].id]
+                flux_tallies[i]._read_results()
+                if rr_tallies:
+                    rr_tallies[i] = sp.tallies[rr_tallies[i].id]
+                    rr_tallies[i]._read_results()
 
-    # Get flux values and make energy groups last dimension
-    flux = flux_tally.get_reshaped_data()  # (domains, groups, 1, 1)
-    flux = np.moveaxis(flux, 1, -1)  # (domains, 1, 1, groups)
+    # Concatenate results across all domain filters
+    fluxes = []
+    all_flux_arrays = []
+    for flux_tally in flux_tallies:
+        # Get flux values and make energy groups last dimension
+        flux = flux_tally.get_reshaped_data()  # (domains, groups, 1, 1)
+        flux = np.moveaxis(flux, 1, -1)  # (domains, 1, 1, groups)
+        all_flux_arrays.append(flux)
+        fluxes.extend(flux.squeeze((1, 2)))
 
-    # Create list where each item corresponds to one domain
-    fluxes = list(flux.squeeze((1, 2)))
+    # If we built reaction-rate tallies, compute microscopic cross sections
+    if rr_tallies:
+        direct_micros = []
+        for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies):
+            flux = flux_arr
+            # Get reaction rates and make energy groups last dimension
+            reaction_rates = rr_tally.get_reshaped_data()  # (domains, groups, nuclides, reactions)
+            reaction_rates = np.moveaxis(reaction_rates, 1, -1)  # (domains, nuclides, reactions, groups)
 
-    # If we built a reaction-rate tally, compute microscopic cross sections
-    if rr_tally is not None:
-        # Get reaction rates
-        reaction_rates = rr_tally.get_reshaped_data()  # (domains, groups, nuclides, reactions)
+            # If RR is 1-group, sum flux over groups
+            if reaction_rate_mode == "flux":
+                flux = flux.sum(axis=-1, keepdims=True)
 
-        # Make energy groups last dimension
-        reaction_rates = np.moveaxis(reaction_rates, 1, -1)  # (domains, nuclides, reactions, groups)
-
-        # If RR is 1-group, sum flux over groups
-        if reaction_rate_mode == "flux":
-            flux = flux.sum(axis=-1, keepdims=True)  # (domains, 1, 1, 1)
-
-        # Divide RR by flux to get microscopic cross sections. The indexing
-        # ensures that only non-zero flux values are used, and broadcasting is
-        # applied to align the shapes of reaction_rates and flux for division.
-        xs = np.zeros_like(reaction_rates)  # (domains, nuclides, reactions, groups)
-        d, _, _, g = np.nonzero(flux)
-        xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
-
-        # Create lists where each item corresponds to one domain
-        direct_micros = [MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs]
+            xs = np.zeros_like(reaction_rates)
+            d, _, _, g = np.nonzero(flux)
+            xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
+            direct_micros.extend(
+                MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs)
 
     # If using flux mode, compute flux-collapsed microscopic XS
     if reaction_rate_mode == 'flux':
@@ -273,9 +290,9 @@ def get_microxs_and_flux(
         ) for flux_i in fluxes]
 
     # Decide which micros to use and merge if needed
-    if reaction_rate_mode == 'flux' and rr_tally is not None:
+    if reaction_rate_mode == 'flux' and rr_tallies:
         micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)]
-    elif rr_tally is not None:
+    elif rr_tallies:
         micros = direct_micros
     else:
         micros = flux_micros

openmc/deplete/pool.py

@@ -42,14 +42,15 @@ def _distribute(items):
         j += chunk_size
 
 def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
-            transfer_rates=None, external_source_rates=None, *matrix_args):
+            transfer_rates=None, external_source_rates=None, substeps=1,
+            *matrix_args):
     """Deplete materials using given reaction rates for a specified time
 
     Parameters
     ----------
     func : callable
         Function to use to get new compositions. Expected to have the signature
-        ``func(A, n0, t) -> n1``
+        ``func(A, n0, t, substeps=1) -> n1``.
     chain : openmc.deplete.Chain
         Depletion chain
     n : list of numpy.ndarray
@@ -74,6 +75,8 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
         External source rates for continuous removal/feed.
 
         .. versionadded:: 0.15.3
+    substeps : int, optional
+        Number of substeps to pass to solvers that support substepping.
     matrix_args: Any, optional
         Additional arguments passed to matrix_func
 
@@ -164,7 +167,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
 
                 # Concatenate vectors of nuclides in one
                 n_multi = np.concatenate(n)
-                n_result = func(matrix, n_multi, dt)
+                n_result = func(matrix, n_multi, dt, substeps)
 
                 # Split back the nuclide vector result into the original form
                 n_result = np.split(n_result, np.cumsum([len(i) for i in n])[:-1])
@@ -198,7 +201,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
                 matrix.resize(matrix.shape[1], matrix.shape[1])
                 n[i] = np.append(n[i], 1.0)
 
-    inputs = zip(matrices, n, repeat(dt))
+    inputs = zip(matrices, n, repeat(dt), repeat(substeps))
 
     if USE_MULTIPROCESSING:
         with Pool(NUM_PROCESSES) as pool:

openmc/deplete/r2s.py

@@ -13,11 +13,11 @@ from .results import Results
 from ..checkvalue import PathLike
 from ..mpi import comm
 from openmc.lib import TemporarySession
-from openmc.utility_funcs import change_directory
 
 
 def get_activation_materials(
-    model: openmc.Model, mmv: openmc.MeshMaterialVolumes
+    model: openmc.Model,
+    mmv_list: list[openmc.MeshMaterialVolumes]
 ) -> openmc.Materials:
     """Get a list of activation materials for each mesh element/material.
 
@@ -31,35 +31,35 @@ def get_activation_materials(
     ----------
     model : openmc.Model
         The full model containing the geometry and materials.
-    mmv : openmc.MeshMaterialVolumes
-        The mesh material volumes object containing the materials and their
-        volumes for each mesh element.
+    mmv_list : list of openmc.MeshMaterialVolumes
+        List of mesh material volumes objects, one per mesh, containing the
+        materials and their volumes for each mesh element.
 
     Returns
     -------
     openmc.Materials
         A list of materials, each corresponding to a unique mesh element and
-        material combination.
+        material combination across all meshes.
 
     """
-    # Get the material ID, volume, and element index for each element-material
-    # combination
-    mat_ids = mmv._materials[mmv._materials > -1]
-    volumes = mmv._volumes[mmv._materials > -1]
-    elems, _ = np.where(mmv._materials > -1)
-
     # Get all materials in the model
     material_dict = model._get_all_materials()
 
     # Create a new activation material for each element-material combination
+    # across all meshes
     materials = openmc.Materials()
-    for elem, mat_id, vol in zip(elems, mat_ids, volumes):
-        mat = material_dict[mat_id]
-        new_mat = mat.clone()
-        new_mat.depletable = True
-        new_mat.name = f'Element {elem}, Material {mat_id}'
-        new_mat.volume = vol
-        materials.append(new_mat)
+    for mesh_idx, mmv in enumerate(mmv_list):
+        mat_ids = mmv._materials[mmv._materials > -1]
+        volumes = mmv._volumes[mmv._materials > -1]
+        elems, _ = np.where(mmv._materials > -1)
+
+        for elem, mat_id, vol in zip(elems, mat_ids, volumes):
+            mat = material_dict[mat_id]
+            new_mat = mat.clone()
+            new_mat.depletable = True
+            new_mat.name = f'Mesh {mesh_idx}, Element {elem}, Material {mat_id}'
+            new_mat.volume = vol
+            materials.append(new_mat)
 
     return materials
 
@@ -70,7 +70,9 @@ class R2SManager:
     This class is responsible for managing the materials and sources needed for
     mesh-based or cell-based R2S calculations. It provides methods to get
     activation materials and decay photon sources based on the mesh/cells and
-    materials in the OpenMC model.
+    materials in the OpenMC model. Multiple meshes can be specified as domains,
+    in which case each element--material combination of each mesh is treated as
+    an activation region (meshes are assumed to be non-overlapping).
 
     This class supports the use of a different models for the neutron and photon
     transport calculation. However, for cell-based calculations, it assumes that
@@ -83,17 +85,20 @@ class R2SManager:
     ----------
     neutron_model : openmc.Model
         The OpenMC model to use for neutron transport.
-    domains : openmc.MeshBase or Sequence[openmc.Cell]
-        The mesh or a sequence of cells that represent the spatial units over
-        which the R2S calculation will be performed.
+    domains : openmc.MeshBase or Sequence[openmc.MeshBase] or Sequence[openmc.Cell]
+        The mesh(es) or a sequence of cells that represent the spatial units
+        over which the R2S calculation will be performed. When a single
+        :class:`~openmc.MeshBase` or a sequence of meshes is given, each
+        element--material combination across all meshes is treated as an
+        activation region.
     photon_model : openmc.Model, optional
         The OpenMC model to use for photon transport calculations. If None, a
         shallow copy of the neutron_model will be created and used.
 
     Attributes
     ----------
-    domains : openmc.MeshBase or Sequence[openmc.Cell]
-        The mesh or a sequence of cells that represent the spatial units over
+    domains : list of openmc.MeshBase or Sequence[openmc.Cell]
+        The meshes or a sequence of cells that represent the spatial units over
         which the R2S calculation will be performed.
     neutron_model : openmc.Model
         The OpenMC model used for neutron transport.
@@ -101,7 +106,7 @@ class R2SManager:
         The OpenMC model used for photon transport calculations.
     method : {'mesh-based', 'cell-based'}
         Indicates whether the R2S calculation uses mesh elements ('mesh-based')
-        as the spatial discetization or a list of a cells ('cell-based').
+        as the spatial discretization or a list of cells ('cell-based').
     results : dict
         A dictionary that stores results from the R2S calculation.
 
@@ -109,7 +114,7 @@ class R2SManager:
     def __init__(
         self,
         neutron_model: openmc.Model,
-        domains: openmc.MeshBase | Sequence[openmc.Cell],
+        domains: openmc.MeshBase | Sequence[openmc.MeshBase] | Sequence[openmc.Cell],
         photon_model: openmc.Model | None = None,
     ):
         self.neutron_model = neutron_model
@@ -126,9 +131,14 @@ class R2SManager:
             self.photon_model = photon_model
         if isinstance(domains, openmc.MeshBase):
             self.method = 'mesh-based'
+            self.domains = [domains]
+        elif isinstance(domains, Sequence) and len(domains) > 0 and \
+                isinstance(domains[0], openmc.MeshBase):
+            self.method = 'mesh-based'
+            self.domains = list(domains)
         else:
             self.method = 'cell-based'
-        self.domains = domains
+            self.domains = list(domains)
         self.results = {}
 
     def run(
@@ -243,11 +253,13 @@ class R2SManager:
     ):
         """Run the neutron transport step.
 
-        This step computes the material volume fractions on the mesh, creates a
-        mesh-material filter, and retrieves the fluxes and microscopic cross
-        sections for each mesh/material combination. This step will populate the
-        'fluxes' and 'micros' keys in the results dictionary. For a mesh-based
-        calculation, it will also populate the 'mesh_material_volumes' key.
+        This step computes the material volume fractions on each mesh, creates
+        mesh-material filters, and retrieves the fluxes and microscopic cross
+        sections for each mesh/material combination via a single transport
+        solve. This step will populate the 'fluxes' and 'micros' keys in the
+        results dictionary. For a mesh-based calculation, it will also populate
+        the 'mesh_material_volumes' key (a list of
+        :class:`~openmc.MeshMaterialVolumes`, one per mesh).
 
         Parameters
         ----------
@@ -266,19 +278,28 @@ class R2SManager:
         output_dir.mkdir(parents=True, exist_ok=True)
 
         if self.method == 'mesh-based':
-            # Compute material volume fractions on the mesh
+            # Compute material volume fractions on each mesh
             if mat_vol_kwargs is None:
                 mat_vol_kwargs = {}
             mat_vol_kwargs.setdefault('bounding_boxes', True)
-            self.results['mesh_material_volumes'] = mmv = comm.bcast(
-                self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
 
-            # Save results to file
-            if comm.rank == 0:
-                mmv.save(output_dir / 'mesh_material_volumes.npz')
+            mmv_list = []
+            domain_filters = []
+            for i, mesh in enumerate(self.domains):
+                mmv = comm.bcast(
+                    mesh.material_volumes(self.neutron_model, **mat_vol_kwargs))
+                mmv_list.append(mmv)
 
-            # Create mesh-material filter based on what combos were found
-            domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv)
+                # Save results to file
+                if comm.rank == 0:
+                    mmv.save(output_dir / f'mesh_material_volumes_{i}.npz')
+
+                # Create mesh-material filter for this mesh
+                domain_filters.append(
+                    openmc.MeshMaterialFilter.from_volumes(mesh, mmv))
+
+            self.results['mesh_material_volumes'] = mmv_list
+            domains = domain_filters
         else:
             domains: Sequence[openmc.Cell] = self.domains
 
@@ -357,8 +378,9 @@ class R2SManager:
 
         if self.method == 'mesh-based':
             # Get unique material for each (mesh, material) combination
-            mmv = self.results['mesh_material_volumes']
-            self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv)
+            mmv_list = self.results['mesh_material_volumes']
+            self.results['activation_materials'] = get_activation_materials(
+                self.neutron_model, mmv_list)
         else:
             # Create unique material for each cell
             activation_mats = openmc.Materials()
@@ -468,12 +490,20 @@ class R2SManager:
         # photon model if it is different from the neutron model to account for
         # potential material changes
         if self.method == 'mesh-based' and different_photon_model:
-            self.results['mesh_material_volumes_photon'] = photon_mmv = comm.bcast(
-                self.domains.material_volumes(self.photon_model, **mat_vol_kwargs))
+            if mat_vol_kwargs is None:
+                mat_vol_kwargs = {}
+            photon_mmv_list = []
+            for i, mesh in enumerate(self.domains):
+                photon_mmv = comm.bcast(
+                    mesh.material_volumes(self.photon_model, **mat_vol_kwargs))
+                photon_mmv_list.append(photon_mmv)
 
-            # Save photon MMV results to file
-            if comm.rank == 0:
-                photon_mmv.save(output_dir / 'mesh_material_volumes.npz')
+                # Save photon MMV results to file
+                if comm.rank == 0:
+                    photon_mmv.save(
+                        output_dir / f'mesh_material_volumes_{i}.npz')
+
+            self.results['mesh_material_volumes_photon'] = photon_mmv_list
 
         if comm.rank == 0:
             tally_ids = [tally.id for tally in self.photon_model.tallies]
@@ -543,7 +573,7 @@ class R2SManager:
     ) -> list[openmc.IndependentSource]:
         """Create decay photon source for a mesh-based calculation.
 
-        For each mesh element-material combination, an
+        For each mesh element-material combination across all meshes, an
         :class:`~openmc.IndependentSource` is created with a
         :class:`~openmc.stats.Box` spatial distribution based on the bounding
         box of the material within the mesh element. A material constraint is
@@ -575,52 +605,56 @@ class R2SManager:
         index_mat = 0
 
         # Get various results from previous steps
-        mat_vols = self.results['mesh_material_volumes']
+        mmv_list = self.results['mesh_material_volumes']
         materials = self.results['activation_materials']
         results = self.results['depletion_results']
-        photon_mat_vols = self.results.get('mesh_material_volumes_photon')
+        photon_mmv_list = self.results.get('mesh_material_volumes_photon')
 
-        # Total number of mesh elements
-        n_elements = mat_vols.num_elements
+        for mesh_idx, mat_vols in enumerate(mmv_list):
+            photon_mat_vols = photon_mmv_list[mesh_idx] \
+                if photon_mmv_list is not None else None
 
-        for index_elem in range(n_elements):
-            # Determine which materials exist in the photon model for this element
-            if photon_mat_vols is not None:
-                photon_materials = {
-                    mat_id
-                    for mat_id, _ in photon_mat_vols.by_element(index_elem)
-                    if mat_id is not None
-                }
+            # Total number of mesh elements for this mesh
+            n_elements = mat_vols.num_elements
 
-            for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
-                # Skip void volume
-                if mat_id is None:
-                    continue
+            for index_elem in range(n_elements):
+                # Determine which materials exist in the photon model for this element
+                if photon_mat_vols is not None:
+                    photon_materials = {
+                        mat_id
+                        for mat_id, _ in photon_mat_vols.by_element(index_elem)
+                        if mat_id is not None
+                    }
 
-                # Skip if this material doesn't exist in photon model
-                if photon_mat_vols is not None and mat_id not in photon_materials:
+                for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
+                    # Skip void volume
+                    if mat_id is None:
+                        continue
+
+                    # Skip if this material doesn't exist in photon model
+                    if photon_mat_vols is not None and mat_id not in photon_materials:
+                        index_mat += 1
+                        continue
+
+                    # Get activated material composition
+                    original_mat = materials[index_mat]
+                    activated_mat = results[time_index].get_material(str(original_mat.id))
+
+                    # Create decay photon source
+                    energy = activated_mat.get_decay_photon_energy()
+                    if energy is not None:
+                        strength = energy.integral()
+                        space = openmc.stats.Box(*bbox)
+                        sources.append(openmc.IndependentSource(
+                            space=space,
+                            energy=energy,
+                            particle='photon',
+                            strength=strength,
+                            constraints={'domains': [mat_dict[mat_id]]}
+                        ))
+
+                    # Increment index of activated material
                     index_mat += 1
-                    continue
-
-                # Get activated material composition
-                original_mat = materials[index_mat]
-                activated_mat = results[time_index].get_material(str(original_mat.id))
-
-                # Create decay photon source
-                energy = activated_mat.get_decay_photon_energy()
-                if energy is not None:
-                    strength = energy.integral()
-                    space = openmc.stats.Box(*bbox)
-                    sources.append(openmc.IndependentSource(
-                        space=space,
-                        energy=energy,
-                        particle='photon',
-                        strength=strength,
-                        constraints={'domains': [mat_dict[mat_id]]}
-                    ))
-
-                # Increment index of activated material
-                index_mat += 1
 
         return sources
 
@@ -638,10 +672,13 @@ class R2SManager:
         # Load neutron transport results
         neutron_dir = path / 'neutron_transport'
         if self.method == 'mesh-based':
-            mmv_file = neutron_dir / 'mesh_material_volumes.npz'
-            if mmv_file.exists():
-                self.results['mesh_material_volumes'] = \
-                    openmc.MeshMaterialVolumes.from_npz(mmv_file)
+            mmv_files = sorted(neutron_dir.glob('mesh_material_volumes*.npz'),
+                               key=lambda p: int(p.stem.split('_')[-1])
+                               if p.stem[-1].isdigit() else 0)
+            if mmv_files:
+                self.results['mesh_material_volumes'] = [
+                    openmc.MeshMaterialVolumes.from_npz(f) for f in mmv_files
+                ]
         fluxes_file = neutron_dir / 'fluxes.npy'
         if fluxes_file.exists():
             self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True))
@@ -665,10 +702,15 @@ class R2SManager:
 
         # Load photon mesh material volumes if they exist (for mesh-based calculations)
         if self.method == 'mesh-based':
-            photon_mmv_file = photon_dir / 'mesh_material_volumes.npz'
-            if photon_mmv_file.exists():
-                self.results['mesh_material_volumes_photon'] = \
-                    openmc.MeshMaterialVolumes.from_npz(photon_mmv_file)
+            photon_mmv_files = sorted(
+                photon_dir.glob('mesh_material_volumes*.npz'),
+                key=lambda p: int(p.stem.split('_')[-1])
+                if p.stem[-1].isdigit() else 0)
+            if photon_mmv_files:
+                self.results['mesh_material_volumes_photon'] = [
+                    openmc.MeshMaterialVolumes.from_npz(f)
+                    for f in photon_mmv_files
+                ]
 
         # Load tally IDs from JSON file
         tally_ids_path = photon_dir / 'tally_ids.json'

openmc/deplete/results.py

@@ -113,7 +113,7 @@ class Results(list):
         ----------
         mat : openmc.Material, str
             Material object or material id to evaluate
-        units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'}
+        units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'}
             Specifies the type of activity to return, options include total
             activity [Bq], specific [Bq/g, Bq/kg] or volumetric activity [Bq/cm3].
         by_nuclide : bool
@@ -231,7 +231,7 @@ class Results(list):
         ----------
         mat : openmc.Material, str
             Material object or material id to evaluate.
-        units : {'W', 'W/g', 'W/kg', 'W/cm3'}
+        units : {'W', 'W/g', 'W/kg', 'W/cm3', 'W/m3'}
             Specifies the units of decay heat to return. Options include total
             heat [W], specific [W/g, W/kg] or volumetric heat [W/cm3].
         by_nuclide : bool

openmc/deplete/stepresult.py

@@ -12,11 +12,12 @@ import h5py
 import numpy as np
 
 import openmc
-from openmc.mpi import comm, MPI
 from openmc.checkvalue import PathLike
+from openmc.mpi import MPI, comm
+
 from .reaction_rates import ReactionRates
 
-VERSION_RESULTS = (1, 2)
+VERSION_RESULTS = (1, 3)
 
 
 __all__ = ["StepResult"]
@@ -57,6 +58,8 @@ class StepResult:
     proc_time : int
         Average time spent depleting a material across all
         materials and processes
+    keff_search_root : float
+        The root returned by the keff search control.
 
     """
     def __init__(self):
@@ -73,6 +76,7 @@ class StepResult:
         self.name_list = None
 
         self.data = None
+        self.keff_search_root = None
 
     def __repr__(self):
         t = self.time[0]
@@ -153,14 +157,14 @@ class StepResult:
         full_burn_list : list of str
             List of all burnable material IDs
         name_list : list of str, optional
-            Material names corresponding to materials in burn_list
+            Material names corresponding to materials in full_burn_list
 
         """
         self.volume = copy.deepcopy(volume)
         self.index_nuc = {nuc: i for i, nuc in enumerate(nuc_list)}
         self.index_mat = {mat: i for i, mat in enumerate(burn_list)}
         self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
-        self.mat_to_name = dict(zip(burn_list, name_list)) if name_list is not None else {}
+        self.mat_to_name = dict(zip(full_burn_list, name_list)) if name_list is not None else {}
 
         # Create storage array
         self.data = np.zeros((self.n_mat, self.n_nuc))
@@ -196,15 +200,15 @@ class StepResult:
         new.rates = self.rates[ranges]
         return new
 
-    def get_material(self, mat_id):
+    def get_material(self, mat_id: str | int) -> openmc.Material:
         """Return material object for given depleted composition
 
         .. versionadded:: 0.13.2
 
         Parameters
         ----------
-        mat_id : str
-            Material ID as a string
+        mat_id : str or int
+            Material ID as a string or integer
 
         Returns
         -------
@@ -217,6 +221,9 @@ class StepResult:
             If specified material ID is not found in the StepResult
 
         """
+        # Coerce to str since internal dictionaries use str keys
+        mat_id = str(mat_id)
+
         with warnings.catch_warnings():
             warnings.simplefilter('ignore', openmc.IDWarning)
             material = openmc.Material(material_id=int(mat_id))
@@ -364,6 +371,10 @@ class StepResult:
             "depletion time", (1,), maxshape=(None,),
             dtype="float64")
 
+        handle.create_dataset(
+            "keff_search_root", (1,), maxshape=(None,),
+            dtype="float64")
+
     def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
         """Converts results object into an hdf5 object.
 
@@ -396,6 +407,7 @@ class StepResult:
         time_dset = handle["/time"]
         source_rate_dset = handle["/source_rate"]
         proc_time_dset = handle["/depletion time"]
+        keff_search_root_dset = handle["/keff_search_root"]
 
         # Get number of results stored
         number_shape = list(number_dset.shape)
@@ -429,6 +441,10 @@ class StepResult:
             proc_shape[0] = new_shape
             proc_time_dset.resize(proc_shape)
 
+            keff_search_root_shape = list(keff_search_root_dset.shape)
+            keff_search_root_shape[0] = new_shape
+            keff_search_root_dset.resize(keff_search_root_shape)
+
         # If nothing to write, just return
         if len(self.index_mat) == 0:
             return
@@ -448,6 +464,7 @@ class StepResult:
                 proc_time_dset[index] = (
                     self.proc_time / (comm.size * self.n_hdf5_mats)
                 )
+            keff_search_root_dset[index] = self.keff_search_root
 
     @classmethod
     def from_hdf5(cls, handle, step):
@@ -496,6 +513,10 @@ class StepResult:
             if step < proc_time_dset.shape[0]:
                 results.proc_time = proc_time_dset[step]
 
+        if "keff_search_root" in handle:
+            keff_search_root_dset = handle["/keff_search_root"]
+            results.keff_search_root = keff_search_root_dset[step]
+
         if results.proc_time is None:
             results.proc_time = np.array([np.nan])
 
@@ -550,6 +571,7 @@ class StepResult:
         step_ind,
         proc_time=None,
         write_rates: bool = False,
+        keff_search_root=None,
         path: PathLike = "depletion_results.h5"
     ):
         """Creates and writes depletion results to disk
@@ -574,6 +596,8 @@ class StepResult:
             processes.
         write_rates : bool, optional
             Whether reaction rates should be written to the results file.
+        keff_search_root : float
+            The root returned by the keff search control.
         path : PathLike
             Path to file to write. Defaults to 'depletion_results.h5'.
 
@@ -601,6 +625,7 @@ class StepResult:
         results.proc_time = proc_time
         if results.proc_time is not None:
             results.proc_time = comm.reduce(proc_time, op=MPI.SUM)
+        results.keff_search_root = keff_search_root
 
         if not Path(path).is_file():
             Path(path).parent.mkdir(parents=True, exist_ok=True)

openmc/examples.py

@@ -1310,3 +1310,312 @@ def random_ray_three_region_cube() -> openmc.Model:
     model.tallies = tallies
 
     return model
+
+def random_ray_three_region_cube_with_detectors() -> openmc.Model:
+    """Create a three region cube model with two external tally regions.
+
+    This is an adaptation of the simple monoenergetic problem of a cube with 
+    three concentric cubic regions. The innermost region is near void (with 
+    Sigma_t around 10^-5) and contains an external isotropic source term, the 
+    middle region is a mild scatterer (with Sigma_t around 10^-3), and the 
+    outer region of the cube is a scatterer and absorber (with Sigma_t around 
+    1).
+
+    Two cubic "detector" regions are found outside this geometry, one along the 
+    y-axis near z=0, and the other in the upper right corner of the system. 
+    The size of each detector is scaled to be equal to that of the source 
+    region. The model returned by this function contains cell tallies on each 
+    detector.
+
+    Returns
+    -------
+    model : openmc.Model
+        A three region cube model
+
+    """
+
+    model = openmc.Model()
+
+    ###########################################################################
+    # Helper function creates a 3 region cube with different fills in each region
+    def fill_cube(N, n_1, n_2, fill_1, fill_2, fill_3):
+        cube = [[[0 for _ in range(N)] for _ in range(N)] for _ in range(N)]
+        for i in range(N):
+            for j in range(N):
+                for k in range(N):
+                    if i < n_1 and j >= (N-n_1) and k < n_1:
+                        cube[i][j][k] = fill_1
+                    elif i < n_2 and j >= (N-n_2) and k < n_2:
+                        cube[i][j][k] = fill_2
+                    else:
+                        cube[i][j][k] = fill_3
+        return cube
+
+    ###########################################################################
+    # Create multigroup data
+
+    # Instantiate the energy group data
+    ebins = [1e-5, 20.0e6]
+    groups = openmc.mgxs.EnergyGroups(group_edges=ebins)
+
+    cavity_sigma_a = 4.0e-5
+    cavity_sigma_s = 3.0e-3
+    cavity_mat_data = openmc.XSdata('cavity', groups)
+    cavity_mat_data.order = 0
+    cavity_mat_data.set_total([cavity_sigma_a + cavity_sigma_s])
+    cavity_mat_data.set_absorption([cavity_sigma_a])
+    cavity_mat_data.set_scatter_matrix(
+        np.rollaxis(np.array([[[cavity_sigma_s]]]), 0, 3))
+
+    absorber_sigma_a = 0.50
+    absorber_sigma_s = 0.50
+    absorber_mat_data = openmc.XSdata('absorber', groups)
+    absorber_mat_data.order = 0
+    absorber_mat_data.set_total([absorber_sigma_a + absorber_sigma_s])
+    absorber_mat_data.set_absorption([absorber_sigma_a])
+    absorber_mat_data.set_scatter_matrix(
+        np.rollaxis(np.array([[[absorber_sigma_s]]]), 0, 3))
+
+    multiplier = 0.01
+    source_sigma_a = cavity_sigma_a * multiplier
+    source_sigma_s = cavity_sigma_s * multiplier
+    source_mat_data = openmc.XSdata('source', groups)
+    source_mat_data.order = 0
+    source_mat_data.set_total([source_sigma_a + source_sigma_s])
+    source_mat_data.set_absorption([source_sigma_a])
+    source_mat_data.set_scatter_matrix(
+        np.rollaxis(np.array([[[source_sigma_s]]]), 0, 3))
+
+    mg_cross_sections_file = openmc.MGXSLibrary(groups)
+    mg_cross_sections_file.add_xsdatas(
+        [source_mat_data, cavity_mat_data, absorber_mat_data])
+    mg_cross_sections_file.export_to_hdf5()
+
+    ###########################################################################
+    # Create materials for the problem
+
+    # Instantiate some Macroscopic Data
+    source_data = openmc.Macroscopic('source')
+    cavity_data = openmc.Macroscopic('cavity')
+    absorber_data = openmc.Macroscopic('absorber')
+
+    # Instantiate some Materials and register the appropriate Macroscopic objects
+    source_mat = openmc.Material(name='source')
+    source_mat.set_density('macro', 1.0)
+    source_mat.add_macroscopic(source_data)
+
+    cavity_mat = openmc.Material(name='cavity')
+    cavity_mat.set_density('macro', 1.0)
+    cavity_mat.add_macroscopic(cavity_data)
+
+    absorber_mat = openmc.Material(name='absorber')
+    absorber_mat.set_density('macro', 1.0)
+    absorber_mat.add_macroscopic(absorber_data)
+
+    # Instantiate a Materials collection
+    materials_file = openmc.Materials([source_mat, cavity_mat, absorber_mat])
+    materials_file.cross_sections = "mgxs.h5"
+
+    ###########################################################################
+    # Define problem geometry
+
+    source_cell = openmc.Cell(fill=source_mat, name='infinite source region')
+    cavity_cell = openmc.Cell(fill=cavity_mat, name='cube cavity region')
+    absorber_cell = openmc.Cell(
+        fill=absorber_mat, name='absorber region')
+
+    source_universe = openmc.Universe(name='source universe')
+    source_universe.add_cells([source_cell])
+
+    cavity_universe = openmc.Universe()
+    cavity_universe.add_cells([cavity_cell])
+
+    absorber_universe = openmc.Universe()
+    absorber_universe.add_cells([absorber_cell])
+
+    absorber_width = 30.0
+    n_base = 6
+
+    # This variable can be increased above 1 to refine the FSR mesh resolution further
+    refinement_level = 2
+
+    n = n_base * refinement_level
+    pitch = absorber_width / n
+
+    pattern = fill_cube(n, 1*refinement_level, 5*refinement_level,
+                        source_universe, cavity_universe, absorber_universe)
+
+    lattice = openmc.RectLattice()
+    lattice.lower_left = [0.0, 0.0, 0.0]
+    lattice.pitch = [pitch, pitch, pitch]
+    lattice.universes = pattern
+
+    lattice_cell = openmc.Cell(fill=lattice)
+
+    lattice_uni = openmc.Universe()
+    lattice_uni.add_cells([lattice_cell])
+
+    x_low = openmc.XPlane(x0=0.0, boundary_type='reflective')
+    x_high = openmc.XPlane(x0=absorber_width)
+
+    y_low = openmc.YPlane(y0=0.0, boundary_type='reflective')
+    y_high = openmc.YPlane(y0=absorber_width)
+
+    z_low = openmc.ZPlane(z0=0.0, boundary_type='reflective')
+    z_high = openmc.ZPlane(z0=absorber_width)
+
+    cube_domain = openmc.Cell(fill=lattice_uni, region=+x_low & -
+                              x_high & +y_low & -y_high & +z_low & -z_high, name='full domain')
+
+    detect_width = absorber_width / n_base
+    outer_width = absorber_width + detect_width
+
+    x_outer = openmc.XPlane(x0=outer_width, boundary_type='vacuum')
+    y_outer = openmc.YPlane(y0=outer_width, boundary_type='vacuum')
+    z_outer = openmc.ZPlane(z0=outer_width, boundary_type='vacuum')
+
+    detector1_right = openmc.XPlane(x0=detect_width)
+    detector1_top = openmc.ZPlane(z0=detect_width)
+
+    detector1_region = (
+        +x_low & -detector1_right &
+        +y_high & -y_outer &
+        +z_low & -detector1_top
+    )
+    detector1 = openmc.Cell(
+        name='detector 1',
+        fill=absorber_mat,
+        region=detector1_region
+    )
+
+    detector2_region = (
+        +x_high & -x_outer &
+        +y_high & -y_outer &
+        +z_high & -z_outer
+    )
+    detector2 = openmc.Cell(
+        name='detector 2',
+        fill=absorber_mat,
+        region=detector2_region
+    )
+    
+    external_x = (
+        +x_high & +y_low & +z_low & -x_outer & 
+        ((-y_outer & -z_high) | (-y_high & +z_high & -z_outer))
+    )
+    external_y = (
+        +y_high & -y_outer & 
+        (
+            (+detector1_right & -x_high & +z_low & -z_outer) | 
+            (-detector1_right & +x_low & +detector1_top & -z_outer) | 
+            (+x_high & -x_outer & +z_low & -z_high)
+        )
+    )
+    external_z = (
+        +x_low & +y_low & +z_high & -z_outer & 
+        ((-y_outer & -x_high) | (-y_high & +x_high & -x_outer))
+    )
+    external_cell = openmc.Cell(fill=cavity_mat, 
+                                region=(external_x | external_y | external_z), 
+                                name='outside cube')
+
+    root = openmc.Universe(
+        name='root universe', 
+        cells=[cube_domain, detector1, detector2, external_cell]
+    )
+
+    # Create a geometry with the two cells and export to XML
+    geometry = openmc.Geometry(root)
+
+    ###########################################################################
+    # Define problem settings
+
+    # Instantiate a Settings object, set all runtime parameters, and export to XML
+    settings = openmc.Settings()
+    settings.energy_mode = "multi-group"
+    settings.inactive = 5
+    settings.batches = 10
+    settings.particles = 500
+    settings.run_mode = 'fixed source'
+
+    # Create an initial uniform spatial source for ray integration
+    lower_left_ray = [0.0, 0.0, 0.0]
+    upper_right_ray = [outer_width, outer_width, outer_width]
+    uniform_dist_ray = openmc.stats.Box(
+        lower_left_ray, upper_right_ray, only_fissionable=False)
+    rr_source = openmc.IndependentSource(space=uniform_dist_ray)
+
+    settings.random_ray['distance_active'] = 800.0
+    settings.random_ray['distance_inactive'] = 100.0
+    settings.random_ray['ray_source'] = rr_source
+    settings.random_ray['volume_normalized_flux_tallies'] = True
+
+    # Create a rectilinear source region mesh
+    sr_mesh = openmc.RegularMesh()
+    sr_mesh.dimension = (14, 14, 14)
+    sr_mesh.lower_left = (0.0, 0.0, 0.0)
+    sr_mesh.upper_right = (outer_width, outer_width, outer_width)
+    settings.random_ray['source_region_meshes'] = [(sr_mesh, [root])]
+
+    # Create the neutron source in the bottom right of the moderator
+    # Good - fast group appears largest (besides most thermal)
+    strengths = [1.0]
+    midpoints = [100.0]
+    energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
+
+    source = openmc.IndependentSource(energy=energy_distribution, constraints={
+                                      'domains': [source_universe]}, strength=3.14)
+
+    settings.source = [source]
+
+    ###########################################################################
+    # Define tallies
+
+    estimator = 'tracklength'
+
+    detector1_filter = openmc.CellFilter(detector1)
+    detector1_tally = openmc.Tally(name="Detector 1 Tally")
+    detector1_tally.filters = [detector1_filter]
+    detector1_tally.scores = ['flux']
+    detector1_tally.estimator = estimator
+
+    detector2_filter = openmc.CellFilter(detector2)
+    detector2_tally = openmc.Tally(name="Detector 2 Tally")
+    detector2_tally.filters = [detector2_filter]
+    detector2_tally.scores = ['flux']
+    detector2_tally.estimator = estimator
+
+    absorber_filter = openmc.MaterialFilter(absorber_mat)
+    absorber_tally = openmc.Tally(name="Absorber Tally")
+    absorber_tally.filters = [absorber_filter]
+    absorber_tally.scores = ['flux']
+    absorber_tally.estimator = estimator
+
+    cavity_filter = openmc.MaterialFilter(cavity_mat)
+    cavity_tally = openmc.Tally(name="Cavity Tally")
+    cavity_tally.filters = [cavity_filter]
+    cavity_tally.scores = ['flux']
+    cavity_tally.estimator = estimator
+
+    source_filter = openmc.MaterialFilter(source_mat)
+    source_tally = openmc.Tally(name="Source Tally")
+    source_tally.filters = [source_filter]
+    source_tally.scores = ['flux']
+    source_tally.estimator = estimator
+
+    # Instantiate a Tallies collection and export to XML
+    tallies = openmc.Tallies([detector1_tally, 
+                              detector2_tally, 
+                              absorber_tally,
+                              cavity_tally,
+                              source_tally])
+
+    ###########################################################################
+    # Assmble Model
+
+    model.geometry = geometry
+    model.materials = materials_file
+    model.settings = settings
+    model.tallies = tallies
+
+    return model

openmc/executor.py

@@ -103,6 +103,7 @@ def _run(args, output, cwd):
         # If OpenMC is finished, break loop
         line = p.stdout.readline()
         if not line and p.poll() is not None:
+            p.stdout.close()
             break
 
         lines.append(line)

openmc/lib/core.py

@@ -33,7 +33,6 @@ class _SourceSite(Structure):
                 ('parent_id', c_int64),
                 ('progeny_id', c_int64)]
 
-
 # Define input type for numpy arrays that will be passed into C++ functions
 # Must be an int or double array, with single dimension that is contiguous
 _array_1d_int = np.ctypeslib.ndpointer(dtype=np.int32, ndim=1,
@@ -494,8 +493,9 @@ def run_random_ray(output=True):
 
 def sample_external_source(
         n_samples: int = 1000,
-        prn_seed: int | None = None
-) -> openmc.ParticleList:
+        prn_seed: int | None = None,
+        as_array: bool = False
+) -> openmc.ParticleList | np.ndarray:
     """Sample external source and return source particles.
 
     .. versionadded:: 0.13.1
@@ -507,11 +507,20 @@ def sample_external_source(
     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`.  The array has fields ``'r'`` (float64,
+        shape 3), ``'u'`` (float64, shape 3), ``'E'`` (float64), ``'time'``
+        (float64), ``'wgt'`` (float64), ``'delayed_group'`` (int32),
+        ``'surf_id'`` (int32), and ``'particle'`` (int32).  This avoids the
+        overhead of constructing individual :class:`~openmc.SourceParticle`
+        objects and is substantially faster for large sample counts.
 
     Returns
     -------
-    openmc.ParticleList
-        List of sampled source particles
+    openmc.ParticleList or numpy.ndarray
+        List of sampled source particles, or a structured array when
+        *as_array* is True.
 
     """
     if n_samples <= 0:
@@ -519,18 +528,28 @@ def sample_external_source(
     if prn_seed is None:
         prn_seed = getrandbits(63)
 
-    # Call into C API to sample source
-    sites_array = (_SourceSite * n_samples)()
-    _dll.openmc_sample_external_source(c_size_t(n_samples), c_uint64(prn_seed), sites_array)
+    # Pre-allocate output array and sample all particles in a single C call
+    result = np.empty(n_samples, dtype=_SourceSite)
+    sites_array = (_SourceSite * n_samples).from_buffer(result)
+    _dll.openmc_sample_external_source(
+        c_size_t(n_samples),
+        c_uint64(prn_seed),
+        sites_array,
+    )
 
-    # Convert to list of SourceParticle and return
-    return openmc.ParticleList([openmc.SourceParticle(
-            r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt,
-            delayed_group=site.delayed_group, surf_id=site.surf_id,
-            particle=openmc.ParticleType(site.particle)
+    if as_array:
+        return result
+
+    particles = [
+        openmc.SourceParticle(
+            r=site.r, u=site.u, E=site.E, time=site.time,
+            wgt=site.wgt, delayed_group=site.delayed_group,
+            surf_id=site.surf_id,
+            particle=openmc.ParticleType(site.particle),
         )
         for site in sites_array
-    ])
+    ]
+    return openmc.ParticleList(particles)
 
 
 def simulation_init():
@@ -674,8 +693,8 @@ class TemporarySession:
         self.model = model
 
         # Determine MPI intercommunicator
-        self.init_kwargs.setdefault('intracomm', comm)
-        self.comm = self.init_kwargs['intracomm']
+        self.comm = self.init_kwargs.get('intracomm') or comm
+        self.init_kwargs['intracomm'] = self.comm
 
     def __enter__(self):
         """Initialize the OpenMC library in a temporary directory."""

openmc/material.py

@@ -349,7 +349,7 @@ class Material(IDManagerMixin):
         clip_tolerance : float
             Maximum fraction of :math:`\sum_i x_i p_i` for discrete distributions
             that will be discarded.
-        units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'}
+        units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'}
             Specifies the units on the integral of the distribution.
         volume : float, optional
             Volume of the material. If not passed, defaults to using the
@@ -367,7 +367,7 @@ class Material(IDManagerMixin):
             the total intensity of the photon source in the requested units.
 
         """
-        cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'})
+        cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3'})
 
         if exclude_nuclides is not None and include_nuclides is not None:
             raise ValueError("Cannot specify both exclude_nuclides and include_nuclides")
@@ -378,6 +378,8 @@ class Material(IDManagerMixin):
                 raise ValueError("volume must be specified if units='Bq'")
         elif units == 'Bq/cm3':
             multiplier = 1
+        elif units == 'Bq/m3':
+            multiplier = 1e6
         elif units == 'Bq/g':
             multiplier = 1.0 / self.get_mass_density()
         elif units == 'Bq/kg':
@@ -1383,16 +1385,16 @@ class Material(IDManagerMixin):
 
     def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
                      volume: float | None = None) -> dict[str, float] | float:
-        """Returns the activity of the material or of each nuclide within.
+        """Return the activity of the material or each nuclide within.
 
         .. versionadded:: 0.13.1
 
         Parameters
         ----------
-        units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Ci', 'Ci/m3'}
+        units : {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3', 'Ci', 'Ci/m3'}
             Specifies the type of activity to return, options include total
             activity [Bq,Ci], specific [Bq/g, Bq/kg] or volumetric activity
-            [Bq/cm3,Ci/m3]. Default is volumetric activity [Bq/cm3].
+            [Bq/cm3, Bq/m3, Ci/m3]. Default is volumetric activity [Bq/cm3].
         by_nuclide : bool
             Specifies if the activity should be returned for the material as a
             whole or per nuclide. Default is False.
@@ -1410,7 +1412,7 @@ class Material(IDManagerMixin):
             of the material is returned as a float.
         """
 
-        cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Ci', 'Ci/m3'})
+        cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3', 'Bq/m3', 'Ci', 'Ci/m3'})
         cv.check_type('by_nuclide', by_nuclide, bool)
 
         if volume is None:
@@ -1420,6 +1422,8 @@ class Material(IDManagerMixin):
             multiplier = volume
         elif units == 'Bq/cm3':
             multiplier = 1
+        elif units == 'Bq/m3':
+            multiplier = 1e6
         elif units == 'Bq/g':
             multiplier = 1.0 / self.get_mass_density()
         elif units == 'Bq/kg':
@@ -1438,16 +1442,15 @@ class Material(IDManagerMixin):
 
     def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False,
                        volume: float | None = None) -> dict[str, float] | float:
-        """Returns the decay heat of the material or for each nuclide in the
-        material in units of [W], [W/g], [W/kg] or [W/cm3].
+        """Return the decay heat of the material or each nuclide within.
 
         .. versionadded:: 0.13.3
 
         Parameters
         ----------
-        units : {'W', 'W/g', 'W/kg', 'W/cm3'}
+        units : {'W', 'W/g', 'W/kg', 'W/cm3', 'W/m3'}
             Specifies the units of decay heat to return. Options include total
-            heat [W], specific [W/g, W/kg] or volumetric heat [W/cm3].
+            heat [W], specific [W/g, W/kg] or volumetric heat [W/cm3, W/m3].
             Default is total heat [W].
         by_nuclide : bool
             Specifies if the decay heat should be returned for the material as a
@@ -1466,13 +1469,15 @@ class Material(IDManagerMixin):
             of the material is returned as a float.
         """
 
-        cv.check_value('units', units, {'W', 'W/g', 'W/kg', 'W/cm3'})
+        cv.check_value('units', units, {'W', 'W/g', 'W/kg', 'W/cm3', 'W/m3'})
         cv.check_type('by_nuclide', by_nuclide, bool)
 
         if units == 'W':
             multiplier = volume if volume is not None else self.volume
         elif units == 'W/cm3':
             multiplier = 1
+        elif units == 'W/m3':
+            multiplier = 1e6
         elif units == 'W/g':
             multiplier = 1.0 / self.get_mass_density()
         elif units == 'W/kg':

openmc/mesh.py

@@ -936,6 +936,87 @@ class StructuredMesh(MeshBase):
                 f"with dimensions {self.dimension}"
             )
 
+    @classmethod
+    def from_domain(
+        cls,
+        domain: HasBoundingBox | BoundingBox,
+        dimension: Sequence[int] | int | None = None,
+        mesh_id: int | None = None,
+        name: str = '',
+        **kwargs
+    ) -> StructuredMesh:
+        """Create a structured mesh from a domain using its bounding box.
+
+        Parameters
+        ----------
+        domain : HasBoundingBox | openmc.BoundingBox
+            Object used as a template for the mesh extents. If ``domain`` has a
+            ``bounding_box`` attribute, that bounding box is used directly.
+        dimension : Iterable of int or int, optional
+            Number of mesh cells. When omitted, the subclass-specific default is
+            used. If provided as a single integer, subclasses that support it
+            interpret it as a target total number of mesh cells.
+        mesh_id : int, optional
+            Unique identifier for the mesh.
+        name : str, optional
+            Name of the mesh.
+        **kwargs
+            Additional keyword arguments forwarded to
+            :meth:`from_bounding_box`.
+
+        Returns
+        -------
+        openmc.StructuredMesh
+            Structured mesh instance.
+        """
+        if isinstance(domain, BoundingBox):
+            bbox = domain
+        elif hasattr(domain, 'bounding_box'):
+            bbox = domain.bounding_box
+        else:
+            raise TypeError("Domain must be a BoundingBox or have a "
+                            "bounding_box property")
+
+        if dimension is None:
+            return cls.from_bounding_box(
+                bbox, mesh_id=mesh_id, name=name, **kwargs)
+
+        return cls.from_bounding_box(
+            bbox, dimension=dimension, mesh_id=mesh_id, name=name, **kwargs)
+
+    @classmethod
+    @abstractmethod
+    def from_bounding_box(
+        cls,
+        bbox: openmc.BoundingBox,
+        dimension: Sequence[int] | int,
+        mesh_id: int | None = None,
+        name: str = '',
+        **kwargs
+    ) -> StructuredMesh:
+        """Create a structured mesh from a bounding box.
+
+        Parameters
+        ----------
+        bbox : openmc.BoundingBox
+            Bounding box used to define the mesh extents.
+        dimension : Iterable of int or int
+            Number of mesh cells. The interpretation and any default value are
+            defined by the concrete mesh type.
+        mesh_id : int, optional
+            Unique identifier for the mesh.
+        name : str, optional
+            Name of the mesh.
+        **kwargs
+            Additional keyword arguments accepted by specific subclasses.
+
+        Returns
+        -------
+        openmc.StructuredMesh
+            Structured mesh instance.
+        """
+        pass
+
 
 class HasBoundingBox(Protocol):
     """Object that has a ``bounding_box`` attribute."""
@@ -1190,62 +1271,47 @@ class RegularMesh(StructuredMesh):
         return mesh
 
     @classmethod
-    def from_domain(
+    def from_bounding_box(
         cls,
-        domain: HasBoundingBox | BoundingBox,
+        bbox: openmc.BoundingBox,
         dimension: Sequence[int] | int = 1000,
         mesh_id: int | None = None,
-        name: str = ''
-    ):
-        """Create RegularMesh from a domain using its bounding box.
+        name: str = '',
+    ) -> RegularMesh:
+        """Create a RegularMesh from a bounding box.
 
         Parameters
         ----------
-        domain : HasBoundingBox | openmc.BoundingBox
-            The object passed in will be used as a template for this mesh. The
-            bounding box of the property of the object passed will be used to
-            set the lower_left and upper_right and of the mesh instance.
-            Alternatively, a :class:`openmc.BoundingBox` can be passed
-            directly.
-        dimension : Iterable of int | int
-            The number of mesh cells in total or number of mesh cells in each
-            direction (x, y, z). If a single integer is provided, the domain
-            will will be divided into that many mesh cells with roughly equal
-            lengths in each direction (cubes).
-        mesh_id : int
-            Unique identifier for the mesh
-        name : str
-            Name of the mesh
+        bbox : openmc.BoundingBox
+            Bounding box used to set the mesh extents.
+        dimension : Iterable of int or int, optional
+            The number of mesh cells in each direction (x, y, z). If a single
+            integer is provided, the total number of cells is distributed
+            across directions to produce cells with roughly equal widths.
+        mesh_id : int, optional
+            Unique identifier for the mesh.
+        name : str, optional
+            Name of the mesh.
 
         Returns
         -------
         openmc.RegularMesh
-            RegularMesh instance
-
+            RegularMesh instance.
         """
-        if isinstance(domain, BoundingBox):
-            bb = domain
-        elif hasattr(domain, 'bounding_box'):
-            bb = domain.bounding_box
-        else:
-            raise TypeError("Domain must be a BoundingBox or have a "
-                            "bounding_box property")
-
         mesh = cls(mesh_id=mesh_id, name=name)
-        mesh.lower_left = bb[0]
-        mesh.upper_right = bb[1]
+        mesh.lower_left = bbox[0]
+        mesh.upper_right = bbox[1]
         if isinstance(dimension, int):
             cv.check_greater_than("dimension", dimension, 1, equality=True)
             # If a single integer is provided, divide the domain into that many
             # mesh cells with roughly equal lengths in each direction
-            ideal_cube_volume = bb.volume / dimension
+            ideal_cube_volume = bbox.volume / dimension
             ideal_cube_size = ideal_cube_volume ** (1 / 3)
             dimension = [
                 max(1, int(round(side / ideal_cube_size)))
-                for side in bb.width
+                for side in bbox.width
             ]
         mesh.dimension = dimension
-
         return mesh
 
     def to_xml_element(self):
@@ -1688,10 +1754,89 @@ class RectilinearMesh(StructuredMesh):
 
         return element
 
-    def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
-        raise NotImplementedError(
-            "get_indices_at_coords is not yet implemented for RectilinearMesh"
-        )
+    def get_indices_at_coords(self, coords: Sequence[float]) -> tuple[int, int, int]:
+        """Find the mesh cell indices containing the specified coordinates.
+
+        .. versionadded:: 0.15.4
+
+        Parameters
+        ----------
+        coords : Sequence[float]
+            Cartesian coordinates of the point as (x, y, z).
+
+        Returns
+        -------
+        tuple[int, int, int]
+            Mesh indices (ix, iy, iz).
+
+        Raises
+        ------
+        ValueError
+            If coords does not contain exactly 3 values, or if a coordinate is
+            outside the mesh grid boundaries.
+        """
+        if len(coords) != 3:
+            raise ValueError(
+                f"coords must contain exactly 3 values for a rectilinear mesh, "
+                f"got {len(coords)}"
+            )
+
+        grids = (self.x_grid, self.y_grid, self.z_grid)
+        indices = []
+
+        for grid, value in zip(grids, coords):
+            if value < grid[0] or value > grid[-1]:
+                raise ValueError(
+                    f"Coordinate value {value} is outside the mesh grid boundaries: "
+                    f"[{grid[0]}, {grid[-1]}]"
+                )
+
+            idx = np.searchsorted(grid, value, side="right") - 1
+            indices.append(int(min(idx, len(grid) - 2)))
+
+        return tuple(indices)
+
+    @classmethod
+    def from_bounding_box(
+        cls,
+        bbox: openmc.BoundingBox,
+        dimension: Sequence[int] | int = 1000,
+        mesh_id: int | None = None,
+        name: str = '',
+    ) -> RectilinearMesh:
+        """Create a RectilinearMesh from a bounding box with uniform grids.
+
+        Parameters
+        ----------
+        bbox : openmc.BoundingBox
+            Bounding box used to set the mesh extents.
+        dimension : Iterable of int or int, optional
+            The number of mesh cells in each direction (x, y, z). If a single
+            integer is provided, the total number of cells is distributed across
+            the three directions proportionally to the side lengths.
+        mesh_id : int, optional
+            Unique identifier for the mesh.
+        name : str, optional
+            Name of the mesh.
+
+        Returns
+        -------
+        openmc.RectilinearMesh
+            RectilinearMesh instance with uniform grids along each axis.
+        """
+        if isinstance(dimension, int):
+            cv.check_greater_than("dimension", dimension, 1, equality=True)
+            ideal_cube_volume = bbox.volume / dimension
+            ideal_cube_size = ideal_cube_volume ** (1 / 3)
+            dimension = [
+                max(1, int(round(side / ideal_cube_size)))
+                for side in bbox.width
+            ]
+        mesh = cls(mesh_id=mesh_id, name=name)
+        mesh.x_grid = np.linspace(bbox[0][0], bbox[1][0], num=dimension[0] + 1)
+        mesh.y_grid = np.linspace(bbox[0][1], bbox[1][1], num=dimension[1] + 1)
+        mesh.z_grid = np.linspace(bbox[0][2], bbox[1][2], num=dimension[2] + 1)
+        return mesh
 
 
 class CylindricalMesh(StructuredMesh):
@@ -1959,34 +2104,31 @@ class CylindricalMesh(StructuredMesh):
         return mesh
 
     @classmethod
-    def from_domain(
+    def from_bounding_box(
         cls,
-        domain: HasBoundingBox | BoundingBox,
+        bbox: openmc.BoundingBox,
         dimension: Sequence[int] = (10, 10, 10),
         mesh_id: int | None = None,
-        phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
         name: str = '',
-        enclose_domain: bool = False
-    ):
-        """Create CylindricalMesh from a domain using its bounding box.
+        phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
+        enclose_domain: bool = False,
+    ) -> CylindricalMesh:
+        """Create CylindricalMesh from a bounding box.
 
         Parameters
         ----------
-        domain : HasBoundingBox | openmc.BoundingBox
-            The object passed in will be used as a template for this mesh. The
-            bounding box of the property of the object passed will be used to
-            set the r_grid, z_grid ranges. Alternatively, a
-            :class:`openmc.BoundingBox` can be passed directly.
+        bbox : openmc.BoundingBox
+            Bounding box used to set the r_grid and z_grid ranges.
         dimension : Iterable of int
             The number of equally spaced mesh cells in each direction (r_grid,
             phi_grid, z_grid)
-        mesh_id : int
+        mesh_id : int, optional
             Unique identifier for the mesh
+        name : str, optional
+            Name of the mesh
         phi_grid_bounds : numpy.ndarray
             Mesh bounds points along the phi-axis in radians. The default value
             is (0, 2π), i.e., the full phi range.
-        name : str
-            Name of the mesh
         enclose_domain : bool
             If True, the mesh will encompass the bounding box of the domain. If
             False, the mesh will be inscribed within the domain's bounding box.
@@ -1997,40 +2139,28 @@ class CylindricalMesh(StructuredMesh):
             CylindricalMesh instance
 
         """
-        if isinstance(domain, BoundingBox):
-            cached_bb = domain
-        elif hasattr(domain, 'bounding_box'):
-            cached_bb = domain.bounding_box
-        else:
-            raise TypeError("Domain must be a BoundingBox or have a "
-                            "bounding_box property")
-
         if enclose_domain:
-            outer_radius = 0.5 * np.linalg.norm(cached_bb.width[:2])
+            outer_radius = 0.5 * np.linalg.norm(bbox.width[:2])
         else:
-            outer_radius = 0.5 * min(cached_bb.width[:2])
+            outer_radius = 0.5 * min(bbox.width[:2])
 
-        r_grid = np.linspace(
-            0,
-            outer_radius,
-            num=dimension[0]+1
-        )
+        r_grid = np.linspace(0, outer_radius, num=dimension[0]+1)
         phi_grid = np.linspace(
             phi_grid_bounds[0],
             phi_grid_bounds[1],
             num=dimension[1]+1
         )
         z_grid = np.linspace(
-            cached_bb[0][2],
-            cached_bb[1][2],
+            bbox[0][2],
+            bbox[1][2],
             num=dimension[2]+1
         )
-        origin = (cached_bb.center[0], cached_bb.center[1], z_grid[0])
+        origin = (bbox.center[0], bbox.center[1], z_grid[0])
 
         # make z-grid relative to the origin
         z_grid -= origin[2]
 
-        mesh = cls(
+        return cls(
             r_grid=r_grid,
             z_grid=z_grid,
             phi_grid=phi_grid,
@@ -2039,8 +2169,6 @@ class CylindricalMesh(StructuredMesh):
             origin=origin
         )
 
-        return mesh
-
     def to_xml_element(self):
         """Return XML representation of the mesh
 
@@ -2348,39 +2476,36 @@ class SphericalMesh(StructuredMesh):
         return mesh
 
     @classmethod
-    def from_domain(
+    def from_bounding_box(
         cls,
-        domain: HasBoundingBox | BoundingBox,
+        bbox: openmc.BoundingBox,
         dimension: Sequence[int] = (10, 10, 10),
         mesh_id: int | None = None,
+        name: str = '',
         phi_grid_bounds: Sequence[float] = (0.0, 2*pi),
         theta_grid_bounds: Sequence[float] = (0.0, pi),
-        name: str = '',
-        enclose_domain: bool = False
-    ):
-        """Create SphericalMesh from a domain using its bounding box.
+        enclose_domain: bool = False,
+    ) -> SphericalMesh:
+        """Create SphericalMesh from a bounding box.
 
         Parameters
         ----------
-        domain : HasBoundingBox | openmc.BoundingBox
-            The object passed in will be used as a template for this mesh. The
-            bounding box of the property of the object passed will be used to
-            set the r_grid, phi_grid, and theta_grid ranges. Alternatively, a
-            :class:`openmc.BoundingBox` can be passed directly.
+        bbox : openmc.BoundingBox
+            Bounding box used to set the r_grid, phi_grid, and theta_grid ranges.
         dimension : Iterable of int
             The number of equally spaced mesh cells in each direction (r_grid,
             phi_grid, theta_grid). Spacing is in angular space (radians) for
             phi and theta, and in absolute space for r.
-        mesh_id : int
+        mesh_id : int, optional
             Unique identifier for the mesh
+        name : str, optional
+            Name of the mesh
         phi_grid_bounds : numpy.ndarray
             Mesh bounds points along the phi-axis in radians. The default value
             is (0, 2π), i.e., the full phi range.
         theta_grid_bounds : numpy.ndarray
             Mesh bounds points along the theta-axis in radians. The default value
             is (0, π), i.e., the full theta range.
-        name : str
-            Name of the mesh
         enclose_domain : bool
             If True, the mesh will encompass the bounding box of the domain. If
             False, the mesh will be inscribed within the domain's bounding box.
@@ -2391,18 +2516,10 @@ class SphericalMesh(StructuredMesh):
             SphericalMesh instance
 
         """
-        if isinstance(domain, BoundingBox):
-            cached_bb = domain
-        elif hasattr(domain, 'bounding_box'):
-            cached_bb = domain.bounding_box
-        else:
-            raise TypeError("Domain must be a BoundingBox or have a "
-                            "bounding_box property")
-
         if enclose_domain:
-            outer_radius = 0.5 * np.linalg.norm(cached_bb.width)
+            outer_radius = 0.5 * np.linalg.norm(bbox.width)
         else:
-            outer_radius = 0.5 * min(cached_bb.width)
+            outer_radius = 0.5 * min(bbox.width)
 
         r_grid = np.linspace(0, outer_radius, num=dimension[0] + 1)
         theta_grid = np.linspace(
@@ -2415,8 +2532,7 @@ class SphericalMesh(StructuredMesh):
             phi_grid_bounds[1],
             num=dimension[2]+1
         )
-        origin = np.array([
-            cached_bb.center[0], cached_bb.center[1], cached_bb.center[2]])
+        origin = np.array([bbox.center[0], bbox.center[1], bbox.center[2]])
 
         return cls(r_grid=r_grid, phi_grid=phi_grid, theta_grid=theta_grid,
                    origin=origin, mesh_id=mesh_id, name=name)

openmc/mgxs/groups.py

@@ -78,6 +78,7 @@ class EnergyGroups:
     @group_edges.setter
     def group_edges(self, edges):
         cv.check_type('group edges', edges, Iterable, Real)
+        cv.check_increasing('group edges', edges)
         cv.check_greater_than('number of group edges', len(edges), 1)
         self._group_edges = np.array(edges)
 

openmc/model/model.py

@@ -265,6 +265,46 @@ class Model:
             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 == 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(
@@ -576,6 +616,7 @@ class Model:
         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)
 
@@ -634,6 +675,9 @@ class Model:
                           "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)
@@ -1294,8 +1338,9 @@ class Model:
         self,
         n_samples: int = 1000,
         prn_seed: int | None = None,
+        as_array: bool = False,
         **init_kwargs
-    ) -> openmc.ParticleList:
+    ) -> openmc.ParticleList | np.ndarray:
         """Sample external source and return source particles.
 
         .. versionadded:: 0.15.1
@@ -1307,13 +1352,17 @@ class Model:
         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
-            List of samples source particles
+        openmc.ParticleList or numpy.ndarray
+            List of sampled source particles, or a structured array when
+            *as_array* is True.
         """
         import openmc.lib
 
@@ -1324,7 +1373,7 @@ class Model:
 
         with openmc.lib.TemporarySession(self, **init_kwargs):
             return openmc.lib.sample_external_source(
-                n_samples=n_samples, prn_seed=prn_seed
+                n_samples=n_samples, prn_seed=prn_seed, as_array=as_array
             )
 
     def apply_tally_results(self, statepoint: PathLike | openmc.StatePoint):
@@ -2515,7 +2564,7 @@ class Model:
     def convert_to_multigroup(
         self,
         method: str = "material_wise",
-        groups: str = "CASMO-2",
+        groups: str | Sequence[float] | openmc.mgxs.EnergyGroups = "CASMO-2",
         nparticles: int = 2000,
         overwrite_mgxs_library: bool = False,
         mgxs_path: PathLike = "mgxs.h5",
@@ -2533,9 +2582,13 @@ class Model:
         ----------
         method : {"material_wise", "stochastic_slab", "infinite_medium"}, optional
             Method to generate the MGXS.
-        groups : openmc.mgxs.EnergyGroups or str, optional
-            Energy group structure for the MGXS or the name of the group
-            structure (based on keys from openmc.mgxs.GROUP_STRUCTURES).
+        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.
         overwrite_mgxs_library : bool, optional
@@ -2572,7 +2625,7 @@ class Model:
             Valid entries for temperature_settings are the same as the valid
             entries in openmc.Settings.temperature_settings.
         """
-        if isinstance(groups, str):
+        if not isinstance(groups, openmc.mgxs.EnergyGroups):
             groups = openmc.mgxs.EnergyGroups(groups)
 
         # Do all work (including MGXS generation) in a temporary directory
@@ -2588,7 +2641,7 @@ class Model:
                     # This mode doesn't require
                     # valid transport settings like particles/batches
                     original_run_mode = self.settings.run_mode
-                    self.settings.run_mode = 'volume' 
+                    self.settings.run_mode = 'volume'
                     self.init_lib(directory=tmpdir)
                     self.sync_dagmc_universes()
                     self.finalize_lib()

openmc/plots.py

@@ -197,13 +197,13 @@ _PLOT_PARAMS = dedent("""\
             Assigns colors to specific materials or cells. Keys are instances of
             :class:`Cell` or :class:`Material` and values are RGB 3-tuples, RGBA
             4-tuples, or strings indicating SVG color names. Red, green, blue,
-            and alpha should all be floats in the range [0.0, 1.0], for example:
+            and alpha should all be integers in the range [0, 255], for example:
 
             .. code-block:: python
 
                 # Make water blue
                 water = openmc.Cell(fill=h2o)
-                universe.plot(..., colors={water: (0., 0., 1.))
+                universe.plot(..., colors={water: (0, 0, 255))
         seed : int
             Seed for the random number generator
         openmc_exec : str

openmc/settings.py

@@ -41,6 +41,10 @@ class Settings:
 
     Attributes
     ----------
+    atomic_relaxation : bool
+        Whether to simulate the atomic relaxation cascade (fluorescence photons
+        and Auger electrons) following photoelectric and incoherent scattering
+        interactions.
     batches : int
         Number of batches to simulate
     confidence_intervals : bool
@@ -179,6 +183,9 @@ class Settings:
        Initial seed for randomly generated plot colors.
     ptables : bool
         Determine whether probability tables are used.
+    properties_file : PathLike
+        Location of the properties file to load cell temperatures/densities
+        and material densities
     random_ray : dict
         Options for configuring the random ray solver. Acceptable keys are:
 
@@ -229,6 +236,9 @@ class Settings:
             stabilization, which may be desirable as stronger diagonal stabilization
             also tends to dampen the convergence rate of the solver, thus requiring
             more iterations to converge.
+        :adjoint_source:
+            Source object used to define localized adjoint source/detector response 
+            function.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -402,8 +412,10 @@ class Settings:
         self._confidence_intervals = None
         self._electron_treatment = None
         self._photon_transport = None
+        self._atomic_relaxation = None
         self._plot_seed = None
         self._ptables = None
+        self._properties_file = None
         self._uniform_source_sampling = None
         self._seed = None
         self._stride = None
@@ -663,6 +675,15 @@ class Settings:
                        electron_treatment, ['led', 'ttb'])
         self._electron_treatment = electron_treatment
 
+    @property
+    def atomic_relaxation(self) -> bool:
+        return self._atomic_relaxation
+
+    @atomic_relaxation.setter
+    def atomic_relaxation(self, atomic_relaxation: bool):
+        cv.check_type('atomic relaxation', atomic_relaxation, bool)
+        self._atomic_relaxation = atomic_relaxation
+
     @property
     def ptables(self) -> bool:
         return self._ptables
@@ -1053,6 +1074,18 @@ class Settings:
 
         self._temperature = temperature
 
+    @property
+    def properties_file(self) -> PathLike | None:
+        return self._properties_file
+
+    @properties_file.setter
+    def properties_file(self, value: PathLike | None):
+        if value is None:
+            self._properties_file = None
+        else:
+            cv.check_type('properties file', value, PathLike)
+            self._properties_file = input_path(value)
+
     @property
     def trace(self) -> Iterable:
         return self._trace
@@ -1391,6 +1424,14 @@ class Settings:
                 cv.check_type('diagonal stabilization rho', value, Real)
                 cv.check_greater_than('diagonal stabilization rho',
                                       value, 0.0, True)
+            elif key == 'adjoint_source':
+                if not isinstance(value, MutableSequence):
+                    value = [value]
+                for source in value:
+                    if not isinstance(source, SourceBase):
+                        raise ValueError(
+                            f'Invalid adjoint source type: {type(source)}. '
+                            'Expected openmc.SourceBase.')
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1631,6 +1672,11 @@ class Settings:
             element = ET.SubElement(root, "electron_treatment")
             element.text = str(self._electron_treatment)
 
+    def _create_atomic_relaxation_subelement(self, root):
+        if self._atomic_relaxation is not None:
+            element = ET.SubElement(root, "atomic_relaxation")
+            element.text = str(self._atomic_relaxation).lower()
+
     def _create_photon_transport_subelement(self, root):
         if self._photon_transport is not None:
             element = ET.SubElement(root, "photon_transport")
@@ -1753,6 +1799,12 @@ class Settings:
                 else:
                     element.text = str(value)
 
+    def _create_properties_file_element(self, root):
+        if self.properties_file is not None:
+            element = ET.Element("properties_file")
+            element.text = str(self.properties_file)
+            root.append(element)
+
     def _create_trace_subelement(self, root):
         if self._trace is not None:
             element = ET.SubElement(root, "trace")
@@ -1932,11 +1984,12 @@ class Settings:
             element = ET.SubElement(root, "random_ray")
             for key, value in self._random_ray.items():
                 if key == 'ray_source' and isinstance(value, SourceBase):
+                    subelement = ET.SubElement(element, 'ray_source')
                     source_element = value.to_xml_element()
                     if source_element.find('bias') is not None:
                         raise RuntimeError(
                             "Ray source distributions should not be biased.")
-                    element.append(source_element)
+                    subelement.append(source_element)
 
                 elif key == 'source_region_meshes':
                     subelement = ET.SubElement(element, 'source_region_meshes')
@@ -1954,8 +2007,20 @@ class Settings:
                         path = f"./mesh[@id='{mesh.id}']"
                         if root.find(path) is None:
                             root.append(mesh.to_xml_element())
-                            if mesh_memo is not None:
+                            if mesh_memo is not None:    
                                 mesh_memo.add(mesh.id)
+                elif key == 'adjoint_source':
+                    subelement = ET.SubElement(element, 'adjoint_source')
+                    # Check that all entries are valid SourceBase instances, in case 
+                    # the random_ray setter was not used to populate dict entries.
+                    if not isinstance(value, MutableSequence):
+                        value = [value]
+                    for source in value:
+                        if not isinstance(source, SourceBase):
+                            raise ValueError(
+                                f'Invalid adjoint source type: {type(source)}. '
+                                'Expected openmc.SourceBase.')
+                        subelement.append(source.to_xml_element())
                 elif isinstance(value, bool):
                     subelement = ET.SubElement(element, key)
                     subelement.text = str(value).lower()
@@ -2129,6 +2194,11 @@ class Settings:
         if text is not None:
             self.electron_treatment = text
 
+    def _atomic_relaxation_from_xml_element(self, root):
+        text = get_text(root, 'atomic_relaxation')
+        if text is not None:
+            self.atomic_relaxation = text in ('true', '1')
+
     def _energy_mode_from_xml_element(self, root):
         text = get_text(root, 'energy_mode')
         if text is not None:
@@ -2260,6 +2330,11 @@ class Settings:
         if text is not None:
             self.temperature['multipole'] = text in ('true', '1')
 
+    def _properties_file_from_xml_element(self, root):
+        text = get_text(root, 'properties_file')
+        if text is not None:
+            self.properties_file = text
+
     def _trace_from_xml_element(self, root):
         text = get_elem_list(root, "trace", int)
         if text is not None:
@@ -2392,8 +2467,9 @@ class Settings:
             for child in elem:
                 if child.tag in ('distance_inactive', 'distance_active', 'diagonal_stabilization_rho'):
                     self.random_ray[child.tag] = float(child.text)
-                elif child.tag == 'source':
-                    source = SourceBase.from_xml_element(child)
+                elif child.tag == 'ray_source':
+                    source_element = child.find('source')
+                    source = SourceBase.from_xml_element(source_element)
                     if child.find('bias') is not None:
                         raise RuntimeError(
                             "Ray source distributions should not be biased.")
@@ -2410,6 +2486,12 @@ class Settings:
                     self.random_ray['adjoint'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'adjoint_source':
+                    self.random_ray['adjoint_source'] = []
+                    for subelem in child.findall('source'):
+                        src = SourceBase.from_xml_element(subelem)
+                        # add newly constructed source object to the list
+                        self.random_ray['adjoint_source'].append(src)
                 elif child.tag == 'sample_method':
                     self.random_ray['sample_method'] = child.text
                 elif child.tag == 'source_region_meshes':
@@ -2478,6 +2560,7 @@ class Settings:
         self._create_collision_track_subelement(element)
         self._create_confidence_intervals(element)
         self._create_electron_treatment_subelement(element)
+        self._create_atomic_relaxation_subelement(element)
         self._create_energy_mode_subelement(element)
         self._create_max_order_subelement(element)
         self._create_photon_transport_subelement(element)
@@ -2497,6 +2580,7 @@ class Settings:
         self._create_ifp_n_generation_subelement(element)
         self._create_tabular_legendre_subelements(element)
         self._create_temperature_subelements(element)
+        self._create_properties_file_element(element)
         self._create_trace_subelement(element)
         self._create_track_subelement(element)
         self._create_ufs_mesh_subelement(element, mesh_memo)
@@ -2594,6 +2678,7 @@ class Settings:
         settings._collision_track_from_xml_element(elem)
         settings._confidence_intervals_from_xml_element(elem)
         settings._electron_treatment_from_xml_element(elem)
+        settings._atomic_relaxation_from_xml_element(elem)
         settings._energy_mode_from_xml_element(elem)
         settings._max_order_from_xml_element(elem)
         settings._photon_transport_from_xml_element(elem)
@@ -2613,6 +2698,7 @@ class Settings:
         settings._ifp_n_generation_from_xml_element(elem)
         settings._tabular_legendre_from_xml_element(elem)
         settings._temperature_from_xml_element(elem)
+        settings._properties_file_from_xml_element(elem)
         settings._trace_from_xml_element(elem)
         settings._track_from_xml_element(elem)
         settings._ufs_mesh_from_xml_element(elem, meshes)

openmc/stats/univariate.py

@@ -2,8 +2,7 @@ from __future__ import annotations
 from abc import ABC, abstractmethod
 from collections import defaultdict
 from collections.abc import Iterable, Sequence
-from copy import deepcopy
-from math import sqrt, pi, exp
+from math import sqrt, pi, exp, log
 from numbers import Real
 from warnings import warn
 
@@ -14,6 +13,7 @@ from scipy.special import exprel, hyp1f1, lambertw
 import scipy
 
 import openmc.checkvalue as cv
+from openmc.data import atomic_mass, NEUTRON_MASS
 from .._xml import get_elem_list, get_text
 from ..mixin import EqualityMixin
 
@@ -1277,6 +1277,138 @@ def Muir(*args, **kwargs):
     return muir(*args, **kwargs)
 
 
+def fusion_neutron_spectrum(
+    ion_temp: float,
+    reactants: str = 'DD',
+    bias: Univariate | None = None
+) -> Normal:
+    r"""Return a Gaussian energy distribution for fusion neutron emission.
+
+    Computes the mean energy and spectral width of the neutron energy spectrum
+    from thermonuclear fusion reactions in a plasma with Maxwellian ion velocity
+    distributions. The mean neutron energy is calculated as
+
+    .. math::
+
+        \langle E_n \rangle = E_0 + \Delta E_\text{th}
+
+    where :math:`E_0` is the neutron energy at zero ion temperature and
+    :math:`\Delta E_\text{th}` is the thermal peak shift due to the motion of
+    the reacting ions. The spectral width is characterized by the FWHM:
+
+    .. math::
+
+        W_{1/2} = \omega_0 (1 + \delta_\omega) \sqrt{T_i}
+
+    where :math:`\omega_0` is the width at the :math:`T_i \to 0` limit and
+    :math:`\delta_\omega` is a temperature-dependent correction term. Both
+    :math:`\Delta E_\text{th}` and :math:`\delta_\omega` are evaluated using
+    interpolation formulas from `Ballabio et al.
+    <https://doi.org/10.1088/0029-5515/38/11/310>`_: Table III for :math:`0 <
+    T_i \le 40` keV and Table IV for :math:`40 < T_i < 100` keV. The returned
+    distribution is a normal (Gaussian) approximation to the spectrum.
+
+    .. versionadded:: 0.15.4
+
+    Parameters
+    ----------
+    ion_temp : float
+        Ion temperature of the plasma in [eV].
+    reactants : {'DD', 'DT'}
+        Fusion reactants. 'DD' corresponds to the D(d,n)\ :sup:`3`\ He reaction
+        and 'DT' to the T(d,n)\ :math:`\alpha` reaction.
+    bias : openmc.stats.Univariate, optional
+        Distribution for biased sampling.
+
+    Returns
+    -------
+    openmc.stats.Normal
+        Normal distribution with mean and standard deviation corresponding to
+        the first and second moments of the fusion neutron energy spectrum. Both
+        the mean and standard deviation are in [eV].
+
+    """
+    if ion_temp < 0.0 or ion_temp > 100e3:
+        raise ValueError("Ion temperature must be between 0 and 100 keV.")
+
+    # Formulas from doi:10.1088/0029-5515/38/11/310
+    mn = NEUTRON_MASS
+    md = atomic_mass('H2')
+    ev_per_c2 = 931.49410372*1e6
+    if reactants == 'DD':
+        mhe3 = atomic_mass('He3')
+        Q = (md + md - mhe3 - mn)*ev_per_c2
+        E_n = mhe3/(mhe3 + mn)*Q
+        w0 = 82.542
+
+        # Low-T constants for peak shift (Table III)
+        a1 = 4.69515
+        a2 = -0.040729
+        a3 = 0.47
+        a4 = 0.81844
+
+        # Low-T constants for width correction (Table III)
+        b1 = 1.7013e-3
+        b2 = 0.16888
+        b3 = 0.49
+        b4 = 7.9460e-4
+
+        # High-T constants for peak shift (Table IV)
+        a5 = 18.225
+        a6 = 2.1525
+
+        # High-T constants for width correction (Table IV)
+        b5 = 8.4619e-3
+        b6 = 8.3241e-4
+
+    elif reactants == 'DT':
+        mt = atomic_mass('H3')
+        ma = atomic_mass('He4')
+        Q = (md + mt - ma - mn)*ev_per_c2
+        E_n = ma/(ma + mn)*Q
+        w0 = 177.259
+
+        # Low-T constants for peak shift (Table III)
+        a1 = 5.30509
+        a2 = 2.4736e-3
+        a3 = 1.84
+        a4 = 1.3818
+
+        # Low-T constants for width correction (Table III)
+        b1 = 5.1068e-4
+        b2 = 7.6223e-3
+        b3 = 1.78
+        b4 = 8.7691e-5
+
+        # High-T constants for peak shift (Table IV)
+        a5 = 37.771
+        a6 = 0.92181
+
+        # High-T constants for width correction (Table IV)
+        b5 = 2.0199e-3
+        b6 = 5.9501e-5
+    else:
+        raise ValueError("Invalid reactants specified. Must be 'DD' or 'DT'.")
+
+    # Ion temperature in keV
+    T = ion_temp * 1e-3
+
+    if T <= 40.0:
+        # Low-temperature interpolation (Table III, 0 < T_i <= 40 keV)
+        Delta_E = a1/(1 + a2*T**a3)*T**(2/3) + a4*T
+        delta_w = b1/(1 + b2*T**b3)*T**(2/3) + b4*T
+    else:
+        # High-temperature interpolation (Table IV, 40 < T_i < 100 keV)
+        Delta_E = a5 + a6*T
+        delta_w = b5 + b6*T
+
+    # Calculate FWHM
+    fwhm = (w0*(1 + delta_w) * sqrt(T))*1e3
+
+    sigma = fwhm / (2*sqrt(2*log(2)))
+    return Normal(E_n + Delta_E * 1e3, sigma, bias=bias)
+
+
 class Tabular(Univariate):
     """Piecewise continuous probability distribution.
 

openmc/tallies.py

@@ -16,6 +16,23 @@ from scipy.stats import chi2, norm
 
 import openmc
 import openmc.checkvalue as cv
+from openmc.filter import (
+    Filter, 
+    DistribcellFilter, 
+    EnergyFunctionFilter, 
+    DelayedGroupFilter, 
+    FilterMeta, 
+    MeshFilter,
+    MeshBornFilter,
+)
+from openmc.arithmetic import (
+    CrossFilter, 
+    AggregateFilter,
+    CrossScore,
+    AggregateScore,
+    CrossNuclide,
+    AggregateNuclide,
+)
 from ._sparse_compat import lil_array
 from ._xml import clean_indentation, get_elem_list, get_text
 from .mixin import IDManagerMixin
@@ -31,9 +48,9 @@ _PRODUCT_TYPES = ['tensor', 'entrywise']
 
 # The following indicate acceptable types when setting Tally.scores,
 # Tally.nuclides, and Tally.filters
-_SCORE_CLASSES = (str, openmc.CrossScore, openmc.AggregateScore)
-_NUCLIDE_CLASSES = (str, openmc.CrossNuclide, openmc.AggregateNuclide)
-_FILTER_CLASSES = (openmc.Filter, openmc.CrossFilter, openmc.AggregateFilter)
+_SCORE_CLASSES = (str, CrossScore, AggregateScore)
+_NUCLIDE_CLASSES = (str, CrossNuclide, AggregateNuclide)
+_FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter)
 
 # Valid types of estimators
 ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog'}
@@ -421,7 +438,7 @@ class Tally(IDManagerMixin):
             self._num_realizations = int(group['n_realizations'][()])
 
             for filt in self.filters:
-                if isinstance(filt, openmc.DistribcellFilter):
+                if isinstance(filt, DistribcellFilter):
                     filter_group = f[f'tallies/filters/filter {filt.id}']
                     filt._num_bins = int(filter_group['n_bins'][()])
 
@@ -1089,8 +1106,8 @@ class Tally(IDManagerMixin):
             return False
 
         # Return False if only one tally has a delayed group filter
-        tally1_dg = self.contains_filter(openmc.DelayedGroupFilter)
-        tally2_dg = other.contains_filter(openmc.DelayedGroupFilter)
+        tally1_dg = self.contains_filter(DelayedGroupFilter)
+        tally2_dg = other.contains_filter(DelayedGroupFilter)
         if tally1_dg != tally2_dg:
             return False
 
@@ -1602,7 +1619,7 @@ class Tally(IDManagerMixin):
 
             # Also check to see if the desired filter is wrapped up in an
             # aggregate
-            elif isinstance(test_filter, openmc.AggregateFilter):
+            elif isinstance(test_filter, AggregateFilter):
                 if isinstance(test_filter.aggregate_filter, filter_type):
                     return test_filter
 
@@ -1704,7 +1721,7 @@ class Tally(IDManagerMixin):
 
         """
 
-        cv.check_type('filters', filters, Iterable, openmc.FilterMeta)
+        cv.check_type('filters', filters, Iterable, FilterMeta)
         cv.check_type('filter_bins', filter_bins, Iterable, tuple)
 
         # If user did not specify any specific Filters, use them all
@@ -1787,7 +1804,7 @@ class Tally(IDManagerMixin):
         """
 
         for score in scores:
-            if not isinstance(score, (str, openmc.CrossScore)):
+            if not isinstance(score, (str, CrossScore)):
                 msg = f'Unable to get score indices for score "{score}" in ' \
                       f'ID="{self.id}" since it is not a string or CrossScore ' \
                       'Tally'
@@ -1984,9 +2001,9 @@ class Tally(IDManagerMixin):
             column_name = 'score'
 
             for score in self.scores:
-                if isinstance(score, (str, openmc.CrossScore)):
+                if isinstance(score, (str, CrossScore)):
                     scores.append(str(score))
-                elif isinstance(score, openmc.AggregateScore):
+                elif isinstance(score, AggregateScore):
                     scores.append(score.name)
                     column_name = f'{score.aggregate_op}(score)'
 
@@ -2086,7 +2103,7 @@ class Tally(IDManagerMixin):
         for i, f in enumerate(self.filters):
             if expand_dims:
                 # Mesh filter indices are backwards so we need to flip them
-                if type(f) in {openmc.MeshFilter, openmc.MeshBornFilter}:
+                if type(f) in {MeshFilter, MeshBornFilter}:
                     fshape = f.shape[::-1]
                     new_shape += fshape
                     idx0, idx1 = i, i + len(fshape) - 1
@@ -2273,7 +2290,7 @@ class Tally(IDManagerMixin):
         else:
             all_filters = [self_copy.filters, other_copy.filters]
             for self_filter, other_filter in product(*all_filters):
-                new_filter = openmc.CrossFilter(self_filter, other_filter,
+                new_filter = CrossFilter(self_filter, other_filter,
                                                 binary_op)
                 new_tally.filters.append(new_filter)
 
@@ -2284,7 +2301,7 @@ class Tally(IDManagerMixin):
         else:
             all_nuclides = [self_copy.nuclides, other_copy.nuclides]
             for self_nuclide, other_nuclide in product(*all_nuclides):
-                new_nuclide = openmc.CrossNuclide(self_nuclide, other_nuclide,
+                new_nuclide = CrossNuclide(self_nuclide, other_nuclide,
                                                   binary_op)
                 new_tally.nuclides.append(new_nuclide)
 
@@ -2295,9 +2312,9 @@ class Tally(IDManagerMixin):
                 if score1 == score2:
                     return score1
                 else:
-                    return openmc.CrossScore(score1, score2, binary_op)
+                    return CrossScore(score1, score2, binary_op)
             else:
-                return openmc.CrossScore(score1, score2, binary_op)
+                return CrossScore(score1, score2, binary_op)
 
         # Add scores to the new tally
         if score_product == 'entrywise':
@@ -2506,16 +2523,16 @@ class Tally(IDManagerMixin):
 
         # Construct lists of tuples for the bins in each of the two filters
         filters = [type(filter1), type(filter2)]
-        if isinstance(filter1, openmc.DistribcellFilter):
+        if isinstance(filter1, DistribcellFilter):
             filter1_bins = [b for b in range(filter1.num_bins)]
-        elif isinstance(filter1, openmc.EnergyFunctionFilter):
+        elif isinstance(filter1, EnergyFunctionFilter):
             filter1_bins = [None]
         else:
             filter1_bins = filter1.bins
 
-        if isinstance(filter2, openmc.DistribcellFilter):
+        if isinstance(filter2, DistribcellFilter):
             filter2_bins = [b for b in range(filter2.num_bins)]
-        elif isinstance(filter2, openmc.EnergyFunctionFilter):
+        elif isinstance(filter2, EnergyFunctionFilter):
             filter2_bins = [None]
         else:
             filter2_bins = filter2.bins
@@ -2648,11 +2665,11 @@ class Tally(IDManagerMixin):
             raise ValueError(msg)
 
         # Check that the scores are valid
-        if not isinstance(score1, (str, openmc.CrossScore)):
+        if not isinstance(score1, (str, CrossScore)):
             msg = 'Unable to swap score1 "{}" in Tally ID="{}" since it is ' \
                   'not a string or CrossScore'.format(score1, self.id)
             raise ValueError(msg)
-        elif not isinstance(score2, (str, openmc.CrossScore)):
+        elif not isinstance(score2, (str, CrossScore)):
             msg = 'Unable to swap score2 "{}" in Tally ID="{}" since it is ' \
                   'not a string or CrossScore'.format(score2, self.id)
             raise ValueError(msg)
@@ -3296,7 +3313,7 @@ class Tally(IDManagerMixin):
                 new_filter.bins = [f.bins[i] for i in bin_indices]
 
                 # Set number of bins manually for mesh/distribcell filters
-                if filter_type is openmc.DistribcellFilter:
+                if filter_type is DistribcellFilter:
                     new_filter._num_bins = f._num_bins
 
                 # Replace existing filter with new one
@@ -3362,16 +3379,16 @@ class Tally(IDManagerMixin):
         std_dev = self.get_reshaped_data(value='std_dev')
 
         # Sum across any filter bins specified by the user
-        if isinstance(filter_type, openmc.FilterMeta):
+        if isinstance(filter_type, FilterMeta):
             find_filter = self.find_filter(filter_type)
 
             # If user did not specify filter bins, sum across all bins
             if len(filter_bins) == 0:
                 bin_indices = np.arange(find_filter.num_bins)
 
-                if isinstance(find_filter, openmc.DistribcellFilter):
+                if isinstance(find_filter, DistribcellFilter):
                     filter_bins = np.arange(find_filter.num_bins)
-                elif isinstance(find_filter, openmc.EnergyFunctionFilter):
+                elif isinstance(find_filter, EnergyFunctionFilter):
                     filter_bins = [None]
                 else:
                     filter_bins = find_filter.bins
@@ -3400,7 +3417,7 @@ class Tally(IDManagerMixin):
 
                     # Add AggregateFilter to the tally sum
                     if not remove_filter:
-                        filter_sum = openmc.AggregateFilter(self_filter,
+                        filter_sum = AggregateFilter(self_filter,
                             [tuple(filter_bins)], 'sum')
                         tally_sum.filters.append(filter_sum)
 
@@ -3423,7 +3440,7 @@ class Tally(IDManagerMixin):
             std_dev = np.sqrt(std_dev)
 
             # Add AggregateNuclide to the tally sum
-            nuclide_sum = openmc.AggregateNuclide(nuclides, 'sum')
+            nuclide_sum = AggregateNuclide(nuclides, 'sum')
             tally_sum.nuclides.append(nuclide_sum)
 
         # Add a copy of this tally's nuclides to the tally sum
@@ -3441,7 +3458,7 @@ class Tally(IDManagerMixin):
             std_dev = np.sqrt(std_dev)
 
             # Add AggregateScore to the tally sum
-            score_sum = openmc.AggregateScore(scores, 'sum')
+            score_sum = AggregateScore(scores, 'sum')
             tally_sum.scores.append(score_sum)
 
         # Add a copy of this tally's scores to the tally sum
@@ -3514,16 +3531,16 @@ class Tally(IDManagerMixin):
         std_dev = self.get_reshaped_data(value='std_dev')
 
         # Average across any filter bins specified by the user
-        if isinstance(filter_type, openmc.FilterMeta):
+        if isinstance(filter_type, FilterMeta):
             find_filter = self.find_filter(filter_type)
 
             # If user did not specify filter bins, average across all bins
             if len(filter_bins) == 0:
                 bin_indices = np.arange(find_filter.num_bins)
 
-                if isinstance(find_filter, openmc.DistribcellFilter):
+                if isinstance(find_filter, DistribcellFilter):
                     filter_bins = np.arange(find_filter.num_bins)
-                elif isinstance(find_filter, openmc.EnergyFunctionFilter):
+                elif isinstance(find_filter, EnergyFunctionFilter):
                     filter_bins = [None]
                 else:
                     filter_bins = find_filter.bins
@@ -3553,7 +3570,7 @@ class Tally(IDManagerMixin):
 
                     # Add AggregateFilter to the tally avg
                     if not remove_filter:
-                        filter_sum = openmc.AggregateFilter(self_filter,
+                        filter_sum = AggregateFilter(self_filter,
                             [tuple(filter_bins)], 'avg')
                         tally_avg.filters.append(filter_sum)
 
@@ -3577,7 +3594,7 @@ class Tally(IDManagerMixin):
             std_dev = np.sqrt(std_dev)
 
             # Add AggregateNuclide to the tally avg
-            nuclide_avg = openmc.AggregateNuclide(nuclides, 'avg')
+            nuclide_avg = AggregateNuclide(nuclides, 'avg')
             tally_avg.nuclides.append(nuclide_avg)
 
         # Add a copy of this tally's nuclides to the tally avg
@@ -3596,7 +3613,7 @@ class Tally(IDManagerMixin):
             std_dev = np.sqrt(std_dev)
 
             # Add AggregateScore to the tally avg
-            score_sum = openmc.AggregateScore(scores, 'avg')
+            score_sum = AggregateScore(scores, 'avg')
             tally_avg.scores.append(score_sum)
 
         # Add a copy of this tally's scores to the tally avg
@@ -3786,7 +3803,7 @@ class Tallies(cv.CheckedList):
         already_written = memo if memo else set()
         for tally in self:
             for f in tally.filters:
-                if isinstance(f, openmc.MeshFilter):
+                if isinstance(f, MeshFilter):
                     if f.mesh.id in already_written:
                         continue
                     if len(f.mesh.name) > 0:
@@ -3881,7 +3898,7 @@ class Tallies(cv.CheckedList):
         # Read filter elements
         filters = {}
         for e in elem.findall('filter'):
-            filter = openmc.Filter.from_xml_element(e, meshes=meshes)
+            filter = Filter.from_xml_element(e, meshes=meshes)
             filters[filter.id] = filter
 
         # Read derivative elements

openmc/weight_windows.py

@@ -11,6 +11,7 @@ import h5py
 
 import openmc
 from openmc.mesh import MeshBase, RectilinearMesh, CylindricalMesh, SphericalMesh, UnstructuredMesh
+from openmc.tallies import Tallies
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
 from ._xml import get_elem_list, get_text, clean_indentation
@@ -499,6 +500,8 @@ class WeightWindowGenerator:
         Particle type the weight windows apply to
     method : {'magic', 'fw_cadis'}
         The weight window generation methodology applied during an update.
+    targets : :class:`openmc.Tallies` or iterable of int
+        Target tallies for local variance reduction via FW-CADIS.
     max_realizations : int
         The upper limit for number of tally realizations when generating weight
         windows.
@@ -518,6 +521,8 @@ class WeightWindowGenerator:
         Particle type the weight windows apply to
     method : {'magic', 'fw_cadis'}
         The weight window generation methodology applied during an update.
+    targets : :class:`openmc.Tallies` or numpy.ndarray
+        Target tallies for local variance reduction via FW-CADIS.
     max_realizations : int
         The upper limit for number of tally realizations when generating weight
         windows.
@@ -529,7 +534,7 @@ class WeightWindowGenerator:
         Whether or not to apply weight windows on the fly.
     """
 
-    _MAGIC_PARAMS = {'value': str, 'threshold': float, 'ratio': float}
+    _WWG_PARAMS = {'value': str, 'threshold': float, 'ratio': float}
 
     def __init__(
         self,
@@ -537,6 +542,7 @@ class WeightWindowGenerator:
         energy_bounds: Sequence[float] | None = None,
         particle_type: str | int | openmc.ParticleType = 'neutron',
         method: str = 'magic',
+        targets: openmc.Tallies | Iterable[int] | None = None,
         max_realizations: int = 1,
         update_interval: int = 1,
         on_the_fly: bool = True
@@ -549,6 +555,7 @@ class WeightWindowGenerator:
             self.energy_bounds = energy_bounds
         self.particle_type = particle_type
         self.method = method
+        self.targets = targets
         self.max_realizations = max_realizations
         self.update_interval = update_interval
         self.on_the_fly = on_the_fly
@@ -611,6 +618,22 @@ class WeightWindowGenerator:
                 self._check_update_parameters()
             except (TypeError, KeyError):
                 warnings.warn(f'Update parameters are invalid for the "{m}" method.')
+    
+    @property
+    def targets(self) -> openmc.Tallies:
+        return self._targets
+
+    @targets.setter
+    def targets(self, t):
+        if t is None:
+            self._targets = t
+        else:
+            cv.check_type('Local FW-CADIS target tallies', t, Iterable)
+            cv.check_greater_than('Local FW-CADIS target tallies', len(t), 0)
+            if not isinstance(t, openmc.Tallies):
+                cv.check_iterable_type('Local FW-CADIS target tallies', t, int)
+                t = np.asarray(list(t), dtype=int)
+            self._targets = t
 
     @property
     def max_realizations(self) -> int:
@@ -638,13 +661,13 @@ class WeightWindowGenerator:
 
     def _check_update_parameters(self, params: dict):
         if self.method == 'magic' or self.method == 'fw_cadis':
-            check_params = self._MAGIC_PARAMS
+            check_params = self._WWG_PARAMS
 
         for key, val in params.items():
             if key not in check_params:
                 raise ValueError(f'Invalid param "{key}" for {self.method} '
                                   'weight window generation')
-            cv.check_type(f'weight window generation param: "{key}"', val, self._MAGIC_PARAMS[key])
+            cv.check_type(f'weight window generation param: "{key}"', val, self._WWG_PARAMS[key])
 
     @update_parameters.setter
     def update_parameters(self, params: dict):
@@ -681,7 +704,7 @@ class WeightWindowGenerator:
             The update parameters as-read from the XML node (keys: str, values: str)
         """
         if method == 'magic' or method == 'fw_cadis':
-            check_params = cls._MAGIC_PARAMS
+            check_params = cls._WWG_PARAMS
 
         for param, param_type in check_params.items():
             if param in update_parameters:
@@ -707,6 +730,20 @@ class WeightWindowGenerator:
         otf_elem.text = str(self.on_the_fly).lower()
         method_elem = ET.SubElement(element, 'method')
         method_elem.text = self.method
+        if self.targets is not None:
+            if self.method != 'fw_cadis':
+                raise ValueError(
+                    "FW-CADIS update method is required in order to use " \
+                    "target tallies for WeightWindowGenerator.")
+            elif isinstance(self.targets, openmc.Tallies):
+                raise RuntimeError(
+                    "FW-CADIS target tallies must be checked to ensure they are " \
+                    "present on model.tallies. Use model.export_to_xml() or " \
+                    "model.export_to_model_xml() to link FW-CADIS target tallies.")
+            else:
+                targets_elem = ET.SubElement(element, 'targets')
+                targets_elem.text = ' '.join(str(tally_id) for tally_id in self.targets)
+
         if self.update_parameters is not None:
             self._update_parameters_subelement(element)
 
@@ -733,8 +770,8 @@ class WeightWindowGenerator:
 
         mesh_id = int(get_text(elem, 'mesh'))
         mesh = meshes[mesh_id]
-
-        energy_bounds = get_elem_list(elem, "energy_bounds, float")
+        
+        energy_bounds = get_elem_list(elem, "energy_bounds", float)
         particle_type = get_text(elem, 'particle_type')
 
         wwg = cls(mesh, energy_bounds, particle_type)
@@ -743,6 +780,14 @@ class WeightWindowGenerator:
         wwg.update_interval = int(get_text(elem, 'update_interval'))
         wwg.on_the_fly = bool(get_text(elem, 'on_the_fly'))
         wwg.method = get_text(elem, 'method')
+        targets_elem = elem.find('targets')
+        if targets_elem is not None:
+            if wwg.method != 'fw_cadis':
+                raise ValueError(
+                    "FW-CADIS update method is required in order to use " \
+                    "target tallies for WeightWindowGenerator.")
+            else:
+                wwg.targets = get_elem_list(elem, "targets")
 
         if elem.find('update_parameters') is not None:
             update_parameters = {}

src/chain.cpp

@@ -74,6 +74,13 @@ void DecayPhotonAngleEnergy::sample(
   mu = Uniform(-1., 1.).sample(seed).first;
 }
 
+double DecayPhotonAngleEnergy::sample_energy_and_pdf(
+  double E_in, double mu, double& E_out, uint64_t* seed) const
+{
+  E_out = photon_energy_->sample(seed).first;
+  return 0.5;
+}
+
 //==============================================================================
 // Global variables
 //==============================================================================

src/distribution_angle.cpp

@@ -82,4 +82,19 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
   return mu;
 }
 
+double AngleDistribution::evaluate(double E, double mu) const
+{
+  // Find energy bin and calculate interpolation factor
+  int i;
+  double r;
+  get_energy_index(energy_, E, i, r);
+
+  double pdf = 0.0;
+  if (r > 0.0)
+    pdf += r * distribution_[i + 1]->evaluate(mu);
+  if (r < 1.0)
+    pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
+  return pdf;
+}
+
 } // namespace openmc

src/distribution_multi.cpp

@@ -1,6 +1,6 @@
 #include "openmc/distribution_multi.h"
 
-#include <algorithm> // for move
+#include <algorithm> // for move, clamp
 #include <cmath>     // for sqrt, sin, cos, max
 
 #include "openmc/constants.h"
@@ -44,6 +44,7 @@ UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)
       fatal_error("Angular distribution reference direction must have "
                   "three parameters specified.");
     u_ref_ = Direction(u_ref.data());
+    u_ref_ /= u_ref_.norm();
   }
 }
 
@@ -65,6 +66,7 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
       fatal_error("Angular distribution reference v direction must have "
                   "three parameters specified.");
     v_ref_ = Direction(v_ref.data());
+    v_ref_ /= v_ref_.norm();
   }
   w_ref_ = u_ref_.cross(v_ref_);
   if (check_for_node(node, "mu")) {
@@ -116,6 +118,22 @@ std::pair<Direction, double> PolarAzimuthal::sample_impl(
     weight};
 }
 
+double PolarAzimuthal::evaluate(Direction u) const
+{
+  double mu = std::clamp(u.dot(u_ref_), -1.0, 1.0);
+  double phi = 0.0;
+  double sin_theta_sq = std::max(0.0, 1.0 - mu * mu);
+  if (sin_theta_sq > 0.0) {
+    double sin_theta = std::sqrt(sin_theta_sq);
+    double cos_phi = u.dot(v_ref_) / sin_theta;
+    double sin_phi = u.dot(w_ref_) / sin_theta;
+    phi = std::atan2(sin_phi, cos_phi);
+    if (phi < 0.0)
+      phi += 2.0 * PI;
+  }
+  return mu_->evaluate(mu) * phi_->evaluate(phi);
+}
+
 //==============================================================================
 // Isotropic implementation
 //==============================================================================
@@ -157,6 +175,11 @@ std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const
   }
 }
 
+double Isotropic::evaluate(Direction u) const
+{
+  return 1.0 / (4.0 * PI);
+}
+
 //==============================================================================
 // Monodirectional implementation
 //==============================================================================

src/finalize.cpp

@@ -140,6 +140,7 @@ int openmc_finalize()
   settings::temperature_multipole = false;
   settings::temperature_range = {0.0, 0.0};
   settings::temperature_tolerance = 10.0;
+  settings::properties_file.clear();
   settings::trigger_on = false;
   settings::trigger_predict = false;
   settings::trigger_batch_interval = 1;

src/initialize.cpp

@@ -118,6 +118,10 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
   if (!read_model_xml())
     read_separate_xml_files();
 
+  if (!settings::properties_file.empty()) {
+    openmc_properties_import(settings::properties_file.c_str());
+  }
+
   // Reset locale to previous state
   if (std::setlocale(LC_ALL, prev_locale.c_str()) == NULL) {
     fatal_error("Cannot reset locale.");

src/lattice.cpp

@@ -340,6 +340,26 @@ Position RectLattice::get_local_position(
 
 //==============================================================================
 
+Direction RectLattice::get_normal(
+  const array<int, 3>& i_xyz, bool& is_valid) const
+{
+  is_valid = false;
+  Direction dir = {0.0, 0.0, 0.0};
+  if ((std::abs(i_xyz[0]) == 1) && (i_xyz[1] == 0) && (i_xyz[2] == 0)) {
+    is_valid = true;
+    dir[0] = std::copysign(1.0, i_xyz[0]);
+  } else if ((i_xyz[0] == 0) && (std::abs(i_xyz[1]) == 1) && (i_xyz[2] == 0)) {
+    is_valid = true;
+    dir[1] = std::copysign(1.0, i_xyz[1]);
+  } else if ((i_xyz[0] == 0) && (i_xyz[1] == 0) && (std::abs(i_xyz[2]) == 1)) {
+    is_valid = true;
+    dir[2] = std::copysign(1.0, i_xyz[2]);
+  }
+  return dir;
+}
+
+//==============================================================================
+
 int32_t& RectLattice::offset(int map, const array<int, 3>& i_xyz)
 {
   return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map +
@@ -986,6 +1006,91 @@ Position HexLattice::get_local_position(
 
 //==============================================================================
 
+Direction HexLattice::get_normal(
+  const array<int, 3>& i_xyz, bool& is_valid) const
+{
+  // Short description of the direction vectors used here.  The beta, gamma, and
+  // delta vectors point towards the flat sides of each hexagonal tile.
+  // Y - orientation:
+  //   basis0 = (1, 0)
+  //   basis1 = (-1/sqrt(3), 1)   = +120 degrees from basis0
+  //   beta   = (sqrt(3)/2, 1/2)  = +30 degrees from basis0
+  //   gamma  = (sqrt(3)/2, -1/2) = -60 degrees from beta
+  //   delta  = (0, 1)            = +60 degrees from beta
+  // X - orientation:
+  //   basis0 = (1/sqrt(3), -1)
+  //   basis1 = (0, 1)            = +120 degrees from basis0
+  //   beta   = (1, 0)            = +30 degrees from basis0
+  //   gamma  = (1/2, -sqrt(3)/2) = -60 degrees from beta
+  //   delta  = (1/2, sqrt(3)/2)  = +60 degrees from beta
+
+  is_valid = false;
+  Direction dir = {0.0, 0.0, 0.0};
+  if ((i_xyz[0] == 0) && (i_xyz[1] == 0) && (std::abs(i_xyz[2]) == 1)) {
+    is_valid = true;
+    dir[2] = std::copysign(1.0, i_xyz[2]);
+  } else if ((i_xyz[2] == 0) &&
+             std::max({std::abs(i_xyz[0]), std::abs(i_xyz[1]),
+               std::abs(i_xyz[0] + i_xyz[1])}) == 1) {
+    is_valid = true;
+    // beta direction
+    if ((i_xyz[0] == 1) && (i_xyz[1] == 0)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.5 * std::sqrt(3.0);
+        dir[1] = 0.5;
+      } else {
+        dir[0] = 1.0;
+        dir[1] = 0.0;
+      }
+    } else if ((i_xyz[0] == -1) && (i_xyz[1] == 0)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = -0.5 * std::sqrt(3.0);
+        dir[1] = -0.5;
+      } else {
+        dir[0] = -1.0;
+        dir[1] = 0.0;
+      }
+      // gamma direction
+    } else if ((i_xyz[0] == 1) && (i_xyz[1] == -1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.5 * std::sqrt(3.0);
+        dir[1] = -0.5;
+      } else {
+        dir[0] = 0.5;
+        dir[1] = -0.5 * std::sqrt(3.0);
+      }
+    } else if ((i_xyz[0] == -1) && (i_xyz[1] == 1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = -0.5 * std::sqrt(3.0);
+        dir[1] = 0.5;
+      } else {
+        dir[0] = -0.5;
+        dir[1] = 0.5 * std::sqrt(3.0);
+      }
+      // delta direction
+    } else if ((i_xyz[0] == 0) && (i_xyz[1] == 1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.0;
+        dir[1] = 1.0;
+      } else {
+        dir[0] = 0.5;
+        dir[1] = 0.5 * std::sqrt(3.0);
+      }
+    } else if ((i_xyz[0] == 0) && (i_xyz[1] == -1)) {
+      if (orientation_ == Orientation::y) {
+        dir[0] = 0.0;
+        dir[1] = -1.0;
+      } else {
+        dir[0] = -0.5;
+        dir[1] = -0.5 * std::sqrt(3.0);
+      }
+    }
+  }
+  return dir;
+}
+
+//==============================================================================
+
 bool HexLattice::is_valid_index(int indx) const
 {
   int nx {2 * n_rings_ - 1};