Merge branch 'develop' into particle-ray

This commit is contained in:
GuySten 2026-05-11 10:52:35 +03:00
commit 71da8ed745
165 changed files with 6847 additions and 1137 deletions

View file

@ -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()

View file

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

View file

@ -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()

View file

@ -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()

View file

@ -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()

View file

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

34
.claude/tools/start_server.sh Executable file
View file

@ -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"

View file

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

View file

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

View file

@ -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 }}

View file

@ -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 }}

View file

@ -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 }}

View file

@ -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 }}

View file

@ -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 }}

View file

@ -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 }}

View file

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

View file

@ -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 }}

View file

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

View file

@ -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 }}

View file

@ -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 }}

View file

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

3
.gitignore vendored
View file

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

9
.mcp.json Normal file
View file

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

View file

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

View file

@ -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`

14
CLAUDE.md Normal file
View file

@ -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.

View file

@ -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 \

View file

@ -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.

View file

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

View file

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

View file

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

View file

@ -597,7 +597,7 @@ found in the :ref:`random ray user guide <random_ray>`.
*Default*: None
:source:
:ray_source:
Specifies the starting ray distribution, and follows the format for
:ref:`source_element`. It must be uniform in space and angle and cover the
full domain. It does not represent a physical neutron or photon source -- it
@ -605,6 +605,35 @@ found in the :ref:`random ray user guide <random_ray>`.
*Default*: None
:adjoint_source:
Specifies an adjoint fixed source for adjoint transport simulations, and
follows the format for :ref:`source_element`. The distributions which make
up the adjoint source are subject to the same restrictions as forward
fixed sources in Random Ray mode.
*Default*: None
:adjoint:
Specifies whether to perform adjoint transport. The default is 'False',
corresponding to forward transport.
*Default*: None
:volume_estimator:
Specifies choice of volume estimator for the random ray solver. Options
are 'naive', 'simulation_averaged', or 'hybrid'. The default is 'hybrid'.
*Default*: None
:volume_normalized_flux_tallies:
Specifies whether to normalize flux tallies by volume (bool). The
default is 'False'. When enabled, flux tallies will be reported in units
of cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
region).
*Default*: None
:sample_method:
Specifies the method for sampling the starting ray distribution. This
element can be set to "prng" or "halton".
@ -1029,17 +1058,19 @@ variable and whose sub-elements/attributes are as follows:
:type:
The type of the distribution. Valid options are "uniform", "discrete",
"tabular", "maxwell", "watt", and "mixture". The "uniform" option produces
variates sampled from a uniform distribution over a finite interval. The
"discrete" option produces random variates that can assume a finite number
of values (i.e., a distribution characterized by a probability mass function).
The "tabular" option produces random variates sampled from a tabulated
distribution where the density function is either a histogram or
"tabular", "maxwell", "watt", "mixture", and "decay_spectrum". The "uniform"
option produces variates sampled from a uniform distribution over a finite
interval. The "discrete" option produces random variates that can assume a
finite number of values (i.e., a distribution characterized by a probability
mass function). The "tabular" option produces random variates sampled from a
tabulated distribution where the density function is either a histogram or
linearly-interpolated between tabulated points. The "watt" option produces
random variates is sampled from a Watt fission spectrum (only used for
energies). The "maxwell" option produce variates sampled from a Maxwell
fission spectrum (only used for energies). The "mixture" option produces samples
from univariate sub-distributions with given probabilities.
fission spectrum (only used for energies). The "mixture" option produces
samples from univariate sub-distributions with given probabilities. The
"decay_spectrum" option produces photon energies sampled from decay photon
spectra in a depletion chain (only used for energies).
*Default*: None
@ -1057,6 +1088,10 @@ variable and whose sub-elements/attributes are as follows:
:math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x`
points are given first followed by corresponding :math:`p` points.
For a "decay_spectrum" distribution, ``parameters`` gives the atom densities
in [atom/b-cm] for the nuclides listed in the ``nuclides`` element, in the
same order.
For a "watt" distribution, ``parameters`` should be given as two real numbers
:math:`a` and :math:`b` that parameterize the distribution :math:`p(x) dx = c
e^{-x/a} \sinh \sqrt{b \, x} dx`.
@ -1086,6 +1121,21 @@ variable and whose sub-elements/attributes are as follows:
This sub-element of a ``pair`` element provides information on the
corresponding univariate distribution.
:volume:
For a "decay_spectrum" distribution, this attribute specifies the source
region volume in cm\ :sup:`3`. It is used together with atom densities to
determine the absolute photon emission rate. When a source uses a
"decay_spectrum" energy distribution, the source strength is set from this
emission rate.
:nuclides:
For a "decay_spectrum" distribution, this element specifies a
whitespace-separated list of nuclide names contributing to the decay photon
source. The atom densities for these nuclides are given by the ``parameters``
element in the same order. Nuclides are resolved against the depletion chain,
and nuclides without decay photon spectra do not contribute to the
distribution.
:bias:
This optional element specifies a biased distribution for importance sampling.
For continuous distributions, the ``bias`` element should contain another
@ -1271,7 +1321,7 @@ The ``<surface_grazing_cutoff>`` element specifies the surface flux cosine cutof
``<surface_grazing_ratio>`` Element
-----------------------------------
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
The ``<surface_grazing_ratio>`` element specifies the surface flux cosine
substitution ratio.
*Default*: 0.5
@ -1696,6 +1746,14 @@ mesh-based weight windows.
The ratio of the lower to upper weight window bounds.
*Default*: 5.0
For FW-CADIS:
:targets:
A sequence of IDs corresponding to the tallies which cover phase
space regions of interest for local variance reduction.
*Default*: None
---------------------------------------
``<weight_window_checkpoints>`` Element

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

@ -101,6 +101,8 @@ extern vector<unique_ptr<ChainNuclide>> chain_nuclides;
void read_chain_file_xml();
void free_memory_chain();
} // namespace openmc
#endif // OPENMC_CHAIN_H

View file

@ -407,6 +407,71 @@ private:
double integral_; //!< Integral of distribution
};
//==============================================================================
// DecaySpectrum — non-owning mixture of decay photon distributions
//==============================================================================
//! Energy distribution formed by mixing multiple decay photon spectra.
//!
//! Unlike the general Mixture distribution, this class holds non-owning
//! pointers to the component distributions (which live in
//! data::chain_nuclides). Each component is weighted by the activity
//! (atoms * decay_constant) of the corresponding nuclide.
class DecaySpectrum : public Distribution {
public:
//============================================================================
// Types, aliases
struct Sample {
double energy;
double weight;
int parent_nuclide;
};
//============================================================================
// Constructors
//! Construct from an XML node containing nuclide names and atom densities.
//!
//! Reads child ``<nuclide>`` elements with ``name`` and ``density``
//! attributes, resolves them against the loaded depletion chain, and
//! constructs the mixed distribution.
explicit DecaySpectrum(pugi::xml_node node);
//============================================================================
// Methods
//! Sample a value from the distribution and return the parent nuclide index
//! \param seed Pseudorandom number seed pointer
//! \return (Sampled energy, sample weight, chain nuclide index)
Sample sample_with_parent(uint64_t* seed) const;
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return (sampled value, sample weight)
std::pair<double, double> sample(uint64_t* seed) const override;
double integral() const override;
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
//! Initialize decay spectrum sampling data
//! \param nuclide_indices Indices of decay photon emitters in
//! data::chain_nuclides
//! \param atoms Number of atoms for each component.
void init(vector<int> nuclide_indices, const vector<double>& atoms);
vector<int> nuclide_indices_; //!< Indices of emitting nuclides in the chain
DiscreteIndex di_; //!< Discrete index for component selection
double integral_; //!< Total photon emission rate
};
} // namespace openmc
#endif // OPENMC_DISTRIBUTION_H

View file

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

View file

@ -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 {

View file

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

View file

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

View file

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

View file

@ -43,6 +43,7 @@ class Source;
namespace model {
extern vector<unique_ptr<Source>> external_sources;
extern vector<unique_ptr<Source>> adjoint_sources;
// Probability distribution for selecting external sources
extern DiscreteIndex external_sources_probability;

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

@ -1,5 +1,6 @@
from __future__ import annotations
from collections.abc import Sequence
from contextlib import nullcontext
import copy
from datetime import datetime
import json
@ -8,16 +9,17 @@ from pathlib import Path
import numpy as np
import openmc
from . import IndependentOperator, PredictorIntegrator
from .chain import Chain
from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5
from .results import Results
from ..checkvalue import PathLike
from ..mpi import comm
from openmc.lib import TemporarySession
from openmc.utility_funcs import change_directory
def get_activation_materials(
model: openmc.Model, mmv: openmc.MeshMaterialVolumes
model: openmc.Model,
mmv_list: list[openmc.MeshMaterialVolumes]
) -> openmc.Materials:
"""Get a list of activation materials for each mesh element/material.
@ -31,35 +33,35 @@ def get_activation_materials(
----------
model : openmc.Model
The full model containing the geometry and materials.
mmv : openmc.MeshMaterialVolumes
The mesh material volumes object containing the materials and their
volumes for each mesh element.
mmv_list : list of openmc.MeshMaterialVolumes
List of mesh material volumes objects, one per mesh, containing the
materials and their volumes for each mesh element.
Returns
-------
openmc.Materials
A list of materials, each corresponding to a unique mesh element and
material combination.
material combination across all meshes.
"""
# Get the material ID, volume, and element index for each element-material
# combination
mat_ids = mmv._materials[mmv._materials > -1]
volumes = mmv._volumes[mmv._materials > -1]
elems, _ = np.where(mmv._materials > -1)
# Get all materials in the model
material_dict = model._get_all_materials()
# Create a new activation material for each element-material combination
# across all meshes
materials = openmc.Materials()
for elem, mat_id, vol in zip(elems, mat_ids, volumes):
mat = material_dict[mat_id]
new_mat = mat.clone()
new_mat.depletable = True
new_mat.name = f'Element {elem}, Material {mat_id}'
new_mat.volume = vol
materials.append(new_mat)
for mesh_idx, mmv in enumerate(mmv_list):
mat_ids = mmv._materials[mmv._materials > -1]
volumes = mmv._volumes[mmv._materials > -1]
elems, _ = np.where(mmv._materials > -1)
for elem, mat_id, vol in zip(elems, mat_ids, volumes):
mat = material_dict[mat_id]
new_mat = mat.clone()
new_mat.depletable = True
new_mat.name = f'Mesh {mesh_idx}, Element {elem}, Material {mat_id}'
new_mat.volume = vol
materials.append(new_mat)
return materials
@ -70,7 +72,9 @@ class R2SManager:
This class is responsible for managing the materials and sources needed for
mesh-based or cell-based R2S calculations. It provides methods to get
activation materials and decay photon sources based on the mesh/cells and
materials in the OpenMC model.
materials in the OpenMC model. Multiple meshes can be specified as domains,
in which case each element--material combination of each mesh is treated as
an activation region (meshes are assumed to be non-overlapping).
This class supports the use of a different models for the neutron and photon
transport calculation. However, for cell-based calculations, it assumes that
@ -83,17 +87,20 @@ class R2SManager:
----------
neutron_model : openmc.Model
The OpenMC model to use for neutron transport.
domains : openmc.MeshBase or Sequence[openmc.Cell]
The mesh or a sequence of cells that represent the spatial units over
which the R2S calculation will be performed.
domains : openmc.MeshBase or Sequence[openmc.MeshBase] or Sequence[openmc.Cell]
The mesh(es) or a sequence of cells that represent the spatial units
over which the R2S calculation will be performed. When a single
:class:`~openmc.MeshBase` or a sequence of meshes is given, each
element--material combination across all meshes is treated as an
activation region.
photon_model : openmc.Model, optional
The OpenMC model to use for photon transport calculations. If None, a
shallow copy of the neutron_model will be created and used.
Attributes
----------
domains : openmc.MeshBase or Sequence[openmc.Cell]
The mesh or a sequence of cells that represent the spatial units over
domains : list of openmc.MeshBase or Sequence[openmc.Cell]
The meshes or a sequence of cells that represent the spatial units over
which the R2S calculation will be performed.
neutron_model : openmc.Model
The OpenMC model used for neutron transport.
@ -101,7 +108,7 @@ class R2SManager:
The OpenMC model used for photon transport calculations.
method : {'mesh-based', 'cell-based'}
Indicates whether the R2S calculation uses mesh elements ('mesh-based')
as the spatial discetization or a list of a cells ('cell-based').
as the spatial discretization or a list of cells ('cell-based').
results : dict
A dictionary that stores results from the R2S calculation.
@ -109,7 +116,7 @@ class R2SManager:
def __init__(
self,
neutron_model: openmc.Model,
domains: openmc.MeshBase | Sequence[openmc.Cell],
domains: openmc.MeshBase | Sequence[openmc.MeshBase] | Sequence[openmc.Cell],
photon_model: openmc.Model | None = None,
):
self.neutron_model = neutron_model
@ -126,9 +133,14 @@ class R2SManager:
self.photon_model = photon_model
if isinstance(domains, openmc.MeshBase):
self.method = 'mesh-based'
self.domains = [domains]
elif isinstance(domains, Sequence) and len(domains) > 0 and \
isinstance(domains[0], openmc.MeshBase):
self.method = 'mesh-based'
self.domains = list(domains)
else:
self.method = 'cell-based'
self.domains = domains
self.domains = list(domains)
self.results = {}
def run(
@ -139,7 +151,7 @@ class R2SManager:
photon_time_indices: Sequence[int] | None = None,
output_dir: PathLike | None = None,
bounding_boxes: dict[int, openmc.BoundingBox] | None = None,
chain_file: PathLike | None = None,
chain_file: PathLike | Chain | None = None,
micro_kwargs: dict | None = None,
mat_vol_kwargs: dict | None = None,
run_kwargs: dict | None = None,
@ -177,9 +189,10 @@ class R2SManager:
Dictionary mapping cell IDs to bounding boxes used for spatial
source sampling in cell-based R2S calculations. Required if method
is 'cell-based'.
chain_file : PathLike, optional
Path to the depletion chain XML file to use during activation. If
not provided, the default configured chain file will be used.
chain_file : PathLike or openmc.deplete.Chain, optional
Path to the depletion chain XML file or depletion chain object to
use during activation. If not provided, the default configured
chain file will be used.
micro_kwargs : dict, optional
Additional keyword arguments passed to
:func:`openmc.deplete.get_microxs_and_flux` during the neutron
@ -206,6 +219,8 @@ class R2SManager:
# consistency (different ranks may have slightly different times)
stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S')
output_dir = Path(comm.bcast(f'r2s_{stamp}'))
else:
output_dir = Path(output_dir)
# Set run_kwargs for the neutron transport step
if micro_kwargs is None:
@ -216,22 +231,35 @@ class R2SManager:
operator_kwargs = {}
run_kwargs.setdefault('output', False)
micro_kwargs.setdefault('run_kwargs', run_kwargs)
# If a chain file is provided, prefer it for steps 1 and 2
if chain_file is not None:
micro_kwargs.setdefault('chain_file', chain_file)
operator_kwargs.setdefault('chain_file', chain_file)
self.step1_neutron_transport(
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
)
self.step2_activation(
timesteps, source_rates, timestep_units, output_dir / 'activation',
operator_kwargs=operator_kwargs
)
self.step3_photon_transport(
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
# DecaySpectrum distributions are resolved in the C++ solver using
# OPENMC_CHAIN_FILE. If a Chain object was passed, write an XML
# representation alongside the R2S outputs.
if isinstance(chain_file, Chain):
output_dir.mkdir(parents=True, exist_ok=True)
chain_path = output_dir / 'chain.xml'
if comm.rank == 0:
chain_file.export_to_xml(chain_path)
comm.barrier()
else:
chain_path = chain_file
chain_context = (
openmc.config.patch('chain_file', chain_path)
if chain_path is not None else nullcontext()
)
with chain_context:
self.step1_neutron_transport(
output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs
)
self.step2_activation(
timesteps, source_rates, timestep_units,
output_dir / 'activation', operator_kwargs=operator_kwargs
)
self.step3_photon_transport(
photon_time_indices, bounding_boxes, output_dir / 'photon_transport',
mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs
)
return output_dir
@ -243,11 +271,13 @@ class R2SManager:
):
"""Run the neutron transport step.
This step computes the material volume fractions on the mesh, creates a
mesh-material filter, and retrieves the fluxes and microscopic cross
sections for each mesh/material combination. This step will populate the
'fluxes' and 'micros' keys in the results dictionary. For a mesh-based
calculation, it will also populate the 'mesh_material_volumes' key.
This step computes the material volume fractions on each mesh, creates
mesh-material filters, and retrieves the fluxes and microscopic cross
sections for each mesh/material combination via a single transport
solve. This step will populate the 'fluxes' and 'micros' keys in the
results dictionary. For a mesh-based calculation, it will also populate
the 'mesh_material_volumes' key (a list of
:class:`~openmc.MeshMaterialVolumes`, one per mesh).
Parameters
----------
@ -266,19 +296,28 @@ class R2SManager:
output_dir.mkdir(parents=True, exist_ok=True)
if self.method == 'mesh-based':
# Compute material volume fractions on the mesh
# Compute material volume fractions on each mesh
if mat_vol_kwargs is None:
mat_vol_kwargs = {}
mat_vol_kwargs.setdefault('bounding_boxes', True)
self.results['mesh_material_volumes'] = mmv = comm.bcast(
self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
# Save results to file
if comm.rank == 0:
mmv.save(output_dir / 'mesh_material_volumes.npz')
mmv_list = []
domain_filters = []
for i, mesh in enumerate(self.domains):
mmv = comm.bcast(
mesh.material_volumes(self.neutron_model, **mat_vol_kwargs))
mmv_list.append(mmv)
# Create mesh-material filter based on what combos were found
domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv)
# Save results to file
if comm.rank == 0:
mmv.save(output_dir / f'mesh_material_volumes_{i}.npz')
# Create mesh-material filter for this mesh
domain_filters.append(
openmc.MeshMaterialFilter.from_volumes(mesh, mmv))
self.results['mesh_material_volumes'] = mmv_list
domains = domain_filters
else:
domains: Sequence[openmc.Cell] = self.domains
@ -357,8 +396,9 @@ class R2SManager:
if self.method == 'mesh-based':
# Get unique material for each (mesh, material) combination
mmv = self.results['mesh_material_volumes']
self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv)
mmv_list = self.results['mesh_material_volumes']
self.results['activation_materials'] = get_activation_materials(
self.neutron_model, mmv_list)
else:
# Create unique material for each cell
activation_mats = openmc.Materials()
@ -468,12 +508,20 @@ class R2SManager:
# photon model if it is different from the neutron model to account for
# potential material changes
if self.method == 'mesh-based' and different_photon_model:
self.results['mesh_material_volumes_photon'] = photon_mmv = comm.bcast(
self.domains.material_volumes(self.photon_model, **mat_vol_kwargs))
if mat_vol_kwargs is None:
mat_vol_kwargs = {}
photon_mmv_list = []
for i, mesh in enumerate(self.domains):
photon_mmv = comm.bcast(
mesh.material_volumes(self.photon_model, **mat_vol_kwargs))
photon_mmv_list.append(photon_mmv)
# Save photon MMV results to file
if comm.rank == 0:
photon_mmv.save(output_dir / 'mesh_material_volumes.npz')
# Save photon MMV results to file
if comm.rank == 0:
photon_mmv.save(
output_dir / f'mesh_material_volumes_{i}.npz')
self.results['mesh_material_volumes_photon'] = photon_mmv_list
if comm.rank == 0:
tally_ids = [tally.id for tally in self.photon_model.tallies]
@ -486,45 +534,30 @@ class R2SManager:
if different_photon_model:
photon_cells = self.photon_model.geometry.get_all_cells()
for time_index in time_indices:
# Create decay photon source
if self.method == 'mesh-based':
self.photon_model.settings.source = \
self.get_decay_photon_source_mesh(time_index)
else:
sources = []
results = self.results['depletion_results']
for cell, original_mat in zip(self.domains, self.results['activation_materials']):
# Skip if the cell is not in the photon model or the
# material has changed
if different_photon_model:
if cell.id not in photon_cells or \
# Determine eligible work items upfront (independent of time index).
if self.method == 'mesh-based':
work_items = self._get_mesh_work_items()
else:
work_items = []
for cell, original_mat in zip(
self.domains, self.results['activation_materials']):
if different_photon_model:
if cell.id not in photon_cells or \
cell.fill.id != photon_cells[cell.id].fill.id:
continue
continue
work_items.append((cell, original_mat, bounding_boxes[cell.id]))
# Get bounding box for the cell
bounding_box = bounding_boxes[cell.id]
# Get activated material composition
activated_mat = results[time_index].get_material(str(original_mat.id))
# Create decay photon source source
space = openmc.stats.Box(*bounding_box)
energy = activated_mat.get_decay_photon_energy()
strength = energy.integral() if energy is not None else 0.0
source = openmc.IndependentSource(
space=space,
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [cell]}
)
sources.append(source)
self.photon_model.settings.source = sources
# Ensure photon transport is enabled in settings
self.photon_model.settings.photon_transport = True
for time_index in time_indices:
# Convert time_index (which may be negative) to a normal index
if time_index < 0:
time_index = len(self.results['depletion_results']) + time_index
time_index += len(self.results['depletion_results'])
# Build decay photon sources and assign to the photon model
sources = self._create_photon_sources(time_index, work_items)
self.photon_model.settings.source = sources
# Run photon transport calculation
photon_dir = Path(output_dir) / f'time_{time_index}'
@ -537,90 +570,107 @@ class R2SManager:
sp.tallies[tally.id] for tally in self.photon_model.tallies
]
def get_decay_photon_source_mesh(
self,
time_index: int = -1
) -> list[openmc.IndependentSource]:
"""Create decay photon source for a mesh-based calculation.
def _get_mesh_work_items(self):
"""Enumerate mesh-based work items across all meshes.
For each mesh element-material combination, an
:class:`~openmc.IndependentSource` is created with a
:class:`~openmc.stats.Box` spatial distribution based on the bounding
box of the material within the mesh element. A material constraint is
also applied so that sampled source sites are limited to the correct
region.
Returns a list of (index_mat, mat_id, bbox) tuples for each eligible
mesh element--material combination, where index_mat is the index into
the activation materials list, mat_id is the material ID, and bbox is
the bounding box for that mesh element--material combination.
When the photon transport model is different from the neutron model, the
photon MeshMaterialVolumes is used to determine whether an (element,
material) combination exists in the photon model.
Returns
-------
list of tuple
Each tuple is (index_mat, mat_id, bbox).
"""
mmv_list = self.results['mesh_material_volumes']
photon_mmv_list = self.results.get('mesh_material_volumes_photon')
work_items = []
index_mat = 0
for mesh_idx, mat_vols in enumerate(mmv_list):
photon_mat_vols = photon_mmv_list[mesh_idx] \
if photon_mmv_list is not None else None
n_elements = mat_vols.num_elements
for index_elem in range(n_elements):
if photon_mat_vols is not None:
photon_materials = {
mat_id
for mat_id, _ in photon_mat_vols.by_element(index_elem)
if mat_id is not None
}
for mat_id, _, bbox in mat_vols.by_element(
index_elem, include_bboxes=True):
if mat_id is None:
continue
if photon_mat_vols is not None \
and mat_id not in photon_materials:
index_mat += 1
continue
work_items.append((index_mat, mat_id, bbox))
index_mat += 1
return work_items
def _create_photon_sources(self, time_index, work_items):
"""Create decay photon sources for a set of regions.
Builds :class:`openmc.IndependentSource` objects with
:class:`openmc.stats.DecaySpectrum` energy distributions that will be
serialized to XML and resolved against the depletion chain by the C++
solver.
Parameters
----------
time_index : int, optional
Time index for the decay photon source. Default is -1 (last time).
time_index : int
Index into depletion results.
work_items : list of tuple
For mesh-based: list of (index_mat, mat_id, bbox).
For cell-based: list of (cell, original_mat, bbox).
Returns
-------
list of openmc.IndependentSource
A list of IndependentSource objects for the decay photons, one for
each mesh element-material combination with non-zero source strength.
Photon sources for each activated region.
"""
mat_dict = self.neutron_model._get_all_materials()
# List to hold all sources
sources = []
# Index in the overall list of activated materials
index_mat = 0
# Get various results from previous steps
mat_vols = self.results['mesh_material_volumes']
step_result = self.results['depletion_results'][time_index]
materials = self.results['activation_materials']
results = self.results['depletion_results']
photon_mat_vols = self.results.get('mesh_material_volumes_photon')
mesh_based = self.method == 'mesh-based'
if mesh_based:
mat_dict = self.neutron_model._get_all_materials()
# Total number of mesh elements
n_elements = mat_vols.num_elements
for index_elem in range(n_elements):
# Determine which materials exist in the photon model for this element
if photon_mat_vols is not None:
photon_materials = {
mat_id
for mat_id, _ in photon_mat_vols.by_element(index_elem)
if mat_id is not None
}
for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
# Skip void volume
if mat_id is None:
continue
# Skip if this material doesn't exist in photon model
if photon_mat_vols is not None and mat_id not in photon_materials:
index_mat += 1
continue
# Get activated material composition
sources = []
for item in work_items:
if mesh_based:
index_mat, domain_id, bbox = item
original_mat = materials[index_mat]
activated_mat = results[time_index].get_material(str(original_mat.id))
domain = mat_dict[domain_id]
else:
cell, original_mat, bbox = item
domain = cell
# Create decay photon source
energy = activated_mat.get_decay_photon_energy()
if energy is not None:
strength = energy.integral()
space = openmc.stats.Box(*bbox)
sources.append(openmc.IndependentSource(
space=space,
energy=energy,
particle='photon',
strength=strength,
constraints={'domains': [mat_dict[mat_id]]}
))
activated_mat = step_result.get_material(str(original_mat.id))
nuclides = activated_mat.get_nuclide_atom_densities()
if not nuclides:
continue
# Increment index of activated material
index_mat += 1
# Eliminate nuclides with zero density
nuclides = {nuclide: density for nuclide, density in nuclides.items()
if density > 0}
energy = openmc.stats.DecaySpectrum(nuclides, activated_mat.volume)
energy.clip(inplace=True)
if not energy.nuclides:
continue
sources.append(openmc.IndependentSource(
space=openmc.stats.Box(bbox.lower_left, bbox.upper_right),
energy=energy,
particle='photon',
constraints={'domains': [domain]},
))
return sources
@ -638,10 +688,13 @@ class R2SManager:
# Load neutron transport results
neutron_dir = path / 'neutron_transport'
if self.method == 'mesh-based':
mmv_file = neutron_dir / 'mesh_material_volumes.npz'
if mmv_file.exists():
self.results['mesh_material_volumes'] = \
openmc.MeshMaterialVolumes.from_npz(mmv_file)
mmv_files = sorted(neutron_dir.glob('mesh_material_volumes*.npz'),
key=lambda p: int(p.stem.split('_')[-1])
if p.stem[-1].isdigit() else 0)
if mmv_files:
self.results['mesh_material_volumes'] = [
openmc.MeshMaterialVolumes.from_npz(f) for f in mmv_files
]
fluxes_file = neutron_dir / 'fluxes.npy'
if fluxes_file.exists():
self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True))
@ -665,10 +718,15 @@ class R2SManager:
# Load photon mesh material volumes if they exist (for mesh-based calculations)
if self.method == 'mesh-based':
photon_mmv_file = photon_dir / 'mesh_material_volumes.npz'
if photon_mmv_file.exists():
self.results['mesh_material_volumes_photon'] = \
openmc.MeshMaterialVolumes.from_npz(photon_mmv_file)
photon_mmv_files = sorted(
photon_dir.glob('mesh_material_volumes*.npz'),
key=lambda p: int(p.stem.split('_')[-1])
if p.stem[-1].isdigit() else 0)
if photon_mmv_files:
self.results['mesh_material_volumes_photon'] = [
openmc.MeshMaterialVolumes.from_npz(f)
for f in photon_mmv_files
]
# Load tally IDs from JSON file
tally_ids_path = photon_dir / 'tally_ids.json'

View file

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

View file

@ -17,7 +17,7 @@ from openmc.mpi import MPI, comm
from .reaction_rates import ReactionRates
VERSION_RESULTS = (1, 2)
VERSION_RESULTS = (1, 3)
__all__ = ["StepResult"]
@ -58,6 +58,8 @@ class StepResult:
proc_time : int
Average time spent depleting a material across all
materials and processes
keff_search_root : float
The root returned by the keff search control.
"""
def __init__(self):
@ -74,6 +76,7 @@ class StepResult:
self.name_list = None
self.data = None
self.keff_search_root = None
def __repr__(self):
t = self.time[0]
@ -154,14 +157,14 @@ class StepResult:
full_burn_list : list of str
List of all burnable material IDs
name_list : list of str, optional
Material names corresponding to materials in burn_list
Material names corresponding to materials in full_burn_list
"""
self.volume = copy.deepcopy(volume)
self.index_nuc = {nuc: i for i, nuc in enumerate(nuc_list)}
self.index_mat = {mat: i for i, mat in enumerate(burn_list)}
self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
self.mat_to_name = dict(zip(burn_list, name_list)) if name_list is not None else {}
self.mat_to_name = dict(zip(full_burn_list, name_list)) if name_list is not None else {}
# Create storage array
self.data = np.zeros((self.n_mat, self.n_nuc))
@ -368,6 +371,10 @@ class StepResult:
"depletion time", (1,), maxshape=(None,),
dtype="float64")
handle.create_dataset(
"keff_search_root", (1,), maxshape=(None,),
dtype="float64")
def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
"""Converts results object into an hdf5 object.
@ -400,6 +407,7 @@ class StepResult:
time_dset = handle["/time"]
source_rate_dset = handle["/source_rate"]
proc_time_dset = handle["/depletion time"]
keff_search_root_dset = handle["/keff_search_root"]
# Get number of results stored
number_shape = list(number_dset.shape)
@ -433,6 +441,10 @@ class StepResult:
proc_shape[0] = new_shape
proc_time_dset.resize(proc_shape)
keff_search_root_shape = list(keff_search_root_dset.shape)
keff_search_root_shape[0] = new_shape
keff_search_root_dset.resize(keff_search_root_shape)
# If nothing to write, just return
if len(self.index_mat) == 0:
return
@ -452,6 +464,7 @@ class StepResult:
proc_time_dset[index] = (
self.proc_time / (comm.size * self.n_hdf5_mats)
)
keff_search_root_dset[index] = self.keff_search_root
@classmethod
def from_hdf5(cls, handle, step):
@ -500,6 +513,10 @@ class StepResult:
if step < proc_time_dset.shape[0]:
results.proc_time = proc_time_dset[step]
if "keff_search_root" in handle:
keff_search_root_dset = handle["/keff_search_root"]
results.keff_search_root = keff_search_root_dset[step]
if results.proc_time is None:
results.proc_time = np.array([np.nan])
@ -554,6 +571,7 @@ class StepResult:
step_ind,
proc_time=None,
write_rates: bool = False,
keff_search_root=None,
path: PathLike = "depletion_results.h5"
):
"""Creates and writes depletion results to disk
@ -578,6 +596,8 @@ class StepResult:
processes.
write_rates : bool, optional
Whether reaction rates should be written to the results file.
keff_search_root : float
The root returned by the keff search control.
path : PathLike
Path to file to write. Defaults to 'depletion_results.h5'.
@ -605,6 +625,7 @@ class StepResult:
results.proc_time = proc_time
if results.proc_time is not None:
results.proc_time = comm.reduce(proc_time, op=MPI.SUM)
results.keff_search_root = keff_search_root
if not Path(path).is_file():
Path(path).parent.mkdir(parents=True, exist_ok=True)

View file

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

View file

@ -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':

View file

@ -268,10 +268,11 @@ class MeshBase(IDManagerMixin, ABC):
return string
def _volume_dim_check(self):
if self.n_dimension != 3 or \
any([d == 0 for d in self.dimension]):
raise RuntimeError(f'Mesh {self.id} is not 3D. '
'Volumes cannot be provided.')
if any(d == 0 for d in self.dimension):
raise RuntimeError(
f'Mesh {self.id} has a zero-size dimension. '
'Volumes cannot be provided.'
)
@classmethod
def from_hdf5(cls, group: h5py.Group):
@ -936,6 +937,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 +1272,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 +1797,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 +2105,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 +2140,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 +2170,6 @@ class CylindricalMesh(StructuredMesh):
origin=origin
)
return mesh
def to_xml_element(self):
"""Return XML representation of the mesh
@ -2385,39 +2477,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 +2517,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 +2533,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)

View file

@ -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)

View file

@ -265,6 +265,46 @@ class Model:
denom_tally = openmc.Tally(name='IFP denominator')
denom_tally.scores = ['ifp-denominator']
self.tallies.append(denom_tally)
# TODO: This should also be incorporated into lower-level calls in
# settings.py, but it requires information about the tallies currently
# on the active Model
def _assign_fw_cadis_tally_IDs(self):
# Verify that all tallies assigned as targets on WeightWindowGenerators
# exist within model.tallies. If this is the case, convert the .targets
# attribute of each WeightWindowGenerator to a sequence of tally IDs.
if len(self.settings.weight_window_generators) == 0:
return
# List of valid tally IDs
reference_tally_ids = np.asarray([tal.id for tal in self.tallies])
for wwg in self.settings.weight_window_generators:
# Only proceeds if the "targets" attribute is an openmc.Tallies,
# which means it hasn't been checked against model.tallies.
if isinstance(wwg.targets, openmc.Tallies):
id_vec = []
for tal in wwg.targets:
# check against model tallies for equivalence
id_next = None
for reference_tal in self.tallies:
if tal == reference_tal:
id_next = reference_tal.id
break
if id_next == None:
raise RuntimeError(
f'Local FW-CADIS target tally {tal.id} not found on model.tallies!')
else:
id_vec.append(id_next)
wwg.targets = id_vec
elif isinstance(wwg.targets, np.ndarray):
invalid = wwg.targets[~np.isin(wwg.targets, reference_tally_ids)]
if len(invalid) > 0:
raise RuntimeError(
f'Local FW-CADIS target tally IDs {invalid} not found on model.tallies!')
@classmethod
def from_xml(
@ -576,6 +616,7 @@ class Model:
if not d.is_dir():
d.mkdir(parents=True, exist_ok=True)
self._assign_fw_cadis_tally_IDs()
self.settings.export_to_xml(d)
self.geometry.export_to_xml(d, remove_surfs=remove_surfs)
@ -634,6 +675,9 @@ class Model:
"set the Geometry.merge_surfaces attribute instead.")
self.geometry.merge_surfaces = True
# Link FW-CADIS WeightWindowGenerator target tallies, if present
self._assign_fw_cadis_tally_IDs()
# provide a memo to track which meshes have been written
mesh_memo = set()
settings_element = self.settings.to_xml_element(mesh_memo)

View file

@ -236,6 +236,9 @@ class Settings:
stabilization, which may be desirable as stronger diagonal stabilization
also tends to dampen the convergence rate of the solver, thus requiring
more iterations to converge.
:adjoint_source:
Source object used to define localized adjoint source/detector response
function.
.. versionadded:: 0.15.0
resonance_scattering : dict
@ -1421,6 +1424,14 @@ class Settings:
cv.check_type('diagonal stabilization rho', value, Real)
cv.check_greater_than('diagonal stabilization rho',
value, 0.0, True)
elif key == 'adjoint_source':
if not isinstance(value, MutableSequence):
value = [value]
for source in value:
if not isinstance(source, SourceBase):
raise ValueError(
f'Invalid adjoint source type: {type(source)}. '
'Expected openmc.SourceBase.')
else:
raise ValueError(f'Unable to set random ray to "{key}" which is '
'unsupported by OpenMC')
@ -1973,11 +1984,12 @@ class Settings:
element = ET.SubElement(root, "random_ray")
for key, value in self._random_ray.items():
if key == 'ray_source' and isinstance(value, SourceBase):
subelement = ET.SubElement(element, 'ray_source')
source_element = value.to_xml_element()
if source_element.find('bias') is not None:
raise RuntimeError(
"Ray source distributions should not be biased.")
element.append(source_element)
subelement.append(source_element)
elif key == 'source_region_meshes':
subelement = ET.SubElement(element, 'source_region_meshes')
@ -1995,8 +2007,20 @@ class Settings:
path = f"./mesh[@id='{mesh.id}']"
if root.find(path) is None:
root.append(mesh.to_xml_element())
if mesh_memo is not None:
if mesh_memo is not None:
mesh_memo.add(mesh.id)
elif key == 'adjoint_source':
subelement = ET.SubElement(element, 'adjoint_source')
# Check that all entries are valid SourceBase instances, in case
# the random_ray setter was not used to populate dict entries.
if not isinstance(value, MutableSequence):
value = [value]
for source in value:
if not isinstance(source, SourceBase):
raise ValueError(
f'Invalid adjoint source type: {type(source)}. '
'Expected openmc.SourceBase.')
subelement.append(source.to_xml_element())
elif isinstance(value, bool):
subelement = ET.SubElement(element, key)
subelement.text = str(value).lower()
@ -2443,8 +2467,9 @@ class Settings:
for child in elem:
if child.tag in ('distance_inactive', 'distance_active', 'diagonal_stabilization_rho'):
self.random_ray[child.tag] = float(child.text)
elif child.tag == 'source':
source = SourceBase.from_xml_element(child)
elif child.tag == 'ray_source':
source_element = child.find('source')
source = SourceBase.from_xml_element(source_element)
if child.find('bias') is not None:
raise RuntimeError(
"Ray source distributions should not be biased.")
@ -2461,6 +2486,12 @@ class Settings:
self.random_ray['adjoint'] = (
child.text in ('true', '1')
)
elif child.tag == 'adjoint_source':
self.random_ray['adjoint_source'] = []
for subelem in child.findall('source'):
src = SourceBase.from_xml_element(subelem)
# add newly constructed source object to the list
self.random_ray['adjoint_source'].append(src)
elif child.tag == 'sample_method':
self.random_ray['sample_method'] = child.text
elif child.tag == 'source_region_meshes':

View file

@ -2,9 +2,10 @@ from __future__ import annotations
from abc import ABC, abstractmethod
from collections import defaultdict
from collections.abc import Iterable, Sequence
from copy import deepcopy
from math import sqrt, pi, exp
from functools import cache
from math import sqrt, pi, exp, log
from numbers import Real
from pathlib import Path
from warnings import warn
import lxml.etree as ET
@ -14,6 +15,8 @@ from scipy.special import exprel, hyp1f1, lambertw
import scipy
import openmc.checkvalue as cv
from openmc.data import atomic_mass, NEUTRON_MASS
import openmc.data
from .._xml import get_elem_list, get_text
from ..mixin import EqualityMixin
@ -124,6 +127,8 @@ class Univariate(EqualityMixin, ABC):
return Legendre.from_xml_element(elem)
elif distribution == 'mixture':
return Mixture.from_xml_element(elem)
elif distribution == 'decay_spectrum':
return DecaySpectrum.from_xml_element(elem)
@abstractmethod
def _sample_unbiased(self, n_samples: int = 1, seed: int | None = None):
@ -1277,6 +1282,138 @@ def Muir(*args, **kwargs):
return muir(*args, **kwargs)
def fusion_neutron_spectrum(
ion_temp: float,
reactants: str = 'DD',
bias: Univariate | None = None
) -> Normal:
r"""Return a Gaussian energy distribution for fusion neutron emission.
Computes the mean energy and spectral width of the neutron energy spectrum
from thermonuclear fusion reactions in a plasma with Maxwellian ion velocity
distributions. The mean neutron energy is calculated as
.. math::
\langle E_n \rangle = E_0 + \Delta E_\text{th}
where :math:`E_0` is the neutron energy at zero ion temperature and
:math:`\Delta E_\text{th}` is the thermal peak shift due to the motion of
the reacting ions. The spectral width is characterized by the FWHM:
.. math::
W_{1/2} = \omega_0 (1 + \delta_\omega) \sqrt{T_i}
where :math:`\omega_0` is the width at the :math:`T_i \to 0` limit and
:math:`\delta_\omega` is a temperature-dependent correction term. Both
:math:`\Delta E_\text{th}` and :math:`\delta_\omega` are evaluated using
interpolation formulas from `Ballabio et al.
<https://doi.org/10.1088/0029-5515/38/11/310>`_: Table III for :math:`0 <
T_i \le 40` keV and Table IV for :math:`40 < T_i < 100` keV. The returned
distribution is a normal (Gaussian) approximation to the spectrum.
.. versionadded:: 0.15.4
Parameters
----------
ion_temp : float
Ion temperature of the plasma in [eV].
reactants : {'DD', 'DT'}
Fusion reactants. 'DD' corresponds to the D(d,n)\ :sup:`3`\ He reaction
and 'DT' to the T(d,n)\ :math:`\alpha` reaction.
bias : openmc.stats.Univariate, optional
Distribution for biased sampling.
Returns
-------
openmc.stats.Normal
Normal distribution with mean and standard deviation corresponding to
the first and second moments of the fusion neutron energy spectrum. Both
the mean and standard deviation are in [eV].
"""
if ion_temp < 0.0 or ion_temp > 100e3:
raise ValueError("Ion temperature must be between 0 and 100 keV.")
# Formulas from doi:10.1088/0029-5515/38/11/310
mn = NEUTRON_MASS
md = atomic_mass('H2')
ev_per_c2 = 931.49410372*1e6
if reactants == 'DD':
mhe3 = atomic_mass('He3')
Q = (md + md - mhe3 - mn)*ev_per_c2
E_n = mhe3/(mhe3 + mn)*Q
w0 = 82.542
# Low-T constants for peak shift (Table III)
a1 = 4.69515
a2 = -0.040729
a3 = 0.47
a4 = 0.81844
# Low-T constants for width correction (Table III)
b1 = 1.7013e-3
b2 = 0.16888
b3 = 0.49
b4 = 7.9460e-4
# High-T constants for peak shift (Table IV)
a5 = 18.225
a6 = 2.1525
# High-T constants for width correction (Table IV)
b5 = 8.4619e-3
b6 = 8.3241e-4
elif reactants == 'DT':
mt = atomic_mass('H3')
ma = atomic_mass('He4')
Q = (md + mt - ma - mn)*ev_per_c2
E_n = ma/(ma + mn)*Q
w0 = 177.259
# Low-T constants for peak shift (Table III)
a1 = 5.30509
a2 = 2.4736e-3
a3 = 1.84
a4 = 1.3818
# Low-T constants for width correction (Table III)
b1 = 5.1068e-4
b2 = 7.6223e-3
b3 = 1.78
b4 = 8.7691e-5
# High-T constants for peak shift (Table IV)
a5 = 37.771
a6 = 0.92181
# High-T constants for width correction (Table IV)
b5 = 2.0199e-3
b6 = 5.9501e-5
else:
raise ValueError("Invalid reactants specified. Must be 'DD' or 'DT'.")
# Ion temperature in keV
T = ion_temp * 1e-3
if T <= 40.0:
# Low-temperature interpolation (Table III, 0 < T_i <= 40 keV)
Delta_E = a1/(1 + a2*T**a3)*T**(2/3) + a4*T
delta_w = b1/(1 + b2*T**b3)*T**(2/3) + b4*T
else:
# High-temperature interpolation (Table IV, 40 < T_i < 100 keV)
Delta_E = a5 + a6*T
delta_w = b5 + b6*T
# Calculate FWHM
fwhm = (w0*(1 + delta_w) * sqrt(T))*1e3
sigma = fwhm / (2*sqrt(2*log(2)))
return Normal(E_n + Delta_E * 1e3, sigma, bias=bias)
class Tabular(Univariate):
"""Piecewise continuous probability distribution.
@ -2064,6 +2201,310 @@ class Mixture(Univariate):
return new_dist
class DecaySpectrum(Univariate):
"""Energy distribution from decay photon spectra of a mixture of nuclides.
This distribution stores nuclide names, their atom densities, and the volume
of the region. When written to XML and read by the C++ solver, the nuclide
names are resolved against the depletion chain to obtain the decay photon
energy spectra and decay constants. The resulting distribution is a mixture
of per-nuclide photon spectra weighted by absolute activity. The volume is
necessary so that the C++ solver can compute the total photon emission rate
in [photons/s], which is used as the source strength.
.. versionadded:: 0.15.4
Parameters
----------
nuclides : dict
Dictionary mapping nuclide name (str) to atom density (float) in units
of [atom/b-cm].
volume : float
Volume of the source region in [cm³]. Used together with atom densities
to compute the absolute photon emission rate.
Attributes
----------
nuclides : dict
Dictionary mapping nuclide name to atom density in [atom/b-cm].
volume : float
Volume of the source region in [cm³].
"""
def __init__(self, nuclides: dict[str, float], volume: float):
super().__init__(bias=None)
self._dist_cache = None
self._dist_cache_key = None
self.nuclides = nuclides
self.volume = volume
def __len__(self):
return len(self.nuclides)
@property
def nuclides(self):
return self._nuclides
@nuclides.setter
def nuclides(self, nuclides):
cv.check_type('nuclides', nuclides, dict)
for name, density in nuclides.items():
cv.check_type('nuclide name', name, str)
cv.check_type(f'atom density for {name}', density, Real)
cv.check_greater_than(f'atom density for {name}', density, 0.0, True)
self._nuclides = dict(nuclides)
self._dist_cache = None
self._dist_cache_key = None
@property
def volume(self):
return self._volume
@volume.setter
def volume(self, volume):
cv.check_type('volume', volume, Real)
cv.check_greater_than('volume', volume, 0.0)
self._volume = float(volume)
self._dist_cache = None
self._dist_cache_key = None
@staticmethod
def _chain_file_cache_key():
"""Return a hashable key for the active depletion chain."""
chain_file = openmc.config.get('chain_file')
if chain_file is None:
return None
path = Path(chain_file).resolve()
try:
stat = path.stat()
except OSError:
return (path, None, None)
return (path, stat.st_mtime, stat.st_size)
def to_distribution(self):
"""Convert to a concrete distribution using decay chain data.
Builds a combined photon energy distribution by looking up each nuclide
in the depletion chain via :func:`openmc.data.decay_photon_energy` and
weighting by absolute atom count (``density * 1e24 * volume``). The
result is cached on the object; the cache is invalidated automatically
when :attr:`nuclides` or :attr:`volume` are reassigned.
Requires ``openmc.config['chain_file']`` to be set.
Returns
-------
openmc.stats.Univariate or None
Combined photon energy distribution, or ``None`` if no nuclide in
:attr:`nuclides` has a photon source in the chain.
"""
chain_key = self._chain_file_cache_key()
if self._dist_cache is not None and self._dist_cache_key == chain_key:
return self._dist_cache
dists = []
weights = []
for name, density in self.nuclides.items():
dist = openmc.data.decay_photon_energy(name)
if dist is not None:
dists.append(dist)
weights.append(density * 1e24 * self.volume)
if not dists:
return None
self._dist_cache = combine_distributions(dists, weights)
self._dist_cache_key = chain_key
return self._dist_cache
def to_xml_element(self, element_name: str):
"""Return XML representation of the decay photon distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : lxml.etree._Element
XML element containing decay photon distribution data
"""
element = ET.Element(element_name)
element.set("type", "decay_spectrum")
element.set("volume", str(self.volume))
nuclides = ET.SubElement(element, "nuclides")
nuclides.text = ' '.join(self.nuclides)
parameters = ET.SubElement(element, "parameters")
parameters.text = ' '.join(str(density) for density in self.nuclides.values())
return element
@classmethod
def from_xml_element(cls, elem: ET.Element):
"""Generate decay photon distribution from an XML element
Parameters
----------
elem : lxml.etree._Element
XML element
Returns
-------
openmc.stats.DecaySpectrum
Decay photon distribution generated from XML element
"""
volume = float(elem.get('volume'))
names = get_elem_list(elem, 'nuclides', str)
densities = get_elem_list(elem, 'parameters', float)
nuclides = dict(zip(names, densities))
return cls(nuclides, volume)
def _sample_unbiased(self, n_samples=1, seed=None):
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum._sample_unbiased requires chain data but none "
"was found. Ensure openmc.config['chain_file'] is set and the "
"chain contains photon sources for the nuclides present."
)
return dist.sample(n_samples, seed)[0]
def integral(self):
"""Return integral of the distribution
Returns the total photon emission rate in [photons/s] by delegating to
:meth:`to_distribution`. Returns ``0.0`` when no chain data is
available (e.g., ``openmc.config['chain_file']`` is not set).
Returns
-------
float
Total photon emission rate in [photons/s], or ``0.0`` if chain
data is unavailable.
"""
try:
dist = self.to_distribution()
except Exception:
return 0.0
if dist is None:
return 0.0
return dist.integral()
@staticmethod
@cache
def _photon_integral(nuclide: str, chain_key) -> float | None:
"""Return the per-atom photon emission integral for a nuclide"""
dist = openmc.data.decay_photon_energy(nuclide)
return dist.integral() if dist is not None else None
def clip(self, tolerance: float = 1e-9, inplace: bool = False):
"""Remove nuclides with negligible contribution to photon emission.
Nuclides that are stable or have no photon source in the depletion
chain are removed unconditionally. The remaining nuclides are ranked
by their photon emission rate (proportional to
``atom_density * decay_constant * photon_yield``) and the least
important are discarded until the cumulative discarded fraction of the
total emission rate exceeds *tolerance*.
Requires ``openmc.config['chain_file']`` to be set.
Parameters
----------
tolerance : float
Maximum fraction of total photon emission rate that may be
discarded.
inplace : bool
Whether to modify the current object in-place or return a new one.
Returns
-------
openmc.stats.DecaySpectrum
Distribution with negligible nuclides removed.
"""
# Compute per-nuclide emission rate; drop non-emitters
emitting_names = []
emitting_densities = []
rates = []
chain_key = self._chain_file_cache_key()
for name, density in self.nuclides.items():
integral = DecaySpectrum._photon_integral(name, chain_key)
if integral is None:
continue
emitting_names.append(name)
emitting_densities.append(density)
rates.append(density * self.volume * integral)
if not emitting_names:
new_nuclides = {}
else:
indices = _intensity_clip(rates, tolerance=tolerance)
new_nuclides = {
emitting_names[i]: emitting_densities[i] for i in indices
}
if inplace:
self._nuclides = new_nuclides
self._dist_cache = None
self._dist_cache_key = None
return self
return type(self)(new_nuclides, self.volume)
@property
def support(self):
return (0.0, np.inf)
def evaluate(self, x):
"""Evaluate the probability density at a given value.
Delegates to the combined distribution built from chain data. Raises
``NotImplementedError`` if the combined distribution is a
:class:`~openmc.stats.Mixture` (which does not support
``evaluate()``).
Parameters
----------
x : float
Value at which to evaluate the PDF.
Returns
-------
float
Probability density at *x*.
"""
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum.evaluate requires chain data. Ensure "
"openmc.config['chain_file'] is set."
)
return dist.evaluate(x)
def mean(self):
"""Return the mean of the distribution.
Delegates to the combined distribution built from chain data.
Returns
-------
float
Mean photon energy in [eV].
"""
dist = self.to_distribution()
if dist is None:
raise RuntimeError(
"DecaySpectrum.mean requires chain data. Ensure "
"openmc.config['chain_file'] is set."
)
return dist.mean()
def combine_distributions(
dists: Sequence[Discrete | Tabular | Mixture],
probs: Sequence[float]

View file

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

View file

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

View file

@ -69,7 +69,7 @@ include = ['openmc*']
exclude = ['tests*']
[tool.setuptools.package-data]
"openmc.data.dose" = ["**/*.txt"]
"openmc.data.dose" = ["**/*.txt", "*.h5"]
"openmc.data" = ["*.txt", "*.DAT", "*.json", "*.h5"]
"openmc.lib" = ["libopenmc.dylib", "libopenmc.so"]

View file

@ -98,6 +98,8 @@ vector<unique_ptr<ChainNuclide>> chain_nuclides;
void read_chain_file_xml()
{
free_memory_chain();
char* chain_file_path = std::getenv("OPENMC_CHAIN_FILE");
if (!chain_file_path) {
return;
@ -120,4 +122,10 @@ void read_chain_file_xml()
}
}
void free_memory_chain()
{
data::chain_nuclides.clear();
data::chain_nuclide_map.clear();
}
} // namespace openmc

View file

@ -7,7 +7,9 @@
#include <numeric> // for accumulate
#include <stdexcept> // for runtime_error
#include <string> // for string, stod
#include <unordered_set>
#include "openmc/chain.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
@ -15,6 +17,10 @@
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"
namespace {
std::unordered_set<std::string> decay_spectrum_missing_chain_nuclides;
}
namespace openmc {
//==============================================================================
@ -758,6 +764,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
dist = UPtrDist {new Tabular(node)};
} else if (type == "mixture") {
dist = UPtrDist {new Mixture(node)};
} else if (type == "decay_spectrum") {
dist = UPtrDist {new DecaySpectrum(node)};
} else if (type == "muir") {
openmc::fatal_error(
"'muir' distributions are now specified using the openmc.stats.muir() "
@ -768,4 +776,120 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
return dist;
}
//==============================================================================
// DecaySpectrum implementation
//==============================================================================
DecaySpectrum::DecaySpectrum(pugi::xml_node node)
{
// Read the region volume [cm^3] needed for absolute emission rate
if (!check_for_node(node, "volume"))
fatal_error("DecaySpectrum: 'volume' attribute is required.");
double volume = std::stod(get_node_value(node, "volume"));
// Read nuclide names and atom densities from XML
vector<int> nuclide_indices;
vector<double> atoms;
auto names = get_node_array<std::string>(node, "nuclides");
auto densities = get_node_array<double>(node, "parameters");
if (names.size() != densities.size()) {
fatal_error("DecaySpectrum nuclides and parameters must have the same "
"length.");
}
for (size_t i = 0; i < names.size(); ++i) {
const auto& name = names[i];
double density = densities[i];
// Look up nuclide in the depletion chain
auto it = data::chain_nuclide_map.find(name);
if (it == data::chain_nuclide_map.end()) {
if (decay_spectrum_missing_chain_nuclides.insert(name).second) {
warning("Nuclide '" + name +
"' appears in a DecaySpectrum source but is not present in "
"the depletion chain; it will be ignored.");
}
continue;
}
int nuclide_index = it->second;
const auto& chain_nuc = data::chain_nuclides[nuclide_index];
const Distribution* photon_dist = chain_nuc->photon_energy();
if (!photon_dist)
continue;
// Skip non-positive densities and warn if negative
if (density <= 0.0) {
if (density < 0.0) {
warning("Nuclide '" + name +
"' has a negative density in a DecaySpectrum source; it will "
"be ignored.");
}
continue;
}
// atoms = density [atom/b-cm] * 1e24 [b/cm^2] * volume [cm^3]
double atoms_i = density * 1.0e24 * volume;
nuclide_indices.push_back(nuclide_index);
atoms.push_back(atoms_i);
}
init(std::move(nuclide_indices), atoms);
}
void DecaySpectrum::init(
vector<int> nuclide_indices, const vector<double>& atoms)
{
if (nuclide_indices.size() != atoms.size()) {
fatal_error("DecaySpectrum nuclide index and atoms arrays must have "
"the same length.");
}
vector<double> probs;
probs.reserve(nuclide_indices.size());
for (size_t i = 0; i < nuclide_indices.size(); ++i) {
// Distribution integral is in [photons/s/atom]; multiplying by atoms gives
// the total emission rate [photons/s] for this nuclide.
const auto* dist =
data::chain_nuclides[nuclide_indices[i]]->photon_energy();
probs.push_back(atoms[i] * dist->integral());
}
nuclide_indices_ = std::move(nuclide_indices);
integral_ = std::accumulate(probs.begin(), probs.end(), 0.0);
if (nuclide_indices_.empty() || integral_ <= 0.0) {
fatal_error("DecaySpectrum source did not resolve any nuclides with decay "
"photon spectra and positive atom densities. Ensure "
"OPENMC_CHAIN_FILE is set and matches the nuclides in the "
"source definition.");
}
di_.assign(probs);
}
DecaySpectrum::Sample DecaySpectrum::sample_with_parent(uint64_t* seed) const
{
size_t idx = di_.sample(seed);
int parent_nuclide = nuclide_indices_[idx];
const auto* dist = data::chain_nuclides[parent_nuclide]->photon_energy();
auto [energy, weight] = dist->sample(seed);
return {energy, weight, parent_nuclide};
}
std::pair<double, double> DecaySpectrum::sample(uint64_t* seed) const
{
auto sample = sample_with_parent(seed);
return {sample.energy, sample.weight};
}
double DecaySpectrum::integral() const
{
return integral_;
}
double DecaySpectrum::sample_unbiased(uint64_t* seed) const
{
return sample_with_parent(seed).energy;
}
} // namespace openmc

View file

@ -2,6 +2,7 @@
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/chain.h"
#include "openmc/cmfd_solver.h"
#include "openmc/collision_track.h"
#include "openmc/constants.h"
@ -44,6 +45,7 @@ void free_memory()
free_memory_photon();
free_memory_settings();
free_memory_thermal();
free_memory_chain();
library_clear();
nuclides_clear();
free_memory_source();

View file

@ -405,6 +405,10 @@ bool read_model_xml()
write_message(
fmt::format("Reading model XML file '{}' ...", model_filename), 5);
// Read chain data before settings so DecaySpectrum source distributions can
// resolve nuclides while sources are constructed.
read_chain_file_xml();
read_settings_xml(settings_root);
// If other XML files are present, display warning
@ -420,9 +424,6 @@ bool read_model_xml()
}
}
// Read data from chain file
read_chain_file_xml();
// Read materials and cross sections
if (!check_for_node(root, "materials")) {
fatal_error(fmt::format(
@ -475,14 +476,15 @@ bool read_model_xml()
void read_separate_xml_files()
{
// Read chain data before settings so DecaySpectrum source distributions can
// resolve nuclides while sources are constructed.
read_chain_file_xml();
read_settings_xml();
if (settings::run_mode != RunMode::PLOTTING) {
read_cross_sections_xml();
}
// Read data from chain file
read_chain_file_xml();
read_materials_xml();
read_geometry_xml();

View file

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

View file

@ -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) {

View file

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

View file

@ -321,7 +321,6 @@ void print_build_info()
std::string png(n);
std::string profiling(n);
std::string coverage(n);
std::string mcpl(n);
std::string uwuw(n);
std::string strict_fp(n);
@ -337,9 +336,6 @@ void print_build_info()
#ifdef OPENMC_LIBMESH_ENABLED
libmesh = y;
#endif
#ifdef OPENMC_MCPL
mcpl = y;
#endif
#ifdef USE_LIBPNG
png = y;
#endif
@ -369,7 +365,6 @@ void print_build_info()
fmt::print("PNG support: {}\n", png);
fmt::print("DAGMC support: {}\n", dagmc);
fmt::print("libMesh support: {}\n", libmesh);
fmt::print("MCPL support: {}\n", mcpl);
fmt::print("Coverage testing: {}\n", coverage);
fmt::print("Profiling flags: {}\n", profiling);
fmt::print("UWUW support: {}\n", uwuw);

View file

@ -14,6 +14,7 @@
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/hdf5_interface.h"
#include "openmc/lattice.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
@ -47,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;
}
}
@ -302,8 +310,6 @@ void Particle::event_cross_surface()
surface() = boundary().surface();
n_coord() = boundary().coord_level();
const auto& surf {*model::surfaces[surface_index()].get()};
if (boundary().lattice_translation()[0] != 0 ||
boundary().lattice_translation()[1] != 0 ||
boundary().lattice_translation()[2] != 0) {
@ -312,7 +318,23 @@ void Particle::event_cross_surface()
bool verbose = settings::verbosity >= 10 || trace();
cross_lattice(*this, boundary(), verbose);
event() = TallyEvent::LATTICE;
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
auto& lat {*model::lattices[lowest_coord().lattice()]};
bool is_valid;
Direction normal =
lat.get_normal(boundary().lattice_translation(), is_valid);
if (is_valid) {
normal /= normal.norm();
score_surface_tally(*this, model::active_surface_tallies, normal);
}
}
} else {
const auto& surf {*model::surfaces[surface_index()].get()};
// Particle crosses surface
// If BC, add particle to surface source before crossing surface
if (surf.surf_source_ && surf.bc_) {
@ -327,10 +349,13 @@ void Particle::event_cross_surface()
apply_weight_windows(*this);
}
event() = TallyEvent::SURFACE;
}
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
score_surface_tally(*this, model::active_surface_tallies, surf);
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
Direction normal = surf.normal(r());
normal /= normal.norm();
score_surface_tally(*this, model::active_surface_tallies, normal);
}
}
}
@ -681,7 +706,9 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
// with a mesh boundary
if (!model::active_surface_tallies.empty()) {
score_surface_tally(*this, model::active_surface_tallies, surf);
Direction normal = surf.normal(r());
normal /= normal.norm();
score_surface_tally(*this, model::active_surface_tallies, normal);
}
if (!model::active_meshsurf_tallies.empty()) {

View file

@ -31,9 +31,11 @@ RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ {
RandomRayVolumeEstimator::HYBRID};
bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};
bool FlatSourceDomain::adjoint_ {false};
bool FlatSourceDomain::fw_cadis_local_ {false};
double FlatSourceDomain::diagonal_stabilization_rho_ {1.0};
std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
FlatSourceDomain::mesh_domain_map_;
std::vector<size_t> FlatSourceDomain::fw_cadis_local_targets_;
FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_)
{
@ -1000,7 +1002,9 @@ void FlatSourceDomain::output_to_vtk() const
void FlatSourceDomain::apply_external_source_to_source_region(
int src_idx, SourceRegionHandle& srh)
{
auto s = model::external_sources[src_idx].get();
auto s = (adjoint_ && !model::adjoint_sources.empty())
? model::adjoint_sources[src_idx].get()
: model::external_sources[src_idx].get();
auto is = dynamic_cast<IndependentSource*>(s);
auto discrete = dynamic_cast<Discrete*>(is->energy());
double strength_factor = is->strength();
@ -1071,13 +1075,17 @@ void FlatSourceDomain::count_external_source_regions()
}
}
void FlatSourceDomain::convert_external_sources()
void FlatSourceDomain::convert_external_sources(bool use_adjoint_sources)
{
// Determine whether forward or (local) adjoint sources are desired
const auto& sources =
use_adjoint_sources ? model::adjoint_sources : model::external_sources;
// Loop over external sources
for (int es = 0; es < model::external_sources.size(); es++) {
for (int es = 0; es < sources.size(); es++) {
// Extract source information
Source* s = model::external_sources[es].get();
Source* s = sources[es].get();
IndependentSource* is = dynamic_cast<IndependentSource*>(s);
Discrete* energy = dynamic_cast<Discrete*>(is->energy());
const std::unordered_set<int32_t>& domain_ids = is->domain_ids();
@ -1223,7 +1231,7 @@ void FlatSourceDomain::flatten_xs()
}
}
void FlatSourceDomain::set_adjoint_sources()
void FlatSourceDomain::set_fw_adjoint_sources()
{
// Set the adjoint external source to 1/forward_flux. If the forward flux is
// negative, zero, or extremely close to zero, set the adjoint source to zero,
@ -1252,6 +1260,10 @@ void FlatSourceDomain::set_adjoint_sources()
source_regions_.external_source(sr, g) = 0.0;
} else {
source_regions_.external_source(sr, g) = 1.0 / flux;
if (!std::isfinite(source_regions_.external_source(sr, g))) {
// If the flux is NaN or Inf, set the adjoint source to zero
source_regions_.external_source(sr, g) = 0.0;
}
}
if (flux > 0.0) {
source_regions_.external_source_present(sr) = 1;
@ -1283,6 +1295,7 @@ void FlatSourceDomain::set_adjoint_sources()
source_regions_.external_source_present(sr) = 0;
}
}
// Divide the fixed source term by sigma t (to save time when applying each
// iteration)
#pragma omp parallel for
@ -1297,8 +1310,86 @@ void FlatSourceDomain::set_adjoint_sources()
sigma_t_[(material * ntemperature_ + temp) * negroups_ + g] *
source_regions_.density_mult(sr);
source_regions_.external_source(sr, g) /= sigma_t;
if (!std::isfinite(source_regions_.external_source(sr, g))) {
// If the flux is NaN or Inf, set the adjoint source to zero
source_regions_.external_source(sr, g) = 0.0;
}
}
}
if (fw_cadis_local_) {
// Only external sources that have a non-mesh type tally task should remain
// non-zero. Everything else gets zero'd out.
#pragma omp parallel for
for (int64_t sr = 0; sr < n_source_regions(); sr++) {
// If there is already no external source, don't need to do anything
if (source_regions_.external_source_present(sr) == 0) {
continue;
}
// If there is an adjoint source term here, then we need to check it.
// We will track if ANY group has a valid local FW-CADIS source term
bool has_any_sources = false;
// Now, loop over groups
for (int g = 0; g < negroups_; g++) {
// If there are no tally tasks associated with this source element
// then it is not a local FW-CADIS source, so we continue to the next
// group
if (source_regions_.tally_task(sr, g).empty()) {
source_regions_.external_source(sr, g) = 0.0;
continue;
}
// If there are tally tasks, we can through them and check if
// any of them are local FW-CADIS targets.
// We track if ANY of the tasks are local FW-CADIS target tallies
bool local_fw_cadis_target_region = false;
// Now we loop through
for (const auto& task : source_regions_.tally_task(sr, g)) {
Tally& tally {*model::tallies[task.tally_idx]};
const auto t_id = tally.id();
// Search for target tallies
if (std::find(fw_cadis_local_targets_.begin(),
fw_cadis_local_targets_.end(),
t_id) != fw_cadis_local_targets_.end()) {
local_fw_cadis_target_region = true;
break;
}
}
// If ANY of the tasks is a local FW-CADIS target,
// Then we keep the source term and set that this
// source region has a valid FW-CADIS source term.
// Otherwise, we zero out the source term.
if (local_fw_cadis_target_region) {
has_any_sources = true;
} else {
source_regions_.external_source(sr, g) = 0.0;
}
} // End loop over groups
// If there were any valid FW-CADIS source terms for any
// of the groups, then the SR as a whole counts as a source
if (has_any_sources) {
source_regions_.external_source_present(sr) = 1;
} else {
source_regions_.external_source_present(sr) = 0;
}
} // End loop over source regions
} // End local FW-CADIS logic
}
void FlatSourceDomain::set_local_adjoint_sources()
{
// Set the external source to user-specified adjoint sources.
convert_external_sources(true);
}
void FlatSourceDomain::transpose_scattering_matrix()

View file

@ -168,14 +168,12 @@ void validate_random_ray_inputs()
"constrained by domain id (cell, material, or universe) in "
"random ray mode.");
} else if (is->domain_ids().size() > 0 && sp) {
// If both a domain constraint and a non-default point source location
// are specified, notify user that domain constraint takes precedence.
if (sp->r().x == 0.0 && sp->r().y == 0.0 && sp->r().z == 0.0) {
warning("Fixed source has both a domain constraint and a point "
"type spatial distribution. The domain constraint takes "
"precedence in random ray mode -- point source coordinate "
"will be ignored.");
}
// If both a domain constraint and a point source location are
// specified, notify user that domain constraint takes precedence.
warning("Fixed source has both a domain constraint and a point "
"type spatial distribution. The domain constraint takes "
"precedence in random ray mode -- point source coordinate "
"will be ignored.");
}
// Check that a discrete energy distribution was used
@ -189,6 +187,56 @@ void validate_random_ray_inputs()
}
}
// Validate adjoint sources
///////////////////////////////////////////////////////////////////
if (FlatSourceDomain::adjoint_ && !model::adjoint_sources.empty()) {
for (int i = 0; i < model::adjoint_sources.size(); i++) {
Source* s = model::adjoint_sources[i].get();
// Check for independent source
IndependentSource* is = dynamic_cast<IndependentSource*>(s);
if (!is) {
fatal_error(
"Only IndependentSource adjoint source types are allowed in "
"random ray mode");
}
// Check for isotropic source
UnitSphereDistribution* angle_dist = is->angle();
Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
if (!id) {
fatal_error(
"Invalid source definition -- only isotropic adjoint sources are "
"allowed in random ray mode.");
}
// Validate that a domain ID was specified OR that it is a point source
auto sp = dynamic_cast<SpatialPoint*>(is->space());
if (is->domain_ids().size() == 0 && !sp) {
fatal_error("Adjoint sources must be point source or spatially "
"constrained by domain id (cell, material, or universe) in "
"random ray mode.");
} else if (is->domain_ids().size() > 0 && sp) {
// If both a domain constraint and a point source location are
// specified, notify user that domain constraint takes precedence.
warning("Adjoint source has both a domain constraint and a point "
"type spatial distribution. The domain constraint takes "
"precedence in random ray mode -- point source coordinate "
"will be ignored.");
}
// Check that a discrete energy distribution was used
Distribution* d = is->energy();
Discrete* dd = dynamic_cast<Discrete*>(d);
if (!dd) {
fatal_error(
"Only discrete (multigroup) energy distributions are allowed for "
"adjoint sources in random ray mode.");
}
}
}
// Validate plotting files
///////////////////////////////////////////////////////////////////
for (int p = 0; p < model::plots.size(); p++) {
@ -232,20 +280,22 @@ void validate_random_ray_inputs()
warning(
"Linear sources may result in negative fluxes in small source regions "
"generated by mesh subdivision. Negative sources may result in low "
"quality FW-CADIS weight windows. We recommend you use flat source mode "
"when generating weight windows with an overlaid mesh tally.");
"quality FW-CADIS weight windows. We recommend you use flat source "
"mode when generating weight windows with an overlaid mesh tally.");
}
}
void print_adjoint_header()
void openmc_finalize_random_ray()
{
if (!FlatSourceDomain::adjoint_)
// If we're going to do an adjoint simulation afterwards, report that this
// is the initial forward flux solve.
header("FORWARD FLUX SOLVE", 3);
else
// Otherwise report that we are doing the adjoint simulation
header("ADJOINT FLUX SOLVE", 3);
FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID;
FlatSourceDomain::volume_normalized_flux_tallies_ = false;
FlatSourceDomain::adjoint_ = false;
FlatSourceDomain::fw_cadis_local_ = false;
FlatSourceDomain::fw_cadis_local_targets_.clear();
FlatSourceDomain::mesh_domain_map_.clear();
RandomRay::ray_source_.reset();
RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
}
//==============================================================================
@ -278,16 +328,6 @@ RandomRaySimulation::RandomRaySimulation()
// Convert OpenMC native MGXS into a more efficient format
// internal to the random ray solver
domain_->flatten_xs();
// Check if adjoint calculation is needed. If it is, we will run the forward
// calculation first and then the adjoint calculation later.
adjoint_needed_ = FlatSourceDomain::adjoint_;
// Adjoint is always false for the forward calculation
FlatSourceDomain::adjoint_ = false;
// The first simulation is run after initialization
is_first_simulation_ = true;
}
void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
@ -295,30 +335,52 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
domain_->apply_meshes();
if (settings::run_mode == RunMode::FIXED_SOURCE) {
// Transfer external source user inputs onto random ray source regions
domain_->convert_external_sources();
domain_->convert_external_sources(false);
domain_->count_external_source_regions();
}
}
void RandomRaySimulation::prepare_fixed_sources_adjoint()
void RandomRaySimulation::prepare_fw_fixed_sources_adjoint()
{
// Prepare adjoint fixed sources using forward flux
domain_->source_regions_.adjoint_reset();
if (settings::run_mode == RunMode::FIXED_SOURCE) {
domain_->set_adjoint_sources();
domain_->set_fw_adjoint_sources();
}
}
void RandomRaySimulation::prepare_adjoint_simulation()
void RandomRaySimulation::prepare_local_fixed_sources_adjoint()
{
// Configure the domain for adjoint simulation
FlatSourceDomain::adjoint_ = true;
if (settings::run_mode == RunMode::FIXED_SOURCE) {
domain_->set_local_adjoint_sources();
}
}
void RandomRaySimulation::prepare_adjoint_simulation(bool fw_adjoint)
{
reset_timers();
if (mpi::master)
header("ADJOINT FLUX SOLVE", 3);
if (fw_adjoint) {
// Forward simulation has already been run;
// Configure the domain for adjoint simulation and
// re-initialize OpenMC general data structures
FlatSourceDomain::adjoint_ = true;
openmc_simulation_init();
prepare_fw_fixed_sources_adjoint();
} else {
// Initialize adjoint fixed sources
domain_->apply_meshes();
prepare_local_fixed_sources_adjoint();
domain_->count_external_source_regions();
}
// Reset k-eff
domain_->k_eff_ = 1.0;
// Initialize adjoint fixed sources, if present
prepare_fixed_sources_adjoint();
// Transpose scattering matrix
domain_->transpose_scattering_matrix();
@ -328,18 +390,6 @@ void RandomRaySimulation::prepare_adjoint_simulation()
void RandomRaySimulation::simulate()
{
if (!is_first_simulation_) {
if (mpi::master && adjoint_needed_)
openmc::print_adjoint_header();
// Reset the timers and reinitialize the general OpenMC datastructures if
// this is after the first simulation
reset_timers();
// Initialize OpenMC general data structures
openmc_simulation_init();
}
// Begin main simulation timer
simulation::time_total.start();
@ -435,7 +485,7 @@ void RandomRaySimulation::simulate()
// End main simulation timer
simulation::time_total.stop();
// Normalize and save the final flux
// Normalize and save the final forward flux
double source_normalization_factor =
domain_->compute_fixed_source_normalization_factor() /
(settings::n_batches - settings::n_inactive);
@ -451,11 +501,6 @@ void RandomRaySimulation::simulate()
// Output all simulation results
output_simulation_results();
// Toggle that the simulation object has been initialized after the first
// simulation
if (is_first_simulation_)
is_first_simulation_ = false;
}
void RandomRaySimulation::output_simulation_results() const
@ -622,17 +667,6 @@ void RandomRaySimulation::print_results_random_ray(
}
}
void openmc_finalize_random_ray()
{
FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID;
FlatSourceDomain::volume_normalized_flux_tallies_ = false;
FlatSourceDomain::adjoint_ = false;
FlatSourceDomain::mesh_domain_map_.clear();
RandomRay::ray_source_.reset();
RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
}
} // namespace openmc
//==============================================================================
@ -645,12 +679,25 @@ void openmc_run_random_ray()
// Run forward simulation
//////////////////////////////////////////////////////////
if (openmc::mpi::master) {
if (openmc::FlatSourceDomain::adjoint_) {
openmc::FlatSourceDomain::adjoint_ = false;
openmc::print_adjoint_header();
openmc::FlatSourceDomain::adjoint_ = true;
}
// Check if adjoint calculation is needed, and if local adjoint source(s)
// are present. If an adjoint calculation is needed and no sources are
// specified, we will run a forward calculation first to calculate adjoint
// sources for global variance reduction, then perform an adjoint
// calculation later.
bool adjoint_needed = openmc::FlatSourceDomain::adjoint_;
bool fw_adjoint = openmc::model::adjoint_sources.empty() && adjoint_needed;
// If we're going to do an adjoint simulation with forward-weighted adjoint
// sources afterwards, report that this is the initial forward flux solve.
if (!adjoint_needed || fw_adjoint) {
// Configure the domain for forward simulation
openmc::FlatSourceDomain::adjoint_ = false;
if (adjoint_needed && openmc::mpi::master)
openmc::header("FORWARD FLUX SOLVE", 3);
} else {
// Configure domain for adjoint simulation (later)
openmc::FlatSourceDomain::adjoint_ = true;
}
// Initialize OpenMC general data structures
@ -663,21 +710,25 @@ void openmc_run_random_ray()
// Initialize Random Ray Simulation Object
openmc::RandomRaySimulation sim;
// Initialize fixed sources, if present
sim.apply_fixed_sources_and_mesh_domains();
if (!adjoint_needed || fw_adjoint) {
// Initialize fixed sources, if present
sim.apply_fixed_sources_and_mesh_domains();
// Run initial random ray simulation
sim.simulate();
// Execute random ray simulation
sim.simulate();
}
//////////////////////////////////////////////////////////
// Run adjoint simulation (if enabled)
//////////////////////////////////////////////////////////
if (sim.adjoint_needed_) {
// Setup for adjoint simulation
sim.prepare_adjoint_simulation();
// Run adjoint simulation
sim.simulate();
if (!adjoint_needed) {
return;
}
// Setup for adjoint simulation
sim.prepare_adjoint_simulation(fw_adjoint);
// Execute random ray simulation
sim.simulate();
}

View file

@ -280,8 +280,9 @@ void get_run_parameters(pugi::xml_node node_base)
} else {
fatal_error("Specify random ray inactive distance in settings XML");
}
if (check_for_node(random_ray_node, "source")) {
xml_node source_node = random_ray_node.child("source");
if (check_for_node(random_ray_node, "ray_source")) {
xml_node ray_source_node = random_ray_node.child("ray_source");
xml_node source_node = ray_source_node.child("source");
// Get point to list of <source> elements and make sure there is at least
// one
RandomRay::ray_source_ = Source::create(source_node);
@ -369,6 +370,13 @@ void get_run_parameters(pugi::xml_node node_base)
"between 0 and 1");
}
}
if (check_for_node(random_ray_node, "adjoint_source")) {
pugi::xml_node adj_source_node = random_ray_node.child("adjoint_source");
for (pugi::xml_node source_node : adj_source_node.children("source")) {
// Find any local adjoint sources
model::adjoint_sources.push_back(Source::create(source_node));
}
}
}
}
@ -1276,6 +1284,16 @@ void read_settings_xml(pugi::xml_node root)
break;
}
}
// If any weight window generators have local FW-CADIS target tallies,
// user-defined adjoint sources cannot be used at the same time.
if (!model::adjoint_sources.empty()) {
for (const auto& wwg : variance_reduction::weight_windows_generators) {
if (!wwg->targets_.empty()) {
fatal_error("Cannot use both user-defined adjoint sources and "
"FW-CADIS target tallies at the same time.");
}
}
}
}
// Set up weight window checkpoints

View file

@ -62,6 +62,8 @@ namespace model {
vector<unique_ptr<Source>> external_sources;
vector<unique_ptr<Source>> adjoint_sources;
DiscreteIndex external_sources_probability;
} // namespace model
@ -352,6 +354,19 @@ IndependentSource::IndependentSource(pugi::xml_node node) : Source(node)
if (check_for_node(node, "energy")) {
pugi::xml_node node_dist = node.child("energy");
energy_ = distribution_from_xml(node_dist);
// For decay photon sources, use the absolute photon emission rate in
// [photons/s] as the source strength
if (dynamic_cast<DecaySpectrum*>(energy_.get())) {
if (strength_ != 1.0) {
warning(fmt::format(
"Source strength of {} is ignored because the source uses a "
"DecaySpectrum energy distribution. The source strength will be "
"set from the DecaySpectrum emission rate.",
strength_));
}
strength_ = energy_->integral();
}
} else {
// Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
energy_ = UPtrDist {new Watt(0.988e6, 2.249e-6)};
@ -412,6 +427,7 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
// Check for monoenergetic source above maximum particle energy
auto p = particle_.transport_index();
auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
auto decay_spectrum = dynamic_cast<DecaySpectrum*>(energy_.get());
if (energy_ptr) {
auto energies =
tensor::Tensor<double>(energy_ptr->x().data(), energy_ptr->x().size());
@ -422,10 +438,18 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
}
while (true) {
// Sample energy spectrum
auto [E, E_wgt_temp] = energy_->sample(seed);
site.E = E;
E_wgt = E_wgt_temp;
// Sample energy spectrum. For decay photon sources, also get the parent
// nuclide index to store in the source site for tallying purposes.
if (decay_spectrum) {
auto sample = decay_spectrum->sample_with_parent(seed);
site.E = sample.energy;
E_wgt = sample.weight;
site.parent_nuclide = sample.parent_nuclide;
} else {
auto [E, E_wgt_temp] = energy_->sample(seed);
site.E = E;
E_wgt = E_wgt_temp;
}
// Resample if energy falls above maximum particle energy
if (site.E < data::energy_max[p] &&
@ -710,6 +734,7 @@ SourceSite sample_external_source(uint64_t* seed)
void free_memory_source()
{
model::external_sources.clear();
model::adjoint_sources.clear();
reset_source_rejection_counters();
}

View file

@ -592,8 +592,16 @@ void write_source_point(std::string filename, span<SourceSite> source_bank,
const vector<int64_t>& bank_index, bool use_mcpl)
{
std::string ext = use_mcpl ? "mcpl" : "h5";
int total_surf_particles = source_bank.size();
#ifdef OPENMC_MPI
int num_particles = source_bank.size();
MPI_Allreduce(
&num_particles, &total_surf_particles, 1, MPI_INT, MPI_SUM, mpi::intracomm);
#endif
write_message("Creating source file {}.{} with {} particles ...", filename,
ext, source_bank.size(), 5);
ext, total_surf_particles, 5);
// Dispatch to appropriate function based on file type
if (use_mcpl) {

View file

@ -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;
}
}
@ -2656,19 +2654,19 @@ void score_meshsurface_tally(Particle& p, const vector<int>& tallies)
}
void score_surface_tally(
Particle& p, const vector<int>& tallies, const Surface& surf)
Particle& p, const vector<int>& tallies, const Direction& normal)
{
double wgt = p.wgt_last();
double mu = std::clamp(p.u().dot(normal), -1.0, 1.0);
// Sign for net current: +1 if crossing outward (in direction of normal),
// -1 if crossing inward
double current_sign = (p.surface() > 0) ? 1.0 : -1.0;
double current_sign = std::copysign(1.0, mu);
// Determine absolute cosine of angle between particle direction and surface
// normal, needed for the surface-crossing flux estimator.
auto n = surf.normal(p.r());
n /= n.norm();
double abs_mu = std::min(std::abs(p.u().dot(n)), 1.0);
double abs_mu = std::abs(mu);
if (abs_mu < settings::surface_grazing_cutoff)
abs_mu = settings::surface_grazing_ratio * settings::surface_grazing_cutoff;

View file

@ -621,7 +621,7 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value,
}
}
} else {
// For FW-CADIS, weight windows are inversely proportional to the adjoint
// For (FW-)CADIS, weight windows are inversely proportional to the adjoint
// fluxes. We normalize the weight windows across all energy groups.
#pragma omp parallel for collapse(2) schedule(static)
for (int e = 0; e < e_bins; e++) {
@ -801,6 +801,13 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node)
fatal_error("FW-CADIS can only be run in random ray solver mode.");
}
FlatSourceDomain::adjoint_ = true;
if (check_for_node(node, "targets")) {
FlatSourceDomain::fw_cadis_local_ = true;
targets_ = get_node_array<size_t>(node, "targets");
FlatSourceDomain::fw_cadis_local_targets_.insert(
std::end(FlatSourceDomain::fw_cadis_local_targets_),
std::begin(targets_), std::end(targets_));
}
} else {
fatal_error(fmt::format(
"Unknown weight window update method '{}' specified", method_string));

View file

@ -24,7 +24,7 @@ DepletionSolutionTuple = namedtuple(
predictor_solution = DepletionSolutionTuple(
PredictorIntegrator, np.array([1.0, 2.46847546272295, 4.11525874568034]),
np.array([1.0, 0.986431226850467, -0.0581692232513460]))
np.array([1.0, 0.986431226850467, 0.0]))
cecm_solution = DepletionSolutionTuple(

View file

@ -0,0 +1,140 @@
""" Tests for KeffSearchControl class """
from hmac import new
from pathlib import Path
import shutil
import sys
import pytest
import numpy as np
import openmc
import openmc.lib
from openmc.deplete import CoupledOperator
from tests.regression_tests import config
@pytest.fixture
def model():
f = openmc.Material(name='f')
f.set_density('g/cm3', 10.29769)
f.add_element('U', 1., enrichment=2.4)
f.add_element('O', 2.)
h = openmc.Material(name='h')
h.set_density('g/cm3', 0.001598)
h.add_element('He', 2.4044e-4)
w = openmc.Material(name='w')
w.set_density('g/cm3', 0.740582)
w.add_element('H', 2)
w.add_element('O', 1)
# Define overall material
materials = openmc.Materials([f, h, w])
# Define surfaces
radii = [0.5, 0.8, 1]
height = 80
surf_in = openmc.ZCylinder(r=radii[0])
surf_mid = openmc.ZCylinder(r=radii[1])
surf_out = openmc.ZCylinder(r=radii[2], boundary_type='reflective')
surf_top = openmc.ZPlane(z0=height/2, boundary_type='vacuum')
surf_bot = openmc.ZPlane(z0=-height/2, boundary_type='vacuum')
surf_trans = openmc.ZPlane(z0=0)
surf_rot1 = openmc.XPlane(x0=0)
surf_rot2 = openmc.YPlane(y0=0)
# Define cells
cell_f = openmc.Cell(name='fuel_cell', fill=f,
region=-surf_in & -surf_top & +surf_bot)
cell_g = openmc.Cell(fill=h,
region = +surf_in & -surf_mid & -surf_top & +surf_bot & +surf_rot2)
# Define unbounded cells for rotation universe
cell_w = openmc.Cell(fill=w, region = -surf_rot1)
cell_h = openmc.Cell(fill=h, region = +surf_rot1)
universe_rot = openmc.Universe(cells=(cell_w, cell_h))
cell_rot = openmc.Cell(name="rot_cell", fill=universe_rot,
region = +surf_in & -surf_mid & -surf_top & +surf_bot & -surf_rot2)
# Define unbounded cells for translation universe
cell_w = openmc.Cell(fill=w, region=+surf_in & -surf_trans )
cell_h = openmc.Cell(fill=h, region=+surf_in & +surf_trans)
universe_trans = openmc.Universe(cells=(cell_w, cell_h))
cell_trans = openmc.Cell(name="trans_cell", fill=universe_trans,
region=+surf_mid & -surf_out & -surf_top & +surf_bot)
# Define overall geometry
geometry = openmc.Geometry([cell_f, cell_g, cell_rot, cell_trans])
# Set material volume for depletion fuel.
f.volume = np.pi * radii[0]**2 * height
settings = openmc.Settings()
settings.particles = 1000
settings.inactive = 10
settings.batches = 50
return openmc.Model(geometry, materials, settings)
def translate_cell(position):
cell_trans = [cell for cell in openmc.lib.cells.values() if cell.name == 'trans_cell'][0]
openmc.lib.cells[cell_trans.id].translation = [0, 0, position]
def rotate_cell(angle):
cell_rot = [cell for cell in openmc.lib.cells.values() if cell.name == 'rot_cell'][0]
openmc.lib.cells[cell_rot.id].rotation = [0, 0, angle]
def set_u235_density(u235_density):
fuel = [material for material in openmc.lib.materials.values()
if material.name == 'f'][0]
nuclides = openmc.lib.materials[fuel.id].nuclides
densities = openmc.lib.materials[fuel.id].densities
nuc_idx = nuclides.index('U235')
densities[nuc_idx] = u235_density
openmc.lib.materials[fuel.id].set_densities(nuclides, densities)
@pytest.mark.parametrize("function, x0, x1, bracket, ref_result", [
(translate_cell, -11, -5, (-15, 0), 'depletion_with_translation'),
(rotate_cell, -80, -50, (-90, 0), 'depletion_with_rotation'),
(set_u235_density, 2e-4, 1e-3, (1e-4, 2e-3), 'depletion_with_refuel')
])
def test_keff_search_control(run_in_tmpdir, model, function, x0, x1, bracket, ref_result):
chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
model.settings.verbosity = 1
op = CoupledOperator(model, chain_file)
integrator = openmc.deplete.PredictorIntegrator(
op, [1], 174., timestep_units = 'd')
integrator.add_keff_search_control(
function=function,
x0=x0,
x1=x1,
bracket=bracket,
output=True,
k_tol=0.1,
sigma_final=5e-2)
integrator.integrate()
# Get path to test and reference results
path_test = op.output_dir / 'depletion_results.h5'
path_reference = Path(__file__).with_name(f'ref_{ref_result}.h5')
# If updating results, do so and return
if config['update']:
shutil.copyfile(str(path_test), str(path_reference))
return
# Load the reference/test results
res_test = openmc.deplete.Results(path_test)
res_ref = openmc.deplete.Results(path_reference)
# Use high tolerance here
assert res_test[0].keff_search_root == pytest.approx(res_ref[0].keff_search_root, rel=2)

View file

@ -207,11 +207,13 @@
<random_ray>
<distance_active>500.0</distance_active>
<distance_inactive>100.0</distance_inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
<ray_source>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>0.0 0.0 0.0 30.0 30.0 30.0</parameters>
</space>
</source>
</ray_source>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<adjoint>true</adjoint>
<volume_estimator>naive</volume_estimator>

View file

@ -80,11 +80,13 @@
<random_ray>
<distance_active>100.0</distance_active>
<distance_inactive>20.0</distance_inactive>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<ray_source>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
</ray_source>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<adjoint>true</adjoint>
</random_ray>

View file

@ -0,0 +1,293 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="source">
<density value="1.0" units="macro"/>
<macroscopic name="source"/>
</material>
<material id="2" name="cavity">
<density value="1.0" units="macro"/>
<macroscopic name="cavity"/>
</material>
<material id="3" name="absorber">
<density value="1.0" units="macro"/>
<macroscopic name="absorber"/>
</material>
</materials>
<geometry>
<cell id="1" name="infinite source region" material="1" universe="1"/>
<cell id="2" name="cube cavity region" material="2" universe="2"/>
<cell id="3" name="absorber region" material="3" universe="3"/>
<cell id="4" fill="4" universe="5"/>
<cell id="5" name="full domain" fill="5" region="1 -2 3 -4 5 -6" universe="6"/>
<cell id="6" name="detector 1" material="3" region="1 -10 4 -8 5 -11" universe="6"/>
<cell id="7" name="detector 2" material="3" region="2 -7 4 -8 6 -9" universe="6"/>
<cell id="8" name="outside cube" material="2" region="(2 3 5 -7 ((-8 -6) | (-4 6 -9))) | (4 -8 ((10 -2 5 -9) | (-10 1 11 -9) | (2 -7 5 -6))) | (1 3 6 -9 ((-8 -2) | (-4 2 -7)))" universe="6"/>
<lattice id="4">
<pitch>2.5 2.5 2.5</pitch>
<dimension>12 12 12</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<universes>
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
1 1 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
2 2 2 2 2 2 2 2 2 2 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 </universes>
</lattice>
<surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
<surface id="2" type="x-plane" coeffs="30.0"/>
<surface id="3" type="y-plane" boundary="reflective" coeffs="0.0"/>
<surface id="4" type="y-plane" coeffs="30.0"/>
<surface id="5" type="z-plane" boundary="reflective" coeffs="0.0"/>
<surface id="6" type="z-plane" coeffs="30.0"/>
<surface id="7" type="x-plane" boundary="vacuum" coeffs="35.0"/>
<surface id="8" type="y-plane" boundary="vacuum" coeffs="35.0"/>
<surface id="9" type="z-plane" boundary="vacuum" coeffs="35.0"/>
<surface id="10" type="x-plane" coeffs="5.0"/>
<surface id="11" type="z-plane" coeffs="5.0"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>500</particles>
<batches>10</batches>
<inactive>5</inactive>
<source type="independent" strength="3.14" particle="neutron">
<energy type="discrete">
<parameters>100.0 1.0</parameters>
</energy>
<constraints>
<domain_type>universe</domain_type>
<domain_ids>1</domain_ids>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<random_ray>
<distance_active>800.0</distance_active>
<distance_inactive>100.0</distance_inactive>
<ray_source>
<source type="independent" strength="1.0" particle="neutron">
<space type="box">
<parameters>0.0 0.0 0.0 35.0 35.0 35.0</parameters>
</space>
</source>
</ray_source>
<volume_normalized_flux_tallies>true</volume_normalized_flux_tallies>
<source_region_meshes>
<mesh id="1">
<domain id="6" type="universe"/>
</mesh>
</source_region_meshes>
<adjoint>true</adjoint>
<adjoint_source>
<source type="independent" strength="3.14" particle="neutron">
<energy type="discrete">
<parameters>100.0 1.0</parameters>
</energy>
<constraints>
<domain_type>cell</domain_type>
<domain_ids>6 7</domain_ids>
</constraints>
</source>
</adjoint_source>
<volume_estimator>naive</volume_estimator>
</random_ray>
<mesh id="1">
<dimension>14 14 14</dimension>
<lower_left>0.0 0.0 0.0</lower_left>
<upper_right>35.0 35.0 35.0</upper_right>
</mesh>
</settings>
<tallies>
<filter id="1" type="cell">
<bins>6</bins>
</filter>
<filter id="2" type="cell">
<bins>7</bins>
</filter>
<filter id="3" type="material">
<bins>3</bins>
</filter>
<filter id="4" type="material">
<bins>2</bins>
</filter>
<filter id="5" type="material">
<bins>1</bins>
</filter>
<tally id="1" name="Detector 1 Tally">
<filters>1</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="2" name="Detector 2 Tally">
<filters>2</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="3" name="Absorber Tally">
<filters>3</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="4" name="Cavity Tally">
<filters>4</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
<tally id="5" name="Source Tally">
<filters>5</filters>
<scores>flux</scores>
<estimator>tracklength</estimator>
</tally>
</tallies>
</model>

View file

@ -0,0 +1,15 @@
tally 1:
2.215273E+01
9.815738E+01
tally 2:
1.873933E+01
7.023420E+01
tally 3:
4.802282E-01
4.612707E-02
tally 4:
2.516720E-01
1.271063E-02
tally 5:
1.169938E-02
3.277334E-05

View file

@ -0,0 +1,35 @@
import os
import openmc
from openmc.examples import random_ray_three_region_cube_with_detectors
from tests.testing_harness import TolerantPyAPITestHarness
class MGXSTestHarness(TolerantPyAPITestHarness):
def _cleanup(self):
super()._cleanup()
f = 'mgxs.h5'
if os.path.exists(f):
os.remove(f)
def test_random_ray_adjoint_local():
model = random_ray_three_region_cube_with_detectors()
detector1_cells = model.geometry.get_cells_by_name("detector 1")
detector2_cells = model.geometry.get_cells_by_name("detector 2")
detector_cells = detector1_cells + detector2_cells
strengths = [1.0]
midpoints = [100.0]
energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
adj_source = openmc.IndependentSource(energy=energy_distribution, constraints={
'domains': detector_cells}, strength=3.14)
model.settings.random_ray['adjoint'] = True
model.settings.random_ray['adjoint_source'] = adj_source
model.settings.random_ray['volume_estimator'] = 'naive'
harness = MGXSTestHarness('statepoint.10.h5', model)
harness.main()

Some files were not shown because too many files have changed in this diff Show more