Feature/northroj/icsbep

icsbep/README.md

@@ -0,0 +1,18 @@
+# Instructions for MC/DC's ICSBEP validation suite
+
+This continuous energy criticality validation test suite is taken from OpenMC's benchmark repository: https://github.com/mit-crpg/benchmarks/tree/master/icsbep
+
+It contains several hundred criticality problems from the international criticality safety benchmark evaluation project (ICSBEP). Currently only 118 of them translate and run with the current state of the translator.
+
+
+## Running problems:
+
+Inside the icsbep directory, the icsbep_translator.ipynb file is a slightly modified version of the OpenMC to MC/DC xml translator found in the CEMeNT organization's repositories (may be private).
+
+There is also an openmc directory which houses all the folders for the input decks and the bash scripts to run them. Inside each problem directory, an openmc directory has subdirectories for each case, and an mcdc directory places all the translated input files in a single place.
+
+Inside the base icsbeb/openmc directory, there are two bash scripts: run_mcdc.sh and run_openmc.sh. These should scan all subdirectories, create slurm submission scripts for each input deck, and submit those decks to the pbatch partition. Please note that this has been set up for LLNL's Dane machine, and each of these input scripts will need to be customized to the user.
+
+The python script collect_results.py will scan those subdirectories and scrape the keff results for both codes.
+
+Also note that OpenMC's default submission scripts have been left in, but they will not work for MC/DC or on LLNL's systems.

icsbep/openmc/collect_results.py

@@ -0,0 +1,146 @@
+import os
+import h5py
+import pandas as pd
+import re
+import numpy as np
+import matplotlib.pyplot as plt
+
+def extract_mcdc_results(mcdc_dir):
+    results = {}
+    for file in os.listdir(mcdc_dir):
+        if file.endswith(".h5"):
+            file_path = os.path.join(mcdc_dir, file)
+            print("Extracting MC/DC results from:", file_path)
+            try:
+                with h5py.File(file_path, 'r') as f:
+                    k_mean = np.array(f["k_mean"])
+                    k_sdev = np.array(f["k_sdev"])
+                    print(k_mean)
+                    if k_mean is not None and k_sdev is not None:
+                        problem_name = file.replace(".h5", "")
+                        results[problem_name] = {"mcdc_k_mean": k_mean, "mcdc_k_sdev": k_sdev}
+            except Exception as e:
+                print(f"Error reading {file_path}: {e}")
+    return results
+
+def extract_openmc_results(openmc_dir):
+    results = {}
+    parent_dir = os.path.basename(os.path.dirname(openmc_dir))
+    case_dirs = [os.path.join(openmc_dir, d) for d in os.listdir(openmc_dir) if d.startswith("case-") and os.path.isdir(os.path.join(openmc_dir, d))]
+
+    if not case_dirs:
+        case_dirs = [openmc_dir]  # If no case-# subdirectories, use openmc_dir itself
+
+    for case_dir in case_dirs:
+        stdout_file = os.path.join(case_dir, "openmc.stdout")
+        if os.path.exists(stdout_file):
+            print("Extracting openmc results from:", case_dir)
+            try:
+                with open(stdout_file, 'r') as f:
+                    for line in f:
+                        match = re.search(r'Combined k-effective\s+=\s+([\d\.]+)\s+\+/-\s+([\d\.]+)', line)
+                        if match:
+                            k_mean = float(match.group(1))
+                            k_sdev = float(match.group(2))
+                            case_name = os.path.basename(case_dir) if "case-" in case_dir else ""
+                            problem_name = f"{parent_dir}_{case_name}" if case_name else parent_dir
+                            results[problem_name] = {"openmc_k_mean": k_mean, "openmc_k_sdev": k_sdev}
+                            break  # Stop after finding the first match
+            except Exception as e:
+                print(f"Error reading {stdout_file}: {e}")
+    return results
+
+def collect_results1(root_dir):
+    dataset = {}
+
+    for dirpath, dirnames, _ in os.walk(root_dir):
+        if "mcdc" in dirnames:
+            mcdc_results = extract_mcdc_results(os.path.join(dirpath, "mcdc"))
+            dataset.update(mcdc_results)
+        if "openmc" in dirnames:
+            openmc_results = extract_openmc_results(os.path.join(dirpath, "openmc"))
+            dataset.update(openmc_results)
+    df = pd.DataFrame.from_dict(dataset, orient='index').reset_index().rename(columns={'index': 'problem_name'})
+    df.to_csv("results_summary.csv", index=False)
+    print("Results collected and saved to results_summary.csv")
+
+def collect_results(root_dir):
+    dataset = {}
+
+    for dirpath, dirnames, _ in os.walk(root_dir):
+        if "mcdc" in dirnames:
+            mcdc_results = extract_mcdc_results(os.path.join(dirpath, "mcdc"))
+            for problem_name, values in mcdc_results.items():
+                if problem_name not in dataset:
+                    dataset[problem_name] = {}
+                dataset[problem_name].update(values)
+
+        if "openmc" in dirnames:
+            openmc_results = extract_openmc_results(os.path.join(dirpath, "openmc"))
+            for problem_name, values in openmc_results.items():
+                if problem_name not in dataset:
+                    dataset[problem_name] = {}
+                dataset[problem_name].update(values)
+
+    # Convert to DataFrame and ensure all columns exist
+    df = pd.DataFrame.from_dict(dataset, orient='index').reset_index().rename(columns={'index': 'problem_name'})
+
+    # Ensure all columns exist, filling missing values with NaN
+    expected_columns = ["problem_name", "openmc_k_mean", "openmc_k_sdev", "mcdc_k_mean", "mcdc_k_sdev"]
+    df = df.reindex(columns=expected_columns)
+
+    # Save to CSV
+    df.to_csv("results_summary.csv", index=False)
+    print("Results collected and saved to results_summary.csv")
+
+#root_directory = os.getcwd()  # Start search from current directory
+#collect_results(root_directory)
+
+
+
+
+### plotting
+plotting_toggle = 1
+
+def plot_icsbep_comparison(df, title):
+
+    # Set up the plot
+    plt.figure(figsize=(12, 6))
+    x_labels = df["problem_name"]
+    x = np.arange(len(x_labels))  # Create x positions for labels
+
+    # Plot OpenMC and MCDC k_mean values
+    plt.scatter(x, df["openmc_k_mean"], marker='o',  label="OpenMC" , color='b')
+    plt.errorbar(x, df["openmc_k_mean"], yerr=df["openmc_k_sdev"], fmt='o', color='b', capsize=3)
+    plt.scatter(x, df["mcdc_k_mean"], marker='s', label="MC/DC", color='r')
+    plt.errorbar(x, df["mcdc_k_mean"], yerr=df["mcdc_k_sdev"], fmt='o', color='r', capsize=3)
+
+    # Format the x-axis for better readability
+    plt.xticks(x, x_labels, rotation=45, ha="right", fontsize=8)  # Rotate labels
+    #plt.xlabel("Problem Name")
+    plt.ylabel("keff")
+    plt.title(title)
+    plt.legend()
+    plt.grid(axis="y", linestyle="--", alpha=0.7)
+
+# Load the CSV file
+df = pd.read_csv("results_summary.csv")
+df = df[:118]
+
+df_heu_met_fast = df[:25]
+df_heu_sol_therm = df[25:44]
+df_misc = pd.concat([df[44:48], df[56:57], df[117:118]])
+df_leu_sol_therm = df[48:56]
+df_pu_met_fast = df[57:64]
+df_pu_sol_therm = df[64:108]
+df_u233_met_fast = df[108:117]
+
+if plotting_toggle == 1:
+    plot_icsbep_comparison(df_heu_met_fast, "ICSBEP: HEU-MET-FAST")
+    plot_icsbep_comparison(df_heu_sol_therm, "ICSBEP: HEU-SOL-THERM")
+    plot_icsbep_comparison(df_misc, "ICSBEP: Miscellaneous")
+    plot_icsbep_comparison(df_leu_sol_therm, "ICSBEP: LEU-SOL-THERM")
+    plot_icsbep_comparison(df_pu_met_fast, "ICSBEP: PU-MET-FAST")
+    plot_icsbep_comparison(df_pu_sol_therm, "ICSBEP: PU-SOL-THERM")
+    plot_icsbep_comparison(df_u233_met_fast, "ICSBEP: U233-MET-FAST")
+

icsbep/openmc/heu-comp-inter-003/mcdc/heu-comp-inter-003_case-1.py

@@ -0,0 +1,218 @@
+import mcdc
+import numpy as np
+
+
+##############################
+#__________Materials__________
+##############################
+
+# Material Name: Aluminum
+m1 = mcdc.material(nuclides=[
+	['Al27',0.060239]])
+# Material Name: Depleted Uranium
+# Depletable
+m2 = mcdc.material(nuclides=[
+	['U235',9.736e-05],
+	['U238',0.047969]])
+# Material Name: SAE 1020 Steel
+m3 = mcdc.material(nuclides=[
+	['Fe54',0.00404865615],
+	['Fe56',0.06355524318],
+	['Fe57',0.00146776773],
+	['Fe58',0.00019533293999999998],
+	['C0',0.00054977],
+	['Mn55',0.00016969]])
+# Material Name: Fe
+m4 = mcdc.material(nuclides=[
+	['Fe54',0.0050169973],
+	['Fe56',0.07875612836],
+	['Fe57',0.00181882246],
+	['Fe58',0.00024205187999999997]])
+# Material Name: HEU, Can I
+# Depletable
+m11 = mcdc.material(nuclides=[
+	['U235',0.02345],
+	['U238',0.0013421],
+	['H1',0.07385949533046],
+	['H2',1.1504669540000001e-05],
+	['O16',0.00085520575413],
+	['O17',3.2424587e-07],
+	['C0',0.00014562],
+	['N14',0.0013543208841],
+	['N15',4.9791159e-06],
+	['Fe54',2.14377065e-06],
+	['Fe56',3.3652614580000005e-05],
+	['Fe57',7.7718563e-07],
+	['Fe58',1.0342914000000001e-07],
+	['Au197',0.0003398]])
+# Material Name: HEU, Can II
+# Depletable
+m12 = mcdc.material(nuclides=[
+	['U235',0.023409],
+	['U238',0.0013371],
+	['H1',0.07538425782696001],
+	['H2',1.174217304e-05],
+	['O16',0.00056681509563],
+	['O17',2.1490437000000002e-07],
+	['C0',0.00019893],
+	['N14',0.00019626842563],
+	['N15',7.2157437e-07],
+	['Fe54',2.2366477e-06],
+	['Fe56',3.5110585640000004e-05],
+	['Fe57',8.108565399999999e-07],
+	['Fe58',1.0791012000000002e-07],
+	['Au197',0.00033787]])
+# Material Name: HEU, Can III
+# Depletable
+m13 = mcdc.material(nuclides=[
+	['U235',0.02355],
+	['U238',0.0013584],
+	['H1',0.06990211174938],
+	['H2',1.088825062e-05],
+	['O16',0.00048405647303999993],
+	['O17',1.8352696e-07],
+	['C0',0.00015288],
+	['N14',0.00020568381028],
+	['N15',7.5618972e-07],
+	['Fe54',6.058926999999999e-06],
+	['Fe56',9.511219639999999e-05],
+	['Fe57',2.1965554e-06],
+	['Fe58',2.9232119999999997e-07],
+	['Au197',0.00034242]])
+# Material Name: HEU, Can IV
+# Depletable
+m14 = mcdc.material(nuclides=[
+	['U235',0.022939],
+	['U238',0.0013117],
+	['H1',0.0727136738085],
+	['H2',1.13261915e-05],
+	['O16',0.00063688852773],
+	['O17',2.4147227e-07],
+	['C0',0.00015539],
+	['N14',0.00024310622800000002],
+	['N15',8.93772e-07],
+	['Fe54',1.19442575e-06],
+	['Fe56',1.8749929900000002e-05],
+	['Fe57',4.3301765e-07],
+	['Fe58',5.76267e-08],
+	['Au197',0.00034196]])
+# Material Name: HEU, Can A
+# Depletable
+m15 = mcdc.material(nuclides=[
+	['U235',0.02303],
+	['U238',0.0013124],
+	['H1',0.07027905303540001],
+	['H2',1.09469646e-05],
+	['O16',0.00060008248251],
+	['O17',2.2751749e-07],
+	['C0',9.6792e-05],
+	['N14',0.00026208644785000006],
+	['N15',9.6355215e-07],
+	['Fe54',5.25202475e-06],
+	['Fe56',8.24455567e-05],
+	['Fe57',1.9040274499999997e-06],
+	['Fe58',2.533911e-07],
+	['Au197',0.00032844]])
+# Material Name: HEU, Can B
+# Depletable
+m16 = mcdc.material(nuclides=[
+	['U235',0.022951],
+	['U238',0.0013058],
+	['H1',0.0732285936024],
+	['H2',1.14063976e-05],
+	['O16',0.0004934129256],
+	['O17',1.870744e-07],
+	['C0',0.00012804],
+	['N14',0.00022934681403],
+	['N15',8.4318597e-07],
+	['Fe54',3.9574157e-06],
+	['Fe56',6.212296324000001e-05],
+	['Fe57',1.4346901400000002e-06],
+	['Fe58',1.9093092e-07],
+	['Au197',0.0003014]])
+mvoid = mcdc.material(nuclides=[['N14',1e-10]])
+##############################
+#__________Geometry__________
+##############################
+
+# Surface(s)
+s1 = mcdc.surface('plane-z', z=-1.27, bc='vacuum')
+s2 = mcdc.surface('plane-z', z=0.0, bc='interface')
+s3 = mcdc.surface('plane-z', z=7.42, bc='interface')
+s4 = mcdc.surface('plane-z', z=7.5788, bc='interface')
+s5 = mcdc.surface('plane-z', z=7.5998, bc='interface')
+s6 = mcdc.surface('plane-z', z=9.6089, bc='interface')
+s7 = mcdc.surface('plane-z', z=9.6299, bc='interface')
+s8 = mcdc.surface('plane-z', z=9.6553, bc='interface')
+s9 = mcdc.surface('plane-z', z=12.6274, bc='interface')
+s10 = mcdc.surface('plane-z', z=12.6528, bc='interface')
+s11 = mcdc.surface('plane-z', z=14.8959, bc='interface')
+s12 = mcdc.surface('plane-z', z=14.9213, bc='interface')
+s13 = mcdc.surface('plane-z', z=17.9219, bc='interface')
+s14 = mcdc.surface('plane-z', z=17.9473, bc='interface')
+s15 = mcdc.surface('plane-z', z=17.9683, bc='interface')
+s16 = mcdc.surface('plane-z', z=19.9808, bc='interface')
+s17 = mcdc.surface('plane-z', z=20.0018, bc='interface')
+s18 = mcdc.surface('plane-z', z=20.1606, bc='interface')
+s19 = mcdc.surface('plane-z', z=27.5806, bc='interface')
+s20 = mcdc.surface('plane-z', z=30.0685, bc='vacuum')
+s21 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=6.0325, bc='interface')
+s22 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=6.35, bc='interface')
+s23 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=7.505, bc='interface')
+s24 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=7.5438, bc='interface')
+s25 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=7.5489, bc='interface')
+s26 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=7.6759, bc='interface')
+s27 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=7.7013, bc='interface')
+s28 = mcdc.surface('cylinder-z', center=[0.0, 0.0], radius=15.1409, bc='vacuum')
+
+# Material Cell(s)
+c1 = mcdc.cell(+s1 & -s2 & -s28, fill=m1)
+c2 = mcdc.cell(+s2 & -s3 & -s24, fill=m2)
+c3 = mcdc.cell(+s3 & -s4 & -s23, fill=m3)
+c4 = mcdc.cell(+s3 & -s4 & +s23 & -s24, fill=mvoid)
+c5 = mcdc.cell(+s4 & -s5 & -s24, fill=m3)
+c6 = mcdc.cell(+s5 & -s6 & -s23, fill=m15)
+c7 = mcdc.cell(+s5 & -s6 & +s23 & -s24, fill=m3)
+c8 = mcdc.cell(+s6 & -s7 & -s24, fill=m3)
+c9 = mcdc.cell(+s7 & -s8 & -s24, fill=m3)
+c10 = mcdc.cell(+s8 & -s9 & -s23, fill=m12)
+c11 = mcdc.cell(+s8 & -s9 & +s23 & -s24, fill=m3)
+c12 = mcdc.cell(+s9 & -s10 & -s24, fill=m3)
+c13 = mcdc.cell(+s10 & -s11 & -s21, fill=mvoid)
+c14 = mcdc.cell(+s10 & -s11 & +s21 & -s22, fill=m1)
+c15 = mcdc.cell(+s10 & -s11 & +s22 & -s24, fill=mvoid)
+c16 = mcdc.cell(+s11 & -s12 & -s24, fill=m3)
+c17 = mcdc.cell(+s12 & -s13 & -s23, fill=m11)
+c18 = mcdc.cell(+s12 & -s13 & +s23 & -s24, fill=m3)
+c19 = mcdc.cell(+s13 & -s14 & -s24, fill=m3)
+c20 = mcdc.cell(+s14 & -s15 & -s24, fill=m3)
+c21 = mcdc.cell(+s15 & -s16 & -s23, fill=m16)
+c22 = mcdc.cell(+s15 & -s16 & +s23 & -s24, fill=m3)
+c23 = mcdc.cell(+s16 & -s17 & -s24, fill=m3)
+c24 = mcdc.cell(+s17 & -s18 & -s23, fill=m3)
+c25 = mcdc.cell(+s17 & -s18 & +s23 & -s24, fill=mvoid)
+c26 = mcdc.cell(+s18 & -s19 & -s24, fill=m2)
+c27 = mcdc.cell(+s19 & -s20 & -s24, fill=mvoid)
+c28 = mcdc.cell(+s2 & -s20 & +s24 & -s25, fill=mvoid)
+c29 = mcdc.cell(+s2 & -s20 & +s25 & -s26, fill=m4)
+c30 = mcdc.cell(+s2 & -s20 & +s26 & -s27, fill=mvoid)
+c31 = mcdc.cell(+s2 & -s20 & +s27 & -s28, fill=m2)
+# Root Universe Cells List:
+u0_cells = []
+# Material Universe(s)
+
+##############################
+#__________Settings__________
+##############################
+
+# Simulation Parameters
+mcdc.eigenmode(N_inactive=20, N_active=180)
+mcdc.setting(N_particle=10000)
+
+# Source Parameters
+# Particle: neutron
+# Space Type: box
+mcdc.source(x=[-1.0,1.0], y=[-1.0,1.0], z=[9.8803,10.4322], prob=1.0)
+
+mcdc.run()
+