diff --git a/Model/Analysis/Analysis.ipynb b/Model/Analysis/Analysis.ipynb new file mode 100644 index 0000000..5205bc1 --- /dev/null +++ b/Model/Analysis/Analysis.ipynb @@ -0,0 +1,517 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "302f22f9-375f-4773-a32b-c51cb16bbb90", + "metadata": {}, + "source": [ + "# Post Processing, Visualizing $k_{eff}$ Convergence and Shannon Entropy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c2e26fe9-3689-4f85-8e8e-da3867961bb3", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 1) Load statepoint and basic parameters\n", + "# ----------------------------------------------------------------------------\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "n_inactive = sp.n_inactive\n", + "n_batches = sp.n_batches\n", + "n_active = n_batches - n_inactive\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 2) Generation-based data: raw k_generation + Shannon entropy\n", + "# ----------------------------------------------------------------------------\n", + "k_eff_raw = np.array(sp.k_generation) # one k-eff per generation\n", + "H = np.array(sp.entropy) # one Shannon entropy per generation\n", + "\n", + "n_generations = len(k_eff_raw)\n", + "gens = np.arange(1, n_generations + 1)\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 3) Batch-based data: average each batch from k_generation\n", + "# ----------------------------------------------------------------------------\n", + "# Determine how many generations per batch (assumes even division)\n", + "gens_per_batch = n_generations // n_batches\n", + "\n", + "# Average each row => one k value per batch\n", + "k_per_batch = np.mean(k_eff_raw.reshape(n_batches, gens_per_batch), axis=1)\n", + "\n", + "# Compute running average over active batches\n", + "batch_mean = np.full(n_batches, np.nan)\n", + "batch_stderr = np.full(n_batches, np.nan)\n", + "\n", + "for i in range(n_batches):\n", + " if i < n_inactive:\n", + " # Skip inactive\n", + " continue\n", + " else:\n", + " active_slice = k_per_batch[n_inactive : i+1]\n", + " batch_mean[i] = np.mean(active_slice)\n", + " if len(active_slice) > 1:\n", + " std_sample = np.std(active_slice, ddof=1)\n", + " batch_stderr[i] = std_sample / np.sqrt(len(active_slice))\n", + " else:\n", + " batch_stderr[i] = 0.0\n", + "\n", + "# Plot each batch point at the generation index where the batch ends\n", + "batch_x = (np.arange(n_batches) + 1) * gens_per_batch\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 4) Create the figure and plot\n", + "# ----------------------------------------------------------------------------\n", + "fig, ax1 = plt.subplots()\n", + "\n", + "# (a) Plot raw generation-based k (no error bars)\n", + "line1, = ax1.plot(\n", + " gens, k_eff_raw, 'o-', label='Gen-based k', alpha=0.7\n", + ")\n", + "\n", + "# (b) Plot batch-based running average k (with error bars)\n", + "line2 = ax1.errorbar(\n", + " batch_x, batch_mean, yerr=batch_stderr,\n", + " fmt='s-', capsize=3, label='Batch-based k', alpha=0.7\n", + ")\n", + "\n", + "# Left y-axis for k\n", + "ax1.set_xlabel('Generation')\n", + "ax1.set_ylabel('k-effective', color='tab:blue')\n", + "ax1.tick_params(axis='y', labelcolor='tab:blue')\n", + "ax1.grid(True)\n", + "\n", + "# (c) Plot Shannon entropy on a second y-axis\n", + "ax2 = ax1.twinx()\n", + "line3, = ax2.plot(\n", + " gens, H, 'r-', label='Shannon entropy'\n", + ")\n", + "ax2.set_ylabel('Shannon Entropy', color='tab:red')\n", + "ax2.tick_params(axis='y', labelcolor='tab:red')\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 5) Combine legends and place them on the middle right\n", + "# ----------------------------------------------------------------------------\n", + "lines1, labels1 = ax1.get_legend_handles_labels()\n", + "lines2, labels2 = ax2.get_legend_handles_labels()\n", + "\n", + "# Merge them into a single legend\n", + "ax2.legend(lines1 + lines2, labels1 + labels2, loc='center right')\n", + "\n", + "plt.title('Gen-based k, Batch-based k, and Shannon Entropy')\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "bc3342e7-ace6-43d8-ace4-6a250c21336f", + "metadata": {}, + "source": [ + "# Post Processing, Visualizing Flux and Fission" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71db4e8b-d18a-4e02-9789-981b45eee8f7", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from ipywidgets import interact, IntSlider\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Function to get data based on tally and score\n", + "def get_tally_data(tally_name, score='fission'):\n", + " tally = sp.get_tally(name=tally_name)\n", + " return tally.get_values(scores=[score])\n", + "\n", + "# Retrieve the fission data (change tally name and score if needed)\n", + "data_type = 'fission' # Change this to 'flux', 'fission', or any other score available\n", + "plot_title = 'Fission Events' # Adjust title for visualization\n", + "tally_name = 'Fission' # Adjust the tally name if necessary\n", + "data = get_tally_data(tally_name=tally_name, score=data_type)\n", + "\n", + "# Define the mesh size and reshape the data\n", + "total_data_points = 300**3\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "data_reshaped = data.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Determine the data range for consistent color scaling\n", + "data_min = data_reshaped.min()\n", + "data_max = data_reshaped.max()\n", + "\n", + "# Function to display a single slice\n", + "def show_slice(z_slice_index):\n", + " plt.figure(figsize=(8, 6))\n", + " data_2d = data_reshaped[z_slice_index, :, :]\n", + "\n", + " # Generate the heatmap for the x-y slice with fixed color scale\n", + " plt.imshow(data_2d, origin='lower', cmap='nipy_spectral', vmin=data_min, vmax=data_max)\n", + " plt.colorbar(label='Fission Events')\n", + " plt.title(f'Fission Events in the x-y plane at z index {z_slice_index}')\n", + " plt.xlabel('X-axis')\n", + " plt.ylabel('Y-axis')\n", + " plt.show()\n", + "\n", + "# Create an interactive slider\n", + "interact(show_slice, z_slice_index=IntSlider(min=0, max=mesh_dimension_size-1, step=1, value=0));\n", + "\n", + "#Up the contrast of the color bars to make it nicer to look at. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cba77ca7-dbda-44c5-904b-2a0856b751d8", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from ipywidgets import interact, IntSlider\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Function to get data based on tally and score\n", + "def get_tally_data(tally_name, score='flux'):\n", + " tally = sp.get_tally(name=tally_name)\n", + " return tally.get_values(scores=[score])\n", + "\n", + "# Retrieve the fission data (change tally name and score if needed)\n", + "data_type = 'flux' # Change this to 'flux', 'fission', or any other score available\n", + "plot_title = 'Neutron Flux' # Adjust title for visualization\n", + "tally_name = 'Flux' # Adjust the tally name if necessary\n", + "data = get_tally_data(tally_name=tally_name, score=data_type)\n", + "\n", + "# Define the mesh size and reshape the data\n", + "total_data_points = 300**3\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "data_reshaped = data.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Determine the data range for consistent color scaling\n", + "data_min = data_reshaped.min()\n", + "data_max = data_reshaped.max()\n", + "\n", + "# Function to display a single slice\n", + "def show_slice(z_slice_index):\n", + " plt.figure(figsize=(8, 6))\n", + " data_2d = data_reshaped[z_slice_index, :, :]\n", + "\n", + " # Generate the heatmap for the x-y slice with fixed color scale\n", + " plt.imshow(data_2d, origin='lower', cmap='nipy_spectral', vmin=data_min, vmax=data_max)\n", + " plt.colorbar(label='Neutron Flux')\n", + " plt.title(f'{plot_title} in the x-y plane at z index {z_slice_index}')\n", + " plt.xlabel('X-axis')\n", + " plt.ylabel('Y-axis')\n", + " plt.show()\n", + "\n", + "# Create an interactive slider\n", + "interact(show_slice, z_slice_index=IntSlider(min=0, max=mesh_dimension_size-1, step=1, value=0));\n" + ] + }, + { + "cell_type": "markdown", + "id": "d841bbe8-700d-40b5-86ca-c66a66fc978d", + "metadata": {}, + "source": [ + "# Create a gif File Showing the Fission or Flux" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a879a2f-894a-4c89-aefa-cf6e4f2ffd0c", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from matplotlib.animation import FuncAnimation, PillowWriter\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Function to get data based on tally and score\n", + "def get_tally_data(tally_name, score='flux'):\n", + " tally = sp.get_tally(name=tally_name)\n", + " return tally.get_values(scores=[score])\n", + "\n", + "# Retrieve the fission data (change tally name and score if needed)\n", + "data_type = 'flux' # Change this to 'flux', 'fission', or any other score available\n", + "plot_title = 'Neutron Flux' # Adjust title for visualization\n", + "tally_name = 'Flux' # Adjust the tally name if necessary\n", + "data = get_tally_data(tally_name=tally_name, score=data_type)\n", + "\n", + "# Define the mesh size and reshape the data\n", + "total_data_points = 300**3\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "data_reshaped = data.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Determine the data range for consistent color scaling\n", + "data_min = data_reshaped.min()\n", + "data_max = data_reshaped.max()\n", + "\n", + "# Initialize the figure and axis\n", + "fig, ax = plt.subplots()\n", + "\n", + "# Function to update the animation for each frame (z-slice)\n", + "def update(z_slice_index, plot_title):\n", + " ax.clear() # Clear previous frame\n", + " data_2d = data_reshaped[z_slice_index, :, :]\n", + " \n", + " # Generate the heatmap for the x-y slice\n", + " im = ax.imshow(data_2d, origin='lower', cmap='viridis', vmin=data_min, vmax=data_max)\n", + " ax.set_title(f'{plot_title} in the x-y plane at z index {z_slice_index}')\n", + " ax.set_xlabel('X-axis')\n", + " ax.set_ylabel('Y-axis')\n", + " return im,\n", + "\n", + "custom_title = \"Neutron Flux\"\n", + "\n", + "# Create the animation\n", + "num_frames = mesh_dimension_size # One frame per z-slice\n", + "ani = FuncAnimation(fig, update, frames=num_frames, fargs=(custom_title,), blit=True, interval=100)\n", + "\n", + "# Save the animation as a video\n", + "ani.save('neutron_flux_animation.gif', writer='ffmpeg', fps=10)\n", + "\n", + "# Optionally display the animation as HTML\n", + "from IPython.display import HTML\n", + "HTML(ani.to_jshtml())\n", + "\n", + "# Close the plot to clean up\n", + "plt.close()" + ] + }, + { + "cell_type": "markdown", + "id": "af7e2d8d-b0d7-49f9-b684-3ac8321b78cb", + "metadata": {}, + "source": [ + "# Visualizing Errors" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fab13cce-2c0a-4cdb-a0ba-9471043badcf", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Access the flux tally\n", + "flux_tally = sp.get_tally(name='Flux')\n", + "\n", + "# Get the mean and standard deviation for the flux tally\n", + "flux_mean = flux_tally.get_values(scores=['flux'], value='mean')\n", + "flux_std_dev = flux_tally.get_values(scores=['flux'], value='std_dev')\n", + "\n", + "# Total number of data points\n", + "total_data_points = 300**3\n", + "\n", + "# Since the data is from a cubic mesh, we take the cube root to find the size of one dimension\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "\n", + "# Reshape the mean and std_dev data to the mesh dimensions\n", + "flux_mean_reshaped = flux_mean.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "flux_std_dev_reshaped = flux_std_dev.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Compute the relative error\n", + "relative_error = np.zeros_like(flux_mean_reshaped)\n", + "nonzero = flux_mean_reshaped > 0\n", + "relative_error[nonzero] = flux_std_dev_reshaped[nonzero] / flux_mean_reshaped[nonzero]\n", + "\n", + "# Plot the distribution of relative errors\n", + "plt.hist(relative_error[nonzero], bins=500)\n", + "plt.title(\"Distribution of Relative Errors (Flux)\")\n", + "plt.xlabel(\"Relative Error\")\n", + "plt.ylabel(\"Frequency\")\n", + "plt.show()\n", + "\n", + "# Access the fission tally\n", + "fission_tally = sp.get_tally(name='Fission')\n", + "\n", + "# Get the mean and standard deviation for the fission tally\n", + "fission_mean = fission_tally.get_values(scores=['fission'], value='mean')\n", + "fission_std_dev = fission_tally.get_values(scores=['fission'], value='std_dev')\n", + "\n", + "# Reshape the mean and std_dev data to the mesh dimensions\n", + "fission_mean_reshaped = fission_mean.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "fission_std_dev_reshaped = fission_std_dev.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Compute the relative error\n", + "relative_error_fission = np.zeros_like(fission_mean_reshaped)\n", + "nonzero_fission = fission_mean_reshaped > 0\n", + "relative_error_fission[nonzero_fission] = fission_std_dev_reshaped[nonzero_fission] / fission_mean_reshaped[nonzero_fission]\n", + "\n", + "# Plot the distribution of relative errors for fission\n", + "plt.hist(relative_error_fission[nonzero_fission], bins=500)\n", + "plt.title(\"Distribution of Relative Errors (Fission)\")\n", + "plt.xlabel(\"Relative Error\")\n", + "plt.ylabel(\"Frequency\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "df721e47-fc3f-4bd9-9ea7-afcd253636d5", + "metadata": {}, + "source": [ + "# Plotting Heat Map of Errors" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b79d437-5b20-45da-b010-5daf751bcc23", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from ipywidgets import interact, IntSlider\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Access the tallies\n", + "flux_tally = sp.get_tally(name='Flux')\n", + "fission_tally = sp.get_tally(name='Fission')\n", + "\n", + "# Extract mean and standard deviation for flux and fission\n", + "flux_data_mean = flux_tally.get_values(scores=['flux'], value='mean')\n", + "fission_data_mean = fission_tally.get_values(scores=['fission'], value='mean')\n", + "flux_data_std_dev = flux_tally.get_values(scores=['flux'], value='std_dev')\n", + "fission_data_std_dev = fission_tally.get_values(scores=['fission'], value='std_dev')\n", + "\n", + "# Define the mesh dimensions (assumed cubic with 300 points in each dimension)\n", + "data_points = 300\n", + "\n", + "# Reshape the 1D arrays into 3D volumes\n", + "flux_data_mean = flux_data_mean.reshape((data_points, data_points, data_points))\n", + "fission_data_mean = fission_data_mean.reshape((data_points, data_points, data_points))\n", + "flux_data_std_dev = flux_data_std_dev.reshape((data_points, data_points, data_points))\n", + "fission_data_std_dev = fission_data_std_dev.reshape((data_points, data_points, data_points))\n", + "\n", + "# Calculate the relative errors (std_dev / mean)\n", + "Relative_errors_flux = flux_data_std_dev / flux_data_mean\n", + "Relative_errors_fission = fission_data_std_dev / fission_data_mean\n", + "\n", + "# Convert to percentage\n", + "flux_error_percent = Relative_errors_flux * 100\n", + "fission_error_percent = Relative_errors_fission * 100\n", + "\n", + "vmin = 0\n", + "vmax = 100\n", + "\n", + "# Determine a common color scale range for both images\n", + "data_min = min(Relative_errors_flux.min(), Relative_errors_fission.min())\n", + "data_max = max(Relative_errors_flux.max(), Relative_errors_fission.max())\n", + "\n", + "def show_slice(z_slice_index):\n", + " # Create a figure with two subplots side by side\n", + " fig, axes = plt.subplots(1, 2, figsize=(12, 6))\n", + " \n", + " # Rotate the slices to orient them correctly for display\n", + " flux_slice = np.rot90(flux_error_percent[z_slice_index, :, :], 4)\n", + " fission_slice = np.rot90(fission_error_percent[z_slice_index, :, :], 4)\n", + " \n", + " # Plot flux relative error\n", + " im1 = axes[0].imshow(flux_slice, origin='lower', cmap='nipy_spectral', vmin=vmin, vmax=vmax)\n", + " axes[0].set_title(f'Flux Relative Error (Slice {z_slice_index})')\n", + " axes[0].set_xlabel('X-axis')\n", + " axes[0].set_ylabel('Y-axis')\n", + " plt.colorbar(im1, ax=axes[0])\n", + " \n", + " # Plot fission relative error\n", + " im2 = axes[1].imshow(fission_slice, origin='lower', cmap='nipy_spectral', vmin=vmin, vmax=vmax)\n", + " axes[1].set_title(f'Fission Relative Error (Slice {z_slice_index})')\n", + " axes[1].set_xlabel('X-axis')\n", + " axes[1].set_ylabel('Y-axis')\n", + " plt.colorbar(im2, ax=axes[1])\n", + " \n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "# Create an interactive slider to select the z-slice\n", + "interact(show_slice, z_slice_index=IntSlider(min=0, max=data_points-1, step=1, value=data_points//2));\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7d59a68-247a-481a-a38b-33838708545c", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2632fda8-84df-4c67-9c14-da8c19992268", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30f7b040-1f69-4987-9706-21557518e5ba", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f83ccbd-e5b8-4296-9ea7-940f57f98712", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Model/RadCenterSubcriticalAssemblyOSU.ipynb b/Model/RadCenterSubcriticalAssemblyOSU.ipynb new file mode 100644 index 0000000..8ba34f3 --- /dev/null +++ b/Model/RadCenterSubcriticalAssemblyOSU.ipynb @@ -0,0 +1,994 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "44c829d0-579c-4033-9b28-13e10139156e", + "metadata": {}, + "source": [ + "# Create material files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a6fdd011-90d7-4773-b42d-53d9716ca4a5", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "################################################################################\n", + "#Materials\n", + "################################################################################\n", + "#Plutonium (PuBe)\n", + "PuBe = openmc.Material(name = 'PuBe')\n", + "PuBe.add_nuclide('Pu239', 1.0, 'ao') \n", + "PuBe.add_element('Be', 13.0, 'ao') \n", + "PuBe.set_density('g/cm3', 15.0) \n", + "\n", + "################################################################################\n", + "#Aluminum Cladding\n", + "aluminum_cladding = openmc.Material(name = 'Aluminum Cladding')\n", + "aluminum_cladding.add_element('Al', 1.0)\n", + "aluminum_cladding.set_density('g/cm3', 2.7)\n", + "\n", + "################################################################################\n", + "#Fuel (UO2)\n", + "fuel = openmc.Material(name='UO2 Fuel')\n", + "fuel.add_nuclide('U238', 0.99284)\n", + "fuel.add_nuclide('U235', 0.00711)\n", + "fuel.add_nuclide('U234', 0.00005)\n", + "fuel.add_s_alpha_beta('c_O_in_UO2')\n", + "fuel.add_element('O', 2)\n", + "fuel.set_density('g/cm3', 10.5)\n", + "\n", + "################################################################################\n", + "#Graphite Moderator\n", + "graphite = openmc.Material(name='Graphite')\n", + "graphite.add_element('C', 1.0)\n", + "graphite.set_density('g/cm3', 1.7)\n", + "graphite.add_s_alpha_beta('c_Graphite')\n", + "\n", + "################################################################################\n", + "#Air\n", + "air = openmc.Material(name='Air')\n", + "air.add_element('N', 0.78084)\n", + "air.add_element('O', 0.209476)\n", + "air.add_element('Ar', 0.00934)\n", + "air.set_density('g/cm3', 0.001204)\n", + "\n", + "################################################################################\n", + "#Concrete\n", + "concrete = openmc.Material(name='Concrete')\n", + "concrete.add_element('H', 0.149857)\n", + "concrete.add_element('C', 0.074206)\n", + "concrete.add_element('O', 0.526836)\n", + "concrete.add_element('Mg', 0.017713)\n", + "concrete.add_element('Al', 0.023794)\n", + "concrete.add_element('Si', 0.091975)\n", + "concrete.add_element('S', 0.001649)\n", + "concrete.add_element('K', 0.000773)\n", + "concrete.add_element('Ca', 0.109681)\n", + "concrete.add_element('Fe', 0.003516)\n", + "concrete.set_density('g/cm3', 2.35)\n", + "\n", + "################################################################################\n", + "#Export Materials\n", + "materials = openmc.Materials([fuel, PuBe, aluminum_cladding, graphite, air, concrete])\n", + "materials.export_to_xml()" + ] + }, + { + "cell_type": "markdown", + "id": "9d81da67-1b43-4145-9a43-9ffa18e37b84", + "metadata": {}, + "source": [ + "# Create Main Dimensions of Reactor" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17224e24-7392-44c4-a29a-4fe93da80f7a", + "metadata": {}, + "outputs": [], + "source": [ + "################################################################################\n", + "#Brick dimensions in the x-y-z planes\n", + "brick_side_length_x = 15\n", + "brick_side_length_y = 15\n", + "brick_side_length_z = 50.0\n", + "\n", + "################################################################################\n", + "# Radii\n", + "r_air_hole = 0.46\n", + "r_inner_clad = 0.61875\n", + "r_fuel = 1.76125\n", + "r_outer_clad = 1.92\n", + "\n", + "# Lengths\n", + "total_length = 135.0\n", + "air_gap_back = 15.0\n", + "\n", + "#Dimensions in total:\n", + "x_direction = 165.24\n", + "y_direction = 150\n", + "z_direction = 150" + ] + }, + { + "cell_type": "markdown", + "id": "cd3157db-2425-4998-840b-e696b5f1a687", + "metadata": {}, + "source": [ + "# Create Reactor Boundary and Fuel Universe" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92959fa3-1c6d-4eb0-8296-c85f2fac0988", + "metadata": {}, + "outputs": [], + "source": [ + "# Reactor Boundary Conditions\n", + "x_min = openmc.XPlane(x0=0, boundary_type='vacuum')\n", + "x_max = openmc.XPlane(x0=165.24, boundary_type='vacuum')\n", + "y_min = openmc.YPlane(y0=-30, boundary_type='vacuum')\n", + "y_max = openmc.YPlane(y0=150.0, boundary_type='vacuum')\n", + "z_min = openmc.ZPlane(z0=0, boundary_type='vacuum')\n", + "z_max = openmc.ZPlane(z0=150.0, boundary_type='vacuum')\n", + "\n", + "# Combine the boundary planes to form a cuboid region\n", + "reactor_boundary = (+x_min & -x_max & +y_min & -y_max & +z_min & -z_max)\n", + "\n", + "############################################################################\n", + "# Make Fuel Universe\n", + "\n", + "# Define the fuel region with a single cylinder per hole\n", + "z_plane_start = openmc.ZPlane(z0=0)\n", + "z_plane_air = openmc.ZPlane(z0=total_length)\n", + "z_plane_end = openmc.ZPlane(z0=z_direction)\n", + "\n", + "# Define cylinders for air hole, cladding, and fuel\n", + "air_hole = openmc.ZCylinder(r=r_air_hole, boundary_type='transmission')\n", + "inner_clad = openmc.ZCylinder(r=r_inner_clad, boundary_type='transmission')\n", + "fuel_rod = openmc.ZCylinder(r=r_fuel, boundary_type='transmission')\n", + "outer_clad = openmc.ZCylinder(r=r_outer_clad, boundary_type='transmission')\n", + "\n", + "# Define cells for air hole, cladding, and fuel\n", + "air_hole_cell = openmc.Cell(name='Air Hole', fill=air)\n", + "air_hole_cell.region = -air_hole & +z_plane_start & -z_plane_air\n", + "\n", + "inner_clad_cell = openmc.Cell(name='Inner Cladding', fill=aluminum_cladding)\n", + "inner_clad_cell.region = +air_hole & -inner_clad & +z_plane_start & -z_plane_air\n", + "\n", + "fuel_cell = openmc.Cell(name='Fuel Rod', fill=fuel)\n", + "fuel_cell.region = +inner_clad & -fuel_rod & +z_plane_start & -z_plane_air\n", + "\n", + "outer_clad_cell = openmc.Cell(name='Outer Cladding', fill=aluminum_cladding)\n", + "outer_clad_cell.region = +fuel_rod & -outer_clad & +z_plane_start & -z_plane_air\n", + "\n", + "#Don't forget about the 15 cm air gap on the back\n", + "air_cylinder = openmc.ZCylinder(r=r_outer_clad, boundary_type='transmission')\n", + "air_cell = openmc.Cell(name='Air Gap', fill=air)\n", + "air_cell.region = -air_cylinder & +z_plane_air & -z_plane_end \n", + "\n", + "pitch = 15.0\n", + "#Graphite Bricks\n", + "graphite_bricks = openmc.model.rectangular_prism(width=pitch, height=pitch)\n", + "graphite_bricks_cell = openmc.Cell(name='Graphite_bricks')\n", + "graphite_bricks_cell.fill = graphite\n", + "graphite_bricks_cell.region = graphite_bricks & +air_cylinder & +outer_clad & -z_plane_end & +z_min\n", + "\n", + "fuel_rod_universe = openmc.Universe(cells=[air_hole_cell, inner_clad_cell, fuel_cell, outer_clad_cell, air_cell, graphite_bricks_cell])" + ] + }, + { + "cell_type": "markdown", + "id": "5c4d973a-980e-486e-8595-c89a1a5089d8", + "metadata": {}, + "source": [ + "### Visualize Fuel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a932512-5fde-45f5-9a7f-6b3e115e8345", + "metadata": {}, + "outputs": [], + "source": [ + "fuel_rod_universe.plot(width=(20,20), colors = {fuel_cell:'green', air_hole_cell:'blue', graphite_bricks_cell:'black', inner_clad_cell:'grey', outer_clad_cell:'red', air_cell:'yellow'}, basis=('xy'))" + ] + }, + { + "cell_type": "markdown", + "id": "350a1665-e571-4c42-b264-d3eac3afe383", + "metadata": {}, + "source": [ + "# Create Fuel Lattice Structure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c68f8666-3ba7-4a06-8104-f976ee44e9ea", + "metadata": {}, + "outputs": [], + "source": [ + "quarter_pitch = pitch*10\n", + "fuel_universe_lattice = openmc.RectLattice(name='Fuel Assembly')\n", + "fuel_universe_lattice.pitch = (pitch, pitch)\n", + "fuel_universe_lattice.lower_left = (0, 0) # Centered in a 150x150 box\n", + "\n", + "\n", + "# Make a 10x10 lattice for the fuel\n", + "fuel_universe_lattice.universes = [\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe],\n", + " [fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe, fuel_rod_universe]\n", + "]\n", + "\n", + "fuel_assembly_region = openmc.model.rectangular_prism(width=quarter_pitch, height=quarter_pitch, origin=(75, 75))\n", + "fuel_assembly_cell = openmc.Cell(name='quarter assembly cell', fill=fuel_universe_lattice, region=fuel_assembly_region)\n", + "\n", + "quarter_assembly_universe = openmc.Universe(cells=[fuel_assembly_cell])" + ] + }, + { + "cell_type": "markdown", + "id": "937a879f-8afd-4899-a0ca-30cac32a7b19", + "metadata": {}, + "source": [ + "### Visualize the Lattice Structure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c083bee0-7505-497a-a653-f6ab4d234906", + "metadata": {}, + "outputs": [], + "source": [ + "quarter_assembly_universe.plot(width=(200,200), colors = {fuel_cell:'green', air_hole_cell:'blue', graphite_bricks_cell:'black', inner_clad_cell:'grey', outer_clad_cell:'red', air_cell:'yellow'}, origin=(75.0, 75.0, 75.0), basis=('xy'))" + ] + }, + { + "cell_type": "markdown", + "id": "3e0ba83c-696e-4feb-a1ff-37e2ec549155", + "metadata": {}, + "source": [ + "# Create Graphite Floor and Plutonium Source" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d78c17e9-0551-4a41-be70-60398663fd33", + "metadata": {}, + "outputs": [], + "source": [ + "plutonium_region = openmc.Sphere(r=2.0, x0=75.0, y0=-15.0, z0=75.0)\n", + "plutonium_cell = openmc.Cell(fill=PuBe, region=-plutonium_region)\n", + "\n", + "graphite_floor_region = openmc.model.RectangularParallelepiped(\n", + " xmin=0, xmax=150, ymin=-30, ymax=0, zmin=0, zmax=150\n", + ")\n", + "graphite_floor_cell = openmc.Cell(fill=graphite, region=-graphite_floor_region & +plutonium_region)\n", + "\n", + "graphite_floor_universe = openmc.Universe(cells=[plutonium_cell, graphite_floor_cell])" + ] + }, + { + "cell_type": "markdown", + "id": "a70b8c6a-9c00-4dab-adda-659918f674bb", + "metadata": {}, + "source": [ + "### Visualize the Graphite Floor and Plutonium Fuel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acf31a4e-44bb-4b0b-9447-7736aac70e83", + "metadata": {}, + "outputs": [], + "source": [ + "graphite_floor_universe.plot(width=(100,100), colors = {graphite_floor_cell:'black', air_cell:'yellow'}, basis=('xy'))" + ] + }, + { + "cell_type": "markdown", + "id": "f43fb189-44a1-4dd3-854c-a73b17b55391", + "metadata": {}, + "source": [ + "# Create Concrete Wall" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97526e66-44ea-4d6a-b46d-5b4baabb2955", + "metadata": {}, + "outputs": [], + "source": [ + "concrete_wall_region = openmc.model.RectangularParallelepiped(\n", + " xmin=150, xmax=165.24, ymin=-30, ymax=150, zmin=0, zmax=150\n", + ")\n", + "concrete_wall_cell = openmc.Cell(fill=concrete, region=-concrete_wall_region)\n", + "\n", + "concrete_wall_universe = openmc.Universe(cells=[concrete_wall_cell])" + ] + }, + { + "cell_type": "markdown", + "id": "b571239f-8793-4817-8ac1-d2087bc4339c", + "metadata": {}, + "source": [ + "### Visualize the Concrete Wall" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "caa6ac71-0715-4e1f-8c53-52c68310eba7", + "metadata": {}, + "outputs": [], + "source": [ + "concrete_wall_universe.plot(width=(100,100), colors = {concrete_wall_cell:'grey'}, basis=('xy'))" + ] + }, + { + "cell_type": "markdown", + "id": "6ae48851-126d-4ed6-b7e9-1ccc0c1f2404", + "metadata": {}, + "source": [ + "# Visualizing the Whole Geometry and Combining Universes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0239edb-30ab-4b71-a9bc-c2a5c9339e4f", + "metadata": {}, + "outputs": [], + "source": [ + "# Combine All Universes\n", + "universe_cells = [fuel_assembly_cell, plutonium_cell, graphite_floor_cell, concrete_wall_cell]\n", + "root_universe = openmc.Universe(cells=[openmc.Cell(fill=openmc.Universe(cells=universe_cells), region=reactor_boundary)])\n", + "\n", + "geometry = openmc.Geometry(root_universe)\n", + "geometry.export_to_xml()\n", + "\n", + "material_colors = {\n", + " PuBe: 'yellow', # Plutonium-Beryllium in yellow\n", + " aluminum_cladding: 'grey', # Aluminum Cladding in grey\n", + " fuel: 'green', # Fuel (UO2) in green\n", + " graphite: 'black', # Graphite in black\n", + " air: 'blue', # Air in blue\n", + " concrete: 'brown' # Concrete in brown\n", + "}\n", + "\n", + "plot_xy = openmc.Plot()\n", + "plot_xy.basis = 'xy'\n", + "plot_xy.origin = (75.0, 90.0, 75.0)\n", + "plot_xy.width = (300., 300.)\n", + "plot_xy.pixels = (4000, 4000)\n", + "plot_xy.color_by = 'material'\n", + "plot_xy.colors = material_colors\n", + "\n", + "plot_xz = openmc.Plot()\n", + "plot_xz.basis = 'xz'\n", + "plot_xz.origin = (75.0, 90.0, 75.0)\n", + "plot_xz.width = (150., 150.)\n", + "plot_xz.pixels = (4000, 4000)\n", + "plot_xz.color_by = 'material'\n", + "plot_xz.colors = material_colors\n", + "\n", + "plot_yz = openmc.Plot()\n", + "plot_yz.basis = 'yz'\n", + "plot_yz.origin = (75.0, 90.0, 75.0)\n", + "plot_yz.width = (150., 150.)\n", + "plot_yz.pixels = (4000, 4000)\n", + "plot_yz.color_by = 'material'\n", + "plot_yz.colors = material_colors\n", + "\n", + "plots = openmc.Plots([plot_xy, plot_xz, plot_yz])\n", + "plots.export_to_xml()\n", + "openmc.plot_geometry()" + ] + }, + { + "cell_type": "markdown", + "id": "f921c625-886a-4e74-98e0-a1a645583661", + "metadata": {}, + "source": [ + "# Simulation Settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45a00f2d-f870-48d9-b1b3-fe962213d5e0", + "metadata": {}, + "outputs": [], + "source": [ + "# Make the initial guess start uniformly within the fuel lattice\n", + "x_min, x_max = 0.0, 150.0\n", + "y_min, y_max = -30.0, 150.0\n", + "z_min, z_max = 0.0, 150.0\n", + "\n", + "x_distribution = openmc.stats.Uniform(x_min, x_max)\n", + "y_distribution = openmc.stats.Uniform(y_min, y_max)\n", + "z_distribution = openmc.stats.Uniform(z_min, z_max)\n", + "\n", + "spatial_distribution = openmc.stats.CartesianIndependent(x_distribution, y_distribution, z_distribution)\n", + "\n", + "source = openmc.IndependentSource()\n", + "source.space = spatial_distribution\n", + "source.angle = openmc.stats.Isotropic()\n", + "source.strength = 1e8\n", + "\n", + "################################################################################\n", + "# Setup simulation settings\n", + "settings = openmc.Settings()\n", + "settings.source = source\n", + "settings.batches = 15\n", + "settings.inactive = 5\n", + "settings.particles = 10000000\n", + "settings.export_to_xml()\n", + "settings.create_fission_neutrons = True\n", + "################################################################################\n", + "# Define a mesh that spans the x-y-z planes\n", + "mesh = openmc.RegularMesh()\n", + "mesh.dimension = [300, 300, 300]\n", + "mesh.lower_left = [0.0, -30.0, 0.0]\n", + "mesh.upper_right = [150.0, 150.0, 150.0]" + ] + }, + { + "cell_type": "markdown", + "id": "6c9ab7d4-69cb-42bd-b7f9-ddbda2a1c94a", + "metadata": {}, + "source": [ + "# Running the Simulation, Flux, Fission, and Shannon Entropy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9bf6f0a0-a2b6-41d9-8299-b70b71eb4d77", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a mesh filter\n", + "mesh_filter = openmc.MeshFilter(mesh)\n", + "\n", + "################################################################################\n", + "# Define a tally for the flux\n", + "flux_tally = openmc.Tally(name='Flux')\n", + "flux_tally.filters = [mesh_filter]\n", + "flux_tally.scores = ['flux']\n", + "\n", + "################################################################################\n", + "# Define a tally for the fission\n", + "fission_tally = openmc.Tally(name='Fission')\n", + "fission_tally.filters = [mesh_filter]\n", + "fission_tally.scores = ['fission']\n", + "\n", + "################################################################################\n", + "# Add both tallies to a single tallies collection\n", + "tallies = openmc.Tallies()\n", + "tallies.append(flux_tally)\n", + "tallies.append(fission_tally)\n", + "\n", + "# Export tallies to XML\n", + "tallies.export_to_xml()\n", + "\n", + "################################################################################\n", + "# Create shannon entropy mesh\n", + "entropy_mesh = openmc.RegularMesh()\n", + "entropy_mesh.dimension = [8, 8, 8] # Adjust dimensions as needed\n", + "entropy_mesh.lower_left = [x_min, y_min, z_min]\n", + "entropy_mesh.upper_right = [x_max, y_max, z_max]\n", + "settings.entropy_mesh = entropy_mesh\n", + "\n", + "################################################################################\n", + "# Assemble and Run\n", + "model = openmc.model.Model(geometry, materials, settings, tallies=tallies)\n", + "model.run(threads=18, geometry_debug=False)\n" + ] + }, + { + "cell_type": "markdown", + "id": "bb7d585b-2196-470b-acd3-72ba7485a14a", + "metadata": {}, + "source": [ + "# Post Processing, Visualizing $k_{eff}$ Convergence and Shannon Entropy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05fad1ea-e31b-459f-8c4c-3816f640b0a6", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 1) Load statepoint and basic parameters\n", + "# ----------------------------------------------------------------------------\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "n_inactive = sp.n_inactive\n", + "n_batches = sp.n_batches\n", + "n_active = n_batches - n_inactive\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 2) Generation-based data: raw k_generation + Shannon entropy\n", + "# ----------------------------------------------------------------------------\n", + "k_eff_raw = np.array(sp.k_generation) # one k-eff per generation\n", + "H = np.array(sp.entropy) # one Shannon entropy per generation\n", + "\n", + "n_generations = len(k_eff_raw)\n", + "gens = np.arange(1, n_generations + 1)\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 3) Batch-based data: average each batch from k_generation\n", + "# ----------------------------------------------------------------------------\n", + "# Determine how many generations per batch (assumes even division)\n", + "gens_per_batch = n_generations // n_batches\n", + "\n", + "# Average each row => one k value per batch\n", + "k_per_batch = np.mean(k_eff_raw.reshape(n_batches, gens_per_batch), axis=1)\n", + "\n", + "# Compute running average over active batches\n", + "batch_mean = np.full(n_batches, np.nan)\n", + "batch_stderr = np.full(n_batches, np.nan)\n", + "\n", + "for i in range(n_batches):\n", + " if i < n_inactive:\n", + " # Skip inactive\n", + " continue\n", + " else:\n", + " active_slice = k_per_batch[n_inactive : i+1]\n", + " batch_mean[i] = np.mean(active_slice)\n", + " if len(active_slice) > 1:\n", + " std_sample = np.std(active_slice, ddof=1)\n", + " batch_stderr[i] = std_sample / np.sqrt(len(active_slice))\n", + " else:\n", + " batch_stderr[i] = 0.0\n", + "\n", + "# Plot each batch point at the generation index where the batch ends\n", + "batch_x = (np.arange(n_batches) + 1) * gens_per_batch\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 4) Create the figure and plot\n", + "# ----------------------------------------------------------------------------\n", + "fig, ax1 = plt.subplots()\n", + "\n", + "# (a) Plot raw generation-based k (no error bars)\n", + "line1, = ax1.plot(\n", + " gens, k_eff_raw, 'o-', label='Gen-based k', alpha=0.7\n", + ")\n", + "\n", + "# (b) Plot batch-based running average k (with error bars)\n", + "line2 = ax1.errorbar(\n", + " batch_x, batch_mean, yerr=batch_stderr,\n", + " fmt='s-', capsize=3, label='Batch-based k', alpha=0.7\n", + ")\n", + "\n", + "# Left y-axis for k\n", + "ax1.set_xlabel('Generation')\n", + "ax1.set_ylabel('k-effective', color='tab:blue')\n", + "ax1.tick_params(axis='y', labelcolor='tab:blue')\n", + "ax1.grid(True)\n", + "\n", + "# (c) Plot Shannon entropy on a second y-axis\n", + "ax2 = ax1.twinx()\n", + "line3, = ax2.plot(\n", + " gens, H, 'r-', label='Shannon entropy'\n", + ")\n", + "ax2.set_ylabel('Shannon Entropy', color='tab:red')\n", + "ax2.tick_params(axis='y', labelcolor='tab:red')\n", + "\n", + "# ----------------------------------------------------------------------------\n", + "# 5) Combine legends and place them on the middle right\n", + "# ----------------------------------------------------------------------------\n", + "lines1, labels1 = ax1.get_legend_handles_labels()\n", + "lines2, labels2 = ax2.get_legend_handles_labels()\n", + "\n", + "# Merge them into a single legend\n", + "ax2.legend(lines1 + lines2, labels1 + labels2, loc='center right')\n", + "\n", + "plt.title('Gen-based k, Batch-based k, and Shannon Entropy')\n", + "plt.tight_layout()\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "4c944ed9-679b-4ba9-b0db-3e356a1d3cee", + "metadata": {}, + "source": [ + "# Post Processing, Visualizing Flux and Fission" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43973f07-ffbc-411e-a575-53e086bb13c2", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from ipywidgets import interact, IntSlider\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Function to get data based on tally and score\n", + "def get_tally_data(tally_name, score='fission'):\n", + " tally = sp.get_tally(name=tally_name)\n", + " return tally.get_values(scores=[score])\n", + "\n", + "# Retrieve the fission data (change tally name and score if needed)\n", + "data_type = 'fission' # Change this to 'flux', 'fission', or any other score available\n", + "plot_title = 'Fission Events' # Adjust title for visualization\n", + "tally_name = 'Fission' # Adjust the tally name if necessary\n", + "data = get_tally_data(tally_name=tally_name, score=data_type)\n", + "\n", + "# Define the mesh size and reshape the data\n", + "total_data_points = 300**3\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "data_reshaped = data.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Determine the data range for consistent color scaling\n", + "data_min = data_reshaped.min()\n", + "data_max = data_reshaped.max()\n", + "\n", + "# Function to display a single slice\n", + "def show_slice(z_slice_index):\n", + " plt.figure(figsize=(8, 6))\n", + " data_2d = data_reshaped[z_slice_index, :, :]\n", + "\n", + " # Generate the heatmap for the x-y slice with fixed color scale\n", + " plt.imshow(data_2d, origin='lower', cmap='nipy_spectral', vmin=data_min, vmax=data_max)\n", + " plt.colorbar(label='Fission Events')\n", + " plt.title(f'Fission Events in the x-y plane at z index {z_slice_index}')\n", + " plt.xlabel('X-axis')\n", + " plt.ylabel('Y-axis')\n", + " plt.show()\n", + "\n", + "# Create an interactive slider\n", + "interact(show_slice, z_slice_index=IntSlider(min=0, max=mesh_dimension_size-1, step=1, value=0));\n", + "\n", + "#Up the contrast of the color bars to make it nicer to look at. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3ce0cc75-fec2-4ecf-98a1-a7c9e22fbe43", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from ipywidgets import interact, IntSlider\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Function to get data based on tally and score\n", + "def get_tally_data(tally_name, score='flux'):\n", + " tally = sp.get_tally(name=tally_name)\n", + " return tally.get_values(scores=[score])\n", + "\n", + "# Retrieve the fission data (change tally name and score if needed)\n", + "data_type = 'flux' # Change this to 'flux', 'fission', or any other score available\n", + "plot_title = 'Neutron Flux' # Adjust title for visualization\n", + "tally_name = 'Flux' # Adjust the tally name if necessary\n", + "data = get_tally_data(tally_name=tally_name, score=data_type)\n", + "\n", + "# Define the mesh size and reshape the data\n", + "total_data_points = 300**3\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "data_reshaped = data.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Determine the data range for consistent color scaling\n", + "data_min = data_reshaped.min()\n", + "data_max = data_reshaped.max()\n", + "\n", + "# Function to display a single slice\n", + "def show_slice(z_slice_index):\n", + " plt.figure(figsize=(8, 6))\n", + " data_2d = data_reshaped[z_slice_index, :, :]\n", + "\n", + " # Generate the heatmap for the x-y slice with fixed color scale\n", + " plt.imshow(data_2d, origin='lower', cmap='nipy_spectral', vmin=data_min, vmax=data_max)\n", + " plt.colorbar(label='Neutron Flux')\n", + " plt.title(f'{plot_title} in the x-y plane at z index {z_slice_index}')\n", + " plt.xlabel('X-axis')\n", + " plt.ylabel('Y-axis')\n", + " plt.show()\n", + "\n", + "# Create an interactive slider\n", + "interact(show_slice, z_slice_index=IntSlider(min=0, max=mesh_dimension_size-1, step=1, value=0));\n" + ] + }, + { + "cell_type": "markdown", + "id": "63859286-5698-47b2-b20c-c187f4d52496", + "metadata": {}, + "source": [ + "# Create a gif File Showing the Fission or Flux" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "376b34e7-e55b-474f-a7d3-7d40d6b3ad15", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from matplotlib.animation import FuncAnimation, PillowWriter\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Function to get data based on tally and score\n", + "def get_tally_data(tally_name, score='flux'):\n", + " tally = sp.get_tally(name=tally_name)\n", + " return tally.get_values(scores=[score])\n", + "\n", + "# Retrieve the fission data (change tally name and score if needed)\n", + "data_type = 'flux' # Change this to 'flux', 'fission', or any other score available\n", + "plot_title = 'Neutron Flux' # Adjust title for visualization\n", + "tally_name = 'Flux' # Adjust the tally name if necessary\n", + "data = get_tally_data(tally_name=tally_name, score=data_type)\n", + "\n", + "# Define the mesh size and reshape the data\n", + "total_data_points = 300**3\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "data_reshaped = data.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Determine the data range for consistent color scaling\n", + "data_min = data_reshaped.min()\n", + "data_max = data_reshaped.max()\n", + "\n", + "# Initialize the figure and axis\n", + "fig, ax = plt.subplots()\n", + "\n", + "# Function to update the animation for each frame (z-slice)\n", + "def update(z_slice_index, plot_title):\n", + " ax.clear() # Clear previous frame\n", + " data_2d = data_reshaped[z_slice_index, :, :]\n", + " \n", + " # Generate the heatmap for the x-y slice\n", + " im = ax.imshow(data_2d, origin='lower', cmap='viridis', vmin=data_min, vmax=data_max)\n", + " ax.set_title(f'{plot_title} in the x-y plane at z index {z_slice_index}')\n", + " ax.set_xlabel('X-axis')\n", + " ax.set_ylabel('Y-axis')\n", + " return im,\n", + "\n", + "custom_title = \"Neutron Flux\"\n", + "\n", + "# Create the animation\n", + "num_frames = mesh_dimension_size # One frame per z-slice\n", + "ani = FuncAnimation(fig, update, frames=num_frames, fargs=(custom_title,), blit=True, interval=100)\n", + "\n", + "# Save the animation as a video\n", + "ani.save('neutron_flux_animation.gif', writer='ffmpeg', fps=10)\n", + "\n", + "# Optionally display the animation as HTML\n", + "from IPython.display import HTML\n", + "HTML(ani.to_jshtml())\n", + "\n", + "# Close the plot to clean up\n", + "plt.close()" + ] + }, + { + "cell_type": "markdown", + "id": "57fcf7cf-1be5-4cf5-9c3d-f4a2dfca99e1", + "metadata": {}, + "source": [ + "# Visualizing Errors" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "569e8033-1b7f-48a3-8de0-d3f61aeada83", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Access the flux tally\n", + "flux_tally = sp.get_tally(name='Flux')\n", + "\n", + "# Get the mean and standard deviation for the flux tally\n", + "flux_mean = flux_tally.get_values(scores=['flux'], value='mean')\n", + "flux_std_dev = flux_tally.get_values(scores=['flux'], value='std_dev')\n", + "\n", + "# Total number of data points\n", + "total_data_points = 300**3\n", + "\n", + "# Since the data is from a cubic mesh, we take the cube root to find the size of one dimension\n", + "mesh_dimension_size = np.round(np.cbrt(total_data_points)).astype(int)\n", + "\n", + "# Reshape the mean and std_dev data to the mesh dimensions\n", + "flux_mean_reshaped = flux_mean.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "flux_std_dev_reshaped = flux_std_dev.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Compute the relative error\n", + "relative_error = np.zeros_like(flux_mean_reshaped)\n", + "nonzero = flux_mean_reshaped > 0\n", + "relative_error[nonzero] = flux_std_dev_reshaped[nonzero] / flux_mean_reshaped[nonzero]\n", + "\n", + "# Plot the distribution of relative errors\n", + "plt.hist(relative_error[nonzero], bins=500)\n", + "plt.title(\"Distribution of Relative Errors (Flux)\")\n", + "plt.xlabel(\"Relative Error\")\n", + "plt.ylabel(\"Frequency\")\n", + "plt.show()\n", + "\n", + "# Access the fission tally\n", + "fission_tally = sp.get_tally(name='Fission')\n", + "\n", + "# Get the mean and standard deviation for the fission tally\n", + "fission_mean = fission_tally.get_values(scores=['fission'], value='mean')\n", + "fission_std_dev = fission_tally.get_values(scores=['fission'], value='std_dev')\n", + "\n", + "# Reshape the mean and std_dev data to the mesh dimensions\n", + "fission_mean_reshaped = fission_mean.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "fission_std_dev_reshaped = fission_std_dev.reshape((mesh_dimension_size, mesh_dimension_size, mesh_dimension_size))\n", + "\n", + "# Compute the relative error\n", + "relative_error_fission = np.zeros_like(fission_mean_reshaped)\n", + "nonzero_fission = fission_mean_reshaped > 0\n", + "relative_error_fission[nonzero_fission] = fission_std_dev_reshaped[nonzero_fission] / fission_mean_reshaped[nonzero_fission]\n", + "\n", + "# Plot the distribution of relative errors for fission\n", + "plt.hist(relative_error_fission[nonzero_fission], bins=500)\n", + "plt.title(\"Distribution of Relative Errors (Fission)\")\n", + "plt.xlabel(\"Relative Error\")\n", + "plt.ylabel(\"Frequency\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "f9fd4892-f958-4396-84d5-3111bde4fdf1", + "metadata": {}, + "source": [ + "# Plotting Heat Map of Errors" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "65c9224c-9ebc-4579-8147-af7a5c72f042", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from ipywidgets import interact, IntSlider\n", + "\n", + "# Load the statepoint file\n", + "sp = openmc.StatePoint('statepoint.15.h5')\n", + "\n", + "# Access the tallies\n", + "flux_tally = sp.get_tally(name='Flux')\n", + "fission_tally = sp.get_tally(name='Fission')\n", + "\n", + "# Extract mean and standard deviation for flux and fission\n", + "flux_data_mean = flux_tally.get_values(scores=['flux'], value='mean')\n", + "fission_data_mean = fission_tally.get_values(scores=['fission'], value='mean')\n", + "flux_data_std_dev = flux_tally.get_values(scores=['flux'], value='std_dev')\n", + "fission_data_std_dev = fission_tally.get_values(scores=['fission'], value='std_dev')\n", + "\n", + "# Define the mesh dimensions (assumed cubic with 300 points in each dimension)\n", + "data_points = 300\n", + "\n", + "# Reshape the 1D arrays into 3D volumes\n", + "flux_data_mean = flux_data_mean.reshape((data_points, data_points, data_points))\n", + "fission_data_mean = fission_data_mean.reshape((data_points, data_points, data_points))\n", + "flux_data_std_dev = flux_data_std_dev.reshape((data_points, data_points, data_points))\n", + "fission_data_std_dev = fission_data_std_dev.reshape((data_points, data_points, data_points))\n", + "\n", + "# Calculate the relative errors (std_dev / mean)\n", + "Relative_errors_flux = flux_data_std_dev / flux_data_mean\n", + "Relative_errors_fission = fission_data_std_dev / fission_data_mean\n", + "\n", + "# Convert to percentage\n", + "flux_error_percent = Relative_errors_flux * 100\n", + "fission_error_percent = Relative_errors_fission * 100\n", + "\n", + "vmin = 0\n", + "vmax = 100\n", + "\n", + "# Determine a common color scale range for both images\n", + "data_min = min(Relative_errors_flux.min(), Relative_errors_fission.min())\n", + "data_max = max(Relative_errors_flux.max(), Relative_errors_fission.max())\n", + "\n", + "def show_slice(z_slice_index):\n", + " # Create a figure with two subplots side by side\n", + " fig, axes = plt.subplots(1, 2, figsize=(12, 6))\n", + " \n", + " # Rotate the slices to orient them correctly for display\n", + " flux_slice = np.rot90(flux_error_percent[z_slice_index, :, :], 4)\n", + " fission_slice = np.rot90(fission_error_percent[z_slice_index, :, :], 4)\n", + " \n", + " # Plot flux relative error\n", + " im1 = axes[0].imshow(flux_slice, origin='lower', cmap='nipy_spectral', vmin=vmin, vmax=vmax)\n", + " axes[0].set_title(f'Flux Relative Error (Slice {z_slice_index})')\n", + " axes[0].set_xlabel('X-axis')\n", + " axes[0].set_ylabel('Y-axis')\n", + " plt.colorbar(im1, ax=axes[0])\n", + " \n", + " # Plot fission relative error\n", + " im2 = axes[1].imshow(fission_slice, origin='lower', cmap='nipy_spectral', vmin=vmin, vmax=vmax)\n", + " axes[1].set_title(f'Fission Relative Error (Slice {z_slice_index})')\n", + " axes[1].set_xlabel('X-axis')\n", + " axes[1].set_ylabel('Y-axis')\n", + " plt.colorbar(im2, ax=axes[1])\n", + " \n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "# Create an interactive slider to select the z-slice\n", + "interact(show_slice, z_slice_index=IntSlider(min=0, max=data_points-1, step=1, value=data_points//2));\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/README.md b/README.md index 4edf0e4..0d35457 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,67 @@ -# OSU-subcriticalassembly -This holds an openmc model of the subcritical assembly at OSU's radiation center +# OSU Subcritical Assembly + +This repository contains an OpenMC model of the subcritical assembly at Oregon State University's Radiation Center. The assembly is composed of graphite bricks, UO_2 fuel rods, a concrete reflector, and a plutonium–beryllium neutron source. + +## Running the Python Script +- Running the “RadCenterSubcriticalAssemblyOSU.ipynb” will require OpenMC to be installed correctly on your computer. When running the script, look at the “model.run” command at the very end of the script; currently, it is being multithreaded to 20 threads. Change this number to whatever you feel is necessary if you have enough threads in your CPU. To see how many CPU threads you have, look at Task Manager and click on the Performance Tab, there you will see on the bottom-left a list of technical specifications on your CPU, look to Logical Processors, that will be the number of threads you have. +- Take note as well to “settings.particles” and “mesh.dimension”, these values can be adjusted by the user as needed. Currently, the code will be very demanding on 20 cores due to the high values and will take a long time to compile. +- To run any analysis on the “Analysis.ipynb” script, make sure the outputs from the “RadCenterSubcriticalAssemblyOSU.ipynb” such as the “tallies.xml”, “summary.xml”, and the statepoint file are placed in the analysis folder. The script will not run at all if they are not present. +- If you change any parameters on the “RadCenterSubcriticalAssemblyOSU.ipynb” such as mesh dimension size, or the number of batches, some parameters on the “Analysis.ipynb” will need to be changed as a result to match. Running the script will tell you what you need to change. + +## Overview + +The subcritical assembly is designed to teach students about nuclear reactions and how to take measurements for reactors. The main components of the assembly are as follows: + +- **Graphite Bricks:** Serve as both structural support and neutron moderator. +- **UO_2 Fuel:** Annular fuel rods with natural uranium enrichment (0.072% U-235) and aluminum cladding. +- **Concrete Structural Material:** Enhances neutron economy by reflecting neutrons back into the assembly. +- **Plutonium–Beryllium Neutron Source:** Initiates the neutron population required for the assembly's operation. + +## Assembly Dimensions + +### UO₂ Fuel Rods +- **Type:** Annular rods with aluminum cladding. +- **Enrichment:** Natural uranium (0.072% U-235). +- **Dimensions:** + - **Fuel Length in z-direction:** 45 cm + - **Innermost Air Hole Radii:** 0.46 cm + - **Clad Inner Radii:** 0.61875 cm + - **Radii of the UO_2 Fuel:** 1.76125 cm + - **Radii of Outermost Clad:** 1.92 cm + - **Pitch of the Fuel Rods:**: 15 cm + +### Graphite Bricks +- **Dimensions:** + - **Length in x-direction:** 50 cm + - **Height in y-direction:** 15 cm + - **Width in z-direction:** 15 cm +- **Brick Arrangement:** + This assembly is made up of 300 graphite bricks with a total dimension size of 100 cm x 100 cm by 150 cm in the y-z-x direction respectively. Each brick with the dimensions listed above. +- **Fuel Arrangement:** + Each brick features a central hole for fuel placement that is the same radii as the outermost cladding of the fuel. Along the z-axis of the brick, three UO_2 fuel rods are inserted end to end. The combined length of these rods is 135 cm, with an additional 15 cm gap filled by air. + +### Concrete Reflector +- **Dimensions:** + - **Length in x-direction:** 15.24 cm + - **Height in y-direction:** 180 cm + - **Width in z-direction:** 150 cm +- **Placement:** + Positioned to the right of the assembly, the concrete wall functions as a reflector to bounce neutrons back into the reactor core. This concrete reflector's leftmost side is touching the rightmost side of the assembly. + +### Graphite Floor +- **Dimensions:** + - **Length in x-direction:** 150 cm + - **Height in y-direction:** 30 cm + - **Width in z-direction:** 150 cm + +### Neutron Source +- **Type:** Plutonium–Beryllium. +- **Location:** + Centrally located on the graphite floor at coordinates (75, -15, 75) with a raii of 2 cm within the reactor framework. + +## Documentation +- All documentation for how OpenMC works can be found on the OpenMC website linked below: +- https://docs.openmc.org/en/stable/ + +## acknowledgements +- Thank you to Dr. Steve Reese for providing necessary information to model this subcritical assembly such as the dimensions of the overall assembly, the x-y locations of each brick, the composition of the fuel and concrete, and the dimensions of the PuBe neutron source.