diff --git a/.gitignore b/.gitignore index 704b9fb04..2447591aa 100644 --- a/.gitignore +++ b/.gitignore @@ -33,6 +33,9 @@ output*.h5 # Ignore all png and giffiles except those in the docs/source/images directory !docs/source/images/**/*.png !docs/source/images/**/*.gif +!mcdc-vv/Al/results/*.png +!mcdc-vv/Al/results/old/*.png + # Cluster job *.out @@ -68,4 +71,5 @@ source_particles.h5 # Container artifacts *.sif *_sandbox -*.tar \ No newline at end of file +*.tar + diff --git a/mcdc-vv/Al/E00050keV/th000deg/input.py b/mcdc-vv/Al/E00050keV/th000deg/input.py new file mode 100644 index 000000000..b6ec55af8 --- /dev/null +++ b/mcdc-vv/Al/E00050keV/th000deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 1e5 # eV +CSDA_RANGE = 0.00604 # g/cm2 +ANGLE = 0.0 + +# MCDC Simulation Parameters +N_PARTICLES = 1000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() diff --git a/mcdc-vv/Al/E00050keV/th000deg/lw_Al_5e+04eV_1e3p_20260418_1528.h5 b/mcdc-vv/Al/E00050keV/th000deg/lw_Al_5e+04eV_1e3p_20260418_1528.h5 new file mode 100644 index 000000000..243b69474 Binary files /dev/null and b/mcdc-vv/Al/E00050keV/th000deg/lw_Al_5e+04eV_1e3p_20260418_1528.h5 differ diff --git a/mcdc-vv/Al/E00050keV/th000deg/process.py b/mcdc-vv/Al/E00050keV/th000deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E00050keV/th000deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E00050keV/th000deg/reference.npz b/mcdc-vv/Al/E00050keV/th000deg/reference.npz new file mode 100644 index 000000000..0b63943df Binary files /dev/null and b/mcdc-vv/Al/E00050keV/th000deg/reference.npz differ diff --git a/mcdc-vv/Al/E00100keV/th000deg/input.py b/mcdc-vv/Al/E00100keV/th000deg/input.py new file mode 100644 index 000000000..961a6db7d --- /dev/null +++ b/mcdc-vv/Al/E00100keV/th000deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 1e5 # eV +CSDA_RANGE = 0.0203 # g/cm2 +ANGLE = 0.0 + +# MCDC Simulation Parameters +N_PARTICLES = 1000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/mcdc-vv/Al/E00100keV/th000deg/lw_Al_1e+05eV_1e3p_20260418_1529.h5 b/mcdc-vv/Al/E00100keV/th000deg/lw_Al_1e+05eV_1e3p_20260418_1529.h5 new file mode 100644 index 000000000..e8c394df0 Binary files /dev/null and b/mcdc-vv/Al/E00100keV/th000deg/lw_Al_1e+05eV_1e3p_20260418_1529.h5 differ diff --git a/mcdc-vv/Al/E00100keV/th000deg/process.py b/mcdc-vv/Al/E00100keV/th000deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E00100keV/th000deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E00100keV/th000deg/reference.npz b/mcdc-vv/Al/E00100keV/th000deg/reference.npz new file mode 100644 index 000000000..30738e86b Binary files /dev/null and b/mcdc-vv/Al/E00100keV/th000deg/reference.npz differ diff --git a/mcdc-vv/Al/E00300keV/th000deg/input.py b/mcdc-vv/Al/E00300keV/th000deg/input.py new file mode 100644 index 000000000..6138ebc14 --- /dev/null +++ b/mcdc-vv/Al/E00300keV/th000deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 3e5 # eV +CSDA_RANGE = 0.113 # g/cm2 +ANGLE = 0.0 + +# MCDC Simulation Parameters +N_PARTICLES = 1000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/mcdc-vv/Al/E00300keV/th000deg/lw_Al_3e+05eV_1e3p_20260418_1530.h5 b/mcdc-vv/Al/E00300keV/th000deg/lw_Al_3e+05eV_1e3p_20260418_1530.h5 new file mode 100644 index 000000000..0fcd6b7bc Binary files /dev/null and b/mcdc-vv/Al/E00300keV/th000deg/lw_Al_3e+05eV_1e3p_20260418_1530.h5 differ diff --git a/mcdc-vv/Al/E00300keV/th000deg/process.py b/mcdc-vv/Al/E00300keV/th000deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E00300keV/th000deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E00300keV/th000deg/reference.npz b/mcdc-vv/Al/E00300keV/th000deg/reference.npz new file mode 100644 index 000000000..08e6d02da Binary files /dev/null and b/mcdc-vv/Al/E00300keV/th000deg/reference.npz differ diff --git a/mcdc-vv/Al/E00300keV/th060deg/input.py b/mcdc-vv/Al/E00300keV/th060deg/input.py new file mode 100644 index 000000000..5eb27daa3 --- /dev/null +++ b/mcdc-vv/Al/E00300keV/th060deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 3e5 # eV +CSDA_RANGE = 0.113 # g/cm2 +ANGLE = 60.0 + +# MCDC Simulation Parameters +N_PARTICLES = 1000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/mcdc-vv/Al/E00300keV/th060deg/lw_Al_3e+05eV_1e3p_20260418_1531.h5 b/mcdc-vv/Al/E00300keV/th060deg/lw_Al_3e+05eV_1e3p_20260418_1531.h5 new file mode 100644 index 000000000..b8b501482 Binary files /dev/null and b/mcdc-vv/Al/E00300keV/th060deg/lw_Al_3e+05eV_1e3p_20260418_1531.h5 differ diff --git a/mcdc-vv/Al/E00300keV/th060deg/process.py b/mcdc-vv/Al/E00300keV/th060deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E00300keV/th060deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E00300keV/th060deg/reference.npz b/mcdc-vv/Al/E00300keV/th060deg/reference.npz new file mode 100644 index 000000000..95b6be040 Binary files /dev/null and b/mcdc-vv/Al/E00300keV/th060deg/reference.npz differ diff --git a/mcdc-vv/Al/E00500keV/th000deg/input.py b/mcdc-vv/Al/E00500keV/th000deg/input.py new file mode 100644 index 000000000..c4bd9c84e --- /dev/null +++ b/mcdc-vv/Al/E00500keV/th000deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 5e5 # eV +CSDA_RANGE = 0.234 # g/cm2 +ANGLE = 0.0 + +# MCDC Simulation Parameters +N_PARTICLES = 1000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1328.h5 b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1328.h5 new file mode 100644 index 000000000..523f28720 Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1328.h5 differ diff --git a/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1344.h5 b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1344.h5 new file mode 100644 index 000000000..b99d596ca Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1344.h5 differ diff --git a/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1347.h5 b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1347.h5 new file mode 100644 index 000000000..1e089cb85 Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1347.h5 differ diff --git a/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1526.h5 b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1526.h5 new file mode 100644 index 000000000..be509edea Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1526.h5 differ diff --git a/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1532.h5 b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1532.h5 new file mode 100644 index 000000000..de32bbd86 Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e3p_20260418_1532.h5 differ diff --git a/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e4p_20260418_1349.h5 b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e4p_20260418_1349.h5 new file mode 100644 index 000000000..e9baaf473 Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th000deg/lw_Al_5e+05eV_1e4p_20260418_1349.h5 differ diff --git a/mcdc-vv/Al/E00500keV/th000deg/process.py b/mcdc-vv/Al/E00500keV/th000deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E00500keV/th000deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E00500keV/th000deg/reference.npz b/mcdc-vv/Al/E00500keV/th000deg/reference.npz new file mode 100644 index 000000000..3cc934e5d Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th000deg/reference.npz differ diff --git a/mcdc-vv/Al/E00500keV/th060deg/input.py b/mcdc-vv/Al/E00500keV/th060deg/input.py new file mode 100644 index 000000000..104c8ab2a --- /dev/null +++ b/mcdc-vv/Al/E00500keV/th060deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 5e5 # eV +CSDA_RANGE = 0.234 # g/cm2 +ANGLE = 60.0 + +# MCDC Simulation Parameters +N_PARTICLES = 1000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/mcdc-vv/Al/E00500keV/th060deg/lw_Al_5e+05eV_1e3p_20260418_1533.h5 b/mcdc-vv/Al/E00500keV/th060deg/lw_Al_5e+05eV_1e3p_20260418_1533.h5 new file mode 100644 index 000000000..2477f2682 Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th060deg/lw_Al_5e+05eV_1e3p_20260418_1533.h5 differ diff --git a/mcdc-vv/Al/E00500keV/th060deg/process.py b/mcdc-vv/Al/E00500keV/th060deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E00500keV/th060deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E00500keV/th060deg/reference.npz b/mcdc-vv/Al/E00500keV/th060deg/reference.npz new file mode 100644 index 000000000..3772b6f71 Binary files /dev/null and b/mcdc-vv/Al/E00500keV/th060deg/reference.npz differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/answer.h5 b/mcdc-vv/Al/E01000keV/th000deg/answer.h5 new file mode 100644 index 000000000..a10234351 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/answer.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/input.py b/mcdc-vv/Al/E01000keV/th000deg/input.py new file mode 100644 index 000000000..26bd3c7d8 --- /dev/null +++ b/mcdc-vv/Al/E01000keV/th000deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 1e6 # eV +CSDA_RANGE = 0.569 # g/cm2 +ANGLE = 0.0 + +# MCDC Simulation Parameters +N_PARTICLES = 10000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1424.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1424.h5 new file mode 100644 index 000000000..06b1a2a93 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1424.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1459.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1459.h5 new file mode 100644 index 000000000..86c9ddf2d Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1459.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1502.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1502.h5 new file mode 100644 index 000000000..1f9a88f05 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1502.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1510.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1510.h5 new file mode 100644 index 000000000..0a2a2f101 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1510.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1512.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1512.h5 new file mode 100644 index 000000000..951869961 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e3p_20260418_1512.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1450.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1450.h5 new file mode 100644 index 000000000..3f1d35ac6 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1450.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1515.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1515.h5 new file mode 100644 index 000000000..114c9a63b Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1515.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1534.h5 b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1534.h5 new file mode 100644 index 000000000..a457ee256 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/lw_Al_1e+06eV_1e4p_20260418_1534.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th000deg/process.py b/mcdc-vv/Al/E01000keV/th000deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E01000keV/th000deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E01000keV/th000deg/reference.npz b/mcdc-vv/Al/E01000keV/th000deg/reference.npz new file mode 100644 index 000000000..79d69223e Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th000deg/reference.npz differ diff --git a/mcdc-vv/Al/E01000keV/th060deg/input.py b/mcdc-vv/Al/E01000keV/th060deg/input.py new file mode 100644 index 000000000..87866865d --- /dev/null +++ b/mcdc-vv/Al/E01000keV/th060deg/input.py @@ -0,0 +1,119 @@ +import numpy as np +import os +import math +import mcdc +from datetime import datetime + +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +DATA_LIBRARY_DIR = os.path.abspath( + os.path.join( + os.path.dirname(__file__), + "..", + "..", + "..", + "..", + "..", + "..", + "..", + "electron-vv-data", + ) +) + +# Set the XS library directory +os.environ["MCDC_LIB"] = os.environ.get(PROCESS_DATA_LIBRARY_ENV, DATA_LIBRARY_DIR) + +# ============================================================================= +# Set problem parameters +# ============================================================================= +# Energy and Angle Parameters +MATERIAL_SYMBOL = "Al" +ENERGY = 1e6 # eV +CSDA_RANGE = 0.569 # g/cm2 +ANGLE = 60.0 + +# MCDC Simulation Parameters +N_PARTICLES = 1000 +z0 = 0.0 # Starting source position + +# Material Properties +RHO_G_CM3 = 2.70 # g/cm3 +ATOMIC_WEIGHT_G_MOL = 26.7497084 # g/mol +AREAL_DENSITY_G_CM2 = 5.05e-3 #g/cm2 + +# Standard Calculations +dz = AREAL_DENSITY_G_CM2 / RHO_G_CM3 +AVAGADRO_NUMBER = 6.02214076e23 # atoms/mol +MAT_DENSITY_ATOMS_PER_BARN_CM = AVAGADRO_NUMBER / ATOMIC_WEIGHT_G_MOL * RHO_G_CM3 / 1e24 # atoms/barn-cm +TINY = 1e-30 +L = CSDA_RANGE / RHO_G_CM3 # cm +N_LAYERS = int(L / dz) +THETA = math.radians(ANGLE) + +# Output variables for naming +np_name = len(str(N_PARTICLES)) - 1 +e_name = f"{ENERGY:.2g}" + +# ============================================================================= +# Debug prints +# ============================================================================= +print(f"[DEBUG] Material Symbol = {MATERIAL_SYMBOL}") +print(f"[DEBUG] CSDA Range = {CSDA_RANGE:.6e} g/cm2") +print(f"[DEBUG] Density = {RHO_G_CM3:.6e} g/cm3") +print(f"[DEBUG] Atomic Weight = {ATOMIC_WEIGHT_G_MOL:.6e} g/mol") +print(f"[DEBUG] Areal density per layer = {AREAL_DENSITY_G_CM2:.6e} g/cm2") +print(f"[DEBUG] Material density = {MAT_DENSITY_ATOMS_PER_BARN_CM:.6e} atoms/barn-cm") +print(f"[DEBUG] Layer thickness = {dz:.6e} cm") +print(f"[DEBUG] Total thickness = {L:.6e} cm") +print(f"[DEBUG] Number of layers = {N_LAYERS}") + +# ============================================================================= +# Set materials +# ============================================================================= +mat = mcdc.Material( + element_composition={MATERIAL_SYMBOL: MAT_DENSITY_ATOMS_PER_BARN_CM} +) + +# ============================================================================= +# Set geometry (surfaces and cells) +# ============================================================================= +# Z-direction surfaces for layers + +s1 = mcdc.Surface.PlaneZ(z=0.0, boundary_condition="vacuum") +s2 = mcdc.Surface.PlaneZ(z=L, boundary_condition="vacuum") + +mcdc.Cell(region=+s1 & -s2, fill=mat) + +# ============================================================================= +# Set source +# ============================================================================= +# Parallel beam of 1 MeV electrons entering at z=0 + +mcdc.Source( + z=[z0 + TINY, z0 + TINY], + particle_type='electron', + energy=np.array([[ENERGY - 1, ENERGY + 1], [0.5, 0.5]]), + direction=[math.sin(THETA), 0.0+TINY, math.cos(THETA)] +) + +# ============================================================================= +# Set tally +# ============================================================================= +# Energy deposition tally along z-axis +z_bins = np.linspace(0.0, L, N_LAYERS + 1) +mesh = mcdc.MeshStructured(z=z_bins) + +mcdc.Tally(name="edep", mesh=mesh, scores=["energy_deposition"]) +mcdc.Tally(name="flux", scores=["flux"], mesh=mesh) +mcdc.Tally(name="s1_current", surface=s1, scores=["net-current"]) +mcdc.Tally(name="s2_current", surface=s2, scores=["net-current"]) + +# ============================================================================= +# Settings and run +# ============================================================================= +mcdc.settings.set_transported_particles(["electron"]) +mcdc.settings.N_particle = N_PARTICLES +mcdc.settings.active_bank_buffer = N_PARTICLES * 100 +mcdc.settings.output_name = f"lw_{MATERIAL_SYMBOL}_{e_name}eV_1e{np_name}p_{datetime.now():%Y%m%d_%H%M%S}" +mcdc.settings.use_progress_bar = True + +mcdc.run() \ No newline at end of file diff --git a/mcdc-vv/Al/E01000keV/th060deg/lw_Al_1e+06eV_1e3p_20260418_1545.h5 b/mcdc-vv/Al/E01000keV/th060deg/lw_Al_1e+06eV_1e3p_20260418_1545.h5 new file mode 100644 index 000000000..1f68f00d0 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th060deg/lw_Al_1e+06eV_1e3p_20260418_1545.h5 differ diff --git a/mcdc-vv/Al/E01000keV/th060deg/process.py b/mcdc-vv/Al/E01000keV/th060deg/process.py new file mode 100644 index 000000000..8f71e15a9 --- /dev/null +++ b/mcdc-vv/Al/E01000keV/th060deg/process.py @@ -0,0 +1,319 @@ +import argparse +import os +import re +import subprocess +import sys +from pathlib import Path + + +DEFAULT_DATA_LIBRARY_NAME = "electron-vv-data" +PROCESS_DATA_LIBRARY_DIR = None +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" + + +def _find_repo_root(start_path): + path = Path(start_path).resolve() + for parent in (path.parent, *path.parents): + if parent.name == "MCDC-VV-electron": + return parent + raise FileNotFoundError("Could not locate the MCDC-VV-electron repository root.") + + +def resolve_process_data_library_dir(process_path=__file__): + repo_root = _find_repo_root(process_path) + + if PROCESS_DATA_LIBRARY_DIR is None: + return str((repo_root / DEFAULT_DATA_LIBRARY_NAME).resolve()) + + candidate = Path(PROCESS_DATA_LIBRARY_DIR).expanduser() + if not candidate.is_absolute(): + candidate = repo_root / candidate + return str(candidate.resolve()) + + +def _format_energy_dir(energy_keV): + return f"E{int(round(float(energy_keV))):05d}keV" + + +def _format_angle_dir(angle_deg): + return f"th{int(round(float(angle_deg))):03d}deg" + + +def resolve_case_directory(problem, material, energy_keV, angle_deg, process_path=__file__): + repo_root = _find_repo_root(process_path) + case_dir = ( + repo_root + / "verification" + / "benchmark" + / "continuous_energy" + / problem + / material + / _format_energy_dir(energy_keV) + / _format_angle_dir(angle_deg) + ) + if not case_dir.is_dir(): + raise FileNotFoundError(f"Case directory not found: {case_dir}") + return case_dir + + +def _find_input_value(name, text): + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case_dir): + input_path = Path(case_dir) / "input.py" + if not input_path.exists(): + raise FileNotFoundError(f"input.py not found in {case_dir}") + + text = input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "energy_eV": float(eval(_find_input_value("ENERGY", text), {}, {})), + "angle_deg": float(eval(_find_input_value("ANGLE", text), {}, {})), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": float(eval(_find_input_value("N_PARTICLES", text), {}, {})), + } + + +def resolve_h5_output(case_dir): + answer_path = Path(case_dir) / "answer.h5" + if answer_path.exists(): + return answer_path + + candidates = [ + path + for path in Path(case_dir).iterdir() + if path.is_file() and path.suffix == ".h5" + ] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + return max(candidates, key=lambda path: path.stat().st_mtime) + + +def run_case_input(case_dir, process_path=__file__): + data_library_dir = resolve_process_data_library_dir(process_path) + env = os.environ.copy() + env[PROCESS_DATA_LIBRARY_ENV] = data_library_dir + env["MCDC_LIB"] = data_library_dir + + print(f"Running input: {Path(case_dir) / 'input.py'}") + print(f"Using process data library: {data_library_dir}") + subprocess.run([sys.executable, "input.py"], cwd=case_dir, check=True, env=env) + + +def process_case(case_dir): + import h5py + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np + + case_dir = Path(case_dir).resolve() + ref_path = case_dir / "reference.npz" + if not ref_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case_dir}") + + h5_path = resolve_h5_output(case_dir) + if h5_path.name != "answer.h5": + print(f"Using HDF5 output: {h5_path.name}") + + metadata = load_case_metadata(case_dir) + + ref = np.load(ref_path) + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + + exp_fmr_A = exp_edep_A = None + exp_fmr_B = exp_edep_B = None + exp_fmr = exp_edep = None + + if "fmr_exp_lw_A" in ref: + exp_fmr_A = ref["fmr_exp_lw_A"].astype(float) + exp_edep_A = ref["edep_exp_lw_A"].astype(float) + + if "fmr_exp_lw_B" in ref: + exp_fmr_B = ref["fmr_exp_lw_B"].astype(float) + exp_edep_B = ref["edep_exp_lw_B"].astype(float) + else: + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = file["tallies/edep/grid/z"][:].astype(float) + dz = z[1:] - z[:-1] + edep = np.atleast_1d(file["tallies/edep/energy_deposition/mean"][()]).astype( + float + ) + + L = metadata["csda_range"] / metadata["rho_g_cm3"] + z_centers = 0.5 * (z[:-1] + z[1:]) + mcdc_fmr = z_centers / L + edep_mcdc = edep / metadata["rho_g_cm3"] / dz / 1e6 + + plt.figure(figsize=(14, 9), constrained_layout=True) + mpl.rcParams.update( + { + "font.size": 20, + "axes.titlesize": 29, + "axes.labelsize": 25, + "legend.fontsize": 22, + "xtick.labelsize": 25, + "ytick.labelsize": 25, + } + ) + + plt.minorticks_on() + plt.grid( + True, which="major", linestyle="-", linewidth=0.8, color="#3B3A3AFF", alpha=0.5 + ) + plt.grid( + True, which="minor", linestyle=":", linewidth=0.6, color="#8A8686", alpha=0.7 + ) + + plt.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=9, + linewidth=3.0, + ) + plt.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=9, + linewidth=3.0, + ) + + if exp_fmr_A is not None: + plt.plot( + exp_fmr_A, + exp_edep_A, + label="Lockwood Experimental A", + marker="D", + markersize=9, + linewidth=3.0, + ) + if exp_fmr_B is not None: + plt.plot( + exp_fmr_B, + exp_edep_B, + label="Lockwood Experimental B", + marker="*", + markersize=9, + linewidth=3.0, + ) + if exp_fmr is not None: + plt.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=9, + linewidth=3.0, + ) + + energy_MeV = metadata["energy_eV"] / 1e6 + angle_deg = metadata["angle_deg"] + + plt.xlim(0, 1) + plt.xlabel("Fraction of Mean Range", labelpad=15) + plt.ylabel("Energy Deposition (MeV/g/cm²)", labelpad=15) + plt.title( + f"Energy Deposition of {energy_MeV:g} MeV Electrons in {metadata['material_symbol']} at {angle_deg:g}° Incidence", + pad=20, + ) + plt.legend() + + out_dir = case_dir.parents[1] / "results" + out_dir.mkdir(exist_ok=True) + + out_name = ( + f"fig_{metadata['material_symbol']}_{energy_MeV:g}MeV_" + f"th{int(round(angle_deg))}_{metadata['n_particles']:.0e}".replace("+", "") + + ".png" + ) + out_path = out_dir / out_name + + plt.savefig(out_path, dpi=400, bbox_inches="tight") + print(f"Saved: {out_path}") + + +def build_parser(): + parser = argparse.ArgumentParser( + description="Run and/or post-process VV electron benchmark cases." + ) + parser.add_argument("-p", "--problem", help="Problem name, e.g. lockwood") + parser.add_argument("-m", "--material", help="Material name, e.g. Al") + parser.add_argument( + "-e", "--energy", type=float, help="Incident energy in keV, e.g. 1000" + ) + parser.add_argument( + "-a", "--angle", type=float, help="Incident angle in degrees, e.g. 30" + ) + parser.add_argument( + "--input-only", + action="store_true", + help="Run the selected input but skip post-processing.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip the input run and only post-process the selected case.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Print the resolved case and data paths without running anything.", + ) + return parser + + +def main(argv=None, process_path=__file__): + parser = build_parser() + args = parser.parse_args(argv) + + has_selector = any( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ) + if has_selector and not all( + value is not None for value in (args.problem, args.material, args.energy, args.angle) + ): + parser.error("Provide -p, -m, -e, and -a together.") + if args.input_only and args.process_only: + parser.error("--input-only and --process-only cannot be used together.") + if not has_selector and args.input_only: + parser.error("--input-only requires -p, -m, -e, and -a.") + + if has_selector: + case_dir = resolve_case_directory( + args.problem, args.material, args.energy, args.angle, process_path + ) + data_library_dir = resolve_process_data_library_dir(process_path) + + if args.dry_run: + print(f"Resolved case directory: {case_dir}") + print(f"Resolved process data library: {data_library_dir}") + return + + if not args.process_only: + run_case_input(case_dir, process_path) + if not args.input_only: + process_case(case_dir) + return + + if args.dry_run: + print(f"Resolved case directory: {Path.cwd().resolve()}") + print(f"Resolved process data library: {resolve_process_data_library_dir(process_path)}") + return + + process_case(Path.cwd()) + + +if __name__ == "__main__": + main() diff --git a/mcdc-vv/Al/E01000keV/th060deg/reference.npz b/mcdc-vv/Al/E01000keV/th060deg/reference.npz new file mode 100644 index 000000000..45f1ff650 Binary files /dev/null and b/mcdc-vv/Al/E01000keV/th060deg/reference.npz differ diff --git a/mcdc-vv/Al/results/fig_Al_0.05MeV_th0_1e03.png b/mcdc-vv/Al/results/fig_Al_0.05MeV_th0_1e03.png new file mode 100644 index 000000000..421ea2b77 Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_0.05MeV_th0_1e03.png differ diff --git a/mcdc-vv/Al/results/fig_Al_0.1MeV_th0_1e03.png b/mcdc-vv/Al/results/fig_Al_0.1MeV_th0_1e03.png new file mode 100644 index 000000000..5b513d88c Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_0.1MeV_th0_1e03.png differ diff --git a/mcdc-vv/Al/results/fig_Al_0.3MeV_th0_1e03.png b/mcdc-vv/Al/results/fig_Al_0.3MeV_th0_1e03.png new file mode 100644 index 000000000..b0d473ae4 Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_0.3MeV_th0_1e03.png differ diff --git a/mcdc-vv/Al/results/fig_Al_0.3MeV_th60_1e03.png b/mcdc-vv/Al/results/fig_Al_0.3MeV_th60_1e03.png new file mode 100644 index 000000000..119de2971 Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_0.3MeV_th60_1e03.png differ diff --git a/mcdc-vv/Al/results/fig_Al_0.5MeV_th0_1e03.png b/mcdc-vv/Al/results/fig_Al_0.5MeV_th0_1e03.png new file mode 100644 index 000000000..b79502737 Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_0.5MeV_th0_1e03.png differ diff --git a/mcdc-vv/Al/results/fig_Al_0.5MeV_th60_1e03.png b/mcdc-vv/Al/results/fig_Al_0.5MeV_th60_1e03.png new file mode 100644 index 000000000..c5be4d933 Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_0.5MeV_th60_1e03.png differ diff --git a/mcdc-vv/Al/results/fig_Al_1MeV_th0_1e04.png b/mcdc-vv/Al/results/fig_Al_1MeV_th0_1e04.png new file mode 100644 index 000000000..f95d1f6f5 Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_1MeV_th0_1e04.png differ diff --git a/mcdc-vv/Al/results/fig_Al_1MeV_th60_1e03.png b/mcdc-vv/Al/results/fig_Al_1MeV_th60_1e03.png new file mode 100644 index 000000000..89151f98a Binary files /dev/null and b/mcdc-vv/Al/results/fig_Al_1MeV_th60_1e03.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03-after-oley.png b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03-after-oley.png new file mode 100644 index 000000000..5c881fcfd Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03-after-oley.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03-beforeMCNP.png b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03-beforeMCNP.png new file mode 100644 index 000000000..64e20ebbf Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03-beforeMCNP.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03_mcdc-data.png b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03_mcdc-data.png new file mode 100644 index 000000000..d1891f29c Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e03_mcdc-data.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e04-beforeMCNP.png b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e04-beforeMCNP.png new file mode 100644 index 000000000..464af761f Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_0.5MeV_th0_1e04-beforeMCNP.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-2.png b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-2.png new file mode 100644 index 000000000..28ecf2082 Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-2.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-after.png b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-after.png new file mode 100644 index 000000000..14f899da5 Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-after.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-beforeMCNP.png b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-beforeMCNP.png new file mode 100644 index 000000000..fde65160b Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03-beforeMCNP.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03.png b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03.png new file mode 100644 index 000000000..28ecf2082 Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e03.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04-after.png b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04-after.png new file mode 100644 index 000000000..f95d1f6f5 Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04-after.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04-beforeMCNP.png b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04-beforeMCNP.png new file mode 100644 index 000000000..68ee05344 Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04-beforeMCNP.png differ diff --git a/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04.png b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04.png new file mode 100644 index 000000000..18d41ce09 Binary files /dev/null and b/mcdc-vv/Al/results/old/fig_Al_1MeV_th0_1e04.png differ diff --git a/mcdc-vv/run.py b/mcdc-vv/run.py new file mode 100644 index 000000000..04fa9cbab --- /dev/null +++ b/mcdc-vv/run.py @@ -0,0 +1,751 @@ +#!/usr/bin/env python3 + +from __future__ import annotations + +import argparse +import math +import os +import re +import shutil +import subprocess +import sys +import tempfile +from dataclasses import dataclass +from pathlib import Path + + +CASE_ROOT = Path(__file__).resolve().parent +SCRIPT_MCDC_ROOT = CASE_ROOT.parent.resolve() +WORKSPACE_ROOT = SCRIPT_MCDC_ROOT.parent +DEFAULT_MCDC_ROOT = (WORKSPACE_ROOT / "MCDC").resolve() +MCDC_ROOT = DEFAULT_MCDC_ROOT if DEFAULT_MCDC_ROOT.exists() else SCRIPT_MCDC_ROOT +DEFAULT_DATA_LIBRARY_DIR = (WORKSPACE_ROOT / "PyEEDL" / "mcdc_data").resolve() +PROCESS_DATA_LIBRARY_ENV = "MCDC_VV_PROCESS_DATA_LIBRARY_DIR" +ADDITIONAL_DATA_LIBRARY_ENVS = ("MCDC_VV_DATA_LIBRARY_DIR", "MCDC_LIB") +DEFAULT_RUN_MODE = "numba" +DEFAULT_TARGET = "cpu" +DEFAULT_MPI_PROCS = "auto" + + +@dataclass(frozen=True) +class CaseInfo: + material: str + energy_keV: int + angle_deg: int + case_dir: Path + + @property + def input_path(self) -> Path: + return self.case_dir / "input.py" + + @property + def reference_path(self) -> Path: + return self.case_dir / "reference.npz" + + @property + def energy_eV(self) -> float: + return float(self.energy_keV) * 1.0e3 + + @property + def results_dir(self) -> Path: + return CASE_ROOT / self.material / "results" + + @property + def label(self) -> str: + return f"{self.material} | {self.energy_keV} keV | {self.angle_deg} deg" + + +def parse_positive_int(value: str) -> int: + parsed = float(value) + if parsed <= 0 or not parsed.is_integer(): + raise argparse.ArgumentTypeError("Value must be a positive integer.") + return int(parsed) + + +def parse_mpi_procs(value: str) -> str | int: + text = value.strip().lower() + if text == "auto": + return "auto" + return parse_positive_int(value) + + +def format_scientific_int(value: int) -> str: + return f"{float(value):.0e}".replace("+", "") + + +def discover_cases() -> list[CaseInfo]: + cases: list[CaseInfo] = [] + pattern = re.compile(r"^E(\d+)keV$") + angle_pattern = re.compile(r"^th(\d+)deg$") + + for input_path in sorted(CASE_ROOT.rglob("input.py")): + rel_parts = input_path.relative_to(CASE_ROOT).parts + if len(rel_parts) != 4: + continue + + material, energy_dir, angle_dir, _ = rel_parts + energy_match = pattern.match(energy_dir) + angle_match = angle_pattern.match(angle_dir) + if not energy_match or not angle_match: + continue + + cases.append( + CaseInfo( + material=material, + energy_keV=int(energy_match.group(1)), + angle_deg=int(angle_match.group(1)), + case_dir=input_path.parent, + ) + ) + + return cases + + +def filter_cases( + cases: list[CaseInfo], + material: str | None = None, + energy_keV: float | None = None, + angle_deg: float | None = None, +) -> list[CaseInfo]: + selected = cases + + if material is not None: + selected = [case for case in selected if case.material.lower() == material.lower()] + + if energy_keV is not None: + selected = [ + case + for case in selected + if math.isclose(case.energy_keV, energy_keV, rel_tol=0.0, abs_tol=1.0e-9) + ] + + if angle_deg is not None: + selected = [ + case + for case in selected + if math.isclose(case.angle_deg, angle_deg, rel_tol=0.0, abs_tol=1.0e-9) + ] + + return selected + + +def _find_input_value(name: str, text: str) -> str: + match = re.search(rf"^\s*{name}\s*=\s*(.+?)\s*(#.*)?$", text, re.M) + if not match: + raise ValueError(f"{name} not found in input.py") + return match.group(1).strip() + + +def load_case_metadata(case: CaseInfo) -> dict[str, float | str]: + text = case.input_path.read_text(encoding="utf-8") + return { + "material_symbol": _find_input_value("MATERIAL_SYMBOL", text).strip("\"' "), + "csda_range": float(eval(_find_input_value("CSDA_RANGE", text), {}, {})), + "rho_g_cm3": float(eval(_find_input_value("RHO_G_CM3", text), {}, {})), + "n_particles": int(float(eval(_find_input_value("N_PARTICLES", text), {}, {}))), + } + + +def resolve_data_library_dir(cli_value: str | None) -> Path: + candidates: list[Path] = [] + + if cli_value: + candidates.append(Path(cli_value).expanduser()) + + candidates.append(DEFAULT_DATA_LIBRARY_DIR) + + for env_name in (PROCESS_DATA_LIBRARY_ENV, *ADDITIONAL_DATA_LIBRARY_ENVS): + env_value = os.environ.get(env_name) + if env_value: + candidates.append(Path(env_value).expanduser()) + + candidates.extend( + [ + MCDC_ROOT / "mcdc-regression_test_data", + ] + ) + + for candidate in candidates: + resolved = candidate.resolve() + if (resolved / "Al.h5").exists(): + return resolved + + tried = "\n".join(f" - {path.resolve()}" for path in candidates) or " - " + raise FileNotFoundError( + "Could not resolve the MCDC electron data library directory.\n" + "Set --data-library or one of the supported environment variables.\n" + f"Tried:\n{tried}" + ) + + +def build_runtime_env(data_library_dir: Path) -> dict[str, str]: + env = os.environ.copy() + repo_python_path = str(MCDC_ROOT) + current_python_path = env.get("PYTHONPATH") + env["PYTHONPATH"] = ( + f"{repo_python_path}{os.pathsep}{current_python_path}" + if current_python_path + else repo_python_path + ) + env[PROCESS_DATA_LIBRARY_ENV] = str(data_library_dir) + for env_name in ADDITIONAL_DATA_LIBRARY_ENVS: + env[env_name] = str(data_library_dir) + return env + + +def resolve_mpi_exec(cli_value: str | None) -> str | None: + candidates: list[str] = [] + + if cli_value: + candidates.append(cli_value) + + candidates.extend(["mpiexec", "mpirun"]) + + for candidate in candidates: + resolved = shutil.which(candidate) + if resolved: + return resolved + + return None + + +def resolve_mpi_procs(requested: str | int, particle_count: int, mpi_exec: str | None) -> int: + if requested == "auto": + if mpi_exec is None: + return 1 + + cpu_count = os.cpu_count() or 1 + return max(1, min(cpu_count, particle_count)) + + requested_int = int(requested) + if requested_int > 1 and mpi_exec is None: + raise FileNotFoundError( + "MPI execution was requested but neither `mpiexec` nor `mpirun` was found." + ) + return max(1, min(requested_int, particle_count)) + + +def build_mcdc_command( + python_executable: str, + input_name: str, + mode: str, + target: str, + caching: bool, + clear_cache: bool, + progress_bar: bool, + runtime_output: bool, + mpi_exec: str | None, + mpi_procs: int, +) -> list[str]: + command: list[str] = [] + + if mpi_procs > 1: + if mpi_exec is None: + raise FileNotFoundError( + "MPI execution was requested but neither `mpiexec` nor `mpirun` was found." + ) + command.extend([mpi_exec, "-n", str(mpi_procs)]) + + command.extend( + [ + python_executable, + input_name, + "--mode", + mode, + "--target", + target, + ] + ) + + command.append("--caching" if caching else "--no_caching") + if clear_cache: + command.append("--clear_cache") + if progress_bar: + command.append("--progress_bar") + else: + command.append("--no-progress_bar") + if runtime_output: + command.append("--runtime_output") + + return command + + +def replace_assignment(text: str, name: str, value: str) -> str: + pattern = re.compile(rf"^(\s*{name}\s*=\s*)(.+?)(\s*(#.*)?$)", re.M) + + def repl(match: re.Match[str]) -> str: + return f"{match.group(1)}{value}{match.group(3)}" + + updated, count = pattern.subn(repl, text, count=1) + if count != 1: + raise ValueError(f"Could not replace assignment for {name}") + return updated + + +def build_runtime_input(case: CaseInfo, particle_override: int | None) -> str: + text = case.input_path.read_text(encoding="utf-8") + text = replace_assignment(text, "MATERIAL_SYMBOL", repr(case.material)) + text = replace_assignment(text, "ENERGY", f"{case.energy_eV:.12g}") + text = replace_assignment(text, "ANGLE", f"{float(case.angle_deg):.12g}") + + if particle_override is not None: + text = replace_assignment(text, "N_PARTICLES", str(particle_override)) + + return text + + +def snapshot_h5_files(case_dir: Path) -> dict[str, int]: + return { + path.name: path.stat().st_mtime_ns + for path in case_dir.glob("*.h5") + if path.is_file() + } + + +def resolve_h5_output(case_dir: Path, before_run: dict[str, int] | None = None) -> Path: + candidates = [path for path in case_dir.glob("*.h5") if path.is_file()] + if not candidates: + raise FileNotFoundError(f"No HDF5 output found in {case_dir}") + + if before_run is not None: + fresh = [ + path + for path in candidates + if before_run.get(path.name) != path.stat().st_mtime_ns + ] + if fresh: + return max(fresh, key=lambda path: path.stat().st_mtime_ns) + + generated = [path for path in candidates if path.name != "answer.h5"] + if generated: + return max(generated, key=lambda path: path.stat().st_mtime_ns) + + return max(candidates, key=lambda path: path.stat().st_mtime_ns) + + +def run_case( + case: CaseInfo, + particle_override: int | None, + data_library_dir: Path, + *, + mode: str, + target: str, + caching: bool, + clear_cache: bool, + progress_bar: bool, + runtime_output: bool, + mpi_exec: str | None, + mpi_procs_request: str | int, +) -> Path: + runtime_input = build_runtime_input(case, particle_override) + env = build_runtime_env(data_library_dir) + before_run = snapshot_h5_files(case.case_dir) + metadata = load_case_metadata(case) + particle_count = particle_override or int(metadata["n_particles"]) + mpi_procs = resolve_mpi_procs(mpi_procs_request, particle_count, mpi_exec) + + with tempfile.NamedTemporaryFile( + mode="w", + suffix="_runtime_input.py", + prefix="tmp_", + dir=case.case_dir, + delete=False, + encoding="utf-8", + ) as handle: + handle.write(runtime_input) + temp_input_path = Path(handle.name) + + try: + command = build_mcdc_command( + python_executable=sys.executable, + input_name=temp_input_path.name, + mode=mode, + target=target, + caching=caching, + clear_cache=clear_cache, + progress_bar=progress_bar, + runtime_output=runtime_output, + mpi_exec=mpi_exec, + mpi_procs=mpi_procs, + ) + subprocess.run( + command, + cwd=case.case_dir, + env=env, + check=True, + ) + finally: + temp_input_path.unlink(missing_ok=True) + + return resolve_h5_output(case.case_dir, before_run=before_run) + + +def read_h5_dataset(file_handle, paths: list[str]): + import numpy as np + + for path in paths: + if path in file_handle: + return np.atleast_1d(file_handle[path][()]).astype(float) + joined = ", ".join(paths) + raise KeyError(f"Could not find any of the datasets: {joined}") + + +def read_h5_scalar_int(file_handle, paths: list[str]) -> int | None: + for path in paths: + if path in file_handle: + return int(file_handle[path][()]) + return None + + +def process_case( + case: CaseInfo, + h5_path: Path | None = None, + particle_override: int | None = None, +) -> Path: + try: + import h5py + import matplotlib + import numpy as np + except ModuleNotFoundError as exc: + raise ModuleNotFoundError( + "Figure generation requires `numpy`, `h5py`, and `matplotlib` " + "in the active Python environment." + ) from exc + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + if not case.reference_path.exists(): + raise FileNotFoundError(f"reference.npz not found in {case.case_dir}") + + if h5_path is None: + h5_path = resolve_h5_output(case.case_dir) + + metadata = load_case_metadata(case) + effective_particles = particle_override or metadata["n_particles"] + + with np.load(case.reference_path) as ref: + theo_fmr = ref["fmr_theo_tiger"].astype(float) + theo_edep = ref["edep_theo_tiger"].astype(float) + exp_fmr = ref["fmr_exp_lw"].astype(float) + exp_edep = ref["edep_exp_lw"].astype(float) + + with h5py.File(h5_path, "r") as file: + z = read_h5_dataset( + file, + ["tallies/edep/grid/z", "tallies/mesh_tally_0/grid/z"], + ) + edep = read_h5_dataset( + file, + [ + "tallies/edep/energy_deposition/mean", + "tallies/mesh_tally_0/edep/mean", + ], + ) + h5_particles = read_h5_scalar_int(file, ["input_deck/setting/N_particle"]) + if h5_particles is not None: + effective_particles = h5_particles + + dz = z[1:] - z[:-1] + z_centers = 0.5 * (z[:-1] + z[1:]) + total_thickness = float(metadata["csda_range"]) / float(metadata["rho_g_cm3"]) + mcdc_fmr = z_centers / total_thickness + edep_mcdc = edep / float(metadata["rho_g_cm3"]) / dz / 1.0e6 + + case.results_dir.mkdir(parents=True, exist_ok=True) + energy_mev = case.energy_keV / 1000.0 + out_name = ( + f"fig_{case.material}_{energy_mev:g}MeV_" + f"th{case.angle_deg}_{format_scientific_int(int(effective_particles))}.png" + ) + out_path = case.results_dir / out_name + + plt.rcParams.update( + { + "font.size": 16, + "axes.titlesize": 20, + "axes.labelsize": 18, + "legend.fontsize": 14, + "xtick.labelsize": 14, + "ytick.labelsize": 14, + } + ) + + fig, ax = plt.subplots(figsize=(12, 7), constrained_layout=True) + ax.minorticks_on() + ax.grid(True, which="major", linestyle="-", linewidth=0.8, alpha=0.4) + ax.grid(True, which="minor", linestyle=":", linewidth=0.6, alpha=0.5) + ax.plot( + mcdc_fmr, + edep_mcdc, + label="MCDC Simulation", + marker="o", + markersize=7, + linewidth=2.5, + ) + ax.plot( + theo_fmr, + theo_edep, + label="TIGER Theoretical", + marker="s", + markersize=7, + linewidth=2.5, + ) + ax.plot( + exp_fmr, + exp_edep, + label="Lockwood Experimental", + marker="D", + markersize=7, + linewidth=2.5, + ) + ax.set_xlim(0.0, 1.0) + ax.set_xlabel("Fraction of Mean Range") + ax.set_ylabel("Energy Deposition (MeV/g/cm²)") + ax.set_title( + f"Energy Deposition of {energy_mev:g} MeV Electrons in " + f"{case.material} at {case.angle_deg:g}° Incidence" + ) + ax.legend() + fig.savefig(out_path, dpi=400, bbox_inches="tight") + plt.close(fig) + + return out_path + + +def print_case_table(cases: list[CaseInfo]) -> None: + for case in cases: + print(f"- {case.label} -> {case.case_dir.relative_to(CASE_ROOT)}") + + +def build_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + description="Run MCDC VV benchmark cases and generate result figures.", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument( + "--mode", + choices=["python", "numba", "numba_debug"], + default=DEFAULT_RUN_MODE, + help="MCDC execution mode.", + ) + parser.add_argument( + "--target", + choices=["cpu", "gpu"], + default=DEFAULT_TARGET, + help="MCDC execution target.", + ) + parser.add_argument( + "--mpi-procs", + type=parse_mpi_procs, + default=DEFAULT_MPI_PROCS, + help="MPI process count. Use `auto` to use all available CPU cores up to the particle count.", + ) + parser.add_argument( + "--mpi-exec", + help="MPI launcher to use, e.g. mpiexec or mpirun.", + ) + parser.add_argument( + "--caching", + dest="caching", + action="store_true", + help="Enable MCDC/Numba caching for faster repeated runs.", + ) + parser.add_argument( + "--no-caching", + dest="caching", + action="store_false", + help="Disable MCDC/Numba caching.", + ) + parser.set_defaults(caching=True) + parser.add_argument( + "--clear-cache", + action="store_true", + help="Clear MCDC/Numba caches before the run.", + ) + parser.add_argument( + "--progress-bar", + dest="progress_bar", + action="store_true", + help="Show the MCDC progress bar during transport.", + ) + parser.add_argument( + "--no-progress-bar", + dest="progress_bar", + action="store_false", + help="Disable the MCDC progress bar for slightly cleaner and lighter runs.", + ) + parser.set_defaults(progress_bar=False) + parser.add_argument( + "--runtime-output", + action="store_true", + help="Ask MCDC to emit extra runtime datasets.", + ) + parser.add_argument( + "--material", + help="Filter cases by material name, e.g. Al.", + ) + parser.add_argument( + "--energy", + type=float, + help="Filter cases by incident energy in keV, e.g. 500 or 1000.", + ) + parser.add_argument( + "--angle", + type=float, + help="Filter cases by incident angle in degrees, e.g. 0 or 60.", + ) + parser.add_argument( + "--particles", + type=parse_positive_int, + help="Override the particle count for selected runs. Scientific notation is allowed.", + ) + parser.add_argument( + "--data-library", + help="Path to the processed MCDC electron data library directory.", + ) + parser.add_argument( + "--list", + action="store_true", + help="List matching cases and exit.", + ) + parser.add_argument( + "--run-only", + action="store_true", + help="Run selected cases but skip figure generation.", + ) + parser.add_argument( + "--process-only", + action="store_true", + help="Skip simulation and only generate figures from existing HDF5 outputs.", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Show resolved configuration without running anything.", + ) + return parser + + +def main(argv: list[str] | None = None) -> int: + parser = build_parser() + args = parser.parse_args(argv) + + if args.run_only and args.process_only: + parser.error("--run-only and --process-only cannot be used together.") + + cases = discover_cases() + selected = filter_cases( + cases, + material=args.material, + energy_keV=args.energy, + angle_deg=args.angle, + ) + + if not selected: + print("No matching cases found.\nAvailable cases:") + print_case_table(cases) + return 1 + + data_library_dir = None + mpi_exec = resolve_mpi_exec(args.mpi_exec) + if args.dry_run: + data_library_dir = None if args.process_only else resolve_data_library_dir(args.data_library) + + if args.list or args.dry_run: + resolved_mpi_procs = None + if not args.process_only: + first_case = selected[0] + first_metadata = load_case_metadata(first_case) + first_particles = args.particles or int(first_metadata["n_particles"]) + resolved_mpi_procs = resolve_mpi_procs( + args.mpi_procs, + first_particles, + mpi_exec, + ) + + print("Selected cases:") + print_case_table(selected) + print(f"\nMCDC root: {MCDC_ROOT}") + print(f"Run mode: {args.mode}") + print(f"Target: {args.target}") + print(f"MPI launcher: {mpi_exec or ''}") + print(f"MPI processes: {args.mpi_procs}") + if resolved_mpi_procs is not None: + print(f"Resolved MPI processes: {resolved_mpi_procs}") + print(f"Caching: {args.caching}") + print(f"Progress bar: {args.progress_bar}") + if data_library_dir is not None: + print(f"Data library: {data_library_dir}") + if args.particles is not None: + print(f"Particle override: {args.particles}") + if args.dry_run: + return 0 + if args.list and not args.dry_run: + return 0 + + figure_paths: list[Path] = [] + h5_paths: list[Path] = [] + + if not args.process_only and data_library_dir is None: + data_library_dir = resolve_data_library_dir(args.data_library) + + for index, case in enumerate(selected, start=1): + print(f"\n[{index}/{len(selected)}] {case.label}") + + h5_path = None + if not args.process_only: + metadata = load_case_metadata(case) + particle_count = args.particles or int(metadata["n_particles"]) + mpi_procs = resolve_mpi_procs(args.mpi_procs, particle_count, mpi_exec) + print(f"Using MCDC root: {MCDC_ROOT}") + print(f"Using data library: {data_library_dir}") + print(f"Using mode: {args.mode}") + print(f"Using target: {args.target}") + print(f"Using MPI launcher: {mpi_exec or ''}") + print(f"Using MPI processes: {mpi_procs}") + print(f"Using caching: {args.caching}") + print(f"Running case in {case.case_dir}") + h5_path = run_case( + case, + args.particles, + data_library_dir, + mode=args.mode, + target=args.target, + caching=args.caching, + clear_cache=args.clear_cache, + progress_bar=args.progress_bar, + runtime_output=args.runtime_output, + mpi_exec=mpi_exec, + mpi_procs_request=args.mpi_procs, + ) + h5_paths.append(h5_path) + print(f"Produced HDF5: {h5_path}") + + if not args.run_only: + process_input = h5_path if h5_path is not None else resolve_h5_output(case.case_dir) + print(f"Generating figure from {process_input}") + figure_path = process_case(case, h5_path=process_input, particle_override=args.particles) + figure_paths.append(figure_path) + print(f"Saved figure: {figure_path}") + + print("\nCompleted.") + if h5_paths: + print(f"HDF5 files: {len(h5_paths)}") + if figure_paths: + print(f"Figures: {len(figure_paths)}") + return 0 + + +if __name__ == "__main__": + try: + raise SystemExit(main()) + except ( + FileNotFoundError, + KeyError, + ModuleNotFoundError, + subprocess.CalledProcessError, + ValueError, + ) as exc: + print(f"Error: {exc}", file=sys.stderr) + raise SystemExit(1) diff --git a/mcdc/object_/distribution.py b/mcdc/object_/distribution.py index 5585c29c2..652f29ded 100644 --- a/mcdc/object_/distribution.py +++ b/mcdc/object_/distribution.py @@ -155,16 +155,24 @@ class DistributionMultiTable(DistributionBase): pdf: NDArray[float64] cdf: NDArray[float64] - def __init__(self, grid, offset, value, pdf): + def __init__(self, grid, offset, value, pdf=None, cdf=None): type_ = DISTRIBUTION_MULTITABLE super().__init__(type_) self.grid = grid self.offset = offset self.value = value - self.pdf = pdf - self.pdf, self.cdf = multi_cdf_from_pdf(offset, value, pdf) + if pdf is not None and cdf is not None: + self.pdf = pdf + self.cdf = cdf + elif cdf is not None: + # Direct-CDF mode for electron data + self.pdf = np.array([], dtype=float64) + self.cdf = cdf + else: + self.pdf = pdf + self.pdf, self.cdf = multi_cdf_from_pdf(offset, value, pdf) def __repr__(self): text = super().__repr__() diff --git a/mcdc/object_/electron_reaction.py b/mcdc/object_/electron_reaction.py index 0feeb8fff..68b59b05b 100644 --- a/mcdc/object_/electron_reaction.py +++ b/mcdc/object_/electron_reaction.py @@ -68,6 +68,35 @@ def decode_reference_frame(type_): return "Center of mass" +def load_multi_table_distribution(h5_group): + if "PDF" in h5_group and "CDF" in h5_group: + return DistributionMultiTable( + h5_group["energy_grid"][()], + h5_group["energy_offset"][()], + h5_group["value"][()], + pdf=h5_group["PDF"][()], + cdf=h5_group["CDF"][()], + ) + + if "PDF" in h5_group: + return DistributionMultiTable( + h5_group["energy_grid"][()], + h5_group["energy_offset"][()], + h5_group["value"][()], + h5_group["PDF"][()], + ) + + if "CDF" in h5_group: + return DistributionMultiTable( + h5_group["energy_grid"][()], + h5_group["energy_offset"][()], + h5_group["value"][()], + cdf=h5_group["CDF"][()], + ) + + raise KeyError("Expected either PDF or CDF dataset in multi-table distribution") + + # ====================================================================================== # Electron ionization # ====================================================================================== @@ -117,14 +146,7 @@ def from_h5_group(cls, h5_group): # Secondary electron energy distribution product = subshell["product"] - subshell_product.append( - DistributionMultiTable( - product["energy_grid"][()], - product["energy_offset"][()], - product["value"][()], - product["PDF"][()], - ) - ) + subshell_product.append(load_multi_table_distribution(product)) return cls( MT, @@ -180,14 +202,7 @@ def from_h5_group(cls, h5_group): large_angle = h5_group["large_angle"] xs_large = DataTable(large_angle["xs_energy"][()], large_angle["xs"][()]) - - mu_group = large_angle["scattering_cosine"] - mu = DistributionMultiTable( - mu_group["energy_grid"][()], - mu_group["energy_offset"][()], - mu_group["value"][()], - mu_group["PDF"][()], - ) + mu = load_multi_table_distribution(large_angle["scattering_cosine"]) return cls(MT, xs, xs_offset, reference_frame, xs_large, mu) diff --git a/mcdc/object_/element.py b/mcdc/object_/element.py index acc1270a0..9af4963ee 100644 --- a/mcdc/object_/element.py +++ b/mcdc/object_/element.py @@ -83,29 +83,41 @@ def set_electron_data(self): # ========================================================================== self.electron_xs_energy_grid = file["electron_reactions/xs_energy_grid"][()] - self.electron_total_xs = np.zeros_like(self.electron_xs_energy_grid) - self.electron_elastic_xs = np.zeros_like(self.electron_xs_energy_grid) - self.electron_excitation_xs = np.zeros_like(self.electron_xs_energy_grid) - self.electron_bremsstrahlung_xs = np.zeros_like(self.electron_xs_energy_grid) - self.electron_ionization_xs = np.zeros_like(self.electron_xs_energy_grid) - - xs_containers = [ - self.electron_elastic_xs, - self.electron_excitation_xs, - self.electron_bremsstrahlung_xs, - self.electron_ionization_xs, - ] - for xs_container, rx_name in list(zip(xs_containers, rx_names)): + def load_xs_dataset(dataset): + values = np.zeros_like(self.electron_xs_energy_grid) + offset = int(dataset.attrs.get("offset", 0)) + data = dataset[()] + values[offset : offset + len(data)] = data + return values + + def load_reaction_total_xs(rx_name): + total_path = f"electron_reactions/{rx_name}/total/xs" + if total_path in file: + return load_xs_dataset(file[total_path]) + + values = np.zeros_like(self.electron_xs_energy_grid) for MT in MTs[rx_name]: xs = file[f"electron_reactions/{rx_name}/{MT}/xs"] - xs_container[xs.attrs["offset"] :] += xs[()] - - self.electron_total_xs = ( - self.electron_elastic_xs - + self.electron_excitation_xs - + self.electron_bremsstrahlung_xs - + self.electron_ionization_xs - ) + offset = int(xs.attrs.get("offset", 0)) + data = xs[()] + values[offset : offset + len(data)] += data + return values + + self.electron_elastic_xs = load_reaction_total_xs("elastic_scattering") + self.electron_excitation_xs = load_reaction_total_xs("excitation") + self.electron_bremsstrahlung_xs = load_reaction_total_xs("bremsstrahlung") + self.electron_ionization_xs = load_reaction_total_xs("ionization") + + total_path = "electron_reactions/total/xs" + if total_path in file: + self.electron_total_xs = load_xs_dataset(file[total_path]) + else: + self.electron_total_xs = ( + self.electron_elastic_xs + + self.electron_excitation_xs + + self.electron_bremsstrahlung_xs + + self.electron_ionization_xs + ) # ========================================================================== # The reactions diff --git a/mcdc/transport/distribution.py b/mcdc/transport/distribution.py index a2bd93e87..1a3e09d97 100644 --- a/mcdc/transport/distribution.py +++ b/mcdc/transport/distribution.py @@ -307,27 +307,47 @@ def _sample_multi_table(E, rng_state, multi_table, data, scale): # The CDF offset = multi_table["cdf_offset"] cdf = data[start + offset : start + offset + size] - # Above is equivalent to: cdf = mcdc_get.multi_table_distribution.cdf_chunk(start, size, multi_table, data) - # Generate random numbers + # Generate random number xi = rng.lcg(rng_state) - # Sample bin index - idx = find_bin(xi, cdf) - c = cdf[idx] + # Degenerate table + if size == 1: + sample = mcdc_get.multi_table_distribution.value(start, multi_table, data) - # Get the other values - idx += start # Apply the offset as these are not chunk-extracted like the cdf - p0 = mcdc_get.multi_table_distribution.pdf(idx, multi_table, data) - p1 = mcdc_get.multi_table_distribution.pdf(idx + 1, multi_table, data) - val0 = mcdc_get.multi_table_distribution.value(idx, multi_table, data) - val1 = mcdc_get.multi_table_distribution.value(idx + 1, multi_table, data) - - m = (p1 - p0) / (val1 - val0) - if m == 0.0: - sample = val0 + (xi - c) / p0 else: - sample = val0 + 1.0 / m * (math.sqrt(p0**2 + 2 * m * (xi - c)) - p0) + # Sample bin index from the tabulated CDF + idx = find_bin(xi, cdf) + + # Direct-CDF mode: used by electron data written from ACE files + if multi_table["pdf_length"] == 0: + + c0 = cdf[idx] + c1 = cdf[idx + 1] + + idx += start + val0 = mcdc_get.multi_table_distribution.value(idx, multi_table, data) + val1 = mcdc_get.multi_table_distribution.value(idx + 1, multi_table, data) + + frac = (xi - c0) / (c1 - c0) + if frac < 0.0: + frac = 0.0 + elif frac > 1.0: + frac = 1.0 + sample = val0 + frac * (val1 - val0) + + # PDF-based mode: used by neutron data and PyEPICS libraries + else: + c = cdf[idx] + + idx += start # Apply the offset + p0 = mcdc_get.multi_table_distribution.pdf(idx, multi_table, data) + p1 = mcdc_get.multi_table_distribution.pdf(idx + 1, multi_table, data) + val0 = mcdc_get.multi_table_distribution.value(idx, multi_table, data) + val1 = mcdc_get.multi_table_distribution.value(idx + 1, multi_table, data) + + m = (p1 - p0) / (val1 - val0) + sample = val0 + 1.0 / m * (math.sqrt(p0**2 + 2 * m * (xi - c)) - p0) if not scale: return sample diff --git a/mcdc/transport/physics/electron/native.py b/mcdc/transport/physics/electron/native.py index 4c296ceb5..068b91c0d 100644 --- a/mcdc/transport/physics/electron/native.py +++ b/mcdc/transport/physics/electron/native.py @@ -26,13 +26,12 @@ from mcdc.transport.data import evaluate_data from mcdc.transport.distribution import ( sample_distribution, - sample_multi_table, ) from mcdc.transport.physics.util import ( evaluate_electron_xs_energy_grid, scatter_direction, ) -from mcdc.transport.util import linear_interpolation +from mcdc.transport.util import find_bin, linear_interpolation # ====================================================================================== # Particle attributes @@ -268,55 +267,284 @@ def elastic_scattering(reaction, particle_container, element, simulation, data): # Current energy E = particle["E"] + Z = int(element["atomic_number"]) + multi_table = simulation["multi_table_distributions"][reaction["mu_ID"]] + mu0 = sample_coupled_elastic_mu(E, Z, particle_container, multi_table, data) - # ------------------------------------------------------------------------- - # Total elastic xs - # ------------------------------------------------------------------------- + # Update direction + azi = 2.0 * PI * rng.lcg(particle_container) + ux = particle["ux"] + uy = particle["uy"] + uz = particle["uz"] + ux_new, uy_new, uz_new = scatter_direction(ux, uy, uz, mu0, azi) + + particle["ux"] = ux_new + particle["uy"] = uy_new + particle["uz"] = uz_new - reaction_base_ID = int(reaction["parent_ID"]) - reaction_base = simulation["electron_reactions"][reaction_base_ID] - xs_total = reaction_micro_xs(E, reaction_base, element, data) - # If large-angle, xs from data table - xs_large = elastic_large_xs(E, reaction, simulation, data) +@njit +def sample_coupled_elastic_mu(E, Z, rng_state, multi_table, data): + idx0, idx1, frac = elastic_table_energy_interval(E, multi_table, data) + fu0, mu_n0, eta0 = elastic_forward_branch_probability(idx0, Z, multi_table, data) - # Important to check because of numerical issues - if xs_large < 0.0: - xs_large = 0.0 - if xs_large > xs_total: - xs_large = xs_total + if idx1 == idx0: + fu = fu0 + else: + fu1, mu_n1, eta1 = elastic_forward_branch_probability( + idx1, Z, multi_table, data + ) + fu = fu0 + frac * (fu1 - fu0) + + if fu < 0.0: + fu = 0.0 + elif fu > 1.0: + fu = 1.0 + + if rng.lcg(rng_state) < fu: + u = rng.lcg(rng_state) + if idx1 == idx0: + mu = sample_forward_peak_mu(eta0, mu_n0, u) + else: + mu0 = sample_forward_peak_mu(eta0, mu_n0, u) + mu1 = sample_forward_peak_mu(eta1, mu_n1, u) + mu = mu0 + frac * (mu1 - mu0) + else: + u = rng.lcg(rng_state) + if idx1 == idx0: + mu = sample_multi_table_fixed_index(idx0, u, multi_table, data) + else: + mu = sample_multi_table_log_cdf(E, idx0, idx1, u, multi_table, data) - prob_large = xs_large / xs_total - mu_cut = float(reaction["mu_cut"]) + if mu < -1.0: + return -1.0 + if mu > 1.0: + return 1.0 + return mu - xi = rng.lcg(particle_container) - if xi < prob_large: - # --------------------------------------------------------------------- - # Large-angle elastic scattering - # --------------------------------------------------------------------- +@njit +def elastic_table_energy_interval(E, multi_table, data): + grid = mcdc_get.multi_table_distribution.grid_all(multi_table, data) + if E <= grid[0]: + return 0, 0, 0.0 + if E >= grid[-1]: + idx = len(grid) - 1 + return idx, idx, 0.0 + + idx0 = find_bin(E, grid) + idx1 = idx0 + 1 + E0 = grid[idx0] + E1 = grid[idx1] + frac = (E - E0) / (E1 - E0) + return idx0, idx1, frac - multi_table = simulation["multi_table_distributions"][reaction["mu_ID"]] - mu0 = sample_multi_table(E, particle_container, multi_table, data) +@njit +def elastic_table_bounds(index, multi_table, data): + start = int(mcdc_get.multi_table_distribution.offset(index, multi_table, data)) + if index + 1 == multi_table["grid_length"]: + end = multi_table["value_length"] else: - # --------------------------------------------------------------------- - # Small-angle elastic scattering (Coulomb tail sampling) - # --------------------------------------------------------------------- + end = int(mcdc_get.multi_table_distribution.offset(index + 1, multi_table, data)) + return start, end - Z = int(element["atomic_number"]) - mu0 = sample_small_angle_mu_coulomb(E, Z, particle_container, mu_cut) - # Update direction - azi = 2.0 * PI * rng.lcg(particle_container) - ux = particle["ux"] - uy = particle["uy"] - uz = particle["uz"] - ux_new, uy_new, uz_new = scatter_direction(ux, uy, uz, mu0, azi) +@njit +def elastic_table_area(index, multi_table, data): + start, end = elastic_table_bounds(index, multi_table, data) + size = end - start + if size < 2: + return 0.0 - particle["ux"] = ux_new - particle["uy"] = uy_new - particle["uz"] = uz_new + if multi_table["pdf_length"] == 0: + c0 = mcdc_get.multi_table_distribution.cdf(start, multi_table, data) + c1 = mcdc_get.multi_table_distribution.cdf(end - 1, multi_table, data) + return c1 - c0 + + total = 0.0 + for i in range(start, end - 1): + mu0 = mcdc_get.multi_table_distribution.value(i, multi_table, data) + mu1 = mcdc_get.multi_table_distribution.value(i + 1, multi_table, data) + p0 = mcdc_get.multi_table_distribution.pdf(i, multi_table, data) + p1 = mcdc_get.multi_table_distribution.pdf(i + 1, multi_table, data) + total += 0.5 * (p0 + p1) * (mu1 - mu0) + return total + + +@njit +def elastic_table_tail_density(index, multi_table, data): + start, end = elastic_table_bounds(index, multi_table, data) + size = end - start + if size < 2: + return 0.0 + + mu_prev = mcdc_get.multi_table_distribution.value(end - 2, multi_table, data) + mu_n = mcdc_get.multi_table_distribution.value(end - 1, multi_table, data) + dmu = mu_n - mu_prev + if dmu <= 0.0: + return 0.0 + + if multi_table["pdf_length"] == 0: + c_prev = mcdc_get.multi_table_distribution.cdf(end - 2, multi_table, data) + c_last = mcdc_get.multi_table_distribution.cdf(end - 1, multi_table, data) + return (c_last - c_prev) / dmu + + return mcdc_get.multi_table_distribution.pdf(end - 1, multi_table, data) + + +@njit +def elastic_forward_branch_probability(index, Z, multi_table, data): + energy = mcdc_get.multi_table_distribution.grid(index, multi_table, data) + eta = compute_scattering_eta(energy, Z) + + start, end = elastic_table_bounds(index, multi_table, data) + size = end - start + if size < 2: + mu_n = mcdc_get.multi_table_distribution.value(end - 1, multi_table, data) + return 0.0, mu_n, eta + + mu_n = mcdc_get.multi_table_distribution.value(end - 1, multi_table, data) + p_n = elastic_table_tail_density(index, multi_table, data) + t_t = elastic_table_area(index, multi_table, data) + + x_n = eta + 1.0 - mu_n + if p_n <= 0.0 or x_n <= 0.0 or eta <= 0.0: + return 0.0, mu_n, eta + + A = p_n * x_n * x_n + t_fu = A * ((1.0 / eta) - (1.0 / x_n)) + denom = t_fu + t_t + if denom <= 0.0: + return 0.0, mu_n, eta + + return t_fu / denom, mu_n, eta + + +@njit +def sample_forward_peak_mu(eta, mu_n, u): + x_n = eta + 1.0 - mu_n + if x_n <= 0.0 or eta <= 0.0: + return mu_n + + inv_n = 1.0 / x_n + inv = inv_n + u * ((1.0 / eta) - inv_n) + if inv <= 0.0: + return mu_n + + return 1.0 + eta - (1.0 / inv) + + +@njit +def sample_multi_table_fixed_index(index, xi, multi_table, data): + start, end = elastic_table_bounds(index, multi_table, data) + size = end - start + if size <= 1: + return mcdc_get.multi_table_distribution.value(start, multi_table, data) + + cdf_offset = multi_table["cdf_offset"] + cdf = data[start + cdf_offset : start + cdf_offset + size] + idx = find_bin(xi, cdf) + c0 = cdf[idx] + c1 = cdf[idx + 1] + if c1 <= c0: + return mcdc_get.multi_table_distribution.value(start + idx, multi_table, data) + + val0 = mcdc_get.multi_table_distribution.value(start + idx, multi_table, data) + val1 = mcdc_get.multi_table_distribution.value(start + idx + 1, multi_table, data) + frac = (xi - c0) / (c1 - c0) + if frac < 0.0: + frac = 0.0 + elif frac > 1.0: + frac = 1.0 + return val0 + frac * (val1 - val0) + + +@njit +def evaluate_multi_table_cdf_at_mu(index, mu, multi_table, data): + start, end = elastic_table_bounds(index, multi_table, data) + size = end - start + if size <= 0: + return 0.0 + + value_offset = multi_table["value_offset"] + values = data[start + value_offset : start + value_offset + size] + if mu <= values[0]: + return 0.0 + if mu >= values[-1]: + return 1.0 + + cdf_offset = multi_table["cdf_offset"] + cdf = data[start + cdf_offset : start + cdf_offset + size] + idx = find_bin(mu, values) + mu0 = values[idx] + mu1 = values[idx + 1] + c0 = cdf[idx] + c1 = cdf[idx + 1] + + if mu1 <= mu0: + return c0 + return c0 + (mu - mu0) * (c1 - c0) / (mu1 - mu0) + + +@njit +def elastic_log_energy_fraction(E, E0, E1): + if E <= 0.0 or E0 <= 0.0 or E1 <= 0.0: + return (E - E0) / (E1 - E0) + return (math.log(E) - math.log(E0)) / (math.log(E1) - math.log(E0)) + + +@njit +def log_interpolate_cdf_value(E, E0, E1, c0, c1): + eps = 1.0e-300 + if c0 <= 0.0 and c1 <= 0.0: + return 0.0 + if c0 >= 1.0 and c1 >= 1.0: + return 1.0 + + y0 = max(min(c0, 1.0), eps) + y1 = max(min(c1, 1.0), eps) + frac = elastic_log_energy_fraction(E, E0, E1) + if frac < 0.0: + frac = 0.0 + elif frac > 1.0: + frac = 1.0 + return math.exp(math.log(y0) + frac * (math.log(y1) - math.log(y0))) + + +@njit +def sample_multi_table_log_cdf(E, idx0, idx1, xi, multi_table, data): + if idx0 == idx1: + return sample_multi_table_fixed_index(idx0, xi, multi_table, data) + + E0 = mcdc_get.multi_table_distribution.grid(idx0, multi_table, data) + E1 = mcdc_get.multi_table_distribution.grid(idx1, multi_table, data) + + start0, end0 = elastic_table_bounds(idx0, multi_table, data) + start1, end1 = elastic_table_bounds(idx1, multi_table, data) + lo0 = mcdc_get.multi_table_distribution.value(start0, multi_table, data) + lo1 = mcdc_get.multi_table_distribution.value(start1, multi_table, data) + hi0 = mcdc_get.multi_table_distribution.value(end0 - 1, multi_table, data) + hi1 = mcdc_get.multi_table_distribution.value(end1 - 1, multi_table, data) + + lo = min(lo0, lo1) + hi = max(hi0, hi1) + if xi <= 0.0: + return lo + if xi >= 1.0: + return hi + + for _ in range(64): + mid = 0.5 * (lo + hi) + c0 = evaluate_multi_table_cdf_at_mu(idx0, mid, multi_table, data) + c1 = evaluate_multi_table_cdf_at_mu(idx1, mid, multi_table, data) + c = log_interpolate_cdf_value(E, E0, E1, c0, c1) + if c < xi: + lo = mid + else: + hi = mid + + return 0.5 * (lo + hi) @njit @@ -348,12 +576,6 @@ def sample_small_angle_mu_coulomb(E, Z, rng_state, mu_cut): return 1.0 - x -@njit -def elastic_large_xs(E, reaction, simulation, data): - data_base = simulation["data"][int(reaction["xs_large_ID"])] - return evaluate_data(E, data_base, simulation, data) - - # ====================================================================================== # Excitation # ======================================================================================