Fix M_PI being unavailable where _USE_MATH_DEFINES comes too late (#4023)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

.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/pull_request_template.md

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

.github/workflows/ci.yml

@@ -1,4 +1,4 @@
-name: CI
+name: Tests and Coverage
 
 on:
   # allows us to run workflows manually
@@ -21,49 +21,64 @@ env:
   GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
 
 jobs:
+  filter-changes:
+    runs-on: ubuntu-latest
+    outputs:
+      source_changed: ${{ steps.filter.outputs.source_changed }}
+    steps:
+      - name: Check out the repository
+        uses: actions/checkout@v6
+      - name: Examine changed files
+        id: filter
+        uses: dorny/paths-filter@v4
+        with:
+          filters: |
+            source_changed:
+              - '!docs/**'
+              - '!**/*.md'
+          predicate-quantifier: 'every'
   main:
+    needs: filter-changes
+    if: ${{ needs.filter-changes.outputs.source_changed == 'true' }}
     runs-on: ubuntu-22.04
     strategy:
       matrix:
-        python-version: ["3.11"]
+        python-version: ["3.12"]
         mpi: [n, y]
         omp: [n, y]
         dagmc: [n]
         libmesh: [n]
         event: [n]
-        vectfit: [n]
 
         include:
-          - python-version: "3.12"
-            omp: n
-            mpi: n
           - python-version: "3.13"
             omp: n
             mpi: n
+          - python-version: "3.14"
+            omp: n
+            mpi: n
+          - python-version: "3.14t"
+            omp: n
+            mpi: n
           - dagmc: y
-            python-version: "3.11"
+            python-version: "3.12"
             mpi: y
             omp: y
           - libmesh: y
-            python-version: "3.11"
+            python-version: "3.12"
             mpi: y
             omp: y
           - libmesh: y
-            python-version: "3.11"
+            python-version: "3.12"
             mpi: n
             omp: y
           - event: y
-            python-version: "3.11"
+            python-version: "3.12"
             omp: y
             mpi: n
-          - vectfit: y
-            python-version: "3.11"
-            omp: n
-            mpi: y
     name: "Python ${{ matrix.python-version }} (omp=${{ matrix.omp }},
       mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }},
-      libmesh=${{ matrix.libmesh }}, event=${{ matrix.event }}
-      vectfit=${{ matrix.vectfit }})"
+      libmesh=${{ matrix.libmesh }}, event=${{ matrix.event }}"
 
     env:
       MPI: ${{ matrix.mpi }}
@@ -71,7 +86,6 @@ jobs:
       OMP: ${{ matrix.omp }}
       DAGMC: ${{ matrix.dagmc }}
       EVENT: ${{ matrix.event }}
-      VECTFIT: ${{ matrix.vectfit }}
       LIBMESH: ${{ matrix.libmesh }}
       NPY_DISABLE_CPU_FEATURES: "AVX512F AVX512_SKX"
       OPENBLAS_NUM_THREADS: 1
@@ -88,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 }}
 
@@ -132,11 +146,6 @@ jobs:
           sudo update-alternatives --set mpirun /usr/bin/mpirun.mpich
           sudo update-alternatives --set mpi-x86_64-linux-gnu /usr/include/x86_64-linux-gnu/mpich
 
-      - name: Optional apt dependencies for vectfit
-        shell: bash
-        if: ${{ matrix.vectfit == 'y' }}
-        run: sudo apt install -y libblas-dev liblapack-dev
-
       - name: install
         shell: bash
         run: |
@@ -149,12 +158,12 @@ jobs:
           openmc -v
 
       - name: cache-xs
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: |
             ~/nndc_hdf5
             ~/endf-b-vii.1
-          key: ${{ runner.os }}-build-xs-cache
+          key: ${{ runner.os }}-build-xs-cache-${{ hashFiles(format('{0}/tools/ci/download-xs.sh', github.workspace)) }}
 
       - name: before
         shell: bash
@@ -168,7 +177,7 @@ jobs:
 
       - name: Setup tmate debug session
         continue-on-error: true
-        if: ${{ contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
+        if: ${{ failure() && contains(env.COMMIT_MESSAGE, '[gha-debug]') }}
         uses: mxschmitt/action-tmate@v3
         timeout-minutes: 10
 
@@ -186,6 +195,7 @@ jobs:
             --gcov-ignore-errors source_not_found \
             --gcov-ignore-errors output_error \
             --gcov-ignore-parse-errors suspicious_hits.warn \
+            --merge-mode-functions=separate \
             --print-summary \
             --lcov -o coverage-cpp.lcov || true
 
@@ -203,13 +213,36 @@ jobs:
           parallel: true
           flag-name: C++ and Python
           path-to-lcov: coverage.lcov
+          fail-on-error: false
 
-  finish:
-    needs: main
+  coverage:
+    needs: [filter-changes, main]
+    if: ${{ always() }}
     runs-on: ubuntu-latest
     steps:
       - name: Coveralls Finished
+        if: ${{ needs.filter-changes.outputs.source_changed == 'true' }}
         uses: coverallsapp/github-action@v2
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           parallel-finished: true
+          fail-on-error: false
+
+  ci-pass:
+    needs: [filter-changes, main, coverage]
+    name: Check CI status
+    if: ${{ always() }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Check CI status
+        run: |
+          if [[ "${{ needs.filter-changes.outputs.source_changed }}" == "false" ]]; then
+            echo "Documentation-only change - CI skipped successfully"
+            exit 0
+          fi
+          if [[ "${{ needs.main.result }}" == "success" && "${{ needs.coverage.result }}" == "success" ]]; then
+            echo "CI passed"
+            exit 0
+          fi
+          echo "CI failed"
+          exit 1

.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

@@ -5,6 +5,12 @@ on:
   workflow_dispatch:
 
   pull_request:
+    types:
+      - opened
+      - synchronize
+      - reopened
+      - labeled
+      - unlabeled
     branches:
       - develop
       - master
@@ -12,8 +18,11 @@ on:
 jobs:
   cpp-linter:
     runs-on: ubuntu-latest
+    permissions:
+      contents: read
+      pull-requests: write
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
       - uses: cpp-linter/cpp-linter-action@v2
         id: linter
         env:
@@ -22,11 +31,30 @@ jobs:
           style: file
           files-changed-only: true
           tidy-checks: '-*'
-          version: '15' # clang-format version
+          version: '18' # clang-format version
+          format-review: ${{ github.event_name == 'pull_request' && contains(github.event.pull_request.labels.*.name, 'cpp-format-suggest') }}
+          passive-reviews: ${{ github.event_name == 'pull_request' && contains(github.event.pull_request.labels.*.name, 'cpp-format-suggest') }}
           file-annotations: true
           step-summary: true
           extensions: 'cpp,h'
 
+      - name: Comment with suggestion instructions
+        if: steps.linter.outputs.checks-failed > 0 && !contains(github.event.pull_request.labels.*.name, 'cpp-format-suggest')
+        uses: actions/github-script@v7
+        with:
+          script: |
+            const {owner, repo} = context.repo;
+            const issue_number = context.payload.pull_request.number;
+            await github.rest.issues.createComment({
+              owner,
+              repo,
+              issue_number,
+              body: "C++ formatting checks failed. Add the `cpp-format-suggest` label to this PR for inline formatting suggestions on the next run."
+            });
+
       - name: Failure Check
         if: steps.linter.outputs.checks-failed > 0
-        run: echo "Some files failed the formatting check! See job summary and file annotations for more info" && exit 1
+        run: |
+          echo "Some files failed the formatting check."
+          echo "See job summary and file annotations for details."
+          exit 1

.gitignore

@@ -25,12 +25,13 @@ examples/**/*.xml
 
 # Documentation builds
 docs/build
+docs/doxygen/xml
 docs/source/_images/*.pdf
 docs/source/_images/*.aux
 docs/source/pythonapi/generated/
 
 # Source build
-build
+build*/
 
 # build from src/utils/setup.py
 src/utils/build
@@ -104,5 +105,8 @@ CMakeSettings.json
 # Visual Studio Code configuration files
 .vscode/
 
+# Claude Code agent tools (cached/generated artifacts)
+.claude/cache/
+
 # Python pickle files
 *.pkl

.gitmodules

@@ -1,12 +1,6 @@
 [submodule "vendor/pugixml"]
 	path = vendor/pugixml
 	url = https://github.com/zeux/pugixml.git
-[submodule "vendor/xtensor"]
-	path = vendor/xtensor
-	url = https://github.com/xtensor-stack/xtensor.git
-[submodule "vendor/xtl"]
-	path = vendor/xtl
-	url = https://github.com/xtensor-stack/xtl.git
 [submodule "vendor/fmt"]
 	path = vendor/fmt
 	url = https://github.com/fmtlib/fmt.git

.mcp.json

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

.readthedocs.yaml

@@ -7,9 +7,14 @@ build:
   jobs:
     post_checkout:
       - git fetch --unshallow || true
+      - cd docs/doxygen && doxygen && cd -
+
 sphinx:
   configuration: docs/source/conf.py
 
+formats:
+  - pdf
+
 python:
   install:
     - method: pip

AGENTS.md

@@ -0,0 +1,348 @@
+# OpenMC AI Coding Agent Instructions
+
+## Project Overview
+
+OpenMC is a Monte Carlo particle transport code for simulating nuclear reactors,
+fusion devices, or other systems with neutron/photon radiation. It's a hybrid
+C++17/Python codebase where:
+- **C++ core** (`src/`, `include/openmc/`) handles the computationally intensive transport simulation
+- **Python API** (`openmc/`) provides user-facing model building, post-processing, and depletion capabilities
+- **C API bindings** (`openmc/lib/`) wrap the C++ library via ctypes for runtime control
+
+## Architecture & Key Components
+
+### C++ Component Structure
+- **Global vectors of unique_ptrs**: Core objects like `model::cells`, `model::universes`, `nuclides` are stored as `vector<unique_ptr<T>>` in nested namespaces (`openmc::model`, `openmc::simulation`, `openmc::settings`, `openmc::data`)
+- **Custom container types**: OpenMC provides its own `vector`, `array`, `unique_ptr`, and `make_unique` in the `openmc::` namespace (defined in `vector.h`, `array.h`, `memory.h`). These are currently typedefs to `std::` equivalents but may become custom implementations for accelerator support. Always use `openmc::vector`, not `std::vector`.
+- **Geometry systems**:
+  - **CSG (default)**: Arbitrarily complex Constructive Solid Geometry using `Surface`, `Region`, `Cell`, `Universe`, `Lattice`
+  - **DAGMC**: CAD-based geometry via Direct Accelerated Geometry Monte Carlo (optional, requires `OPENMC_USE_DAGMC`)
+  - **Unstructured mesh**: libMesh-based geometry (optional, requires `OPENMC_USE_LIBMESH`)
+- **Particle tracking**: `Particle` class with `GeometryState` manages particle transport through geometry
+- **Tallies**: Score quantities during simulation via `Filter` and `Tally` objects
+- **Random ray solver**: Alternative deterministic method in `src/random_ray/`
+- **Optional features**: DAGMC (CAD geometry), libMesh (unstructured mesh), MPI, all controlled by `#ifdef OPENMC_MPI`, etc.
+
+### Python Component Structure
+- **ID management**: All geometry objects (Cell, Surface, Material, etc.) inherit from `IDManagerMixin` which auto-assigns unique integer IDs and tracks them via class-level `used_ids` and `next_id`
+- **Input validation**: Extensive use of `openmc.checkvalue` module functions (`check_type`, `check_value`, `check_length`) for all setters
+- **XML I/O**: Most classes implement `to_xml_element()` and `from_xml_element()` for serialization to OpenMC's XML input format
+- **HDF5 output**: Post-simulation data in statepoint files read via `openmc.StatePoint`
+- **Depletion**: `openmc.deplete` implements burnup via operator-splitting with various integrators (Predictor, CECM, etc.)
+- **Nuclear Data**: `openmc.data` provides programmatic access to nuclear data files (ENDF, ACE, HDF5)
+
+## Git Branching Workflow
+
+OpenMC uses a git flow branching model with two primary branches:
+
+- **`develop` branch**: The main development branch where all ongoing development takes place. This is the **primary branch against which pull requests are submitted and merged**. This branch is not guaranteed to be stable and may contain work-in-progress features.
+- **`master` branch**: The stable release branch containing the latest stable release of OpenMC. This branch only receives merges from `develop` when the development team decides a release should occur.
+
+### Instructions for Code Review
+
+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
+
+1. Create a feature/bugfix branch off `develop`
+2. Make changes and commit to the feature branch
+3. Open a pull request to merge the feature branch into `develop`
+4. A committer reviews and merges the PR into `develop`
+
+## Critical Build & Test Workflows
+
+### Build Dependencies
+- **C++17 compiler**: GCC, Clang, or Intel
+- **CMake** (3.16+): Required for configuring and building the C++ library
+- **HDF5**: Required for cross section data and output file formats
+- **libpng**: Used for generating visualization when OpenMC is run in plotting mode
+
+Without CMake and HDF5, OpenMC cannot be compiled.
+
+### Building the C++ Library
+```bash
+# Configure with CMake (from build/ directory)
+cmake .. -DOPENMC_USE_MPI=ON -DOPENMC_USE_OPENMP=ON -DCMAKE_BUILD_TYPE=RelWithDebInfo
+
+# Available CMake options (all default OFF except OPENMC_USE_OPENMP and OPENMC_BUILD_TESTS):
+# -DOPENMC_USE_OPENMP=ON/OFF        # OpenMP parallelism
+# -DOPENMC_USE_MPI=ON/OFF           # MPI support
+# -DOPENMC_USE_DAGMC=ON/OFF         # CAD geometry support
+# -DOPENMC_USE_LIBMESH=ON/OFF       # Unstructured mesh
+# -DOPENMC_ENABLE_PROFILE=ON/OFF    # Profiling flags
+# -DOPENMC_ENABLE_COVERAGE=ON/OFF   # Coverage analysis
+
+# Build
+make -j
+
+# C++ unit tests (uses Catch2)
+ctest
+```
+
+### Python Development
+```bash
+# Install in development mode (requires building C++ library first)
+pip install -e .
+
+# Python tests (uses pytest)
+pytest tests/unit_tests/           # Fast unit tests
+pytest tests/regression_tests/     # Full regression suite (requires nuclear data)
+```
+
+### Nuclear Data Setup (CRITICAL for Running OpenMC)
+Most tests require the NNDC HDF5 nuclear cross-section library.
+
+**Important**: Check if `OPENMC_CROSS_SECTIONS` is already set in the user's
+environment before downloading, as many users already have nuclear data
+installed. Though do note that if this variable is present that it may point to
+different cross section data and that the NNDC data is required for tests to
+pass.
+
+**If not already configured, download and setup:**
+```bash
+# Download NNDC HDF5 cross section library (~800 MB compressed)
+wget -q -O - https://anl.box.com/shared/static/teaup95cqv8s9nn56hfn7ku8mmelr95p.xz | tar -C $HOME -xJ
+
+# Set environment variable (add to ~/.bashrc or ~/.zshrc for persistence)
+export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
+```
+
+**Alternative**: Use the provided download script (checks if data exists before downloading):
+```bash
+bash tools/ci/download-xs.sh  # Downloads both NNDC HDF5 and ENDF/B-VII.1 data
+```
+
+Without this data, regression tests will fail with "No cross_sections.xml file
+found" errors, or, in the case that alternative cross section data is configured
+the tests will execute but will not pass. The `cross_sections.xml` file is an
+index listing paths to individual HDF5 nuclear data files for each nuclide.
+
+## Testing Expectations
+
+### Environment Requirements
+
+  - **Data**: As described above, OpenMC's test suite requires OpenMC to be configured with NNDC data.
+  - **OpenMP Settings**: OpenMC's tests may fail is more than two OpenMP threads are used. The environment variable `OMP_NUM_THREADS=2` should be set to avoid sporadic test failures.
+  - **Executable configuration**: The OpenMC executable should compiled with debug symbols enabled.
+
+### C++ Tests
+Located in `tests/cpp_unit_tests/`, use Catch2 framework. Run via `ctest` after building with `-DOPENMC_BUILD_TESTS=ON`.
+
+### Python Unit Tests
+Located in `tests/unit_tests/`, these are fast, standalone tests that verify Python API functionality without running full simulations. Use standard pytest patterns:
+
+**Categories**:
+- **API validation**: Test object creation, property setters/getters, XML serialization (e.g., `test_material.py`, `test_cell.py`, `test_source.py`)
+- **Data processing**: Test nuclear data handling, cross sections, depletion chains (e.g., `test_data_neutron.py`, `test_deplete_chain.py`)
+- **Library bindings**: Test `openmc.lib` ctypes interface with `model.init_lib()`/`model.finalize_lib()` (e.g., `test_lib.py`)
+- **Geometry operations**: Test bounding boxes, containment, lattice generation (e.g., `test_bounding_box.py`, `test_lattice.py`)
+
+**Common patterns**:
+- Use fixtures from `tests/unit_tests/conftest.py` (e.g., `uo2`, `water`, `sphere_model`)
+- Test invalid inputs with `pytest.raises(ValueError)` or `pytest.raises(TypeError)`
+- Use `run_in_tmpdir` fixture for tests that create files
+- Tests with `openmc.lib` require calling `model.init_lib()` in try/finally with `model.finalize_lib()`
+
+**Example**:
+```python
+def test_material_properties():
+    m = openmc.Material()
+    m.add_nuclide('U235', 1.0)
+    assert 'U235' in m.nuclides
+
+    with pytest.raises(TypeError):
+        m.add_nuclide('H1', '1.0')  # Invalid type
+```
+
+Unit tests should be fast. For tests requiring simulation output, use regression tests instead.
+
+### Python Regression Tests
+Regression tests compare OpenMC output against reference data. **Prefer using existing models from `openmc.examples` or those found in tests/unit_tests/conftest.py** (like `pwr_pin_cell()`, `pwr_assembly()`, `slab_mg()`) rather than building from scratch.
+
+**Test Harness Types** (in `tests/testing_harness.py`):
+- **PyAPITestHarness**: Standard harness for Python API tests. Compares `inputs_true.dat` (XML hash) and `results_true.dat` (statepoint k-eff and tally values). Requires `model.xml` generation.
+- **HashedPyAPITestHarness**: Like PyAPITestHarness but hashes the results for compact comparison
+- **TolerantPyAPITestHarness**: For tests with floating-point non-associativity (e.g., random ray solver with single precision). Uses relative tolerance comparisons.
+- **WeightWindowPyAPITestHarness**: Compares weight window bounds from `weight_windows.h5`
+- **CollisionTrackTestHarness**: Compares collision track data from `collision_track.h5` against `collision_track_true.h5`
+- **TestHarness**: Base harness for XML-based tests (no Python model building)
+- **PlotTestHarness**: Compares plot output files (PNG or voxel HDF5)
+- **CMFDTestHarness**: Specialized for CMFD acceleration tests
+- **ParticleRestartTestHarness**: Tests particle restart functionality
+
+Almost all cases use either `PyAPITestHarness` or `HashedPyAPITestHarness`
+
+**Example Test**:
+```python
+from openmc.examples import pwr_pin_cell
+from tests.testing_harness import PyAPITestHarness
+
+def test_my_feature():
+    model = pwr_pin_cell()
+    model.settings.particles = 1000  # Modify to exercise feature
+    harness = PyAPITestHarness('statepoint.10.h5', model)
+    harness.main()
+```
+
+**Workflow**: Create `test.py` and `__init__.py` in `tests/regression_tests/my_test/`, run `pytest --update` to generate reference files (`inputs_true.dat`, `results_true.dat`, etc.), then verify with `pytest` without `--update`. Test results should be generated with `-DOPENMC_ENABLE_STRICT_FP=on` to ensure reproducibility across platforms and optimization levels.
+
+**Critical**: When modifying OpenMC code, regenerate affected test references with `pytest --update` and commit updated reference files.
+
+### Test Configuration
+
+`pytest.ini` sets: `python_files = test*.py`, `python_classes = NoThanks` (disables class-based test collection).
+
+### Testing Options
+
+For builds of OpenMC with MPI enabled, the `--mpi` flag should be passed to the test suite to ensure that appropriate tests are executed using two MPI processes.
+
+The entire test suite can be executed with OpenMC running in event-based mode (instead of the default history-based mode) by providing the `--event` flag to the `pytest` command.
+
+## Cross-Language Boundaries
+
+The C API (defined in `include/openmc/capi.h`) exposes C++ functionality to Python via ctypes bindings in `openmc/lib/`. Example:
+```cpp
+// C++ API in capi.h
+extern "C" int openmc_run();
+
+// Python binding in openmc/lib/core.py
+_dll.openmc_run.restype = c_int
+def run():
+    _dll.openmc_run()
+```
+
+When modifying C++ public APIs, update corresponding ctypes signatures in `openmc/lib/*.py`.
+
+## Code Style & Conventions
+
+### C++ Style (enforced by .clang-format)
+  OpenMC generally tries to follow C++ core guidelines where possible
+  (https://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines) and follow
+  modern C++ practices (e.g. RAII) whenever possible.
+
+- **Naming**:
+  - Classes: `CamelCase` (e.g., `HexLattice`)
+  - Functions/methods: `snake_case` (e.g., `get_indices`)
+  - Variables: `snake_case` with trailing underscore for class members (e.g., `n_particles_`, `energy_`)
+  - Constants: `UPPER_SNAKE_CASE` (e.g., `SQRT_PI`)
+- **Namespaces**: All code in `openmc::` namespace, global state in sub-namespaces
+- **Include order**: Related header first, then C/C++ stdlib, third-party libs, local headers
+- **Comments**: C++-style (`//`) only, never C-style (`/* */`)
+- **Standard**: C++17 features allowed
+- **Formatting**: Run `clang-format` (version 18) before committing; install via `tools/dev/install-commit-hooks.sh`
+
+### Python Style
+- **PEP8** compliant
+- **Docstrings**: numpydoc format for all public functions/methods
+- **Type hints**: Use sparingly, primarily for complex signatures
+- **Path handling**: Use `pathlib.Path` for filesystem operations, accept `str | os.PathLike` in function arguments
+- **Dependencies**: Core dependencies only (numpy, scipy, h5py, pandas, matplotlib, lxml, ipython, uncertainties, endf). Other packages must be optional
+- **Python version**: Minimum 3.11 (as of Nov 2025)
+
+### ID Management Pattern (Python)
+When creating geometry objects, IDs can be auto-assigned or explicit:
+```python
+# Auto-assigned ID
+cell = openmc.Cell()  # Gets next available ID
+
+# Explicit ID
+cell = openmc.Cell(id=10)  # Warning if ID already used
+
+# Reset all IDs (useful in test fixtures)
+openmc.reset_auto_ids()
+```
+
+### Input Validation Pattern (Python)
+All setters use checkvalue functions:
+```python
+import openmc.checkvalue as cv
+
+@property
+def temperature(self):
+    return self._temperature
+
+@temperature.setter
+def temperature(self, temp):
+    cv.check_type('temperature', temp, Real)
+    cv.check_greater_than('temperature', temp, 0.0)
+    self._temperature = temp
+```
+
+### Working with HDF5 Files
+C++ uses custom HDF5 wrappers in `src/hdf5_interface.cpp`. Python uses h5py directly. Statepoint format version is `VERSION_STATEPOINT` in `include/openmc/constants.h`.
+
+### Conditional Compilation
+Check for optional features:
+```cpp
+#ifdef OPENMC_MPI
+  // MPI-specific code
+#endif
+
+#ifdef OPENMC_DAGMC
+  // DAGMC-specific code
+#endif
+```
+
+## Documentation
+
+- **User docs**: Sphinx documentation in `docs/source/` hosted at https://docs.openmc.org
+- **C++ docs**: Doxygen-style comments with `\brief`, `\param` tags
+- **Python docs**: numpydoc format docstrings
+
+## Common Pitfalls
+
+1. **Forgetting nuclear data**: Tests fail without `OPENMC_CROSS_SECTIONS` environment variable
+2. **ID conflicts**: Python objects with duplicate IDs trigger `IDWarning`, use `reset_auto_ids()` between tests
+3. **MPI builds**: Code must work with and without MPI; use `#ifdef OPENMC_MPI` guards
+4. **Path handling**: Use `pathlib.Path` in new Python code, not `os.path`
+5. **Clang-format version**: CI uses version 18; other versions may produce different formatting

CITATION.cff

@@ -1,9 +1,43 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+title: OpenMC
+authors:
+- family-names: Romano
+  given-names: Paul K.
+  orcid: "https://orcid.org/0000-0002-1147-045X"
+- family-names: Shriwise
+  given-names: Patrick C.
+  orcid: "https://orcid.org/0000-0002-3979-7665"
+- family-names: Shimwell
+  given-names: Jonathan
+  orcid: "https://orcid.org/0000-0001-6909-0946"
+- family-names: Harper
+  given-names: Sterling
+- family-names: Boyd
+  given-names: Will
+- family-names: Nelson
+  given-names: Adam G.
+  orcid: "https://orcid.org/0000-0002-3614-0676"
+- family-names: Tramm
+  given-names: John R.
+  orcid: "https://orcid.org/0000-0002-5397-4402"
+- family-names: Ridley
+  given-names: Gavin
+  orcid: "https://orcid.org/0000-0003-1635-8042"
+- family-names: Johnson
+  given-names: Andrew
+  orcid: "https://orcid.org/0000-0003-2125-8775"
+- family-names: Peterson
+  given-names: Ethan E.
+  orcid: "https://orcid.org/0000-0002-5694-7194"
+- family-names: Herman
+  given-names: Bryan R.
 preferred-citation:
   authors:
   - family-names: Romano
     given-names: Paul K.
     orcid: "https://orcid.org/0000-0002-1147-045X"
-  - final-names: Horelik
+  - family-names: Horelik
     given-names: Nicholas E.
   - family-names: Herman
     given-names: Bryan R.

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)
@@ -38,6 +43,7 @@ option(OPENMC_USE_LIBMESH     "Enable support for libMesh unstructured mesh tall
 option(OPENMC_USE_MPI         "Enable MPI"                                           OFF)
 option(OPENMC_USE_UWUW        "Enable UWUW"                                          OFF)
 option(OPENMC_FORCE_VENDORED_LIBS "Explicitly use submodules defined in 'vendor'"    OFF)
+option(OPENMC_ENABLE_STRICT_FP "Enable strict FP flags to improve test portability"  OFF)
 
 message(STATUS "OPENMC_USE_OPENMP ${OPENMC_USE_OPENMP}")
 message(STATUS "OPENMC_BUILD_TESTS ${OPENMC_BUILD_TESTS}")
@@ -48,6 +54,7 @@ message(STATUS "OPENMC_USE_LIBMESH ${OPENMC_USE_LIBMESH}")
 message(STATUS "OPENMC_USE_MPI ${OPENMC_USE_MPI}")
 message(STATUS "OPENMC_USE_UWUW ${OPENMC_USE_UWUW}")
 message(STATUS "OPENMC_FORCE_VENDORED_LIBS ${OPENMC_FORCE_VENDORED_LIBS}")
+message(STATUS "OPENMC_ENABLE_STRICT_FP ${OPENMC_ENABLE_STRICT_FP}")
 
 # Warnings for deprecated options
 foreach(OLD_OPT IN ITEMS "openmp" "profile" "coverage" "dagmc" "libmesh")
@@ -89,6 +96,19 @@ if(NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE RelWithDebInfo CACHE STRING "Choose the type of build" FORCE)
 endif()
 
+#===============================================================================
+# When STRICT_FP is enabled, remove NDEBUG from RelWithDebInfo flags so that
+# assert() remains active. CMake normally adds -DNDEBUG for both Release and
+# RelWithDebInfo, which disables C/C++ assert() statements.
+#===============================================================================
+
+if(OPENMC_ENABLE_STRICT_FP)
+  foreach(FLAG_VAR CMAKE_CXX_FLAGS_RELWITHDEBINFO CMAKE_C_FLAGS_RELWITHDEBINFO)
+    string(REPLACE "-DNDEBUG" "" ${FLAG_VAR} "${${FLAG_VAR}}")
+    string(REPLACE "/DNDEBUG" "" ${FLAG_VAR} "${${FLAG_VAR}}")
+  endforeach()
+endif()
+
 #===============================================================================
 # OpenMP for shared-memory parallelism (and GPU support some day!)
 #===============================================================================
@@ -193,6 +213,26 @@ endif()
 # Set compile/link flags based on which compiler is being used
 #===============================================================================
 
+# When OPENMC_ENABLE_STRICT_FP is enabled, disable compiler optimizations that change
+# floating-point results relative to -O0, improving cross-platform and
+# cross-optimization-level reproducibility for regression testing:
+#   -ffp-contract=off    Prevents FMA contraction (fused multiply-add changes rounding)
+#   -fno-builtin         Prevents replacing math function calls (pow, exp, log, etc.)
+#                        with builtin versions that may differ from libm
+# By default (OFF), the compiler is free to use all optimizations for best
+# performance.
+if(OPENMC_ENABLE_STRICT_FP)
+  include(CheckCXXCompilerFlag)
+  check_cxx_compiler_flag(-ffp-contract=off SUPPORTS_FP_CONTRACT_OFF)
+  if(SUPPORTS_FP_CONTRACT_OFF)
+    list(APPEND cxxflags -ffp-contract=off)
+  endif()
+  check_cxx_compiler_flag(-fno-builtin SUPPORTS_NO_BUILTIN)
+  if(SUPPORTS_NO_BUILTIN)
+    list(APPEND cxxflags -fno-builtin)
+  endif()
+endif()
+
 # Skip for Visual Studio which has its own configurations through GUI
 if(NOT MSVC)
 
@@ -266,23 +306,6 @@ else()
   endif()
 endif()
 
-#===============================================================================
-# xtensor header-only library
-#===============================================================================
-
-if(OPENMC_FORCE_VENDORED_LIBS)
-  add_subdirectory(vendor/xtl)
-  set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl)
-  add_subdirectory(vendor/xtensor)
-else()
-  find_package_write_status(xtensor)
-  if (NOT xtensor_FOUND)
-    add_subdirectory(vendor/xtl)
-    set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl)
-    add_subdirectory(vendor/xtensor)
-  endif()
-endif()
-
 #===============================================================================
 # Catch2 library
 #===============================================================================
@@ -332,6 +355,7 @@ endif()
 #===============================================================================
 
 list(APPEND libopenmc_SOURCES
+  src/atomic_mass.cpp
   src/bank.cpp
   src/boundary_condition.cpp
   src/bremsstrahlung.cpp
@@ -372,6 +396,7 @@ list(APPEND libopenmc_SOURCES
   src/particle.cpp
   src/particle_data.cpp
   src/particle_restart.cpp
+  src/particle_type.cpp
   src/photon.cpp
   src/physics.cpp
   src/physics_common.cpp
@@ -387,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
@@ -425,7 +451,9 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_musurface.cpp
   src/tallies/filter_parent_nuclide.cpp
   src/tallies/filter_particle.cpp
+  src/tallies/filter_particle_production.cpp
   src/tallies/filter_polar.cpp
+  src/tallies/filter_reaction.cpp
   src/tallies/filter_sph_harm.cpp
   src/tallies/filter_sptl_legendre.cpp
   src/tallies/filter_surface.cpp
@@ -495,7 +523,7 @@ endif()
 # target_link_libraries treats any arguments starting with - but not -l as
 # linker flags. Thus, we can pass both linker flags and libraries together.
 target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}
-                      xtensor fmt::fmt ${CMAKE_DL_LIBS})
+                      fmt::fmt ${CMAKE_DL_LIBS})
 
 if(TARGET pugixml::pugixml)
   target_link_libraries(libopenmc pugixml::pugixml)
@@ -504,7 +532,7 @@ else()
 endif()
 
 if(OPENMC_USE_DAGMC)
-  target_compile_definitions(libopenmc PRIVATE OPENMC_DAGMC_ENABLED)
+  target_compile_definitions(libopenmc PUBLIC OPENMC_DAGMC_ENABLED)
   target_link_libraries(libopenmc dagmc-shared)
 
   if(OPENMC_USE_UWUW)
@@ -552,6 +580,9 @@ endif()
 if (OPENMC_ENABLE_COVERAGE)
   target_compile_definitions(libopenmc PRIVATE COVERAGEBUILD)
 endif()
+if (OPENMC_ENABLE_STRICT_FP)
+  target_compile_definitions(libopenmc PRIVATE OPENMC_ENABLE_STRICT_FP)
+endif()
 
 #===============================================================================
 # openmc executable
@@ -580,9 +611,7 @@ add_custom_command(TARGET libopenmc POST_BUILD
 #===============================================================================
 # Install executable, scripts, manpage, license
 #===============================================================================
-
-configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY)
-configure_file(cmake/OpenMCConfigVersion.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake" @ONLY)
+include(CMakePackageConfigHelpers)
 
 set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC)
 install(TARGETS openmc libopenmc
@@ -596,10 +625,24 @@ install(EXPORT openmc-targets
   NAMESPACE OpenMC::
   DESTINATION ${INSTALL_CONFIGDIR})
 
+configure_package_config_file(
+    "cmake/OpenMCConfig.cmake.in"
+    "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
+    INSTALL_DESTINATION ${INSTALL_CONFIGDIR}
+)
+
+write_basic_package_version_file(
+    "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
+    VERSION ${OPENMC_VERSION}
+    COMPATIBILITY AnyNewerVersion
+)
+
 install(FILES
-  "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
-  "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
-  DESTINATION ${INSTALL_CONFIGDIR})
+    "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake"
+    "${CMAKE_BINARY_DIR}/${CMAKE_FILES_DIRECTORY}/OpenMCConfigVersion.cmake"
+    DESTINATION "${INSTALL_CONFIGDIR}"
+)
+
 install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1)
 install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright)
 install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

Dockerfile

@@ -33,11 +33,6 @@ ARG build_libmesh
 # Set default value of HOME to /root
 ENV HOME=/root
 
-# Embree variables
-ENV EMBREE_TAG='v4.3.1'
-ENV EMBREE_REPO='https://github.com/embree/embree'
-ENV EMBREE_INSTALL_DIR=$HOME/EMBREE/
-
 # MOAB variables
 ENV MOAB_TAG='5.5.1'
 ENV MOAB_REPO='https://bitbucket.org/fathomteam/moab/'
@@ -58,10 +53,11 @@ 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
-ENV LD_LIBRARY_PATH=${DAGMC_INSTALL_DIR}/lib:$LD_LIBRARY_PATH \
+ENV LD_LIBRARY_PATH=${DAGMC_INSTALL_DIR}/lib:${LD_LIBRARY_PATH:-} \
     OPENMC_ENDF_DATA=/root/endf-b-vii.1 \
     DEBIAN_FRONTEND=noninteractive
 
@@ -71,7 +67,7 @@ RUN apt-get update -y && \
     apt-get install -y \
         python3-pip python-is-python3 wget git build-essential cmake \
         mpich libmpich-dev libhdf5-serial-dev libhdf5-mpich-dev \
-        libpng-dev python3-venv && \
+        libpng-dev libpugixml-dev libfmt-dev catch2 python3-venv && \
     apt-get autoremove
 
 # create virtual enviroment to avoid externally managed environment error
@@ -83,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 \
@@ -94,22 +90,12 @@ RUN cd $HOME \
 
 RUN if [ "$build_dagmc" = "on" ]; then \
         # Install addition packages required for DAGMC
-        apt-get -y install libeigen3-dev libnetcdf-dev libtbb-dev libglfw3-dev \
+        apt-get -y install \
+            libeigen3-dev libnetcdf-dev libtbb-dev libglfw3-dev libembree-dev \
         && pip install --upgrade numpy \
         && pip install --no-cache-dir setuptools cython \
-        # Clone and install EMBREE
-        && mkdir -p $HOME/EMBREE && cd $HOME/EMBREE \
-        && git clone --single-branch -b ${EMBREE_TAG} --depth 1 ${EMBREE_REPO} \
-        && mkdir build && cd build \
-        && cmake ../embree \
-                    -DCMAKE_INSTALL_PREFIX=${EMBREE_INSTALL_DIR} \
-                    -DEMBREE_MAX_ISA=NONE \
-                    -DEMBREE_ISA_SSE42=ON \
-                    -DEMBREE_ISPC_SUPPORT=OFF \
-        && make 2>/dev/null -j${compile_cores} install \
-        && rm -rf ${EMBREE_INSTALL_DIR}/build ${EMBREE_INSTALL_DIR}/embree ; \
         # Clone and install MOAB
-        mkdir -p $HOME/MOAB && cd $HOME/MOAB \
+        && mkdir -p $HOME/MOAB && cd $HOME/MOAB \
         && git clone  --single-branch -b ${MOAB_TAG} --depth 1 ${MOAB_REPO} \
         && mkdir build && cd build \
         && cmake ../moab -DCMAKE_BUILD_TYPE=Release \
@@ -118,6 +104,7 @@ RUN if [ "$build_dagmc" = "on" ]; then \
                       -DBUILD_SHARED_LIBS=OFF \
                       -DENABLE_FORTRAN=OFF \
                       -DENABLE_BLASLAPACK=OFF \
+                      -DENABLE_TESTING=OFF \
         && make 2>/dev/null -j${compile_cores} install \
         && cmake ../moab \
                     -DENABLE_PYMOAB=ON \
@@ -133,7 +120,7 @@ RUN if [ "$build_dagmc" = "on" ]; then \
         && mkdir build && cd build \
         && cmake ../double-down -DCMAKE_INSTALL_PREFIX=${DD_INSTALL_DIR} \
                              -DMOAB_DIR=/usr/local \
-                             -DEMBREE_DIR=${EMBREE_INSTALL_DIR} \
+                             -DEMBREE_DIR=/usr \
         && make 2>/dev/null -j${compile_cores} install \
         && rm -rf ${DD_INSTALL_DIR}/build ${DD_INSTALL_DIR}/double-down ; \
         # Clone and install DAGMC
@@ -147,6 +134,7 @@ RUN if [ "$build_dagmc" = "on" ]; then \
                        -DDOUBLE_DOWN_DIR=${DD_INSTALL_DIR} \
                        -DCMAKE_PREFIX_PATH=${DD_INSTALL_DIR}/lib \
                        -DBUILD_STATIC_LIBS=OFF \
+                       -DBUILD_TESTS=OFF \
         && make 2>/dev/null -j${compile_cores} install \
         && rm -rf ${DAGMC_INSTALL_DIR}/DAGMC ${DAGMC_INSTALL_DIR}/build ; \
     fi

LICENSE

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

cmake/OpenMCConfig.cmake.in

@@ -1,14 +1,18 @@
-get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
+@PACKAGE_INIT@
 
-# Compute the install prefix from this file's location
-get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE)
+include("${CMAKE_CURRENT_LIST_DIR}/OpenMCConfigVersion.cmake")
+include(CMakeFindDependencyMacro)
+
+# Explicitly calculate prefix if it was not generated above
+if(NOT DEFINED PACKAGE_PREFIX_DIR)
+  get_filename_component(PACKAGE_PREFIX_DIR "${CMAKE_CURRENT_LIST_DIR}/../../.." ABSOLUTE)
+endif()
+
+find_dependency(fmt CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR})
+find_dependency(pugixml CONFIG REQUIRED HINTS ${PACKAGE_PREFIX_DIR})
 
-find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
-find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
-find_package(xtl CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
-find_package(xtensor CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
 if(@OPENMC_USE_DAGMC@)
-  find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@)
+  find_dependency(DAGMC REQUIRED HINTS @DAGMC_DIR@)
 endif()
 
 if(@OPENMC_USE_LIBMESH@)
@@ -18,20 +22,24 @@ if(@OPENMC_USE_LIBMESH@)
   pkg_check_modules(LIBMESH REQUIRED @LIBMESH_PC_FILE@>=1.7.0 IMPORTED_TARGET)
 endif()
 
-find_package(PNG)
-
-if(NOT TARGET OpenMC::libopenmc)
-  include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")
+if("@PNG_FOUND@")
+  find_dependency(PNG)
 endif()
 
 if(@OPENMC_USE_MPI@)
-  find_package(MPI REQUIRED)
+  find_dependency(MPI REQUIRED)
 endif()
 
 if(@OPENMC_USE_OPENMP@)
-  find_package(OpenMP REQUIRED)
+  find_dependency(OpenMP REQUIRED)
 endif()
 
 if(@OPENMC_USE_UWUW@ AND NOT ${DAGMC_BUILD_UWUW})
   message(FATAL_ERROR "UWUW is enabled in OpenMC but the DAGMC installation discovered was not configured with UWUW.")
 endif()
+
+include("${CMAKE_CURRENT_LIST_DIR}/OpenMCTargets.cmake")
+
+if(NOT OpenMC_FIND_QUIETLY)
+    message(STATUS "Found OpenMC: ${PACKAGE_VERSION} (found in ${PACKAGE_PREFIX_DIR})")
+endif()

cmake/OpenMCConfigVersion.cmake.in

@@ -1,11 +0,0 @@
-set(PACKAGE_VERSION "@OPENMC_VERSION@")
-
-# Check whether the requested PACKAGE_FIND_VERSION is compatible
-if("${PACKAGE_VERSION}" VERSION_LESS "${PACKAGE_FIND_VERSION}")
-  set(PACKAGE_VERSION_COMPATIBLE FALSE)
-else()
-  set(PACKAGE_VERSION_COMPATIBLE TRUE)
-  if ("${PACKAGE_VERSION}" VERSION_EQUAL "${PACKAGE_FIND_VERSION}")
-    set(PACKAGE_VERSION_EXACT TRUE)
-  endif()
-endif()

docs/Makefile

@@ -45,6 +45,7 @@ help:
 clean:
 	-rm -rf $(BUILDDIR)/*
 	-rm -rf source/pythonapi/generated/
+	-rm -rf doxygen/xml
 
 html:
 	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html

docs/doxygen/Doxyfile

@@ -0,0 +1,13 @@
+# Doxyfile 1.9.1
+
+# This file describes the settings to be used by the documentation system
+# doxygen (www.doxygen.org) for a project.
+
+# Difference with default Doxyfile 1.9.1
+PROJECT_NAME           = OpenMC
+QUIET                  = YES
+WARN_IF_UNDOCUMENTED   = NO
+INPUT                  = ../../include/openmc/capi.h
+GENERATE_HTML          = NO
+GENERATE_LATEX         = NO
+GENERATE_XML           = YES

docs/source/capi/index.rst

@@ -46,43 +46,20 @@ Type Definitions
 Functions
 ---------
 
-.. c:function:: int openmc_calculate_volumes()
+..
+   Once documentation is complete in capi.h, use:
+   .. doxygenfile:: capi.h
+   to populate this documentation without using
+   .. doxygenfunction::
+   for every function.
 
-   Run a stochastic volume calculation
+.. doxygenfunction:: openmc_calculate_volumes
 
-   :return: Return status (negative if an error occurred)
-   :rtype: int
+.. doxygenfunction:: openmc_cell_get_fill
 
-.. c:function:: int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
+.. doxygenfunction:: openmc_cell_get_id
 
-   Get the fill for a cell
-
-   :param int32_t index: Index in the cells array
-   :param int* type: Type of the fill
-   :param int32_t** indices: Array of material indices for cell
-   :param int32_t* n: Length of indices array
-   :return: Return status (negative if an error occurred)
-   :rtype: int
-
-.. c:function:: int openmc_cell_get_id(int32_t index, int32_t* id)
-
-   Get the ID of a cell
-
-   :param int32_t index: Index in the cells array
-   :param int32_t* id: ID of the cell
-   :return: Return status (negative if an error occurred)
-   :rtype: int
-
-.. c:function:: int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T)
-
-   Get the temperature of a cell
-
-   :param int32_t index: Index in the cells array
-   :param int32_t* instance: Which instance of the cell. If a null pointer is passed, the temperature
-                             of the first instance is returned.
-   :param double* T: temperature of the cell
-   :return: Return status (negative if an error occurred)
-   :rtype: int
+.. doxygenfunction:: openmc_cell_get_temperature
 
 .. c:function:: int openmc_cell_get_density(int32_t index, const int32_t* instance, double* density)
 
@@ -580,6 +557,279 @@ Functions
    :return: Return status (negative if an error occurs)
    :rtype: int
 
+.. c:function:: int openmc_get_plot_index(int32_t id, int32_t* index)
+
+   Get the index in the plots array for a plot with a given ID.
+
+   :param int32_t id: Plot ID
+   :param int32_t* index: Index in the plots array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_plot_get_id(int32_t index, int32_t* id)
+
+   Get the ID of a plot.
+
+   :param int32_t index: Index in the plots array
+   :param int32_t* id: Plot ID
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_plot_set_id(int32_t index, int32_t id)
+
+   Set the ID of a plot.
+
+   :param int32_t index: Index in the plots array
+   :param int32_t id: Plot ID
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: size_t openmc_plots_size()
+
+   Number of plots currently allocated.
+
+   :return: Number of plots in the plots array
+   :rtype: size_t
+
+.. c:function:: int openmc_solidraytrace_plot_create(int32_t* index)
+
+   Create a new solid raytrace plot.
+
+   :param int32_t* index: Index of the newly created plot
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_pixels(int32_t index, int32_t* width, int32_t* height)
+
+   Get output pixel dimensions for a solid raytrace plot.
+
+   :param int32_t index: Index in the plots array
+   :param int32_t* width: Image width in pixels
+   :param int32_t* height: Image height in pixels
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_pixels(int32_t index, int32_t width, int32_t height)
+
+   Set output pixel dimensions for a solid raytrace plot.
+
+   :param int32_t index: Index in the plots array
+   :param int32_t width: Image width in pixels
+   :param int32_t height: Image height in pixels
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_color_by(int32_t index, int32_t* color_by)
+
+   Get the domain type used for coloring (0=materials, 1=cells).
+
+   :param int32_t index: Index in the plots array
+   :param int32_t* color_by: Coloring mode
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_color_by(int32_t index, int32_t color_by)
+
+   Set the domain type used for coloring (0=materials, 1=cells).
+
+   :param int32_t index: Index in the plots array
+   :param int32_t color_by: Coloring mode
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_default_colors(int32_t index)
+
+   Set default random colors for the current ``color_by`` mode.
+
+   :param int32_t index: Index in the plots array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_all_opaque(int32_t index)
+
+   Mark all domains in the current ``color_by`` mode as opaque.
+
+   :param int32_t index: Index in the plots array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_opaque(int32_t index, int32_t id, bool visible)
+
+   Set whether a specific domain ID is opaque (visible) in the rendered image.
+
+   :param int32_t index: Index in the plots array
+   :param int32_t id: Cell/material ID (based on ``color_by``)
+   :param bool visible: Whether the domain is opaque
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_color(int32_t index, int32_t id, uint8_t r, uint8_t g, uint8_t b)
+
+   Set RGB color for a specific domain ID.
+
+   :param int32_t index: Index in the plots array
+   :param int32_t id: Cell/material ID (based on ``color_by``)
+   :param uint8_t r: Red channel
+   :param uint8_t g: Green channel
+   :param uint8_t b: Blue channel
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_color(int32_t index, int32_t id, uint8_t* r, uint8_t* g, uint8_t* b)
+
+   Get RGB color for a specific domain ID.
+
+   :param int32_t index: Index in the plots array
+   :param int32_t id: Cell/material ID (based on ``color_by``)
+   :param uint8_t* r: Red channel
+   :param uint8_t* g: Green channel
+   :param uint8_t* b: Blue channel
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_camera_position(int32_t index, double* x, double* y, double* z)
+
+   Get camera position.
+
+   :param int32_t index: Index in the plots array
+   :param double* x: X coordinate
+   :param double* y: Y coordinate
+   :param double* z: Z coordinate
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_camera_position(int32_t index, double x, double y, double z)
+
+   Set camera position.
+
+   :param int32_t index: Index in the plots array
+   :param double x: X coordinate
+   :param double y: Y coordinate
+   :param double z: Z coordinate
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_look_at(int32_t index, double* x, double* y, double* z)
+
+   Get camera target point.
+
+   :param int32_t index: Index in the plots array
+   :param double* x: X coordinate
+   :param double* y: Y coordinate
+   :param double* z: Z coordinate
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_look_at(int32_t index, double x, double y, double z)
+
+   Set camera target point.
+
+   :param int32_t index: Index in the plots array
+   :param double x: X coordinate
+   :param double y: Y coordinate
+   :param double z: Z coordinate
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_up(int32_t index, double* x, double* y, double* z)
+
+   Get the camera up vector.
+
+   :param int32_t index: Index in the plots array
+   :param double* x: X component
+   :param double* y: Y component
+   :param double* z: Z component
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_up(int32_t index, double x, double y, double z)
+
+   Set the camera up vector.
+
+   :param int32_t index: Index in the plots array
+   :param double x: X component
+   :param double y: Y component
+   :param double z: Z component
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_light_position(int32_t index, double* x, double* y, double* z)
+
+   Get light source position.
+
+   :param int32_t index: Index in the plots array
+   :param double* x: X coordinate
+   :param double* y: Y coordinate
+   :param double* z: Z coordinate
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_light_position(int32_t index, double x, double y, double z)
+
+   Set light source position.
+
+   :param int32_t index: Index in the plots array
+   :param double x: X coordinate
+   :param double y: Y coordinate
+   :param double z: Z coordinate
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_fov(int32_t index, double* fov)
+
+   Get horizontal field of view in degrees.
+
+   :param int32_t index: Index in the plots array
+   :param double* fov: Field of view in degrees
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_fov(int32_t index, double fov)
+
+   Set horizontal field of view in degrees.
+
+   :param int32_t index: Index in the plots array
+   :param double fov: Field of view in degrees
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_get_diffuse_fraction(int32_t index, double* diffuse_fraction)
+
+   Get diffuse-light fraction.
+
+   :param int32_t index: Index in the plots array
+   :param double* diffuse_fraction: Diffuse fraction in [0, 1]
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_set_diffuse_fraction(int32_t index, double diffuse_fraction)
+
+   Set diffuse-light fraction.
+
+   :param int32_t index: Index in the plots array
+   :param double diffuse_fraction: Diffuse fraction in [0, 1]
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_update_view(int32_t index)
+
+   Recompute internal camera/view transforms after camera changes.
+
+   :param int32_t index: Index in the plots array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_solidraytrace_plot_create_image(int32_t index, uint8_t* data_out, int32_t width, int32_t height)
+
+   Render the plot to an RGB image buffer.
+
+   :param int32_t index: Index in the plots array
+   :param uint8_t* data_out: Output buffer of shape ``height*width*3``
+   :param int32_t width: Image width in pixels
+   :param int32_t height: Image height in pixels
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_reset()
 
    Resets all tally scores
@@ -601,6 +851,10 @@ Functions
    :return: Return status (negative if an error occurs)
    :rtype: int
 
+.. c:function:: void openmc_run_random_ray()
+
+   Run a random ray simulation
+
 .. c:function:: int openmc_set_n_batches(int32_t n_batches, bool set_max_batches, bool add_statepoint_batch)
 
    Set number of batches and number of max batches

docs/source/conf.py

@@ -11,7 +11,10 @@
 # All configuration values have a default; values that are commented out
 # serve to show the default.
 
-import sys, os
+import os
+from pathlib import Path
+import subprocess
+import sys
 
 # Determine if we're on Read the Docs server
 on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
@@ -37,6 +40,7 @@ sys.path.insert(0, os.path.abspath('../..'))
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
 extensions = [
+    'breathe',
     'sphinx.ext.autodoc',
     'sphinx.ext.napoleon',
     'sphinx.ext.autosummary',
@@ -47,6 +51,8 @@ extensions = [
 ]
 if not on_rtd:
     extensions.append('sphinxcontrib.rsvgconverter')
+    doxygen_dir = Path(__file__).parents[1] / 'doxygen'
+    subprocess.run(['doxygen'], cwd=doxygen_dir, check=True)
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
@@ -62,7 +68,7 @@ master_doc = 'index'
 
 # General information about the project.
 project = 'OpenMC'
-copyright = '2011-2025, Massachusetts Institute of Technology, UChicago Argonne LLC, and OpenMC contributors'
+copyright = '2011-2026, Massachusetts Institute of Technology, UChicago Argonne LLC, and OpenMC contributors'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -117,6 +123,11 @@ pygments_style = 'tango'
 # A list of ignored prefixes for module index sorting.
 #modindex_common_prefix = []
 
+# -- Options breathe + doxygen -------------------------------------------------
+
+breathe_projects = {"OpenMC": "../doxygen/xml"}
+breathe_default_project = "OpenMC"
+breathe_domain_by_file_pattern = {"*capi.h": "c"}
 
 # -- Options for HTML output ---------------------------------------------------
 

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/contributing.rst

@@ -111,7 +111,8 @@ The TC consists of the following individuals:
 - `Paul Romano <https://github.com/paulromano>`_
 - `Patrick Shriwise <https://github.com/pshriwise>`_
 - `Adam Nelson <https://github.com/nelsonag>`_
-- `Benoit Forget <https://github.com/bforget>`_
+- `Jonathan Shimwell <https://github.com/shimwell>`_
+- `John Tramm <https://github.com/jtramm>`_
 
 The Project Lead is Paul Romano.
 

docs/source/devguide/docbuild.rst

@@ -14,6 +14,11 @@ Python API. That is, from the root directory of the OpenMC repository:
 
     python -m pip install ".[docs]"
 
+The OpenMC documentation also uses Doxygen to automatically generate its
+C/C++ API documentation directly from the docstrings available in the source
+code. You will need to have a working installation of Doxygen to generate the
+documentation locally.
+
 -----------------------------------
 Building Documentation as a Webpage
 -----------------------------------

docs/source/devguide/index.rst

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

docs/source/devguide/policies.rst

@@ -21,8 +21,8 @@ C++ code in OpenMC must conform to the most recent C++ standard that is fully
 supported in the `version of the gcc compiler
 <https://gcc.gnu.org/projects/cxx-status.html>`_ that is distributed with the
 oldest version of Ubuntu that is still within its `standard support period
-<https://ubuntu.com/about/release-cycle>`_. Ubuntu 20.04 LTS will be supported
-through April 2025 and is distributed with gcc 9.3.0, which fully supports the
+<https://ubuntu.com/about/release-cycle>`_. Ubuntu 22.04 LTS will be supported
+through April 2027 and is distributed with gcc 11.4.0, which fully supports the
 C++17 standard.
 
 --------------------
@@ -31,5 +31,5 @@ CMake Version Policy
 
 Similar to the C++ standard policy, the minimum supported version of CMake
 corresponds to whatever version is distributed with the oldest version of Ubuntu
-still within its standard support period. Ubuntu 20.04 LTS is distributed with
-CMake 3.16.
+still within its standard support period. Ubuntu 22.04 LTS is distributed with
+CMake 3.22.

docs/source/devguide/styleguide.rst

@@ -30,7 +30,7 @@ whenever a file is saved. For example, `Visual Studio Code
 support for running clang-format.
 
 .. note::
-    OpenMC's CI uses `clang-format` version 15. A different version of `clang-format`
+    OpenMC's CI uses `clang-format` version 18. A different version of `clang-format`
     may produce different line changes and as a result fail the CI test.
 
 Miscellaneous

docs/source/devguide/tests.rst

@@ -37,6 +37,9 @@ Prerequisites
 - Some tests require `NJOY <https://www.njoy21.io/NJOY2016>`_ to preprocess
   cross section data. The test suite assumes that you have an ``njoy``
   executable available on your :envvar:`PATH`.
+- OpenMC should be compiled with ``-DOPENMC_ENABLE_STRICT_FP=on`` to ensure
+  reproducible floating-point results across platforms and optimization levels.
+  Without this flag, regression tests may not match reference values.
 
 Running Tests
 -------------
@@ -67,9 +70,11 @@ make sure you have satisfied all the prerequisites above. After you have done
 that, consider the following:
 
 - When building OpenMC, make sure you run CMake with
-  ``-DCMAKE_BUILD_TYPE=Debug``. Building with a release build will result in
-  some test failures due to differences in which compiler optimizations are
-  used.
+  ``-DOPENMC_ENABLE_STRICT_FP=on``. This prevents the compiler from applying
+  floating-point optimizations (such as replacing math library calls with
+  builtins or contracting multiply-add into FMA instructions) that can produce
+  bit-level differences across platforms and optimization levels. Any
+  ``CMAKE_BUILD_TYPE`` can be used.
 - Because tallies involve the sum of many floating point numbers, the
   non-associativity of floating point numbers can result in different answers
   especially when the number of threads is high (different order of operations).

docs/source/io_formats/collision_track.rst

@@ -10,7 +10,7 @@ may also be written after each batch when multiple files are requested
 (``collision_track.N.h5``) or when the run is performed in parallel. The file
 contains the information needed to reconstruct each recorded collision.
 
-The current revision of the collision track file format is 1.0.
+The current revision of the collision track file format is 1.2.
 
 **/**
 
@@ -33,13 +33,13 @@ The current revision of the collision track file format is 1.0.
               - ``event_mt`` (*int*) -- ENDF MT number identifying the reaction.
               - ``delayed_group`` (*int*) -- Delayed neutron group index (non-zero for delayed events).
               - ``cell_id`` (*int*) -- ID of the cell in which the collision occurred.
-              - ``nuclide_id`` (*int*) -- ZA identifier of the nuclide (ZZZAAAM format).
+              - ``nuclide_id`` (*int*) -- PDG number of the nuclide (100ZZZAAAM).
               - ``material_id`` (*int*) -- ID of the material containing the collision site.
               - ``universe_id`` (*int*) -- ID of the universe containing the collision site.
               - ``n_collision`` (*int*) -- Collision counter for the particle history.
-              - ``particle`` (*int*) -- Particle type (0=neutron, 1=photon, 2=electron, 3=positron).
-              - ``parent_id`` (*int64*) -- Unique ID of the parent particle.
-              - ``progeny_id`` (*int64*) -- Progeny ID of the particle.
+              - ``particle`` (*int32_t*) -- Particle type (PDG number).
+              - ``parent_id`` (*int64_t*) -- Unique ID of the parent particle.
+              - ``progeny_id`` (*int64_t*) -- Progeny ID of the particle.
 
 In an MPI run, OpenMC writes the combined dataset by gathering collision-track
 entries from all ranks before flushing them to disk, so the final file appears

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,11 +29,14 @@ 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>/**
 
 :Attributes: - **index** (*int*) -- Index used in results for this material
              - **volume** (*double*) -- Volume of this material in [cm^3]
+             - **name** (*char[]*) -- Name of this material
 
 **/nuclides/<name>/**
 

docs/source/io_formats/geometry.rst

@@ -38,11 +38,9 @@ Each ``<surface>`` element can have the following attributes or sub-elements:
 
   :boundary:
      The boundary condition for the surface. This can be "transmission",
-     "vacuum", "reflective", or "periodic". Periodic boundary conditions can
-     only be applied to x-, y-, and z-planes. Only axis-aligned periodicity is
-     supported, i.e., x-planes can only be paired with x-planes. Specify which
-     planes are periodic and the code will automatically identify which planes
-     are paired together.
+     "vacuum", "reflective", or "periodic". Specify which planes are
+     periodic and the code will automatically identify which planes are
+     paired together.
 
     *Default*: "transmission"
 
@@ -318,9 +316,10 @@ the following attributes or sub-elements:
     *Default*: None
 
   :orientation:
-    The orientation of the hexagonal lattice. The string "x" indicates that two
-    sides of the lattice are parallel to the x-axis, whereas the string "y"
-    indicates that two sides are parallel to the y-axis.
+    The orientation of the hexagonal lattice. The string "x" indicates that each
+    lattice element has two faces that are perpendicular to the x-axis, whereas
+    the string "y" indicates that each lattice element has two faces that are
+    perpendicular to the y-axis.
 
     *Default*: "y"
 
@@ -407,24 +406,55 @@ Each ``<dagmc_universe>`` element can have the following attributes or sub-eleme
 
     *Default*: None
 
-  :material_overrides:
-    This element contains information on material overrides to be applied to the
-    DAGMC universe. It has the following attributes and sub-elements:
+  :cell:
+    Zero or more ``<cell>`` sub-elements may appear to override properties of
+    individual DAGMC volumes. Each ``<cell>`` element supports the following
+    attributes and sub-elements:
 
-    :cell:
-      Material override information for a single cell. It contains the following
-      attributes and sub-elements:
+    :id:
+      The integer cell ID in the DAGMC geometry to override. Required.
 
-      :id:
-        The cell ID in the DAGMC geometry for which the material override will
-        apply.
+    :name:
+      An optional string label for the cell.
 
-      :materials:
-        A list of material IDs that will apply to instances of the cell. If the
-        list contains only one ID, it will replace the original material
-        assignment of all instances of the DAGMC cell. If the list contains more
-        than one material, each material ID of the list will be assigned to the
-        various instances of the DAGMC cell.
+      *Default*: None
+
+    :material:
+      The material ID to assign to this cell. Use ``void`` for vacuum. Multiple
+      space-separated IDs may be given to specify a distribmat (distributed
+      material) assignment. Required.
+
+    :temperature:
+      Temperature(s) in [K] to assign to the cell. Must be greater than or equal
+      to 0. Multiple space-separated values may be given.
+
+      *Default*: None
+
+    :density:
+      Density in [g/cm³] to assign to the cell. Must be greater than 0. Requires a non-void
+      material fill. Multiple space-separated values may be given.
+
+      *Default*: None
+
+    :volume:
+      Volume of the cell in [cm³].
+
+      .. note:: DAGMC can compute cell volumes exactly from the triangulated
+                mesh surfaces. Specifying a manual volume risks inconsistency
+                with that capability.
+
+      *Default*: None
+
+    The following standard ``<cell>`` attributes are **not** supported inside
+    ``<dagmc_universe>`` and will raise an error if present: ``region``,
+    ``fill``, ``universe``, ``translation``, ``rotation``.
+
+    .. deprecated::
+      The ``<material_overrides>`` sub-element (containing ``<cell_override>``
+      children with ``<material_ids>``) is deprecated. A deprecation warning is
+      emitted and the overrides are converted to the ``<cell>`` format at parse
+      time. It is an error to specify both ``<material_overrides>`` and
+      ``<cell>`` sub-elements on the same ``<dagmc_universe>``.
 
     *Default*: None
 

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/particle_restart.rst

@@ -4,7 +4,7 @@
 Particle Restart File Format
 ============================
 
-The current version of the particle restart file format is 2.0.
+The current version of the particle restart file format is 2.1.
 
 **/**
 
@@ -26,8 +26,7 @@ The current version of the particle restart file format is 2.0.
            - **run_mode** (*char[]*) -- Run mode used, either 'fixed source',
              'eigenvalue', or 'particle restart'.
            - **id** (*int8_t*) -- Unique identifier of the particle.
-           - **type** (*int*) -- Particle type (0=neutron, 1=photon, 2=electron,
-             3=positron)
+           - **type** (*int32_t*) -- Particle type (PDG number)
            - **weight** (*double*) -- Weight of the particle.
            - **energy** (*double*) -- Energy of the particle in eV for
              continuous-energy mode, or the energy group of the particle for

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
 ---------------------
@@ -85,6 +98,11 @@ sub-elements:
     A list of strings representing the nuclide, to define specific
     define specific target nuclide collisions to be banked.
 
+    .. note::
+      Electron and positron collision-track events are not associated with
+      a specific nuclide. If a ``nuclides`` entry is specified, these events
+      are omitted.
+
     *Default*: None
 
   :reactions:
@@ -277,6 +295,15 @@ ignored for all run modes other than "eigenvalue".
 
   *Default*: 1
 
+------------------------------
+``<ifp_n_generation>`` Element
+------------------------------
+
+The ``<ifp_n_generation>`` element indicates the number of generations to
+consider for the Iterated Fission Probability method.
+
+  *Default*: 10
+
 ----------------------
 ``<inactive>`` Element
 ----------------------
@@ -402,7 +429,25 @@ then, OpenMC will only use up to the :math:`P_1` data.
 ``<max_history_splits>`` Element
 --------------------------------
 
-The ``<max_history_splits>`` element indicates the number of times a particle can split during a history.
+The ``<max_history_splits>`` element indicates the number of times a particle
+can split during a history.
+
+  *Default*: 1000
+
+-----------------------------
+``<max_secondaries>`` Element
+-----------------------------
+
+The ``<max_secondaries>`` element indicates the maximum secondary bank size.
+
+  *Default*: 10000
+
+------------------------
+``<max_tracks>`` Element
+------------------------
+
+The ``<max_tracks>`` element indicates the maximum number of tracks written to a
+track file (per MPI process).
 
   *Default*: 1000
 
@@ -515,6 +560,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
 ---------------------
@@ -545,7 +602,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
@@ -553,6 +610,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".
@@ -660,14 +746,16 @@ pseudo-random number generator.
 
   *Default*: 1
 
---------------------
-``<stride>`` Element
---------------------
+-----------------------------------
+``<shared_secondary_bank>`` Element
+-----------------------------------
 
-The ``stride`` element is used to specify how many random numbers are allocated
-for each source particle history.
-
-  *Default*: 152,917
+  The ``shared_secondary_bank`` element indicates whether to use a shared
+  secondary particle bank. When enabled, secondary particles are collected into
+  a global bank, sorted for reproducibility, and load-balanced across MPI ranks
+  between generations. If not specified, the shared secondary bank is enabled
+  automatically for fixed-source simulations with weight windows active, and
+  disabled otherwise.
 
 .. _source_element:
 
@@ -689,12 +777,15 @@ attributes/sub-elements:
     *Default*: 1.0
 
   :type:
-    Indicator of source type. One of ``independent``, ``file``, ``compiled``, or
-    ``mesh``. The type of the source will be determined by this attribute if it
-    is present.
+    Indicator of source type. One of ``independent``, ``file``, ``compiled``,
+    ``mesh``, or ``tokamak``. The type of the source will be determined by this
+    attribute if it is present.
 
   :particle:
-    The source particle type, either ``neutron`` or ``photon``.
+    The source particle type, specified as a PDG number or a string alias (e.g.,
+    ``neutron``/``n``, ``photon``/``gamma``, ``electron``, ``positron``,
+    ``proton``/``p``, ``deuteron``/``d``, ``triton``/``t``, ``alpha``, or GNDS
+    nuclide names like ``Fe57``).
 
     *Default*: neutron
 
@@ -784,6 +875,7 @@ attributes/sub-elements:
 
       For a "cylindrical" distribution, no parameters are specified. Instead,
       the ``r``, ``phi``, ``z``, and ``origin`` elements must be specified.
+      Optionally, the ``r_dir`` and ``z_dir`` elements could be specified.
 
       For a "spherical" distribution, no parameters are specified. Instead,
       the ``r``, ``theta``, ``phi``, and ``origin`` elements must be specified.
@@ -815,6 +907,10 @@ attributes/sub-elements:
       of a univariate probability distribution (see the description in
       :ref:`univariate`).
 
+    :r_dir:
+      For "cylindrical" distributions, this element specifies the direction
+      of the cylinder r-axis at phi=0. Defaults to (1.0, 0.0, 0.0).
+
     :theta:
       For a "spherical" distribution, this element specifies the distribution
       of theta-coordinates. The necessary sub-elements/attributes are those of a
@@ -827,6 +923,10 @@ attributes/sub-elements:
       sub-elements/attributes are those of a univariate probability
       distribution (see the description in :ref:`univariate`).
 
+    :z_dir:
+      For "cylindrical" distributions, this element specifies the direction
+      of the cylinder z-axis. Defaults to (0.0, 0.0, 1.0).
+
     :origin:
       For "cylindrical and "spherical" distributions, this element specifies
       the coordinates for the origin of the coordinate system.
@@ -844,13 +944,18 @@ attributes/sub-elements:
       relative source strength of each mesh element or each point in the cloud.
 
     :volume_normalized:
-      For "mesh" spatial distrubtions, this optional boolean element specifies
+      For "mesh" spatial distributions, this optional boolean element specifies
       whether the vector of relative strengths should be multiplied by the mesh
       element volume. This is most common if the strengths represent a source
       per unit volume.
 
       *Default*: false
 
+    :bias:
+      For "mesh" and "cloud" spatial distributions, this optional element
+      specifies floating point values corresponding to alternative probabilities
+      for each value/component to use for biased sampling.
+
   :angle:
     An element specifying the angular distribution of source sites. This element
     has the following attributes:
@@ -883,6 +988,10 @@ attributes/sub-elements:
       are those of a univariate probability distribution (see the description in
       :ref:`univariate`).
 
+    :bias:
+      For "isotropic" angular distributions, this optional element specifies a
+      "mu-phi" angular distribution used for biased sampling.
+
   :energy:
     An element specifying the energy distribution of source sites. The necessary
     sub-elements/attributes are those of a univariate probability distribution
@@ -906,6 +1015,84 @@ attributes/sub-elements:
     mesh element and follows the format for :ref:`source_element`. The number of
     ``<source>`` sub-elements should correspond to the number of mesh elements.
 
+  For a source with ``type="tokamak"``, the spatial distribution is described by
+  a Miller-style flux-surface parameterization and the following sub-elements
+  are used instead of the ``space`` element:
+
+  :major_radius:
+    The major radius :math:`R_0` of the plasma in [cm].
+
+  :minor_radius:
+    The minor radius :math:`a` of the plasma in [cm]. Must be smaller than
+    ``major_radius``.
+
+  :elongation:
+    The plasma elongation :math:`\kappa` (must be > 0).
+
+  :triangularity:
+    The plasma triangularity :math:`\delta` (must be in [-1, 1]). Negative
+    values describe negative-triangularity plasmas.
+
+  :shafranov_shift:
+    The Shafranov shift :math:`\Delta` in [cm] (must be >= 0 and less than
+    ``minor_radius``/2).
+
+  :r_over_a:
+    A list of normalized minor-radius grid points :math:`r/a`. Must be strictly
+    increasing, start at 0, and end at 1.
+
+  :emission_density:
+    A list of neutron emission densities :math:`S(r)` evaluated at each
+    ``r_over_a`` grid point (arbitrary units, must be non-negative). Only the
+    shape matters, since the profile is normalized internally. Values are
+    interpolated linearly between grid points and the profile is refined on an
+    internal grid for radial sampling. Must have the same length as
+    ``r_over_a`` and contain at least one positive value.
+
+  :phi_start:
+    The starting toroidal angle in [rad].
+
+    *Default*: 0.0
+
+  :phi_extent:
+    The toroidal angle extent in [rad]. The source is sampled uniformly in
+    :math:`[\phi_\text{start},\ \phi_\text{start} + \phi_\text{extent}]`.
+
+    *Default*: :math:`2\pi`
+
+  :n_alpha:
+    The number of poloidal-angle grid points used to build the sampling CDFs
+    (must be > 2). Larger values reduce discretization bias; values below 51
+    produce a warning.
+
+    *Default*: 101
+
+  :vertical_shift:
+    A vertical shift of the plasma center in [cm].
+
+    *Default*: 0.0
+
+  :energy:
+    For a tokamak source, one or more ``energy`` sub-elements specify the
+    neutron energy distribution(s). Either a single distribution is given (used
+    at all radii) or exactly one distribution per ``r_over_a`` grid point is
+    given, in which case the energy is sampled from one of the two
+    distributions bracketing the sampled radius, selected stochastically with
+    probability proportional to the proximity of the radius to each grid point
+    (stochastic interpolation). Each follows the format of a univariate
+    probability distribution (see :ref:`univariate`).
+
+  :time:
+    An optional ``time`` sub-element specifying the time distribution of source
+    particles, following the format of a univariate probability distribution
+    (see :ref:`univariate`).
+
+    *Default*: particles are born at :math:`t=0`
+
+  .. note:: Biased sampling can be applied to the spatial and energy distributions
+            of a source by using the ``<bias>`` sub-element (see
+            :ref:`univariate` for details on how to specify bias distributions).
+
   :constraints:
     This sub-element indicates the presence of constraints on sampled source
     sites (see :ref:`usersguide_source_constraints` for details). It may have
@@ -952,17 +1139,19 @@ variable and whose sub-elements/attributes are as follows:
 
 :type:
   The type of the distribution. Valid options are "uniform", "discrete",
-  "tabular", "maxwell", "watt", and "mixture". The "uniform" option produces
-  variates sampled from a uniform distribution over a finite interval. The
-  "discrete" option produces random variates that can assume a finite number
-  of values (i.e., a distribution characterized by a probability mass function).
-  The "tabular" option produces random variates sampled from a tabulated
-  distribution where the density function is either a histogram or
+  "tabular", "maxwell", "watt", "mixture", and "decay_spectrum". The "uniform"
+  option produces variates sampled from a uniform distribution over a finite
+  interval. The "discrete" option produces random variates that can assume a
+  finite number of values (i.e., a distribution characterized by a probability
+  mass function). The "tabular" option produces random variates sampled from a
+  tabulated distribution where the density function is either a histogram or
   linearly-interpolated between tabulated points. The "watt" option produces
   random variates is sampled from a Watt fission spectrum (only used for
   energies). The "maxwell" option produce variates sampled from a Maxwell
-  fission spectrum (only used for energies). The "mixture" option produces samples
-  from univariate sub-distributions with given probabilities.
+  fission spectrum (only used for energies). The "mixture" option produces
+  samples from univariate sub-distributions with given probabilities. The
+  "decay_spectrum" option produces photon energies sampled from decay photon
+  spectra in a depletion chain (only used for energies).
 
   *Default*: None
 
@@ -980,6 +1169,10 @@ variable and whose sub-elements/attributes are as follows:
   :math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x`
   points are given first followed by corresponding :math:`p` points.
 
+  For a "decay_spectrum" distribution, ``parameters`` gives the atom densities
+  in [atom/b-cm] for the nuclides listed in the ``nuclides`` element, in the
+  same order.
+
   For a "watt" distribution, ``parameters`` should be given as two real numbers
   :math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c
   e^{-x/a} \sinh \sqrt{b \, x} dx`.
@@ -998,13 +1191,41 @@ variable and whose sub-elements/attributes are as follows:
   *Default*: histogram
 
 :pair:
-  For a "mixture" distribution, this element provides a distribution and its corresponding probability.
+  For a "mixture" distribution, this element provides a distribution and its
+  corresponding probability.
 
   :probability:
-    An attribute or ``pair`` that provides the probability of a univariate distribution within a "mixture" distribution.
+    An attribute or ``pair`` that provides the probability of a univariate
+    distribution within a "mixture" distribution.
 
   :dist:
-    This sub-element of a ``pair`` element provides information on the corresponding univariate distribution.
+    This sub-element of a ``pair`` element provides information on the
+    corresponding univariate distribution.
+
+:volume:
+  For a "decay_spectrum" distribution, this attribute specifies the source
+  region volume in cm\ :sup:`3`. It is used together with atom densities to
+  determine the absolute photon emission rate. When a source uses a
+  "decay_spectrum" energy distribution, the source strength is set from this
+  emission rate.
+
+:nuclides:
+  For a "decay_spectrum" distribution, this element specifies a
+  whitespace-separated list of nuclide names contributing to the decay photon
+  source. The atom densities for these nuclides are given by the ``parameters``
+  element in the same order. Nuclides are resolved against the depletion chain,
+  and nuclides without decay photon spectra do not contribute to the
+  distribution.
+
+:bias:
+  This optional element specifies a biased distribution for importance sampling.
+  For continuous distributions, the ``bias`` element should contain another
+  univariate distribution with the same support (interval) as the parent
+  distribution. For discrete distributions, the ``bias`` element should contain
+  floating point values corresponding to alternative probabilities for each
+  value/component to be used for biased sampling.
+
+  *Default*: None
 
 ---------------------------------------
 ``<source_rejection_fraction>`` Element
@@ -1016,23 +1237,6 @@ based on constraints.
 
    *Default*: 0.05
 
--------------------------
-``<state_point>`` Element
--------------------------
-
-The ``<state_point>`` element indicates at what batches a state point file
-should be written. A state point file can be used to restart a run or to get
-tally results at any batch. The default behavior when using this tag is to
-write out the source bank in the state_point file. This behavior can be
-customized by using the ``<source_point>`` element. This element has the
-following attributes/sub-elements:
-
-  :batches:
-    A list of integers separated by spaces indicating at what batches a state
-    point file should be written.
-
-    *Default*: Last batch only
-
 --------------------------
 ``<source_point>`` Element
 --------------------------
@@ -1082,6 +1286,32 @@ attributes/sub-elements:
 
     *Default*: false
 
+-------------------------
+``<state_point>`` Element
+-------------------------
+
+The ``<state_point>`` element indicates at what batches a state point file
+should be written. A state point file can be used to restart a run or to get
+tally results at any batch. The default behavior when using this tag is to
+write out the source bank in the state_point file. This behavior can be
+customized by using the ``<source_point>`` element. This element has the
+following attributes/sub-elements:
+
+  :batches:
+    A list of integers separated by spaces indicating at what batches a state
+    point file should be written.
+
+    *Default*: Last batch only
+
+--------------------
+``<stride>`` Element
+--------------------
+
+The ``stride`` element is used to specify how many random numbers are allocated
+for each source particle history.
+
+  *Default*: 152,917
+
 ------------------------------
 ``<surf_source_read>`` Element
 ------------------------------
@@ -1169,6 +1399,23 @@ attributes/sub-elements:
           are not eligible to store any particles when using ``cell``, ``cellfrom``
           or ``cellto`` attributes. It is recommended to use surface IDs instead.
 
+------------------------------------
+``<surface_grazing_cutoff>`` Element
+------------------------------------
+
+The ``<surface_grazing_cutoff>`` element specifies the surface flux cosine cutoff.
+
+  *Default*: 0.001
+
+-----------------------------------
+``<surface_grazing_ratio>`` Element
+-----------------------------------
+
+The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
+substitution ratio.
+
+  *Default*: 0.5
+
 ------------------------------
 ``<survival_biasing>`` Element
 ------------------------------
@@ -1352,6 +1599,15 @@ has the following attributes/sub-elements:
               for fixed source and small criticality calculations, but is very
               optimistic for highly coupled full-core reactor problems.
 
+-------------------------------------
+``<uniform_source_sampling>`` Element
+-------------------------------------
+
+The ``<uniform_source_sampling>`` element indicates whether to sample among
+multiple sources uniformly, applying their strengths as weights to sampled
+particles.
+
+  *Default*: False
 
 ------------------------
 ``<ufs_mesh>`` Element
@@ -1364,6 +1620,16 @@ Agency Monte Carlo Performance Benchmark Problem," Proceedings of *Physor 2012*,
 Knoxville, TN (2012). The mesh should cover all possible fissionable materials
 in the problem and is specified using a :ref:`mesh_element`.
 
+-------------------------------
+``<use_decay_photons>`` Element
+-------------------------------
+
+The ``<use_decay_photons>`` element indicates whether to produce decay photons
+from neutron reactions instead of prompt photons. This is used in conjunction
+with the direct 1-step method for shutdown dose rate calculations.
+
+  *Default*: False
+
 .. _verbosity:
 
 -----------------------
@@ -1465,7 +1731,8 @@ sub-elements/attributes:
     *Default*: None
 
   :particle_type:
-    The particle that the weight windows will apply to (e.g., 'neutron')
+    The particle that the weight windows will apply to, specified as a PDG
+    code or string (e.g., ``neutron``).
 
     *Default*: 'neutron'
 
@@ -1525,7 +1792,8 @@ mesh-based weight windows.
     *Default*: None
 
   :particle_type:
-    The particle that the weight windows will apply to (e.g., 'neutron')
+    The particle that the weight windows will apply to, specified as a PDG
+    code or string (e.g., ``neutron``).
 
     *Default*: neutron
 
@@ -1569,6 +1837,14 @@ mesh-based weight windows.
 
         *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
 ---------------------------------------
@@ -1594,3 +1870,21 @@ following sub-elements/attributes:
 
   The ``weight_windows_file`` element has no attributes and contains the path to
   a weight windows HDF5 file to load during simulation initialization.
+
+-------------------------------
+``<weight_windows_on>`` Element
+-------------------------------
+
+  The ``weight_windows_on`` element indicates whether weight windows are
+  enabled.
+
+  *Default*: False
+
+----------------------------------
+``<write_initial_source>`` Element
+----------------------------------
+
+  The ``write_initial_source`` element indicates whether to write the initial
+  source distribution to file.
+
+  *Default*: False

docs/source/io_formats/source.rst

@@ -15,6 +15,8 @@ following the same format.
 **/**
 
 :Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
+             - **version** (*int[2]*) -- Major and minor version of the source
+               file format.
 
 :Datasets:
 
@@ -22,5 +24,5 @@ following the same format.
              particle. The compound type has fields ``r``, ``u``, ``E``,
              ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``,
              which represent the position, direction, energy, time, weight,
-             delayed group, surface ID, and particle type (0=neutron, 1=photon,
-             2=electron, 3=positron), respectively.
+             delayed group, surface ID, and particle type (PDG number),
+             respectively.

docs/source/io_formats/statepoint.rst

@@ -4,7 +4,7 @@
 State Point File Format
 =======================
 
-The current version of the statepoint file format is 18.1.
+The current version of the statepoint file format is 18.2.
 
 **/**
 
@@ -56,8 +56,8 @@ The current version of the statepoint file format is 18.1.
              ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and
              ``particle``, which represent the position, direction, energy,
              time, weight, delayed group, surface ID, and particle type
-             (0=neutron, 1=photon, 2=electron, 3=positron), respectively. Only
-             present when `run_mode` is 'eigenvalue'.
+             (PDG number), respectively. Only present when `run_mode` is
+             'eigenvalue'.
 
 **/tallies/**
 

docs/source/io_formats/tallies.rst

@@ -142,9 +142,9 @@ attributes/sub-elements:
 
   :type:
     The type of the filter. Accepted options are "cell", "cellfrom",
-    "cellborn", "surface", "material", "universe", "energy", "energyout", "mu",
-    "polar", "azimuthal", "mesh", "distribcell", "delayedgroup",
-    "energyfunction", and "particle".
+    "cellborn", "surface", "material", "universe", "energy", "energyout",
+    "mu", "polar", "azimuthal", "mesh", "distribcell", "delayedgroup",
+    "energyfunction", "particle", and "particleproduction".
 
   :bins:
      A description of the bins for each type of filter can be found in
@@ -318,8 +318,34 @@ should be set to:
   they use ``energy`` and ``y``.
 
 :particle:
-  A list of integers indicating the type of particles to tally ('neutron' = 1,
-  'photon' = 2, 'electron' = 3, 'positron' = 4).
+  A list of particle identifiers to tally, specified as strings (e.g.,
+  ``neutron``, ``photon``, ``He4``) or as integer PDG numbers.
+
+:particleproduction:
+  This filter tallies secondary particles produced in reactions, binned by
+  particle type and, optionally, by energy. Unlike other energy filters, the
+  weight applied is the weight of the secondary particle. To obtain secondary
+  particle production rates, use this filter with the ``events`` score.
+
+  The filter uses the following sub-elements instead of ``bins``:
+
+  :particles:
+    A space-separated list of secondary particle types to tally (e.g.,
+    ``photon``, ``neutron``, ``electron``).
+
+  :energies:
+    An optional monotonically increasing list of energy boundaries in [eV]
+    for binning the secondary particle energies. If omitted, total production
+    is tallied without energy binning.
+
+  For example, to tally photon and neutron production in three energy groups:
+
+  .. code-block:: xml
+
+      <filter id="1" type="particleproduction">
+        <particles>photon neutron</particles>
+        <energies>0.0 1.0e5 1.0e6 20.0e6</energies>
+      </filter>
 
 ------------------
 ``<mesh>`` Element

docs/source/io_formats/track.rst

@@ -4,7 +4,7 @@
 Track File Format
 =================
 
-The current revision of the particle track file format is 3.0.
+The current revision of the particle track file format is 3.1.
 
 **/**
 
@@ -32,6 +32,5 @@ The current revision of the particle track file format is 3.0.
                             the array for each primary/secondary particle. The
                             last offset should match the total size of the
                             array.
-                          - **particles** (*int[]*) -- Particle type for each
-                            primary/secondary particle (0=neutron, 1=photon,
-                            2=electron, 3=positron).
+                          - **particles** (*int32_t[]*) -- Particle type for
+                            each primary/secondary particle (PDG number).

docs/source/license.rst

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

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/tallies.rst

@@ -205,7 +205,71 @@ had a collision at every event. Thus, for tallies with outgoing-energy filters
 or for tallies of scattering moments (which require the scattering cosine of
 the change-in-angle), we must use an analog estimator.
 
-.. TODO: Add description of surface current tallies
+-----------------------------------
+Surface-Integrated Flux and Current
+-----------------------------------
+
+Surface tallies allow you to measure particle behavior as they cross specific
+boundaries in your geometry. Unlike volume tallies, which integrate over a
+volumetric region, surface tallies capture the current or flux passing through a
+surface. Surface tallies are estimated using an analog estimator.
+
+Current Score
+-------------
+
+When tallying the current across a surface, we simply count the weight of
+particles that cross the surface of interest:
+
+
+.. math::
+    :label: analog-current-estimator
+
+    J = \frac{1}{W} \sum_{i \in S} w_i.
+
+where :math:`J` is the area-integrated current passing through surface
+:math:`S`, :math:`W` is the total starting weight of the particles, and
+:math:`w_i` is the weight of the particle as it crosses the surface :math:`S`.
+
+Flux Score
+----------
+
+When tallying flux over a surface, we use the relationship between current and
+flux:
+
+
+.. math::
+    :label: surface-flux-estimator
+
+    \phi_S = \frac{1}{W} \sum_{i \in S} \frac{w_i}{|\mu|}.
+
+where :math:`\phi_S` is the area-integrated flux over surface :math:`S`,
+:math:`W` is the total starting weight of the particles, :math:`w_i` is the
+weight of the particle as it crosses the surface :math:`S` and :math:`\mu` is
+the cosine of angle between the particle direction and the surface normal.
+
+This equation diverges when the particle crossing the surface is nearly parallel
+to it (that is, as :math:`\mu` approaches zero). To remove this divergence,
+OpenMC scores:
+
+.. math::
+    :label: modified-surface-flux-estimator
+
+    \phi_S = \frac{1}{W} \sum_{i \in S} w_i f(\mu).
+
+and the function :math:`f` is defined by:
+
+.. math::
+     f(\mu) = \begin{cases}
+      \frac{1}{|\mu|} & |\mu| > \mu_\text{cut} \\
+      \frac{1}{c\mu_\text{cut}} & |\mu| \le \mu_\text{cut}
+     \end{cases}
+
+where :math:`\mu_\text{cut}` is the grazing cosine cutoff and :math:`c` is the
+cosine substitution ratio. The parameters :math:`\mu_\text{cut}` and :math:`c`
+can be set by the user via the :attr:`openmc.Settings.surface_grazing_cutoff`
+and :attr:`openmc.Settings.surface_grazing_ratio` attributes, respectively. The
+default values for these parameters are 0.001 and 0.5 as recommended by
+`Favorite, Thomas, and Booth <https://doi.org/10.13182/NSE09-72>`_.
 
 .. _tallies_statistics:
 

docs/source/methods/variance_reduction.rst

@@ -22,12 +22,14 @@ not experience a single scoring event, even after billions of analog histories.
 Variance reduction techniques aim to either flatten the global uncertainty
 distribution, such that all regions of phase space have a fairly similar
 uncertainty, or to reduce the uncertainty in specific locations (such as a
-detector). There are two strategies available in OpenMC for variance reduction:
-the Monte Carlo MAGIC method and the FW-CADIS method. Both strategies work by
-developing a weight window mesh that can be utilized by subsequent Monte Carlo
-solves to split particles heading towards areas of lower flux densities while
-terminating particles in higher flux regions---all while maintaining a fair
-game.
+detector). There are three strategies available in OpenMC for variance
+reduction: weight windows generated via the MAGIC method or the FW-CADIS method,
+and source biasing. Both weight windowing strategies work by developing a mesh
+that can be utilized by subsequent Monte Carlo solves to split particles heading
+towards areas of lower flux densities while terminating particles in higher flux
+regions. In contrast, source biasing modifies source site sampling behavior to
+preferentially track particles more likely to reach phase space regions of
+interest.
 
 ------------
 MAGIC Method
@@ -80,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::
@@ -132,3 +134,83 @@ aware of this.
     :label: variance_fom
 
     \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:
+
+--------------
+Source Biasing
+--------------
+
+In contrast to the previous two methods that introduce population controls
+during transport, source biasing modifies the sampling of the external source
+distribution. The basic premise of the technique is that for each spatial,
+angular, energy, or time distribution of a source, an additional distribution
+can be specified provided that the two share a common support (set of points
+where the distribution is nonzero). Samples are then drawn from this "bias"
+distribution, which can be chosen to preferentially direct particles towards
+phase space regions of interest. In order to avoid biasing the tally results,
+however, a weight adjustment is applied to each sampled site as described below.
+
+Assume that the unbiased probability density function of a random variable
+:math:`X:x \rightarrow \mathbb{R}` is given by :math:`f(x)`, but that using the
+biased distribution :math:`g(x)` will result in a greater number of particle
+trajectories reaching some phase space region of interest. Then a sample
+:math:`x_0` may be drawn from :math:`g(x)` while maintaining a fair game,
+provided that its weight is adjusted as:
+
+.. math::
+    :label: source_bias
+
+    w = w_0 \times \frac{f(x_0)}{g(x_0)}
+
+where :math:`w_0` is the weight of an unbiased sample from :math:`f(x)`,
+typically unity.
+
+Returning now to Equation :eq:`source_bias`, the requirement for common support
+becomes evident. If :math:`\mathrm{supp} (g)` fully contains but is not
+identical to :math:`\mathrm{supp} (f)`, then some samples from :math:`g(x)` will
+correspond to points where :math:`f(x) = 0`. Thus these source sites would be
+assigned a starting weight of 0, meaning the particles would be killed
+immediately upon transport, effectively wasting computation time. Conversely, if
+:math:`\mathrm{supp} (g)` is fully contained by but not identical to
+:math:`\mathrm{supp} (f)`, the contributions of some regions outside
+:math:`\mathrm{supp} (g)` will not be counted towards the integral, potentially
+biasing the tally. The weight assigned to such points would be undefined since
+:math:`g(x) = \mathbf{0}` at these points.
+
+When an independent source is sampled in OpenMC, the particle's coordinate in
+each variable of phase space :math:`(\mathbf{r},\mathbf{\Omega},E,t)` is
+successively drawn from an independent probability distribution. Multiple
+variables can be biased, in which case the resultant weight :math:`w` applied to
+the particle is the product of the weights assigned from all sampled
+distributions: space, angle, energy, and time, as shown in Equation
+:eq:`tot_wgt`.
+
+.. math::
+    :label: tot_wgt
+
+    w = w_r \times w_{\Omega} \times w_E \times w_t
+
+Finally, source biasing and weight windows serve different purposes. Source
+biasing changes how particles are born, allowing the initial source sites to be
+sampled preferentially from important regions of phase space (space, angle,
+energy, and time) with an accompanying weight adjustment. Weight windows, by
+contrast, apply population control during transport (splitting and Russian
+roulette) to help particles reach and contribute in important regions as they
+move through the system. Because particle transport proceeds as usual after a
+biased source is sampled, particle attenuation in optically thick regions
+outside the source volume will not be affected by source biasing; in such
+scenarios, transport biasing techniques such as weight windows are often more
+effective.

docs/source/pythonapi/base.rst

@@ -26,6 +26,7 @@ Simulation Settings
    openmc.FileSource
    openmc.CompiledSource
    openmc.MeshSource
+   openmc.TokamakSource
    openmc.SourceParticle
    openmc.VolumeCalculation
    openmc.Settings
@@ -132,6 +133,7 @@ Constructing Tallies
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter
+   openmc.ParticleProductionFilter
    openmc.MuFilter
    openmc.MuSurfaceFilter
    openmc.PolarFilter
@@ -148,6 +150,7 @@ Constructing Tallies
    openmc.ZernikeRadialFilter
    openmc.ParentNuclideFilter
    openmc.ParticleFilter
+   openmc.ReactionFilter
    openmc.MeshMaterialVolumes
    openmc.Trigger
    openmc.TallyDerivative
@@ -176,7 +179,8 @@ Geometry Plotting
    :nosignatures:
    :template: myclass.rst
 
-   openmc.Plot
+   openmc.SlicePlot
+   openmc.VoxelPlot
    openmc.WireframeRayTracePlot
    openmc.SolidRayTracePlot
    openmc.Plots

docs/source/pythonapi/capi.rst

@@ -40,6 +40,7 @@ Functions
    reset_timers
    run
    run_in_memory
+   run_random_ray
    sample_external_source
    simulation_finalize
    simulation_init
@@ -81,12 +82,15 @@ Classes
    Nuclide
    ParentNuclideFilter
    ParticleFilter
+   ParticleProductionFilter
    PolarFilter
+   ReactionFilter
    RectilinearMesh
    RegularMesh
    SpatialLegendreFilter
    SphericalHarmonicsFilter
    SphericalMesh
+   SolidRayTracePlot
    SurfaceFilter
    Tally
    TemporarySession
@@ -124,6 +128,12 @@ Data
 
    :type: dict
 
+.. data:: plots
+
+   Mapping of plot ID to :class:`openmc.lib.SolidRayTracePlot` instances.
+
+   :type: dict
+
 .. data:: nuclides
 
    Mapping of nuclide name to :class:`openmc.lib.Nuclide` instances.

docs/source/pythonapi/data.rst

@@ -71,6 +71,8 @@ Core Functions
     isotopes
     kalbach_slope
     linearize
+    mass_attenuation_coefficient
+    mass_energy_absorption_coefficient
     thin
     water_density
     zam

docs/source/pythonapi/mgxs.rst

@@ -11,6 +11,16 @@ Module Variables
 .. autodata:: openmc.mgxs.GROUP_STRUCTURES
     :annotation:
 
+Functions 
++++++++++
+
+.. autosummary::
+    :toctree: generated
+    :nosignatures:
+    :template: myfunction.rst
+
+    openmc.mgxs.convert_flux_groups
+
 Classes
 +++++++
 

docs/source/pythonapi/stats.rst

@@ -22,6 +22,7 @@ Univariate Probability Distributions
    openmc.stats.Legendre
    openmc.stats.Mixture
    openmc.stats.Normal
+   openmc.stats.DecaySpectrum
 
 .. autosummary::
    :toctree: generated
@@ -29,6 +30,7 @@ Univariate Probability Distributions
    :template: myfunction.rst
 
    openmc.stats.delta_function
+   openmc.stats.fusion_neutron_spectrum
    openmc.stats.muir
 
 Angular Distributions
@@ -67,3 +69,4 @@ Spatial Distributions
    :template: myfunction.rst
 
    openmc.stats.spherical_uniform
+   openmc.stats.cylindrical_uniform

docs/source/quickinstall.rst

@@ -119,7 +119,7 @@ packages should be installed, for example in Homebrew via:
 
 .. code-block:: sh
 
-   brew install llvm cmake xtensor hdf5 python libomp libpng
+   brew install llvm cmake hdf5 python libomp libpng
 
 The compiler provided by the above LLVM package should be used in place of the
 one provisioned by XCode, which does not support the multithreading library used

docs/source/usersguide/decay_sources.rst

@@ -132,8 +132,9 @@ can be run::
     r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes)
 
 If not specified otherwise, a photon transport calculation is run at each time
-in the depletion schedule. That means in the case above, we would see three
-photon transport calculations. To specify specific times at which photon
+in the depletion schedule for which a decay photon source exists. Times without
+a decay photon source, such as the initial state of a model containing only
+stable nuclides, are omitted. To specify particular times at which photon
 transport calculations should be run, pass the ``photon_time_indices`` argument.
 For example, if we wanted to run a photon transport calculation only on the last
 time (after the 5 hour decay), we would run::
@@ -141,6 +142,19 @@ time (after the 5 hour decay), we would run::
     r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
             photon_time_indices=[2])
 
+To attribute photon tally results to their parent radionuclides, set
+``by_parent_nuclide=True``. This automatically adds a
+:class:`openmc.ParentNuclideFilter` to every photon tally that does not already
+have one. The filter bins are the union of radionuclides contributing to the
+prepared decay photon sources. The resulting bins can be used directly when
+inspecting the tally results::
+
+    r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes,
+            photon_time_indices=[2], by_parent_nuclide=True)
+
+    photon_tally = r2s.results['photon_tallies'][2][0]
+    tally_by_parent = photon_tally.get_pandas_dataframe()
+
 After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager`
 instance will have a ``results`` dictionary that allows you to directly access
 results from each of the steps. It will also write out all the output files into
@@ -148,12 +162,13 @@ a directory that is named "r2s_<timestamp>/". The ``output_dir`` argument to the
 :meth:`~openmc.deplete.R2SManager.run` method enables you to override the
 default output directory name if desired.
 
-The :meth:`~openmc.deplete.R2SManager.run` method actually runs three
+The :meth:`~openmc.deplete.R2SManager.run` method actually runs four
 lower-level methods under the hood::
 
     r2s.step1_neutron_transport(...)
     r2s.step2_activation(...)
-    r2s.step3_photon_transport(...)
+    r2s.step3_photon_source(...)
+    r2s.step4_photon_transport(...)
 
 For users looking for more control over the calculation, these lower-level
 methods can be used in lieu of the :meth:`openmc.deplete.R2SManager.run` method.
@@ -190,6 +205,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
 ================================
 
@@ -236,4 +270,3 @@ relevant tallies. This can be done with the aid of the
 
     # Apply time correction factors
     tally = d1s.apply_time_correction(dose_tally, factors, time_index)
-

docs/source/usersguide/depletion.rst

@@ -449,3 +449,42 @@ to transfer xenon from one material to another, you'd use::
     ...
 
     integrator.add_transfer_rate(mat1, ['Xe'], 0.1, destination_material=mat2)
+
+Comparing to Other Codes
+========================
+
+Comparing depletion results from OpenMC with those from another code, such as
+MCNP or Serpent, requires more than constructing equivalent transport models.
+At each depletion step, differences in the transport solution, nuclear data,
+reaction rate normalization, and numerical integration can all affect the
+result. Small differences can also accumulate over successive depletion steps.
+
+For a meaningful comparison, align as many of the following inputs and methods
+as possible:
+
+- Geometry and material definitions and associated physical properties such as
+  temperature
+- Neutron cross section library (e.g., ENDF/B-VIII.0)
+- Treatment of thermal scattering and unresolved resonance probability tables
+- Neutron reactions accounted for in the depletion chain
+- Decay data in the depletion chain
+- Isomeric branching ratios for reactions in the depletion chain
+- Fission product yields in the depletion chain
+- Fission product yield interpolation method
+  (``CoupledOperator(fission_yield_mode=...)``)
+- Reaction rate normalization, including fission Q values
+  (``CoupledOperator(fission_q=...)``)
+- Depletion integration method (``PredictorIntegrator``, ``CECMIntegrator``,
+  etc.) and time-step sizes
+
+When comparing to codes that use ACE format cross sections, it is recommended to
+directly convert the ACE files to HDF5 format using functionality from the
+:mod:`openmc.data` module (see :ref:`create_xs_library`). Some of the
+LANL-distributed ACE libraries used with MCNP have also been converted to HDF5
+format and are available for download at https://openmc.org/data.
+
+Even after these choices have been aligned, exact agreement should not be
+expected. Codes may use different approximations or numerical methods that
+cannot be configured identically. When investigating a discrepancy, first
+compare transport results and one-group reaction rates at the initial time, then
+compare changes over subsequent timesteps.

docs/source/usersguide/geometry.rst

@@ -182,6 +182,8 @@ boundary condition.
 
 Periodic boundary conditions can be applied to pairs of planar surfaces.
 If there are only two periodic surfaces they will be matched automatically.
+
+
 Otherwise it is necessary to specify pairs explicitly using the
 :attr:`Surface.periodic_surface` attribute as in the following example::
 
@@ -192,7 +194,7 @@ Otherwise it is necessary to specify pairs explicitly using the
 Both rotational and translational periodic boundary conditions are specified in
 the same fashion. If both planes have the same normal vector, a translational
 periodicity is assumed; rotational periodicity is assumed otherwise. Currently,
-only rotations about the :math:`z`-axis are supported.
+rotations must be about the :math:`x`-, :math:`y`-, or :math:`z`-axis.
 
 For a rotational periodic BC, the normal vectors of each surface must point
 inwards---towards the valid geometry. For example, a :class:`XPlane` and
@@ -246,6 +248,28 @@ The classes :class:`Halfspace`, :class:`Intersection`, :class:`Union`, and
 :class:`Complement` and all instances of :class:`openmc.Region` and can be
 assigned to the :attr:`Cell.region` attribute.
 
+Cells also contain :attr:`Cell.temperature` and :attr:`Cell.density`
+attributes which override the temperature and density of the fill. These can
+be quite useful when temperatures and densities are spatially varying, as the
+alternative would be to add a unique :class:`Material` for each permutation of
+temperature, density, and composition. You can set the temperature (K) and
+density (g/cc) of a cell like so::
+
+  fuel.temperature = 800.0
+  fuel.density = 10.0
+
+The real utility of cell temperatures and densities occurs when a cell is
+replicated across the geometry, such as when a cell is the root geometric element
+in a replicated :ref:`universe<usersguide_universes>` or :ref:`lattice
+<usersguide_lattices>`. In those cases, you can provide a list of temperatures
+and densities to apply a temperature/density field to all of the distributed cells::
+
+  fuel.temperature = [800.0, 900.0, 800.0, 900.0]
+  fuel.density = [10.0, 9.0, 10.0, 9.0]
+
+In this example, the fuel cell is distributed four times in the geometry. Each
+distributed instance then receives its own temperature and density.
+
 .. _usersguide_universes:
 
 ---------
@@ -413,11 +437,11 @@ to help figure out how to place universes::
 
 
 Note that by default, hexagonal lattices are positioned such that each lattice
-element has two faces that are parallel to the :math:`y` axis. As one example,
-to create a three-ring lattice centered at the origin with a pitch of 10 cm
-where all the lattice elements centered along the :math:`y` axis are filled with
-universe ``u`` and the remainder are filled with universe ``q``, the following
-code would work::
+element has two faces that are perpendicular to the :math:`y` axis. As one
+example, to create a three-ring lattice centered at the origin with a pitch of
+10 cm where all the lattice elements centered along the :math:`y` axis are
+filled with universe ``u`` and the remainder are filled with universe ``q``, the
+following code would work::
 
   hexlat = openmc.HexLattice()
   hexlat.center = (0, 0)
@@ -530,6 +554,89 @@ UWUW and OpenMC material ID space will cause an error. To automatically resolve
 these ID overlaps, ``auto_ids`` can be set to ``True`` to append the UWUW
 material IDs to the OpenMC material ID space.
 
+
+Material overrides and differentiation
+--------------------------------------
+
+Programmatic access to DAGMC cell information for material overrides
+and differentiation requires synchronization of the DAGMC universe
+representation across Python and C-API::
+
+  model.init_lib()
+  model.sync_dagmc_universes()
+  model.finalize_lib()
+
+Upon completion of these steps, the :attr:`DAGMCUniverse.cells` attribute will
+be populated with :class:`DAGMCCell` proxy objects that represent the cells
+defined in the DAGMC model. The :class:`DAGMCCell` objects  will have
+:class:`openmc.Material`'s' applied according to the assignments upon
+initialization of the model. These materials can be replaced in the same manner
+as :class:`openmc.Cell` objects to override material assignments in the DAGMC
+model.
+
+Depletion with DAGMC geometry
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The synchronization of :class:`openmc.DAGMCUniverse`'s is important for
+depletion calculations using DAGMC geometry when materials need to be
+differentiated to perform material burnup independently in each DAGMC cell. See
+:meth:`openmc.model.Model.differentiate_mats`.
+
+Material overrides
+~~~~~~~~~~~~~~~~~~
+
+OpenMC supports overriding material assignments defined inside a DAGMC HDF5
+model so that CAD-assigned materials can be replaced by :class:`openmc.Material`
+objects. This is useful when the CAD geometry provides the shape but OpenMC
+materials (specific nuclide content, densities, or depletion behavior) are
+required.
+
+
+Replacing materials by name
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+If a DAGMC file includes material name tags, you can replace all cells that
+reference a particular name with an :class:`openmc.Material` using
+:meth:`~openmc.DAGMCUniverse.replace_material_assignment`::
+
+  import openmc
+
+  dag_univ = openmc.DAGMCUniverse('dagmc.h5m')
+
+  fuel = openmc.Material(name='fuel')
+  fuel.add_nuclide('U235', 0.05)
+  fuel.add_nuclide('U238', 0.95)
+  fuel.set_density('g/cm3', 10.5)
+
+  dag_univ.replace_material_assignment('Fuel', fuel)
+
+This lets you keep CAD geometry while adopting OpenMC material definitions.
+
+Per-cell material overrides
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+To assign overrides without initializing :class:`openmc.Model`, the
+:meth:`openmc.DAGMCUniverse.add_material_override` method can be used to assign
+materials to particular DAGMC cells. The method accepts either an integer cell
+ID::
+
+  dag_univ = openmc.DAGMCUniverse('dagmc.h5m')
+
+  enriched = openmc.Material(name='fuel_enriched')
+  enriched.add_nuclide('U235', 0.10)
+  enriched.add_nuclide('U238', 0.90)
+  enriched.set_density('g/cm3', 10.5)
+
+  dag_univ.add_material_override(1, enriched)
+
+In the case that the :class:`openmc.DAGMCUniverse` has already been synchronized,
+a :class:`openmc.DAGMCCell` object can also be provide to assign the material.
+
+Overrides are written to the `<material_overrides>` element of the
+:ref:`<dagmc_universe> <dagmc_element>` XML element so the C++ core can apply
+them on initialization.
+
+
 .. _Direct Accelerated Geometry Monte Carlo: https://svalinn.github.io/DAGMC/
 .. _University of Wisconsin Unified Workflow: https://svalinn.github.io/DAGMC/usersguide/uw2.html
 

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
 
@@ -383,6 +452,20 @@ OPENMC_USE_MPI
   options, please see the `FindMPI.cmake documentation
   <https://cmake.org/cmake/help/latest/module/FindMPI.html>`_.
 
+.. _cmake_strict_fp:
+
+OPENMC_ENABLE_STRICT_FP
+  Disables compiler optimizations that change floating-point results relative to
+  unoptimized builds, improving cross-platform and cross-optimization-level
+  reproducibility. This disables FMA contraction (``-ffp-contract=off``) and
+  compiler builtin replacements of math functions like ``pow``, ``exp``, ``log``
+  (``-fno-builtin``). It also keeps C/C++ assertions active by removing the
+  ``-DNDEBUG`` flag from ``RelWithDebInfo`` builds. Without this flag, these
+  optimizations can produce bit-level differences across platforms, compilers,
+  and optimization levels. This option should be used when running the test
+  suite. By default (off), the compiler is free to use all optimizations for
+  best performance. (Default: off)
+
 OPENMC_FORCE_VENDORED_LIBS
   Forces OpenMC to use the submodules located in the vendor directory, as
   opposed to searching the system for already installed versions of those
@@ -415,7 +498,10 @@ Release
 
 RelWithDebInfo
   (Default if no type is specified.) Enable optimization and debug. On most
-  platforms/compilers, this is equivalent to `-O2 -g`.
+  platforms/compilers, this is equivalent to `-O2 -g`. When
+  :ref:`OPENMC_ENABLE_STRICT_FP <cmake_strict_fp>` is enabled, OpenMC removes the
+  ``-DNDEBUG`` flag that CMake normally adds for this build type, so that
+  C/C++ assertions remain active.
 
 Example of configuring for Debug mode:
 

docs/source/usersguide/plots.rst

@@ -6,13 +6,14 @@ Geometry Visualization
 
 .. currentmodule:: openmc
 
-OpenMC is capable of producing two-dimensional slice plots of a geometry as well
-as three-dimensional voxel plots using the geometry plotting :ref:`run mode
-<usersguide_run_modes>`. The geometry plotting mode relies on the presence of a
-:ref:`plots.xml <io_plots>` file that indicates what plots should be created. To
-create this file, one needs to create one or more :class:`openmc.Plot`
-instances, add them to a :class:`openmc.Plots` collection, and then use the
-:class:`Plots.export_to_xml` method to write the ``plots.xml`` file.
+OpenMC is capable of producing two-dimensional slice plots of a geometry,
+three-dimensional voxel plots, and three-dimensional raytrace plots using the
+geometry plotting :ref:`run mode <usersguide_run_modes>`. The geometry plotting
+mode relies on the presence of a :ref:`plots.xml <io_plots>` file that indicates
+what plots should be created. To create this file, one needs to create one or
+more instances of the various plot classes described below, add them to a
+:class:`openmc.Plots` collection, and then use the :class:`Plots.export_to_xml`
+method to write the ``plots.xml`` file.
 
 -----------
 Slice Plots
@@ -21,15 +22,14 @@ Slice Plots
 .. image:: ../_images/atr.png
    :width: 300px
 
-By default, when an instance of :class:`openmc.Plot` is created, it indicates
-that a 2D slice plot should be made. You can specify the origin of the plot
-(:attr:`Plot.origin`), the width of the plot in each direction
-(:attr:`Plot.width`), the number of pixels to use in each direction
-(:attr:`Plot.pixels`), and the basis directions for the plot. For example, to
-create a :math:`x` - :math:`z` plot centered at (5.0, 2.0, 3.0) with a width of
-(50., 50.)  and 400x400 pixels::
+The :class:`openmc.SlicePlot` class indicates that a 2D slice plot should be
+made. You can specify the origin of the plot (:attr:`SlicePlot.origin`), the
+width of the plot in each direction (:attr:`SlicePlot.width`), the number of
+pixels to use in each direction (:attr:`SlicePlot.pixels`), and the basis
+directions for the plot. For example, to create a :math:`x` - :math:`z` plot
+centered at (5.0, 2.0, 3.0) with a width of (50., 50.)  and 400x400 pixels::
 
-  plot = openmc.Plot()
+  plot = openmc.SlicePlot()
   plot.basis = 'xz'
   plot.origin = (5.0, 2.0, 3.0)
   plot.width = (50., 50.)
@@ -47,7 +47,7 @@ that location.
 
 By default, a unique color will be assigned to each cell in the geometry. If you
 want your plot to be colored by material instead, change the
-:attr:`Plot.color_by` attribute::
+:attr:`SlicePlot.color_by` attribute::
 
   plot.color_by = 'material'
 
@@ -68,8 +68,8 @@ particular cells/materials should be given colors of your choosing::
 Note that colors can be given as RGB tuples or by a string indicating a valid
 `SVG color <https://www.w3.org/TR/SVG11/types.html#ColorKeywords>`_.
 
-When you're done creating your :class:`openmc.Plot` instances, you need to then
-assign them to a :class:`openmc.Plots` collection and export it to XML::
+When you're done creating your :class:`openmc.SlicePlot` instances, you need to
+then assign them to a :class:`openmc.Plots` collection and export it to XML::
 
   plots = openmc.Plots([plot1, plot2, plot3])
   plots.export_to_xml()
@@ -97,13 +97,11 @@ Voxel Plots
 .. image:: ../_images/3dba.png
    :width: 200px
 
-The :class:`openmc.Plot` class can also be told to generate a 3D voxel plot
-instead of a 2D slice plot. Simply change the :attr:`Plot.type` attribute to
-'voxel'. In this case, the :attr:`Plot.width` and :attr:`Plot.pixels` attributes
-should be three items long, e.g.::
+The :class:`openmc.VoxelPlot` class enables the generation of a 3D voxel plot
+instead of a 2D slice plot. In this case, the :attr:`VoxelPlot.width` and
+:attr:`VoxelPlot.pixels` attributes should be three items long, e.g.::
 
-  vox_plot = openmc.Plot()
-  vox_plot.type = 'voxel'
+  vox_plot = openmc.VoxelPlot()
   vox_plot.width = (100., 100., 50.)
   vox_plot.pixels = (400, 400, 200)
 

docs/source/usersguide/processing.rst

@@ -97,7 +97,15 @@ VTK Mesh File Generation
 ------------------------
 
 VTK files of OpenMC meshes can be created using the
-:meth:`openmc.Mesh.write_data_to_vtk`  method. Data can be applied to the
+:meth:`openmc.Mesh.write_data_to_vtk` method. This method supports several VTK
+formats depending on the mesh type. Structured meshes
+(:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`,
+:class:`~openmc.CylindricalMesh`, and :class:`~openmc.SphericalMesh`) can be
+exported to legacy VTK format (``.vtk``). The :class:`~openmc.UnstructuredMesh`
+class supports VTK unstructured grid formats (``.vtu``) as well as an HDF5-based
+format (``.vtkhdf``) that does not require the ``vtk`` module to write.
+
+Data can be applied to the
 elements of the resulting mesh from mesh filter objects. This data can be
 provided either as a flat array or, in the case of structured meshes
 (:class:`~openmc.RegularMesh`, :class:`~openmc.RectilinearMesh`,

docs/source/usersguide/random_ray.rst

@@ -513,6 +513,7 @@ Supported scores:
     - total
     - fission
     - nu-fission
+    - kappa-fission
     - events
 
 Supported Estimators:
@@ -644,7 +645,10 @@ model to use these multigroup cross sections. An example is given below::
       nparticles=2000,
       overwrite_mgxs_library=False,
       mgxs_path="mgxs.h5",
-      correction=None
+      correction=None,
+      source_energy=None,
+      temperatures=None,
+      temperature_settings=None
   )
 
 The most important parameter to set is the ``method`` parameter, which can be
@@ -706,6 +710,45 @@ generation and use an existing library file.
     with a :math:`\rho` default value of 1.0, which can be adjusted with the
     ``settings.random_ray['diagonal_stabilization_rho']`` parameter.
 
+When generating MGXS data with either the ``stochastic_slab`` or
+``infinite_medium`` methods, by default the simulation will use a uniform source
+distribution spread evenly over all energy groups. This ensures that all energy
+groups receive tallies and therefore produce non-zero total multigroup cross
+sections. Additionally, the function will convert any sources in the model into
+simplified spatial sources that retain the original energy distributions.  If
+sources are present, they will be used 99% of the time to sample source energies
+during MGXS generation. The other 1% of the time, energies will be sampled
+uniformly over all energy groups to ensure that all groups receive some tallies.
+However, the user may wish to specify a different source energy spectrum (for
+instance, if they are using a FileSource, such that the energy distribution
+cannot be extracted from the python source object). This can be done by
+providing a :class:`openmc.stats.Univariate` distribution as the
+``source_energy`` parameter of the :meth:`openmc.Model.convert_to_multigroup`
+method. If provided, it will override any sources present in the model and will
+be used 99% of the time to sample source energies during MGXS generation. The
+other 1% of the time, energies will be sampled uniformly over all energy groups
+to ensure that all groups receive some tallies.
+
+For instance, a D-D fusion simulation may involve a complex file source. In this
+case, the user may wish to provide a discrete 2.45 MeV energy source
+distribution for MGXS generation as::
+
+  source_energy = openmc.stats.delta_function(2.45e6)
+
+The ``temperatures`` parameter can be provided if temperature-dependent
+multi-group cross sections are desired for multi-physics simulations. An
+individual cross section generation calculation is run for each temperature
+provided, where the materials in the model are set to the temperature. The
+temperature settings used during cross section generation can be specified with the
+``temperature_settings`` parameter. If no ``temperature_settings`` are provided,
+the settings contained in the model will be used. The valid keys and values in the
+``temperature_settings`` dictionary are identical to
+:attr:`openmc.Settings.temperature_settings`; more information can be found in
+:class:`openmc.Settings` . This approach yields isothermal cross section interpolation
+tables, which can be inaccurate for systems with large differences between temperatures
+in each material (often the case in fission reactors). If a more sophisticated
+temperature-dependence is required, we recommend generating cross sections manually.
+
 Ultimately, the methods described above are all just approximations.
 Approximations in the generated MGXS data will fundamentally limit the potential
 accuracy of the random ray solver. However, the methods described above are all
@@ -901,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
 ---------------------------------
@@ -1030,22 +1075,52 @@ 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. When an initial forward solve is performed (i.e., when no
+user-specified adjoint source is present), its output files are also written to
+disk with a ``forward`` infix, so they are not overwritten by the subsequent
+adjoint solve. This applies to the statepoint, ``tallies.out``, and any voxel
+plots, e.g., ``statepoint.forward.N.h5`` and ``tallies.forward.out``; the
+adjoint solve keeps the usual file names. This allows analyses requiring both
+the forward and adjoint solutions to be performed from a single run. When
+generating FW-CADIS weight windows, no weight window file is written for the
+forward solve, as only the final adjoint-derived weight windows are meaningful.
 
 .. note::
-    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
@@ -1105,11 +1180,10 @@ given below:
     tallies.export_to_xml()
 
     # Create voxel plot
-    plot = openmc.Plot()
+    plot = openmc.VoxelPlot()
     plot.origin = [0, 0, 0]
     plot.width = [2*pitch, 2*pitch, 1]
     plot.pixels = [1000, 1000, 1]
-    plot.type = 'voxel'
 
     # Instantiate a Plots collection and export to XML
     plots = openmc.Plots([plot])
@@ -1189,11 +1263,10 @@ given below:
     tallies.export_to_xml()
 
     # Create voxel plot
-    plot = openmc.Plot()
+    plot = openmc.VoxelPlot()
     plot.origin = [0, 0, 0]
     plot.width = [2*pitch, 2*pitch, 1]
     plot.pixels = [1000, 1000, 1]
-    plot.type = 'voxel'
 
     # Instantiate a Plots collection and export to XML
     plots = openmc.Plots([plot])

docs/source/usersguide/scripts.rst

@@ -48,6 +48,7 @@ flags:
                        restart file
 -s, --threads N        Run with *N* OpenMP threads
 -t, --track            Write tracks for all particles (up to max_tracks)
+-q, --verbosity V      Set the output verbosity to *V*
 -v, --version          Show version information
 -h, --help             Show help message
 

docs/source/usersguide/settings.rst

@@ -272,6 +272,12 @@ option::
   settings.source = [src1, src2]
   settings.uniform_source_sampling = True
 
+Additionally, sampling from an :class:`openmc.IndependentSource` may be biased
+for local or global variance reduction by modifying the
+:attr:`~openmc.IndependentSource.bias` attribute of each of its four main
+distributions. Further discussion of source biasing can be found in
+:ref:`source_biasing`.
+
 Finally, the :attr:`IndependentSource.particle` attribute can be used to
 indicate the source should be composed of particles other than neutrons. For
 example, the following would generate a photon source::
@@ -285,6 +291,53 @@ example, the following would generate a photon source::
 For a full list of all classes related to statistical distributions, see
 :ref:`pythonapi_stats`.
 
+Tokamak Plasma Sources
+----------------------
+
+For fusion applications, the :class:`openmc.TokamakSource` class provides a
+native parametric neutron source for tokamak plasmas. Rather than specifying
+spatial, angular, and energy distributions separately, the source is defined by
+the plasma geometry (using a `Miller-style flux-surface parameterization
+<https://doi.org/10.1063/1.872666>`_) and a radial emission profile. Source
+sites are sampled directly from the plasma volume without rejection.
+
+The plasma shape is described by the major radius :math:`R_0`, minor radius
+:math:`a`, elongation :math:`\kappa`, triangularity :math:`\delta`, and
+Shafranov shift :math:`\Delta`. The neutron birth profile is given as an
+emission density :math:`S(r/a)` tabulated on a normalized minor-radius grid that
+runs from 0 (magnetic axis) to 1 (last closed flux surface); only the shape of
+the profile matters, since it is normalized internally. The emission density is
+linearly interpolated between the supplied points and refined internally for
+radial sampling. For example::
+
+  import numpy as np
+
+  r_over_a = np.linspace(0.0, 1.0, 50)
+  emission = (1.0 - r_over_a**2)**2  # peaked on-axis profile
+
+  source = openmc.TokamakSource(
+      major_radius=620.0,    # cm
+      minor_radius=200.0,    # cm
+      elongation=1.8,
+      triangularity=0.45,
+      shafranov_shift=10.0,  # cm
+      r_over_a=r_over_a,
+      emission_density=emission,
+      energy=openmc.stats.muir(e0=14.08e6, m_rat=5.0, kt=20000.0),
+  )
+
+  settings.source = source
+
+The ``energy`` argument accepts either a single
+:class:`~openmc.stats.Univariate` distribution applied at all radii, or a
+sequence with one distribution per ``r_over_a`` grid point to model a
+radially-varying neutron spectrum (energies are then sampled by stochastic
+interpolation between the two distributions bracketing the sampled radius). A
+time distribution can be given with the ``time`` argument; by default, particles
+are born at :math:`t=0`. The toroidal extent can be restricted with
+``phi_start`` and ``phi_extent`` to model a sector of the plasma, and
+``vertical_shift`` translates the plasma center along the z-axis.
+
 File-based Sources
 ------------------
 
@@ -394,7 +447,7 @@ below.
     {
       openmc::SourceSite particle;
       // weight
-      particle.particle = openmc::ParticleType::neutron;
+      particle.particle = openmc::ParticleType::neutron();
       particle.wgt = 1.0;
       // position
       double angle = 2.0 * M_PI * openmc::prn(seed);
@@ -471,7 +524,7 @@ parameters to the source class when it is created:
     {
       openmc::SourceSite particle;
       // weight
-      particle.particle = openmc::ParticleType::neutron;
+      particle.particle = openmc::ParticleType::neutron();
       particle.wgt = 1.0;
       // position
       particle.r.x = 0.0;
@@ -598,6 +651,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
@@ -779,6 +839,11 @@ collision_track.h5 file at the end of the simulation. The file contains
 300 recorded collisions that occurred in materials with IDs 1 or 2, involving
 fission or (n,2n) reactions on the nuclides U-238 or O-16, within cells
 with IDs 5 and 12.
+
+.. note::
+   Electron and positron collision-track events are not associated with a
+   specific nuclide. If a ``nuclides`` entry is specified, these events are omitted.
+
 The file can be read using :func:`openmc.read_collision_track_file`.
 The example below shows how to extract the data from the collision_track
 feature and displays the fields stored in the file:

docs/source/usersguide/tallies.rst

@@ -105,109 +105,124 @@ The following tables show all valid scores:
 
 .. table:: **Reaction scores: units are reactions per source particle.**
 
-    +----------------------+---------------------------------------------------+
-    |Score                 | Description                                       |
-    +======================+===================================================+
-    |absorption            |Total absorption rate. For incident neutrons, this |
-    |                      |accounts for all reactions that do not produce     |
-    |                      |secondary neutrons as well as fission. For incident|
-    |                      |photons, this includes photoelectric and pair      |
-    |                      |production.                                        |
-    +----------------------+---------------------------------------------------+
-    |elastic               |Elastic scattering reaction rate.                  |
-    +----------------------+---------------------------------------------------+
-    |fission               |Total fission reaction rate.                       |
-    +----------------------+---------------------------------------------------+
-    |scatter               |Total scattering rate.                             |
-    +----------------------+---------------------------------------------------+
-    |total                 |Total reaction rate.                               |
-    +----------------------+---------------------------------------------------+
-    |(n,2nd)               |(n,2nd) reaction rate.                             |
-    +----------------------+---------------------------------------------------+
-    |(n,2n)                |(n,2n) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,3n)                |(n,3n) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,na)                |(n,n\ :math:`\alpha`\ ) reaction rate.             |
-    +----------------------+---------------------------------------------------+
-    |(n,n3a)               |(n,n3\ :math:`\alpha`\ ) reaction rate.            |
-    +----------------------+---------------------------------------------------+
-    |(n,2na)               |(n,2n\ :math:`\alpha`\ ) reaction rate.            |
-    +----------------------+---------------------------------------------------+
-    |(n,3na)               |(n,3n\ :math:`\alpha`\ ) reaction rate.            |
-    +----------------------+---------------------------------------------------+
-    |(n,np)                |(n,np) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,n2a)               |(n,n2\ :math:`\alpha`\ ) reaction rate.            |
-    +----------------------+---------------------------------------------------+
-    |(n,2n2a)              |(n,2n2\ :math:`\alpha`\ ) reaction rate.           |
-    +----------------------+---------------------------------------------------+
-    |(n,nd)                |(n,nd) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,nt)                |(n,nt) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,n3He)              |(n,n\ :sup:`3`\ He) reaction rate.                 |
-    +----------------------+---------------------------------------------------+
-    |(n,nd2a)              |(n,nd2\ :math:`\alpha`\ ) reaction rate.           |
-    +----------------------+---------------------------------------------------+
-    |(n,nt2a)              |(n,nt2\ :math:`\alpha`\ ) reaction rate.           |
-    +----------------------+---------------------------------------------------+
-    |(n,4n)                |(n,4n) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,2np)               |(n,2np) reaction rate.                             |
-    +----------------------+---------------------------------------------------+
-    |(n,3np)               |(n,3np) reaction rate.                             |
-    +----------------------+---------------------------------------------------+
-    |(n,n2p)               |(n,n2p) reaction rate.                             |
-    +----------------------+---------------------------------------------------+
-    |(n,n*X*)              |Level inelastic scattering reaction rate. The *X*  |
-    |                      |indicates what which inelastic level, e.g., (n,n3) |
-    |                      |is third-level inelastic scattering.               |
-    +----------------------+---------------------------------------------------+
-    |(n,nc)                |Continuum level inelastic scattering reaction rate.|
-    +----------------------+---------------------------------------------------+
-    |(n,gamma)             |Radiative capture reaction rate.                   |
-    +----------------------+---------------------------------------------------+
-    |(n,p)                 |(n,p) reaction rate.                               |
-    +----------------------+---------------------------------------------------+
-    |(n,d)                 |(n,d) reaction rate.                               |
-    +----------------------+---------------------------------------------------+
-    |(n,t)                 |(n,t) reaction rate.                               |
-    +----------------------+---------------------------------------------------+
-    |(n,3He)               |(n,\ :sup:`3`\ He) reaction rate.                  |
-    +----------------------+---------------------------------------------------+
-    |(n,a)                 |(n,\ :math:`\alpha`\ ) reaction rate.              |
-    +----------------------+---------------------------------------------------+
-    |(n,2a)                |(n,2\ :math:`\alpha`\ ) reaction rate.             |
-    +----------------------+---------------------------------------------------+
-    |(n,3a)                |(n,3\ :math:`\alpha`\ ) reaction rate.             |
-    +----------------------+---------------------------------------------------+
-    |(n,2p)                |(n,2p) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,pa)                |(n,p\ :math:`\alpha`\ ) reaction rate.             |
-    +----------------------+---------------------------------------------------+
-    |(n,t2a)               |(n,t2\ :math:`\alpha`\ ) reaction rate.            |
-    +----------------------+---------------------------------------------------+
-    |(n,d2a)               |(n,d2\ :math:`\alpha`\ ) reaction rate.            |
-    +----------------------+---------------------------------------------------+
-    |(n,pd)                |(n,pd) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,pt)                |(n,pt) reaction rate.                              |
-    +----------------------+---------------------------------------------------+
-    |(n,da)                |(n,d\ :math:`\alpha`\ ) reaction rate.             |
-    +----------------------+---------------------------------------------------+
-    |coherent-scatter      |Coherent (Rayleigh) scattering reaction rate.      |
-    +----------------------+---------------------------------------------------+
-    |incoherent-scatter    |Incoherent (Compton) scattering reaction rate.     |
-    +----------------------+---------------------------------------------------+
-    |photoelectric         |Photoelectric absorption reaction rate.            |
-    +----------------------+---------------------------------------------------+
-    |pair-production       |Pair production reaction rate.                     |
-    +----------------------+---------------------------------------------------+
-    |*Arbitrary integer*   |An arbitrary integer is interpreted to mean the    |
-    |                      |reaction rate for a reaction with a given ENDF MT  |
-    |                      |number.                                            |
-    +----------------------+---------------------------------------------------+
+    +------------------------+-------------------------------------------------+
+    |Score                   |Description                                      |
+    +========================+=================================================+
+    |absorption              |Total absorption rate. For incident neutrons,    |
+    |                        |this accounts for all reactions that do not      |
+    |                        |produce secondary neutrons as well as fission.   |
+    |                        |For incident photons, this includes              |
+    |                        |photoelectric and pair production.               |
+    +------------------------+-------------------------------------------------+
+    |elastic                 |Elastic scattering reaction rate.                |
+    +------------------------+-------------------------------------------------+
+    |fission                 |Total fission reaction rate.                     |
+    +------------------------+-------------------------------------------------+
+    |scatter                 |Total scattering rate.                           |
+    +------------------------+-------------------------------------------------+
+    |total                   |Total reaction rate.                             |
+    +------------------------+-------------------------------------------------+
+    |(n,2nd)                 |(n,2nd) reaction rate.                           |
+    +------------------------+-------------------------------------------------+
+    |(n,2n)                  |(n,2n) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,3n)                  |(n,3n) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,na)                  |(n,n\ :math:`\alpha`\ ) reaction rate.           |
+    +------------------------+-------------------------------------------------+
+    |(n,n3a)                 |(n,n3\ :math:`\alpha`\ ) reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |(n,2na)                 |(n,2n\ :math:`\alpha`\ ) reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |(n,3na)                 |(n,3n\ :math:`\alpha`\ ) reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |(n,np)                  |(n,np) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,n2a)                 |(n,n2\ :math:`\alpha`\ ) reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |(n,2n2a)                |(n,2n2\ :math:`\alpha`\ ) reaction rate.         |
+    +------------------------+-------------------------------------------------+
+    |(n,nd)                  |(n,nd) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,nt)                  |(n,nt) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,n3He)                |(n,n\ :sup:`3`\ He) reaction rate.               |
+    +------------------------+-------------------------------------------------+
+    |(n,nd2a)                |(n,nd2\ :math:`\alpha`\ ) reaction rate.         |
+    +------------------------+-------------------------------------------------+
+    |(n,nt2a)                |(n,nt2\ :math:`\alpha`\ ) reaction rate.         |
+    +------------------------+-------------------------------------------------+
+    |(n,4n)                  |(n,4n) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,2np)                 |(n,2np) reaction rate.                           |
+    +------------------------+-------------------------------------------------+
+    |(n,3np)                 |(n,3np) reaction rate.                           |
+    +------------------------+-------------------------------------------------+
+    |(n,n2p)                 |(n,n2p) reaction rate.                           |
+    +------------------------+-------------------------------------------------+
+    |(n,npa)                 |(n,np\ :math:`\alpha`\ ) reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |(n,n*X*)                |Level inelastic scattering reaction rate. The    |
+    |                        |*X* indicates which inelastic level, e.g.,       |
+    |                        |(n,n3) is third-level inelastic scattering.      |
+    +------------------------+-------------------------------------------------+
+    |(n,nc)                  |Continuum level inelastic scattering             |
+    |                        |reaction rate.                                   |
+    +------------------------+-------------------------------------------------+
+    |(n,gamma)               |Radiative capture reaction rate.                 |
+    +------------------------+-------------------------------------------------+
+    |(n,p)                   |(n,p) reaction rate.                             |
+    +------------------------+-------------------------------------------------+
+    |(n,d)                   |(n,d) reaction rate.                             |
+    +------------------------+-------------------------------------------------+
+    |(n,t)                   |(n,t) reaction rate.                             |
+    +------------------------+-------------------------------------------------+
+    |(n,3He)                 |(n,\ :sup:`3`\ He) reaction rate.                |
+    +------------------------+-------------------------------------------------+
+    |(n,a)                   |(n,\ :math:`\alpha`\ ) reaction rate.            |
+    +------------------------+-------------------------------------------------+
+    |(n,2a)                  |(n,2\ :math:`\alpha`\ ) reaction rate.           |
+    +------------------------+-------------------------------------------------+
+    |(n,3a)                  |(n,3\ :math:`\alpha`\ ) reaction rate.           |
+    +------------------------+-------------------------------------------------+
+    |(n,2p)                  |(n,2p) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,pa)                  |(n,p\ :math:`\alpha`\ ) reaction rate.           |
+    +------------------------+-------------------------------------------------+
+    |(n,t2a)                 |(n,t2\ :math:`\alpha`\ ) reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |(n,d2a)                 |(n,d2\ :math:`\alpha`\ ) reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |(n,pd)                  |(n,pd) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,pt)                  |(n,pt) reaction rate.                            |
+    +------------------------+-------------------------------------------------+
+    |(n,da)                  |(n,d\ :math:`\alpha`\ ) reaction rate.           |
+    +------------------------+-------------------------------------------------+
+    |photon-total            |Total photo-atomic reaction rate.                |
+    +------------------------+-------------------------------------------------+
+    |coherent-scatter        |Coherent (Rayleigh) scattering reaction rate.    |
+    +------------------------+-------------------------------------------------+
+    |incoherent-scatter      |Incoherent (Compton) scattering reaction rate.   |
+    +------------------------+-------------------------------------------------+
+    |photoelectric           |Photoelectric absorption reaction rate.          |
+    +------------------------+-------------------------------------------------+
+    |photoelectric-*S*       |Subshell photoelectric absorption rate for the   |
+    |                        |*S* shell. For example, "photoelectric-N3" is the|
+    |                        |rate for the N3 subshell.                        |
+    +------------------------+-------------------------------------------------+
+    |pair-production         |Pair production reaction rate (total).           |
+    +------------------------+-------------------------------------------------+
+    |pair-production-electron|Pair production reaction rate in the electron    |
+    |                        |field.                                           |
+    +------------------------+-------------------------------------------------+
+    |pair-production-nuclear |Pair production reaction rate in the nuclear     |
+    |                        |field.                                           |
+    +------------------------+-------------------------------------------------+
+    |*Arbitrary integer*     |An arbitrary integer is interpreted to mean the  |
+    |                        |reaction rate for a reaction with a given ENDF   |
+    |                        |MT number.                                       |
+    +------------------------+-------------------------------------------------+
 
 .. table:: **Particle production scores: units are particles produced per
            source particles.**
@@ -246,18 +261,21 @@ The following tables show all valid scores:
     +----------------------+---------------------------------------------------+
     |Score                 | Description                                       |
     +======================+===================================================+
-    |current               |Used in combination with a meshsurface filter:     |
+    |current               |It may not be used in conjunction with any other   |
+    |                      |score except flux.                                 |
+    |                      |                                                   |
+    |                      |When used in combination with a meshsurface filter:|
     |                      |Partial currents on the boundaries of each cell in |
-    |                      |a mesh. It may not be used in conjunction with any |
-    |                      |other score. Only energy and mesh filters may be   |
-    |                      |used.                                              |
-    |                      |Used in combination with a surface filter:         |
+    |                      |a mesh.                                            |
+    |                      |                                                   |
+    |                      |When used in combination with a surface filter:    |
     |                      |Net currents on any surface previously defined in  |
     |                      |the geometry. It may be used along with any other  |
     |                      |filter, except meshsurface filters.                |
     |                      |Surfaces can alternatively be defined with cell    |
     |                      |from and cell filters thereby resulting in tallying|
     |                      |partial currents.                                  |
+    |                      |                                                   |
     |                      |Units are particles per source particle.           |
     +----------------------+---------------------------------------------------+
     |events                |Number of scoring events. Units are events per     |

docs/source/usersguide/variance_reduction.rst

@@ -4,24 +4,27 @@
 Variance Reduction
 ==================
 
-Global variance reduction in OpenMC is accomplished by weight windowing
-techniques. 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 this section, we break down the
-steps required to both generate and then apply weight windows.
+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`
@@ -69,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::
@@ -88,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::
@@ -146,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
@@ -172,3 +226,148 @@ Weight window mesh information is embedded into the weight window file, so the
 mesh does not need to be redefined. Monte Carlo solves that load a weight window
 file as above will utilize weight windows to reduce the variance of the
 simulation.
+
+.. _source_biasing:
+
+--------------
+Source Biasing
+--------------
+
+In fixed source problems, source biasing provides a means to reduce the variance
+on global or localized responses, depending on the biasing scheme. In either
+case, the premise of the method is to sample source sites from a biased
+distribution that directs a larger fraction of the simulated histories towards
+phase space regions of interest than would be found there under analog sampling.
+In order to preserve an unbiased estimate of the tally mean, the weight of these
+with analog sampling, divided by the probability assigned by the biased
+distribution. While the assignment of statistical weights is outlined in the
+:ref:`methods section <methods_source_biasing>`, this section demonstrates the
+implementation of source biasing to problems in OpenMC.
+
+Source biasing in OpenMC is accomplished by applying a distribution to the
+:attr:`bias` attribute of one or more of the univariate or independent
+multivariate distributions which make up an :class:`~openmc.IndependentSource`
+instance as follows::
+
+    # First create the biased distribution
+    biased_dist = openmc.stats.PowerLaw(a=0, b=3, n=3)
+
+    # Construct a new distribution with the bias applied
+    dist = openmc.stats.PowerLaw(a=0, b=3, n=2, bias=biased_dist)
+
+    # The bias attribute can also be set on an existing "analog" distribution:
+    sphere_dist = openmc.stats.spherical_uniform(r_outer=3)
+    sphere_dist.r.bias = biased_dist
+
+Univariate distributions may be sampled via the Python API, returning the
+sample(s) along with the associated weight(s)::
+
+    sample_vec, wgt_vec = dist.sample(n_samples=100)
+
+Here, if the distribution is unbiased, the weight of each sample will be unity.
+Finally, :class:`~openmc.IndependentSource` instances can be constructed with
+biased distributions::
+
+    # Create a source with a biased spatial distribution
+    source = openmc.IndependentSource(space=sphere_dist)
+
+During the simulation, source sites are then sampled using the biased
+distributions where available and given starting statistical weights
+corresponding to the cumulative product of the weights assigned by each
+distribution in the source object. Hence multiple source variables (e.g.,
+direction and energy) may be biased and the resulting source sites will have
+their weights adjusted accordingly.
+
+.. note::
+    Combining source biasing with weight windows can be a powerful variance
+    reduction technique if each is constructed appropriately for the response
+    of interest. For example, if a source biasing scheme is devised for
+    variance reduction of a specific localized response, the user may be able
+    to specify their own weight window structure that results in more efficient
+    transport than if weight windows were generated by either of OpenMC's
+    automatic weight window generators, which are intended for global variance
+    reduction.
+
+Biased distributions that could result in degenerate weight mappings are not
+recommended; this is most commonly seen when biasing the :math:`\phi`-coordinate
+of spherical or cylindrical independent multivariate distributions. In such
+cases degenerate behavior will be observed at the pole about which :math:`\phi`
+is measured, with all values of :math:`\phi` (hence many possible statistical
+weights) mapping to the same point for :math:`r=0` or :math:`\mu=0`, and large
+weight gradients in the vicinity. In most cases requiring a spherical
+independent source, it would be preferable to reorient the reference vector of
+the distribution such that biasing could be applied to the
+:math:`\mu`-coordinate instead.
+
+When biasing a distribution, care should also be taken to ensure that both the
+unbiased and biased distribution share a common support---that is, every region
+of phase space mapped to a nonzero probability density by the unbiased
+distribution should likewise map to nonzero probability under the biased
+distribution, and vice versa. In OpenMC, this places restrictions on the set of
+compatible distributions that may be used to bias sampling of each distribution
+type. The following table summarizes the method for each distribution in OpenMC
+that permits biased sampling.
+
+.. list-table:: **Distributions that support biased sampling**
+   :header-rows: 1
+   :widths: 35 65
+
+   * - Discrete Univariate PDFs
+     - Biasing Method
+   * - :class:`openmc.stats.Discrete`
+     - Apply a vector of alternative probabilities to the :attr:`bias`
+       attribute
+
+.. list-table::
+   :header-rows: 1
+   :widths: 35 65
+
+   * - Continuous Univariate PDFs
+     - Biasing Method
+   * - :class:`openmc.stats.Uniform`,
+       :class:`openmc.stats.PowerLaw`,
+       :class:`openmc.stats.Maxwell`,
+       :class:`openmc.stats.Watt`,
+       :class:`openmc.stats.Normal`,
+       :class:`openmc.stats.Tabular`
+     - Apply a second, unbiased continous univariate PDF to the :attr:`bias`
+       attribute, ensuring that the :attr:`support` attribute of each
+       distribution is the same
+
+.. list-table::
+   :header-rows: 1
+   :widths: 35 65
+
+   * - Mixed Univariate PDFs
+     - Biasing Method
+   * - :class:`openmc.stats.Mixture`
+     - May be constructed from multiple biased univariate distributions, or a
+       second, unbiased continous univariate PDF may be applied to the
+       :attr:`bias` attribute
+
+.. list-table::
+   :header-rows: 1
+   :widths: 35 65
+
+   * - Discrete Multivariate PDFs
+     - Biasing Method
+   * - :class:`openmc.stats.PointCloud`,
+       :class:`openmc.stats.MeshSpatial`
+     - Apply a vector of the new relative probabilities of each point or mesh
+       element under biased sampling to the :attr:`bias` attribute
+
+.. list-table::
+   :header-rows: 1
+   :widths: 35 65
+
+   * - Continuous Multivariate PDFs
+     - Biasing Method
+   * - :class:`openmc.stats.CartesianIndependent`,
+       :class:`openmc.stats.CylindricalIndependent`,
+       :class:`openmc.stats.SphericalIndependent`,
+       :class:`openmc.stats.PolarAzimuthal`
+     - Construct from biased univariate distributions for :attr:`x`, :attr:`y`,
+       :attr:`z`, etc.
+   * - :class:`openmc.stats.Isotropic`
+     - Apply an unbiased :class:`openmc.stats.PolarAzimuthal` to the
+       :attr:`bias` attribute

examples/custom_source/source_ring.cpp

@@ -1,17 +1,17 @@
-#include <cmath> // for M_PI
+#define _USE_MATH_DEFINES
+#include <cmath>  // for M_PI
 #include <memory> // for unique_ptr
 
+#include "openmc/particle.h"
 #include "openmc/random_lcg.h"
 #include "openmc/source.h"
-#include "openmc/particle.h"
 
-class RingSource : public openmc::Source
-{
+class RingSource : public openmc::Source {
   openmc::SourceSite sample(uint64_t* seed) const
   {
     openmc::SourceSite particle;
     // particle type
-    particle.particle = openmc::ParticleType::neutron;
+    particle.particle = openmc::ParticleType::neutron();
     // position
     double angle = 2.0 * M_PI * openmc::prn(seed);
     double radius = 3.0;
@@ -25,10 +25,11 @@ class RingSource : public openmc::Source
   }
 };
 
-// A function to create a unique pointer to an instance of this class when generated
-// via a plugin call using dlopen/dlsym.
-// You must have external C linkage here otherwise dlopen will not find the file
-extern "C" std::unique_ptr<RingSource> openmc_create_source(std::string parameters)
+// A function to create a unique pointer to an instance of this class when
+// generated via a plugin call using dlopen/dlsym. You must have external C
+// linkage here otherwise dlopen will not find the file
+extern "C" std::unique_ptr<RingSource> openmc_create_source(
+  std::string parameters)
 {
   return std::make_unique<RingSource>();
 }

examples/lattice/hexagonal/build_xml.py

@@ -128,14 +128,14 @@ settings_file.export_to_xml()
 #                   Exporting to OpenMC plots.xml file
 ###############################################################################
 
-plot_xy = openmc.Plot(plot_id=1)
+plot_xy = openmc.SlicePlot(plot_id=1)
 plot_xy.filename = 'plot_xy'
 plot_xy.origin = [0, 0, 0]
 plot_xy.width = [6, 6]
 plot_xy.pixels = [400, 400]
 plot_xy.color_by = 'material'
 
-plot_yz = openmc.Plot(plot_id=2)
+plot_yz = openmc.SlicePlot(plot_id=2)
 plot_yz.filename = 'plot_yz'
 plot_yz.basis = 'yz'
 plot_yz.origin = [0, 0, 0]

examples/lattice/nested/build_xml.py

@@ -135,7 +135,7 @@ settings_file.export_to_xml()
 #                   Exporting to OpenMC plots.xml file
 ###############################################################################
 
-plot = openmc.Plot(plot_id=1)
+plot = openmc.SlicePlot(plot_id=1)
 plot.origin = [0, 0, 0]
 plot.width = [4, 4]
 plot.pixels = [400, 400]

examples/lattice/simple/build_xml.py

@@ -128,7 +128,7 @@ settings_file.export_to_xml()
 #                   Exporting to OpenMC plots.xml file
 ###############################################################################
 
-plot = openmc.Plot(plot_id=1)
+plot = openmc.SlicePlot(plot_id=1)
 plot.origin = [0, 0, 0]
 plot.width = [4, 4]
 plot.pixels = [400, 400]

examples/parameterized_custom_source/parameterized_source_ring.cpp

@@ -1,63 +1,66 @@
-#include <cmath> // for M_PI
-#include <memory> // for unique_ptr
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <memory>
 #include <unordered_map>
 
+#include "openmc/particle.h"
 #include "openmc/random_lcg.h"
 #include "openmc/source.h"
-#include "openmc/particle.h"
 
 class RingSource : public openmc::Source {
-  public:
-    RingSource(double radius, double energy) : radius_(radius), energy_(energy) { }
+public:
+  RingSource(double radius, double energy) : radius_(radius), energy_(energy) {}
 
-    // Defines a function that can create a unique pointer to a new instance of this class
-    // by extracting the parameters from the provided string.
-    static std::unique_ptr<RingSource> from_string(std::string parameters)
-    {
-      std::unordered_map<std::string, std::string> parameter_mapping;
+  // Defines a function that can create a unique pointer to a new instance of
+  // this class by extracting the parameters from the provided string.
+  static std::unique_ptr<RingSource> from_string(std::string parameters)
+  {
+    std::unordered_map<std::string, std::string> parameter_mapping;
 
-      std::stringstream ss(parameters);
-      std::string parameter;
-      while (std::getline(ss, parameter, ',')) {
-        parameter.erase(0, parameter.find_first_not_of(' '));
-        std::string key = parameter.substr(0, parameter.find_first_of('='));
-        std::string value = parameter.substr(parameter.find_first_of('=') + 1, parameter.length());
-        parameter_mapping[key] = value;
-      }
-
-      double radius = std::stod(parameter_mapping["radius"]);
-      double energy = std::stod(parameter_mapping["energy"]);
-      return std::make_unique<RingSource>(radius, energy);
+    std::stringstream ss(parameters);
+    std::string parameter;
+    while (std::getline(ss, parameter, ',')) {
+      parameter.erase(0, parameter.find_first_not_of(' '));
+      std::string key = parameter.substr(0, parameter.find_first_of('='));
+      std::string value =
+        parameter.substr(parameter.find_first_of('=') + 1, parameter.length());
+      parameter_mapping[key] = value;
     }
 
-    // Samples from an instance of this class.
-    openmc::SourceSite sample(uint64_t* seed) const
-    {
-      openmc::SourceSite particle;
-      // particle type
-      particle.particle = openmc::ParticleType::neutron;
-      // position
-      double angle = 2.0 * M_PI * openmc::prn(seed);
-      double radius = this->radius_;
-      particle.r.x = radius * std::cos(angle);
-      particle.r.y = radius * std::sin(angle);
-      particle.r.z = 0.0;
-      // angle
-      particle.u = {1.0, 0.0, 0.0};
-      particle.E = this->energy_;
+    double radius = std::stod(parameter_mapping["radius"]);
+    double energy = std::stod(parameter_mapping["energy"]);
+    return std::make_unique<RingSource>(radius, energy);
+  }
 
-      return particle;
-    }
+  // Samples from an instance of this class.
+  openmc::SourceSite sample(uint64_t* seed) const
+  {
+    openmc::SourceSite particle;
+    // particle type
+    particle.particle = openmc::ParticleType::neutron();
+    // position
+    double angle = 2.0 * M_PI * openmc::prn(seed);
+    double radius = this->radius_;
+    particle.r.x = radius * std::cos(angle);
+    particle.r.y = radius * std::sin(angle);
+    particle.r.z = 0.0;
+    // angle
+    particle.u = {1.0, 0.0, 0.0};
+    particle.E = this->energy_;
 
-  private:
-    double radius_;
-    double energy_;
+    return particle;
+  }
+
+private:
+  double radius_;
+  double energy_;
 };
 
-// A function to create a unique pointer to an instance of this class when generated
-// via a plugin call using dlopen/dlsym.
-// You must have external C linkage here otherwise dlopen will not find the file
-extern "C" std::unique_ptr<RingSource> openmc_create_source(std::string parameters)
+// A function to create a unique pointer to an instance of this class when
+// generated via a plugin call using dlopen/dlsym. You must have external C
+// linkage here otherwise dlopen will not find the file
+extern "C" std::unique_ptr<RingSource> openmc_create_source(
+  std::string parameters)
 {
   return RingSource::from_string(parameters);
 }

examples/pincell_random_ray/build_xml.py

@@ -192,11 +192,10 @@ tallies.export_to_xml()
 #                   Exporting to OpenMC plots.xml file
 ###############################################################################
 
-plot = openmc.Plot()
+plot = openmc.VoxelPlot()
 plot.origin = [0, 0, 0]
 plot.width = [pitch, pitch, pitch]
 plot.pixels = [1000, 1000, 1]
-plot.type = 'voxel'
 
 # Instantiate a Plots collection and export to XML
 plots = openmc.Plots([plot])

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/atomic_mass.h

@@ -0,0 +1,28 @@
+//==============================================================================
+// atomic masses definitions
+//==============================================================================
+
+#ifndef OPENMC_ATOMIC_MASS_H
+#define OPENMC_ATOMIC_MASS_H
+
+#include <cstdint>
+#include <unordered_map>
+
+namespace openmc {
+
+// Values here are from the Committee on Data for Science and Technology
+// (CODATA) 2018 recommendation (https://physics.nist.gov/cuu/Constants/).
+
+// Physical constants
+constexpr double MASS_ELECTRON {5.48579909065e-4}; // mass of an electron in amu
+constexpr double MASS_NEUTRON {1.00866491595};     // mass of a neutron in amu
+constexpr double MASS_PROTON {1.007276466621};     // mass of a proton in amu
+constexpr double MASS_DEUTRON {2.013553212745};    // mass of a deutron in amu
+constexpr double MASS_HELION {3.014932247175};     // mass of a helion in amu
+constexpr double MASS_ALPHA {4.001506179127};      // mass of an alpha in amu
+
+extern std::unordered_map<int32_t, double> ATOMIC_MASS;
+
+} // namespace openmc
+
+#endif //  OPENMC_ATOMIC_MASS_H

include/openmc/bank.h

@@ -34,18 +34,24 @@ extern vector<vector<double>> ifp_fission_lifetime_bank;
 
 extern vector<int64_t> progeny_per_particle;
 
+extern SharedArray<SourceSite> shared_secondary_bank_read;
+extern SharedArray<SourceSite> shared_secondary_bank_write;
+
 } // namespace simulation
 
 //==============================================================================
 // Non-member functions
 //==============================================================================
 
-void sort_fission_bank();
+void sort_bank(SharedArray<SourceSite>& bank, bool is_fission_bank);
 
 void free_memory_bank();
 
 void init_fission_bank(int64_t max);
 
+int64_t synchronize_global_secondary_bank(
+  SharedArray<SourceSite>& shared_secondary_bank);
+
 } // namespace openmc
 
 #endif // OPENMC_BANK_H

include/openmc/bank_io.h

@@ -16,8 +16,9 @@
 namespace openmc {
 
 template<typename SiteType>
-void write_bank_dataset(const char* dataset_name, hid_t group_id,
-  span<SiteType> bank, const vector<int64_t>& bank_index, hid_t banktype
+void write_bank_dataset(
+  const char* dataset_name, hid_t group_id, span<SiteType> bank,
+  const vector<int64_t>& bank_index, hid_t membanktype, hid_t filebanktype
 #ifdef OPENMC_MPI
   ,
   MPI_Datatype mpi_dtype
@@ -30,8 +31,8 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
 #ifdef PHDF5
   hsize_t dims[] {static_cast<hsize_t>(dims_size)};
   hid_t dspace = H5Screate_simple(1, dims, nullptr);
-  hid_t dset = H5Dcreate(group_id, dataset_name, banktype, dspace, H5P_DEFAULT,
-    H5P_DEFAULT, H5P_DEFAULT);
+  hid_t dset = H5Dcreate(group_id, dataset_name, filebanktype, dspace,
+    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
   hsize_t count[] {static_cast<hsize_t>(count_size)};
   hid_t memspace = H5Screate_simple(1, count, nullptr);
@@ -42,7 +43,7 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
   hid_t plist = H5Pcreate(H5P_DATASET_XFER);
   H5Pset_dxpl_mpio(plist, H5FD_MPIO_COLLECTIVE);
 
-  H5Dwrite(dset, banktype, memspace, dspace, plist, bank.data());
+  H5Dwrite(dset, membanktype, memspace, dspace, plist, bank.data());
 
   H5Sclose(dspace);
   H5Sclose(memspace);
@@ -52,7 +53,7 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
   if (mpi::master) {
     hsize_t dims[] {static_cast<hsize_t>(dims_size)};
     hid_t dspace = H5Screate_simple(1, dims, nullptr);
-    hid_t dset = H5Dcreate(group_id, dataset_name, banktype, dspace,
+    hid_t dset = H5Dcreate(group_id, dataset_name, filebanktype, dspace,
       H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 
 #ifdef OPENMC_MPI
@@ -75,7 +76,8 @@ void write_bank_dataset(const char* dataset_name, hid_t group_id,
       H5Sselect_hyperslab(
         dspace_rank, H5S_SELECT_SET, start, nullptr, count, nullptr);
 
-      H5Dwrite(dset, banktype, memspace, dspace_rank, H5P_DEFAULT, bank.data());
+      H5Dwrite(
+        dset, membanktype, memspace, dspace_rank, H5P_DEFAULT, bank.data());
 
       H5Sclose(memspace);
       H5Sclose(dspace_rank);

include/openmc/boundary_condition.h

@@ -138,18 +138,28 @@ protected:
 //==============================================================================
 //! A BC that rotates particles about a global axis.
 //
-//! Currently only rotations about the z-axis are supported.
+//! Only rotations about the x, y, and z axes are supported.
 //==============================================================================
 
 class RotationalPeriodicBC : public PeriodicBC {
 public:
-  RotationalPeriodicBC(int i_surf, int j_surf);
-
+  enum PeriodicAxis { x, y, z };
+  RotationalPeriodicBC(int i_surf, int j_surf, PeriodicAxis axis);
+  double compute_periodic_rotation(
+    double rise_1, double run_1, double rise_2, double run_2) const;
   void handle_particle(Particle& p, const Surface& surf) const override;
 
 protected:
   //! Angle about the axis by which particle coordinates will be rotated
   double angle_;
+  //! Do we need to flip surfaces senses when applying the transformation?
+  bool flip_sense_;
+  //! Ensure that choice of axes is right handed. axis_1_idx_ corresponds to the
+  //! independent axis and axis_2_idx_ corresponds to the dependent axis in the
+  //! 2D plane  perpendicular to the planes' axis of rotation
+  int zero_axis_idx_;
+  int axis_1_idx_;
+  int axis_2_idx_;
 };
 
 } // namespace openmc

include/openmc/bounding_box.h

@@ -13,12 +13,19 @@ namespace openmc {
 //==============================================================================
 
 struct BoundingBox {
-  double xmin = -INFTY;
-  double xmax = INFTY;
-  double ymin = -INFTY;
-  double ymax = INFTY;
-  double zmin = -INFTY;
-  double zmax = INFTY;
+  Position min = {-INFTY, -INFTY, -INFTY};
+  Position max = {INFTY, INFTY, INFTY};
+
+  // Constructors
+  BoundingBox() = default;
+  BoundingBox(Position min_, Position max_) : min {min_}, max {max_} {}
+
+  // Static factory methods
+  static BoundingBox infinite() { return {}; }
+  static BoundingBox inverted()
+  {
+    return {{INFTY, INFTY, INFTY}, {-INFTY, -INFTY, -INFTY}};
+  }
 
   inline BoundingBox operator&(const BoundingBox& other)
   {
@@ -35,29 +42,26 @@ struct BoundingBox {
   // intersect operator
   inline BoundingBox& operator&=(const BoundingBox& other)
   {
-    xmin = std::max(xmin, other.xmin);
-    xmax = std::min(xmax, other.xmax);
-    ymin = std::max(ymin, other.ymin);
-    ymax = std::min(ymax, other.ymax);
-    zmin = std::max(zmin, other.zmin);
-    zmax = std::min(zmax, other.zmax);
+    min.x = std::max(min.x, other.min.x);
+    min.y = std::max(min.y, other.min.y);
+    min.z = std::max(min.z, other.min.z);
+    max.x = std::min(max.x, other.max.x);
+    max.y = std::min(max.y, other.max.y);
+    max.z = std::min(max.z, other.max.z);
     return *this;
   }
 
   // union operator
   inline BoundingBox& operator|=(const BoundingBox& other)
   {
-    xmin = std::min(xmin, other.xmin);
-    xmax = std::max(xmax, other.xmax);
-    ymin = std::min(ymin, other.ymin);
-    ymax = std::max(ymax, other.ymax);
-    zmin = std::min(zmin, other.zmin);
-    zmax = std::max(zmax, other.zmax);
+    min.x = std::min(min.x, other.min.x);
+    min.y = std::min(min.y, other.min.y);
+    min.z = std::min(min.z, other.min.z);
+    max.x = std::max(max.x, other.max.x);
+    max.y = std::max(max.y, other.max.y);
+    max.z = std::max(max.z, other.max.z);
     return *this;
   }
-
-  inline Position min() const { return {xmin, ymin, zmin}; }
-  inline Position max() const { return {xmax, ymax, zmax}; }
 };
 
 } // namespace openmc

include/openmc/bremsstrahlung.h

@@ -3,7 +3,7 @@
 
 #include "openmc/particle.h"
 
-#include "xtensor/xtensor.hpp"
+#include "openmc/tensor.h"
 
 namespace openmc {
 
@@ -14,9 +14,9 @@ namespace openmc {
 class BremsstrahlungData {
 public:
   // Data
-  xt::xtensor<double, 2> pdf;   //!< Bremsstrahlung energy PDF
-  xt::xtensor<double, 2> cdf;   //!< Bremsstrahlung energy CDF
-  xt::xtensor<double, 1> yield; //!< Photon yield
+  tensor::Tensor<double> pdf;   //!< Bremsstrahlung energy PDF
+  tensor::Tensor<double> cdf;   //!< Bremsstrahlung energy CDF
+  tensor::Tensor<double> yield; //!< Photon yield
 };
 
 class Bremsstrahlung {
@@ -32,9 +32,9 @@ public:
 
 namespace data {
 
-extern xt::xtensor<double, 1>
+extern tensor::Tensor<double>
   ttb_e_grid; //! energy T of incident electron in [eV]
-extern xt::xtensor<double, 1>
+extern tensor::Tensor<double>
   ttb_k_grid; //! reduced energy W/T of emitted photon
 
 } // namespace data

include/openmc/capi.h

@@ -9,14 +9,41 @@
 extern "C" {
 #endif
 
+//! Run a stochastic volume calculation
+//
+//! \return Status (negative if an error occurred)
 int openmc_calculate_volumes();
+
 int openmc_cell_filter_get_bins(
   int32_t index, const int32_t** cells, int32_t* n);
+
+//! Get the fill for a cell
+//
+//! \param index Index in the cells array
+//! \param type Type of the fill
+//! \param indices Array of material indices for cell
+//! \param n Length of indices array
+//! \return Status (negative if an error occurred)
 int openmc_cell_get_fill(
   int32_t index, int* type, int32_t** indices, int32_t* n);
+
+//! Get the ID of a cell
+//
+//! \param index Index in the cells array
+//! \param id ID of the cell
+//! \return Status (negative if an error occurred)
 int openmc_cell_get_id(int32_t index, int32_t* id);
+
+//! Get the temperature of a cell
+//
+//! \param index Index in the cells array
+//! \param instance Which instance of the cell. If a null pointer is
+//!                 passed, the temperature of the first instance is returned.
+//! \param T temperature of the cell
+//!\return Status (negative if an error occurred)
 int openmc_cell_get_temperature(
   int32_t index, const int32_t* instance, double* T);
+
 int openmc_cell_get_density(
   int32_t index, const int32_t* instance, double* rho);
 int openmc_cell_get_translation(int32_t index, double xyz[]);
@@ -116,15 +143,63 @@ int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
 int openmc_mesh_get_volumes(int32_t index, double* volumes);
 int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz,
-  int max_mats, int32_t* materials, double* volumes);
+  int max_mats, int32_t* materials, double* volumes, double* bboxes);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_new_filter(const char* type, int32_t* index);
 int openmc_next_batch(int* status);
 int openmc_nuclide_name(int index, const char** name);
 int openmc_plot_geometry();
+// Deprecated; use openmc_slice_data.
 int openmc_id_map(const void* slice, int32_t* data_out);
+// Deprecated; use openmc_slice_data.
 int openmc_property_map(const void* slice, double* data_out);
+int openmc_slice_data(const double origin[3], const double u_span[3],
+  const double v_span[3], const size_t pixels[2], bool show_overlaps, int level,
+  int32_t filter_index, int32_t* geom_data, double* property_data);
+int openmc_get_plot_index(int32_t id, int32_t* index);
+int openmc_plot_get_id(int32_t index, int32_t* id);
+int openmc_plot_set_id(int32_t index, int32_t id);
+int openmc_solidraytrace_plot_create(int32_t* index);
+int openmc_solidraytrace_plot_get_pixels(
+  int32_t index, int32_t* width, int32_t* height);
+int openmc_solidraytrace_plot_set_pixels(
+  int32_t index, int32_t width, int32_t height);
+int openmc_solidraytrace_plot_get_color_by(int32_t index, int32_t* color_by);
+int openmc_solidraytrace_plot_set_color_by(int32_t index, int32_t color_by);
+int openmc_solidraytrace_plot_set_default_colors(int32_t index);
+int openmc_solidraytrace_plot_set_all_opaque(int32_t index);
+int openmc_solidraytrace_plot_set_opaque(
+  int32_t index, int32_t id, bool visible);
+int openmc_solidraytrace_plot_set_color(
+  int32_t index, int32_t id, uint8_t r, uint8_t g, uint8_t b);
+int openmc_solidraytrace_plot_get_camera_position(
+  int32_t index, double* x, double* y, double* z);
+int openmc_solidraytrace_plot_set_camera_position(
+  int32_t index, double x, double y, double z);
+int openmc_solidraytrace_plot_get_look_at(
+  int32_t index, double* x, double* y, double* z);
+int openmc_solidraytrace_plot_set_look_at(
+  int32_t index, double x, double y, double z);
+int openmc_solidraytrace_plot_get_up(
+  int32_t index, double* x, double* y, double* z);
+int openmc_solidraytrace_plot_set_up(
+  int32_t index, double x, double y, double z);
+int openmc_solidraytrace_plot_get_light_position(
+  int32_t index, double* x, double* y, double* z);
+int openmc_solidraytrace_plot_set_light_position(
+  int32_t index, double x, double y, double z);
+int openmc_solidraytrace_plot_get_fov(int32_t index, double* fov);
+int openmc_solidraytrace_plot_set_fov(int32_t index, double fov);
+int openmc_solidraytrace_plot_update_view(int32_t index);
+int openmc_solidraytrace_plot_create_image(
+  int32_t index, uint8_t* data_out, int32_t width, int32_t height);
+int openmc_solidraytrace_plot_get_color(
+  int32_t index, int32_t id, uint8_t* r, uint8_t* g, uint8_t* b);
+int openmc_solidraytrace_plot_get_diffuse_fraction(
+  int32_t index, double* diffuse_fraction);
+int openmc_solidraytrace_plot_set_diffuse_fraction(
+  int32_t index, double diffuse_fraction);
 int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx,
   double** grid_y, int* ny, double** grid_z, int* nz);
 int openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x,
@@ -140,6 +215,7 @@ int openmc_remove_tally(int32_t index);
 int openmc_reset();
 int openmc_reset_timers();
 int openmc_run();
+void openmc_run_random_ray();
 int openmc_sample_external_source(size_t n, uint64_t* seed, void* sites);
 void openmc_set_seed(int64_t new_seed);
 void openmc_set_stride(uint64_t new_stride);
@@ -201,8 +277,8 @@ int openmc_weight_windows_set_energy_bounds(
   int32_t index, double* e_bounds, size_t e_bounds_size);
 int openmc_weight_windows_get_energy_bounds(
   int32_t index, const double** e_bounds, size_t* e_bounds_size);
-int openmc_weight_windows_set_particle(int32_t index, int particle);
-int openmc_weight_windows_get_particle(int32_t index, int* particle);
+int openmc_weight_windows_set_particle(int32_t index, int32_t particle);
+int openmc_weight_windows_get_particle(int32_t index, int32_t* particle);
 int openmc_weight_windows_get_bounds(int32_t index, const double** lower_bounds,
   const double** upper_bounds, size_t* size);
 int openmc_weight_windows_set_bounds(int32_t index, const double* lower_bounds,
@@ -218,6 +294,7 @@ int openmc_weight_windows_set_weight_cutoff(int32_t index, double cutoff);
 int openmc_weight_windows_get_max_split(int32_t index, int* max_split);
 int openmc_weight_windows_set_max_split(int32_t index, int max_split);
 size_t openmc_weight_windows_size();
+size_t openmc_plots_size();
 int openmc_weight_windows_export(const char* filename = nullptr);
 int openmc_weight_windows_import(const char* filename = nullptr);
 int openmc_zernike_filter_get_order(int32_t index, int* order);
@@ -227,10 +304,10 @@ int openmc_zernike_filter_set_order(int32_t index, int order);
 int openmc_zernike_filter_set_params(
   int32_t index, const double* x, const double* y, const double* r);
 
-int openmc_particle_filter_get_bins(int32_t idx, int bins[]);
+int openmc_particle_filter_get_bins(int32_t idx, int32_t bins[]);
 
 //! Sets the mesh and energy grid for CMFD reweight
-//! \param[in] meshtyally_id id of CMFD Mesh Tally
+//! \param[in] meshtally_id id of CMFD Mesh Tally
 //! \param[in] cmfd_indices indices storing spatial and energy dimensions of
 //! CMFD problem \param[in] norm CMFD normalization factor
 void openmc_initialize_mesh_egrid(

include/openmc/cell.h

@@ -123,11 +123,11 @@ private:
   //! BoundingBox if the particle is in a complex cell.
   BoundingBox bounding_box_complex(vector<int32_t> postfix) const;
 
-  //! Enfource precedence: Parenthases, Complement, Intersection, Union
-  void add_precedence();
+  //! Enforce precedence between intersections and unions
+  void enforce_precedence();
 
   //! Add parenthesis to enforce precedence
-  int64_t add_parentheses(int64_t start);
+  void add_parentheses(int64_t start);
 
   //! Remove complement operators from the expression
   void remove_complement_ops();
@@ -145,6 +145,30 @@ private:
   bool simple_; //!< Does the region contain only intersections?
 };
 
+//==============================================================================
+// XML parsing helpers for <cell> nodes
+//==============================================================================
+
+//! Parse material IDs from a <cell> XML node.
+//! \param node XML node containing a "material" attribute or child element
+//! \param cell_id Cell ID used in error messages
+//! \return Vector of material IDs (MATERIAL_VOID for "void")
+vector<int32_t> parse_cell_material_xml(pugi::xml_node node, int32_t cell_id);
+
+//! Parse temperatures in [K] from a <cell> XML node.
+//! Validates that all values are non-negative and the list is non-empty.
+//! \param node XML node containing a "temperature" attribute or child element
+//! \param cell_id Cell ID used in error messages
+//! \return Vector of temperatures in [K]
+vector<double> parse_cell_temperature_xml(pugi::xml_node node, int32_t cell_id);
+
+//! Parse densities in [g/cm³] from a <cell> XML node.
+//! Validates that all values are positive and the list is non-empty.
+//! \param node XML node containing a "density" attribute or child element
+//! \param cell_id Cell ID used in error messages
+//! \return Vector of densities in [g/cm³]
+vector<double> parse_cell_density_xml(pugi::xml_node node, int32_t cell_id);
+
 //==============================================================================
 
 class Cell {

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_;
 };
@@ -92,6 +101,8 @@ extern vector<unique_ptr<ChainNuclide>> chain_nuclides;
 
 void read_chain_file_xml();
 
+void free_memory_chain();
+
 } // namespace openmc
 
 #endif // OPENMC_CHAIN_H

include/openmc/constants.h

@@ -9,6 +9,7 @@
 #include <limits>
 
 #include "openmc/array.h"
+#include "openmc/atomic_mass.h"
 #include "openmc/vector.h"
 #include "openmc/version.h"
 
@@ -25,16 +26,16 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
 constexpr int HDF5_VERSION[] {3, 0};
 
 // Version numbers for binary files
-constexpr array<int, 2> VERSION_STATEPOINT {18, 1};
-constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
-constexpr array<int, 2> VERSION_TRACK {3, 0};
+constexpr array<int, 2> VERSION_STATEPOINT {18, 2};
+constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 1};
+constexpr array<int, 2> VERSION_TRACK {3, 1};
 constexpr array<int, 2> VERSION_SUMMARY {6, 1};
 constexpr array<int, 2> VERSION_VOLUME {1, 0};
 constexpr array<int, 2> VERSION_VOXEL {2, 0};
 constexpr array<int, 2> VERSION_MGXS_LIBRARY {1, 0};
 constexpr array<int, 2> VERSION_PROPERTIES {1, 1};
 constexpr array<int, 2> VERSION_WEIGHT_WINDOWS {1, 0};
-constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 0};
+constexpr array<int, 2> VERSION_COLLISION_TRACK {1, 2};
 
 // ============================================================================
 // ADJUSTABLE PARAMETERS
@@ -86,10 +87,10 @@ constexpr double INFTY {std::numeric_limits<double>::max()};
 // (CODATA) 2018 recommendation (https://physics.nist.gov/cuu/Constants/).
 
 // Physical constants
-constexpr double MASS_NEUTRON {1.00866491595}; // mass of a neutron in amu
+constexpr double AMU_EV {
+  9.3149410242e8}; // atomic mass unit energy equivalent in eV/c^2
 constexpr double MASS_NEUTRON_EV {
-  939.56542052e6};                             // mass of a neutron in eV/c^2
-constexpr double MASS_PROTON {1.007276466621}; // mass of a proton in amu
+  939.56542052e6}; // neutron mass energy equivalent in eV/c^2
 constexpr double MASS_ELECTRON_EV {
   0.51099895000e6}; // electron mass energy equivalent in eV/c^2
 constexpr double FINE_STRUCTURE {
@@ -226,6 +227,7 @@ enum ReactionType {
   N_XA = 207,
   HEATING = 301,
   DAMAGE_ENERGY = 444,
+  PHOTON_TOTAL = 501,
   COHERENT = 502,
   INCOHERENT = 504,
   PAIR_PROD_ELEC = 515,
@@ -301,7 +303,7 @@ enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION };
 enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB };
 
 // Tally score type -- if you change these, make sure you also update the
-// _SCORES dictionary in openmc/capi/tally.py
+// _SCORES dictionary in openmc/lib/tally.py
 //
 // These are kept as a normal enum and made negative, since variables which
 // store one of these enum values usually also may be responsible for storing
@@ -365,7 +367,8 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY };
 
 enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
 enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
-enum class RandomRaySampleMethod { PRNG, HALTON };
+enum class RandomRaySampleMethod { PRNG, HALTON, S2 };
+enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT };
 
 //==============================================================================
 // Geometry Constants

include/openmc/dagmc.h

@@ -94,6 +94,10 @@ private:
 class DAGUniverse : public Universe {
 
 public:
+  using MaterialOverrides = std::unordered_map<int32_t, vector<int32_t>>;
+  using TemperatureOverrides = std::unordered_map<int32_t, vector<double>>;
+  using DensityOverrides = std::unordered_map<int32_t, vector<double>>;
+
   explicit DAGUniverse(pugi::xml_node node);
 
   //! Create a new DAGMC universe
@@ -112,6 +116,9 @@ public:
   //! Initialize the DAGMC accel. data structures, indices, material
   //! assignments, etc.
   void initialize();
+  void initialize(const MaterialOverrides& material_overrides,
+    const TemperatureOverrides& temperature_overrides,
+    const DensityOverrides& density_overrides = {});
 
   //! Reads UWUW materials and returns an ID map
   void read_uwuw_materials();
@@ -146,7 +153,8 @@ public:
 
   //! Assign a material overriding normal assignement to a cell
   //! \param[in] c The OpenMC cell to which the material is assigned
-  void override_assign_material(std::unique_ptr<DAGCell>& c) const;
+  void override_assign_material(std::unique_ptr<DAGCell>& c,
+    const MaterialOverrides& material_overrides) const;
 
   //! Return the index into the model cells vector for a given DAGMC volume
   //! handle in the universe
@@ -187,7 +195,9 @@ private:
   void set_id();        //!< Deduce the universe id from model::universes
   void init_dagmc();    //!< Create and initialise DAGMC pointer
   void init_metadata(); //!< Create and initialise dagmcMetaData pointer
-  void init_geometry(); //!< Create cells and surfaces from DAGMC entities
+  void init_geometry(const MaterialOverrides& material_overrides,
+    const TemperatureOverrides& temperature_overrides,
+    const DensityOverrides& density_overrides);
 
   std::string
     filename_; //!< Name of the DAGMC file used to create this universe
@@ -201,11 +211,6 @@ private:
                              //!< generate new material IDs for the universe
   bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard"
                        //!< volume
-  std::unordered_map<int32_t, vector<int32_t>>
-    material_overrides_; //!< Map of material overrides
-                         //!< keys correspond to the DAGMCCell id
-                         //!< values are a list of material ids used
-                         //!< for the override
 };
 
 //==============================================================================

include/openmc/distribution.h

@@ -15,6 +15,17 @@
 
 namespace openmc {
 
+//==============================================================================
+// Helper function for computing importance weights from biased sampling
+//==============================================================================
+
+//! Compute importance weights for biased sampling
+//! \param p Unnormalized original probability vector
+//! \param b Unnormalized bias probability vector
+//! \return Vector of importance weights (p_norm[i] / b_norm[i])
+vector<double> compute_importance_weights(
+  const vector<double>& p, const vector<double>& b);
+
 //==============================================================================
 //! Abstract class representing a univariate probability distribution
 //==============================================================================
@@ -22,11 +33,41 @@ namespace openmc {
 class Distribution {
 public:
   virtual ~Distribution() = default;
-  virtual double sample(uint64_t* seed) const = 0;
+
+  //! Sample a value from the distribution, handling biasing automatically
+  //! \param seed Pseudorandom number seed pointer
+  //! \return (sampled value, importance weight)
+  virtual std::pair<double, double> sample(uint64_t* seed) const;
+
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x)
+  virtual double evaluate(double x) const;
 
   //! Return integral of distribution
   //! \return Integral of distribution
   virtual double integral() const { return 1.0; };
+
+  //! Set bias distribution
+  virtual void set_bias(std::unique_ptr<Distribution> bias)
+  {
+    bias_ = std::move(bias);
+  }
+
+  const Distribution* bias() const { return bias_.get(); }
+
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  virtual double sample_unbiased(uint64_t* seed) const = 0;
+
+  //! Read bias distribution from XML
+  //! \param node XML node that may contain a bias child element
+  void read_bias_from_xml(pugi::xml_node node);
+
+  // Biasing distribution
+  unique_ptr<Distribution> bias_;
 };
 
 using UPtrDist = unique_ptr<Distribution>;
@@ -50,7 +91,7 @@ public:
 
   //! Sample a value from the distribution
   //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
+  //! \return Sampled index
   size_t sample(uint64_t* seed) const;
 
   // Properties
@@ -67,7 +108,7 @@ private:
   //! Normalize distribution so that probabilities sum to unity
   void normalize();
 
-  //! Initialize alias tables for distribution
+  //! Initialize alias table for sampling
   void init_alias();
 };
 
@@ -82,20 +123,30 @@ public:
 
   //! Sample a value from the distribution
   //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! \return (sampled value, sample weight)
+  std::pair<double, double> sample(uint64_t* seed) const override;
 
   double integral() const override { return di_.integral(); };
 
+  //! Override set_bias as no-op (bias handled in constructor)
+  void set_bias(std::unique_ptr<Distribution> bias) override {}
+
   // Properties
   const vector<double>& x() const { return x_; }
   const vector<double>& prob() const { return di_.prob(); }
   const vector<size_t>& alias() const { return di_.alias(); }
+  const vector<double>& weight() const { return weight_; }
+
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
 
 private:
-  vector<double> x_; //!< Possible outcomes
-  DiscreteIndex di_; //!< discrete probability distribution of
-                     //!< outcome indices
+  vector<double> x_;      //!< Possible outcomes
+  vector<double> weight_; //!< Importance weights (empty if unbiased)
+  DiscreteIndex di_; //!< Discrete probability distribution of outcome indices
 };
 
 //==============================================================================
@@ -107,14 +158,20 @@ public:
   explicit Uniform(pugi::xml_node node);
   Uniform(double a, double b) : a_ {a}, b_ {b} {};
 
-  //! Sample a value from the distribution
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x)
+  double evaluate(double x) const override;
 
   double a() const { return a_; }
   double b() const { return b_; }
 
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
 private:
   double a_; //!< Lower bound of distribution
   double b_; //!< Upper bound of distribution
@@ -131,15 +188,21 @@ public:
     : offset_ {std::pow(a, n + 1)}, span_ {std::pow(b, n + 1) - offset_},
       ninv_ {1 / (n + 1)} {};
 
-  //! Sample a value from the distribution
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x)
+  double evaluate(double x) const override;
 
   double a() const { return std::pow(offset_, ninv_); }
   double b() const { return std::pow(offset_ + span_, ninv_); }
   double n() const { return 1 / ninv_ - 1; }
 
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
 private:
   //! Store processed values in object to allow for faster sampling
   double offset_; //!< a^(n+1)
@@ -156,13 +219,19 @@ public:
   explicit Maxwell(pugi::xml_node node);
   Maxwell(double theta) : theta_ {theta} {};
 
-  //! Sample a value from the distribution
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x)
+  double evaluate(double x) const override;
 
   double theta() const { return theta_; }
 
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
 private:
   double theta_; //!< Factor in exponential [eV]
 };
@@ -176,41 +245,66 @@ public:
   explicit Watt(pugi::xml_node node);
   Watt(double a, double b) : a_ {a}, b_ {b} {};
 
-  //! Sample a value from the distribution
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x)
+  double evaluate(double x) const override;
 
   double a() const { return a_; }
   double b() const { return b_; }
 
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
 private:
   double a_; //!< Factor in exponential [eV]
   double b_; //!< Factor in square root [1/eV]
 };
 
 //==============================================================================
-//! Normal distributions with form 1/2*std_dev*sqrt(pi) exp
-//! (-(e-E0)/2*std_dev)^2
+//! Normal distribution with optional truncation bounds.
+//!
+//! The standard normal PDF is 1/(sqrt(2*pi)*sigma) *
+//! exp(-(x-mu)^2/(2*sigma^2)). When truncated to [lower, upper], the PDF is
+//! renormalized so that it integrates to 1 over the truncation interval.
 //==============================================================================
 
 class Normal : public Distribution {
 public:
   explicit Normal(pugi::xml_node node);
-  Normal(double mean_value, double std_dev)
-    : mean_value_ {mean_value}, std_dev_ {std_dev} {};
+  Normal(double mean_value, double std_dev, double lower = -INFTY,
+    double upper = INFTY);
 
-  //! Sample a value from the distribution
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x), accounting for truncation normalization
+  double evaluate(double x) const override;
 
   double mean_value() const { return mean_value_; }
   double std_dev() const { return std_dev_; }
+  double lower() const { return lower_; }
+  double upper() const { return upper_; }
+  bool is_truncated() const { return is_truncated_; }
+
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
 
 private:
-  double mean_value_; //!< middle of distribution [eV]
-  double std_dev_;    //!< standard deviation [eV]
+  double mean_value_;  //!< Mean of distribution
+  double std_dev_;     //!< Standard deviation
+  double lower_;       //!< Lower truncation bound (default: -INFTY)
+  double upper_;       //!< Upper truncation bound (default: +INFTY)
+  bool is_truncated_;  //!< True if bounds are finite
+  double norm_factor_; //!< Normalization factor for truncated distribution
+
+  //! Compute normalization factor for truncated distribution
+  void compute_normalization();
 };
 
 //==============================================================================
@@ -223,10 +317,10 @@ public:
   Tabular(const double* x, const double* p, int n, Interpolation interp,
     const double* c = nullptr);
 
-  //! Sample a value from the distribution
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x)
+  double evaluate(double x) const override;
 
   // properties
   vector<double>& x() { return x_; }
@@ -235,6 +329,12 @@ public:
   Interpolation interp() const { return interp_; }
   double integral() const override { return integral_; };
 
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
 private:
   vector<double> x_;     //!< tabulated independent variable
   vector<double> p_;     //!< tabulated probability density
@@ -259,13 +359,19 @@ public:
   explicit Equiprobable(pugi::xml_node node);
   Equiprobable(const double* x, int n) : x_ {x, x + n} {};
 
-  //! Sample a value from the distribution
-  //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! Evaluate probability density, f(x), at a point
+  //! \param x Point to evaluate f(x)
+  //! \return f(x)
+  double evaluate(double x) const override;
 
   const vector<double>& x() const { return x_; }
 
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
 private:
   vector<double> x_; //! Possible outcomes
 };
@@ -280,18 +386,90 @@ public:
 
   //! Sample a value from the distribution
   //! \param seed Pseudorandom number seed pointer
-  //! \return Sampled value
-  double sample(uint64_t* seed) const override;
+  //! \return (sampled value, sample weight)
+  std::pair<double, double> sample(uint64_t* seed) const override;
 
   double integral() const override { return integral_; }
 
-private:
-  // Storrage for probability + distribution
-  using DistPair = std::pair<double, UPtrDist>;
+  //! Override set_bias as no-op (bias handled in constructor)
+  void set_bias(std::unique_ptr<Distribution> bias) override {}
 
-  vector<DistPair>
-    distribution_;  //!< sub-distributions + cummulative probabilities
-  double integral_; //!< integral of distribution
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
+private:
+  vector<UPtrDist> distribution_; //!< Sub-distributions
+  vector<double> weight_; //!< Importance weights for component selection
+  DiscreteIndex di_;      //!< Discrete probability distribution of indices
+  double integral_;       //!< Integral of distribution
+};
+
+//==============================================================================
+// DecaySpectrum — non-owning mixture of decay photon distributions
+//==============================================================================
+
+//! Energy distribution formed by mixing multiple decay photon spectra.
+//!
+//! Unlike the general Mixture distribution, this class holds non-owning
+//! pointers to the component distributions (which live in
+//! data::chain_nuclides). Each component is weighted by the activity
+//! (atoms * decay_constant) of the corresponding nuclide.
+
+class DecaySpectrum : public Distribution {
+public:
+  //============================================================================
+  // Types, aliases
+
+  struct Sample {
+    double energy;
+    double weight;
+    int parent_nuclide;
+  };
+
+  //============================================================================
+  // Constructors
+
+  //! Construct from an XML node containing nuclide names and atom densities.
+  //!
+  //! Reads child ``<nuclide>`` elements with ``name`` and ``density``
+  //! attributes, resolves them against the loaded depletion chain, and
+  //! constructs the mixed distribution.
+  explicit DecaySpectrum(pugi::xml_node node);
+
+  //============================================================================
+  // Methods
+
+  //! Sample a value from the distribution and return the parent nuclide index
+  //! \param seed Pseudorandom number seed pointer
+  //! \return (Sampled energy, sample weight, chain nuclide index)
+  Sample sample_with_parent(uint64_t* seed) const;
+
+  //! Sample a value from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return (sampled value, sample weight)
+  std::pair<double, double> sample(uint64_t* seed) const override;
+
+  double integral() const override;
+
+protected:
+  //! Sample a value (unbiased) from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled value
+  double sample_unbiased(uint64_t* seed) const override;
+
+private:
+  //! Initialize decay spectrum sampling data
+  //! \param nuclide_indices  Indices of decay photon emitters in
+  //!   data::chain_nuclides
+  //! \param atoms  Number of atoms for each component.
+  void init(vector<int> nuclide_indices, const vector<double>& atoms);
+
+  vector<int> nuclide_indices_; //!< Indices of emitting nuclides in the chain
+  DiscreteIndex di_;            //!< Discrete index for component selection
+  double integral_;             //!< Total photon emission rate
 };
 
 } // namespace openmc

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(); }