Merge remote-tracking branch 'upstream/develop' into covariance_storage

.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/workflows/ci.yml

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

.github/workflows/dockerhub-publish.yml

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

.github/workflows/format-check.yml

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

.gitignore

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

.mcp.json

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

.readthedocs.yaml

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

AGENTS.md

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

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.

Dockerfile

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

docs/source/devguide/agentic-tools.rst

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

docs/source/devguide/index.rst

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

docs/source/io_formats/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/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

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

openmc/deplete/abc.py

@@ -201,7 +201,7 @@ class TransportOperator(ABC):
         Returns
         -------
         volume : dict of str to float
-            Volumes corresponding to materials in burn_list
+            Volumes corresponding to materials in full_burn_list
         nuc_list : list of str
             A list of all nuclide names. Used for sorting the simulation.
         burn_list : list of int
@@ -210,7 +210,7 @@ class TransportOperator(ABC):
         full_burn_list : list of int
             All burnable materials in the geometry.
         name_list : list of str
-            Material names corresponding to materials in burn_list
+            Material names corresponding to materials in full_burn_list
         """
 
     def finalize(self):

openmc/deplete/chain.py

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

openmc/deplete/results.py

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

openmc/deplete/stepresult.py

@@ -154,14 +154,14 @@ class StepResult:
         full_burn_list : list of str
             List of all burnable material IDs
         name_list : list of str, optional
-            Material names corresponding to materials in burn_list
+            Material names corresponding to materials in full_burn_list
 
         """
         self.volume = copy.deepcopy(volume)
         self.index_nuc = {nuc: i for i, nuc in enumerate(nuc_list)}
         self.index_mat = {mat: i for i, mat in enumerate(burn_list)}
         self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
-        self.mat_to_name = dict(zip(burn_list, name_list)) if name_list is not None else {}
+        self.mat_to_name = dict(zip(full_burn_list, name_list)) if name_list is not None else {}
 
         # Create storage array
         self.data = np.zeros((self.n_mat, self.n_nuc))

openmc/material.py

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

openmc/mesh.py

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

openmc/mgxs/groups.py

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

src/mgxs.cpp

@@ -510,6 +510,8 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
     break;
   case MgxsType::INVERSE_VELOCITY:
     val = xs_t->inverse_velocity(a, gin);
+    if (!(val > 0))
+      val = data::mg.default_inverse_velocity_[gin];
     break;
   case MgxsType::DECAY_RATE:
     if (dg != nullptr) {

src/mgxs_interface.cpp

@@ -237,6 +237,19 @@ void MgxsInterface::read_header(const std::string& path_cross_sections)
                 "library file!");
   }
 
+  // Calculate approximate default inverse velocity data
+  for (int i = 0; i < energy_bins_.size() - 1; ++i) {
+    double e_min = std::max(energy_bins_[i + 1], 1e-5);
+    double e_max = energy_bins_[i];
+    double alpha = 1.0 / (C_LIGHT * std::log(e_max / e_min));
+    double k_max = std::sqrt(1 + 2.0 * MASS_NEUTRON_EV / e_max);
+    double k_min = std::sqrt(1 + 2.0 * MASS_NEUTRON_EV / e_min);
+    double inv_v =
+      alpha * (2.0 * (std::atanh(1.0 / k_max) - std::atanh(1.0 / k_min)) -
+                (k_max - k_min));
+    default_inverse_velocity_.push_back(inv_v);
+  }
+
   // Close MGXS HDF5 file
   file_close(file_id);
 }

src/particle.cpp

@@ -48,26 +48,33 @@ double Particle::speed() const
 {
   if (settings::run_CE) {
     // Determine mass in eV/c^2
-    double mass;
-    switch (type().pdg_number()) {
-    case PDG_NEUTRON:
-      mass = MASS_NEUTRON_EV;
-    case PDG_ELECTRON:
-    case PDG_POSITRON:
-      mass = MASS_ELECTRON_EV;
-    default:
-      mass = this->type().mass() * AMU_EV;
-    }
+    double mass = this->mass();
 
     // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
     return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
            (this->E() + mass);
   } else {
-    auto& macro_xs = data::mg.macro_xs_[this->material()];
+    auto mat = this->material();
+    if (mat == MATERIAL_VOID)
+      return 1.0 / data::mg.default_inverse_velocity_[this->g()];
+    auto& macro_xs = data::mg.macro_xs_[mat];
     int macro_t = this->mg_xs_cache().t;
     int macro_a = macro_xs.get_angle_index(this->u());
-    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
-                   nullptr, nullptr, macro_t, macro_a);
+    return 1.0 / macro_xs.get_xs(
+                   MgxsType::INVERSE_VELOCITY, this->g(), macro_t, macro_a);
+  }
+}
+
+double Particle::mass() const
+{
+  switch (type().pdg_number()) {
+  case PDG_NEUTRON:
+    return MASS_NEUTRON_EV;
+  case PDG_ELECTRON:
+  case PDG_POSITRON:
+    return MASS_ELECTRON_EV;
+  default:
+    return this->type().mass() * AMU_EV;
   }
 }
 

src/tallies/tally_scoring.cpp

@@ -2436,24 +2436,22 @@ void score_tracklength_tally_general(
           if (p.material() != MATERIAL_VOID) {
             const auto& mat = model::materials[p.material()];
             auto j = mat->mat_nuclide_index_[i_nuclide];
-            if (j == C_NONE) {
-              // Determine log union grid index
-              if (i_log_union == C_NONE) {
-                int neutron = ParticleType::neutron().transport_index();
-                i_log_union = std::log(p.E() / data::energy_min[neutron]) /
-                              simulation::log_spacing;
-              }
-
-              // Update micro xs cache
-              if (!tally.multiply_density()) {
-                p.update_neutron_xs(i_nuclide, i_log_union);
-                atom_density = 1.0;
-              }
-            } else {
-              atom_density = tally.multiply_density()
-                               ? mat->atom_density(j, p.density_mult())
-                               : 1.0;
+            if (j != C_NONE)
+              atom_density = mat->atom_density(j, p.density_mult());
+          }
+          if (atom_density > 0) {
+            if (!tally.multiply_density())
+              atom_density = 1.0;
+          } else if (!tally.multiply_density()) {
+            // Determine log union grid index
+            if (i_log_union == C_NONE) {
+              int neutron = ParticleType::neutron().transport_index();
+              i_log_union = std::log(p.E() / data::energy_min[neutron]) /
+                            simulation::log_spacing;
             }
+            // Update micro xs cache
+            p.update_neutron_xs(i_nuclide, i_log_union);
+            atom_density = 1.0;
           }
         }
 
@@ -2565,25 +2563,25 @@ void score_collision_tally(Particle& p)
 
         double atom_density = 0.;
         if (i_nuclide >= 0) {
-          const auto& mat = model::materials[p.material()];
-          auto j = mat->mat_nuclide_index_[i_nuclide];
-          if (j == C_NONE) {
+          if (p.material() != MATERIAL_VOID) {
+            const auto& mat = model::materials[p.material()];
+            auto j = mat->mat_nuclide_index_[i_nuclide];
+            if (j != C_NONE)
+              atom_density = mat->atom_density(j, p.density_mult());
+          }
+          if (atom_density > 0) {
+            if (!tally.multiply_density())
+              atom_density = 1.0;
+          } else if (!tally.multiply_density()) {
             // Determine log union grid index
             if (i_log_union == C_NONE) {
               int neutron = ParticleType::neutron().transport_index();
               i_log_union = std::log(p.E() / data::energy_min[neutron]) /
                             simulation::log_spacing;
             }
-
             // Update micro xs cache
-            if (!tally.multiply_density()) {
-              p.update_neutron_xs(i_nuclide, i_log_union);
-              atom_density = 1.0;
-            }
-          } else {
-            atom_density = tally.multiply_density()
-                             ? mat->atom_density(j, p.density_mult())
-                             : 1.0;
+            p.update_neutron_xs(i_nuclide, i_log_union);
+            atom_density = 1.0;
           }
         }
 

tests/unit_tests/test_deplete_chain.py

@@ -246,6 +246,29 @@ def test_form_matrix(simple_chain):
         assert new_mat[r, c] == mat[r, c]
 
 
+def test_decay_matrix(simple_chain):
+    """Test that decay_matrix contains only radioactive decay terms."""
+    # Nuclide order: H1(0), A(1), B(2), C(3)
+    decay_A = log(2) / 2.36520E+04
+    decay_B = log(2) / 3.29040E+04
+
+    expected = np.zeros((4, 4))
+    expected[1, 1] = -decay_A           # Loss: A decays
+    expected[2, 1] = decay_A * 0.6      # A -> B (branching ratio 0.6)
+    expected[3, 1] = decay_A * 0.4      # A -> C (branching ratio 0.4)
+    expected[1, 2] = decay_B            # B -> A (branching ratio 1.0)
+    expected[2, 2] = -decay_B           # Loss: B decays
+
+    assert np.allclose(expected, simple_chain.decay_matrix.toarray())
+
+
+def test_decay_matrix_cached(simple_chain):
+    """Test that decay_matrix is lazily computed and returns the same object."""
+    m1 = simple_chain.decay_matrix
+    m2 = simple_chain.decay_matrix
+    assert m1 is m2
+
+
 def test_getitem():
     """Test nuc_by_ind converter function."""
     chain = Chain()

tests/unit_tests/test_material.py

@@ -594,6 +594,8 @@ def test_get_activity():
     assert pytest.approx(m4.get_activity(units='Bq/g', by_nuclide=True)["H3"]) == 355978108155965.94  # [Bq/g]
     assert pytest.approx(m4.get_activity(units='Bq/cm3')) == 355978108155965.94*3/2 # [Bq/cc]
     assert pytest.approx(m4.get_activity(units='Bq/cm3', by_nuclide=True)["H3"]) == 355978108155965.94*3/2 # [Bq/cc]
+    assert pytest.approx(m4.get_activity(units='Bq/m3')) == 355978108155965.94*3/2*1e6 # [Bq/m3]
+    assert pytest.approx(m4.get_activity(units='Bq/m3', by_nuclide=True)["H3"]) == 355978108155965.94*3/2*1e6 # [Bq/m3]
     # volume is required to calculate total activity
     m4.volume = 10.
     assert pytest.approx(m4.get_activity(units='Bq')) == 355978108155965.94*3/2*10 # [Bq]
@@ -650,6 +652,8 @@ def test_get_decay_heat():
     assert pytest.approx(m4.get_decay_heat(units='W/g', by_nuclide=True)["I135"]) == 40175.15720273193 # [W/g]
     assert pytest.approx(m4.get_decay_heat(units='W/cm3')) == 40175.15720273193*3/2 # [W/cc]
     assert pytest.approx(m4.get_decay_heat(units='W/cm3', by_nuclide=True)["I135"]) == 40175.15720273193*3/2 #[W/cc]
+    assert pytest.approx(m4.get_decay_heat(units='W/m3')) == 40175.15720273193*3/2*1e6 # [W/m3]
+    assert pytest.approx(m4.get_decay_heat(units='W/m3', by_nuclide=True)["I135"]) == 40175.15720273193*3/2*1e6 # [W/m3]
     # volume is required to calculate total decay heat
     m4.volume = 10.
     assert pytest.approx(m4.get_decay_heat(units='W')) == 40175.15720273193*3/2*10 # [W]
@@ -680,6 +684,8 @@ def test_decay_photon_energy():
     src_per_bqg = m.get_decay_photon_energy(units='Bq/g')
     src_per_bqkg = m.get_decay_photon_energy(units='Bq/kg')
     assert pytest.approx(src_per_bqg.integral()) == src_per_bqkg.integral() / 1000.
+    src_per_bqm3 = m.get_decay_photon_energy(units='Bq/m3')
+    assert pytest.approx(src_per_bqm3.integral()) == src_per_cm3.integral() * 1e6
 
     # If we add Xe135 (which has a tabular distribution), the photon source
     # should be a mixture distribution

tests/unit_tests/test_mesh_from_domain.py

@@ -27,6 +27,20 @@ def test_reg_mesh_from_bounding_box():
     assert np.array_equal(mesh.upper_right, bb[1])
 
 
+def test_rectilinear_mesh_from_bounding_box():
+    """Tests a RectilinearMesh can be made from a BoundingBox directly."""
+    bb = openmc.BoundingBox([-8, -7, -5], [12, 13, 15])
+
+    mesh = openmc.RectilinearMesh.from_bounding_box(bb, dimension=[2, 4, 5])
+    assert isinstance(mesh, openmc.RectilinearMesh)
+    assert np.array_equal(mesh.dimension, (2, 4, 5))
+    assert np.array_equal(mesh.lower_left, bb[0])
+    assert np.array_equal(mesh.upper_right, bb[1])
+    assert np.array_equal(mesh.x_grid, [-8., 2., 12.])
+    assert np.array_equal(mesh.y_grid, [-7., -2., 3., 8., 13.])
+    assert np.array_equal(mesh.z_grid, [-5., -1., 3., 7., 11., 15.])
+
+
 def test_cylindrical_mesh_from_cell():
     """Tests a CylindricalMesh can be made from a Cell and the specified
     dimensions are propagated through."""

tests/unit_tests/test_mg_inverse_velocity.py

@@ -0,0 +1,59 @@
+import openmc
+import numpy as np
+import pytest
+
+@pytest.fixture
+def one_group_lib():
+    groups = openmc.mgxs.EnergyGroups([0.0, 20.0e6])
+    xsdata = openmc.XSdata('slab_mat', groups)
+    xsdata.order = 0
+    xsdata.set_total([0.0])               
+    xsdata.set_absorption([0.0])
+    xsdata.set_scatter_matrix([[[0.0]]])
+    
+    mg_library = openmc.MGXSLibrary(groups)
+    mg_library.add_xsdata(xsdata)
+    name = 'mgxs.h5'
+    mg_library.export_to_hdf5(name)
+    yield name
+
+@pytest.fixture
+def slab_model(one_group_lib):
+    model = openmc.Model()
+    mat = openmc.Material(name='slab_material')
+    mat.set_density('macro', 1.0)
+    mat.add_macroscopic('slab_mat')
+
+    model.materials = openmc.Materials([mat])
+    model.materials.cross_sections = one_group_lib
+    
+    x_min = openmc.XPlane(x0=0.0, boundary_type='vacuum')
+    x_max = openmc.XPlane(x0=10.0, boundary_type='vacuum')
+
+    y_min = openmc.YPlane(y0=-10.0, boundary_type='vacuum')
+    y_max = openmc.YPlane(y0=10.0, boundary_type='vacuum')
+    z_min = openmc.ZPlane(z0=-10.0, boundary_type='vacuum')
+    z_max = openmc.ZPlane(z0=19.0, boundary_type='vacuum')
+
+    cell = openmc.Cell(fill=mat, region=+z_min & -x_max & +y_min & -y_max & +z_min & -z_max)
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings = openmc.Settings()
+    model.settings.energy_mode = 'multi-group'
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 3
+    model.settings.particles = 10
+
+    source = openmc.IndependentSource()
+    source.space = openmc.stats.Point((5.0, 0.0, 0.0))
+    model.settings.source = source
+    return model
+
+def test_inverse_velocity(run_in_tmpdir, slab_model):
+    tally = openmc.Tally()
+    tally.scores = ['flux','inverse-velocity']
+    slab_model.tallies = [tally]
+    slab_model.run(apply_tally_results=True)
+    inverse_velocity = tally.mean.squeeze()[1]/tally.mean.squeeze()[0]
+
+    assert inverse_velocity == pytest.approx(1.6144e-5, rel=1e-4)

tests/unit_tests/test_void.py

@@ -0,0 +1,42 @@
+import numpy as np
+import openmc
+import pytest
+
+
+@pytest.fixture
+def empty_sphere():
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+    surf = openmc.Sphere(r=10, boundary_type='vacuum')
+    cell = openmc.Cell(region=-surf)
+    model.geometry = openmc.Geometry([cell])
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 3
+    model.settings.particles = 1000
+
+    tally = openmc.Tally()
+    tally.scores = ['total', 'elastic']
+    tally.nuclides = ['U235']
+    tally.multiply_density = False
+    model.tallies.append(tally)
+
+    return model
+
+
+def test_equivalent_microxs(empty_sphere, run_in_tmpdir):
+    sp_file = empty_sphere.run()
+    with openmc.StatePoint(sp_file) as sp:
+        tally1 = sp.tallies[1]
+
+    mat = openmc.Material()
+    mat.add_nuclide('H1', 1e-16)
+
+    empty_sphere.geometry.get_all_cells()[1].fill = mat
+
+    sp_file = empty_sphere.run()
+    with openmc.StatePoint(sp_file) as sp:
+        tally2 = sp.tallies[1]
+
+    assert np.isclose(tally1.mean.sum(), tally2.mean.sum(), rtol=1e-10, atol=0)
+    assert tally1.mean.sum() > 0

tools/ci/gha-install-njoy.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 set -ex
 cd $HOME
-git clone https://github.com/njoy/NJOY2016
+git clone -b 2016.78 https://github.com/njoy/NJOY2016
 cd NJOY2016
 mkdir build && cd build
 cmake -Dstatic=on .. && make 2>/dev/null && sudo make install