Merge branch 'develop' into IFP-inactive

.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

@@ -27,10 +27,10 @@ jobs:
       source_changed: ${{ steps.filter.outputs.source_changed }}
     steps:
       - name: Check out the repository
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
       - name: Examine changed files
         id: filter
-        uses: dorny/paths-filter@668c092af3649c4b664c54e4b704aa46782f6f7c # latest master commit, not released yet
+        uses: dorny/paths-filter@v4
         with:
           filters: |
             source_changed:
@@ -43,45 +43,42 @@ jobs:
     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 }}
@@ -89,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
@@ -106,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 }}
 
@@ -150,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: |
@@ -167,7 +158,7 @@ jobs:
           openmc -v
 
       - name: cache-xs
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: |
             ~/nndc_hdf5
@@ -222,6 +213,7 @@ jobs:
           parallel: true
           flag-name: C++ and Python
           path-to-lcov: coverage.lcov
+          fail-on-error: false
 
   coverage:
     needs: [filter-changes, main]
@@ -234,6 +226,7 @@ jobs:
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           parallel-finished: true
+          fail-on-error: false
 
   ci-pass:
     needs: [filter-changes, main, coverage]

.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

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

.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

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

AGENTS.md

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

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
@@ -388,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
@@ -426,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
@@ -496,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)
@@ -553,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

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

@@ -5,8 +5,6 @@ get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE)
 
 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@)
 endif()

docs/source/capi/index.rst

@@ -580,6 +580,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 +874,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

@@ -62,7 +62,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

docs/source/devguide/agentic-tools.rst

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

docs/source/devguide/index.rst

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

docs/source/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/depletion_results.rst

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

docs/source/io_formats/mgxs_library.rst

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

docs/source/io_formats/settings.rst

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

docs/source/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
@@ -321,6 +321,32 @@ should be set to:
   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/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

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

docs/source/pythonapi/base.rst

@@ -132,6 +132,7 @@ Constructing Tallies
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter
+   openmc.ParticleProductionFilter
    openmc.MuFilter
    openmc.MuSurfaceFilter
    openmc.PolarFilter
@@ -148,6 +149,7 @@ Constructing Tallies
    openmc.ZernikeRadialFilter
    openmc.ParentNuclideFilter
    openmc.ParticleFilter
+   openmc.ReactionFilter
    openmc.MeshMaterialVolumes
    openmc.Trigger
    openmc.TallyDerivative

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

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

docs/source/usersguide/geometry.rst

@@ -248,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:
 
 ---------

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

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

docs/source/usersguide/settings.rst

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

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

examples/parameterized_custom_source/parameterized_source_ring.cpp

@@ -1,5 +1,5 @@
-#include <cmath>  // for M_PI
-#include <memory> // for unique_ptr
+#include <cmath>
+#include <memory>
 #include <unordered_map>
 
 #include "openmc/particle.h"

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

@@ -125,6 +125,49 @@ int openmc_nuclide_name(int index, const char** name);
 int openmc_plot_geometry();
 int openmc_id_map(const void* slice, int32_t* data_out);
 int openmc_property_map(const void* slice, double* data_out);
+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 +183,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);
@@ -218,6 +262,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);

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"
 
@@ -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,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY };
 
 enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };
 enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
-enum class RandomRaySampleMethod { PRNG, HALTON };
+enum class RandomRaySampleMethod { PRNG, HALTON, S2 };
 
 //==============================================================================
 // Geometry Constants

include/openmc/distribution.h

@@ -265,23 +265,29 @@ private:
 };
 
 //==============================================================================
-//! 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);
 
   //! Evaluate probability density, f(x), at a point
   //! \param x Point to evaluate f(x)
-  //! \return 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
@@ -290,8 +296,15 @@ protected:
   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();
 };
 
 //==============================================================================
@@ -394,6 +407,71 @@ private:
   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
 
 #endif // OPENMC_DISTRIBUTION_H

include/openmc/distribution_angle.h

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

include/openmc/distribution_energy.h

@@ -5,7 +5,7 @@
 #define OPENMC_DISTRIBUTION_ENERGY_H
 
 #include "hdf5.h"
-#include "xtensor/xtensor.hpp"
+#include "openmc/tensor.h"
 
 #include "openmc/constants.h"
 #include "openmc/endf.h"
@@ -86,9 +86,9 @@ private:
   struct CTTable {
     Interpolation interpolation;  //!< Interpolation law
     int n_discrete;               //!< Number of of discrete energies
-    xt::xtensor<double, 1> e_out; //!< Outgoing energies in [eV]
-    xt::xtensor<double, 1> p;     //!< Probability density
-    xt::xtensor<double, 1> c;     //!< Cumulative distribution
+    tensor::Tensor<double> e_out; //!< Outgoing energies in [eV]
+    tensor::Tensor<double> p;     //!< Probability density
+    tensor::Tensor<double> c;     //!< Cumulative distribution
   };
 
   int n_region_;                        //!< Number of inteprolation regions

include/openmc/distribution_multi.h

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

include/openmc/distribution_spatial.h

@@ -67,12 +67,18 @@ public:
   Distribution* phi() const { return phi_.get(); }
   Distribution* z() const { return z_.get(); }
   Position origin() const { return origin_; }
+  Direction r_dir() const { return r_dir_; }
+  Direction phi_dir() const { return phi_dir_; }
+  Direction z_dir() const { return z_dir_; }
 
 private:
-  UPtrDist r_;      //!< Distribution of r coordinates
-  UPtrDist phi_;    //!< Distribution of phi coordinates
-  UPtrDist z_;      //!< Distribution of z coordinates
-  Position origin_; //!< Cartesian coordinates of the cylinder center
+  UPtrDist r_;        //!< Distribution of r coordinates
+  UPtrDist phi_;      //!< Distribution of phi coordinates
+  UPtrDist z_;        //!< Distribution of z coordinates
+  Position origin_;   //!< Cartesian coordinates of the cylinder center
+  Direction r_dir_;   //!< Direction of r-axis at phi=0
+  Direction phi_dir_; //!< Direction of phi-axis at phi=0
+  Direction z_dir_;   //!< Direction of z-axis
 };
 
 //==============================================================================

include/openmc/eigenvalue.h

@@ -6,7 +6,7 @@
 
 #include <cstdint> // for int64_t
 
-#include "xtensor/xtensor.hpp"
+#include "openmc/tensor.h"
 #include <hdf5.h>
 
 #include "openmc/array.h"
@@ -24,7 +24,7 @@ namespace simulation {
 extern double keff_generation; //!<  Single-generation k on each processor
 extern array<double, 2> k_sum; //!< Used to reduce sum and sum_sq
 extern vector<double> entropy; //!< Shannon entropy at each generation
-extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
+extern tensor::Tensor<double> source_frac; //!< Source fraction for UFS
 
 } // namespace simulation
 

include/openmc/endf.h

@@ -33,6 +33,13 @@ bool is_disappearance(int MT);
 //! \return Whether corresponding reaction is an inelastic scattering reaction
 bool is_inelastic_scatter(int MT);
 
+//! Determine whether an MT number matches a target MT, considering that the
+//! target may be a summation reaction.
+//! \param[in] event_mt MT number of the actual event
+//! \param[in] target_mt MT number to check against
+//! \return Whether event_mt is a component of target_mt (or equal to it)
+bool mt_matches(int event_mt, int target_mt);
+
 //==============================================================================
 //! Abstract one-dimensional function
 //==============================================================================

include/openmc/event.h

@@ -112,6 +112,19 @@ void process_collision_events();
 //! \param n_particles The number of particles in the particle buffer
 void process_death_events(int64_t n_particles);
 
+//! Process event queues until all are empty. Each iteration processes the
+//! longest queue first to maximize vectorization efficiency.
+void process_transport_events();
+
+//! Initialize secondary particles from a shared secondary bank for
+//! event-based transport
+//
+//! \param n_particles The number of particles to initialize
+//! \param offset The offset index in the shared secondary bank
+//! \param shared_secondary_bank The shared secondary bank to read from
+void process_init_secondary_events(int64_t n_particles, int64_t offset,
+  const SharedArray<SourceSite>& shared_secondary_bank);
+
 } // namespace openmc
 
 #endif // OPENMC_EVENT_H

include/openmc/hdf5_interface.h

@@ -11,8 +11,7 @@
 
 #include "hdf5.h"
 #include "hdf5_hl.h"
-#include "xtensor/xadapt.hpp"
-#include "xtensor/xarray.hpp"
+#include "openmc/tensor.h"
 
 #include "openmc/array.h"
 #include "openmc/error.h"
@@ -166,24 +165,19 @@ void read_attribute(hid_t obj_id, const char* name, vector<T>& vec)
   read_attr(obj_id, name, H5TypeMap<T>::type_id, vec.data());
 }
 
-// Generic array version
+// Tensor version
 template<typename T>
-void read_attribute(hid_t obj_id, const char* name, xt::xarray<T>& arr)
+void read_attribute(hid_t obj_id, const char* name, tensor::Tensor<T>& tensor)
 {
-  // Get shape of attribute array
+  // Get shape of attribute
   auto shape = attribute_shape(obj_id, name);
 
-  // Allocate new array to read data into
-  std::size_t size = 1;
-  for (const auto x : shape)
-    size *= x;
-  vector<T> buffer(size);
+  // Resize tensor and read data directly
+  vector<size_t> tshape(shape.begin(), shape.end());
+  tensor.resize(tshape);
 
   // Read data from attribute
-  read_attr(obj_id, name, H5TypeMap<T>::type_id, buffer.data());
-
-  // Adapt array into xarray
-  arr = xt::adapt(buffer, shape);
+  read_attr(obj_id, name, H5TypeMap<T>::type_id, tensor.data());
 }
 
 // overload for std::string
@@ -290,63 +284,34 @@ void read_dataset(
 }
 
 template<typename T>
-void read_dataset(hid_t dset, xt::xarray<T>& arr, bool indep = false)
+void read_dataset(hid_t dset, tensor::Tensor<T>& tensor, bool indep = false)
 {
   // Get shape of dataset
   vector<hsize_t> shape = object_shape(dset);
 
-  // Allocate space in the array to read data into
-  std::size_t size = 1;
-  for (const auto x : shape)
-    size *= x;
-  arr.resize(shape);
+  // Resize tensor and read data directly
+  vector<size_t> tshape(shape.begin(), shape.end());
+  tensor.resize(tshape);
 
-  // Read data from attribute
+  // Read data from dataset
   read_dataset_lowlevel(
-    dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep, arr.data());
+    dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep, tensor.data());
 }
 
 template<>
 void read_dataset(
-  hid_t dset, xt::xarray<std::complex<double>>& arr, bool indep);
+  hid_t dset, tensor::Tensor<std::complex<double>>& tensor, bool indep);
 
 template<typename T>
 void read_dataset(
-  hid_t obj_id, const char* name, xt::xarray<T>& arr, bool indep = false)
+  hid_t obj_id, const char* name, tensor::Tensor<T>& tensor, bool indep = false)
 {
-  // Open dataset and read array
+  // Open dataset and read tensor
   hid_t dset = open_dataset(obj_id, name);
-  read_dataset(dset, arr, indep);
+  read_dataset(dset, tensor, indep);
   close_dataset(dset);
 }
 
-template<typename T, std::size_t N>
-void read_dataset(
-  hid_t obj_id, const char* name, xt::xtensor<T, N>& arr, bool indep = false)
-{
-  // Open dataset and read array
-  hid_t dset = open_dataset(obj_id, name);
-
-  // Get shape of dataset
-  vector<hsize_t> hsize_t_shape = object_shape(dset);
-  close_dataset(dset);
-
-  // cast from hsize_t to size_t
-  vector<size_t> shape(hsize_t_shape.size());
-  for (int i = 0; i < shape.size(); i++) {
-    shape[i] = static_cast<size_t>(hsize_t_shape[i]);
-  }
-
-  // Allocate new xarray to read data into
-  xt::xarray<T> xarr(shape);
-
-  // Read data from the dataset
-  read_dataset(obj_id, name, xarr);
-
-  // Copy into xtensor
-  arr = xarr;
-}
-
 // overload for Position
 inline void read_dataset(
   hid_t obj_id, const char* name, Position& r, bool indep = false)
@@ -358,31 +323,22 @@ inline void read_dataset(
   r.z = x[2];
 }
 
-template<typename T, std::size_t N>
+template<typename T>
 inline void read_dataset_as_shape(
-  hid_t obj_id, const char* name, xt::xtensor<T, N>& arr, bool indep = false)
+  hid_t obj_id, const char* name, tensor::Tensor<T>& tensor, bool indep = false)
 {
   hid_t dset = open_dataset(obj_id, name);
 
-  // Allocate new array to read data into
-  std::size_t size = 1;
-  for (const auto x : arr.shape())
-    size *= x;
-  vector<T> buffer(size);
-
-  // Read data from attribute
+  // Read data directly into pre-shaped tensor
   read_dataset_lowlevel(
-    dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep, buffer.data());
-
-  // Adapt into xarray
-  arr = xt::adapt(buffer, arr.shape());
+    dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep, tensor.data());
 
   close_dataset(dset);
 }
 
-template<typename T, std::size_t N>
-inline void read_nd_vector(hid_t obj_id, const char* name,
-  xt::xtensor<T, N>& result, bool must_have = false)
+template<typename T>
+inline void read_nd_tensor(hid_t obj_id, const char* name,
+  tensor::Tensor<T>& result, bool must_have = false)
 {
   if (object_exists(obj_id, name)) {
     read_dataset_as_shape(obj_id, name, result, true);
@@ -496,12 +452,16 @@ inline void write_dataset(
     false, buffer.data());
 }
 
-// Template for xarray, xtensor, etc.
-template<typename D>
-inline void write_dataset(
-  hid_t obj_id, const char* name, const xt::xcontainer<D>& arr)
+// Template for Tensor and StaticTensor2D. A SFINAE guard is used here to
+// prevent this template from matching vector/string types that have their own
+// overloads above. A generic Container parameter avoids duplicating the body
+// for both Tensor<T> and StaticTensor2D<T,R,C>.
+template<typename Container,
+  typename =
+    std::enable_if_t<tensor::is_tensor<std::decay_t<Container>>::value>>
+inline void write_dataset(hid_t obj_id, const char* name, const Container& arr)
 {
-  using T = typename D::value_type;
+  using T = typename std::decay_t<Container>::value_type;
   auto s = arr.shape();
   vector<hsize_t> dims {s.cbegin(), s.cend()};
   write_dataset_lowlevel(obj_id, dims.size(), dims.data(), name,

include/openmc/lattice.h

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

include/openmc/material.h

@@ -5,8 +5,8 @@
 #include <unordered_map>
 
 #include "openmc/span.h"
+#include "openmc/tensor.h"
 #include "pugixml.hpp"
-#include "xtensor/xtensor.hpp"
 #include <hdf5.h>
 
 #include "openmc/bremsstrahlung.h"
@@ -189,7 +189,7 @@ public:
   vector<int> nuclide_;                 //!< Indices in nuclides vector
   vector<int> element_;                 //!< Indices in elements vector
   NCrystalMat ncrystal_mat_;            //!< NCrystal material object
-  xt::xtensor<double, 1> atom_density_; //!< Nuclide atom density in [atom/b-cm]
+  tensor::Tensor<double> atom_density_; //!< Nuclide atom density in [atom/b-cm]
   double density_;                      //!< Total atom density in [atom/b-cm]
   double density_gpcc_;                 //!< Total atom density in [g/cm^3]
   double charge_density_;               //!< Total charge density in [e/b-cm]

include/openmc/math_functions.h

@@ -201,6 +201,18 @@ std::complex<double> faddeeva(std::complex<double> z);
 //! \return Derivative of Faddeeva function evaluated at z
 std::complex<double> w_derivative(std::complex<double> z, int order);
 
+//! Evaluate relative exponential function
+//!
+//! \param x Real argument
+//! \return (exp(x)-1)/x without loss of precision near 0
+double exprel(double x);
+
+//! Evaluate relative logarithm function
+//!
+//! \param x Real argument
+//! \return log(1+x)/x without loss of precision near 0
+double log1prel(double x);
+
 //! Helper function to get index and interpolation function on an incident
 //! energy grid
 //!
@@ -211,5 +223,15 @@ std::complex<double> w_derivative(std::complex<double> z, int order);
 void get_energy_index(
   const vector<double>& energies, double E, int& i, double& f);
 
+//==============================================================================
+//! Calculate the cumulative distribution function of the standard normal
+//! distribution at a given value.
+//!
+//! \param z The value at which to evaluate the CDF
+//! \return Phi(z) = P(X <= z) for X ~ N(0,1)
+//==============================================================================
+
+double standard_normal_cdf(double z);
+
 } // namespace openmc
 #endif // OPENMC_MATH_FUNCTIONS_H

include/openmc/mesh.h

@@ -8,8 +8,8 @@
 #include <unordered_map>
 
 #include "hdf5.h"
+#include "openmc/tensor.h"
 #include "pugixml.hpp"
-#include "xtensor/xtensor.hpp"
 
 #include "openmc/bounding_box.h"
 #include "openmc/error.h"
@@ -284,8 +284,8 @@ public:
   virtual Position upper_right() const = 0;
 
   // Data members
-  xt::xtensor<double, 1> lower_left_;  //!< Lower-left coordinates of mesh
-  xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
+  tensor::Tensor<double> lower_left_;  //!< Lower-left coordinates of mesh
+  tensor::Tensor<double> upper_right_; //!< Upper-right coordinates of mesh
   int id_ {-1};                        //!< Mesh ID
   std::string name_;                   //!< User-specified name
   int n_dimension_ {-1};               //!< Number of dimensions
@@ -348,7 +348,7 @@ public:
   //! \param[in] Pointer to bank sites
   //! \param[in] Number of bank sites
   //! \param[out] Whether any bank sites are outside the mesh
-  xt::xtensor<double, 1> count_sites(
+  tensor::Tensor<double> count_sites(
     const SourceSite* bank, int64_t length, bool* outside) const;
 
   //! Get bin given mesh indices
@@ -419,8 +419,8 @@ public:
   //! Get a label for the mesh bin
   std::string bin_label(int bin) const override;
 
-  //! Get shape as xt::xtensor
-  xt::xtensor<int, 1> get_x_shape() const;
+  //! Get mesh dimensions as a tensor
+  tensor::Tensor<int> get_shape_tensor() const;
 
   double volume(int bin) const override
   {
@@ -515,7 +515,7 @@ public:
   //! \param[in] bank Array of bank sites
   //! \param[out] Whether any bank sites are outside the mesh
   //! \return Array indicating number of sites in each mesh/energy bin
-  xt::xtensor<double, 1> count_sites(
+  tensor::Tensor<double> count_sites(
     const SourceSite* bank, int64_t length, bool* outside) const;
 
   //! Return the volume for a given mesh index
@@ -526,7 +526,7 @@ public:
   // Data members
   double volume_frac_;           //!< Volume fraction of each mesh element
   double element_volume_;        //!< Volume of each mesh element
-  xt::xtensor<double, 1> width_; //!< Width of each mesh element
+  tensor::Tensor<double> width_; //!< Width of each mesh element
 };
 
 class RectilinearMesh : public StructuredMesh {

include/openmc/mgxs.h

@@ -6,7 +6,7 @@
 
 #include <string>
 
-#include "xtensor/xtensor.hpp"
+#include "openmc/tensor.h"
 
 #include "openmc/constants.h"
 #include "openmc/hdf5_interface.h"
@@ -22,7 +22,7 @@ namespace openmc {
 
 class Mgxs {
 private:
-  xt::xtensor<double, 1> kTs; // temperature in eV (k * T)
+  tensor::Tensor<double> kTs; // temperature in eV (k * T)
   AngleDistributionType
     scatter_format;       // flag for if this is legendre, histogram, or tabular
   int num_groups;         // number of energy groups
@@ -113,6 +113,11 @@ public:
     const vector<Mgxs*>& micros, const vector<double>& atom_densities,
     int num_group, int num_delay);
 
+  //! \brief Get the number of temperature data points.
+  //!
+  //! @return The number of temperature data points for this MGXS
+  inline int n_temperature_points() { return kTs.size(); }
+
   //! \brief Provides a cross section value given certain parameters
   //!
   //! @param xstype Type of cross section requested, according to the

include/openmc/mgxs_interface.h

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

include/openmc/nuclide.h

@@ -96,7 +96,7 @@ public:
   // Temperature dependent cross section data
   vector<double> kTs_;                //!< temperatures in eV (k*T)
   vector<EnergyGrid> grid_;           //!< Energy grid at each temperature
-  vector<xt::xtensor<double, 2>> xs_; //!< Cross sections at each temperature
+  vector<tensor::Tensor<double>> xs_; //!< Cross sections at each temperature
 
   // Multipole data
   unique_ptr<WindowedMultipole> multipole_;

include/openmc/output.h

@@ -10,6 +10,8 @@
 
 namespace openmc {
 
+extern "C" const bool STRICT_FP_ENABLED;
+
 //! \brief Display the main title banner as well as information about the
 //! program developers, version, and date/time which the problem was run.
 void title();
@@ -60,7 +62,6 @@ void write_tallies();
 void show_time(const char* label, double secs, int indent_level = 0);
 
 } // namespace openmc
-#endif // OPENMC_OUTPUT_H
 
 //////////////////////////////////////
 // Custom formatters
@@ -87,3 +88,5 @@ struct formatter<std::array<T, 2>> {
 }; // namespace fmt
 
 } // namespace fmt
+
+#endif // OPENMC_OUTPUT_H

include/openmc/particle.h

@@ -39,6 +39,8 @@ public:
 
   double speed() const;
 
+  double mass() const;
+
   //! create a secondary particle
   //
   //! stores the current phase space attributes of the particle in the
@@ -69,7 +71,8 @@ public:
   void event_advance();
   void event_cross_surface();
   void event_collide();
-  void event_revive_from_secondary();
+  void event_revive_from_secondary(const SourceSite& site);
+  void event_check_limit_and_revive();
   void event_death();
 
   //! pulse-height recording

include/openmc/particle_data.h

@@ -50,8 +50,11 @@ struct SourceSite {
 
   // Extra attributes that don't show up in source written to file
   int parent_nuclide {-1};
-  int64_t parent_id;
-  int64_t progeny_id;
+  int64_t parent_id {0};
+  int64_t progeny_id {0};
+  double wgt_born {1.0};
+  double wgt_ww_born {-1.0};
+  int64_t n_split {0};
 };
 
 struct CollisionTrackSite {
@@ -533,9 +536,14 @@ private:
   uint64_t seeds_[N_STREAMS];
   int stream_;
 
-  vector<SourceSite> secondary_bank_;
+  vector<SourceSite> local_secondary_bank_;
 
-  int64_t current_work_;
+  // Keep track of how many secondary particles were created in the collision
+  // and what the starting index is in the secondary bank for this particle
+  int n_secondaries_ {0};
+  int secondary_bank_index_ {0};
+
+  int64_t current_work_ {0};
 
   vector<double> flux_derivs_;
 
@@ -558,7 +566,9 @@ private:
 
   int n_event_ {0};
 
-  int n_split_ {0};
+  int64_t n_tracks_ {0}; //!< number of tracks in this particle history
+
+  int64_t n_split_ {0};
   double ww_factor_ {0.0};
 
   int64_t n_progeny_ {0};
@@ -688,8 +698,23 @@ public:
   int& stream() { return stream_; }
 
   // secondary particle bank
-  SourceSite& secondary_bank(int i) { return secondary_bank_[i]; }
-  decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; }
+  SourceSite& local_secondary_bank(int i) { return local_secondary_bank_[i]; }
+  const SourceSite& local_secondary_bank(int i) const
+  {
+    return local_secondary_bank_[i];
+  }
+  decltype(local_secondary_bank_)& local_secondary_bank()
+  {
+    return local_secondary_bank_;
+  }
+
+  // Number of secondaries created in a collision
+  int& n_secondaries() { return n_secondaries_; }
+  const int& n_secondaries() const { return n_secondaries_; }
+
+  // Starting index in secondary bank for this collision
+  int& secondary_bank_index() { return secondary_bank_index_; }
+  const int& secondary_bank_index() const { return secondary_bank_index_; }
 
   // Current simulation work index
   int64_t& current_work() { return current_work_; }
@@ -729,13 +754,16 @@ public:
   int& n_event() { return n_event_; }
 
   // Number of times variance reduction has caused a particle split
-  int n_split() const { return n_split_; }
-  int& n_split() { return n_split_; }
+  int64_t n_split() const { return n_split_; }
+  int64_t& n_split() { return n_split_; }
 
   // Particle-specific factor for on-the-fly weight window adjustment
   double ww_factor() const { return ww_factor_; }
   double& ww_factor() { return ww_factor_; }
 
+  // Number of tracks in this particle history
+  int64_t& n_tracks() { return n_tracks_; }
+
   // Number of progeny produced by this particle
   int64_t& n_progeny() { return n_progeny_; }
 

include/openmc/particle_type.h

@@ -12,6 +12,7 @@
 #include <type_traits>
 
 #include "openmc/constants.h"
+#include "openmc/error.h"
 
 namespace openmc {
 
@@ -61,6 +62,17 @@ public:
   //----------------------------------------------------------------------------
   // Methods
 
+  // Get particle mass in [u]
+  double mass() const
+  {
+    int32_t p = std::abs(pdg_number_);
+    if (ATOMIC_MASS.count(p)) {
+      return ATOMIC_MASS[p];
+    } else {
+      fatal_error("Unknown mass for particle " + str());
+    }
+  }
+
   // Convert to string representation
   std::string str() const;
 

include/openmc/photon.h

@@ -6,7 +6,7 @@
 #include "openmc/particle.h"
 #include "openmc/vector.h"
 
-#include "xtensor/xtensor.hpp"
+#include "openmc/tensor.h"
 #include <hdf5.h>
 
 #include <string>
@@ -62,14 +62,14 @@ public:
   int64_t index_;    //!< Index in global elements vector
 
   // Microscopic cross sections
-  xt::xtensor<double, 1> energy_;
-  xt::xtensor<double, 1> coherent_;
-  xt::xtensor<double, 1> incoherent_;
-  xt::xtensor<double, 1> photoelectric_total_;
-  xt::xtensor<double, 1> pair_production_total_;
-  xt::xtensor<double, 1> pair_production_electron_;
-  xt::xtensor<double, 1> pair_production_nuclear_;
-  xt::xtensor<double, 1> heating_;
+  tensor::Tensor<double> energy_;
+  tensor::Tensor<double> coherent_;
+  tensor::Tensor<double> incoherent_;
+  tensor::Tensor<double> photoelectric_total_;
+  tensor::Tensor<double> pair_production_total_;
+  tensor::Tensor<double> pair_production_electron_;
+  tensor::Tensor<double> pair_production_nuclear_;
+  tensor::Tensor<double> heating_;
 
   // Form factors
   Tabulated1D incoherent_form_factor_;
@@ -81,27 +81,27 @@ public:
   // stored separately to improve memory access pattern when calculating the
   // total cross section
   vector<ElectronSubshell> shells_;
-  xt::xtensor<double, 2> cross_sections_;
+  tensor::Tensor<double> cross_sections_;
 
   // Compton profile data
-  xt::xtensor<double, 2> profile_pdf_;
-  xt::xtensor<double, 2> profile_cdf_;
-  xt::xtensor<double, 1> binding_energy_;
-  xt::xtensor<double, 1> electron_pdf_;
+  tensor::Tensor<double> profile_pdf_;
+  tensor::Tensor<double> profile_cdf_;
+  tensor::Tensor<double> binding_energy_;
+  tensor::Tensor<double> electron_pdf_;
 
   // Map subshells from Compton profile data obtained from Biggs et al,
   // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic
   // relaxation data
-  xt::xtensor<int, 1> subshell_map_;
+  tensor::Tensor<int> subshell_map_;
 
   // Stopping power data
   double I_; // mean excitation energy
-  xt::xtensor<int, 1> n_electrons_;
-  xt::xtensor<double, 1> ionization_energy_;
-  xt::xtensor<double, 1> stopping_power_radiative_;
+  tensor::Tensor<int> n_electrons_;
+  tensor::Tensor<double> ionization_energy_;
+  tensor::Tensor<double> stopping_power_radiative_;
 
   // Bremsstrahlung scaled DCS
-  xt::xtensor<double, 2> dcs_;
+  tensor::Tensor<double> dcs_;
 
   // Whether atomic relaxation data is present
   bool has_atomic_relaxation_ {false};
@@ -137,7 +137,7 @@ void free_memory_photon();
 
 namespace data {
 
-extern xt::xtensor<double, 1>
+extern tensor::Tensor<double>
   compton_profile_pz; //! Compton profile momentum grid
 
 //! Photon interaction data for each element

include/openmc/physics_common.h

@@ -13,5 +13,9 @@ namespace openmc {
 //! \param[in] weight_survive Weight assigned to particles that survive
 void russian_roulette(Particle& p, double weight_survive);
 
+//! \brief Performs the global russian roulette operation on a particle
+//! \param[in,out] p  Particle object
+void apply_russian_roulette(Particle& p);
+
 } // namespace openmc
 #endif // OPENMC_PHYSICS_COMMON_H

include/openmc/plot.h

@@ -6,8 +6,8 @@
 #include <unordered_map>
 #include <unordered_set>
 
+#include "openmc/tensor.h"
 #include "pugixml.hpp"
-#include "xtensor/xarray.hpp"
 
 #include "hdf5.h"
 #include "openmc/cell.h"
@@ -17,6 +17,7 @@
 #include "openmc/particle.h"
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
+#include "openmc/ray.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -83,19 +84,27 @@ const RGBColor BLACK {0, 0, 0};
  * \class PlottableInterface
  * \brief Interface for plottable objects.
  *
- * PlottableInterface classes must have a unique ID in the plots.xml file.
- * They guarantee the ability to create output in some form. This interface
- * is designed to be implemented by classes that produce plot-relevant data
- * which can be visualized.
+ * PlottableInterface classes must have unique IDs. If no ID (or -1) is
+ * provided, the next available ID is assigned automatically. They guarantee
+ * the ability to create output in some form. This interface is designed to be
+ * implemented by classes that produce plot-relevant data which can be
+ * visualized.
  */
+
+typedef tensor::Tensor<RGBColor> ImageData;
 class PlottableInterface {
+public:
+  PlottableInterface() = default;
+
+  void set_default_colors();
+
 private:
   void set_id(pugi::xml_node plot_node);
-  int id_; // unique plot ID
+  int id_ {C_NONE}; // unique plot ID
 
   void set_bg_color(pugi::xml_node plot_node);
   void set_universe(pugi::xml_node plot_node);
-  void set_default_colors(pugi::xml_node plot_node);
+  void set_color_by(pugi::xml_node plot_node);
   void set_user_colors(pugi::xml_node plot_node);
   void set_overlap_color(pugi::xml_node plot_node);
   void set_mask(pugi::xml_node plot_node);
@@ -107,30 +116,36 @@ protected:
 public:
   enum class PlotColorBy { cells = 0, mats = 1 };
 
+  // Generates image data based on plot parameters and returns it
+  virtual ImageData create_image() const = 0;
+
   // Creates the output image named path_plot_
   virtual void create_output() const = 0;
 
+  // Write populated image data to file
+  void write_image(const ImageData& data) const;
+
   // Print useful info to the terminal
   virtual void print_info() const = 0;
 
   const std::string& path_plot() const { return path_plot_; }
   std::string& path_plot() { return path_plot_; }
   int id() const { return id_; }
+  void set_id(int id = C_NONE);
   int level() const { return level_; }
+  PlotColorBy color_by() const { return color_by_; }
 
   // Public color-related data
   PlottableInterface(pugi::xml_node plot_node);
   virtual ~PlottableInterface() = default;
-  int level_;                    // Universe level to plot
-  bool color_overlaps_;          // Show overlapping cells?
-  PlotColorBy color_by_;         // Plot coloring (cell/material)
-  RGBColor not_found_ {WHITE};   // Plot background color
-  RGBColor overlap_color_ {RED}; // Plot overlap color
-  vector<RGBColor> colors_;      // Plot colors
+  int level_ {-1};                           // Universe level to plot
+  bool color_overlaps_ {false};              // Show overlapping cells?
+  PlotColorBy color_by_ {PlotColorBy::mats}; // Plot coloring (cell/material)
+  RGBColor not_found_ {WHITE};               // Plot background color
+  RGBColor overlap_color_ {RED};             // Plot overlap color
+  vector<RGBColor> colors_;                  // Plot colors
 };
 
-typedef xt::xtensor<RGBColor, 2> ImageData;
-
 struct IdData {
   // Constructor
   IdData(size_t h_res, size_t v_res);
@@ -140,7 +155,7 @@ struct IdData {
   void set_overlap(size_t y, size_t x);
 
   // Members
-  xt::xtensor<int32_t, 3> data_; //!< 2D array of cell & material ids
+  tensor::Tensor<int32_t> data_; //!< 2D array of cell & material ids
 };
 
 struct PropertyData {
@@ -152,7 +167,7 @@ struct PropertyData {
   void set_overlap(size_t y, size_t x);
 
   // Members
-  xt::xtensor<double, 3> data_; //!< 2D array of temperature & density data
+  tensor::Tensor<double> data_; //!< 2D array of temperature & density data
 };
 
 //===============================================================================
@@ -166,6 +181,11 @@ public:
 
   enum class PlotBasis { xy = 1, xz = 2, yz = 3 };
 
+  // Accessors
+
+  const std::array<size_t, 3>& pixels() const { return pixels_; }
+  std::array<size_t, 3>& pixels() { return pixels_; }
+
   // Members
 public:
   Position origin_;           //!< Plot origin in geometry
@@ -270,11 +290,11 @@ private:
 public:
   // Add mesh lines to ImageData
   void draw_mesh_lines(ImageData& data) const;
-  void create_image() const;
+  ImageData create_image() const override;
   void create_voxel() const;
 
-  virtual void create_output() const;
-  virtual void print_info() const;
+  void create_output() const override;
+  void print_info() const override;
 
   PlotType type_;                 //!< Plot type (Slice/Voxel)
   int meshlines_width_;           //!< Width of lines added to the plot
@@ -294,17 +314,32 @@ public:
  */
 class RayTracePlot : public PlottableInterface {
 public:
+  RayTracePlot() = default;
   RayTracePlot(pugi::xml_node plot);
 
   // Standard getters. No setting since it's done from XML.
   const Position& camera_position() const { return camera_position_; }
+  Position& camera_position() { return camera_position_; }
   const Position& look_at() const { return look_at_; }
+  Position& look_at() { return look_at_; }
+
   const double& horizontal_field_of_view() const
   {
     return horizontal_field_of_view_;
   }
+  double& horizontal_field_of_view() { return horizontal_field_of_view_; }
 
-  virtual void print_info() const;
+  void print_info() const override;
+
+  const std::array<int, 2>& pixels() const { return pixels_; }
+  std::array<int, 2>& pixels() { return pixels_; }
+
+  const Direction& up() const { return up_; }
+  Direction& up() { return up_; }
+
+  //! brief Updates the cached camera-to-model matrix after changes to
+  //! camera parameters.
+  void update_view();
 
 protected:
   Direction camera_x_axis() const
@@ -330,8 +365,6 @@ protected:
    */
   std::pair<Position, Direction> get_pixel_ray(int horiz, int vert) const;
 
-  std::array<int, 2> pixels_; // pixel dimension of resulting image
-
 private:
   void set_look_at(pugi::xml_node node);
   void set_camera_position(pugi::xml_node node);
@@ -341,9 +374,9 @@ private:
 
   double horizontal_field_of_view_ {70.0}; // horiz. f.o.v. in degrees
   Position camera_position_;               // where camera is
-  Position look_at_; // point camera is centered looking at
-
-  Direction up_ {0.0, 0.0, 1.0}; // which way is up
+  Position look_at_;                     // point camera is centered looking at
+  std::array<int, 2> pixels_ {100, 100}; // pixel dimension of resulting image
+  Direction up_ {0.0, 0.0, 1.0};         // which way is up
 
   /* The horizontal thickness, if using an orthographic projection.
    * If set to zero, we assume using a perspective projection.
@@ -377,8 +410,9 @@ class WireframeRayTracePlot : public RayTracePlot {
 public:
   WireframeRayTracePlot(pugi::xml_node plot);
 
-  virtual void create_output() const;
-  virtual void print_info() const;
+  ImageData create_image() const override;
+  void create_output() const override;
+  void print_info() const override;
 
 private:
   void set_opacities(pugi::xml_node node);
@@ -434,10 +468,22 @@ class SolidRayTracePlot : public RayTracePlot {
   friend class PhongRay;
 
 public:
+  SolidRayTracePlot() = default;
+
   SolidRayTracePlot(pugi::xml_node plot);
 
-  virtual void create_output() const;
-  virtual void print_info() const;
+  ImageData create_image() const override;
+  void create_output() const override;
+  void print_info() const override;
+
+  const std::unordered_set<int>& opaque_ids() const { return opaque_ids_; }
+  std::unordered_set<int>& opaque_ids() { return opaque_ids_; }
+
+  const Position& light_location() const { return light_location_; }
+  Position& light_location() { return light_location_; }
+
+  const double& diffuse_fraction() const { return diffuse_fraction_; }
+  double& diffuse_fraction() { return diffuse_fraction_; }
 
 private:
   void set_opaque_ids(pugi::xml_node node);
@@ -452,43 +498,6 @@ private:
   Position light_location_;
 };
 
-// Base class that implements ray tracing logic, not necessarily through
-// defined regions of the geometry but also outside of it.
-class Ray : public GeometryState {
-
-public:
-  Ray(Position r, Direction u) { init_from_r_u(r, u); }
-
-  // Called at every surface intersection within the model
-  virtual void on_intersection() = 0;
-
-  /*
-   * Traces the ray through the geometry, calling on_intersection
-   * at every surface boundary.
-   */
-  void trace();
-
-  // Stops the ray and exits tracing when called from on_intersection
-  void stop() { stop_ = true; }
-
-  // Sets the dist_ variable
-  void compute_distance();
-
-protected:
-  // Records how far the ray has traveled
-  double traversal_distance_ {0.0};
-
-private:
-  // Max intersections before we assume ray tracing is caught in an infinite
-  // loop:
-  static const int MAX_INTERSECTIONS = 1000000;
-
-  bool hit_something_ {false};
-  bool stop_ {false};
-
-  unsigned event_counter_ {0};
-};
-
 class ProjectionRay : public Ray {
 public:
   ProjectionRay(Position r, Direction u, const WireframeRayTracePlot& plot,
@@ -496,7 +505,7 @@ public:
     : Ray(r, u), plot_(plot), line_segments_(line_segments)
   {}
 
-  virtual void on_intersection() override;
+  void on_intersection() override;
 
 private:
   /* Store a reference to the plot object which is running this ray, in order
@@ -519,7 +528,7 @@ public:
     result_color_ = plot_.not_found_;
   }
 
-  virtual void on_intersection() override;
+  void on_intersection() override;
 
   const RGBColor& result_color() { return result_color_; }
 

include/openmc/random_ray/flat_source_domain.h

@@ -40,9 +40,10 @@ public:
   void random_ray_tally();
   virtual void accumulate_iteration_flux();
   void output_to_vtk() const;
-  void convert_external_sources();
+  void convert_external_sources(bool use_adjoint_sources);
   void count_external_source_regions();
-  void set_adjoint_sources();
+  void set_fw_adjoint_sources();
+  void set_local_adjoint_sources();
   void flux_swap();
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
@@ -76,6 +77,7 @@ public:
   // Static Data members
   static bool volume_normalized_flux_tallies_;
   static bool adjoint_; // If the user wants outputs based on the adjoint flux
+  static bool fw_cadis_local_;
   static double
     diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization
                                  // for transport corrected MGXS data
@@ -84,6 +86,8 @@ public:
   static std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
     mesh_domain_map_;
 
+  static std::vector<size_t> fw_cadis_local_targets_;
+
   //----------------------------------------------------------------------------
   // Static data members
   static RandomRayVolumeEstimator volume_estimator_;
@@ -100,17 +104,18 @@ public:
   // in model::cells
   vector<int64_t> source_region_offsets_;
 
-  // 2D arrays stored in 1D representing values for all materials x energy
-  // groups
+  // 3D arrays stored in 1D representing values for all materials x temperature
+  // points x energy groups
   int n_materials_;
+  int ntemperature_;
   vector<double> sigma_t_;
   vector<double> nu_sigma_f_;
   vector<double> sigma_f_;
   vector<double> chi_;
   vector<double> kappa_fission_;
 
-  // 3D arrays stored in 1D representing values for all materials x energy
-  // groups x energy groups
+  // 4D arrays stored in 1D representing values for all materials x temperature
+  // points x energy groups x energy groups
   vector<double> sigma_s_;
 
   // The abstract container holding all source region-specific data
@@ -177,10 +182,10 @@ protected:
   // Volumes for each tally and bin/score combination. This intermediate data
   // structure is used when tallying quantities that must be normalized by
   // volume (i.e., flux). The vector is index by tally index, while the inner 2D
-  // xtensor is indexed by bin index and score index in a similar manner to the
+  // tensor is indexed by bin index and score index in a similar manner to the
   // results tensor in the Tally class, though without the third dimension, as
   // SUM and SUM_SQ do not need to be tracked.
-  vector<xt::xtensor<double, 2>> tally_volumes_;
+  vector<tensor::Tensor<double>> tally_volumes_;
 
 }; // class FlatSourceDomain
 

include/openmc/random_ray/random_ray.h

@@ -41,6 +41,7 @@ public:
   uint64_t transport_history_based_single_ray();
   SourceSite sample_prng();
   SourceSite sample_halton();
+  SourceSite sample_s2();
 
   //----------------------------------------------------------------------------
   // Static data members
@@ -65,6 +66,7 @@ private:
   vector<double> mesh_fractional_lengths_;
 
   int negroups_;
+  int ntemperature_;
   FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source
                                        // data needed for ray transport
   double distance_travelled_ {0};

include/openmc/random_ray/random_ray_simulation.h

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

include/openmc/random_ray/source_region.h

@@ -146,6 +146,7 @@ public:
 
   // Scalar fields
   int* material_;
+  int* temperature_idx_;
   double* density_mult_;
   int* is_small_;
   int* n_hits_;
@@ -199,6 +200,9 @@ public:
   double& density_mult() { return *density_mult_; }
   const double density_mult() const { return *density_mult_; }
 
+  int& temperature_idx() { return *temperature_idx_; }
+  const int temperature_idx() const { return *temperature_idx_; }
+
   int& is_small() { return *is_small_; }
   const int is_small() const { return *is_small_; }
 
@@ -319,8 +323,9 @@ public:
 
   //---------------------------------------
   // Scalar fields
-
-  int material_ {0};          //!< Index in openmc::model::materials array
+  int material_ {0}; //!< Index in openmc::model::materials array
+  int temperature_idx_ {
+    0}; //!< Index into the MGXS array representing temperature
   double density_mult_ {1.0}; //!< A density multiplier queried from the cell
                               //!< corresponding to the source region.
   OpenMPMutex lock_;
@@ -400,6 +405,9 @@ public:
   int& material(int64_t sr) { return material_[sr]; }
   const int material(int64_t sr) const { return material_[sr]; }
 
+  int& temperature_idx(int64_t sr) { return temperature_idx_[sr]; }
+  const int temperature_idx(int64_t sr) const { return temperature_idx_[sr]; }
+
   double& density_mult(int64_t sr) { return density_mult_[sr]; }
   const double density_mult(int64_t sr) const { return density_mult_[sr]; }
 
@@ -634,6 +642,7 @@ private:
 
   // SoA storage for scalar fields (one item per source region)
   vector<int> material_;
+  vector<int> temperature_idx_;
   vector<double> density_mult_;
   vector<int> is_small_;
   vector<int> n_hits_;

include/openmc/ray.h

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

include/openmc/reaction.h

@@ -76,11 +76,21 @@ public:
 //! \return Name of the corresponding reaction
 std::string reaction_name(int mt);
 
-//! Return reaction type (MT value) given a reaction name
+//! Return MT value for given a reaction name (including special tally MT
+//! values)
 //
 //! \param[in] name  Reaction name
-//! \return Corresponding reaction type (MT value)
-int reaction_type(std::string name);
+//! \return Corresponding MT number or special tally score
+int reaction_tally_mt(std::string name);
+
+//! Return ENDF MT number given a reaction name
+//
+//! Unlike reaction_tally_mt(), this function always returns a positive ENDF MT
+//! number and never a special negative tally score.
+//!
+//! \param[in] name  Reaction name
+//! \return Corresponding ENDF MT number
+int reaction_mt(const std::string& name);
 
 } // namespace openmc
 

include/openmc/reaction_product.h

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