diff --git a/.github/workflows/black_lint.yml b/.github/workflows/black_lint.yml new file mode 100644 index 0000000..f33f91d --- /dev/null +++ b/.github/workflows/black_lint.yml @@ -0,0 +1,27 @@ +name: Black Style Check + +on: + push: + pull_request: + +jobs: + black-lint: + name: Lint with Black + runs-on: ubuntu-latest + + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.14" + + - name: Install Black + run: | + pip install --upgrade pip + pip install black + + - name: Run Black (check mode) + run: black --check . diff --git a/.gitignore b/.gitignore index 51891a7..a3529f8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,19 +1,29 @@ -configs/launch_config.py -configs/user_config.py +# Local configuration +/configs/launch_config.py +/configs/user_config.py -results +# Python caches +__pycache__/ -tallies.out -statepoint* -tmp* -submit.pbs -*png -__pycache__ +# Verification results +/results/ + +# Simulation output output* -*.core -*xml -*DS_Store +statepoint* +tallies.out summary.h5 mgxs.h5 -*gif +*.xml + +# Generated plots +*.png +*.gif + +# Batch and temporary files +submit.pbs +tmp* +*.core +# Operating-system files +.DS_Store diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..2407cc3 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,6 @@ +repos: + - repo: https://github.com/psf/black-pre-commit-mirror + rev: 26.1.0 + hooks: + - id: black + language_version: python3.14 diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..ce1d4ba --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,36 @@ +# Changelog + +All notable changes to this project will be documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/2.0.0/), and this project uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html) as a guide. +MC/DC-VVP release numbers align with the corresponding compatible MC/DC release. + +## [Unreleased] + +Target release: 0.16.0, the first versioned MC/DC-VVP release. + +### Added + +- Add the analytical neutron $k$-eigenvalue suite with subcritical and supercritical SHEM-361 cases, analytical matrix-eigenvalue references, active-cycle convergence studies, and uncertainty plots, from [@ilhamv] +- Add energy-dependent weight-window and time-census variants of the infinite homogeneous SHEM-361 problem, from [@ilhamv] +- Add AZURV1 variants for basic variance-reduction techniques, analytical spatial weight windows, time censuses, and census-based tallies, from [@ilhamv] +- Add top-level and suite-level READMEs describing layouts, configuration, launching, processing, cases, and references, from [@ilhamv] +- Add shared platform, user, and launch configuration for local and HPC campaigns, from [@ilhamv] + +### Changed + +- Migration to Maestro-based launch, from [@ilhamv] +- Update analytical fixed-source cases for the simulation-owned MC/DC interface and unified material model, from [@ilhamv] +- Standardize **suite** and **case** as the VVP repository's organizational terminology, from [@ilhamv] +- Organize fixed-source cases around consistent input, reference, processing, and optional plotting scripts, from [@ilhamv] + +### Deprecated + +### Removed + +### Fixed + +### Security + +[Unreleased]: https://github.com/ilhamv/mcdc-vvp/tree/master +[@ilhamv]: https://github.com/ilhamv diff --git a/README.md b/README.md index bb9f2e2..e673e0a 100644 --- a/README.md +++ b/README.md @@ -6,69 +6,93 @@ A collection of verification, validation, and performance (VVP) test suites for [MC/DC](https://github.com/mcdc-project/mcdc). -The repository provides a unified framework for launching, processing, and organizing MC/DC verification campaigns on local workstations and HPC platforms. Each verification, validation, and performance suite is self-contained and can be executed independently, while the top-level workflow enables reproducible campaign-wide execution. On supported HPC systems, campaign orchestration is performed using [Maestro](https://github.com/llnl/maestrowf). +The repository provides a unified framework for launching, processing, and organizing MC/DC-VVP campaigns on local workstations and HPC platforms. +Each suite is self-contained and can be executed independently, while the top-level workflow enables reproducible campaign-wide execution. +Workflow orchestration is performed using [Maestro](https://github.com/llnl/maestrowf). -## Repository organization +## Directory layout ```text -configs/ Shared platform, user, and launch configurations -verification/ Verification test suites -validation/ Validation test suites -performance/ Performance test suites -results/ Processed results from completed campaigns - -launch.py Launch all enabled suites -process.py Process all enabled suites +configs/ Shared platform, user, and launch configurations +verification/ Verification suites and their cases +results/metadata.yaml Append-only launch and processing history +results// Processed results for one recorded launch + +launch.py Launch all enabled suites +process.py Process one recorded launch ``` -Each suite contains its own launch and processing workflow and may also be executed manually without the top-level scripts. +MC/DC-VVP uses **suite** and **case** as standard terms for its two organizational levels: + +- A **suite** is a self-contained collection of related VVP cases with a shared launch and processing workflow. +- A **case** is one individual problem definition and its inputs, reference solution or data, and processing logic. -## Quick start +The top-level workflow launches and processes enabled suites, while each suite workflow runs and processes its cases. +Every integrated suite provides a README that describes its layout, configuration, workflow, and cases. -1. Create a user launch configuration: +## Configuration + +Create the local launch configuration: ```bash cp configs/launch_config.py.template configs/launch_config.py ``` -2. Edit `configs/launch_config.py` to enable the desired suites and select the target platform. +Edit `configs/launch_config.py` to enable the desired suites and set their platform and launch options. +Use `platform=None` for local execution or a name from `configs/platform_config.py` for HPC execution. + +For HPC execution, also create `configs/user_config.py` from its template and provide the account and optional queue, reservation, and Python paths for the target platform. -3. Launch the enabled suites: +## Launching and processing + +Launch locally enabled suites configured with `platform=None`: ```bash -# Local execution python launch.py +``` + +Launch enabled suites configured for a specific HPC platform: -# HPC execution +```bash python launch.py --platform tuolumne ``` -4. After all jobs have completed, process the results: +The `--platform` option selects suites with a matching configured platform. + +After the latest launch has completed, process all suites submitted by that launch: ```bash python process.py ``` -Processed figures, metadata, and summary results are written to the `results/` directory. +Process a specific recorded launch by passing the launch ID printed by `launch.py` and stored in `results/metadata.yaml`: -## Verification suites +```bash +python process.py 20260817T120000123456Z +``` + +Each launch is processed into its own `results//` directory, which contains a metadata snapshot and the suite result hierarchy. +Reprocessing one launch replaces only that launch's subfolder and does not affect results from other launches. + +## Suites ### Analytical verification -Analytical verification demonstrates the expected statistical convergence of MC/DC by comparing numerical solutions against analytical reference solutions as the number of source particles is increased. +Analytical verification demonstrates the expected statistical convergence of MC/DC by comparing numerical solutions against analytical reference solutions as the sampling effort is increased. | Physics | Suite | Description | | :------ | :---- | :---------- | -| Neutron transport | Analytical fixed-source | Multigroup steady-state and transient fixed-source verification problems, including Reed's problem, AZURV1 variants, and infinite SHEM-361 benchmarks. | +| Neutron transport | [Fixed-source](verification/analytical/neutron/fixed_source/README.md) | Multigroup steady-state and transient fixed-source cases, including Reed's problem, AZURV1 variants, and infinite SHEM-361 benchmarks. | +| Neutron transport | [$k$-eigenvalue](verification/analytical/neutron/k_eigenvalue/README.md) | Infinite homogeneous SHEM-361 criticality cases with analytical eigenvalue and energy-spectrum references. | ### Benchmark verification -Benchmark verification compares MC/DC against established reference Monte Carlo codes on problems without analytical solutions. +Benchmark verification compares MC/DC against established reference Monte Carlo codes for cases without analytical solutions. | Physics | Suite | Description | | :------ | :---- | :---------- | -| Neutron transport (multigroup) | Benchmark multigroup | Time-dependent benchmark problems, including the Kobayashi Dog-Leg and C5G7 transient benchmarks. | -| Neutron transport (continuous energy) | Benchmark continuous energy | Continuous-energy benchmark problems for representative reactor systems. | +| Neutron transport (multigroup) | Benchmark multigroup | Time-dependent benchmark cases, including the Kobayashi Dog-Leg and C5G7 transient benchmarks. | +| Neutron transport (continuous energy) | [Benchmark continuous energy](verification/benchmark/neutron/continuous_energy/README.md) | Continuous-energy benchmark cases for representative reactor systems. | ## Validation @@ -84,11 +108,4 @@ Performance suites evaluate computational performance, scalability, and efficien ## Documentation -Comprehensive user and developer documentation will be available on Read the Docs. It will include: - -- Framework architecture -- Launch and processing workflow -- Platform and user configuration -- Adding new suites -- Adding new verification cases -- Developer guidelines +The top-level and suite READMEs provide the repository-specific documentation for MC/DC-VVP. diff --git a/configs/launch_config.py.template b/configs/launch_config.py.template index 084b3be..808614b 100644 --- a/configs/launch_config.py.template +++ b/configs/launch_config.py.template @@ -6,4 +6,11 @@ LAUNCH_CONFIG = { "walltime": None, "rewrite": False, }, + "verification/analytical/neutron/k_eigenvalue": { + "enabled": False, + "platform": None, + "mpi": False, + "walltime": None, + "rewrite": False, + }, } diff --git a/configs/user_config.py b/configs/user_config.py deleted file mode 100644 index 34d9e38..0000000 --- a/configs/user_config.py +++ /dev/null @@ -1,15 +0,0 @@ -USER_CONFIG = { - "dane": { - "account": "orsu", - "queue": "pbatch", - "reservation": None, - "maestro_python": None, - }, - "tuolumne": { - "account": "orsu", - "queue": "pbatch", - "reservation": None, - "maestro_python": "/g/g92/variansyah1/venv/tuolumne/maestrowf/bin/python", - "mcdc_python": "/g/g92/variansyah1/venv/tuolumne/mcdc/bin/python", - }, -} diff --git a/configs/user_config.py.template b/configs/user_config.py.template index 577492d..519d0b2 100644 --- a/configs/user_config.py.template +++ b/configs/user_config.py.template @@ -4,6 +4,6 @@ USER_CONFIG = { "queue": "the_queue_here", "reservation": None, "maestro_python": None, - "mcdc_python": None, + "mcdc_python": None, }, } diff --git a/launch.py b/launch.py index dced030..da1db70 100644 --- a/launch.py +++ b/launch.py @@ -1,3 +1,5 @@ +"""Launch configured VVP suites and record enough context to reproduce each run.""" + import argparse import datetime import importlib.metadata @@ -7,13 +9,13 @@ import yaml - # ====================================================================================== # Bootstrap VVP imports # ====================================================================================== REPO_DIR = Path(__file__).resolve().parent +# Support launching this script from any working directory without installing VVP. if str(REPO_DIR) not in sys.path: sys.path.insert(0, str(REPO_DIR)) @@ -22,14 +24,15 @@ # Load launch configuration # ====================================================================================== -from configs.launch_config import LAUNCH_CONFIG # noqa: E402 - +from configs.launch_config import LAUNCH_CONFIG # ====================================================================================== # Helper functions # ====================================================================================== + def get_git_hash(repo_dir): + """Return the commit that identifies the VVP source used for the launch.""" return subprocess.check_output( ["git", "-C", str(repo_dir), "rev-parse", "HEAD"], stderr=subprocess.DEVNULL, @@ -38,6 +41,7 @@ def get_git_hash(repo_dir): def is_git_dirty(repo_dir): + """Report whether the launch includes changes not captured by the commit hash.""" result = subprocess.run( ["git", "-C", str(repo_dir), "status", "--porcelain"], stdout=subprocess.PIPE, @@ -49,12 +53,19 @@ def is_git_dirty(repo_dir): def get_mcdc_version(): + """Return the installed MC/DC version when package metadata is available.""" try: return importlib.metadata.version("mcdc") except importlib.metadata.PackageNotFoundError: return None +def write_metadata(metadata_file, metadata): + """Persist launch metadata after each submission milestone.""" + with metadata_file.open("w") as f: + yaml.dump(metadata, f, sort_keys=False) + + # ====================================================================================== # Command-line arguments # ====================================================================================== @@ -63,7 +74,7 @@ def get_mcdc_version(): parser.add_argument( "--platform", default=None, - help="Active platform. Use None/omit for local suites.", + help="Select suites configured for PLATFORM; omit for local run.", ) args = parser.parse_args() @@ -87,24 +98,27 @@ def get_mcdc_version(): metadata.setdefault("launches", []) -metadata["launches"].append( - { - "launched_at": datetime.datetime.now(datetime.UTC).isoformat(), - "active_platform": active_platform, - "mcdc_version": get_mcdc_version(), - "mcdc_vvp_hash": get_git_hash(REPO_DIR), - "mcdc_vvp_dirty": is_git_dirty(REPO_DIR), - "launch_config": LAUNCH_CONFIG, - } -) - -with metadata_file.open("w") as f: - yaml.dump(metadata, f, sort_keys=False) +# Keep an append-only launch history so reruns do not erase provenance. +launch_time = datetime.datetime.now(datetime.UTC) +launch_record = { + "launch_id": launch_time.strftime("%Y%m%dT%H%M%S%fZ"), + "launched_at": launch_time.isoformat(), + "active_platform": active_platform, + "mcdc_version": get_mcdc_version(), + "mcdc_vvp_hash": get_git_hash(REPO_DIR), + "mcdc_vvp_dirty": is_git_dirty(REPO_DIR), + "launch_config": LAUNCH_CONFIG, + "suite_runs": {}, +} +metadata["launches"].append(launch_record) + +write_metadata(metadata_file, metadata) print("=" * 80) print("Prepared VVP results metadata") print(f"Results directory : {results_dir}") print(f"Metadata file : {metadata_file}") +print(f"Launch ID : {launch_record['launch_id']}") print(f"Active platform : {active_platform}") print("=" * 80) @@ -120,6 +134,7 @@ def get_mcdc_version(): suite_platform = options.get("platform") + # Launch only suites explicitly assigned to this invocation's platform. if suite_platform != active_platform: print( f"Skip platform mismatch: {suite} " @@ -136,11 +151,13 @@ def get_mcdc_version(): if not launcher.is_file(): raise FileNotFoundError(f"Suite launcher not found: {launcher}") + # Reuse the active interpreter so suite launchers inherit this environment. command = [ sys.executable, str(launcher), ] + # Forward only options that are meaningful for the configured suite. if suite_platform is not None: command.extend(["--platform", suite_platform]) @@ -158,12 +175,30 @@ def get_mcdc_version(): print("Command:", " ".join(command)) print("=" * 80) + maestro_runs_before = set(suite_dir.glob("maestro_run_*")) subprocess.run(command, cwd=suite_dir, check=True) + # Associate this campaign with the exact generated suite launch. + maestro_runs_after = set(suite_dir.glob("maestro_run_*")) + new_maestro_runs = maestro_runs_after - maestro_runs_before + + if not new_maestro_runs: + raise RuntimeError(f"Suite did not create a new Maestro run: {suite}") + + maestro_run = max(new_maestro_runs, key=lambda path: path.stat().st_mtime) + launch_record["suite_runs"][suite] = str(maestro_run.relative_to(REPO_DIR)) + write_metadata(metadata_file, metadata) + # ====================================================================================== # Summary # ====================================================================================== +launch_record["submission_completed_at"] = datetime.datetime.now( + datetime.UTC +).isoformat() +write_metadata(metadata_file, metadata) + print() +print(f"Launch ID: {launch_record['launch_id']}") print("Launch complete.") diff --git a/process.py b/process.py index 91d2516..8ebd68b 100644 --- a/process.py +++ b/process.py @@ -1,47 +1,117 @@ +"""Process one recorded VVP launch into an isolated results directory.""" + +import argparse +import datetime import shutil import subprocess import sys from pathlib import Path +import yaml + +REPO_DIR = Path(__file__).resolve().parent + # ====================================================================================== -# Bootstrap VVP imports +# Helper functions # ====================================================================================== -REPO_DIR = Path(__file__).resolve().parent -if str(REPO_DIR) not in sys.path: - sys.path.insert(0, str(REPO_DIR)) +def get_launch_id(launch): + """Return the stored launch ID or derive one for legacy metadata.""" + if "launch_id" in launch: + return launch["launch_id"] + + launched_at = datetime.datetime.fromisoformat(launch["launched_at"]) + return launched_at.astimezone(datetime.UTC).strftime("%Y%m%dT%H%M%SZ") + + +def write_metadata(metadata_file, metadata): + """Persist updated processing provenance.""" + with metadata_file.open("w") as f: + yaml.dump(metadata, f, sort_keys=False) # ====================================================================================== -# Load launch configuration +# Command-line arguments # ====================================================================================== -from configs.launch_config import LAUNCH_CONFIG # noqa: E402 - +parser = argparse.ArgumentParser(description="Process a recorded MC/DC VVP launch.") +parser.add_argument( + "launch_id", + nargs="?", + default=None, + help="Metadata launch ID to process. Defaults to the latest launch.", +) +args = parser.parse_args() # ====================================================================================== # Set up results # ====================================================================================== results_dir = REPO_DIR / "results" +metadata_file = results_dir / "metadata.yaml" -if not results_dir.is_dir(): +if not metadata_file.is_file(): raise FileNotFoundError( - "Top-level results directory not found. Run launch.py before process.py." + "Top-level launch metadata not found. Run launch.py before process.py." + ) + +with metadata_file.open("r") as f: + metadata = yaml.safe_load(f) or {} + +launches = metadata.get("launches", []) +if not launches: + raise ValueError(f"No launches are recorded in {metadata_file}") + +if args.launch_id is None: + launch = launches[-1] +else: + matches = [item for item in launches if get_launch_id(item) == args.launch_id] + if not matches: + available = ", ".join(get_launch_id(item) for item in launches) + raise ValueError( + f"Launch ID '{args.launch_id}' was not found. Available launches: {available}" + ) + launch = matches[-1] + +launch_id = get_launch_id(launch) +launch_config = launch["launch_config"] +suite_runs = launch.get("suite_runs", {}) +has_recorded_suite_runs = "suite_runs" in launch +launch_results = results_dir / launch_id +launch_results.mkdir(parents=True, exist_ok=True) + +if not has_recorded_suite_runs: + print( + "Warning: this legacy launch does not record exact suite runs; " + "each suite will process its latest Maestro run." ) +print("=" * 80) +print(f"Processing launch : {launch_id}") +print(f"Launched at : {launch['launched_at']}") +print(f"Results directory : {launch_results}") +print("=" * 80) + # ====================================================================================== # Process enabled suites # ====================================================================================== -for suite, options in LAUNCH_CONFIG.items(): +for suite, options in launch_config.items(): if not options.get("enabled", False): print(f"Skip disabled suite: {suite}") continue + if options.get("platform") != launch.get("active_platform"): + print(f"Skip suite not submitted by this launch: {suite}") + continue + + if has_recorded_suite_runs and suite not in suite_runs: + print(f"Skip suite without a recorded Maestro run: {suite}") + continue + suite_dir = REPO_DIR / suite processor = suite_dir / "process.py" @@ -55,19 +125,26 @@ print(f"Processing suite: {suite}") print("=" * 80) - subprocess.run( - [sys.executable, str(processor)], - cwd=suite_dir, - check=True, - ) - suite_results = suite_dir / "results" + # Prevent figures from a previous suite-processing run leaking into this launch. + if suite_results.exists(): + shutil.rmtree(suite_results) + + command = [sys.executable, str(processor)] + if suite in suite_runs: + maestro_run = REPO_DIR / suite_runs[suite] + if not maestro_run.is_dir(): + raise FileNotFoundError(f"Recorded Maestro run not found: {maestro_run}") + command.append(str(maestro_run)) + + subprocess.run(command, cwd=suite_dir, check=True) + if not suite_results.is_dir(): print(f"No suite results found: {suite_results}") continue - destination = results_dir / suite + destination = launch_results / suite destination.parent.mkdir(parents=True, exist_ok=True) if destination.exists(): @@ -76,10 +153,24 @@ shutil.move(str(suite_results), str(destination)) +# ====================================================================================== +# Store processing metadata +# ====================================================================================== + +launch["launch_id"] = launch_id +launch["processed_at"] = datetime.datetime.now(datetime.UTC).isoformat() +launch["results_directory"] = str(launch_results.relative_to(REPO_DIR)) +write_metadata(metadata_file, metadata) + +with (launch_results / "metadata.yaml").open("w") as f: + yaml.dump(launch, f, sort_keys=False) + + # ====================================================================================== # Summary # ====================================================================================== print() -print(f"Results : {results_dir}") +print(f"Launch : {launch_id}") +print(f"Results : {launch_results}") print("Processing complete.") diff --git a/verification/analytical/neutron/fixed_source/README.md b/verification/analytical/neutron/fixed_source/README.md index 75033a5..afb4e99 100644 --- a/verification/analytical/neutron/fixed_source/README.md +++ b/verification/analytical/neutron/fixed_source/README.md @@ -1,61 +1,94 @@ # Analytical Neutron Fixed-Source Verification -This suite provides the foundation for verifying the neutron transport capabilities of MC/DC using multigroup fixed-source problems with analytical reference solutions. Beyond verifying the multigroup fixed-source mode itself, the suite exercises the majority of the transport algorithms shared across other physics modes, including geometry tracking, surface crossings, material lookup, source sampling, particle transport, and tallying. Verification is performed by observing the expected (N^{-1/2}) statistical convergence of the error as the number of source particles increases. +This suite verifies MC/DC neutron transport using multigroup fixed-source cases with analytical reference solutions. +It exercises geometry tracking, surface crossings, material lookup, source sampling, particle transport, and tallying. +Verification checks for the expected $N^{-1/2}$ statistical convergence as the number of source particles increases. The suite can be executed independently or as part of the top-level MC/DC-VVP workflow. -## Cases +## Directory layout -| Case | Description | -| :--- | :---------- | -| `slab_absorbium` | Steady-state flux distribution in a purely absorbing multilayer slab. | -| `reed` | Reed's classic one-dimensional transport benchmark. | -| `slab_isobeam_td` | Time-dependent flux propagation from an isotropic planar source. | -| `azurv1` | AZURV1 transient benchmark. | -| `azurv1_sub` | Subcritical variant of AZURV1. | -| `azurv1_super` | Supercritical variant of AZURV1. | -| `azurv1-census` | AZURV1 with census-based time discretization. | -| `azurv1-census-tally` | AZURV1 using census tallying. | -| `inf_shem361` | Infinite homogeneous SHEM-361 multigroup spectrum evolution. | -| `inf_shem361_td` | Time-dependent SHEM-361 spectrum evolution. | -| `inf_shem361_td-census` | Time-dependent SHEM-361 spectrum evolution using census discretization. | - -## Running the suite - -The suite may be launched directly: +```text +cases/ Verification case definitions and processing scripts +data/ Shared multigroup cross-section data +maestro_run_*/ Generated Maestro workflow directories +results/ Generated figures from processed cases + +task.yaml Configure the particle-count study for each case +study.yaml Generated Maestro study definition + +launch.py Build and launch the Maestro study +run_case.py Run one case over its particle-count study +process.py Process a completed Maestro study + +cleanup.py Remove generated outputs and figures +util.py Provide shared processing and plotting utilities +``` + +## Configuration + +The `task.yaml` file selects the cases and defines the particle-count range and number of tasks for each case. +Edit this file to change the study without modifying the launch or processing scripts. + +HPC runs use the shared platform settings in the repository's `configs/platform_config.py` and the user-specific settings in `configs/user_config.py`. + +## Launching and processing + +From this suite directory, launch the study locally: ```bash python launch.py ``` -or on a supported HPC platform: +Launch the study on a supported HPC platform: ```bash python launch.py --platform tuolumne --mpi ``` -After all jobs have completed: +Use `--walltime HOURS` to limit the requested walltime and `--rewrite` to replace existing case output. + +After all jobs have completed, process the latest Maestro run: ```bash python process.py ``` -## Directory layout +Pass a Maestro run directory to process a specific run: -```text -cases/ Verification cases -launch.py Launch the suite -process.py Process the suite results -run_case.py Execute a single verification case -cleanup.py Remove generated outputs and processed figures -task.yaml Particle-count study definition -util.py Shared processing and plotting utilities -maestro_run_*/ Maestro workflow outputs -results/ Processed figures generated by this suite +```bash +python process.py maestro_run_ ``` +Processed figures are written to this suite's `results/` directory. +The top-level `process.py` collects these figures under the repository's `results/` directory. + +## Cases + +| Case | Description | +| :--- | :---------- | +| [`slab_absorbium`](cases/slab_absorbium/) | Steady-state flux distribution in a purely absorbing multilayer slab. | +| [`slab_isobeam_td`](cases/slab_isobeam_td/) | Time-dependent flux propagation from an isotropic planar source. | +| [`reed`](cases/reed/) | Reed's classic one-dimensional transport benchmark. | +| [`azurv1`](cases/azurv1/) | AZURV1 transient benchmark. | +| [`azurv1_sub`](cases/azurv1_sub/) | Subcritical variant of AZURV1. | +| [`azurv1_super`](cases/azurv1_super/) | Supercritical variant of AZURV1. | +| [`azurv1-census`](cases/azurv1-census/) | AZURV1 with time censuses. | +| [`azurv1-census-tally`](cases/azurv1-census-tally/) | AZURV1 using time-census-based tallying. | +| [`azurv1-basic_techniques`](cases/azurv1-basic_techniques/) | AZURV1 with implicit capture, weighted emission, and global weight roulette. | +| [`azurv1-weight_windows`](cases/azurv1-weight_windows/) | AZURV1 with spatial weight windows derived from the time-averaged analytical solution. | +| [`inf_shem361`](cases/inf_shem361/) | Infinite homogeneous SHEM-361 steady-state spectrum. | +| [`inf_shem361_td`](cases/inf_shem361_td/) | Time-dependent SHEM-361 spectrum evolution. | +| [`inf_shem361_td-census`](cases/inf_shem361_td-census/) | Time-dependent SHEM-361 spectrum evolution with time censuses. | +| [`inf_shem361-weight_windows`](cases/inf_shem361-weight_windows/) | Infinite homogeneous SHEM-361 with energy-dependent weight windows derived from the analytical spectrum. | + +### SHEM-361 multigroup data + +The SHEM-361 cases use `data/SHEM-361.npz` as their multigroup cross-section dataset. +The dataset was generated by homogenizing a continuous-energy calculation of an infinite lattice of borated PWR fuel pin cells onto the SHEM-361 energy-group structure. + ## References - W. H. Reed, *New Difference Schemes for the Neutron Transport Equation*, Nuclear Science and Engineering, 1971. -- AZURV1 benchmark. -- SHEM-361 multigroup library. +- B. D. Ganapol, R. S. Baker, J. A. Dahl, and R. E. Alcouffe, [*Homogeneous Infinite Media Time-Dependent Analytical Benchmarks*](https://www.osti.gov/biblio/975281), Los Alamos National Laboratory, LA-UR-01-1854, 2001. +- A. Hébert and A. Santamarina, *Refinement of the Santamarina-Hfaiedh Energy Mesh Between 22.5 eV and 11.4 keV*, PHYSOR 2008, vol. 2, pp. 929–938, 2008. diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/input.py b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/input.py new file mode 100644 index 0000000..0b2eb7f --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/input.py @@ -0,0 +1,64 @@ +import numpy as np +import mcdc + +simulation = mcdc.Simulation("AZURV1 basic techniques") + +# ====================================================================================== +# Set model +# ====================================================================================== +# Infinite medium with isotropic plane surface at the center +# Based on Ganapol LA-UR-01-1854 (AZURV1 benchmark) +# Effective scattering ratio c = 1.0 + +# Set materials +m = mcdc.Material.multigroup( + capture=np.array([1.0 / 3.0]), + scatter=np.array([[1.0 / 3.0]]), + fission=np.array([1.0 / 3.0]), + nu_p=np.array([2.0]), +) + +# Set surfaces +s1 = mcdc.Surface.PlaneX(x=-1e10, boundary_condition="reflective") +s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") + +# Set cells +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) + +# ====================================================================================== +# Set source +# ====================================================================================== +# Isotropic pulse at x=t=0 + +source = mcdc.Source( + position=[0.0, 0.0, 0.0], + isotropic=True, + energy=0, + time=0.0, +) +simulation.set_sources([source]) + +# ====================================================================================== +# Set tallies, settings, techniques, and run MC/DC +# ====================================================================================== + +# Tallies +mesh = mcdc.MeshStructured(x=np.linspace(-20.5, 20.5, 202)) +tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +simulation.set_tallies([tally]) + +# Settings +simulation.settings.N_particle = 100 +simulation.settings.N_batch = 2 + +# Techniques +simulation.technique.implicit_capture() +simulation.technique.weighted_emission(weight_target=0.5) +simulation.technique.global_weight_roulette( + weight_threshold=0.75, + weight_target=1.0, +) + +# Run +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/plot.py b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/plot.py new file mode 100644 index 0000000..0e57ebe --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/plot.py @@ -0,0 +1,59 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import h5py +import sys + +# Reference solution +data = np.load("reference.npz") +phi_ref = data["phi"] + +# Get results +output = sys.argv[1] +with h5py.File(output, "r") as f: + x = f["tallies/tracklength_tally_0/grid/x"][:] + dx = x[1:] - x[:-1] + x_mid = 0.5 * (x[:-1] + x[1:]) + t = f["tallies/tracklength_tally_0/grid/time"][:] + dt = t[1:] - t[:-1] + K = len(t) - 1 + + phi = f["tallies/tracklength_tally_0/flux/mean"][:] + phi_sd = f["tallies/tracklength_tally_0/flux/sdev"][:] + + # Normalize + for k in range(K): + phi[k] /= dx * dt[k] + phi_sd[k] /= dx * dt[k] + +# Flux - average +fig = plt.figure() +ax = plt.axes( + xlim=(-21.889999999999997, 21.89), ylim=(-0.042992644459595206, 0.9028455336514992) +) +ax.grid() +ax.set_xlabel(r"$x$") +ax.set_ylabel(r"Flux") +ax.set_title(r"$\bar{\phi}_{k,j}$") +(line1,) = ax.plot([], [], "-b", label="MC") +(line2,) = ax.plot([], [], "--r", label="Ref.") +fb = ax.fill_between([], [], [], [], alpha=0.2, color="b") +text = ax.text(0.02, 0.9, "", transform=ax.transAxes) +ax.legend() + + +def animate(k): + global fb + fb.remove() + line1.set_data(x_mid, phi[k, :]) + fb = ax.fill_between( + x_mid, phi[k, :] - phi_sd[k, :], phi[k, :] + phi_sd[k, :], alpha=0.2, color="b" + ) + line2.set_data(x_mid, phi_ref[k, :]) + text.set_text(r"$t \in [%.1f,%.1f]$ s" % (t[k], t[k + 1])) + return line1, line2, text + + +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) +plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/process.py b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/process.py new file mode 100644 index 0000000..a7d6d03 --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/process.py @@ -0,0 +1,46 @@ +import numpy as np +import h5py +import sys + +# Get util +sys.path.append("../../") +import util + +# Particle counts +N_min = int(sys.argv[1]) +N_max = int(sys.argv[2]) +N = int(sys.argv[3]) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) + +# Reference solution +data = np.load("reference.npz") +phi_ref = data["phi"] + +# Error containers +error = np.zeros(len(N_particle_list)) + +error_max = np.zeros(len(N_particle_list)) + +# Calculate error +for i, N_particle in enumerate(N_particle_list): + # Get results + with h5py.File("output_%i.h5" % (int(N_particle)), "r") as f: + x = f["tallies/tracklength_tally_0/grid/x"][:] + dx = x[1:] - x[:-1] + t = f["tallies/tracklength_tally_0/grid/time"][:] + dt = t[1:] - t[:-1] + K = len(t) - 1 + + phi = f["tallies/tracklength_tally_0/flux/mean"][:] + + # Normalize + for k in range(K): + phi[k] /= dx * dt[k] + + # Get error + error[i] = util.error(phi, phi_ref) + + error_max[i] = util.error_max(phi, phi_ref) + +# Plot +util.plot_convergence("flux", N_particle_list, error, error_max) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/reference.npz b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/reference.npz new file mode 100644 index 0000000..f34183c Binary files /dev/null and b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/reference.npz differ diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/reference.py b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/reference.py new file mode 100644 index 0000000..d595570 --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-basic_techniques/reference.py @@ -0,0 +1,65 @@ +import numpy as np +from scipy.integrate import quad + +# ============================================================================= +# Reference solution generator +# ============================================================================= + +# Scattering ratio +c = 1.0 +i = complex(0, 1) + +# Spatial grid +J = 201 +x = np.linspace(-20.5, 20.5, J + 1) + +# Time grid +K = 20 +t = np.linspace(0.0, 20.0, K + 1) + + +def integrand(u, eta, t): + q = (1 + eta) / (1 - eta) + xi = (np.log(q) + i * u) / (eta + i * np.tan(u / 2)) + return ( + 1.0 + / (np.cos(u / 2)) ** 2 + * (xi**2 * np.e ** (c * t / 2 * (1 - eta**2) * xi)).real + ) + + +def phi(x, t): + if t == 0.0 or abs(x) >= t: + return 0.0 + eta = x / t + integral = quad(integrand, 0.0, np.pi, args=(eta, t))[0] + return np.e**-t / 2 / t * (1 + c * t / 4 / np.pi * (1 - eta**2) * integral) + + +def phi_t(t, x): + if t == 0.0 or abs(x) >= t: + return 0.0 + eta = x / t + integral = quad(integrand, 0.0, np.pi, args=(eta, t))[0] + return np.e**-t / 2 / t * (1 + c * t / 4 / np.pi * (1 - eta**2) * integral) + + +def phiX(x, t0, t1): + return quad(phi_t, t0, t1, args=(x,))[0] + + +phi_avg = np.zeros([K, J]) + +for k in range(K): + for j in range(J): + x0 = x[j] + x1 = x[j + 1] + dx = x1 - x0 + t0 = t[k] + t1 = t[k + 1] + dt = t1 - t0 + phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt + +phi_avg = np.nan_to_num(phi_avg) + +np.savez("reference.npz", x=x, t=t, phi=phi_avg) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/input.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/input.py index 7751fcb..22282c2 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/input.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/input.py @@ -1,6 +1,8 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("AZURV1 census tally") + # ====================================================================================== # Set model # ====================================================================================== @@ -9,7 +11,7 @@ # Effective scattering ratio c = 1.0 # Set materials -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=np.array([1.0 / 3.0]), scatter=np.array([[1.0 / 3.0]]), fission=np.array([1.0 / 3.0]), @@ -21,19 +23,21 @@ s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m) +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) # ====================================================================================== # Set source # ====================================================================================== # Isotropic pulse at x=t=0 -mcdc.Source( +source = mcdc.Source( position=[0.0, 0.0, 0.0], isotropic=True, - energy_group=0, + energy=0, time=0.0, ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, techniques, and run MC/DC @@ -41,20 +45,18 @@ # Tallies mesh = mcdc.MeshStructured(x=np.linspace(-20.5, 20.5, 202)) -mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 50 -mcdc.settings.N_batch = 2 -mcdc.settings.census_bank_buffer_ratio = 5.0 -mcdc.settings.source_bank_buffer_ratio = 5.0 -mcdc.settings.set_time_census(np.linspace(0.0, 20.0, 5)[1:], tally_frequency=5) +simulation.settings.N_particle = 50 +simulation.settings.N_batch = 2 +simulation.settings.census_bank_buffer_ratio = 5.0 +simulation.settings.source_bank_buffer_ratio = 5.0 +simulation.settings.set_time_census(np.linspace(0.0, 20.0, 5)[1:], tally_frequency=5) -# Tecniques -mcdc.simulation.population_control() +# Techniques +simulation.technique.population_control() # Run -mcdc.run() - -# Post-processing -mcdc.recombine_tallies() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/plot.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/plot.py index 02ac4e1..0e57ebe 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/plot.py @@ -2,42 +2,34 @@ import matplotlib.pyplot as plt import matplotlib.animation as animation import h5py +import sys # Reference solution data = np.load("reference.npz") phi_ref = data["phi"] -# The grids -x = np.linspace(-20.5, 20.5, 202) -t = np.linspace(0.0, 20.0, 21) -dx = x[1:] - x[:-1] -x_mid = 0.5 * (x[:-1] + x[1:]) -dt = t[1:] - t[:-1] -J = len(dx) -K = len(dt) +# Get results +output = sys.argv[1] +with h5py.File(output, "r") as f: + x = f["tallies/tracklength_tally_0/grid/x"][:] + dx = x[1:] - x[:-1] + x_mid = 0.5 * (x[:-1] + x[1:]) + t = f["tallies/tracklength_tally_0/grid/time"][:] + dt = t[1:] - t[:-1] + K = len(t) - 1 -# Get the solution -phi = np.zeros((K, J)) -N_census = 4 -N_batch = 10 -for i_census in range(N_census): - for i_batch in range(N_batch): - with h5py.File("output-batch_%i-census_%i.h5" % (i_batch, i_census), "r") as f: - phi[5 * i_census : 5 * i_census + 5, :] += f[ - "tallies/tracklength_tally_0/flux/score" - ][:] - phi[5 * i_census : 5 * i_census + 5] /= N_batch + phi = f["tallies/tracklength_tally_0/flux/mean"][:] + phi_sd = f["tallies/tracklength_tally_0/flux/sdev"][:] -# Normalize -for k in range(K): - phi[k] /= dx * dt[k] + # Normalize + for k in range(K): + phi[k] /= dx * dt[k] + phi_sd[k] /= dx * dt[k] # Flux - average fig = plt.figure() ax = plt.axes( - # xlim=(-21.889999999999997, 21.89), ylim=(-0.042992644459595206, 0.9028455336514992) - xlim=(-21.889999999999997, 21.89), - ylim=(1e-7, 0.9028455336514992), + xlim=(-21.889999999999997, 21.89), ylim=(-0.042992644459595206, 0.9028455336514992) ) ax.grid() ax.set_xlabel(r"$x$") @@ -45,24 +37,23 @@ ax.set_title(r"$\bar{\phi}_{k,j}$") (line1,) = ax.plot([], [], "-b", label="MC") (line2,) = ax.plot([], [], "--r", label="Ref.") -# fb = ax.fill_between([], [], [], [], alpha=0.2, color="b") +fb = ax.fill_between([], [], [], [], alpha=0.2, color="b") text = ax.text(0.02, 0.9, "", transform=ax.transAxes) ax.legend() -ax.set_yscale("log") def animate(k): - # global fb - # fb.remove() + global fb + fb.remove() line1.set_data(x_mid, phi[k, :]) - # fb = ax.fill_between( - # x_mid, phi[k, :] - phi_sd[k, :], phi[k, :] + phi_sd[k, :], alpha=0.2, color="b" - # ) + fb = ax.fill_between( + x_mid, phi[k, :] - phi_sd[k, :], phi[k, :] + phi_sd[k, :], alpha=0.2, color="b" + ) line2.set_data(x_mid, phi_ref[k, :]) text.set_text(r"$t \in [%.1f,%.1f]$ s" % (t[k], t[k + 1])) return line1, line2, text -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/process.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/process.py index 4ffa4f9..a7d6d03 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/process.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/process.py @@ -6,11 +6,11 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution data = np.load("reference.npz") diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/reference.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/reference.py index 261ad68..d595570 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census-tally/reference.py @@ -45,7 +45,7 @@ def phi_t(t, x): def phiX(x, t0, t1): - return quad(phi_t, t0, t1, args=(x))[0] + return quad(phi_t, t0, t1, args=(x,))[0] phi_avg = np.zeros([K, J]) @@ -60,12 +60,6 @@ def phiX(x, t0, t1): dt = t1 - t0 phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt -for j in range(J + 1): - for k in range(K): - t0 = t[k] - t1 = t[k + 1] - dt = t1 - t0 - phi_avg = np.nan_to_num(phi_avg) np.savez("reference.npz", x=x, t=t, phi=phi_avg) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census/input.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census/input.py index a58698c..23bf65f 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census/input.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census/input.py @@ -1,6 +1,8 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("AZURV1 census") + # ====================================================================================== # Set model # ====================================================================================== @@ -9,7 +11,7 @@ # Effective scattering ratio c = 1.0 # Set materials -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=np.array([1.0 / 3.0]), scatter=np.array([[1.0 / 3.0]]), fission=np.array([1.0 / 3.0]), @@ -21,19 +23,21 @@ s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m) +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) # ====================================================================================== # Set source # ====================================================================================== # Isotropic pulse at x=t=0 -mcdc.Source( +source = mcdc.Source( position=[0.0, 0.0, 0.0], isotropic=True, - energy_group=0, + energy=0, time=0.0, ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, techniques, and run MC/DC @@ -41,17 +45,18 @@ # Tallies mesh = mcdc.MeshStructured(x=np.linspace(-20.5, 20.5, 202)) -mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 50 -mcdc.settings.N_batch = 2 -mcdc.settings.census_bank_buffer_ratio = 5.0 -mcdc.settings.source_bank_buffer_ratio = 5.0 -mcdc.settings.set_time_census(np.linspace(0.0, 20.0, 21)[1:-1]) +simulation.settings.N_particle = 50 +simulation.settings.N_batch = 2 +simulation.settings.census_bank_buffer_ratio = 5.0 +simulation.settings.source_bank_buffer_ratio = 5.0 +simulation.settings.set_time_census(np.linspace(0.0, 20.0, 21)[1:-1]) -# Tecniques -mcdc.simulation.population_control() +# Techniques +simulation.technique.population_control() # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census/plot.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census/plot.py index 143aaaf..0e57ebe 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census/plot.py @@ -54,6 +54,6 @@ def animate(k): return line1, line2, text -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census/process.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census/process.py index 4ffa4f9..a7d6d03 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census/process.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census/process.py @@ -6,11 +6,11 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution data = np.load("reference.npz") diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-census/reference.py b/verification/analytical/neutron/fixed_source/cases/azurv1-census/reference.py index 261ad68..d595570 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1-census/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-census/reference.py @@ -45,7 +45,7 @@ def phi_t(t, x): def phiX(x, t0, t1): - return quad(phi_t, t0, t1, args=(x))[0] + return quad(phi_t, t0, t1, args=(x,))[0] phi_avg = np.zeros([K, J]) @@ -60,12 +60,6 @@ def phiX(x, t0, t1): dt = t1 - t0 phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt -for j in range(J + 1): - for k in range(K): - t0 = t[k] - t1 = t[k + 1] - dt = t1 - t0 - phi_avg = np.nan_to_num(phi_avg) np.savez("reference.npz", x=x, t=t, phi=phi_avg) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/input.py b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/input.py new file mode 100644 index 0000000..96e654a --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/input.py @@ -0,0 +1,83 @@ +from pathlib import Path + +import numpy as np +import mcdc + +REFERENCE_DATA = Path(__file__).resolve().parent / "reference.npz" +WINDOW_FLOOR = 1.0e-3 +WINDOW_WIDTH = 2.5 + +simulation = mcdc.Simulation("AZURV1 weight windows") + +# ====================================================================================== +# Set model +# ====================================================================================== +# Infinite medium with isotropic plane surface at the center +# Based on Ganapol LA-UR-01-1854 (AZURV1 benchmark) +# Effective scattering ratio c = 1.0 + +# Set materials +m = mcdc.Material.multigroup( + capture=np.array([1.0 / 3.0]), + scatter=np.array([[1.0 / 3.0]]), + fission=np.array([1.0 / 3.0]), + nu_p=np.array([2.0]), +) + +# Set surfaces +s1 = mcdc.Surface.PlaneX(x=-1e10, boundary_condition="reflective") +s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") + +# Set cells +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) + +# ====================================================================================== +# Set source +# ====================================================================================== +# Isotropic pulse at x=t=0 + +source = mcdc.Source( + position=[0.0, 0.0, 0.0], + isotropic=True, + energy=0, + time=0.0, +) +simulation.set_sources([source]) + +# ====================================================================================== +# Set tallies, settings, techniques, and run MC/DC +# ====================================================================================== + +# Load the analytical grid and time-dependent cell-average flux. +with np.load(REFERENCE_DATA) as data: + x = data["x"] + time = data["t"] + phi = data["phi"] + +# Tallies +mesh = mcdc.MeshStructured(x=x) +tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=time) +simulation.set_tallies([tally]) + +# Settings +simulation.settings.N_particle = 100 +simulation.settings.N_batch = 2 +simulation.settings.active_bank_buffer = 10_000 + +# Average the analytical flux over the complete simulation time. +phi_average = np.average(phi, axis=0, weights=np.diff(time)) + +# Convert the normalized analytical profile into lower, target, and upper weights. +target_weight = phi_average / np.max(phi_average) +target_weight = WINDOW_FLOOR + (1.0 - WINDOW_FLOOR) * target_weight +weight_windows = np.empty((1, len(x) - 1, 1, 1, 3)) +weight_windows[0, :, 0, 0, 0] = target_weight / WINDOW_WIDTH +weight_windows[0, :, 0, 0, 1] = target_weight +weight_windows[0, :, 0, 0, 2] = target_weight * WINDOW_WIDTH + +# Techniques +simulation.technique.weight_windows(weight_windows, mesh=mesh) + +# Run +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/plot.py b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/plot.py new file mode 100644 index 0000000..0e57ebe --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/plot.py @@ -0,0 +1,59 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.animation as animation +import h5py +import sys + +# Reference solution +data = np.load("reference.npz") +phi_ref = data["phi"] + +# Get results +output = sys.argv[1] +with h5py.File(output, "r") as f: + x = f["tallies/tracklength_tally_0/grid/x"][:] + dx = x[1:] - x[:-1] + x_mid = 0.5 * (x[:-1] + x[1:]) + t = f["tallies/tracklength_tally_0/grid/time"][:] + dt = t[1:] - t[:-1] + K = len(t) - 1 + + phi = f["tallies/tracklength_tally_0/flux/mean"][:] + phi_sd = f["tallies/tracklength_tally_0/flux/sdev"][:] + + # Normalize + for k in range(K): + phi[k] /= dx * dt[k] + phi_sd[k] /= dx * dt[k] + +# Flux - average +fig = plt.figure() +ax = plt.axes( + xlim=(-21.889999999999997, 21.89), ylim=(-0.042992644459595206, 0.9028455336514992) +) +ax.grid() +ax.set_xlabel(r"$x$") +ax.set_ylabel(r"Flux") +ax.set_title(r"$\bar{\phi}_{k,j}$") +(line1,) = ax.plot([], [], "-b", label="MC") +(line2,) = ax.plot([], [], "--r", label="Ref.") +fb = ax.fill_between([], [], [], [], alpha=0.2, color="b") +text = ax.text(0.02, 0.9, "", transform=ax.transAxes) +ax.legend() + + +def animate(k): + global fb + fb.remove() + line1.set_data(x_mid, phi[k, :]) + fb = ax.fill_between( + x_mid, phi[k, :] - phi_sd[k, :], phi[k, :] + phi_sd[k, :], alpha=0.2, color="b" + ) + line2.set_data(x_mid, phi_ref[k, :]) + text.set_text(r"$t \in [%.1f,%.1f]$ s" % (t[k], t[k + 1])) + return line1, line2, text + + +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) +plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/process.py b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/process.py new file mode 100644 index 0000000..a7d6d03 --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/process.py @@ -0,0 +1,46 @@ +import numpy as np +import h5py +import sys + +# Get util +sys.path.append("../../") +import util + +# Particle counts +N_min = int(sys.argv[1]) +N_max = int(sys.argv[2]) +N = int(sys.argv[3]) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) + +# Reference solution +data = np.load("reference.npz") +phi_ref = data["phi"] + +# Error containers +error = np.zeros(len(N_particle_list)) + +error_max = np.zeros(len(N_particle_list)) + +# Calculate error +for i, N_particle in enumerate(N_particle_list): + # Get results + with h5py.File("output_%i.h5" % (int(N_particle)), "r") as f: + x = f["tallies/tracklength_tally_0/grid/x"][:] + dx = x[1:] - x[:-1] + t = f["tallies/tracklength_tally_0/grid/time"][:] + dt = t[1:] - t[:-1] + K = len(t) - 1 + + phi = f["tallies/tracklength_tally_0/flux/mean"][:] + + # Normalize + for k in range(K): + phi[k] /= dx * dt[k] + + # Get error + error[i] = util.error(phi, phi_ref) + + error_max[i] = util.error_max(phi, phi_ref) + +# Plot +util.plot_convergence("flux", N_particle_list, error, error_max) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/reference.npz b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/reference.npz new file mode 100644 index 0000000..f34183c Binary files /dev/null and b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/reference.npz differ diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/reference.py b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/reference.py new file mode 100644 index 0000000..d595570 --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/azurv1-weight_windows/reference.py @@ -0,0 +1,65 @@ +import numpy as np +from scipy.integrate import quad + +# ============================================================================= +# Reference solution generator +# ============================================================================= + +# Scattering ratio +c = 1.0 +i = complex(0, 1) + +# Spatial grid +J = 201 +x = np.linspace(-20.5, 20.5, J + 1) + +# Time grid +K = 20 +t = np.linspace(0.0, 20.0, K + 1) + + +def integrand(u, eta, t): + q = (1 + eta) / (1 - eta) + xi = (np.log(q) + i * u) / (eta + i * np.tan(u / 2)) + return ( + 1.0 + / (np.cos(u / 2)) ** 2 + * (xi**2 * np.e ** (c * t / 2 * (1 - eta**2) * xi)).real + ) + + +def phi(x, t): + if t == 0.0 or abs(x) >= t: + return 0.0 + eta = x / t + integral = quad(integrand, 0.0, np.pi, args=(eta, t))[0] + return np.e**-t / 2 / t * (1 + c * t / 4 / np.pi * (1 - eta**2) * integral) + + +def phi_t(t, x): + if t == 0.0 or abs(x) >= t: + return 0.0 + eta = x / t + integral = quad(integrand, 0.0, np.pi, args=(eta, t))[0] + return np.e**-t / 2 / t * (1 + c * t / 4 / np.pi * (1 - eta**2) * integral) + + +def phiX(x, t0, t1): + return quad(phi_t, t0, t1, args=(x,))[0] + + +phi_avg = np.zeros([K, J]) + +for k in range(K): + for j in range(J): + x0 = x[j] + x1 = x[j + 1] + dx = x1 - x0 + t0 = t[k] + t1 = t[k + 1] + dt = t1 - t0 + phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt + +phi_avg = np.nan_to_num(phi_avg) + +np.savez("reference.npz", x=x, t=t, phi=phi_avg) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1/input.py b/verification/analytical/neutron/fixed_source/cases/azurv1/input.py index 1c6aa62..2620e7e 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1/input.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1/input.py @@ -1,6 +1,8 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("AZURV1") + # ====================================================================================== # Set model # ====================================================================================== @@ -9,7 +11,7 @@ # Effective scattering ratio c = 1.0 # Set materials -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=np.array([1.0 / 3.0]), scatter=np.array([[1.0 / 3.0]]), fission=np.array([1.0 / 3.0]), @@ -21,19 +23,21 @@ s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m) +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) # ====================================================================================== # Set source # ====================================================================================== # Isotropic pulse at x=t=0 -mcdc.Source( +source = mcdc.Source( position=[0.0, 0.0, 0.0], isotropic=True, - energy_group=0, + energy=0, time=0.0, ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, and run MC/DC @@ -41,11 +45,12 @@ # Tallies mesh = mcdc.MeshStructured(x=np.linspace(-20.5, 20.5, 202)) -mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 100 -mcdc.settings.N_batch = 2 +simulation.settings.N_particle = 100 +simulation.settings.N_batch = 2 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1/plot.py b/verification/analytical/neutron/fixed_source/cases/azurv1/plot.py index 143aaaf..0e57ebe 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1/plot.py @@ -54,6 +54,6 @@ def animate(k): return line1, line2, text -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1/process.py b/verification/analytical/neutron/fixed_source/cases/azurv1/process.py index 4ffa4f9..a7d6d03 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1/process.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1/process.py @@ -6,11 +6,11 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution data = np.load("reference.npz") diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1/reference.py b/verification/analytical/neutron/fixed_source/cases/azurv1/reference.py index 261ad68..d595570 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1/reference.py @@ -45,7 +45,7 @@ def phi_t(t, x): def phiX(x, t0, t1): - return quad(phi_t, t0, t1, args=(x))[0] + return quad(phi_t, t0, t1, args=(x,))[0] phi_avg = np.zeros([K, J]) @@ -60,12 +60,6 @@ def phiX(x, t0, t1): dt = t1 - t0 phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt -for j in range(J + 1): - for k in range(K): - t0 = t[k] - t1 = t[k + 1] - dt = t1 - t0 - phi_avg = np.nan_to_num(phi_avg) np.savez("reference.npz", x=x, t=t, phi=phi_avg) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/input.py b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/input.py index 42a08a0..e79d07b 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/input.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/input.py @@ -1,6 +1,8 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("AZURV1 subcritical") + # ====================================================================================== # Set model # ====================================================================================== @@ -9,7 +11,7 @@ # Effective scattering ratio c = 0.9 # Set materials -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=np.array([1.0 / 3.0]), scatter=np.array([[1.0 / 3.0]]), fission=np.array([1.0 / 3.0]), @@ -21,19 +23,21 @@ s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m) +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) # ====================================================================================== # Set source # ====================================================================================== # Isotropic pulse at x=t=0 -mcdc.Source( +source = mcdc.Source( position=[0.0, 0.0, 0.0], isotropic=True, - energy_group=0, + energy=0, time=0.0, ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, and run MC/DC @@ -41,11 +45,12 @@ # Tallies mesh = mcdc.MeshStructured(x=np.linspace(-20.5, 20.5, 202)) -mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 100 -mcdc.settings.N_batch = 2 +simulation.settings.N_particle = 100 +simulation.settings.N_batch = 2 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/plot.py b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/plot.py index 143aaaf..0e57ebe 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/plot.py @@ -54,6 +54,6 @@ def animate(k): return line1, line2, text -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/process.py b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/process.py index 4ffa4f9..a7d6d03 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/process.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/process.py @@ -6,11 +6,11 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution data = np.load("reference.npz") diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/reference.py b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/reference.py index ea6f7f1..7fffe9c 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_sub/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_sub/reference.py @@ -45,7 +45,7 @@ def phi_t(t, x): def phiX(x, t0, t1): - return quad(phi_t, t0, t1, args=(x))[0] + return quad(phi_t, t0, t1, args=(x,))[0] phi_avg = np.zeros([K, J]) @@ -60,12 +60,6 @@ def phiX(x, t0, t1): dt = t1 - t0 phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt -for j in range(J + 1): - for k in range(K): - t0 = t[k] - t1 = t[k + 1] - dt = t1 - t0 - phi_avg = np.nan_to_num(phi_avg) np.savez("reference.npz", x=x, t=t, phi=phi_avg) diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_super/input.py b/verification/analytical/neutron/fixed_source/cases/azurv1_super/input.py index 29f332f..10a724b 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_super/input.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_super/input.py @@ -1,6 +1,8 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("AZURV1 supercritical") + # ====================================================================================== # Set model # ====================================================================================== @@ -9,7 +11,7 @@ # Effective scattering ratio c = 1.1 # Set materials -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=np.array([1.0 / 3.0]), scatter=np.array([[1.0 / 3.0]]), fission=np.array([1.0 / 3.0]), @@ -21,19 +23,21 @@ s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m) +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) # ====================================================================================== # Set source # ====================================================================================== # Isotropic pulse at x=t=0 -mcdc.Source( +source = mcdc.Source( position=[0.0, 0.0, 0.0], isotropic=True, - energy_group=0, + energy=0, time=0.0, ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, and run MC/DC @@ -41,11 +45,12 @@ # Tallies mesh = mcdc.MeshStructured(x=np.linspace(-20.5, 20.5, 202)) -mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=np.linspace(0.0, 20.0, 21)) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 100 -mcdc.settings.N_batch = 2 +simulation.settings.N_particle = 100 +simulation.settings.N_batch = 2 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_super/plot.py b/verification/analytical/neutron/fixed_source/cases/azurv1_super/plot.py index 143aaaf..0e57ebe 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_super/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_super/plot.py @@ -54,6 +54,6 @@ def animate(k): return line1, line2, text -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_super/process.py b/verification/analytical/neutron/fixed_source/cases/azurv1_super/process.py index 4ffa4f9..a7d6d03 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_super/process.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_super/process.py @@ -6,11 +6,11 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution data = np.load("reference.npz") diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.npz b/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.npz index f4ccceb..3bfccd0 100644 Binary files a/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.npz and b/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.npz differ diff --git a/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.py b/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.py index 7f7deb4..33a3a20 100644 --- a/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/azurv1_super/reference.py @@ -45,7 +45,7 @@ def phi_t(t, x): def phiX(x, t0, t1): - return quad(phi_t, t0, t1, args=(x))[0] + return quad(phi_t, t0, t1, args=(x,))[0] phi_avg = np.zeros([K, J]) @@ -60,12 +60,6 @@ def phiX(x, t0, t1): dt = t1 - t0 phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt -for j in range(J + 1): - for k in range(K): - t0 = t[k] - t1 = t[k + 1] - dt = t1 - t0 - phi_avg = np.nan_to_num(phi_avg) np.savez("reference.npz", x=x, t=t, phi=phi_avg) diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/input.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/input.py new file mode 100644 index 0000000..7fdc240 --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/input.py @@ -0,0 +1,89 @@ +from pathlib import Path + +import numpy as np +import mcdc + +from reference import reference + +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" +WINDOW_FLOOR = 1.0e-3 +WINDOW_WIDTH = 2.5 + +simulation = mcdc.Simulation("Infinite SHEM-361 weight windows") + +# ====================================================================================== +# Set model +# ====================================================================================== +# The infinite homogeneous medium is modeled with a reflecting slab + +# Load material data +with np.load(SHEM361_DATA) as data: + SigmaC = data["SigmaC"] * 1.5 # /cm + SigmaS = data["SigmaS"] + SigmaF = data["SigmaF"] + nu_p = data["nu_p"] + nu_d = data["nu_d"] + chi_p = data["chi_p"] + chi_d = data["chi_d"] + G = int(data["G"]) + +# Set material +m = mcdc.Material.multigroup( + capture=SigmaC, + scatter=SigmaS, + fission=SigmaF, + nu_p=nu_p, + chi_p=chi_p, + nu_d=nu_d, + chi_d=chi_d, +) + +# Set surfaces +s1 = mcdc.Surface.PlaneX(x=-1e10, boundary_condition="reflective") +s2 = mcdc.Surface.PlaneX(x=1e10, boundary_condition="reflective") + +# Set cells +c = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([c]) + +# ====================================================================================== +# Set source +# ====================================================================================== + +source = mcdc.Source( + position=(0.0, 0.0, 0.0), + isotropic=True, + discrete_energy=np.array([[G - 1], [1.0]]), +) +simulation.set_sources([source]) + +# ====================================================================================== +# Set tallies, settings, techniques, and run MC/DC +# ====================================================================================== + +# Tallies +tally = mcdc.Tally(scores=["flux"], energy="all") +simulation.set_tallies([tally]) + +# Settings +simulation.settings.N_particle = 40 +simulation.settings.N_batch = 2 +simulation.settings.active_bank_buffer = 10_000 + +# Normalize the analytical group spectrum into target weights. +phi = reference() +target_weight = phi / np.max(phi) +target_weight = WINDOW_FLOOR + (1.0 - WINDOW_FLOOR) * target_weight + +# Center each standard multigroup coordinate in its corresponding window bin. +group_bounds = np.arange(G + 1, dtype=float) - 0.5 +weight_windows = np.empty((G, 1, 1, 1, 3)) +weight_windows[:, 0, 0, 0, 0] = target_weight / WINDOW_WIDTH +weight_windows[:, 0, 0, 0, 1] = target_weight +weight_windows[:, 0, 0, 0, 2] = target_weight * WINDOW_WIDTH + +# Techniques +simulation.technique.weight_windows(weight_windows, energy=group_bounds) + +# Run +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/plot.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/plot.py new file mode 100644 index 0000000..5395fcb --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/plot.py @@ -0,0 +1,33 @@ +import numpy as np +import matplotlib.pyplot as plt +import h5py +import sys +from pathlib import Path + +from reference import reference + +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + +# Load results +output = sys.argv[1] +with np.load(SHEM361_DATA) as data: + E = data["E"] + E_mid = 0.5 * (E[1:] + E[:-1]) + dE = E[1:] - E[:-1] +with h5py.File(output, "r") as f: + phi = f["tallies/tracklength_tally_0/flux/mean"][:] / dE * E_mid + phi_sd = f["tallies/tracklength_tally_0/flux/sdev"][:] / dE * E_mid + +# Reference solution +phi_ref = reference() / dE * E_mid + +# Flux +plt.plot(E_mid, phi, "-b", label="MC") +plt.fill_between(E_mid, phi - phi_sd, phi + phi_sd, alpha=0.2, color="b") +plt.plot(E_mid, phi_ref, "--r", label="analytical") +plt.xscale("log") +plt.xlabel(r"$E$, eV") +plt.ylabel(r"$E\phi(E)$") +plt.grid() +plt.legend() +plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/process.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/process.py new file mode 100644 index 0000000..adaf35b --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/process.py @@ -0,0 +1,32 @@ +from reference import reference +import numpy as np +import h5py +import sys + +sys.path.append("../../") +import util + +# Particle counts +N_min = int(sys.argv[1]) +N_max = int(sys.argv[2]) +N = int(sys.argv[3]) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) + +# Reference solution +phi_ref = reference() + +# Error container +error = np.zeros(len(N_particle_list)) +error_max = np.zeros(len(N_particle_list)) + +# Calculate error +for k, N_particle in enumerate(N_particle_list): + # Get results + with h5py.File("output_%i.h5" % (int(N_particle)), "r") as f: + phi = f["tallies/tracklength_tally_0/flux/mean"][:] + + error[k] = util.rerror(phi, phi_ref) + error_max[k] = util.rerror_max(phi, phi_ref) + +# Plot +util.plot_convergence("flux", N_particle_list, error, error_max) diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/reference.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/reference.py new file mode 100644 index 0000000..7f9045e --- /dev/null +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361-weight_windows/reference.py @@ -0,0 +1,23 @@ +from pathlib import Path + +import numpy as np + +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + + +def reference(): + # Load material data + with np.load(SHEM361_DATA) as data: + SigmaT = data["SigmaT"] + SigmaC = data["SigmaC"] + SigmaS = data["SigmaS"] + nuSigmaF = data["nuSigmaF"] + G = data["G"] + SigmaT += SigmaC * 0.5 + + A = np.diag(SigmaT) - SigmaS - nuSigmaF + Q = np.zeros(G) + Q[-1] = 1.0 + + phi = np.linalg.solve(A, Q) + return phi diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361/input.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361/input.py index 2068c3b..92c034b 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361/input.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361/input.py @@ -1,13 +1,19 @@ +from pathlib import Path + import numpy as np import mcdc +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + +simulation = mcdc.Simulation("Infinite SHEM-361") + # ====================================================================================== # Set model # ====================================================================================== -# The infinite homogenous medium is modeled with reflecting slab +# The infinite homogeneous medium is modeled with a reflecting slab # Load material data -with np.load("SHEM-361.npz") as data: +with np.load(SHEM361_DATA) as data: SigmaC = data["SigmaC"] * 1.5 # /cm SigmaS = data["SigmaS"] SigmaF = data["SigmaF"] @@ -15,12 +21,9 @@ nu_d = data["nu_d"] chi_p = data["chi_p"] chi_d = data["chi_d"] - G = data["G"] - speed = data["v"] - lamd = data["lamd"] # Set material -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=SigmaC, scatter=SigmaS, fission=SigmaF, @@ -36,26 +39,31 @@ # Set cells c = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([c]) # ====================================================================================== # Set source # ====================================================================================== -mcdc.Source( - position=(0.0, 0.0, 0.0), isotropic=True, energy_group=np.array([[360], [1.0]]) +source = mcdc.Source( + position=(0.0, 0.0, 0.0), + isotropic=True, + discrete_energy=np.array([[360], [1.0]]), ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, and run MC/DC # ====================================================================================== # Tallies -mcdc.Tally(scores=["flux"], energy="all_groups") +tally = mcdc.Tally(scores=["flux"], energy="all") +simulation.set_tallies([tally]) -# Swttings -mcdc.settings.N_particle = 40 -mcdc.settings.N_batch = 2 -mcdc.settings.active_bank_buffer = 1000 +# Settings +simulation.settings.N_particle = 40 +simulation.settings.N_batch = 2 +simulation.settings.active_bank_buffer = 1000 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361/plot.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361/plot.py index 91e4b69..5395fcb 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361/plot.py @@ -2,14 +2,16 @@ import matplotlib.pyplot as plt import h5py import sys +from pathlib import Path from reference import reference +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + # Load results output = sys.argv[1] -with np.load("../../data/MGXS-SHEM361.npz") as data: +with np.load(SHEM361_DATA) as data: E = data["E"] - G = data["G"] E_mid = 0.5 * (E[1:] + E[:-1]) dE = E[1:] - E[:-1] with h5py.File(output, "r") as f: diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361/process.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361/process.py index ccd6d68..adaf35b 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361/process.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361/process.py @@ -6,11 +6,11 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution phi_ref = reference() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361/reference.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361/reference.py index 5ac84f2..7f9045e 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361/reference.py @@ -1,9 +1,13 @@ +from pathlib import Path + import numpy as np +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + def reference(): # Load material data - with np.load("SHEM-361.npz") as data: + with np.load(SHEM361_DATA) as data: SigmaT = data["SigmaT"] SigmaC = data["SigmaC"] SigmaS = data["SigmaS"] diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/SHEM-361.npz b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/SHEM-361.npz deleted file mode 100644 index 6ba2c1f..0000000 Binary files a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/SHEM-361.npz and /dev/null differ diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/input.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/input.py index 6128d7e..24525ee 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/input.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/input.py @@ -1,13 +1,19 @@ +from pathlib import Path + import numpy as np import mcdc +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + +simulation = mcdc.Simulation("Infinite SHEM-361 time-dependent census") + # ====================================================================================== # Set model # ====================================================================================== -# The infinite homogenous medium is modeled with reflecting slab +# The infinite homogeneous medium is modeled with a reflecting slab # Load material data -with np.load("SHEM-361.npz") as data: +with np.load(SHEM361_DATA) as data: SigmaC = data["SigmaC"] * 1.5 # /cm SigmaS = data["SigmaS"] SigmaF = data["SigmaF"] @@ -15,12 +21,11 @@ nu_d = data["nu_d"] chi_p = data["chi_p"] chi_d = data["chi_d"] - G = data["G"] speed = data["v"] lamd = data["lamd"] # Set material -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=SigmaC, scatter=SigmaS, fission=SigmaF, @@ -38,36 +43,41 @@ # Set cells c = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([c]) # ====================================================================================== # Set source # ====================================================================================== -mcdc.Source( - position=(0.0, 0.0, 0.0), isotropic=True, energy_group=np.array([[360], [1.0]]) +source = mcdc.Source( + position=(0.0, 0.0, 0.0), + isotropic=True, + discrete_energy=np.array([[360], [1.0]]), ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, techniques, and run MC/DC # ====================================================================================== # Tallies -mcdc.Tally( +tally = mcdc.Tally( scores=["flux"], time=np.insert(np.logspace(-8, 1, 100), 0, 0.0), - energy="all_groups", + energy="all", ) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 20 -mcdc.settings.N_batch = 2 -mcdc.settings.set_time_census(np.logspace(-5, 1, 6)) -mcdc.settings.active_bank_buffer = 1000 -mcdc.settings.census_bank_buffer_ratio = 5.0 -mcdc.settings.source_bank_buffer_ratio = 5.0 +simulation.settings.N_particle = 20 +simulation.settings.N_batch = 2 +simulation.settings.set_time_census(np.logspace(-5, 1, 6)) +simulation.settings.active_bank_buffer = 1000 +simulation.settings.census_bank_buffer_ratio = 5.0 +simulation.settings.source_bank_buffer_ratio = 5.0 # Techniques -mcdc.simulation.population_control() +simulation.technique.population_control() # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/plot.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/plot.py index 88ec16f..83ad298 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/plot.py @@ -1,9 +1,13 @@ +from pathlib import Path + import numpy as np import matplotlib.pyplot as plt import h5py import sys import matplotlib.animation as animation +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + # Reference solution data = np.load("reference.npz") phi_ref = data["phi"].T @@ -11,9 +15,8 @@ # Load results output = sys.argv[1] -with np.load("../../data/MGXS-SHEM361.npz") as data: +with np.load(SHEM361_DATA) as data: E = data["E"] - G = data["G"] speed = data["v"] E_mid = 0.5 * (E[1:] + E[:-1]) dE = E[1:] - E[:-1] @@ -87,6 +90,6 @@ def animate(k): return line1, line2 -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/process.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/process.py index f540eb2..289307b 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/process.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/process.py @@ -1,15 +1,19 @@ +from pathlib import Path + import numpy as np import h5py import sys +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution data = np.load("reference.npz") @@ -26,9 +30,8 @@ # Calculate error for i, N_particle in enumerate(N_particle_list): # Get results - with np.load("SHEM-361.npz") as data: + with np.load(SHEM361_DATA) as data: E = data["E"] - G = data["G"] speed = data["v"] E_mid = 0.5 * (E[1:] + E[:-1]) dE = E[1:] - E[:-1] diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/reference.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/reference.py index 1ebacb0..267a0aa 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td-census/reference.py @@ -1,24 +1,25 @@ +from pathlib import Path + import numpy as np from scipy.linalg import expm +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + # Time grid t = np.insert(np.logspace(-8, 1, 100), 0, 0.0) K = len(t) - 1 # Load material data -with np.load("../../data/MGXS-SHEM361.npz") as data: +with np.load(SHEM361_DATA) as data: SigmaT = data["SigmaT"] SigmaC = data["SigmaC"] SigmaS = data["SigmaS"] nuSigmaF_p = data["nuSigmaF_p"] SigmaF = data["SigmaF"] - nu_p = data["nu_p"] nu_d = data["nu_d"] - chi_p = data["chi_p"] chi_d = data["chi_d"] G = data["G"] J = data["J"] - E = data["E"] v = data["v"] lamd = data["lamd"] SigmaT += SigmaC * 0.5 diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/SHEM-361.npz b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/SHEM-361.npz deleted file mode 100644 index 6ba2c1f..0000000 Binary files a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/SHEM-361.npz and /dev/null differ diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/input.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/input.py index 4a2ea50..ca87013 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/input.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/input.py @@ -1,13 +1,19 @@ +from pathlib import Path + import numpy as np import mcdc +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + +simulation = mcdc.Simulation("Infinite SHEM-361 time-dependent") + # ====================================================================================== # Set model # ====================================================================================== -# The infinite homogenous medium is modeled with reflecting slab +# The infinite homogeneous medium is modeled with a reflecting slab # Load material data -with np.load("SHEM-361.npz") as data: +with np.load(SHEM361_DATA) as data: SigmaC = data["SigmaC"] * 1.5 # /cm SigmaS = data["SigmaS"] SigmaF = data["SigmaF"] @@ -15,12 +21,11 @@ nu_d = data["nu_d"] chi_p = data["chi_p"] chi_d = data["chi_d"] - G = data["G"] speed = data["v"] lamd = data["lamd"] # Set material -m = mcdc.MaterialMG( +m = mcdc.Material.multigroup( capture=SigmaC, scatter=SigmaS, fission=SigmaF, @@ -38,30 +43,35 @@ # Set cells c = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([c]) # ====================================================================================== # Set source # ====================================================================================== -mcdc.Source( - position=(0.0, 0.0, 0.0), isotropic=True, energy_group=np.array([[360], [1.0]]) +source = mcdc.Source( + position=(0.0, 0.0, 0.0), + isotropic=True, + discrete_energy=np.array([[360], [1.0]]), ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, and run MC/DC # ====================================================================================== # Tallies -mcdc.Tally( +tally = mcdc.Tally( scores=["flux"], time=np.insert(np.logspace(-8, 1, 100), 0, 0.0), - energy="all_groups", + energy="all", ) +simulation.set_tallies([tally]) -# Swttings -mcdc.settings.N_particle = 50 -mcdc.settings.N_batch = 2 -mcdc.settings.active_bank_buffer = 1000 +# Settings +simulation.settings.N_particle = 50 +simulation.settings.N_batch = 2 +simulation.settings.active_bank_buffer = 1000 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/plot.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/plot.py index 88ec16f..83ad298 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/plot.py @@ -1,9 +1,13 @@ +from pathlib import Path + import numpy as np import matplotlib.pyplot as plt import h5py import sys import matplotlib.animation as animation +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + # Reference solution data = np.load("reference.npz") phi_ref = data["phi"].T @@ -11,9 +15,8 @@ # Load results output = sys.argv[1] -with np.load("../../data/MGXS-SHEM361.npz") as data: +with np.load(SHEM361_DATA) as data: E = data["E"] - G = data["G"] speed = data["v"] E_mid = 0.5 * (E[1:] + E[:-1]) dE = E[1:] - E[:-1] @@ -87,6 +90,6 @@ def animate(k): return line1, line2 -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/process.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/process.py index f540eb2..289307b 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/process.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/process.py @@ -1,15 +1,19 @@ +from pathlib import Path + import numpy as np import h5py import sys +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution data = np.load("reference.npz") @@ -26,9 +30,8 @@ # Calculate error for i, N_particle in enumerate(N_particle_list): # Get results - with np.load("SHEM-361.npz") as data: + with np.load(SHEM361_DATA) as data: E = data["E"] - G = data["G"] speed = data["v"] E_mid = 0.5 * (E[1:] + E[:-1]) dE = E[1:] - E[:-1] diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/reference.py b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/reference.py index 1ebacb0..267a0aa 100644 --- a/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/inf_shem361_td/reference.py @@ -1,24 +1,25 @@ +from pathlib import Path + import numpy as np from scipy.linalg import expm +SHEM361_DATA = Path(__file__).resolve().parents[2] / "data" / "SHEM-361.npz" + # Time grid t = np.insert(np.logspace(-8, 1, 100), 0, 0.0) K = len(t) - 1 # Load material data -with np.load("../../data/MGXS-SHEM361.npz") as data: +with np.load(SHEM361_DATA) as data: SigmaT = data["SigmaT"] SigmaC = data["SigmaC"] SigmaS = data["SigmaS"] nuSigmaF_p = data["nuSigmaF_p"] SigmaF = data["SigmaF"] - nu_p = data["nu_p"] nu_d = data["nu_d"] - chi_p = data["chi_p"] chi_d = data["chi_d"] G = data["G"] J = data["J"] - E = data["E"] v = data["v"] lamd = data["lamd"] SigmaT += SigmaC * 0.5 diff --git a/verification/analytical/neutron/fixed_source/cases/reed/README.md b/verification/analytical/neutron/fixed_source/cases/reed/README.md deleted file mode 100644 index 28ffc0a..0000000 --- a/verification/analytical/neutron/fixed_source/cases/reed/README.md +++ /dev/null @@ -1 +0,0 @@ -Four slab layers with different materials. Based on William H. Reed, NSE (1971), 46:2, 309-314, DOI: 10.13182/NSE46-309. Reference solution is generated from OpenMC run with 30 batches and 1000M particles per batch. diff --git a/verification/analytical/neutron/fixed_source/cases/reed/build-xml.py b/verification/analytical/neutron/fixed_source/cases/reed/build-xml.py deleted file mode 100644 index ef424a4..0000000 --- a/verification/analytical/neutron/fixed_source/cases/reed/build-xml.py +++ /dev/null @@ -1,181 +0,0 @@ -import openmc -import numpy as np -import h5py, sys - -# =========================================================================== -# Set Library -# =========================================================================== - -SigmaC = 1.0 / 3.0 -SigmaF = 1.0 / 3.0 -nu = 2.3 -SigmaA = SigmaC + SigmaF -SigmaS = 1.0 / 3.0 -SigmaT = SigmaA + SigmaS -v = 1.0 - -m1_SigmaT = 50.0 -m2_SigmaT = 5.0 -m3_SigmaT = 0.0 -m4_SigmaT = 1.0 - -m1_SigmaA = 50.0 -m2_SigmaA = 5.0 -m3_SigmaA = 0.0 -m4_SigmaA = 0.1 - -m1_SigmaS = m1_SigmaT - m1_SigmaA -m2_SigmaS = m2_SigmaT - m2_SigmaA -m3_SigmaS = m3_SigmaT - m3_SigmaA -m4_SigmaS = m4_SigmaT - m4_SigmaA - -groups = openmc.mgxs.EnergyGroups([0.0, 2e7]) - -m1 = openmc.XSdata("m1", groups) -m2 = openmc.XSdata("m2", groups) -m3 = openmc.XSdata("m3", groups) -m4 = openmc.XSdata("m4", groups) - -m1.order = 0 -m2.order = 0 -m3.order = 0 -m4.order = 0 - -m1.set_total([m1_SigmaT], temperature=294.0) -m2.set_total([m2_SigmaT], temperature=294.0) -m3.set_total([m3_SigmaT], temperature=294.0) -m4.set_total([m4_SigmaT], temperature=294.0) - -m1.set_absorption([m1_SigmaA], temperature=294.0) -m2.set_absorption([m2_SigmaA], temperature=294.0) -m3.set_absorption([m3_SigmaA], temperature=294.0) -m4.set_absorption([m4_SigmaA], temperature=294.0) - -m1.set_scatter_matrix(np.ones((1, 1, 1)) * m1_SigmaS, temperature=294.0) -m2.set_scatter_matrix(np.ones((1, 1, 1)) * m2_SigmaS, temperature=294.0) -m3.set_scatter_matrix(np.ones((1, 1, 1)) * m3_SigmaS, temperature=294.0) -m4.set_scatter_matrix(np.ones((1, 1, 1)) * m4_SigmaS, temperature=294.0) - -mg_cross_sections_file = openmc.MGXSLibrary(groups) -mg_cross_sections_file.add_xsdata(m1) -mg_cross_sections_file.add_xsdata(m2) -mg_cross_sections_file.add_xsdata(m3) -mg_cross_sections_file.add_xsdata(m4) -mg_cross_sections_file.export_to_hdf5("mgxs.h5") - -# =========================================================================== -# Exporting to OpenMC materials.xml file -# =========================================================================== - -materials = {} -materials["m1"] = openmc.Material(name="m1") -materials["m1"].set_density("macro", 1.0) -materials["m1"].add_macroscopic("m1") - -materials["m2"] = openmc.Material(name="m2") -materials["m2"].set_density("macro", 1.0) -materials["m2"].add_macroscopic("m2") - -materials["m3"] = openmc.Material(name="m3") -materials["m3"].set_density("macro", 1.0) -materials["m3"].add_macroscopic("m3") - -materials["m4"] = openmc.Material(name="m4") -materials["m4"].set_density("macro", 1.0) -materials["m4"].add_macroscopic("m4") - -materials_file = openmc.Materials(materials.values()) -materials_file.cross_sections = "mgxs.h5" -materials_file.export_to_xml() - -# =========================================================================== -# Exporting to OpenMC geometry.xml file -# =========================================================================== - -# Instantiate ZCylinder surfaces -s1 = openmc.XPlane(x0=0.0, boundary_type="reflective") -s2 = openmc.XPlane(x0=2.0) -s3 = openmc.XPlane(x0=3.0) -s4 = openmc.XPlane(x0=5.0) -s5 = openmc.XPlane(x0=8.0, boundary_type="vacuum") - -# Instantiate Cells -c1 = openmc.Cell() -c2 = openmc.Cell() -c3 = openmc.Cell() -c4 = openmc.Cell() -# Use surface half-spaces to define regions -c1.region = +s1 & -s2 -c2.region = +s2 & -s3 -c3.region = +s3 & -s4 -c4.region = +s4 & -s5 -# Register Materials with Cells -c1.fill = materials["m1"] -c2.fill = materials["m2"] -c3.fill = materials["m3"] -c4.fill = materials["m4"] - -# Instantiate Universes -root = openmc.Universe(universe_id=0, name="root universe", cells=[c1, c2, c3, c4]) - -# Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry(root) -geometry.export_to_xml() - -# =========================================================================== -# Exporting to OpenMC settings.xml file -# =========================================================================== - -# Instantiate a Settings object, set all runtime parameters, and export to XML -settings_file = openmc.Settings() -settings_file.run_mode = "fixed source" -settings_file.particles = 1000000000 -settings_file.batches = 30 -settings_file.output = {"tallies": False} -settings_file.energy_mode = "multi-group" - -# Define source -src1 = openmc.IndependentSource() -src2 = openmc.IndependentSource() -# -src1.space = openmc.stats.CartesianIndependent( - x=openmc.stats.Uniform(0.0, 2.0), - y=openmc.stats.Discrete([0.0], [1.0]), - z=openmc.stats.Discrete([0.0], [1.0]), -) -src2.space = openmc.stats.CartesianIndependent( - x=openmc.stats.Uniform(5.0, 6.0), - y=openmc.stats.Discrete([0.0], [1.0]), - z=openmc.stats.Discrete([0.0], [1.0]), -) -# -src1.angle = openmc.stats.Isotropic() -src2.angle = openmc.stats.Isotropic() -# -src1.strength = 50.0 / 50.5 -src2.strength = 0.5 / 50.5 - -settings_file.source = [src1, src2] -settings_file.export_to_xml() - - -# =========================================================================== -# Exporting to OpenMC tallies.xml file -# =========================================================================== - -# Create a mesh filter that can be used in a tally -mesh = openmc.RectilinearMesh() -mesh.x_grid = np.linspace(0.0, 8.0, 81) -mesh.y_grid = np.linspace(-1e15, 1e15, 2) -mesh.z_grid = np.linspace(-1e15, 1e15, 2) -mesh_filter = openmc.MeshFilter(mesh) - -# Now use the mesh filter in a tally and indicate what scores are desired -tally1 = openmc.Tally(name="flux") -tally1.estimator = "tracklength" -tally1.filters = [mesh_filter] -tally1.scores = ["flux"] - -# Instantiate a Tallies collection and export to XML -tallies = openmc.Tallies([tally1]) -tallies.export_to_xml() diff --git a/verification/analytical/neutron/fixed_source/cases/reed/input.py b/verification/analytical/neutron/fixed_source/cases/reed/input.py index adc57bd..12ce9c7 100644 --- a/verification/analytical/neutron/fixed_source/cases/reed/input.py +++ b/verification/analytical/neutron/fixed_source/cases/reed/input.py @@ -1,17 +1,19 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("Reed slab") + # ====================================================================================== # Set model # ====================================================================================== -# Three slab layers with different materials +# Four slab layers with different materials # Based on William H. Reed, NSE (1971), 46:2, 309-314, DOI: 10.13182/NSE46-309 # Set materials -m1 = mcdc.MaterialMG(capture=np.array([50.0])) -m2 = mcdc.MaterialMG(capture=np.array([5.0])) -m3 = mcdc.MaterialMG(capture=np.array([0.0])) # Vacuum -m4 = mcdc.MaterialMG(capture=np.array([0.1]), scatter=np.array([[0.9]])) +m1 = mcdc.Material.multigroup(capture=np.array([50.0])) +m2 = mcdc.Material.multigroup(capture=np.array([5.0])) +m3 = mcdc.Material.multigroup(capture=np.array([0.0])) # Vacuum +m4 = mcdc.Material.multigroup(capture=np.array([0.1]), scatter=np.array([[0.9]])) # Set surfaces s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="reflective") @@ -21,21 +23,23 @@ s5 = mcdc.Surface.PlaneZ(z=8.0, boundary_condition="vacuum") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m1) -mcdc.Cell(region=+s2 & -s3, fill=m2) -mcdc.Cell(region=+s3 & -s4, fill=m3) -mcdc.Cell(region=+s4 & -s5, fill=m4) +cell_1 = mcdc.Cell(region=+s1 & -s2, fill=m1) +cell_2 = mcdc.Cell(region=+s2 & -s3, fill=m2) +cell_3 = mcdc.Cell(region=+s3 & -s4, fill=m3) +cell_4 = mcdc.Cell(region=+s4 & -s5, fill=m4) +simulation.set_model([cell_1, cell_2, cell_3, cell_4]) # ====================================================================================== # Set source # ====================================================================================== # Isotropic source in the absorbing medium -mcdc.Source(z=[0.0, 2.0], isotropic=True, energy_group=0, probability=50.0) +source_1 = mcdc.Source(z=[0.0, 2.0], isotropic=True, energy=0, probability=50.0) # Isotropic source in the first half of the outermost medium, # with 1/100 strength -mcdc.Source(z=[5.0, 6.0], isotropic=True, energy_group=0, probability=0.5) +source_2 = mcdc.Source(z=[5.0, 6.0], isotropic=True, energy=0, probability=0.5) +simulation.set_sources([source_1, source_2]) # ====================================================================================== # Set tallies, settings, and run MC/DC @@ -43,11 +47,12 @@ # Tallies mesh = mcdc.MeshStructured(z=np.linspace(0.0, 8.0, 81)) -mcdc.Tally(mesh=mesh, scores=["flux"]) +tally = mcdc.Tally(mesh=mesh, scores=["flux"]) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 5000 -mcdc.settings.N_batch = 2 +simulation.settings.N_particle = 5000 +simulation.settings.N_batch = 2 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/reed/plot.py b/verification/analytical/neutron/fixed_source/cases/reed/plot.py index 27d24f9..ee883ca 100644 --- a/verification/analytical/neutron/fixed_source/cases/reed/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/reed/plot.py @@ -1,25 +1,28 @@ -import matplotlib.pyplot as plt -import h5py +from pathlib import Path import sys +import h5py +import matplotlib.pyplot as plt + from reference import reference -# Reference solution +# Locate case data independently of the caller's working directory. +CASE_DIR = Path(__file__).resolve().parent output = sys.argv[1] +# Analytical reference solution +z_ref, phi_ref = reference() + # Load results with h5py.File(output, "r") as f: z = f["tallies/tracklength_tally_0/grid/z"][:] dz = z[1:] - z[:-1] z_mid = 0.5 * (z[:-1] + z[1:]) - I = len(z) - 1 - phi = f["tallies/tracklength_tally_0/flux/mean"][:] phi_sd = f["tallies/tracklength_tally_0/flux/sdev"][:] # Normalize phi, phi_sd = phi / dz, phi_sd / dz -z_ref, phi_ref = reference() # Flux - spatial average plt.plot(z_mid, phi, "-b", label="MC") @@ -30,5 +33,5 @@ plt.grid() plt.legend() plt.title(r"$\bar{\phi}_i$") -plt.savefig("scalar_flux.png", dpi=300) +plt.savefig(CASE_DIR / "scalar_flux.png", dpi=300) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/reed/process.py b/verification/analytical/neutron/fixed_source/cases/reed/process.py index 0084279..53d7a72 100644 --- a/verification/analytical/neutron/fixed_source/cases/reed/process.py +++ b/verification/analytical/neutron/fixed_source/cases/reed/process.py @@ -2,18 +2,19 @@ import h5py import sys +from reference import reference + sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) -# Reference solution -data = np.load("reference.npz") -phi_ref = data["phi"] +# Analytical reference solution +_, phi_ref = reference() # Error containers error = np.zeros(len(N_particle_list)) diff --git a/verification/analytical/neutron/fixed_source/cases/reed/reference.npz b/verification/analytical/neutron/fixed_source/cases/reed/reference.npz deleted file mode 100644 index 20f52f7..0000000 Binary files a/verification/analytical/neutron/fixed_source/cases/reed/reference.npz and /dev/null differ diff --git a/verification/analytical/neutron/fixed_source/cases/reed/reference.py b/verification/analytical/neutron/fixed_source/cases/reed/reference.py index c3bccc6..f06dd5d 100644 --- a/verification/analytical/neutron/fixed_source/cases/reed/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/reed/reference.py @@ -3,7 +3,7 @@ # Warsa, J.S., "Analytical SN solutions in heterogeneous slabs, using symbolic algebra computer programs" # Ann. Nucl. Energy, 29(7), 851-874 (2002). DOI: 10.1016/S0306-4549(01)00080-9 -# reproduced reference equations from Warsa, 2002 +# Reproduce the reference equations from Warsa (2002). def reference(): diff --git a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/input.py b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/input.py index e6bb936..4ead42c 100644 --- a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/input.py +++ b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/input.py @@ -1,15 +1,17 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("Slab absorbium") + # ====================================================================================== # Set model # ====================================================================================== # Three slab layers with different purely-absorbing materials # Set materials -m1 = mcdc.MaterialMG(capture=np.array([1.0])) -m2 = mcdc.MaterialMG(capture=np.array([1.5])) -m3 = mcdc.MaterialMG(capture=np.array([2.0])) +m1 = mcdc.Material.multigroup(capture=np.array([1.0])) +m2 = mcdc.Material.multigroup(capture=np.array([1.5])) +m3 = mcdc.Material.multigroup(capture=np.array([2.0])) # Set surfaces s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") @@ -18,29 +20,36 @@ s4 = mcdc.Surface.PlaneZ(z=6.0, boundary_condition="vacuum") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m2) -mcdc.Cell(region=+s2 & -s3, fill=m3) -mcdc.Cell(region=+s3 & -s4, fill=m1) +cell_1 = mcdc.Cell(region=+s1 & -s2, fill=m2) +cell_2 = mcdc.Cell(region=+s2 & -s3, fill=m3) +cell_3 = mcdc.Cell(region=+s3 & -s4, fill=m1) +simulation.set_model([cell_1, cell_2, cell_3]) # ====================================================================================== # Set source # ====================================================================================== # Uniform isotropic source throughout the domain -mcdc.Source(z=[0.0, 6.0], isotropic=True, energy_group=0) +source = mcdc.Source(z=[0.0, 6.0], isotropic=True, energy=0) +simulation.set_sources([source]) # ====================================================================================== -# Set tallies, settings, and run mcdc +# Set tallies, settings, and run MC/DC # ====================================================================================== # Tallies -mcdc.Tally(surface=s4, scores=["net-current"]) +surface_tally = mcdc.Tally(surface=s4, scores=["current-net"]) mesh = mcdc.MeshStructured(z=np.linspace(0.0, 6.0, 61)) -mcdc.Tally(mesh=mesh, mu=np.linspace(-1.0, 1.0, 32 + 1), scores=["flux", "collision"]) +mesh_tally = mcdc.Tally( + mesh=mesh, + mu=np.linspace(-1.0, 1.0, 32 + 1), + scores=["flux", "collision"], +) +simulation.set_tallies([surface_tally, mesh_tally]) # Settings -mcdc.settings.N_particle = 100 -mcdc.settings.N_batch = 2 +simulation.settings.N_particle = 100 +simulation.settings.N_batch = 2 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/plot.py b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/plot.py index f936169..90d9285 100644 --- a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/plot.py @@ -10,7 +10,7 @@ with h5py.File(output, "r") as f: z = f["tallies/tracklength_tally_0/grid/z"][:] mu = f["tallies/tracklength_tally_0/grid/mu"][:] -phi_ref, J_ref, psi_ref = reference(z, mu) +phi_ref, _, psi_ref = reference(z, mu) # Load results with h5py.File(output, "r") as f: diff --git a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/process.py b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/process.py index 9f58660..fc2b4c8 100644 --- a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/process.py +++ b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/process.py @@ -6,17 +6,17 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution with h5py.File("output_%i.h5" % (int(N_particle_list[0])), "r") as f: z = f["tallies/tracklength_tally_0/grid/z"][:] mu = f["tallies/tracklength_tally_0/grid/mu"][:] -phi_ref, J_ref, psi_ref = reference(z, mu) +phi_ref, _, psi_ref = reference(z, mu) # Error containers error = np.zeros(len(N_particle_list)) @@ -34,7 +34,7 @@ mu = f["tallies/tracklength_tally_0/grid/mu"][:] dmu = mu[1:] - mu[:-1] I = len(z) - 1 - N = len(mu) - 1 + N_mu = len(mu) - 1 psi = f["tallies/tracklength_tally_0/flux/mean"][:] psi = np.transpose(psi) @@ -44,10 +44,9 @@ for i in range(I): phi[i] += np.sum(psi[i, :]) - psi_norm = np.zeros(psi.shape) # Normalize phi /= dz - for n in range(N): + for n in range(N_mu): psi[:, n] = psi[:, n] / dz / dmu[n] # Get error @@ -60,6 +59,4 @@ # Plot util.plot_convergence("flux", N_particle_list, error, error_max) -util.plot_convergence( - "flux", N_particle_list, error_psi, error_max_psi -) +util.plot_convergence("angular_flux", N_particle_list, error_psi, error_max_psi) diff --git a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/reference.py b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/reference.py index b1bfaa5..15b6efd 100644 --- a/verification/analytical/neutron/fixed_source/cases/slab_absorbium/reference.py +++ b/verification/analytical/neutron/fixed_source/cases/slab_absorbium/reference.py @@ -48,13 +48,13 @@ def psi3(mu, x): # Flux def phi1(x): - return quad(psi1, -1, 1, args=(x), points=[0.0])[0] + return quad(psi1, -1, 1, args=(x,), points=[0.0])[0] def phi2(x): - return quad(psi2, -1, 1, args=(x), points=[0.0])[0] + return quad(psi2, -1, 1, args=(x,), points=[0.0])[0] def phi3(x): - return quad(psi3, -1, 1, args=(x), points=[0.0])[0] + return quad(psi3, -1, 1, args=(x,), points=[0.0])[0] # Integrands for current def mu_psi1(mu, x): @@ -77,23 +77,23 @@ def mu_psi3(mu, x): # Current def J1(x): - return quad(mu_psi1, -1, 1, args=(x), points=[0.0])[0] + return quad(mu_psi1, -1, 1, args=(x,), points=[0.0])[0] def J2(x): - return quad(mu_psi2, -1, 1, args=(x), points=[0.0])[0] + return quad(mu_psi2, -1, 1, args=(x,), points=[0.0])[0] def J3(x): - return quad(mu_psi3, -1, 1, args=(x), points=[0.0])[0] + return quad(mu_psi3, -1, 1, args=(x,), points=[0.0])[0] # Angular flux def psi1_(x, mu0, mu1): - return quad(psi1, mu0, mu1, args=(x), points=[0.0])[0] + return quad(psi1, mu0, mu1, args=(x,), points=[0.0])[0] def psi2_(x, mu0, mu1): - return quad(psi2, mu0, mu1, args=(x), points=[0.0])[0] + return quad(psi2, mu0, mu1, args=(x,), points=[0.0])[0] def psi3_(x, mu0, mu1): - return quad(psi3, mu0, mu1, args=(x), points=[0.0])[0] + return quad(psi3, mu0, mu1, args=(x,), points=[0.0])[0] phi = np.zeros(I) psi = np.zeros((I, N)) @@ -120,8 +120,4 @@ def psi3_(x, mu0, mu1): mu0 = mu[n] mu1 = mu[n + 1] psi[i, n] = quad(psi3_, x[i], x[i + 1], args=(mu0, mu1))[0] / dx[i] / dmu[n] - for n in range(N): - mu0 = mu[n] - mu1 = mu[n + 1] - return phi, J, psi diff --git a/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/input.py b/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/input.py index 987b59a..805ece8 100644 --- a/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/input.py +++ b/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/input.py @@ -2,32 +2,36 @@ import mcdc +simulation = mcdc.Simulation("Slab isotropic beam time dependent") + # ====================================================================================== # Set model # ====================================================================================== # Finite homogeneous pure-absorbing slab # Set materials -m = mcdc.MaterialMG(capture=np.array([1.0])) +m = mcdc.Material.multigroup(capture=np.array([1.0])) # Set surfaces s1 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="vacuum") s2 = mcdc.Surface.PlaneX(x=5.0, boundary_condition="vacuum") # Set cells -mcdc.Cell(region=+s1 & -s2, fill=m) +cell = mcdc.Cell(region=+s1 & -s2, fill=m) +simulation.set_model([cell]) # ====================================================================================== # Set source # ====================================================================================== # Isotropic beam from left-end -mcdc.Source( +source = mcdc.Source( position=(1e-10, 0.0, 0.0), white_direction=(1.0, 0.0, 0.0), - energy_group=0, + energy=0, time=[0.0, 5.0], ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, and run MC/DC @@ -35,15 +39,16 @@ # Tallies mesh = mcdc.MeshUniform(x=(0.0, 0.1, 50)) -mcdc.Tally( +tally = mcdc.Tally( mesh=mesh, scores=["flux"], time=np.linspace(0.0, 5.0, 51), ) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = 100 -mcdc.settings.N_batch = 2 +simulation.settings.N_particle = 100 +simulation.settings.N_batch = 2 # Run -mcdc.run() +simulation.run() diff --git a/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/plot.py b/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/plot.py index 418a72a..b78c26c 100644 --- a/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/plot.py +++ b/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/plot.py @@ -65,6 +65,6 @@ def animate(k): return line1, line2, text -simulation = animation.FuncAnimation(fig, animate, frames=K) -writervideo = animation.FFMpegWriter(fps=6) +# Uncomment to animate all time bins. +# _simulation = animation.FuncAnimation(fig, animate, frames=K) plt.show() diff --git a/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/process.py b/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/process.py index 2e46516..9273e97 100644 --- a/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/process.py +++ b/verification/analytical/neutron/fixed_source/cases/slab_isobeam_td/process.py @@ -6,11 +6,11 @@ sys.path.append("../../") import util -# Cases run +# Particle counts N_min = int(sys.argv[1]) N_max = int(sys.argv[2]) N = int(sys.argv[3]) -N_particle_list = np.logspace(N_min, N_max, N) +N_particle_list = np.logspace(N_min, N_max, N, dtype=int) # Reference solution with h5py.File("output_%i.h5" % (int(N_particle_list[0])), "r") as f: diff --git a/verification/analytical/neutron/fixed_source/cases/inf_shem361/SHEM-361.npz b/verification/analytical/neutron/fixed_source/data/SHEM-361.npz similarity index 100% rename from verification/analytical/neutron/fixed_source/cases/inf_shem361/SHEM-361.npz rename to verification/analytical/neutron/fixed_source/data/SHEM-361.npz diff --git a/verification/analytical/neutron/fixed_source/launch.py b/verification/analytical/neutron/fixed_source/launch.py index 7beaed9..953c417 100644 --- a/verification/analytical/neutron/fixed_source/launch.py +++ b/verification/analytical/neutron/fixed_source/launch.py @@ -1,3 +1,5 @@ +"""Build and launch the analytical fixed-source verification study with Maestro.""" + import argparse import os import subprocess @@ -5,7 +7,6 @@ import yaml - # ====================================================================================== # Bootstrap VVP imports # ====================================================================================== @@ -14,6 +15,7 @@ REPO_DIR = Path(__file__).resolve().parents[4] +# Support launching this suite directly without installing MC/DC-VVP. if str(REPO_DIR) not in sys.path: sys.path.insert(0, str(REPO_DIR)) @@ -22,10 +24,11 @@ # Load shared VVP configs # ====================================================================================== -from configs.platform_config import PLATFORMS # noqa: E402 +from configs.platform_config import PLATFORMS +# User overrides are optional for local runs. try: - from configs.user_config import USER_CONFIG # noqa: E402 + from configs.user_config import USER_CONFIG except ImportError: USER_CONFIG = {} @@ -66,6 +69,7 @@ mcdc_python = user_platform_config.get("mcdc_python") +# Use the active interpreter unless this platform specifies another MC/DC environment. if mcdc_python is None: mcdc_python = sys.executable else: @@ -76,6 +80,7 @@ scheduler = platform["scheduler"] cpu_cores = platform["cpu_cores_per_node"] + # Never request more walltime than the platform permits. walltime_hours = ( platform["max_walltime_hours"] if args.walltime is None @@ -109,10 +114,12 @@ steps = [] for case_name in tasks: + # Normalize case names into stable Maestro step identifiers. safe_case_name = case_name.replace("-", "_") command = f"{mcdc_python} {run_case} --name {case_name}" + # Maestro replaces LAUNCHER with the scheduler-specific MPI launch command. if args.mpi: command += ' --mpi "$(LAUNCHER)"' @@ -121,6 +128,7 @@ run = {"cmd": command} + # Scheduled studies require explicit resources; local studies run directly. if not local: run["nodes"] = 1 run["walltime"] = walltime @@ -151,6 +159,7 @@ } if not local: + # Attach batch settings only when Maestro submits to a scheduler. batch = { "type": scheduler, "host": platform["host"], @@ -185,11 +194,13 @@ env = os.environ.copy() +# Use the configured Maestro environment when it differs from the active one. if maestro_python is None: maestro_command = ["maestro", "run", "study.yaml"] else: maestro_python = Path(maestro_python).expanduser() maestro_bin = maestro_python.parent + # Keep executables spawned by Maestro in the same configured environment. env["PATH"] = f"{maestro_bin}:{env['PATH']}" maestro_command = [ @@ -207,6 +218,7 @@ # Store launch metadata # ====================================================================================== +# Maestro creates timestamped run directories, so capture the newly generated launch. maestro_runs = sorted( suite_dir.glob("maestro_run_*"), key=lambda path: path.stat().st_mtime, @@ -224,6 +236,7 @@ "rewrite": args.rewrite, } +# Snapshot the effective launch and task configuration with the generated run. with (latest_run / "launch_config.yaml").open("w") as f: yaml.dump(launch_config, f, sort_keys=False) diff --git a/verification/analytical/neutron/fixed_source/process.py b/verification/analytical/neutron/fixed_source/process.py index 851ebb4..3e869a2 100644 --- a/verification/analytical/neutron/fixed_source/process.py +++ b/verification/analytical/neutron/fixed_source/process.py @@ -5,7 +5,6 @@ import yaml - # ====================================================================================== # Command-line arguments # ====================================================================================== @@ -27,7 +26,6 @@ # ====================================================================================== suite_dir = Path(__file__).resolve().parent -task_file = suite_dir / "task.yaml" if args.maestro_run is None: maestro_runs = sorted( @@ -46,10 +44,14 @@ maestro_run = suite_dir / maestro_run launch_config_file = maestro_run / "launch_config.yaml" +task_file = maestro_run / "task.yaml" if not launch_config_file.is_file(): raise FileNotFoundError(f"Launch config not found: {launch_config_file}") +if not task_file.is_file(): + raise FileNotFoundError(f"Task config not found: {task_file}") + # ====================================================================================== # Load configuration and tasks diff --git a/verification/analytical/neutron/fixed_source/run_case.py b/verification/analytical/neutron/fixed_source/run_case.py index dcec0e1..1851abe 100644 --- a/verification/analytical/neutron/fixed_source/run_case.py +++ b/verification/analytical/neutron/fixed_source/run_case.py @@ -6,7 +6,6 @@ import numpy as np import yaml - # ====================================================================================== # Command-line arguments # ====================================================================================== diff --git a/verification/analytical/neutron/fixed_source/task.yaml b/verification/analytical/neutron/fixed_source/task.yaml index 0184999..13ec96f 100644 --- a/verification/analytical/neutron/fixed_source/task.yaml +++ b/verification/analytical/neutron/fixed_source/task.yaml @@ -46,6 +46,16 @@ azurv1: logN_max: 7 N_task: 7 +azurv1-basic_techniques: + logN_min: 5 + logN_max: 7 + N_task: 7 + +azurv1-weight_windows: + logN_min: 4 + logN_max: 6 + N_task: 7 + azurv1-census: logN_min: 5 logN_max: 7 @@ -81,6 +91,11 @@ inf_shem361: logN_max: 7 N_task: 7 +inf_shem361-weight_windows: + logN_min: 3 + logN_max: 5 + N_task: 7 + inf_shem361_td: logN_min: 5 logN_max: 7 diff --git a/verification/analytical/neutron/k_eigenvalue/.gitignore b/verification/analytical/neutron/k_eigenvalue/.gitignore new file mode 100644 index 0000000..1d0875e --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/.gitignore @@ -0,0 +1,4 @@ +study.yaml +maestro_run*/ +results/ + diff --git a/verification/analytical/neutron/k_eigenvalue/README.md b/verification/analytical/neutron/k_eigenvalue/README.md new file mode 100644 index 0000000..d701724 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/README.md @@ -0,0 +1,93 @@ +# Analytical Neutron $k$-Eigenvalue Verification + +This suite verifies MC/DC neutron $k$-eigenvalue calculations using infinite homogeneous multigroup problems with analytical reference solutions. +It compares the estimated multiplication factor and fundamental-mode energy spectrum against the dominant solution of the corresponding generalized matrix eigenvalue problem. +Verification checks for the expected $N_\mathrm{active}^{-1/2}$ statistical convergence as the number of active cycles increases. + +The number of particles per cycle and the number of inactive cycles are fixed in each case input so that the study isolates convergence with active cycles. +The suite can be executed independently or as part of the top-level MC/DC-VVP workflow. + +## Directory layout + +```text +cases/ Verification case definitions and processing scripts +maestro_run_*/ Generated Maestro workflow directories +results/ Generated figures from processed cases + +task.yaml Configure the active-cycle study for each case +study.yaml Generated Maestro study definition + +launch.py Build and launch the Maestro study +run_case.py Run one case over its active-cycle study +process.py Process a completed Maestro study + +cleanup.py Remove generated outputs and figures +util.py Provide shared processing and plotting utilities +``` + +The cases use the SHEM-361 dataset in the neighboring fixed-source suite at `../fixed_source/data/SHEM-361.npz`. + +## Configuration + +The `task.yaml` file selects the cases and defines the minimum and maximum active-cycle counts and the number of tasks for each case. +The active-cycle counts are geometrically spaced so that convergence can be assessed efficiently over a range of tens to hundreds of cycles. + +The fixed particle and inactive-cycle counts are defined in each case's `input.py` and should be tuned with a preliminary source-convergence and variance study before running the active-cycle campaign. +HPC runs use the shared platform settings in the repository's `configs/platform_config.py` and the user-specific settings in `configs/user_config.py`. + +## Launching and processing + +From this suite directory, launch the study locally: + +```bash +python launch.py +``` + +Launch the study on a supported HPC platform: + +```bash +python launch.py --platform tuolumne --mpi +``` + +Use `--walltime HOURS` to limit the requested walltime and `--rewrite` to replace existing case output. + +After all jobs have completed, process the latest Maestro run: + +```bash +python process.py +``` + +Pass a Maestro run directory to process a specific run: + +```bash +python process.py maestro_run_ +``` + +Processed figures are written to this suite's `results/` directory. +The top-level `process.py` collects these figures under the repository's `results/` directory. + +Each case's `plot.py` can inspect one result interactively, showing both the energy spectrum and the cycle-by-cycle multiplication factor: + +```bash +python cases/inf_shem361_subcritical/plot.py cases/inf_shem361_subcritical/output_320.h5 +``` + +## Cases + +| Case | Description | +| :--- | :---------- | +| [`inf_shem361_subcritical`](cases/inf_shem361_subcritical/) | Infinite homogeneous SHEM-361 with the capture cross section increased by 50%. | +| [`inf_shem361_supercritical`](cases/inf_shem361_supercritical/) | Infinite homogeneous SHEM-361 with the original capture cross section. | + +For each case, the analytical multiplication factor and spectrum are obtained from + +$$ +\left[\operatorname{diag}(\Sigma_t)-\Sigma_s\right]\phi += \frac{1}{k}\nu\Sigma_f\phi. +$$ + +SciPy solves the generalized eigenvalue problem, and the dominant real eigenpair supplies $k$ and the normalized fundamental-mode spectrum. + +## References + +- A. Hébert and A. Santamarina, *Refinement of the Santamarina-Hfaiedh Energy Mesh Between 22.5 eV and 11.4 keV*, PHYSOR 2008, vol. 2, pp. 929–938, 2008. diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/input.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/input.py new file mode 100644 index 0000000..0060099 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/input.py @@ -0,0 +1,74 @@ +from pathlib import Path + +import numpy as np +import mcdc + +SHEM361_DATA = ( + Path(__file__).resolve().parents[3] / "fixed_source" / "data" / "SHEM-361.npz" +) + +CAPTURE_FACTOR = 1.5 +N_PARTICLE = 1000000 +N_INACTIVE = 20 +N_ACTIVE = 20 + +simulation = mcdc.Simulation("Infinite SHEM-361 subcritical k-eigenvalue") + +# ====================================================================================== +# Set model +# ====================================================================================== +# The infinite homogeneous medium is modeled with a reflecting slab. + +with np.load(SHEM361_DATA) as data: + SigmaC = data["SigmaC"] * CAPTURE_FACTOR + SigmaS = data["SigmaS"] + SigmaF = data["SigmaF"] + nu_p = data["nu_p"] + nu_d = data["nu_d"] + chi_p = data["chi_p"] + chi_d = data["chi_d"] + G = int(data["G"]) + +material = mcdc.Material.multigroup( + capture=SigmaC, + scatter=SigmaS, + fission=SigmaF, + nu_p=nu_p, + chi_p=chi_p, + nu_d=nu_d, + chi_d=chi_d, +) + +left = mcdc.Surface.PlaneX(x=-1.0e10, boundary_condition="reflective") +right = mcdc.Surface.PlaneX(x=1.0e10, boundary_condition="reflective") +cell = mcdc.Cell(region=+left & -right, fill=material) +simulation.set_model([cell]) + +# ====================================================================================== +# Set initial source +# ====================================================================================== + +source = mcdc.Source( + position=(0.0, 0.0, 0.0), + isotropic=True, + discrete_energy=np.array([[G - 1], [1.0]]), +) +simulation.set_sources([source]) + +# ====================================================================================== +# Set tallies, settings, techniques, and run MC/DC +# ====================================================================================== + +tally = mcdc.Tally(scores=["flux"], energy="all") +simulation.set_tallies([tally]) + +# Particle and inactive-cycle counts stay fixed while the study overrides N_active. +simulation.settings.N_particle = N_PARTICLE +simulation.settings.set_eigenmode(N_inactive=N_INACTIVE, N_active=N_ACTIVE) +simulation.settings.census_bank_buffer_ratio = 4.0 +simulation.settings.source_bank_buffer_ratio = 4.0 + +# Keep exactly N_particle source histories after each eigenvalue cycle. +simulation.technique.population_control() + +simulation.run() diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/plot.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/plot.py new file mode 100644 index 0000000..79a3197 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/plot.py @@ -0,0 +1,88 @@ +import sys +from pathlib import Path + +import h5py +import matplotlib.pyplot as plt +import numpy as np + +from reference import reference + +SHEM361_DATA = ( + Path(__file__).resolve().parents[3] / "fixed_source" / "data" / "SHEM-361.npz" +) + +# Load the energy grid and one MC/DC result. +output_file = sys.argv[1] +with np.load(SHEM361_DATA) as data: + energy = data["E"] + energy_midpoint = 0.5 * (energy[1:] + energy[:-1]) + energy_width = np.diff(energy) + +with h5py.File(output_file, "r") as output: + flux = output["tallies/tracklength_tally_0/flux/mean"][:] + flux_sdev = output["tallies/tracklength_tally_0/flux/sdev"][:] + k_cycle = output["k_cycle"][:] + k_mean = output["k_mean"][()] + k_sdev = output["k_sdev"][()] + N_inactive = int(output["settings/N_inactive"][()]) + +# Normalize the arbitrary eigenvector amplitude before comparing spectrum shapes. +k_reference, flux_reference = reference() +normalization = np.sum(flux) +flux /= normalization +flux_sdev /= normalization + +flux = flux / energy_width * energy_midpoint +flux_sdev = flux_sdev / energy_width * energy_midpoint +flux_reference = flux_reference / energy_width * energy_midpoint + +# Compare the MC/DC and analytical fundamental-mode energy spectra. +plt.plot(energy_midpoint, flux, "-b", label="MC/DC") +plt.fill_between( + energy_midpoint, + flux - flux_sdev, + flux + flux_sdev, + alpha=0.2, + color="b", +) +plt.plot(energy_midpoint, flux_reference, "--r", label="analytical") +plt.xscale("log") +plt.xlabel(r"$E$, eV") +plt.ylabel(r"$E\phi(E)$") +plt.grid() +plt.legend() +plt.show() + +# Show convergence of the cycle estimates and distinguish active statistics. +cycles = np.arange(1, len(k_cycle) + 1) +active_cycles = cycles[N_inactive:] + +plt.plot(cycles, k_cycle, "-b", label="MC/DC cycle") +plt.axhline(k_reference, color="k", linestyle="--", label="analytical") +plt.plot( + active_cycles, + np.full(active_cycles.size, k_mean), + ":r", + label="MC/DC active mean", +) +plt.fill_between( + active_cycles, + k_mean - k_sdev, + k_mean + k_sdev, + alpha=0.2, + color="r", +) + +if N_inactive > 0: + plt.axvline( + N_inactive + 0.5, + color="0.5", + linestyle=":", + label="active cycles begin", + ) + +plt.xlabel("Cycle") +plt.ylabel(r"$k$") +plt.grid() +plt.legend() +plt.show() diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/process.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/process.py new file mode 100644 index 0000000..9d09d97 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/process.py @@ -0,0 +1,43 @@ +import sys +from pathlib import Path + +import h5py +import numpy as np + +from reference import reference + +SUITE_DIR = Path(__file__).resolve().parents[2] +if str(SUITE_DIR) not in sys.path: + sys.path.insert(0, str(SUITE_DIR)) + +import util + +# Reproduce the active-cycle counts used by run_case.py. +N_active_min = int(sys.argv[1]) +N_active_max = int(sys.argv[2]) +N_task = int(sys.argv[3]) +active_cycle_counts = np.rint( + np.geomspace( + N_active_min, + N_active_max, + N_task, + ) +).astype(int) + +# Calculate the dominant analytical eigenpair once for all tasks. +k_reference, flux_reference = reference() +k_error = np.zeros(N_task) +flux_error = np.zeros(N_task) + +for index, N_active in enumerate(active_cycle_counts): + with h5py.File(f"output_{N_active}.h5", "r") as output: + k_effective = output["k_mean"][()] + flux = output["tallies/tracklength_tally_0/flux/mean"][:] + + # Only the spectrum shape matters for this infinite-medium eigenproblem. + flux /= np.sum(flux) + k_error[index] = abs(k_effective - k_reference) / k_reference + flux_error[index] = util.relative_error(flux, flux_reference) + +util.plot_convergence("k-effective", active_cycle_counts, k_error) +util.plot_convergence("flux", active_cycle_counts, flux_error) diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/reference.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/reference.py new file mode 100644 index 0000000..aba82fc --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_subcritical/reference.py @@ -0,0 +1,30 @@ +from pathlib import Path + +import numpy as np +from scipy.linalg import eig + +SHEM361_DATA = ( + Path(__file__).resolve().parents[3] / "fixed_source" / "data" / "SHEM-361.npz" +) +CAPTURE_FACTOR = 1.5 + + +def reference(): + """Return the dominant multiplication factor and normalized spectrum.""" + with np.load(SHEM361_DATA) as data: + total = data["SigmaT"] + (CAPTURE_FACTOR - 1.0) * data["SigmaC"] + loss = np.diag(total) - data["SigmaS"] + production = data["nuSigmaF"] + + eigenvalues, eigenvectors = eig(production, loss) + physical = np.isfinite(eigenvalues) & (np.abs(eigenvalues.imag) < 1.0e-10) + indices = np.flatnonzero(physical) + dominant = indices[np.argmax(eigenvalues[physical].real)] + + k_effective = float(eigenvalues[dominant].real) + flux = eigenvectors[:, dominant].real + if np.sum(flux) < 0.0: + flux *= -1.0 + flux /= np.sum(flux) + + return k_effective, flux diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/input.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/input.py new file mode 100644 index 0000000..7557528 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/input.py @@ -0,0 +1,74 @@ +from pathlib import Path + +import numpy as np +import mcdc + +SHEM361_DATA = ( + Path(__file__).resolve().parents[3] / "fixed_source" / "data" / "SHEM-361.npz" +) + +CAPTURE_FACTOR = 1.0 +N_PARTICLE = 1000000 +N_INACTIVE = 20 +N_ACTIVE = 20 + +simulation = mcdc.Simulation("Infinite SHEM-361 supercritical k-eigenvalue") + +# ====================================================================================== +# Set model +# ====================================================================================== +# The infinite homogeneous medium is modeled with a reflecting slab. + +with np.load(SHEM361_DATA) as data: + SigmaC = data["SigmaC"] * CAPTURE_FACTOR + SigmaS = data["SigmaS"] + SigmaF = data["SigmaF"] + nu_p = data["nu_p"] + nu_d = data["nu_d"] + chi_p = data["chi_p"] + chi_d = data["chi_d"] + G = int(data["G"]) + +material = mcdc.Material.multigroup( + capture=SigmaC, + scatter=SigmaS, + fission=SigmaF, + nu_p=nu_p, + chi_p=chi_p, + nu_d=nu_d, + chi_d=chi_d, +) + +left = mcdc.Surface.PlaneX(x=-1.0e10, boundary_condition="reflective") +right = mcdc.Surface.PlaneX(x=1.0e10, boundary_condition="reflective") +cell = mcdc.Cell(region=+left & -right, fill=material) +simulation.set_model([cell]) + +# ====================================================================================== +# Set initial source +# ====================================================================================== + +source = mcdc.Source( + position=(0.0, 0.0, 0.0), + isotropic=True, + discrete_energy=np.array([[G - 1], [1.0]]), +) +simulation.set_sources([source]) + +# ====================================================================================== +# Set tallies, settings, techniques, and run MC/DC +# ====================================================================================== + +tally = mcdc.Tally(scores=["flux"], energy="all") +simulation.set_tallies([tally]) + +# Particle and inactive-cycle counts stay fixed while the study overrides N_active. +simulation.settings.N_particle = N_PARTICLE +simulation.settings.set_eigenmode(N_inactive=N_INACTIVE, N_active=N_ACTIVE) +simulation.settings.census_bank_buffer_ratio = 4.0 +simulation.settings.source_bank_buffer_ratio = 4.0 + +# Keep exactly N_particle source histories after each eigenvalue cycle. +simulation.technique.population_control() + +simulation.run() diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/plot.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/plot.py new file mode 100644 index 0000000..79a3197 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/plot.py @@ -0,0 +1,88 @@ +import sys +from pathlib import Path + +import h5py +import matplotlib.pyplot as plt +import numpy as np + +from reference import reference + +SHEM361_DATA = ( + Path(__file__).resolve().parents[3] / "fixed_source" / "data" / "SHEM-361.npz" +) + +# Load the energy grid and one MC/DC result. +output_file = sys.argv[1] +with np.load(SHEM361_DATA) as data: + energy = data["E"] + energy_midpoint = 0.5 * (energy[1:] + energy[:-1]) + energy_width = np.diff(energy) + +with h5py.File(output_file, "r") as output: + flux = output["tallies/tracklength_tally_0/flux/mean"][:] + flux_sdev = output["tallies/tracklength_tally_0/flux/sdev"][:] + k_cycle = output["k_cycle"][:] + k_mean = output["k_mean"][()] + k_sdev = output["k_sdev"][()] + N_inactive = int(output["settings/N_inactive"][()]) + +# Normalize the arbitrary eigenvector amplitude before comparing spectrum shapes. +k_reference, flux_reference = reference() +normalization = np.sum(flux) +flux /= normalization +flux_sdev /= normalization + +flux = flux / energy_width * energy_midpoint +flux_sdev = flux_sdev / energy_width * energy_midpoint +flux_reference = flux_reference / energy_width * energy_midpoint + +# Compare the MC/DC and analytical fundamental-mode energy spectra. +plt.plot(energy_midpoint, flux, "-b", label="MC/DC") +plt.fill_between( + energy_midpoint, + flux - flux_sdev, + flux + flux_sdev, + alpha=0.2, + color="b", +) +plt.plot(energy_midpoint, flux_reference, "--r", label="analytical") +plt.xscale("log") +plt.xlabel(r"$E$, eV") +plt.ylabel(r"$E\phi(E)$") +plt.grid() +plt.legend() +plt.show() + +# Show convergence of the cycle estimates and distinguish active statistics. +cycles = np.arange(1, len(k_cycle) + 1) +active_cycles = cycles[N_inactive:] + +plt.plot(cycles, k_cycle, "-b", label="MC/DC cycle") +plt.axhline(k_reference, color="k", linestyle="--", label="analytical") +plt.plot( + active_cycles, + np.full(active_cycles.size, k_mean), + ":r", + label="MC/DC active mean", +) +plt.fill_between( + active_cycles, + k_mean - k_sdev, + k_mean + k_sdev, + alpha=0.2, + color="r", +) + +if N_inactive > 0: + plt.axvline( + N_inactive + 0.5, + color="0.5", + linestyle=":", + label="active cycles begin", + ) + +plt.xlabel("Cycle") +plt.ylabel(r"$k$") +plt.grid() +plt.legend() +plt.show() diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/process.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/process.py new file mode 100644 index 0000000..9d09d97 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/process.py @@ -0,0 +1,43 @@ +import sys +from pathlib import Path + +import h5py +import numpy as np + +from reference import reference + +SUITE_DIR = Path(__file__).resolve().parents[2] +if str(SUITE_DIR) not in sys.path: + sys.path.insert(0, str(SUITE_DIR)) + +import util + +# Reproduce the active-cycle counts used by run_case.py. +N_active_min = int(sys.argv[1]) +N_active_max = int(sys.argv[2]) +N_task = int(sys.argv[3]) +active_cycle_counts = np.rint( + np.geomspace( + N_active_min, + N_active_max, + N_task, + ) +).astype(int) + +# Calculate the dominant analytical eigenpair once for all tasks. +k_reference, flux_reference = reference() +k_error = np.zeros(N_task) +flux_error = np.zeros(N_task) + +for index, N_active in enumerate(active_cycle_counts): + with h5py.File(f"output_{N_active}.h5", "r") as output: + k_effective = output["k_mean"][()] + flux = output["tallies/tracklength_tally_0/flux/mean"][:] + + # Only the spectrum shape matters for this infinite-medium eigenproblem. + flux /= np.sum(flux) + k_error[index] = abs(k_effective - k_reference) / k_reference + flux_error[index] = util.relative_error(flux, flux_reference) + +util.plot_convergence("k-effective", active_cycle_counts, k_error) +util.plot_convergence("flux", active_cycle_counts, flux_error) diff --git a/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/reference.py b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/reference.py new file mode 100644 index 0000000..f2868d0 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cases/inf_shem361_supercritical/reference.py @@ -0,0 +1,30 @@ +from pathlib import Path + +import numpy as np +from scipy.linalg import eig + +SHEM361_DATA = ( + Path(__file__).resolve().parents[3] / "fixed_source" / "data" / "SHEM-361.npz" +) +CAPTURE_FACTOR = 1.0 + + +def reference(): + """Return the dominant multiplication factor and normalized spectrum.""" + with np.load(SHEM361_DATA) as data: + total = data["SigmaT"] + (CAPTURE_FACTOR - 1.0) * data["SigmaC"] + loss = np.diag(total) - data["SigmaS"] + production = data["nuSigmaF"] + + eigenvalues, eigenvectors = eig(production, loss) + physical = np.isfinite(eigenvalues) & (np.abs(eigenvalues.imag) < 1.0e-10) + indices = np.flatnonzero(physical) + dominant = indices[np.argmax(eigenvalues[physical].real)] + + k_effective = float(eigenvalues[dominant].real) + flux = eigenvectors[:, dominant].real + if np.sum(flux) < 0.0: + flux *= -1.0 + flux /= np.sum(flux) + + return k_effective, flux diff --git a/verification/analytical/neutron/k_eigenvalue/cleanup.py b/verification/analytical/neutron/k_eigenvalue/cleanup.py new file mode 100644 index 0000000..963c1f8 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/cleanup.py @@ -0,0 +1,13 @@ +"""Remove generated case outputs and processed figures from this suite.""" + +from pathlib import Path + +suite_dir = Path(__file__).resolve().parent + +for output in (suite_dir / "cases").glob("*/output*.h5"): + output.unlink() + +results_dir = suite_dir / "results" +if results_dir.is_dir(): + for figure in results_dir.glob("*.png"): + figure.unlink() diff --git a/verification/analytical/neutron/k_eigenvalue/launch.py b/verification/analytical/neutron/k_eigenvalue/launch.py new file mode 100644 index 0000000..b1fc3d8 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/launch.py @@ -0,0 +1,230 @@ +"""Build and launch the analytical k-eigenvalue verification study.""" + +import argparse +import os +import subprocess +import sys +from pathlib import Path + +import yaml + +# ====================================================================================== +# Bootstrap VVP imports +# ====================================================================================== + +REPO_DIR = Path(__file__).resolve().parents[4] + +# Support launching this suite directly without installing MC/DC-VVP. +if str(REPO_DIR) not in sys.path: + sys.path.insert(0, str(REPO_DIR)) + + +# ====================================================================================== +# Load shared VVP configs +# ====================================================================================== + +from configs.platform_config import PLATFORMS + +# User overrides are optional for local runs. +try: + from configs.user_config import USER_CONFIG +except ImportError: + USER_CONFIG = {} + + +# ====================================================================================== +# Command-line arguments +# ====================================================================================== + +parser = argparse.ArgumentParser( + description="Launch the MC/DC VVP analytical neutron k-eigenvalue suite." +) +parser.add_argument("--platform", default="local", choices=["local"] + list(PLATFORMS)) +parser.add_argument("--mpi", action="store_true") +parser.add_argument("--walltime", type=int, default=None) +parser.add_argument("--rewrite", action="store_true") +args = parser.parse_args() + +if args.mpi and args.platform == "local": + parser.error("--mpi requires a cluster platform. Specify --platform .") + + +# ====================================================================================== +# Paths and platform settings +# ====================================================================================== + +suite_dir = Path(__file__).resolve().parent +task_file = suite_dir / "task.yaml" +run_case = suite_dir / "run_case.py" +study_file = suite_dir / "study.yaml" + +local = args.platform == "local" +user_platform_config = USER_CONFIG.get(args.platform, {}) + +mcdc_python = user_platform_config.get("mcdc_python") +if mcdc_python is None: + mcdc_python = sys.executable +else: + mcdc_python = str(Path(mcdc_python).expanduser()) + +if not local: + platform = PLATFORMS[args.platform] + scheduler = platform["scheduler"] + cpu_cores = platform["cpu_cores_per_node"] + + # Never request more walltime than the platform permits. + walltime_hours = ( + platform["max_walltime_hours"] + if args.walltime is None + else min(args.walltime, platform["max_walltime_hours"]) + ) + walltime = platform["walltime_format"].format(hours=walltime_hours) + + account = user_platform_config.get("account") + queue = user_platform_config.get("queue") + reservation = user_platform_config.get("reservation") + + if account is None: + raise ValueError( + f"Platform '{args.platform}' requires an account. " + "Create configs/user_config.py from configs/user_config.py.template." + ) + + +# ====================================================================================== +# Build Maestro study +# ====================================================================================== + +with task_file.open("r") as f: + tasks = yaml.safe_load(f) + +steps = [] + +for case_name in tasks: + safe_case_name = case_name.replace("-", "_") + command = f"{mcdc_python} {run_case} --name {case_name}" + + # Maestro replaces LAUNCHER with the scheduler-specific MPI launch command. + if args.mpi: + command += ' --mpi "$(LAUNCHER)"' + + if args.rewrite: + command += " --rewrite" + + run = {"cmd": command} + + # Scheduled studies require explicit resources; local studies run directly. + if not local: + run["nodes"] = 1 + run["walltime"] = walltime + if args.mpi: + run["procs"] = cpu_cores + run["exclusive"] = True + else: + run["procs"] = 1 + + steps.append( + { + "name": safe_case_name, + "description": f"Run case: {case_name}", + "run": run, + } + ) + +study = { + "description": { + "name": "maestro_run", + "description": "MC/DC verification - analytical neutron - k-eigenvalue suite", + }, + "env": {"variables": {}}, + "study": steps, +} + +if not local: + batch = { + "type": scheduler, + "host": platform["host"], + "bank": account, + } + if queue is not None: + batch["queue"] = queue + if reservation is not None: + batch["reservation"] = reservation + study["batch"] = batch + +with study_file.open("w") as f: + yaml.dump(study, f, sort_keys=False) + + +# ====================================================================================== +# Launch Maestro +# ====================================================================================== + +maestro_python = None +if not local: + maestro_python = user_platform_config.get("maestro_python") + +env = os.environ.copy() +if maestro_python is None: + maestro_command = ["maestro", "run", "study.yaml"] +else: + maestro_python = Path(maestro_python).expanduser() + maestro_bin = maestro_python.parent + env["PATH"] = f"{maestro_bin}:{env['PATH']}" + maestro_command = [ + str(maestro_python), + "-m", + "maestrowf.maestro", + "run", + "study.yaml", + ] + +subprocess.run(maestro_command, cwd=suite_dir, check=True, env=env) + + +# ====================================================================================== +# Store launch metadata +# ====================================================================================== + +maestro_runs = sorted( + suite_dir.glob("maestro_run_*"), + key=lambda path: path.stat().st_mtime, +) +if not maestro_runs: + raise RuntimeError("Maestro did not create a maestro_run_* directory.") + +latest_run = maestro_runs[-1] +launch_config = { + "platform": args.platform, + "mpi": args.mpi, + "walltime": args.walltime, + "rewrite": args.rewrite, +} + +with (latest_run / "launch_config.yaml").open("w") as f: + yaml.dump(launch_config, f, sort_keys=False) + +with (latest_run / "task.yaml").open("w") as f: + yaml.dump(tasks, f, sort_keys=False) + + +# ====================================================================================== +# Summary +# ====================================================================================== + +print(f"Platform : {args.platform}") +print(f"MPI : {args.mpi}") +print(f"Rewrite : {args.rewrite}") +print(f"Python : {mcdc_python}") +print(f"Study : {study_file}") + +if not local: + print(f"Scheduler: {scheduler}") + print(f"Account : {account}") + print(f"Queue : {queue}") + print(f"Reserv. : {reservation}") + print(f"Walltime : {walltime}") + print("Nodes : 1") + print(f"Procs : {cpu_cores if args.mpi else 1}") + +print(f"Cases : {len(steps)}") diff --git a/verification/analytical/neutron/k_eigenvalue/process.py b/verification/analytical/neutron/k_eigenvalue/process.py new file mode 100644 index 0000000..86d1b59 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/process.py @@ -0,0 +1,102 @@ +"""Process the analytical neutron k-eigenvalue verification suite.""" + +import argparse +import subprocess +import sys +from pathlib import Path + +import yaml + +# ====================================================================================== +# Command-line arguments +# ====================================================================================== + +parser = argparse.ArgumentParser( + description="Process the MC/DC VVP analytical neutron k-eigenvalue suite." +) +parser.add_argument( + "maestro_run", + nargs="?", + default=None, + help="Maestro run directory to process. Defaults to the latest maestro_run_*.", +) +args = parser.parse_args() + + +# ====================================================================================== +# Resolve the launch and load its configuration +# ====================================================================================== + +suite_dir = Path(__file__).resolve().parent + +if args.maestro_run is None: + maestro_runs = sorted( + suite_dir.glob("maestro_run_*"), + key=lambda path: path.stat().st_mtime, + ) + if not maestro_runs: + raise FileNotFoundError("No maestro_run_* directory found.") + maestro_run = maestro_runs[-1] +else: + maestro_run = Path(args.maestro_run).expanduser() + if not maestro_run.is_absolute(): + maestro_run = suite_dir / maestro_run + +launch_config_file = maestro_run / "launch_config.yaml" +task_file = maestro_run / "task.yaml" +if not launch_config_file.is_file(): + raise FileNotFoundError(f"Launch config not found: {launch_config_file}") +if not task_file.is_file(): + raise FileNotFoundError(f"Task config not found: {task_file}") + +with launch_config_file.open("r") as f: + launch_config = yaml.safe_load(f) + +with task_file.open("r") as f: + tasks = yaml.safe_load(f) + + +# ====================================================================================== +# Process cases and collect their figures +# ====================================================================================== + +results_dir = suite_dir / "results" +results_dir.mkdir(exist_ok=True) + +for case_name, task in tasks.items(): + case_dir = suite_dir / "cases" / case_name + process_script = case_dir / "process.py" + + if not process_script.is_file(): + print(f"Skipping {case_name}: no process.py") + continue + + print(f"Processing {case_name}") + subprocess.run( + [ + sys.executable, + str(process_script), + str(task["N_active_min"]), + str(task["N_active_max"]), + str(task["N_task"]), + ], + cwd=case_dir, + check=True, + ) + + for figure in case_dir.glob("*.png"): + destination = results_dir / f"{case_name}_{figure.name}" + figure.replace(destination) + + +# ====================================================================================== +# Summary +# ====================================================================================== + +print() +print(f"Maestro run: {maestro_run}") +print(f"Platform : {launch_config['platform']}") +print(f"MPI : {launch_config['mpi']}") +print(f"Rewrite : {launch_config['rewrite']}") +print(f"Cases : {len(tasks)}") +print("Processing complete.") diff --git a/verification/analytical/neutron/k_eigenvalue/run_case.py b/verification/analytical/neutron/k_eigenvalue/run_case.py new file mode 100644 index 0000000..17f5d89 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/run_case.py @@ -0,0 +1,95 @@ +"""Run one analytical k-eigenvalue case over its active-cycle study.""" + +import argparse +import subprocess +import sys +from pathlib import Path + +import numpy as np +import yaml + +# ====================================================================================== +# Command-line arguments +# ====================================================================================== + +parser = argparse.ArgumentParser( + description="Run one MC/DC VVP analytical k-eigenvalue verification case." +) +parser.add_argument("--name", required=True, help="Verification case name.") +parser.add_argument("--task-file", default="task.yaml") +parser.add_argument( + "--mpi", + default="", + help="MPI launch command supplied by Maestro, e.g. 'srun -n 112'.", +) +parser.add_argument( + "--rewrite", + action="store_true", + help="Rewrite existing output files instead of skipping them.", +) +args = parser.parse_args() + + +# ====================================================================================== +# Paths and task definition +# ====================================================================================== + +suite_dir = Path(__file__).resolve().parent +case_dir = suite_dir / "cases" / args.name +task_file = suite_dir / args.task_file + +if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + +with task_file.open("r") as f: + tasks = yaml.safe_load(f) + +if args.name not in tasks: + raise ValueError(f"Case '{args.name}' is not listed in {task_file}") + +task = tasks[args.name] + + +# ====================================================================================== +# Run active-cycle tasks +# ====================================================================================== + +active_cycle_counts = np.rint( + np.geomspace( + task["N_active_min"], + task["N_active_max"], + task["N_task"], + ) +).astype(int) + +for N_active in active_cycle_counts: + N_active = int(N_active) + + output = f"output_{N_active}" + output_file = case_dir / f"{output}.h5" + + if output_file.is_file(): + if args.rewrite: + print(f"Rewrite (remove existing output): {args.name}, N_active={N_active}") + output_file.unlink() + else: + print(f"Skip (output exists): {args.name}, N_active={N_active}") + continue + + command = ( + f"{args.mpi} {sys.executable} input.py " + "--mode=numba " + f"--N_active={N_active} " + f"--output={output} " + "--no-progress_bar " + "--caching" + ).strip() + + print("=" * 80) + print(f"Case : {args.name}") + print(f"N_active: {N_active}") + print(f"Python : {sys.executable}") + print(f"Command : {command}") + print("=" * 80) + + subprocess.run(command, shell=True, cwd=case_dir, check=True) diff --git a/verification/analytical/neutron/k_eigenvalue/task.yaml b/verification/analytical/neutron/k_eigenvalue/task.yaml new file mode 100644 index 0000000..8085d85 --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/task.yaml @@ -0,0 +1,26 @@ +# ============================================================================= +# MC/DC Verification Task Definitions +# ============================================================================= +# +# This file defines the analytical neutron k-eigenvalue cases. +# +# Each top-level key corresponds to a case directory under this suite. +# The active-cycle counts are generated with: +# +# numpy.rint( +# numpy.geomspace(N_active_min, N_active_max, N_task) +# ).astype(int) +# +# The fixed particle and inactive-cycle counts are set in each input.py. +# +# ============================================================================= + +inf_shem361_subcritical: + N_active_min: 20 + N_active_max: 2000 + N_task: 5 + +inf_shem361_supercritical: + N_active_min: 20 + N_active_max: 2000 + N_task: 5 diff --git a/verification/analytical/neutron/k_eigenvalue/util.py b/verification/analytical/neutron/k_eigenvalue/util.py new file mode 100644 index 0000000..f0e18af --- /dev/null +++ b/verification/analytical/neutron/k_eigenvalue/util.py @@ -0,0 +1,33 @@ +import matplotlib.pyplot as plt +import numpy as np + + +def relative_error(value, reference): + """Return the relative Euclidean error against a reference array.""" + return np.linalg.norm(value - reference) / np.linalg.norm(reference) + + +def plot_convergence(name, active_cycles, error): + """Plot observed error with an active-cycle inverse-square-root guide.""" + midpoint = len(active_cycles) // 2 + + plt.plot(active_cycles, error, "bo", fillstyle="none", label="MC/DC") + + expected = 1.0 / np.sqrt(active_cycles) + expected *= error[midpoint] / expected[midpoint] + plt.plot( + active_cycles, + expected, + "r--", + label=r"$O(N_\mathrm{active}^{-0.5})$", + ) + + plt.xscale("log") + plt.yscale("log") + plt.xlabel(r"Active cycles, $N_\mathrm{active}$") + plt.ylabel("Relative error") + plt.title(name) + plt.grid() + plt.legend() + plt.savefig(f"{name}.png") + plt.clf() diff --git a/verification/benchmark/neutron/continuous_energy/README.md b/verification/benchmark/neutron/continuous_energy/README.md new file mode 100644 index 0000000..c3203eb --- /dev/null +++ b/verification/benchmark/neutron/continuous_energy/README.md @@ -0,0 +1,14 @@ +# Continuous-energy verification + +The continuous-energy benchmarks require a nuclear-data library generated with the same MC/DC version used to run the cases. + +Set `MCDC_LIB` to the generated HDF5 library before launching a benchmark: + +```bash +export MCDC_LIB=/path/to/mcdc/library +``` + +Generate or refresh the library with MC/DC's `tools/data_library_generator/neutron/generate.py` utility. +The source ACE data location is configured with `MCDC_ACELIB`; see the generator's README in the MC/DC repository for installation and invocation details. + +The benchmark inputs perform a lightweight schema check before compiling so an outdated library produces an actionable error rather than a low-level HDF5 exception. diff --git a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/compare-spectrum.py b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/compare-spectrum.py index 9870a55..1c94e12 100644 --- a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/compare-spectrum.py +++ b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/compare-spectrum.py @@ -13,15 +13,15 @@ sdev_openmc = tally.std_dev.reshape((200, 500)) with h5py.File("mcdc/output.h5", "r") as f: - flux_mcdc = f["tallies/global_tally_0/flux/mean"][()].transpose() - sdev_mcdc = f["tallies/global_tally_0/flux/sdev"][()].transpose() + flux_mcdc = f["tallies/tracklength_tally_0/flux/mean"][()].transpose() + sdev_mcdc = f["tallies/tracklength_tally_0/flux/sdev"][()].transpose() # Grids - t = f["tallies/global_tally_0/grid/time"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] dt = t[1:] - t[:-1] t_mid = 0.5 * (t[1:] + t[:-1]) - E = f["tallies/global_tally_0/grid/energy"][()] + E = f["tallies/tracklength_tally_0/grid/energy"][()] E_mid = 0.5 * (E[1:] + E[:-1]) dE = E[1:] - E[:-1] diff --git a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/mcdc/input.py b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/mcdc/input.py index 303fc0d..876b325 100644 --- a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/mcdc/input.py +++ b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-helium/mcdc/input.py @@ -1,6 +1,8 @@ import mcdc import numpy as np +simulation = mcdc.Simulation("Pulsed UO2-helium pincell") + # Material fuel = mcdc.Material( nuclide_composition={ @@ -19,17 +21,19 @@ y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective") y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective") # -mcdc.Cell(-cylinder, fill=fuel) -mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator) +fuel_cell = mcdc.Cell(-cylinder, fill=fuel) +moderator_cell = mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator) +simulation.set_model([fuel_cell, moderator_cell]) # Source -mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6) +source = mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6) +simulation.set_sources([source]) # Setting -mcdc.settings.N_particle = 10000 -mcdc.settings.N_batch = 30 -mcdc.settings.time_boundary = 1.0 -mcdc.settings.active_bank_buffer = 1000 +simulation.settings.N_particle = 10000 +simulation.settings.N_batch = 30 +simulation.settings.time_boundary = 1.0 +simulation.settings.active_bank_buffer = 1000 # Tally t_grid = np.insert(np.logspace(-9, -4, 200), 0, 0.0) @@ -37,6 +41,7 @@ groups = 500 energies = np.logspace(np.log10(e_min), np.log10(e_max), groups + 1) -mcdc.TallyGlobal(scores=["flux"], time=t_grid, energy=energies) +tally = mcdc.Tally(scores=["flux"], time=t_grid, energy=energies) +simulation.set_tallies([tally]) -mcdc.run() +simulation.run() diff --git a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/compare-spectrum.py b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/compare-spectrum.py index 6ec33ef..03981fd 100644 --- a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/compare-spectrum.py +++ b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/compare-spectrum.py @@ -13,15 +13,15 @@ sdev_openmc = tally.std_dev.reshape((200, 500)) with h5py.File("mcdc/output.h5", "r") as f: - flux_mcdc = f["tallies/global_tally_0/flux/mean"][()].transpose() - sdev_mcdc = f["tallies/global_tally_0/flux/sdev"][()].transpose() + flux_mcdc = f["tallies/tracklength_tally_0/flux/mean"][()].transpose() + sdev_mcdc = f["tallies/tracklength_tally_0/flux/sdev"][()].transpose() # Grids - t = f["tallies/global_tally_0/grid/time"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] dt = t[1:] - t[:-1] t_mid = 0.5 * (t[1:] + t[:-1]) - E = f["tallies/global_tally_0/grid/energy"][()] + E = f["tallies/tracklength_tally_0/grid/energy"][()] E_mid = 0.5 * (E[1:] + E[:-1]) dE = E[1:] - E[:-1] diff --git a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/mcdc/input.py b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/mcdc/input.py index 8f9b5e9..9cbcc77 100644 --- a/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/mcdc/input.py +++ b/verification/benchmark/neutron/continuous_energy/pulsed_pincells/uo2-water/mcdc/input.py @@ -1,6 +1,8 @@ import mcdc import numpy as np +simulation = mcdc.Simulation("Pulsed UO2-water pincell") + # Material fuel = mcdc.Material( nuclide_composition={ @@ -25,17 +27,19 @@ y0 = mcdc.Surface.PlaneY(y=-pitch / 2, boundary_condition="reflective") y1 = mcdc.Surface.PlaneY(y=pitch / 2, boundary_condition="reflective") # -mcdc.Cell(-cylinder, fill=fuel) -mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator) +fuel_cell = mcdc.Cell(-cylinder, fill=fuel) +moderator_cell = mcdc.Cell(+x0 & -x1 & +y0 & -y1 & +cylinder, fill=moderator) +simulation.set_model([fuel_cell, moderator_cell]) # Source -mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6) +source = mcdc.Source(position=[0.0, 0.0, 0.0], isotropic=True, time=0.0, energy=14.1e6) +simulation.set_sources([source]) # Setting -mcdc.settings.N_particle = 10000 -mcdc.settings.N_batch = 30 -mcdc.settings.time_boundary = 1.0 -mcdc.settings.active_bank_buffer = 1000 +simulation.settings.N_particle = 10000 +simulation.settings.N_batch = 30 +simulation.settings.time_boundary = 1.0 +simulation.settings.active_bank_buffer = 1000 # Tally t_grid = np.insert(np.logspace(-9, -4, 200), 0, 0.0) @@ -43,6 +47,7 @@ groups = 500 energies = np.logspace(np.log10(e_min), np.log10(e_max), groups + 1) -mcdc.TallyGlobal(scores=["flux"], time=t_grid, energy=energies) +tally = mcdc.Tally(scores=["flux"], time=t_grid, energy=energies) +simulation.set_tallies([tally]) -mcdc.run() +simulation.run() diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/input.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/input.py index 16c7592..fb2e5f9 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/input.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/input.py @@ -1,5 +1,11 @@ +import math + +import h5py import numpy as np -import h5py, mcdc, math + +import mcdc + +simulation = mcdc.Simulation("C5G7 four-phase transient") # ============================================================================= # Materials @@ -11,7 +17,7 @@ # Setter def set_mat(mat): - return mcdc.MaterialMG( + return mcdc.Material.multigroup( capture=mat["capture"][:], scatter=mat["scatter"][:], fission=mat["fission"][:], @@ -313,9 +319,9 @@ def set_mat(mat): reflector_south = mcdc.Cell(+x0 & -x3 & +y0 & -y1 & +z1 & -z2, mat_mod) reflector_east = mcdc.Cell(+x2 & -x3 & +y1 & -y3 & +z1 & -z2, mat_mod) -# Root universe -mcdc.simulation.set_root_universe( - cells=[ +# Set model +simulation.set_model( + [ assembly_1, assembly_2, assembly_3, @@ -323,7 +329,7 @@ def set_mat(mat): reflector_bottom, reflector_south, reflector_east, - ], + ] ) # ============================================================================= @@ -336,9 +342,11 @@ def set_mat(mat): x=np.array([pitch * 17 * 3 / 2] * 2) + np.array([-pitch / 2, +pitch / 2]), y=np.array([-pitch * 17 * 3 / 2] * 2) + np.array([-pitch / 2, +pitch / 2]), z=[-core_height / 2, core_height / 2], - energy_group=0, + isotropic=True, + energy=0, time=[0.0, 15.0], ) +simulation.set_sources([source]) # ============================================================================= # Set tallies, settings, and run MC/DC @@ -354,12 +362,13 @@ def set_mat(mat): y = np.linspace(-pitch * 17 * 2, 0.0, Ny + 1) z = np.linspace(-core_height / 2, core_height / 2, Nz + 1) mesh = mcdc.MeshStructured(x=x, y=y, z=z) -mcdc.Tally(mesh=mesh, scores=["fission"], time=t) +tally = mcdc.Tally(mesh=mesh, scores=["fission"], time=t) +simulation.set_tallies([tally]) # Settings -mcdc.settings.N_particle = int(1e6) -mcdc.settings.N_batch = 30 -mcdc.settings.active_bank_buffer = 10000 +simulation.settings.N_particle = int(1e6) +simulation.settings.N_batch = 30 +simulation.settings.active_bank_buffer = 10000 # Run -mcdc.run() +simulation.run() diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission-sdev.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission-sdev.py index f2f4338..fa0a82b 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission-sdev.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission-sdev.py @@ -7,13 +7,13 @@ # Get fission rates with h5py.File("output_4.h5", "r") as f: - fissions = f["tallies/mesh_tally_0/fission/mean"][()] - fissions_sd = f["tallies/mesh_tally_0/fission/sdev"][()] + fissions = f["tallies/tracklength_tally_0/fission/mean"][()] + fissions_sd = f["tallies/tracklength_tally_0/fission/sdev"][()] - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/t"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] # The grids t_mid = 0.5 * (t[:-1] + t[1:]) diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission.py index 755442b..3c63e92 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-fission.py @@ -7,11 +7,11 @@ # Get results with h5py.File("output_4.h5", "r") as f: - fissions = f["tallies/mesh_tally_0/fission/mean"][()] - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/t"][()] + fissions = f["tallies/tracklength_tally_0/fission/mean"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] # Total fission fission_total = np.average(fissions, axis=(1, 2, 3)) diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast-sdev.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast-sdev.py index c056f8f..fa48188 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast-sdev.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast-sdev.py @@ -7,12 +7,12 @@ # Get results with h5py.File("output.h5", "r") as f: - fission_total = f["tallies/mesh_tally_0/fission/sdev"][()] - fluxes = f["tallies/mesh_tally_2/flux/sdev"][:, 0, :, :] - x = f["tallies/mesh_tally_2/grid/x"][()] - y = f["tallies/mesh_tally_2/grid/y"][()] - z = f["tallies/mesh_tally_2/grid/z"][()] - t = f["tallies/mesh_tally_2/grid/t"][()] + fission_total = f["tallies/tracklength_tally_0/fission/sdev"][()] + fluxes = f["tallies/tracklength_tally_2/flux/sdev"][:, 0, :, :] + x = f["tallies/tracklength_tally_2/grid/x"][()] + y = f["tallies/tracklength_tally_2/grid/y"][()] + z = f["tallies/tracklength_tally_2/grid/z"][()] + t = f["tallies/tracklength_tally_2/grid/time"][()] # The grids t_mid = 0.5 * (t[:-1] + t[1:]) diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast.py index 3a98722..5d82b4f 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-fast.py @@ -7,12 +7,12 @@ # Get results with h5py.File("output.h5", "r") as f: - fission_total = f["tallies/mesh_tally_0/fission/mean"][()] - fluxes = f["tallies/mesh_tally_2/flux/mean"][:, 0, :, :] - x = f["tallies/mesh_tally_2/grid/x"][()] - y = f["tallies/mesh_tally_2/grid/y"][()] - z = f["tallies/mesh_tally_2/grid/z"][()] - t = f["tallies/mesh_tally_2/grid/t"][()] + fission_total = f["tallies/tracklength_tally_0/fission/mean"][()] + fluxes = f["tallies/tracklength_tally_2/flux/mean"][:, 0, :, :] + x = f["tallies/tracklength_tally_2/grid/x"][()] + y = f["tallies/tracklength_tally_2/grid/y"][()] + z = f["tallies/tracklength_tally_2/grid/z"][()] + t = f["tallies/tracklength_tally_2/grid/time"][()] # The grids t_mid = 0.5 * (t[:-1] + t[1:]) diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal-sdev.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal-sdev.py index d5e27cf..1c2313d 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal-sdev.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal-sdev.py @@ -7,12 +7,12 @@ # Get results with h5py.File("output.h5", "r") as f: - fission_total = f["tallies/mesh_tally_0/fission/sdev"][()] - fluxes = f["tallies/mesh_tally_2/flux/sdev"][:, 0, :, :] - x = f["tallies/mesh_tally_2/grid/x"][()] - y = f["tallies/mesh_tally_2/grid/y"][()] - z = f["tallies/mesh_tally_2/grid/z"][()] - t = f["tallies/mesh_tally_2/grid/t"][()] + fission_total = f["tallies/tracklength_tally_0/fission/sdev"][()] + fluxes = f["tallies/tracklength_tally_2/flux/sdev"][:, 0, :, :] + x = f["tallies/tracklength_tally_2/grid/x"][()] + y = f["tallies/tracklength_tally_2/grid/y"][()] + z = f["tallies/tracklength_tally_2/grid/z"][()] + t = f["tallies/tracklength_tally_2/grid/time"][()] # The grids t_mid = 0.5 * (t[:-1] + t[1:]) diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal.py index cc99a65..90bdb18 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/plot-flux-thermal.py @@ -7,12 +7,12 @@ # Get results with h5py.File("output.h5", "r") as f: - fission_total = f["tallies/mesh_tally_0/fission/mean"][()] - fluxes = f["tallies/mesh_tally_2/flux/mean"][:, 1, :, :] - x = f["tallies/mesh_tally_2/grid/x"][()] - y = f["tallies/mesh_tally_2/grid/y"][()] - z = f["tallies/mesh_tally_2/grid/z"][()] - t = f["tallies/mesh_tally_2/grid/t"][()] + fission_total = f["tallies/tracklength_tally_0/fission/mean"][()] + fluxes = f["tallies/tracklength_tally_2/flux/mean"][:, 1, :, :] + x = f["tallies/tracklength_tally_2/grid/x"][()] + y = f["tallies/tracklength_tally_2/grid/y"][()] + z = f["tallies/tracklength_tally_2/grid/z"][()] + t = f["tallies/tracklength_tally_2/grid/time"][()] # The grids t_mid = 0.5 * (t[:-1] + t[1:]) diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/process.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/process.py deleted file mode 100644 index 90aa070..0000000 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/mcdc/process.py +++ /dev/null @@ -1,4 +0,0 @@ -import mcdc -import sys - -mcdc.recombine_tallies(sys.argv[1]) diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-convergence.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-convergence.py index d095fcf..a50c4f7 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-convergence.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-convergence.py @@ -18,7 +18,7 @@ # Getting the reference with h5py.File("mcdc/output_4.h5", "r") as f: - fission_mcdc = f["tallies/mesh_tally_0/fission/mean"][()] + fission_mcdc = f["tallies/tracklength_tally_0/fission/mean"][()] with openmc.StatePoint("openmc_/output_4.h5") as sp: tally = sp.get_tally(name="pincell fission") @@ -36,7 +36,7 @@ for n in range(NN): # Get results with h5py.File("mcdc/output_%i.h5" % n, "r") as f: - fission_mcdc = f["tallies/mesh_tally_0/fission/mean"][()] + fission_mcdc = f["tallies/tracklength_tally_0/fission/mean"][()] # Get results with openmc.StatePoint("openmc_/output_%i.h5" % n) as sp: diff --git a/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-difference.py b/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-difference.py index a8a429c..649ebfe 100644 --- a/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-difference.py +++ b/verification/benchmark/neutron/multigroup/c5g7-4phase/plot-difference.py @@ -10,11 +10,11 @@ # Get results with h5py.File("mcdc/output_4.h5", "r") as f: - fissions_mcdc = f["tallies/mesh_tally_0/fission/mean"][()] - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/t"][()] + fissions_mcdc = f["tallies/tracklength_tally_0/fission/mean"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] # Get results Nt = 200 diff --git a/verification/benchmark/neutron/multigroup/kobayashi/mcdc/input.py b/verification/benchmark/neutron/multigroup/kobayashi/mcdc/input.py index bd57da0..83ae74b 100644 --- a/verification/benchmark/neutron/multigroup/kobayashi/mcdc/input.py +++ b/verification/benchmark/neutron/multigroup/kobayashi/mcdc/input.py @@ -1,6 +1,8 @@ import numpy as np import mcdc +simulation = mcdc.Simulation("Kobayashi dog-leg transient benchmark") + # ====================================================================================== # Set model # ====================================================================================== @@ -8,8 +10,8 @@ # (PNE 2001, https://doi.org/10.1016/S0149-1970(01)00007-5) # Set materials -m = mcdc.MaterialMG(capture=np.array([0.05]), scatter=np.array([[0.05]])) -m_void = mcdc.MaterialMG(capture=np.array([5e-5]), scatter=np.array([[5e-5]])) +m = mcdc.Material.multigroup(capture=np.array([0.05]), scatter=np.array([[0.05]])) +m_void = mcdc.Material.multigroup(capture=np.array([5e-5]), scatter=np.array([[5e-5]])) # Set surfaces sx1 = mcdc.Surface.PlaneX(x=0.0, boundary_condition="reflective") @@ -41,20 +43,22 @@ # Shield box = +sx1 & -sx5 & +sy1 & -sy5 & +sz1 & -sz5 shield_cell = mcdc.Cell(region=box & ~void_channel, fill=m) +simulation.set_model([source_cell, void_cell, shield_cell]) # ====================================================================================== # Set source # ====================================================================================== # The source pulses in t=[0,5] -mcdc.Source( +source = mcdc.Source( x=[0.0, 10.0], y=[0.0, 10.0], z=[0.0, 10.0], isotropic=True, - energy_group=0, + energy=0, time=[0.0, 50.0], ) +simulation.set_sources([source]) # ====================================================================================== # Set tallies, settings, and run MC/DC @@ -62,15 +66,16 @@ # Tallies time_grid = np.linspace(0.0, 200.0, 101) -#mesh = mcdc.MeshUniform(x=(0.0, 1.0, 60), y=(0.0, 1.0, 100), z=(0.0, 1.0, 60)) +# mesh = mcdc.MeshUniform(x=(0.0, 1.0, 60), y=(0.0, 1.0, 100), z=(0.0, 1.0, 60)) mesh = mcdc.MeshUniform(x=(0.0, 1.0, 60)) -mcdc.Tally(mesh=mesh, scores=["flux"], time=time_grid) -mcdc.Tally(scores=["density"], time=time_grid) +flux_tally = mcdc.Tally(mesh=mesh, scores=["flux"], time=time_grid) +density_tally = mcdc.Tally(scores=["density"], time=time_grid) +simulation.set_tallies([flux_tally, density_tally]) # Settings -#mcdc.settings.N_particle = int(1e9) -mcdc.settings.N_particle = int(1e2) -mcdc.settings.N_batch = 30 +# simulation.settings.N_particle = int(1e9) +simulation.settings.N_particle = int(1e2) +simulation.settings.N_batch = 30 # Run -mcdc.run() +simulation.run() diff --git a/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-flux.py b/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-flux.py index aa429bf..bd3ee6a 100644 --- a/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-flux.py +++ b/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-flux.py @@ -7,12 +7,12 @@ # Get results with h5py.File("output_4.h5", "r") as f: - fluxes = f["tallies/mesh_tally_0/flux/mean"][()] - densities = f["tallies/global_tally_0/density/mean"][()] - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/time"][()] + fluxes = f["tallies/tracklength_tally_0/flux/mean"][()] + densities = f["tallies/tracklength_tally_0/density/mean"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] # The grids t_mid = 0.5 * (t[:-1] + t[1:]) diff --git a/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-sdev.py b/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-sdev.py index 3e52727..98c80ab 100644 --- a/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-sdev.py +++ b/verification/benchmark/neutron/multigroup/kobayashi/mcdc/plot-sdev.py @@ -9,14 +9,14 @@ # Get results with h5py.File("output_4.h5", "r") as f: - fluxes = f["tallies/mesh_tally_0/flux/mean"][()] - densities = f["tallies/global_tally_0/density/mean"][()] - fluxes_sd = f["tallies/mesh_tally_0/flux/sdev"][()] - densities_sd = f["tallies/global_tally_0/density/sdev"][()] - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/t"][()] + fluxes = f["tallies/tracklength_tally_0/flux/mean"][()] + densities = f["tallies/tracklength_tally_0/density/mean"][()] + fluxes_sd = f["tallies/tracklength_tally_0/flux/sdev"][()] + densities_sd = f["tallies/tracklength_tally_0/density/sdev"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] # The grids t_mid = 0.5 * (t[:-1] + t[1:]) diff --git a/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-flux.py b/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-flux.py index a681d32..cb50675 100644 --- a/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-flux.py +++ b/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-flux.py @@ -16,10 +16,10 @@ # The grids with h5py.File("../mcdc/output.h5", "r") as f: - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/t"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] t_mid = 0.5 * (t[:-1] + t[1:]) XY_X, XY_Y = np.meshgrid(x, y, indexing="ij") XZ_X, XZ_Z = np.meshgrid(x, z, indexing="ij") diff --git a/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-sdev.py b/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-sdev.py index f4b47c0..0374730 100644 --- a/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-sdev.py +++ b/verification/benchmark/neutron/multigroup/kobayashi/openmc_/plot-sdev.py @@ -19,10 +19,10 @@ # The grids with h5py.File("../mcdc/output.h5", "r") as f: - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/t"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] t_mid = 0.5 * (t[:-1] + t[1:]) XY_X, XY_Y = np.meshgrid(x, y, indexing="ij") XZ_X, XZ_Z = np.meshgrid(x, z, indexing="ij") diff --git a/verification/benchmark/neutron/multigroup/kobayashi/plot-convergence.py b/verification/benchmark/neutron/multigroup/kobayashi/plot-convergence.py index 5fe549f..23aacde 100644 --- a/verification/benchmark/neutron/multigroup/kobayashi/plot-convergence.py +++ b/verification/benchmark/neutron/multigroup/kobayashi/plot-convergence.py @@ -13,7 +13,7 @@ # Getting the reference with h5py.File("mcdc/output_4.h5", "r") as f: - flux_mcdc = f["tallies/mesh_tally_0/flux/mean"][()] + flux_mcdc = f["tallies/tracklength_tally_0/flux/mean"][()] with openmc.StatePoint("openmc_/output_4.h5") as sp: tally = sp.get_tally(scores=["flux"]) flux_openmc = tally.mean.reshape((100, 60, 100, 60)) @@ -29,7 +29,7 @@ for n in range(NN): # Get results with h5py.File("mcdc/output_%i.h5" % n, "r") as f: - flux_mcdc = f["tallies/mesh_tally_0/flux/mean"][()] + flux_mcdc = f["tallies/tracklength_tally_0/flux/mean"][()] # Get results with openmc.StatePoint("openmc_/output_%i.h5" % n) as sp: diff --git a/verification/benchmark/neutron/multigroup/kobayashi/plot-difference.py b/verification/benchmark/neutron/multigroup/kobayashi/plot-difference.py index 58d655d..4a07163 100644 --- a/verification/benchmark/neutron/multigroup/kobayashi/plot-difference.py +++ b/verification/benchmark/neutron/multigroup/kobayashi/plot-difference.py @@ -10,12 +10,12 @@ # Get results with h5py.File("mcdc/output_4.h5", "r") as f: - fluxes_mcdc = f["tallies/mesh_tally_0/flux/mean"][()] - densities_mcdc = f["tallies/mesh_tally_1/density/mean"][()] - x = f["tallies/mesh_tally_0/grid/x"][()] - y = f["tallies/mesh_tally_0/grid/y"][()] - z = f["tallies/mesh_tally_0/grid/z"][()] - t = f["tallies/mesh_tally_0/grid/t"][()] + fluxes_mcdc = f["tallies/tracklength_tally_0/flux/mean"][()] + densities_mcdc = f["tallies/tracklength_tally_1/density/mean"][()] + x = f["tallies/tracklength_tally_0/grid/x"][()] + y = f["tallies/tracklength_tally_0/grid/y"][()] + z = f["tallies/tracklength_tally_0/grid/z"][()] + t = f["tallies/tracklength_tally_0/grid/time"][()] # Get results with openmc.StatePoint("openmc_/output_4.h5") as sp: