diff --git a/.devcontainer/base.Dockerfile b/.devcontainer/base.Dockerfile index 4a6b6326..50761245 100644 --- a/.devcontainer/base.Dockerfile +++ b/.devcontainer/base.Dockerfile @@ -94,6 +94,7 @@ RUN apt-get install libxinerama-dev -y # python packages from the neutronics workflow RUN pip install neutronics_material_maker[density] \ + stl_to_h5m \ remove_dagmc_tags \ openmc-plasma-source \ openmc-dagmc-wrapper \ @@ -247,7 +248,7 @@ RUN openmc_data_downloader -d nuclear_data -l ENDFB-8.0-NNDC TENDL-2019 -p neutr RUN pip install openmc_data && \ mkdir -p /nuclear_data && \ - download_nndc_chain -d nuclear_data -r b8.0 + download_endf_chain -d nuclear_data -r b8.0 # install WMP nuclear data RUN wget https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/WMP_Library_v1.1.tar.gz && \ @@ -255,4 +256,4 @@ RUN wget https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/WMP_Libr rm WMP_Library_v1.1.tar.gz ENV OPENMC_CROSS_SECTIONS=/nuclear_data/cross_sections.xml -ENV OPENMC_CHAIN_FILE=/nuclear_data/chain-nndc-b8.0.xml +ENV OPENMC_CHAIN_FILE=/nuclear_data/chain-endf-b8.0.xml diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml index 69f66d4d..d32ca73c 100644 --- a/.github/workflows/ci_tests.yml +++ b/.github/workflows/ci_tests.yml @@ -71,21 +71,18 @@ jobs: run: | pytest tests/test_task_12.py -v -# VR times out, reimplement onces cpp VR is ready -# - name: test task 13 -# run: | -# pytest tests/test_task_13.py -v + - name: test task 13 + run: | + pytest tests/test_task_13.py -v - name: test task 14 run: | pytest tests/test_task_14.py -v - - # TODO include when these have been speed up - # - name: test task 15 - # run: | - # pytest tests/test_task_15.py -v - - # - name: test task 15 - # run: | - # pytest tests/test_task_15.py -v - +# TODO add task 15 and 16 + - name: test task 17 + run: | + pytest tests/test_task_14.py -v + + - name: test task 18 + run: | + pytest tests/test_task_14.py -v diff --git a/tasks/task_09_CSG_dose_tallies/shut_down_dose_rate_example.py b/tasks/task_09_CSG_dose_tallies/shut_down_dose_rate_example.py index 28c59b51..0c9bcf4b 100644 --- a/tasks/task_09_CSG_dose_tallies/shut_down_dose_rate_example.py +++ b/tasks/task_09_CSG_dose_tallies/shut_down_dose_rate_example.py @@ -109,7 +109,6 @@ # final_step=False, operator_kwargs={ "normalization_mode": "source-rate", # needed as this is a fixed source simulation - "dilute_initial": 0, # need to avoid adding small amounts of fissle material "chain_file": chain_file }, ) diff --git a/tasks/task_13_variance_reduction/1_shielded_room_survival_biasing.ipynb b/tasks/task_13_variance_reduction/1_shielded_room_survival_biasing.ipynb new file mode 100644 index 00000000..77e9fda7 --- /dev/null +++ b/tasks/task_13_variance_reduction/1_shielded_room_survival_biasing.ipynb @@ -0,0 +1,421 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example simulates a shield room / bunker with corridor entrance a neutron source in the center of the room.\n", + "\n", + "Variance reduction is used to accelerate the simulation.\n", + "\n", + "The variance reduction method used for this simulation (survival_biasing) is not as effective as other variance reduction methods available in OpenMC but is the first task in the variance reduction section as it is the simplest to implement.\n", + "\n", + "Openmc documentation on survival-biasing\n", + "https://docs.openmc.org/en/stable/methods/neutron_physics.html#survival-biasing\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First we import openmc and other packages needed for the example" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib.colors import LogNorm # used for plotting log scale graphs" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create a couple of materials for the simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mat_air = openmc.Material(name=\"Air\")\n", + "mat_air.add_element(\"N\", 0.784431)\n", + "mat_air.add_element(\"O\", 0.210748)\n", + "mat_air.add_element(\"Ar\", 0.0046)\n", + "mat_air.set_density(\"g/cc\", 0.001205)\n", + "\n", + "mat_concrete = openmc.Material()\n", + "mat_concrete.add_element(\"H\",0.168759)\n", + "mat_concrete.add_element(\"C\",0.001416)\n", + "mat_concrete.add_element(\"O\",0.562524)\n", + "mat_concrete.add_element(\"Na\",0.011838)\n", + "mat_concrete.add_element(\"Mg\",0.0014)\n", + "mat_concrete.add_element(\"Al\",0.021354)\n", + "mat_concrete.add_element(\"Si\",0.204115)\n", + "mat_concrete.add_element(\"K\",0.005656)\n", + "mat_concrete.add_element(\"Ca\",0.018674)\n", + "mat_concrete.add_element(\"Fe\",0.00426)\n", + "mat_concrete.set_density(\"g/cm3\", 2.3)\n", + "\n", + "materials = openmc.Materials([mat_air, mat_concrete])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we define and plot the geometry. This geometry is define by parameters for every width and height. The parameters input into the geometry in a stacked manner so they can easily be adjusted to change the geometry without creating overlapping cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "width_a = 100\n", + "width_b = 200\n", + "width_c = 500\n", + "width_d = 250\n", + "width_e = 200\n", + "width_f = 200\n", + "width_g = 100\n", + "\n", + "depth_a = 100\n", + "depth_b = 200\n", + "depth_c = 700\n", + "depth_d = 600\n", + "depth_e = 200\n", + "depth_f = 100\n", + "\n", + "height_j = 100\n", + "height_k = 500\n", + "height_l = 100\n", + "\n", + "xplane_0 = openmc.XPlane(x0=0, boundary_type=\"vacuum\")\n", + "xplane_1 = openmc.XPlane(x0=xplane_0.x0 + width_a)\n", + "xplane_2 = openmc.XPlane(x0=xplane_1.x0 + width_b)\n", + "xplane_3 = openmc.XPlane(x0=xplane_2.x0 + width_c)\n", + "xplane_4 = openmc.XPlane(x0=xplane_3.x0 + width_d)\n", + "xplane_5 = openmc.XPlane(x0=xplane_4.x0 + width_e)\n", + "xplane_6 = openmc.XPlane(x0=xplane_5.x0 + width_f)\n", + "xplane_7 = openmc.XPlane(x0=xplane_6.x0 + width_g, boundary_type=\"vacuum\")\n", + "\n", + "yplane_0 = openmc.YPlane(y0=0, boundary_type=\"vacuum\")\n", + "yplane_1 = openmc.YPlane(y0=yplane_0.y0 + depth_a)\n", + "yplane_2 = openmc.YPlane(y0=yplane_1.y0 + depth_b)\n", + "yplane_3 = openmc.YPlane(y0=yplane_2.y0 + depth_c)\n", + "yplane_4 = openmc.YPlane(y0=yplane_3.y0 + depth_d)\n", + "yplane_5 = openmc.YPlane(y0=yplane_4.y0 + depth_e)\n", + "yplane_6 = openmc.YPlane(y0=yplane_5.y0 + depth_f, boundary_type=\"vacuum\")\n", + "\n", + "zplane_1 = openmc.ZPlane(z0=0, boundary_type=\"vacuum\")\n", + "zplane_2 = openmc.ZPlane(z0=zplane_1.z0 + height_j)\n", + "zplane_3 = openmc.ZPlane(z0=zplane_2.z0 + height_k)\n", + "zplane_4 = openmc.ZPlane(z0=zplane_3.z0 + height_l, boundary_type=\"vacuum\")\n", + "\n", + "outside_left_region = +xplane_0 & -xplane_1 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_4\n", + "wall_left_region = +xplane_1 & -xplane_2 & +yplane_2 & -yplane_4 & +zplane_2 & -zplane_3\n", + "wall_right_region = +xplane_5 & -xplane_6 & +yplane_2 & -yplane_5 & +zplane_2 & -zplane_3\n", + "wall_top_region = +xplane_1 & -xplane_4 & +yplane_4 & -yplane_5 & +zplane_2 & -zplane_3\n", + "outside_top_region = +xplane_0 & -xplane_7 & +yplane_5 & -yplane_6 & +zplane_1 & -zplane_4\n", + "wall_bottom_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_2 & +zplane_2 & -zplane_3\n", + "outside_bottom_region = +xplane_0 & -xplane_7 & +yplane_0 & -yplane_1 & +zplane_1 & -zplane_4\n", + "wall_middle_region = +xplane_3 & -xplane_4 & +yplane_3 & -yplane_4 & +zplane_2 & -zplane_3\n", + "outside_right_region = +xplane_6 & -xplane_7 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_4\n", + "\n", + "room_region = +xplane_2 & -xplane_3 & +yplane_2 & -yplane_4 & +zplane_2 & -zplane_3\n", + "gap_region = +xplane_3 & -xplane_4 & +yplane_2 & -yplane_3 & +zplane_2 & -zplane_3\n", + "corridor_region = +xplane_4 & -xplane_5 & +yplane_2 & -yplane_5 & +zplane_2 & -zplane_3\n", + "\n", + "roof_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_2\n", + "floor_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_5 & +zplane_3 & -zplane_4\n", + "\n", + "outside_left_cell = openmc.Cell(region=outside_left_region, fill=mat_air)\n", + "outside_right_cell = openmc.Cell(region=outside_right_region, fill=mat_air)\n", + "outside_top_cell = openmc.Cell(region=outside_top_region, fill=mat_air)\n", + "outside_bottom_cell = openmc.Cell(region=outside_bottom_region, fill=mat_air)\n", + "wall_left_cell = openmc.Cell(region=wall_left_region, fill=mat_concrete)\n", + "wall_right_cell = openmc.Cell(region=wall_right_region, fill=mat_concrete)\n", + "wall_top_cell = openmc.Cell(region=wall_top_region, fill=mat_concrete)\n", + "wall_bottom_cell = openmc.Cell(region=wall_bottom_region, fill=mat_concrete)\n", + "wall_middle_cell = openmc.Cell(region=wall_middle_region, fill=mat_concrete)\n", + "room_cell = openmc.Cell(region=room_region, fill=mat_air)\n", + "gap_cell = openmc.Cell(region=gap_region, fill=mat_air)\n", + "corridor_cell = openmc.Cell(region=corridor_region, fill=mat_air)\n", + "\n", + "roof_cell = openmc.Cell(region=roof_region, fill=mat_concrete)\n", + "floor_cell = openmc.Cell(region=floor_region, fill=mat_concrete)\n", + "\n", + "geometry = openmc.Geometry(\n", + " [\n", + " outside_bottom_cell,\n", + " outside_top_cell,\n", + " outside_left_cell,\n", + " outside_right_cell,\n", + " wall_left_cell,\n", + " wall_right_cell,\n", + " wall_top_cell,\n", + " wall_bottom_cell,\n", + " wall_middle_cell,\n", + " room_cell,\n", + " gap_cell,\n", + " corridor_cell,\n", + " roof_cell,\n", + " floor_cell,\n", + " ]\n", + ")\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we plot the geometry and color by materials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "geometry.root_universe.plot(basis='xy', color_by='material') \n", + "plt.savefig('geometry_top_down_view.png', bbox_inches=\"tight\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we create a point source, this also uses the same geometry parameters to place in the center of the room regardless of the values of the parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# location of the point source\n", + "source_x = width_a + width_b + width_c * 0.5\n", + "source_y = depth_a + depth_b + depth_c * 0.75\n", + "source_z = height_j + height_k * 0.5\n", + "space = openmc.stats.Point((source_x, source_y, source_z))\n", + "\n", + "angle = openmc.stats.Isotropic()\n", + "\n", + "# all (100%) of source particles are 2.5MeV energy\n", + "energy = openmc.stats.Discrete([2.5e6], [1.0])\n", + "\n", + "source = openmc.Source(space=space, angle=angle, energy=energy)\n", + "source.particle = \"neutron\"\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create a mesh that encompasses the entire geoemtry and scores neutron flux" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mesh = openmc.RegularMesh().from_domain(geometry)\n", + "mesh.dimension = (100, 100, 1)\n", + "\n", + "mesh_filter = openmc.MeshFilter(mesh)\n", + "particle_filter = openmc.ParticleFilter('neutron')\n", + "\n", + "flux_tally = openmc.Tally(name=\"flux tally\")\n", + "flux_tally.filters = [mesh_filter, particle_filter]\n", + "flux_tally.scores = [\"flux\"]\n", + "\n", + "tallies = openmc.Tallies([flux_tally])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Creates the simulation settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "settings = openmc.Settings()\n", + "\n", + "settings.run_mode = \"fixed source\"\n", + "settings.source = source\n", + "settings.particles = 2000\n", + "settings.batches = 5" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Creates the model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "model = openmc.Model(geometry, materials, settings, tallies)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We are going to run this model twice. Once with survival_biasing and once without. This next code bock makes a function so that we can easily run the simulation twice with these two different settings." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def run_and_plot(model, image_filename):\n", + "\n", + " sp_filename = model.run()\n", + "\n", + " with openmc.StatePoint(sp_filename) as sp:\n", + " flux_tally = sp.get_tally(name=\"flux tally\")\n", + "\n", + " mesh_extent = mesh.bounding_box.extent['xy']\n", + "\n", + " # create a plot of the mean flux values\n", + " flux_mean = flux_tally.get_reshaped_data(value='mean', expand_dims=True).squeeze()\n", + " plt.subplot(1, 2, 1)\n", + " plt.imshow(\n", + " flux_mean.T,\n", + " origin=\"lower\",\n", + " extent=mesh_extent,\n", + " norm=LogNorm(),\n", + " )\n", + " plt.title(\"Flux Mean\")\n", + "\n", + "\n", + " plt.subplot(1, 2, 2)\n", + " # create a plot of the flux relative error\n", + " flux_std_dev = flux_tally.get_reshaped_data(value='std_dev', expand_dims=True).squeeze()\n", + " plt.imshow(\n", + " flux_std_dev.T,\n", + " origin=\"lower\",\n", + " extent=mesh_extent,\n", + " norm=LogNorm(),\n", + " )\n", + " plt.title(\"Flux Std. Dev.\")\n", + "\n", + " plt.savefig(image_filename)\n", + " return sp" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we run the model without survival_biasing and save the meesh tally image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "run_and_plot(model, \"no_survival_biasing.png\")\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we run the model again with survival_biasing enabled. This needs to me set on the model settings. The values chosen can be fine tunned to the model but for this example we have chosen some values after a little exploration." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "model.settings.survival_biasing = True\n", + "model.settings.cutoff = {\n", + " \"weight\": 0.3, # value needs to be between 0 and 1\n", + " \"weight_avg\": 0.9, # value needs to be between 0 and 1\n", + "}\n", + "run_and_plot(model, \"yes_survival_biasing.png\")\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Comparing the two simulations we notice that the survival_biasing simulation resulted in slightly improved flux tally and a slightly improved standard deviation.\n", + "\n", + "This model is used again in the next example that makes use of weight windows and yields more impressive acceleration. " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "openmc_dev", + "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.10.10" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/tasks/task_13_variance_reduction/2_shielded_room_single_ww.ipynb b/tasks/task_13_variance_reduction/2_shielded_room_single_ww.ipynb new file mode 100644 index 00000000..839ea997 --- /dev/null +++ b/tasks/task_13_variance_reduction/2_shielded_room_single_ww.ipynb @@ -0,0 +1,551 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "72e0b25d-e541-4e8d-8805-b984374ee53d", + "metadata": {}, + "source": [ + "# Variance Reduction - Weight Windows\n", + "\n", + "## Creating and utilizing a wight window to accelerate deep shielding simulations\n", + "\n", + "This example simulates a shield room / bunker with corridor entrance a neutron source in the center of the room with. This example implements a single step of the Magic method of weight window generation. \n", + "\n", + "In this tutorial we shall focus on generating a weight window to accelerate the simulation of particles through a shield.\n", + "\n", + "Weight Windows are found using the MAGIC method and used to accelerate the simulation.\n", + "\n", + "The variance reduction method used for this simulation is well documented in the OpenMC documentation\n", + "https://docs.openmc.org/en/stable/methods/neutron_physics.html\n", + "\n", + "The MAGIC method is well described in the original publication\n", + "https://scientific-publications.ukaea.uk/wp-content/uploads/Published/INTERN1.pdf\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "f00eddb1-1e3a-4973-ba0c-d8feeb2a6704", + "metadata": {}, + "source": [ + "First we import ```openmc``` including ```openmc.lib``` and other packages needed for the example" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3448d1db-4328-42e3-8960-50d53896f541", + "metadata": {}, + "outputs": [], + "source": [ + "import time # used to time the simulation\n", + "\n", + "import openmc\n", + "import openmc.lib # this example makes use of openmc lib to run the simulations\n", + "\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib.colors import LogNorm # used for plotting log scale graphs" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "8179e89a-304b-4684-91a9-8a908b75e8cd", + "metadata": {}, + "source": [ + "We create a couple of materials for the simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47cb21c7-0afa-446c-9a6b-76c8dfb2f93e", + "metadata": {}, + "outputs": [], + "source": [ + "mat_air = openmc.Material(name=\"Air\")\n", + "mat_air.add_element(\"N\", 0.784431)\n", + "mat_air.add_element(\"O\", 0.210748)\n", + "mat_air.add_element(\"Ar\", 0.0046)\n", + "mat_air.set_density(\"g/cc\", 0.001205)\n", + "\n", + "mat_concrete = openmc.Material()\n", + "mat_concrete.add_element(\"H\",0.168759)\n", + "mat_concrete.add_element(\"C\",0.001416)\n", + "mat_concrete.add_element(\"O\",0.562524)\n", + "mat_concrete.add_element(\"Na\",0.011838)\n", + "mat_concrete.add_element(\"Mg\",0.0014)\n", + "mat_concrete.add_element(\"Al\",0.021354)\n", + "mat_concrete.add_element(\"Si\",0.204115)\n", + "mat_concrete.add_element(\"K\",0.005656)\n", + "mat_concrete.add_element(\"Ca\",0.018674)\n", + "mat_concrete.add_element(\"Fe\",0.00426)\n", + "mat_concrete.set_density(\"g/cm3\", 2.3)\n", + "\n", + "materials = openmc.Materials([mat_air, mat_concrete])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "424b234c-48a2-4d1b-b374-f69430e464db", + "metadata": {}, + "source": [ + "Now we define and plot the geometry. This geometry is define by parameters for every width and height. The parameters input into the geometry in a stacked manner so they can easily be adjusted to change the geometry without creating overlapping cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a72dc5e3-cb8e-41c2-9cfb-f8f7489c7f44", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "width_a = 100\n", + "width_b = 100\n", + "width_c = 500\n", + "width_d = 100\n", + "width_e = 100\n", + "width_f = 100\n", + "width_g = 100\n", + "\n", + "depth_a = 100\n", + "depth_b = 100\n", + "depth_c = 700\n", + "depth_d = 600\n", + "depth_e = 100\n", + "depth_f = 100\n", + "\n", + "height_j = 100\n", + "height_k = 500\n", + "height_l = 100\n", + "\n", + "xplane_0 = openmc.XPlane(x0=0, boundary_type=\"vacuum\")\n", + "xplane_1 = openmc.XPlane(x0=xplane_0.x0 + width_a)\n", + "xplane_2 = openmc.XPlane(x0=xplane_1.x0 + width_b)\n", + "xplane_3 = openmc.XPlane(x0=xplane_2.x0 + width_c)\n", + "xplane_4 = openmc.XPlane(x0=xplane_3.x0 + width_d)\n", + "xplane_5 = openmc.XPlane(x0=xplane_4.x0 + width_e)\n", + "xplane_6 = openmc.XPlane(x0=xplane_5.x0 + width_f)\n", + "xplane_7 = openmc.XPlane(x0=xplane_6.x0 + width_g, boundary_type=\"vacuum\")\n", + "\n", + "yplane_0 = openmc.YPlane(y0=0, boundary_type=\"vacuum\")\n", + "yplane_1 = openmc.YPlane(y0=yplane_0.y0 + depth_a)\n", + "yplane_2 = openmc.YPlane(y0=yplane_1.y0 + depth_b)\n", + "yplane_3 = openmc.YPlane(y0=yplane_2.y0 + depth_c)\n", + "yplane_4 = openmc.YPlane(y0=yplane_3.y0 + depth_d)\n", + "yplane_5 = openmc.YPlane(y0=yplane_4.y0 + depth_e)\n", + "yplane_6 = openmc.YPlane(y0=yplane_5.y0 + depth_f, boundary_type=\"vacuum\")\n", + "\n", + "zplane_1 = openmc.ZPlane(z0=0, boundary_type=\"vacuum\")\n", + "zplane_2 = openmc.ZPlane(z0=zplane_1.z0 + height_j)\n", + "zplane_3 = openmc.ZPlane(z0=zplane_2.z0 + height_k)\n", + "zplane_4 = openmc.ZPlane(z0=zplane_3.z0 + height_l, boundary_type=\"vacuum\")\n", + "\n", + "outside_left_region = +xplane_0 & -xplane_1 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_4\n", + "wall_left_region = +xplane_1 & -xplane_2 & +yplane_2 & -yplane_4 & +zplane_2 & -zplane_3\n", + "wall_right_region = +xplane_5 & -xplane_6 & +yplane_2 & -yplane_5 & +zplane_2 & -zplane_3\n", + "wall_top_region = +xplane_1 & -xplane_4 & +yplane_4 & -yplane_5 & +zplane_2 & -zplane_3\n", + "outside_top_region = +xplane_0 & -xplane_7 & +yplane_5 & -yplane_6 & +zplane_1 & -zplane_4\n", + "wall_bottom_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_2 & +zplane_2 & -zplane_3\n", + "outside_bottom_region = +xplane_0 & -xplane_7 & +yplane_0 & -yplane_1 & +zplane_1 & -zplane_4\n", + "wall_middle_region = +xplane_3 & -xplane_4 & +yplane_3 & -yplane_4 & +zplane_2 & -zplane_3\n", + "outside_right_region = +xplane_6 & -xplane_7 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_4\n", + "\n", + "room_region = +xplane_2 & -xplane_3 & +yplane_2 & -yplane_4 & +zplane_2 & -zplane_3\n", + "gap_region = +xplane_3 & -xplane_4 & +yplane_2 & -yplane_3 & +zplane_2 & -zplane_3\n", + "corridor_region = +xplane_4 & -xplane_5 & +yplane_2 & -yplane_5 & +zplane_2 & -zplane_3\n", + "\n", + "roof_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_2\n", + "floor_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_5 & +zplane_3 & -zplane_4\n", + "\n", + "outside_left_cell = openmc.Cell(region=outside_left_region, fill=mat_air)\n", + "outside_right_cell = openmc.Cell(region=outside_right_region, fill=mat_air)\n", + "outside_top_cell = openmc.Cell(region=outside_top_region, fill=mat_air)\n", + "outside_bottom_cell = openmc.Cell(region=outside_bottom_region, fill=mat_air)\n", + "wall_left_cell = openmc.Cell(region=wall_left_region, fill=mat_concrete)\n", + "wall_right_cell = openmc.Cell(region=wall_right_region, fill=mat_concrete)\n", + "wall_top_cell = openmc.Cell(region=wall_top_region, fill=mat_concrete)\n", + "wall_bottom_cell = openmc.Cell(region=wall_bottom_region, fill=mat_concrete)\n", + "wall_middle_cell = openmc.Cell(region=wall_middle_region, fill=mat_concrete)\n", + "room_cell = openmc.Cell(region=room_region, fill=mat_air)\n", + "gap_cell = openmc.Cell(region=gap_region, fill=mat_air)\n", + "corridor_cell = openmc.Cell(region=corridor_region, fill=mat_air)\n", + "\n", + "roof_cell = openmc.Cell(region=roof_region, fill=mat_concrete)\n", + "floor_cell = openmc.Cell(region=floor_region, fill=mat_concrete)\n", + "\n", + "geometry = openmc.Geometry(\n", + " [\n", + " outside_bottom_cell,\n", + " outside_top_cell,\n", + " outside_left_cell,\n", + " outside_right_cell,\n", + " wall_left_cell,\n", + " wall_right_cell,\n", + " wall_top_cell,\n", + " wall_bottom_cell,\n", + " wall_middle_cell,\n", + " room_cell,\n", + " gap_cell,\n", + " corridor_cell,\n", + " roof_cell,\n", + " floor_cell,\n", + " ]\n", + ")\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "51713fe1", + "metadata": {}, + "source": [ + "Now we plot the geometry and color by materials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "94c41f17", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "geometry.root_universe.plot(basis='xy', color_by='material') \n", + "plt.savefig('geometry_top_down_view.png', bbox_inches=\"tight\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "8d6bd790", + "metadata": {}, + "source": [ + "Next we create a point source, this also uses the same geometry parameters to place in the center of the room regardless of the values of the parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd4986b0", + "metadata": {}, + "outputs": [], + "source": [ + "# location of the point source\n", + "source_x = width_a + width_b + width_c * 0.5\n", + "source_y = depth_a + depth_b + depth_c * 0.75\n", + "source_z = height_j + height_k * 0.5\n", + "space = openmc.stats.Point((source_x, source_y, source_z))\n", + "\n", + "angle = openmc.stats.Isotropic()\n", + "\n", + "# all (100%) of source particles are 2.5MeV energy\n", + "energy = openmc.stats.Discrete([2.5e6], [1.0])\n", + "\n", + "source = openmc.Source(space=space, angle=angle, energy=energy)\n", + "source.particle = \"neutron\"\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "1456ac2f", + "metadata": {}, + "source": [ + "Next we create a mesh that encompasses the entire geometry and scores neutron flux" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d7bbee6", + "metadata": {}, + "outputs": [], + "source": [ + "mesh = openmc.RegularMesh().from_domain(geometry)\n", + "mesh.dimension = (500, 500, 1)\n", + "\n", + "mesh_filter = openmc.MeshFilter(mesh)\n", + "particle_filter = openmc.ParticleFilter('neutron')\n", + "\n", + "flux_tally = openmc.Tally(name=\"flux tally\")\n", + "flux_tally.filters = [mesh_filter, particle_filter]\n", + "flux_tally.scores = [\"flux\"]\n", + "flux_tally.id = 42 # we set the ID because we need to access this later\n", + "\n", + "tallies = openmc.Tallies([flux_tally])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "2bc34eac", + "metadata": {}, + "source": [ + "Creates the simulation settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea32fb73", + "metadata": {}, + "outputs": [], + "source": [ + "settings = openmc.Settings()\n", + "\n", + "settings.run_mode = \"fixed source\"\n", + "settings.source = source\n", + "settings.particles = 80000\n", + "settings.batches = 5\n", + "# no need to write the tallies.out file which saves space and time when large meshes are used\n", + "settings.output = {'tallies': False}" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "9bc4e52b", + "metadata": {}, + "source": [ + "Creates and export the model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d8042a5", + "metadata": {}, + "outputs": [], + "source": [ + "model = openmc.Model(geometry, materials, settings, tallies)\n", + "model.export_to_xml() # this is necessary as openmc.lib loads up the model.xml file" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "488391db", + "metadata": {}, + "source": [ + "Now we make use of openmc.lib to control the simulation. Documentation on openmc.lib is here\n", + "https://docs.openmc.org/en/stable/pythonapi/capi.html\n", + "\n", + "We also time the simulation so that we can perform the same simulation again with weight windows and try to fine tune the particle so that both simulations take the same time and we are making a fair comparision." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9159a336", + "metadata": {}, + "outputs": [], + "source": [ + "# this helps time the simulation\n", + "t0 = time.time()\n", + "\n", + "# first we initialize openmc lib, this reads in the model.xml and material cross sections \n", + "openmc.lib.init()\n", + "\n", + "# This runs openmc with the settings provided earlier\n", + "openmc.lib.run()\n", + "\n", + "# End the connection to openmc lib and write out the statepoint file\n", + "openmc.lib.finalize()\n", + "\n", + "t1 = time.time()\n", + "\n", + "total = t1-t0\n", + "\n", + "print(f'total time without weight windows {total}s')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "79ec2340", + "metadata": {}, + "source": [ + "Now we want to plot the results of the simulation. We want to do this twice to compare the results so I've written this up as a function that we can call." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6f60d98", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_mesh_tally(statepoint_filename, image_filename):\n", + "\n", + " with openmc.StatePoint(statepoint_filename) as sp:\n", + " flux_tally = sp.get_tally(name=\"flux tally\")\n", + "\n", + " mesh_extent = mesh.bounding_box.extent['xy']\n", + "\n", + " # get a slice of mean values on the xy basis mid z axis\n", + " flux_mean = flux_tally.get_reshaped_data(value='mean', expand_dims=True).squeeze()\n", + " plt.subplot(1, 2, 1)\n", + " # create a plot of the mean flux values\n", + " plt.imshow(\n", + " flux_mean.T,\n", + " origin=\"lower\",\n", + " extent=mesh_extent,\n", + " norm=LogNorm(),\n", + " )\n", + " plt.title(\"Flux Mean\")\n", + "\n", + "\n", + " plt.subplot(1, 2, 2)\n", + " # get a slice of std dev values on the xy basis mid z axis\n", + " flux_std_dev = flux_tally.get_reshaped_data(value='std_dev', expand_dims=True).squeeze()\n", + " # create a plot of the flux relative error\n", + " plt.imshow(\n", + " flux_std_dev.T,\n", + " origin=\"lower\",\n", + " extent=mesh_extent,\n", + " norm=LogNorm(),\n", + " )\n", + " plt.title(\"Flux Std. Dev.\")\n", + "\n", + " plt.savefig(image_filename)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b15e186f", + "metadata": {}, + "source": [ + "This next section calls the plotting function and saves an image of the mesh tally as no_ww.png" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8fb8f31e", + "metadata": {}, + "outputs": [], + "source": [ + "plot_mesh_tally(f'statepoint.{settings.batches}.h5', 'no_ww.png')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "2abd8995", + "metadata": {}, + "source": [ + "Now we want to run the simulation twice. Once to generate weight windows and once to make use of the weight windows. We are going to adjust the number of particles to try and spend the same total amount of computer time as the last simulation and make this a fair test. This might require some adaption of the particles as the simulation is set up for my laptop." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d519ac33", + "metadata": {}, + "outputs": [], + "source": [ + "# this helps time the simulation\n", + "t0 = time.time()\n", + "\n", + "# first we initialize openmc lib, this reads in the model.xml and material cross sections \n", + "openmc.lib.init()\n", + "\n", + "# Now we find the mesh tally that will be used for creating the weight windows \n", + "# This ID of 42 matches the tally ID we set earlier\n", + "# At this point the tally is empty.\n", + "tally = openmc.lib.tallies[42]\n", + "\n", + "# We create a weight window object from the empty tally. \n", + "wws = openmc.lib.WeightWindows.from_tally(tally, particle='neutron')\n", + "\n", + "# We are running a fraction of particles of the original simulation as we want\n", + "# to take less simulation time for the first run and leave some simulation time\n", + "# for the run with weight windows.\n", + "openmc.lib.settings.particles = 12000\n", + "\n", + "# This runs openmc with the settings provided earlier\n", + "openmc.lib.run()\n", + "\n", + "# The tally now contains meaningful information.\n", + "# So we can update the weight windows from the tally.\n", + "# This uses the MAGIC method\n", + "wws.update_magic(tally)\n", + "\n", + "# Now we enable weight window usage, previously this was set to False by default\n", + "openmc.lib.settings.weight_windows_on = True\n", + "\n", + "# As the particles are now splitting the number of particles run per second\n", + "# will be significantly lower. Therefore we will reduce the number of particles\n", + "# so that the simulation runs in reasonable amount of time\n", + "openmc.lib.settings.particles = 380\n", + "\n", + "# Now we run the simulation again, this time with the weight windows in use\n", + "openmc.lib.run()\n", + "\n", + "# End the connection to openmc lib and write out the statepoint file\n", + "openmc.lib.finalize()\n", + "\n", + "t1 = time.time()\n", + "\n", + "total = t1-t0\n", + "\n", + "print(f'total time without weight windows {total}s')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fbf35ffc", + "metadata": {}, + "outputs": [], + "source": [ + "plot_mesh_tally(f'statepoint.{settings.batches}.h5', 'ww.png')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b086b6ec", + "metadata": {}, + "source": [ + "On my laptop both simulations took 20 seconds to complete but the resulting flux map from the simulation with weight windows shows neutrons got further through the geometry.\n", + "\n", + "The weight windows allow the computer to spend a larger proportion of time simulating valuable neutrons that have reached regions of low neutron flux." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c56069a4", + "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.10.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tasks/task_13_variance_reduction/3_sphere_iterative_per_run_ww.ipynb b/tasks/task_13_variance_reduction/3_sphere_iterative_per_run_ww.ipynb new file mode 100644 index 00000000..1706741f --- /dev/null +++ b/tasks/task_13_variance_reduction/3_sphere_iterative_per_run_ww.ipynb @@ -0,0 +1,450 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "72e0b25d-e541-4e8d-8805-b984374ee53d", + "metadata": {}, + "source": [ + "# Variance Reduction - Weight Windows\n", + "\n", + "## Iteratively creating and utilizing a wight window to accelerate deep shielding simulations\n", + "\n", + "This example simulates a sphere of material with a neutron source in the center.This example implements the MAGIC method of weight window generation on each simulation run.\n", + "\n", + "In this tutorial we shall focus on generating a weight window to accelerate the simulation of particles through a shield and improving the weight window with each iteration.\n", + "\n", + "Weight Windows are found using the MAGIC method and used to accelerate the simulation.\n", + "\n", + "The variance reduction method used for this simulation is well documented in the OpenMC documentation\n", + "https://docs.openmc.org/en/stable/methods/neutron_physics.html\n", + "\n", + "The MAGIC method is well described in the original publication\n", + "https://scientific-publications.ukaea.uk/wp-content/uploads/Published/INTERN1.pdf\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "f00eddb1-1e3a-4973-ba0c-d8feeb2a6704", + "metadata": {}, + "source": [ + "First we import ```openmc``` including ```openmc.lib``` and other packages needed for the example" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3448d1db-4328-42e3-8960-50d53896f541", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "import openmc.lib # this example makes use of openmc lib to run the simulations\n", + "\n", + "import numpy as np\n", + "\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib.colors import LogNorm # used for plotting log scale graphs" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "8179e89a-304b-4684-91a9-8a908b75e8cd", + "metadata": {}, + "source": [ + "We create a couple of materials for the simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47cb21c7-0afa-446c-9a6b-76c8dfb2f93e", + "metadata": {}, + "outputs": [], + "source": [ + "mat_water = openmc.Material()\n", + "mat_water.add_element(\"H\", 1)\n", + "mat_water.add_element(\"O\", 2)\n", + "mat_water.set_density(\"g/cm3\", 1.0)\n", + "\n", + "my_materials = openmc.Materials([mat_water])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "424b234c-48a2-4d1b-b374-f69430e464db", + "metadata": {}, + "source": [ + "Now we define and plot the spherical geometry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a72dc5e3-cb8e-41c2-9cfb-f8f7489c7f44", + "metadata": {}, + "outputs": [], + "source": [ + "# outer surface at 500 cm\n", + "outer_surface = openmc.model.RectangularParallelepiped(-300, 300, -300, 300, -300, 300, boundary_type=\"vacuum\")\n", + "\n", + "# A single region below the surface\n", + "region_1 = -outer_surface\n", + "\n", + "# A single cell full of water\n", + "cell_1 = openmc.Cell(region=region_1)\n", + "cell_1.fill = mat_water\n", + "\n", + "my_geometry = openmc.Geometry([cell_1])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "51713fe1", + "metadata": {}, + "source": [ + "Now we plot the geometry and color by materials." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "94c41f17", + "metadata": {}, + "outputs": [], + "source": [ + "my_geometry.root_universe.plot(basis='xy', color_by='material') \n", + "plt.savefig('geometry_top_down_view.png', bbox_inches=\"tight\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "8d6bd790", + "metadata": {}, + "source": [ + "Next we create a point source, this also uses the same geometry parameters to place in the center of the room regardless of the values of the parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd4986b0", + "metadata": {}, + "outputs": [], + "source": [ + "# location of the point source\n", + "space = openmc.stats.Point((0, 0, 0))\n", + "angle = openmc.stats.Isotropic()\n", + "\n", + "# all (100%) of source particles are 14MeV energy\n", + "energy = openmc.stats.Discrete([14e6], [1.0])\n", + "\n", + "source = openmc.Source(space=space, angle=angle, energy=energy)\n", + "source.particle = \"neutron\"" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "1456ac2f", + "metadata": {}, + "source": [ + "Next we create a mesh that encompasses the entire geometry and scores neutron flux" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d7bbee6", + "metadata": {}, + "outputs": [], + "source": [ + "mesh = openmc.RegularMesh().from_domain(domain=my_geometry)\n", + "print(mesh)\n", + "# mesh.r_grid = np.linspace(0, outer_surface.r, 100)\n", + "\n", + "mesh_filter = openmc.MeshFilter(mesh)\n", + "\n", + "flux_tally = openmc.Tally(name=\"flux tally\")\n", + "flux_tally.filters = [mesh_filter]\n", + "flux_tally.scores = [\"flux\"]\n", + "flux_tally.id = 55 # we set the ID number here as we need to access it during the openmc lib run\n", + "\n", + "# adds the mesh tally to the model\n", + "my_tallies = openmc.Tallies()\n", + "my_tallies.append(flux_tally)\n", + "\n", + "tallies = openmc.Tallies([flux_tally])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "2bc34eac", + "metadata": {}, + "source": [ + "Creates the simulation settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea32fb73", + "metadata": {}, + "outputs": [], + "source": [ + "my_settings = openmc.Settings()\n", + "my_settings.run_mode = \"fixed source\"\n", + "my_settings.source = source\n", + "my_settings.particles = 120\n", + "my_settings.batches = 10\n", + "\n", + "# no need to write the tallies.out file which saves space and time when large meshes are used\n", + "my_settings.output = {'tallies': False}" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "9bc4e52b", + "metadata": {}, + "source": [ + "Creates and export the model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d8042a5", + "metadata": {}, + "outputs": [], + "source": [ + "model = openmc.Model(my_geometry, my_materials, my_settings, my_tallies)\n", + "\n", + "# deletes old input and output files\n", + "!rm *.xml\n", + "!rm *.h5\n", + "\n", + "model.export_to_xml() # this is necessary as openmc.lib loads up the model.xml file" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "bf5eeb8b", + "metadata": {}, + "source": [ + "Now we want to plot the results of the simulation. We want to do this twice to compare the results so I've written this up as a function that we can call." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc9ec8cb", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_mesh_tally_and_weight_window(statepoint_filename, weight_window_filename, image_filename):\n", + " \n", + " with openmc.StatePoint(statepoint_filename) as sp:\n", + " flux_tally = sp.get_tally(name=\"flux tally\")\n", + "\n", + " tally_mesh = flux_tally.find_filter(openmc.MeshFilter).mesh\n", + " tally_mesh_extent = tally_mesh.bounding_box.extent['xy']\n", + "\n", + " # get a slice of mean values on the xy basis mid z axis\n", + " flux_mean = flux_tally.get_reshaped_data(value='mean', expand_dims=True).squeeze()[:,:,int(mesh.dimension[2]/2)]\n", + " plt.subplot(1, 3, 1)\n", + " # create a plot of the mean flux values\n", + " plt.imshow(\n", + " flux_mean.T,\n", + " origin=\"lower\",\n", + " extent=tally_mesh_extent,\n", + " norm=LogNorm(),\n", + " )\n", + " plt.title(\"Flux Mean\")\n", + "\n", + "\n", + " plt.subplot(1, 3, 2)\n", + " # get a slice of std dev values on the xy basis mid z axis\n", + " flux_std_dev = flux_tally.get_reshaped_data(value='std_dev', expand_dims=True).squeeze()[:,:,int(mesh.dimension[1]/2)]\n", + " # create a plot of the flux relative error\n", + " plt.imshow(\n", + " flux_std_dev.T,\n", + " origin=\"lower\",\n", + " extent=tally_mesh_extent,\n", + " norm=LogNorm(),\n", + " )\n", + " plt.title(\"Flux Std. Dev.\")\n", + "\n", + " wws=openmc.hdf5_to_wws(weight_window_filename)\n", + " ww = wws[0] # get the one and only mesh tally\n", + " ww_mesh = ww.mesh # get the mesh that the weight window is mapped on\n", + " ww_mesh_extent = ww_mesh.bounding_box.extent['xy']\n", + " reshaped_ww_vals = ww.lower_ww_bounds.reshape(mesh.dimension)\n", + " print('reshaped_ww_vals.shape', reshaped_ww_vals.shape)\n", + " # slice on XZ basis, midplane Y axis\n", + " slice_of_ww = reshaped_ww_vals[:,:,int(mesh.dimension[1]/2)]\n", + " plt.subplot(1, 3, 3)\n", + " plt.imshow(\n", + " slice_of_ww.T,\n", + " origin=\"lower\",\n", + " extent=ww_mesh_extent,\n", + " norm=LogNorm(),\n", + " )\n", + " plt.title(\"Weight Window lower bound\")\n", + "\n", + " plt.savefig(image_filename, bbox_inches=\"tight\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "647a81d8", + "metadata": {}, + "source": [ + "Now we make use of openmc.lib to control the simulation. Documentation on openmc.lib is here\n", + "https://docs.openmc.org/en/stable/pythonapi/capi.html\n", + "\n", + "We run 5 iterations with each iteration improving the weight window." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9159a336", + "metadata": {}, + "outputs": [], + "source": [ + "with openmc.lib.run_in_memory():\n", + "\n", + " # loads up a live pointer to the tally with id=55, at this stage the tally is empty\n", + " tally = openmc.lib.tallies[55]\n", + "\n", + " # makes weight windows from the tally, at this stage the values are empty\n", + " wws = openmc.lib.WeightWindows.from_tally(tally, particle='neutron')\n", + "\n", + " # 5 iterations of weight window improvements\n", + " for i in range(5):\n", + "\n", + " # run the simulation\n", + " openmc.lib.run()\n", + "\n", + " # improves the weight window with the latest tally results\n", + " wws.update_magic(tally)\n", + "\n", + " # we write out the weight window map for plotting later\n", + " openmc.lib.export_weight_windows(filename=f'weight_windows{i}.h5')\n", + " # we write out the statepoint so that the tally can be plotted later\n", + " openmc.lib.statepoint_write(filename=f'statepoint_simulation_{i}.h5')\n", + "\n", + " # turns on the weight windows to ensure they are used\n", + " openmc.lib.settings.weight_windows_on = True\n", + "\n", + " # creates a plot of the flux, std_dev and weight window\n", + " plot_mesh_tally_and_weight_window(\n", + " f'statepoint_simulation_{i}.h5',\n", + " f'weight_windows{i}.h5',\n", + " f'plot_{i}.png'\n", + " )" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "e293ce43", + "metadata": {}, + "source": [] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "79ec2340", + "metadata": {}, + "source": [] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "0c6548fc-a538-432f-9694-3092db42b1b4", + "metadata": {}, + "source": [ + "The iterative improvment of the flux / standard deviation / weight windows with each simulation run can be seen when we plot all the simulation results one after each other." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b89dd15-1596-4fb8-9988-9b37d38b9b7c", + "metadata": {}, + "outputs": [], + "source": [ + "from PIL import Image\n", + "images = [Image.open(x) for x in [f'plot_{c}.png' for c in range(5)]]\n", + "widths, heights = zip(*(i.size for i in images))\n", + "\n", + "total_height = sum(heights)\n", + "max_width = max(widths)\n", + "\n", + "new_im = Image.new('RGB', (max_width, total_height))\n", + "y_offset = 0\n", + "for im in images:\n", + " new_im.paste(im, (0,y_offset))\n", + " y_offset += im.size[1]\n", + "new_im.save('flux_std-dev_ww_for_all_simulations_reset.png')\n", + "new_im" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d6daa0ba-5a5b-4701-aeaa-d57db34fb485", + "metadata": {}, + "source": [ + "Learning outcome\n", + "\n", + "Weight windows can be incrementally improved by running subsequent simulations.\n", + "\n", + "Running lots of small simulations where the weight window can incrementally improve the wieght window and can yields better results that a big single simulation to generate weight windows and a single big simulation to make use of the weight windows.\n", + "\n", + "Doing this iteration with openmc.lib means we don't need to reload the nuclear data between simulations which saves time.\n", + "\n", + "Additionally we have access to openmc.lib methods which are nessecary for updating the weight window with the MAGIC method and exporting the weight window to a h5 file.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74b59746-2e58-45b4-baeb-37b69244d145", + "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.8.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tasks/task_13_variance_reduction/4_sphere_iterative_per_batch_ww.ipynb b/tasks/task_13_variance_reduction/4_sphere_iterative_per_batch_ww.ipynb new file mode 100644 index 00000000..334bdf93 --- /dev/null +++ b/tasks/task_13_variance_reduction/4_sphere_iterative_per_batch_ww.ipynb @@ -0,0 +1,399 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "db6c76be-284f-4355-88aa-b4c89b9bbcea", + "metadata": {}, + "source": [ + "This example has a sphere of concrete with a second smaller shell of concrete\n", + "surrounding the sphere.\n", + "\n", + "The first simulation is analog with no variance reduction / weight windows.\n", + "This simulation shows that not many neutrons get to the shell and the \n", + "consequently the neutron spectra on the shell cell is unresolved. Additional\n", + "batches improve the neutron spectra but it is clear that it would take many\n", + "batches to get a reasonable neutron spectra.\n", + "\n", + "The second simulation makes use of a variance reduction method called weight\n", + "windows. The value of the weight windows is assigned using the MAGIC method.\n", + "https://scientific-publications.ukaea.uk/papers/application-of-novel-global-variance-reduction-methods-to-fusion-radiation-transport-problems/\n", + "The value of the weight windows are updated with each simulated batch and as\n", + "the simulation runs for longer the weight windows improve gradually as does\n", + "spectra tally." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecaafb63-bd33-4597-baf3-8e8cccd32f9c", + "metadata": {}, + "outputs": [], + "source": [ + "import openmc\n", + "# Note this example makes use of OpenMC lib which provides python bindings to\n", + "# the C/C++ methods in OpenMC and allows more direct control of the Monte Carlo\n", + "# simulation. In this example we iterate through the batches and access the\n", + "# tally result each time.\n", + "# Link to openmc.lib documentation https://docs.openmc.org/en/stable/pythonapi/capi.html\n", + "import openmc.lib\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "1e8755e5-f1f8-4c4f-b142-eb8b1e6e6f84", + "metadata": {}, + "source": [ + "This makes concrete which is the shielding material for this simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b26aff6f-e542-4c62-af29-9a4f7a793733", + "metadata": {}, + "outputs": [], + "source": [ + "mat_concrete = openmc.Material()\n", + "mat_concrete.add_element(\"H\",0.168759)\n", + "mat_concrete.add_element(\"C\",0.001416)\n", + "mat_concrete.add_element(\"O\",0.562524)\n", + "mat_concrete.add_element(\"Na\",0.011838)\n", + "mat_concrete.add_element(\"Mg\",0.0014)\n", + "mat_concrete.add_element(\"Al\",0.021354)\n", + "mat_concrete.add_element(\"Si\",0.204115)\n", + "mat_concrete.add_element(\"K\",0.005656)\n", + "mat_concrete.add_element(\"Ca\",0.018674)\n", + "mat_concrete.add_element(\"Fe\",0.00426)\n", + "mat_concrete.set_density(\"g/cm3\", 2.3)\n", + "\n", + "my_materials = openmc.Materials([mat_concrete])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "78616dec-bac6-405e-bcfc-b0724670d209", + "metadata": {}, + "source": [ + "This makes the spherical geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cb35d51-3d94-48a3-98fa-804a3195383a", + "metadata": {}, + "outputs": [], + "source": [ + "# surfaces\n", + "surf1 = openmc.Sphere(r=170)\n", + "outer_surface = openmc.Sphere(r=200, boundary_type=\"vacuum\")\n", + "\n", + "# regions\n", + "region_1 = -surf1\n", + "region_2 = -outer_surface & +surf1\n", + "\n", + "# cells\n", + "cell_1 = openmc.Cell(region=region_1)\n", + "cell_1.fill = mat_concrete\n", + "cell_2 = openmc.Cell(region=region_2)\n", + "cell_2.fill = mat_concrete\n", + "\n", + "# settings\n", + "my_settings = openmc.Settings()\n", + "\n", + "my_geometry = openmc.Geometry([cell_1, cell_2])\n", + "\n", + "my_geometry.root_universe.plot(basis='xy', color_by='cell') \n", + "plt.savefig('geometry_top_down_view.png', bbox_inches=\"tight\")" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b8be07cf-a90c-4fe9-ab86-44e26b0e6bfa", + "metadata": {}, + "source": [ + "This code makes a 14MeV neutron point source" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4868792-6ee6-440e-982c-4698a215c65b", + "metadata": {}, + "outputs": [], + "source": [ + "source = openmc.IndependentSource()\n", + "source.space = openmc.stats.Point((0.0, 0.0, 0.0))\n", + "source.angle = openmc.stats.Isotropic()\n", + "source.energy = openmc.stats.Discrete([14e6], [1.0])\n", + "source.particle = \"neutron\"" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "d70dd451-02ff-44d6-9628-84ed2753b578", + "metadata": {}, + "source": [ + "This section makes the simulation settings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f11ceb32-f439-4ecd-a99f-fbf0471d9405", + "metadata": {}, + "outputs": [], + "source": [ + "my_settings = openmc.Settings()\n", + "my_settings.run_mode = \"fixed source\"\n", + "my_settings.source = source\n", + "my_settings.particles = 50000\n", + "my_settings.batches = 5\n", + "# the mesh tallies produce large tallies.out files so this output setting avoids writing the tallies.out and saves time\n", + "my_settings.output = {'tallies': False}" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "e4a965a3-0c7e-4634-9f91-16b0a4325b8d", + "metadata": {}, + "source": [ + "This section makes the tallies.\n", + "\n", + "We have a spherical mesh tally for the getting the flux. This is used to generate the weight windows.\n", + "\n", + "We also have a neutron spectra tally to show how the neutrons energy distribution in the outer shell cell is resolved." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a51cd616-b0ff-4322-82db-4608f3082611", + "metadata": {}, + "outputs": [], + "source": [ + "my_tallies = openmc.Tallies()\n", + "\n", + "# This spherical mesh tally is used for generating the weight windows.\n", + "mesh = openmc.SphericalMesh()\n", + "mesh.r_grid = np.linspace(0, outer_surface.r, 5000)\n", + "mesh_filter = openmc.MeshFilter(mesh)\n", + "flux_tally_for_ww = openmc.Tally(name=\"flux tally\")\n", + "flux_tally_for_ww.filters = [mesh_filter]\n", + "flux_tally_for_ww.scores = [\"flux\"]\n", + "flux_tally_for_ww.id = 42\n", + "my_tallies.append(flux_tally_for_ww)\n", + "\n", + "# This spectrum tally is on the outer shell and shows then energy distribution\n", + "# of neutrons present in the cell.\n", + "energy_filter = openmc.EnergyFilter.from_group_structure('CCFE-709')\n", + "surface_filter = openmc.CellFilter(cell_2)\n", + "outer_surface_spectra_tally = openmc.Tally(name='outer_surface_spectra_tally')\n", + "outer_surface_spectra_tally.scores = ['current']\n", + "outer_surface_spectra_tally.filters = [surface_filter, energy_filter]\n", + "outer_surface_spectra_tally.id = 12\n", + "my_tallies.append(outer_surface_spectra_tally)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "f62aacde-8f5e-4fdb-82bd-58f43037b6ac", + "metadata": {}, + "source": [ + "creates and exports the model to an xml file. When using openmc.lib this\n", + "export is needed as we don't use the normal model.run() method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "540dde33-25a4-45d9-859f-11ab2e222cfe", + "metadata": {}, + "outputs": [], + "source": [ + "model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies)\n", + "model.export_to_xml()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "580415f1-d56e-4a4a-b1b4-1ef30aeb92e4", + "metadata": {}, + "source": [ + "Creates a plotting figure that we will build up with the simulations results. Initially this will be empty but we shall simulate and populate these rows and columns with neutron spectra in the outer shell cell." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec4a2ecf-ec9e-43a3-b2f1-249ef97aae0e", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(my_settings.batches, 2, sharex=True, sharey=True)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "862c2695-7aec-43e9-9461-0689ceed2482", + "metadata": {}, + "source": [ + "We run the model in analog mode batch by batch. Each time we plot the spectra\n", + "tally result. The spectra tally will gradually to get better with each batch\n", + "as the batches combine to continually improve the result." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ddb7274d-0782-4462-acf1-edf8f99895e6", + "metadata": {}, + "outputs": [], + "source": [ + "# this context manager helps close openmc lib when the code indent closes\n", + "with openmc.lib.run_in_memory():\n", + "\n", + " # gets a live pointer to the tally, this updates as the tally is accumulated\n", + " spectra_tally = openmc.lib.tallies[outer_surface_spectra_tally.id]\n", + " \n", + " # simulation_init is needed prior to iter_batches\n", + " openmc.lib.simulation_init()\n", + "\n", + " # loops through each batch getting the latest tally result and plotting it\n", + " for counter, batch in enumerate(openmc.lib.iter_batches()):\n", + "\n", + " axs[counter][0].step(energy_filter.values[:-1], spectra_tally.mean.flatten())\n", + " axs[counter][0].set_title(f'Batch {counter+1}')\n", + " axs[counter][0].set_yscale('log')\n", + " axs[counter][0].set_xscale('log')\n", + "\n", + " openmc.lib.simulation_finalize()" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "e000825a-033c-4663-8e9f-70390ed54b8e", + "metadata": {}, + "source": [ + "Now we simulate the same model but with weight windows that are improved on each batch" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "05c69e62-54f8-45b3-919f-f20bb1d3cc80", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7701329-7fba-480c-98de-661969d743b4", + "metadata": {}, + "outputs": [], + "source": [ + "# originally we had 2000 particles per batch\n", + "# on my computer the analog simulation ran with 12470 particles/second\n", + "# and the weight windows simulation that comes next runs with 87 particles/second\n", + "# therefore we are going to decrease the settings.particles so that both simulations\n", + "# get the same amount of compute time. These numbers might differ for your computer.\n", + "model.settings.particles = int(model.settings.particles*(87/12470))\n", + "\n", + "# export the model again so the particles is updated\n", + "model.export_to_xml()\n", + "with openmc.lib.run_in_memory():\n", + "\n", + " # gets a live pointer to the mesh tally that we use to generate the \n", + " ww_tally = openmc.lib.tallies[flux_tally_for_ww.id]\n", + " # generates a weight window from the tally (which is currently empty)\n", + " wws = openmc.lib.WeightWindows.from_tally(ww_tally, particle='neutron')\n", + " \n", + " # gets a live pointer to the spectra tally that we will plot with each batch\n", + " spectra_tally = openmc.lib.tallies[outer_surface_spectra_tally.id]\n", + "\n", + " # turn the weight windows on\n", + " openmc.lib.settings.weight_windows_on = True\n", + "\n", + " openmc.lib.simulation_init()\n", + " for counter, batch in enumerate(openmc.lib.iter_batches()):\n", + "\n", + " # updates the weight window with the latest mesh tally flux results \n", + " wws.update_magic(ww_tally)\n", + "\n", + " # plots the spectra tally for the batch\n", + " axs[counter][1].step(energy_filter.values[:-1], spectra_tally.mean.flatten())\n", + " axs[counter][1].set_title(f'Batch {counter+1}')\n", + " axs[counter][1].set_yscale('log')\n", + " axs[counter][1].set_xscale('log')\n", + "\n", + " openmc.lib.simulation_finalize()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73db1d34-39f5-4065-9113-fe27cf654c1a", + "metadata": {}, + "outputs": [], + "source": [ + "# sets titles, labels and saves the plot\n", + "axs[0][0].set_title('Analog simulation')\n", + "axs[0][1].set_title('Iterative weight windows simulation')\n", + "axs[4][1].set_xlabel(f'Energy [eV]')\n", + "axs[4][0].set_xlabel(f'Energy [eV]')\n", + "fig.savefig('improving_spectra_with_ww_iteration.jpg', bbox_inches=\"tight\")\n", + "# plt.show()\n", + "from IPython.display import Image\n", + "Image('improving_spectra_with_ww_iteration.jpg')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "b2c74788-8de5-4be7-9263-7f988fc43890", + "metadata": {}, + "source": [ + "In the left hand side colum we see that just running a non varience reduction simulation and only occational neutrons get to the outer shell cell to contribute to the tally. This would take a long time to converge the spectra.\n", + "\n", + "However with weight windows we see the tally results are developing relatively quickly and we can start to see some spectra structure after just a few batches.\n", + "\n", + "We have previously seen flux maps benefit from weight windows but this example showed that spectra simulations can also be improved.\n", + "\n", + "The other difference is that we are improving the weight window each batch instead of each simulation. This is quite fine grain evolution of the weight windows. The finaly tally result ends up getting contributions from several batches where each batch has different weight windows." + ] + } + ], + "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.8.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tasks/task_13_variance_reduction/README.md b/tasks/task_13_variance_reduction/README.md new file mode 100644 index 00000000..48614670 --- /dev/null +++ b/tasks/task_13_variance_reduction/README.md @@ -0,0 +1,13 @@ +Varience reduction examples + +Currently OpenMC supports two types of variance reduction. +A detailed description of each method can be found in the [documentation](https://docs.openmc.org/en/stable/methods/neutron_physics.html?highlight=survival#variance-reduction-techniques). + +The workshop contains the following variance reduction examples: + +| Filename | variance reduction technique | geometry | mesh type | +|---|---|---|---| +| 1_shielded_room_survival_biasing.py | survival_biasing | shielded bunker | RegularMesh | Flux map | air space and concrete | +| 2_shielded_room_single_ww.ipynb | weight windows | sphere | RegularMesh | air space and concrete | +| 3_sphere_iterative_per_run_ww.py | weight windows | cube | RegularMesh | Water | +| 4_sphere_iterative_per_batch_ww.py | weight windows | sphere | SphericalMesh | concrete | \ No newline at end of file diff --git a/tasks/task_13_variance_reduction/generate_iterative_weight_windows.ipynb b/tasks/task_13_variance_reduction/generate_iterative_weight_windows.ipynb deleted file mode 100644 index c3cf5c60..00000000 --- a/tasks/task_13_variance_reduction/generate_iterative_weight_windows.ipynb +++ /dev/null @@ -1,259 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "d4d248c5-e7cd-428f-b099-061a1328853e", - "metadata": {}, - "source": [ - "# variance reduction techniques\n", - "\n", - "## Iteratively create and utilising a wight window to accelerate the deep shielding simulations" - ] - }, - { - "cell_type": "markdown", - "id": "7b02559d-a4f3-4e6f-83f4-1aa030bf59f5", - "metadata": {}, - "source": [ - "This example implements the Magic method of weight window generation. The theory of weight windows and the method are best described by the paper here https://scientific-publications.ukaea.uk/wp-content/uploads/Published/INTERN1.pdf\n", - "\n", - "In this tutorial we shall focus on generating a weight window to accelerate the simulation of particles through a shield. Then using that weight window to generate a new weight window and repeating the process to obtain a final weight window.\n", - "\n", - "First we must make a model. This is kept as simple as possible as the focus of this notebook is on generating and then using a weight window.\n", - "\n", - "The model is a single sphere of 200 cm radius filled with water and containing a 14MeV point source in the center" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "779490fe-0dfe-4523-982a-3881eff46e80", - "metadata": {}, - "outputs": [], - "source": [ - "import openmc\n", - "import openmc_weight_window_generator\n", - "# this is a minimal package that adds some functionality to openmc, namely it adds:\n", - "# - openmc.StatePoint.generate_wws which we use in a previous task\n", - "# - openmc.Model.generate_wws_magic_method which we use in this task\n", - "\n", - "\n", - "# creates a shielding material\n", - "water = openmc.Material(name='Water')\n", - "water.set_density('g/cc', 1.0)\n", - "water.add_element('H', 2)\n", - "water.add_element('O', 1)\n", - "materials = openmc.Materials([water])\n", - "\n", - "sphere1 = openmc.Sphere(r=200, boundary_type='vacuum')\n", - "\n", - "region1 = -sphere1\n", - "\n", - "cell1 = openmc.Cell(fill=water, region=region1)\n", - "\n", - "geometry = openmc.Geometry([cell1])\n", - "\n", - "source = openmc.Source()\n", - "source.space = openmc.stats.Point((0.0, 0.0, 0.0))\n", - "source.angle = openmc.stats.Isotropic()\n", - "source.energy = openmc.stats.Discrete([14e6], [1.0])\n", - "source.particle = 'neutron'\n", - "\n", - "my_settings = openmc.Settings()\n", - "my_settings.run_mode = 'fixed source'\n", - "my_settings.source = source\n", - "my_settings.particles = 10000\n", - "my_settings.batches = 10\n", - "# the mesh tallies produce large tallies.out files so this output setting avoids writing the tallies.out and saves time\n", - "my_settings.output = {'tallies': False}\n", - "\n", - "model = openmc.Model(geometry, materials, my_settings)" - ] - }, - { - "cell_type": "markdown", - "id": "980c9abe-cb26-4bf8-998d-08bc46b0cac7", - "metadata": {}, - "source": [ - "Now we can plot the simple geometry" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "da6728e8-a947-4dee-8e20-04a980019d33", - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "plt.figure(figsize=(10,10))\n", - "model.geometry.root_universe.plot(pixels=(600, 600))\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "c3a3547e-f9d7-4b53-806b-0c84648e4df6", - "metadata": {}, - "source": [ - "Now we shall add a regular mesh tally over the the geometry.\n", - "\n", - "The mesh will be used to record the neutron flux in each mesh voxel" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "44dba5ef-ec9c-4280-8fc1-b9191b5c62ea", - "metadata": {}, - "outputs": [], - "source": [ - "mesh = openmc.RegularMesh()\n", - "mesh.lower_left = model.geometry.bounding_box[0]\n", - "mesh.upper_right = model.geometry.bounding_box[1]\n", - "mesh.dimension = (50, 50, 50)\n", - "\n", - "mesh_filter = openmc.MeshFilter(mesh)\n", - "\n", - "flux_tally = openmc.Tally(name='flux tally')\n", - "flux_tally.filters = [mesh_filter]\n", - "flux_tally.scores = ['flux']\n", - "\n", - "# adds the mesh tally to the model\n", - "model.tallies = [flux_tally]\n", - "\n", - "output_file = model.run()" - ] - }, - { - "cell_type": "markdown", - "id": "a2404758-0b31-442f-84c5-9e5cc2a455c3", - "metadata": {}, - "source": [ - "Now we can plot the flux and the standard deviation of the flux tally to see how far into the shield the neutrons got. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a5997678-95a6-4178-99d2-b66df38cfc1c", - "metadata": {}, - "outputs": [], - "source": [ - "from matplotlib.colors import LogNorm\n", - "\n", - "with openmc.StatePoint(output_file) as sp:\n", - " flux_tally = sp.get_tally(id=flux_tally.id)\n", - "\n", - "llc, urc = model.geometry.bounding_box\n", - "\n", - "fig, (ax1, ax2) = plt.subplots(1,2)\n", - "fig.suptitle('Flux and std. dev. without weight window')\n", - "\n", - "slice_index = int(mesh.dimension[-1]/2)\n", - "\n", - "# create a plot of the mean values\n", - "flux_mean = flux_tally.mean.reshape(*mesh.dimension)\n", - "img1 = ax1.imshow(flux_mean[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]), norm=LogNorm())\n", - "ax1.set_title('Flux Mean')\n", - "plt.colorbar(img1, ax=ax1, fraction=0.046)\n", - "img1.set_clim(vmin=1e-30, vmax=1.0)\n", - "\n", - "# create a plot of the flux relative error\n", - "flux_rel_err = flux_tally.get_values(value='rel_err').reshape(*mesh.dimension)\n", - "img2 = ax2.imshow(flux_rel_err[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]))\n", - "ax2.set_title('Flux Rel. Err.')\n", - "plt.colorbar(img2, ax=ax2, fraction=0.046)\n", - "# ax2.set_colorbar(img2, ax=ax2)\n", - "img2.set_clim(vmin=0.0, vmax=1.0)\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "e8bc03db-f38b-4c85-ae11-f222f4fe80c7", - "metadata": {}, - "source": [ - "Now we shall run the simulation several times, with each iteration we shall use the flux tally to produce a weight window. Then the weight window will be used in the subsequent simulation." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5194f55d-ab70-44d5-a7c7-5dbb3f6efa8f", - "metadata": {}, - "outputs": [], - "source": [ - "# change this to 5 iterations and 50000 max splits to see a nice improvement in weight windows\n", - "model.generate_wws_magic_method(tally=flux_tally, iterations=3, max_split=300)" - ] - }, - { - "cell_type": "markdown", - "id": "76c22db3-6959-4ea7-b819-3ba49ac19eab", - "metadata": {}, - "source": [ - "Now we can plot the flux and the standard deviation of the flux tally to see how far into the shield the neutrons got. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6b81b0c7-6313-43fc-895c-55494870508d", - "metadata": {}, - "outputs": [], - "source": [ - "from matplotlib.colors import LogNorm\n", - "\n", - "with openmc.StatePoint(output_file) as sp:\n", - " flux_tally = sp.get_tally(id=flux_tally.id)\n", - "\n", - "llc, urc = model.geometry.bounding_box\n", - "\n", - "fig, (ax1, ax2) = plt.subplots(1,2)\n", - "fig.suptitle('Flux and std. dev. without weight window')\n", - "\n", - "slice_index = int(mesh.dimension[-1]/2)\n", - "\n", - "# create a plot of the mean values\n", - "flux_mean = flux_tally.mean.reshape(*mesh.dimension)\n", - "img1 = ax1.imshow(flux_mean[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]), norm=LogNorm())\n", - "ax1.set_title('Flux Mean')\n", - "plt.colorbar(img1, ax=ax1, fraction=0.046)\n", - "img1.set_clim(vmin=1e-30, vmax=1.0)\n", - "\n", - "# create a plot of the flux relative error\n", - "flux_rel_err = flux_tally.get_values(value='rel_err').reshape(*mesh.dimension)\n", - "img2 = ax2.imshow(flux_rel_err[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]))\n", - "ax2.set_title('Flux Rel. Err.')\n", - "plt.colorbar(img2, ax=ax2, fraction=0.046)\n", - "# ax2.set_colorbar(img2, ax=ax2)\n", - "img2.set_clim(vmin=0.0, vmax=1.0)\n", - "\n", - "plt.show()" - ] - } - ], - "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.8.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/tasks/task_13_variance_reduction/generate_single_ww_and_apply.ipynb b/tasks/task_13_variance_reduction/generate_single_ww_and_apply.ipynb deleted file mode 100644 index 39914f4c..00000000 --- a/tasks/task_13_variance_reduction/generate_single_ww_and_apply.ipynb +++ /dev/null @@ -1,347 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "72e0b25d-e541-4e8d-8805-b984374ee53d", - "metadata": {}, - "source": [ - "# variance reduction techniques\n", - "\n", - "## Creating and utilising a wight window to accelerate the deep shielding simulations" - ] - }, - { - "cell_type": "markdown", - "id": "f00eddb1-1e3a-4973-ba0c-d8feeb2a6704", - "metadata": {}, - "source": [ - "This example implements a single step of the Magic method of weight window generation. The theory of weight windows and the method are best described by the paper here https://scientific-publications.ukaea.uk/wp-content/uploads/Published/INTERN1.pdf\n", - "\n", - "In this tutorial we shall focus on generating a weight window to accelerate the simulation of particles through a shield.\n", - "\n", - "First we must make a model. This is kept as simple as possible as the focus of this notebook is on generating and then using a weight window.\n", - "\n", - "The model is a single sphere of 200 cm radius filled with water and containing a 14MeV point source in the center" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3448d1db-4328-42e3-8960-50d53896f541", - "metadata": {}, - "outputs": [], - "source": [ - "import openmc\n", - "\n", - "# creates a shielding material\n", - "water = openmc.Material(name='Water')\n", - "water.set_density('g/cc', 1.0)\n", - "water.add_element('H', 2)\n", - "water.add_element('O', 1)\n", - "materials = openmc.Materials([water])\n", - "\n", - "sphere1 = openmc.Sphere(r=200, boundary_type='vacuum')\n", - "\n", - "region1 = -sphere1\n", - "\n", - "cell1 = openmc.Cell(fill=water, region=region1)\n", - "\n", - "geometry = openmc.Geometry([cell1])\n", - "\n", - "source = openmc.Source()\n", - "source.space = openmc.stats.Point((0.0, 0.0, 0.0))\n", - "source.angle = openmc.stats.Isotropic()\n", - "source.energy = openmc.stats.Discrete([14e6], [1.0])\n", - "source.particle = 'neutron'\n", - "\n", - "my_settings = openmc.Settings()\n", - "my_settings.run_mode = 'fixed source'\n", - "my_settings.source = source\n", - "my_settings.particles = 10000\n", - "my_settings.batches = 10\n", - "# the mesh tallies produce large tallies.out files so this output setting avoids writing the tallies.out and saves time\n", - "my_settings.output = {'tallies': False}\n", - "\n", - "model = openmc.model.Model(geometry, materials, my_settings)" - ] - }, - { - "cell_type": "markdown", - "id": "8179e89a-304b-4684-91a9-8a908b75e8cd", - "metadata": {}, - "source": [ - "Now we can plot the simple geometry" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "47cb21c7-0afa-446c-9a6b-76c8dfb2f93e", - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "plt.figure(figsize=(10,10))\n", - "model.geometry.root_universe.plot(pixels=(600, 600))\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "424b234c-48a2-4d1b-b374-f69430e464db", - "metadata": {}, - "source": [ - "Now we shall add a regular mesh tally over the the geometry.\n", - "\n", - "The mesh will be used to record the neutron flux in each mesh voxel" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a72dc5e3-cb8e-41c2-9cfb-f8f7489c7f44", - "metadata": {}, - "outputs": [], - "source": [ - "mesh = openmc.RegularMesh()\n", - "mesh.lower_left = model.geometry.bounding_box[0]\n", - "mesh.upper_right = model.geometry.bounding_box[1]\n", - "mesh.dimension = (50, 50, 50)\n", - "\n", - "mesh_filter = openmc.MeshFilter(mesh)\n", - "\n", - "flux_tally = openmc.Tally(name='flux tally')\n", - "flux_tally.filters = [mesh_filter]\n", - "flux_tally.scores = ['flux']\n", - "\n", - "# adds the mesh tally to the model\n", - "model.tallies = [flux_tally]" - ] - }, - { - "cell_type": "markdown", - "id": "458867ab-5b7b-400a-bbf4-c290988a99da", - "metadata": {}, - "source": [ - "Now we shall run the simulation and record flux in each mesh voxel" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1c5b84c7-0fe5-48d5-b2b7-58ce8cea57f1", - "metadata": {}, - "outputs": [], - "source": [ - "output_file = model.run()" - ] - }, - { - "cell_type": "markdown", - "id": "e54d695d-1f35-477c-8506-fcbca57b179a", - "metadata": {}, - "source": [ - "Now we can plot the flux and the standard deviation of the flux tally to see how far into the shield the neutrons got. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0494f6db-ec78-43f2-a6be-c4c55c7e7d0a", - "metadata": {}, - "outputs": [], - "source": [ - "from matplotlib.colors import LogNorm\n", - "\n", - "with openmc.StatePoint(output_file) as sp:\n", - " flux_tally = sp.get_tally(id=flux_tally.id)\n", - "\n", - "llc, urc = model.geometry.bounding_box\n", - "\n", - "fig, (ax1, ax2) = plt.subplots(1,2)\n", - "fig.suptitle('Flux and std. dev. without weight window')\n", - "\n", - "slice_index = int(mesh.dimension[-1]/2)\n", - "\n", - "# create a plot of the mean values\n", - "flux_mean = flux_tally.mean.reshape(*mesh.dimension)\n", - "img1 = ax1.imshow(flux_mean[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]), norm=LogNorm())\n", - "ax1.set_title('Flux Mean')\n", - "plt.colorbar(img1, ax=ax1, fraction=0.046)\n", - "img1.set_clim(vmin=1e-30, vmax=1.0)\n", - "\n", - "# create a plot of the flux relative error\n", - "flux_rel_err = flux_tally.get_values(value='rel_err').reshape(*mesh.dimension)\n", - "img2 = ax2.imshow(flux_rel_err[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]))\n", - "ax2.set_title('Flux Rel. Err.')\n", - "plt.colorbar(img2, ax=ax2, fraction=0.046)\n", - "# ax2.set_colorbar(img2, ax=ax2)\n", - "img2.set_clim(vmin=0.0, vmax=1.0)\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "f0b0c196-0829-4b36-abec-8fec944c8f2b", - "metadata": {}, - "source": [ - "As this flux map tells us where the neutrons go we can use it to create a wieght window that promotes neutron transport in areas they normally don't reach." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "03cc2b9f-5421-41d3-9f0d-4d43cf82977f", - "metadata": {}, - "outputs": [], - "source": [ - "import openmc_weight_window_generator\n", - "# this is a minimal package that adds some functionality to openmc, namely it adds:\n", - "# - openmc.StatePoint.generate_wws which we use in this task\n", - "# - openmc.Model.generate_wws_magic_method which we use in the next task\n", - "\n", - "sp_file = openmc.StatePoint(output_file)\n", - "# this generates an openmc.WeightWindow object from the flux tally\n", - "weight_windows = sp_file.generate_wws(tally=flux_tally, rel_err_tol=0.7)" - ] - }, - { - "cell_type": "markdown", - "id": "cfb6e00c-aada-4951-9cdf-022477defa6c", - "metadata": {}, - "source": [ - "The weight window generated uses the same mesh as the flux tally, uses the flux to generate lower_ww_bounds, sets the upper_bound_ratio to 5 (rule of thumb used here) and sets the max_split to 1_000_000. These can all be changed to customise the weight window but are reasonable defaults.\n", - "\n", - "We can plot the lower_ww_bounds of the generated weight window to see how it changes over the geometry" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c2cc0474-ab96-4db9-a261-ff0b1fe9256b", - "metadata": {}, - "outputs": [], - "source": [ - "plt.imshow(weight_windows[0].lower_ww_bounds[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]), norm=LogNorm())\n", - "plt.title('lower_ww_bounds')\n", - "plt.colorbar()" - ] - }, - { - "cell_type": "markdown", - "id": "02af110d-9ea2-4782-a355-bbd6822cc394", - "metadata": {}, - "source": [ - "Now we can rerun the simulation but this time using the weight window to push the particles further into the geometry" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4c5f7fda-870f-4c84-8820-e87966de2ae0", - "metadata": {}, - "outputs": [], - "source": [ - "#deletes the old output files\n", - "!rm summary.h5\n", - "!rm statepoint.*.h5\n", - "\n", - "\n", - "model.settings.weight_windows = weight_windows\n", - "model.settings.max_split = 1_000 # might want to increase this to 1_000_000 and get a better result \n", - "model.run()" - ] - }, - { - "cell_type": "markdown", - "id": "4eee750a-8913-4dcc-a528-fcea1b74a854", - "metadata": {}, - "source": [ - "We can now plot the flux and standard devation of the flux for the simulation that used weight windows and see that the particles penetrated further into the shield and that the standard deviation has been reduced." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fb2b523c-f925-4670-927b-22ce695de56e", - "metadata": {}, - "outputs": [], - "source": [ - "with openmc.StatePoint(output_file) as sp:\n", - " flux_tally = sp.get_tally(id=flux_tally.id)\n", - "\n", - "fig, (ax1, ax2) = plt.subplots(1,2)\n", - "fig.tight_layout() \n", - "fig.suptitle('Flux and std. dev. with weight window')\n", - "\n", - "# create a plot of the mean values\n", - "flux_mean = flux_tally.mean.reshape(*mesh.dimension)\n", - "img1 = ax1.imshow(flux_mean[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]), norm=LogNorm())\n", - "ax1.set_title('Flux Mean')\n", - "plt.colorbar(img1, ax=ax1, fraction=0.046)\n", - "img1.set_clim(vmin=1e-30, vmax=1.0)\n", - "\n", - "# create a plot of the flux relative error\n", - "flux_rel_err = flux_tally.get_values(value='rel_err').reshape(*mesh.dimension)\n", - "img2 = ax2.imshow(flux_rel_err[slice_index], origin='lower', extent=(llc[0], urc[0], llc[1], urc[1]))\n", - "ax2.set_title('Flux Rel. Err.')\n", - "plt.colorbar(img2, ax=ax2, fraction=0.046)\n", - "# ax2.set_colorbar(img2, ax=ax2)\n", - "img2.set_clim(vmin=0.0, vmax=1.0)\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "2285d40b-2e8e-4f52-ac14-c742599cf08f", - "metadata": {}, - "source": [ - "Notice that the particles now get further into the shielding and the error has been reduce across the simulation.\n", - "\n", - "This is not exactly a fair comparison as the second simulation takes a little longer to run. To make it fairer we could use a trigger to stop each simulation after the same amount of time. However that would complicate the example." - ] - }, - { - "cell_type": "markdown", - "id": "4e9a1ae6-7ab7-4bdb-be1a-39fa6a8ca4c2", - "metadata": {}, - "source": [ - "Learning Outcomes:\n", - "* Weight windows can be useful for accelerating deep shielding simulations where particles\n", - "* Weight windows can be generated from a neutron flux field\n", - "* The MAGIC method is a popular method of generating weight windows" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "91d6b1bb-8b69-46fb-9a09-91727db67c29", - "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.8.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/tasks/task_13_variance_reduction/iterative_spherical_mesh_weight_windows.py b/tasks/task_13_variance_reduction/iterative_spherical_mesh_weight_windows.py deleted file mode 100644 index a633a90b..00000000 --- a/tasks/task_13_variance_reduction/iterative_spherical_mesh_weight_windows.py +++ /dev/null @@ -1,115 +0,0 @@ -# This example has a sphere of water showing how to increase the depth of -# neutron transport through the water. -# A series of simulations are performed that iteratively improve the weight window values -# The resulting neutron flux can be observed to propagate further through the geometry - -import openmc -import numpy as np -import matplotlib.pyplot as plt -import plotly.graph_objects as go -from pathlib import Path - - -# materials -mat_water = openmc.Material() -mat_water.add_element("H", 1) -mat_water.add_element("O", 2) -mat_water.set_density("g/cm3", 1.0) - -my_materials = openmc.Materials([mat_water]) - -# surfaces -outer_surface = openmc.Sphere(r=50*100, boundary_type="vacuum") - -# regions -region_1 = -outer_surface - -# cells -cell_1 = openmc.Cell(region=region_1) -cell_1.fill = mat_water - -# settings -my_settings = openmc.Settings() - -my_geometry = openmc.Geometry([cell_1]) - -source = openmc.Source() -source.space = openmc.stats.Point((0.0, 0.0, 0.0)) -source.angle = openmc.stats.Isotropic() -source.energy = openmc.stats.Discrete([14e6], [1.0]) -source.particle = "neutron" - -my_settings = openmc.Settings() -my_settings.run_mode = "fixed source" -my_settings.source = source -my_settings.particles = 10000 -my_settings.batches = 10 - -# tally -mesh = openmc.SphericalMesh() -mesh.r_grid = np.linspace(0, outer_surface.r, 100) - -mesh_filter = openmc.MeshFilter(mesh) - -flux_tally = openmc.Tally(name="flux tally") -flux_tally.filters = [mesh_filter] -flux_tally.scores = ["flux"] - -# adds the mesh tally to the model -my_tallies = openmc.Tallies() -my_tallies.append(flux_tally) - -import openmc_weight_window_generator -# import adds generate_wws_magic_method method to the model class - -model = openmc.Model(my_geometry, my_materials, my_settings, my_tallies) - - -all_wws = model.generate_wws_magic_method( - tally=flux_tally, iterations=50, - max_split=500_000, - output_dir="magic_ww", rel_err_tol=0.99 -) - -# plots the flux as a function of radius for each iteration -output_files = [Path("magic_ww") / str(c) / f"statepoint.{my_settings.batches}.h5" for c in range(1, len(all_wws))] -fig = go.Figure() -for i, output_file in enumerate(output_files): - with openmc.StatePoint(output_file) as sp: - flux_tally = sp.get_tally(name="flux tally") - fig.add_trace( - go.Scatter( - x=mesh.r_grid[1:], y=flux_tally.mean.flatten(), name=f"flux tally, iteration {i+1}" - ) - ) -fig.update_yaxes(type="log") -fig.update_layout(xaxis_title="Radius [cm]", yaxis_title="Flux") -fig.show() - -# plots the lower bound of the the weight window values as a function of radius for each iteration -fig = go.Figure() -for i, weight_windows in enumerate(all_wws): - fig.add_trace( - go.Scatter( - x=mesh.r_grid[1:], - y=weight_windows[0].lower_ww_bounds.flatten(), - name=f"lower ww bounds, iteration {i+1}", - ) - ) -fig.update_yaxes(type="log") -fig.update_layout(xaxis_title="Radius [cm]", yaxis_title="weight window lower bound value") -fig.show() - -# plots the upper bound of the the weight window values as a function of radius for each iteration -fig = go.Figure() -for i, weight_windows in enumerate(all_wws): - fig.add_trace( - go.Scatter( - x=mesh.r_grid[1:], - y=weight_windows[0].upper_ww_bounds.flatten(), - name=f"upper ww bounds, iteration {i+1}", - ) - ) -fig.update_yaxes(type="log") -fig.update_layout(xaxis_title="Radius [cm]", yaxis_title="weight window upper bound value") -fig.show() diff --git a/tasks/task_13_variance_reduction/shielded_room_survival_biasing.py b/tasks/task_13_variance_reduction/shielded_room_survival_biasing.py deleted file mode 100644 index bab67163..00000000 --- a/tasks/task_13_variance_reduction/shielded_room_survival_biasing.py +++ /dev/null @@ -1,269 +0,0 @@ -import openmc -import openmc_geometry_plot # adds extra plotting functions to openmc.Geometry object -import numpy as np -from matplotlib import pyplot as plt -from matplotlib.colors import LogNorm - - -air = openmc.Material(name="Air") -air.set_density("g/cc", 0.001205) -air.add_element("N", 0.784431) -air.add_element("O", 0.210748) -air.add_element("Ar", 0.0046) - -concrete = openmc.Material(name="concrete") -concrete.set_density("g/cm3", 7.874) -concrete.add_element("Fe", 1) - -materials = openmc.Materials([air, concrete]) - -width_a = 100 -width_b = 200 -width_c = 500 -width_d = 250 -width_e = 200 -width_f = 200 -width_g = 100 - -depth_a = 100 -depth_b = 200 -depth_c = 700 -depth_d = 600 -depth_e = 200 -depth_f = 100 - -height_j = 100 -height_k = 500 -height_l = 100 - - -xplane_0 = openmc.XPlane(x0=0, boundary_type="vacuum") -xplane_1 = openmc.XPlane(x0=xplane_0.x0 + width_a) -xplane_2 = openmc.XPlane(x0=xplane_1.x0 + width_b) -xplane_3 = openmc.XPlane(x0=xplane_2.x0 + width_c) -xplane_4 = openmc.XPlane(x0=xplane_3.x0 + width_d) -xplane_5 = openmc.XPlane(x0=xplane_4.x0 + width_e) -xplane_6 = openmc.XPlane(x0=xplane_5.x0 + width_f) -xplane_7 = openmc.XPlane(x0=xplane_6.x0 + width_g, boundary_type="vacuum") - -yplane_0 = openmc.YPlane(y0=0, boundary_type="vacuum") -yplane_1 = openmc.YPlane(y0=yplane_0.y0 + depth_a) -yplane_2 = openmc.YPlane(y0=yplane_1.y0 + depth_b) -yplane_3 = openmc.YPlane(y0=yplane_2.y0 + depth_c) -yplane_4 = openmc.YPlane(y0=yplane_3.y0 + depth_d) -yplane_5 = openmc.YPlane(y0=yplane_4.y0 + depth_e) -yplane_6 = openmc.YPlane(y0=yplane_5.y0 + depth_f, boundary_type="vacuum") - -zplane_1 = openmc.ZPlane(z0=0, boundary_type="vacuum") -zplane_2 = openmc.ZPlane(z0=zplane_1.z0 + height_j) -zplane_3 = openmc.ZPlane(z0=zplane_2.z0 + height_k) -zplane_4 = openmc.ZPlane(z0=zplane_3.z0 + height_l, boundary_type="vacuum") - -outside_left_region = +xplane_0 & -xplane_1 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_4 -wall_left_region = +xplane_1 & -xplane_2 & +yplane_2 & -yplane_4 & +zplane_2 & -zplane_3 -wall_right_region = +xplane_5 & -xplane_6 & +yplane_2 & -yplane_5 & +zplane_2 & -zplane_3 -wall_top_region = +xplane_1 & -xplane_4 & +yplane_4 & -yplane_5 & +zplane_2 & -zplane_3 -outside_top_region = +xplane_0 & -xplane_7 & +yplane_5 & -yplane_6 & +zplane_1 & -zplane_4 -wall_bottom_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_2 & +zplane_2 & -zplane_3 -outside_bottom_region = +xplane_0 & -xplane_7 & +yplane_0 & -yplane_1 & +zplane_1 & -zplane_4 -wall_middle_region = +xplane_3 & -xplane_4 & +yplane_3 & -yplane_4 & +zplane_2 & -zplane_3 -outside_right_region = +xplane_6 & -xplane_7 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_4 - -room_region = +xplane_2 & -xplane_3 & +yplane_2 & -yplane_4 & +zplane_2 & -zplane_3 -gap_region = +xplane_3 & -xplane_4 & +yplane_2 & -yplane_3 & +zplane_2 & -zplane_3 -corridor_region = +xplane_4 & -xplane_5 & +yplane_2 & -yplane_5 & +zplane_2 & -zplane_3 - -roof_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_5 & +zplane_1 & -zplane_2 -floor_region = +xplane_1 & -xplane_6 & +yplane_1 & -yplane_5 & +zplane_3 & -zplane_4 - -outside_left_cell = openmc.Cell(region=outside_left_region, fill=air, name="outside_left_cell") -outside_right_cell = openmc.Cell(region=outside_right_region, fill=air, name="outside_right_cell") -outside_top_cell = openmc.Cell(region=outside_top_region, fill=air, name="outside_top_cell") -outside_bottom_cell = openmc.Cell(region=outside_bottom_region, fill=air, name="outside_bottom_cell") -wall_left_cell = openmc.Cell(region=wall_left_region, fill=concrete) -wall_right_cell = openmc.Cell(region=wall_right_region, fill=concrete) -wall_top_cell = openmc.Cell(region=wall_top_region, fill=concrete) -wall_bottom_cell = openmc.Cell(region=wall_bottom_region, fill=concrete) -wall_middle_cell = openmc.Cell(region=wall_middle_region, fill=concrete) -room_cell = openmc.Cell(region=room_region, fill=air, name="room_cell") -gap_cell = openmc.Cell(region=gap_region, fill=air) -corridor_cell = openmc.Cell(region=corridor_region, fill=air) - -roof_cell = openmc.Cell(region=roof_region, fill=concrete) -floor_cell = openmc.Cell(region=floor_region, fill=concrete) - -materials = openmc.Materials([air, concrete]) -geometry = openmc.Geometry( - [ - outside_bottom_cell, - outside_top_cell, - outside_left_cell, - outside_right_cell, - wall_left_cell, - wall_right_cell, - wall_top_cell, - wall_bottom_cell, - wall_middle_cell, - room_cell, - gap_cell, - corridor_cell, - roof_cell, - floor_cell, - ] -) - -model = openmc.Model() -model.geometry = geometry - -# location of the point source -source_x = width_a + width_b + width_c * 0.5 -source_y = depth_a + depth_b + depth_c * 0.75 -source_z = height_j + height_k * 0.5 - -xlabel, ylabel = geometry.get_axis_labels(view_direction="z") -plt.xlabel(xlabel) -plt.ylabel(ylabel) - -plot_extent = geometry.bounding_box.extent['xy'] -data_slice = geometry.get_slice_of_material_ids(view_direction="z") -# plots the materials with randomly assigned colors -plt.imshow( - np.fliplr(data_slice), - extent=plot_extent, -) - -# plots the outline of the cells -plt.contour( - np.fliplr(data_slice), - origin="upper", - colors="k", - linestyles="solid", - linewidths=1, - extent=plot_extent, -) - -# plots the source location -plt.scatter(source_x, source_y, c="red") -plt.savefig('geometry_view.png') - -plt.clf() -plt.cla() - -xlabel, ylabel = geometry.get_axis_labels(view_direction="y") -plt.xlabel(xlabel) -plt.ylabel(ylabel) - - -plot_extent = geometry.bounding_box.extent['xz'] - -data_slice = geometry.get_slice_of_material_ids(view_direction="y") -# plots the materials with randomly assigned colors -plt.imshow( - np.fliplr(data_slice), - extent=plot_extent, -) - -# plots the outline of the cells -plt.contour( - np.fliplr(data_slice), - origin="upper", - colors="k", - linestyles="solid", - linewidths=1, - extent=plot_extent, -) - -# plots the source location -plt.scatter(source_x, source_z, c="red") -plt.savefig('geometry_view_2.png') - -mesh = openmc.RegularMesh().from_domain(geometry) -mesh.dimension = (100, 100, 1) - -mesh_filter = openmc.MeshFilter(mesh) - -flux_tally = openmc.Tally(name="flux tally") -flux_tally.filters = [mesh_filter] -flux_tally.scores = ["flux"] - -model.tallies = [flux_tally] - -space = openmc.stats.Point((source_x, source_y, source_z)) -angle = openmc.stats.Isotropic() -energy = openmc.stats.Discrete([2.5e6], [1.0]) - -source = openmc.Source(space=space, angle=angle, energy=energy) -source.particle = "neutron" -model.settings.run_mode = "fixed source" -model.settings.source = source -model.settings.particles = 2000 -model.settings.batches = 5 - - -def run_and_plot(model, filename, output=True): - - sp_filename = model.run(output=output) - - with openmc.StatePoint(sp_filename) as sp: - flux_tally = sp.get_tally(name="flux tally") - - mesh_extent = mesh.bounding_box.extent['xy'] - - # create a plot of the mean flux values - flux_mean = flux_tally.mean.reshape(100, 100) - plt.subplot(1, 2, 1) - plt.imshow( - flux_mean, - origin="lower", - extent=mesh_extent, - norm=LogNorm(), - ) - plt.title("Flux Mean") - - data_slice = geometry.get_slice_of_material_ids(view_direction="z") - xlabel, ylabel = geometry.get_axis_labels(view_direction="z") - plt.xlabel(xlabel) - plt.ylabel(ylabel) - - plt.contour( - np.fliplr(data_slice), - origin="upper", - colors="k", - linestyles="solid", - linewidths=1, - extent=mesh_extent, - ) - - plt.subplot(1, 2, 2) - # create a plot of the flux relative error - flux_std_dev = flux_tally.get_values(value="std_dev").reshape(*mesh.dimension) - plt.imshow( - flux_std_dev, - origin="lower", - extent=mesh_extent, - norm=LogNorm(), - ) - plt.title("Flux Std. Dev.") - - plt.xlabel(xlabel) - plt.ylabel(ylabel) - plt.contour( - np.fliplr(data_slice), - origin="upper", - colors="k", - linestyles="solid", - linewidths=1, - extent=mesh_extent, - ) - plt.savefig(filename) - return sp - - -run_and_plot(model, "no_survival_biasing.png") - -model.settings.survival_biasing = True -model.settings.cutoff = { - "weight": 0.3, # value needs to be between 0 and 1 - "weight_avg": 0.9, # value needs to be between 0 and 1 -} -run_and_plot(model, "yes_survival_biasing.png") diff --git a/tasks/task_13_variance_reduction/spherical_mesh_weight_windows.py b/tasks/task_13_variance_reduction/spherical_mesh_weight_windows.py deleted file mode 100644 index d262ff34..00000000 --- a/tasks/task_13_variance_reduction/spherical_mesh_weight_windows.py +++ /dev/null @@ -1,160 +0,0 @@ -# This example has a sphere of water showing how to increase the depth of -# neutron transport through the water. -# First a simulation with no weight windows is shown -# Weight windows are found from the flux values obtained with the first simulation -# Secondly a simulation with weight windows is performed to find new flux values deeper into the water sphere -# Another set of weight windows is made from the second flux simulation -# The weight window value as a function of depth is plotted to show how these improve with each simulation iteration - -import openmc -import numpy as np -import matplotlib.pyplot as plt - - -# materials -mat_water = openmc.Material() -mat_water.add_element("H", 1) -mat_water.add_element("O", 2) -mat_water.set_density("g/cm3", 1.0) - -my_materials = openmc.Materials([mat_water]) - -# surfaces -outer_surface = openmc.Sphere(r=500, boundary_type="vacuum") - -# regions -region_1 = -outer_surface - -# cells -cell_1 = openmc.Cell(region=region_1) -cell_1.fill = mat_water - -# settings -my_settings = openmc.Settings() - -my_geometry = openmc.Geometry([cell_1]) - -source = openmc.Source() -source.space = openmc.stats.Point((0.0, 0.0, 0.0)) -source.angle = openmc.stats.Isotropic() -source.energy = openmc.stats.Discrete([14e6], [1.0]) -source.particle = "neutron" - -my_settings = openmc.Settings() -my_settings.run_mode = "fixed source" -my_settings.source = source -my_settings.particles = 10000 -my_settings.batches = 10 - -# tally -mesh = openmc.SphericalMesh() -mesh.r_grid = np.linspace(0, outer_surface.r, 100) - -mesh_filter = openmc.MeshFilter(mesh) - -flux_tally = openmc.Tally(name="flux tally") -flux_tally.filters = [mesh_filter] -flux_tally.scores = ["flux"] - -# adds the mesh tally to the model -my_tallies = openmc.Tallies() -my_tallies.append(flux_tally) - -# model -model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies) -output_file = model.run(cwd="no_ww") - -import openmc_weight_window_generator - -# post process the flux values to create the first batch of weight windows -with openmc.StatePoint(output_file) as sp: - flux_tally = sp.get_tally(id=flux_tally.id) - # weight windows from flux - weight_windows = sp.generate_wws(tally=flux_tally, rel_err_tol=0.7)[0] - - -# plot flux against distance -plt.plot( - mesh.r_grid[1:], flux_tally.mean.flatten(), label="flux tally no ww", color="red" -) -plt.legend() -plt.yscale("log") -plt.xlabel("Radius [cm]") -plt.ylabel("Flux") -plt.show() - - -# plot weight window against distance -plt.plot( - mesh.r_grid[1:], - weight_windows.lower_ww_bounds.flatten(), - label="lower ww bounds", - color="red", -) -plt.plot( - mesh.r_grid[1:], - weight_windows.upper_ww_bounds.flatten(), - label="lower up bounds", - color="lightcoral", -) -plt.legend() -plt.xlabel("Radius [cm]") -plt.ylabel("weight window bound value") -plt.show() - -model.settings.weight_windows = weight_windows -# model update by assigning the weight windows found. - -output_file = model.run(cwd="initial_ww") - -with openmc.StatePoint(output_file) as sp: - flux_tally_with_ww = sp.get_tally(id=flux_tally.id) - # weight windows from flux - weight_windows_2 = sp.generate_wws(tally=flux_tally_with_ww, rel_err_tol=0.7)[0] - - -# plot flux with and with out weight windows against distance -plt.plot( - mesh.r_grid[1:], flux_tally.mean.flatten(), label="flux tally no ww", color="red" -) -plt.plot( - mesh.r_grid[1:], - flux_tally_with_ww.mean.flatten(), - label="flux tally with ww", - color="blue", -) -plt.legend() -plt.yscale("log") -plt.xlabel("Radius [cm]") -plt.ylabel("Flux") -plt.show() - -# plot weight window against distance for the firstand second set of weight windows made -plt.plot( - mesh.r_grid[1:], - weight_windows.lower_ww_bounds.flatten(), - label="lower ww bounds", - color="red", -) -plt.plot( - mesh.r_grid[1:], - weight_windows.upper_ww_bounds.flatten(), - label="lower up bounds", - color="lightcoral", -) -plt.plot( - mesh.r_grid[1:], - weight_windows_2.lower_ww_bounds.flatten(), - label="lower ww bounds iteration 2", - color="blue", -) -plt.plot( - mesh.r_grid[1:], - weight_windows_2.upper_ww_bounds.flatten(), - label="lower up bounds iteration 2", - color="cornflowerblue", -) -plt.legend() -plt.xlabel("Radius [cm]") -plt.ylabel("weight window bound value") -plt.show() diff --git a/tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb b/tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb index a981e978..40781722 100644 --- a/tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb +++ b/tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb @@ -60,8 +60,6 @@ "\n", "# MATERIALS\n", "\n", - "mats = openmc.Materials()\n", - "\n", "# makes a simple material from Silver\n", "my_material = openmc.Material() \n", "my_material.add_element('Ag', 1, percent_type='ao')\n", @@ -148,18 +146,16 @@ "# this file tells openmc the decay paths between isotopes including probabilities of different routes and half lives\n", "# To download this xml file you can run these commands\n", "# pip install openmc_data\n", - "# download_nndc_chain -d nuclear_data -r b8.0\n", + "# download_endf_chain -d nuclear_data -r b8.0\n", "\n", "openmc.config['chain_file'] = '/nuclear_data/chain-nndc-b8.0.xml'\n", "\n", "operator = openmc.deplete.CoupledOperator(\n", " model=model,\n", " normalization_mode=\"source-rate\", # set for fixed source simulation, otherwise defaults to fission simulation\n", - " dilute_initial=0, # set to zero to avoid adding small amounts of isotopes, defaults to adding small amounts of fissionable isotopes\n", " reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny\n", " reduce_chain_level=5,\n", - ")\n", - "\n" + ")" ] }, { @@ -241,7 +237,6 @@ "metadata": {}, "outputs": [], "source": [ - "\n", "integrator.integrate()\n", "\n", "# bash command to show the output files produce\n", @@ -324,7 +319,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3.10.6 ('openmc_plot_dev')", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -338,7 +333,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.10" + "version": "3.8.13" }, "vscode": { "interpreter": { diff --git a/tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb b/tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb index b3c05f18..164d6f61 100644 --- a/tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb +++ b/tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb @@ -16,14 +16,13 @@ "outputs": [], "source": [ "# remove any old files\n", - "!rm settings.xm model.xml materials.xml geometry.xml settings.xml\n", + "!rm settings.xm model.xml materials.xml geometry.xml settings.xml model.xml\n", "\n", "import openmc\n", "import openmc.deplete\n", "import math\n", "\n", "\n", - "\n", "# MATERIALS\n", "\n", "mats = openmc.Materials()\n", @@ -128,13 +127,12 @@ "# this file tells openmc the decay paths between isotopes including probabilities of different routes and half lives\n", "# To download this xml file you can run these commands\n", "# pip install openmc_data\n", - "# download_nndc_chain -d nuclear_data -r b8.0\n", - "openmc.config['chain_file'] = '/nuclear_data/chain-nndc-b8.0.xml'\n", + "# download_endf_chain -d nuclear_data -r b8.0\n", + "# openmc.config['chain_file'] = '/nuclear_data/chain-nndc-b8.0.xml'\n", "\n", "operator = openmc.deplete.CoupledOperator(\n", " model=model,\n", " normalization_mode=\"source-rate\", # set for fixed source simulation, otherwise defaults to fission simulation\n", - " dilute_initial=0, # set to zero to avoid adding small amounts of isotopes, defaults to adding small amounts of fissionable isotopes\n", " reduce_chain=True # reduced to only the isotopes present in depletable materials and their possible progeny\n", ")\n", "\n", @@ -186,14 +184,41 @@ "source": [ "results = openmc.deplete.ResultsList.from_hdf5(\"depletion_results.h5\")\n", "\n", - "times, number_of_n_gamma_reactions = results.get_reaction_rate(breeding_material, 'Li6', '(n,gamma)')\n", - "number_of_n_gamma_reactions" + "times, number_of_n_gamma_reactions = results.get_reaction_rate(breeding_material, 'Li6', '(n,gamma)')" ] + }, + { + "cell_type": "markdown", + "id": "f4e5eaee-f68c-48cc-b8c2-fb79ba4abee5", + "metadata": {}, + "source": [ + "Then we can plot the changing reaction rate as a function of time steps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bba28c5c-21fc-4563-a473-d9638541d56a", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "plt.plot(times, number_of_n_gamma_reactions)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a268d5ff-d90a-4e2c-be05-5f8c6d7f61ea", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3.10.6 ('openmc_plot_dev')", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -207,7 +232,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.10" + "version": "3.8.13" }, "vscode": { "interpreter": { diff --git a/tasks/task_14_activation_transmutation_depletion/activation_from_spectra.py b/tasks/task_14_activation_transmutation_depletion/activation_from_spectra.py new file mode 100644 index 00000000..5704cf2e --- /dev/null +++ b/tasks/task_14_activation_transmutation_depletion/activation_from_spectra.py @@ -0,0 +1,99 @@ + +# spectra_tally.mean +# normalised_spectra = spectra_tally.mean / sum(spectra_tally.mean) + +# normalised_spectra is probability of each enegry , sum of normalised_spectra = 1 + + + +import openmc + +# MATERIALS + +# Due to the hydrogen content water is a very good neutron moderator +my_material = openmc.Material() +my_material.add_element('H', 1, percent_type='ao') +my_material.add_element('O', 2, percent_type='ao') +my_material.set_density('g/cm3', 1) + +my_materials = openmc.Materials([my_material]) + + +# GEOMETRY + +# surfaces +outer_surface = openmc.Sphere(r=500, boundary_type='vacuum') + +# cells +cell_1 = openmc.Cell(region=-outer_surface) +cell_1.fill = my_material + +my_geometry = openmc.Geometry([cell_1]) + + +# SIMULATION SETTINGS + +# Instantiate a Settings object +my_settings = openmc.Settings() +my_settings.batches = 10 +my_settings.particles = 10000 +my_settings.run_mode = 'fixed source' + +# Create a DT point source +my_source = openmc.Source() +my_source.space = openmc.stats.Point((0, 0, 0)) +my_source.angle = openmc.stats.Isotropic() +my_source.energy = openmc.stats.Discrete([14e6], [1]) +my_settings.source = my_source + +#creates an empty tally object +my_tallies = openmc.Tallies() + +# sets up filters for the tallies +neutron_particle_filter = openmc.ParticleFilter(['neutron']) + +# creates an array of 300 linear spaced energy bins from 0MeV to 15MeV +# our source is 14MeV so this should capture all the neutron energies in the simulation +# there is a disadvantage of using a linear group structure which is covered in part 2 of this task +import numpy as np +energy_filter = openmc.EnergyFilter(np.linspace(0, 15e6, 300)) + +# setup the filters for the cell tally +cell_filter = openmc.CellFilter(cell_1) + +# create the tally +cell_spectra_tally = openmc.Tally(name='cell_spectra_tally') +cell_spectra_tally.scores = ['flux'] +cell_spectra_tally.filters = [cell_filter, neutron_particle_filter, energy_filter] +my_tallies.append(cell_spectra_tally) + +# combine all the required parts to make a model +model = openmc.model.Model(my_geometry, my_materials, my_settings, my_tallies) + +# remove old files and runs OpenMC + +results_filename = model.run() + +# open the results file +results = openmc.StatePoint(results_filename) + +#extracts the tally values from the simulation results +cell_tally = results.get_tally(name='cell_spectra_tally') + +# flattens the ndarray into a 1d array +flux = cell_tally.mean.flatten() + +spectrum_probability = flux / sum(flux) + +proberbility_per_ev = spectrum_probability / np.diff(energy_filter.bins).flatten() + +tab = openmc.stats.Tabular( + x = energy_filter.bins, + p=proberbility_per_ev, + interpolation='histogram' +) + +energy_filter.tabular(cell_tally.mean.flatten()) + +source = openmc.Source() +source.energy = tab diff --git a/tests/test_task_15.py b/tests/test_task_15.py new file mode 100644 index 00000000..eb284bfc --- /dev/null +++ b/tests/test_task_15.py @@ -0,0 +1,50 @@ + +""" +tests the create_isotope_plot from plotting_utils in the same way the examples +use the function. +""" + +import os +import sys +import unittest +from pathlib import Path + +import nbformat +from nbconvert.preprocessors import ExecutePreprocessor +from nbconvert.preprocessors.execute import CellExecutionError + + +def _notebook_run(path): + """ + Execute a notebook via nbconvert and collect output. + :returns (parsed nb object, execution errors) + """ + kernel_name = 'python%d' % sys.version_info[0] + this_file_directory = os.path.dirname(__file__) + errors = [] + + with open(path) as f: + nb = nbformat.read(f, as_version=4) + nb.metadata.get('kernelspec', {})['name'] = kernel_name + ep = ExecutePreprocessor(kernel_name=kernel_name, timeout=1900) #, allow_errors=True + + try: + ep.preprocess(nb, {'metadata': {'path': this_file_directory}}) + + except CellExecutionError as e: + if "SKIP" in e.traceback: + print(str(e.traceback).split("\n")[-2]) + else: + raise e + + return nb, errors + + +class test_tasks(unittest.TestCase): + +# No module named 'openmc_model' + def test_task_14(self): + for notebook in Path().rglob("tasks/task_15_*/*.ipynb"): + print(notebook) + nb, errors = _notebook_run(notebook) + assert errors == [] diff --git a/tests/test_task_16.py b/tests/test_task_16.py new file mode 100644 index 00000000..7a479631 --- /dev/null +++ b/tests/test_task_16.py @@ -0,0 +1,50 @@ + +""" +tests the create_isotope_plot from plotting_utils in the same way the examples +use the function. +""" + +import os +import sys +import unittest +from pathlib import Path + +import nbformat +from nbconvert.preprocessors import ExecutePreprocessor +from nbconvert.preprocessors.execute import CellExecutionError + + +def _notebook_run(path): + """ + Execute a notebook via nbconvert and collect output. + :returns (parsed nb object, execution errors) + """ + kernel_name = 'python%d' % sys.version_info[0] + this_file_directory = os.path.dirname(__file__) + errors = [] + + with open(path) as f: + nb = nbformat.read(f, as_version=4) + nb.metadata.get('kernelspec', {})['name'] = kernel_name + ep = ExecutePreprocessor(kernel_name=kernel_name, timeout=1900) #, allow_errors=True + + try: + ep.preprocess(nb, {'metadata': {'path': this_file_directory}}) + + except CellExecutionError as e: + if "SKIP" in e.traceback: + print(str(e.traceback).split("\n")[-2]) + else: + raise e + + return nb, errors + + +class test_tasks(unittest.TestCase): + +# No module named 'openmc_model' + def test_task_14(self): + for notebook in Path().rglob("tasks/task_16_*/*.ipynb"): + print(notebook) + nb, errors = _notebook_run(notebook) + assert errors == [] diff --git a/tests/test_task_17.py b/tests/test_task_17.py new file mode 100644 index 00000000..126824eb --- /dev/null +++ b/tests/test_task_17.py @@ -0,0 +1,50 @@ + +""" +tests the create_isotope_plot from plotting_utils in the same way the examples +use the function. +""" + +import os +import sys +import unittest +from pathlib import Path + +import nbformat +from nbconvert.preprocessors import ExecutePreprocessor +from nbconvert.preprocessors.execute import CellExecutionError + + +def _notebook_run(path): + """ + Execute a notebook via nbconvert and collect output. + :returns (parsed nb object, execution errors) + """ + kernel_name = 'python%d' % sys.version_info[0] + this_file_directory = os.path.dirname(__file__) + errors = [] + + with open(path) as f: + nb = nbformat.read(f, as_version=4) + nb.metadata.get('kernelspec', {})['name'] = kernel_name + ep = ExecutePreprocessor(kernel_name=kernel_name, timeout=1900) #, allow_errors=True + + try: + ep.preprocess(nb, {'metadata': {'path': this_file_directory}}) + + except CellExecutionError as e: + if "SKIP" in e.traceback: + print(str(e.traceback).split("\n")[-2]) + else: + raise e + + return nb, errors + + +class test_tasks(unittest.TestCase): + +# No module named 'openmc_model' + def test_task_14(self): + for notebook in Path().rglob("tasks/task_17_*/*.ipynb"): + print(notebook) + nb, errors = _notebook_run(notebook) + assert errors == [] diff --git a/tests/test_task_18.py b/tests/test_task_18.py new file mode 100644 index 00000000..0c608989 --- /dev/null +++ b/tests/test_task_18.py @@ -0,0 +1,50 @@ + +""" +tests the create_isotope_plot from plotting_utils in the same way the examples +use the function. +""" + +import os +import sys +import unittest +from pathlib import Path + +import nbformat +from nbconvert.preprocessors import ExecutePreprocessor +from nbconvert.preprocessors.execute import CellExecutionError + + +def _notebook_run(path): + """ + Execute a notebook via nbconvert and collect output. + :returns (parsed nb object, execution errors) + """ + kernel_name = 'python%d' % sys.version_info[0] + this_file_directory = os.path.dirname(__file__) + errors = [] + + with open(path) as f: + nb = nbformat.read(f, as_version=4) + nb.metadata.get('kernelspec', {})['name'] = kernel_name + ep = ExecutePreprocessor(kernel_name=kernel_name, timeout=1900) #, allow_errors=True + + try: + ep.preprocess(nb, {'metadata': {'path': this_file_directory}}) + + except CellExecutionError as e: + if "SKIP" in e.traceback: + print(str(e.traceback).split("\n")[-2]) + else: + raise e + + return nb, errors + + +class test_tasks(unittest.TestCase): + +# No module named 'openmc_model' + def test_task_14(self): + for notebook in Path().rglob("tasks/task_18_*/*.ipynb"): + print(notebook) + nb, errors = _notebook_run(notebook) + assert errors == []