This notebook was created for the THE GRE@T-PIONEeR course “Laboratory Experiments at the Training Reactor of BME”.
The purpose of this notebook is to verify the proper installation and usability of OpenMC and the python API in Jupyter Lab environment as well as to introduce the syntax of creating some objects for a geometry through a very simple test case. This excercise is not at all intended as a comprehensive guide or tutorial, however it can provide some basic experience regarding the structure of an OpenMC input script.
TASK: Execute all blocks of this notebook: build a geometry by solving tasks 1 (material definition) and 2 (geometry definition), then run the simulation. Finally, make sure to complete the related Matlab Grader assignment
Tip: Check out the available commads in the Run menu.
Tip: To get familiar with the use of the OpenMC through the python API, it is recommended to examine the contents of each block and execute them one by one.
Tip: Note that a notebook allows non-consecutive execution of the notebook-cells. This is a nice feature, but it can confuse the workspace if parts of intertwined objects are changed or re-created with different identification numbers. It is recommended to re-run the notebook from the beginning if you experience strange behaviour.
STEP 0 - Initialize necessary python modules
numpy for arrays and some other stuff
import numpy as np
matplotlib.pyplot for plotting
from matplotlib import pyplot as plt
OpenMC python API
import openmc
This python file contains fuctions to check proper object definitions in the tasks.
from inspector import *
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 11
8 import openmc
10 # This python file contains fuctions to check proper object definitions in the tasks. —> 11 from inspector import *
ModuleNotFoundError: No module named ‘inspector’
STEP 1 - Defining materials for the geometry
Different types of materials can be created with the openmc.Material() class.
The material objects contain information about e.g.:
composition: .add_nuclide() and .add_element() ,
density: set_density() ,
S(a,b) treatment: add_s_alpha_beta() ,
…
The materials use in the model are bundled together in the openmc.Materials() object This will be converted into an .xml file for the openmc application.
Some of the properties/methods used for defining the material object can be seen below.
No description has been provided for this image No description has been provided for this image No description has been provided for this image
(Source: openmc.Material; accessed on 10/12/2023)
materials = openmc.Materials()
water
MAT_water=openmc.Material(material_id=1,name=“water”) MAT_water.add_nuclide(‘H1’, 0.6672, ‘ao’) # adding isotope MAT_water.add_nuclide(‘O16’, 0.3328, ‘ao’)
MAT_water.set_density(‘g/cm3’, 9.978e-01)
MAT_water.add_s_alpha_beta(‘c_H_in_H2O’)
materials.append(MAT_water)
U fuel
MAT_fuel=openmc.Material(material_id=2,name=“fuel”)
MAT_fuel.add_nuclide(‘U235’,0.25, ‘wo’)
MAT_fuel.add_nuclide(‘U238’,0.75, ‘wo’)
MAT_fuel.set_density(‘g/cm3’, 19.1)
materials.append(MAT_fuel)
graphite
MAT_graphite=openmc.Material(material_id=3,name=“graphite”) MAT_graphite.add_element(‘C’, 1, ‘ao’) # adding element MAT_graphite.set_density(‘g/cm3’, 1.7)
MAT_graphite.add_s_alpha_beta(‘c_Graphite’) materials.append(MAT_graphite)
TASK 1: Create a material object with the following properties!
Variable name
Set the variable name of the object in the python api to ‘MAT_U235’.
Identification
Set the “material_id” of the object to 42.
Composition
The composition of the material sould be pure U235 isotope. (In such a case, the fraction type ‘ao’ or ‘wo’ does not matter.)
Density
Set the density of the material to 19 g/cm^3
Make sure to include the material in the materials list!
Use materials.append() method.
Your solution
import openmc
Initialize the materials collection
materials = openmc.Materials()
Water
MAT_water = openmc.Material(name=“water”)
MAT_water.add_nuclide(‘H1’, 0.6672, ‘ao’)
MAT_water.add_nuclide(‘O16’, 0.3328, ‘ao’)
MAT_water.set_density(‘g/cm3’, 0.9978)
MAT_water.add_s_alpha_beta(‘c_H_in_H2O’)
materials.append(MAT_water)
U fuel
MAT_fuel = openmc.Material(name=“fuel”)
MAT_fuel.add_nuclide(‘U235’, 0.25, ‘wo’)
MAT_fuel.add_nuclide(‘U238’, 0.75, ‘wo’)
MAT_fuel.set_density(‘g/cm3’, 19.1)
materials.append(MAT_fuel)
Graphite
MAT_graphite = openmc.Material(name=“graphite”)
MAT_graphite.add_element(‘C’, 1.0, ‘ao’)
MAT_graphite.set_density(‘g/cm3’, 1.7)
MAT_graphite.add_s_alpha_beta(‘c_Graphite’)
materials.append(MAT_graphite)
MAT_U235 (Task 1 specific)
MAT_U235 = openmc.Material(material_id=42, name=‘MAT_U235’) MAT_U235.add_nuclide(‘U235’, 1.0, ‘ao’) # Pure U-235
MAT_U235.set_density(‘g/cm3’, 19) materials.append(MAT_U235)
Export materials to XML
materials.export_to_xml()
Export materials to xml format file which will be used by openmc
materials.export_to_xml()
DO NOT MODIFY THIS BLOCK
check_material(MAT_U235) check_materials(materials)
Tip: Examine ‘materials.xml’ in the file browser.
At this point, use the “Kernel > Restart Kernel and Run up to Selected Cell…” option to clean up the workspace of the notebook!
If at any point you enconter errors during running a notebook, this may be one of the top solutions, because due to repeated and out-of-sequence cell-executions, the workspace can get messy. Another, more efficient cleanup can be the removal of generated xml, out and h5 files as well.
STEP 2 - Defining the geometry
In OpenMC, complex arbitrary 3D geometrical regions are defined using Constructive Solid Geometry (CSG) .
The geometry of a region is defined along the following steps:
simple surfaces (planes, spheres, etc.) are defined,
then half-spaces are defined using these surfaces,
then complex regions are defined by doing boolean operations with the half-spaces
The different regions are described by Cell objects. Besides the region, these objects may store information about more complex sub-geometry (eg. in the form of a Universe, see below) and/or material properties. Cells are bundled together into Universe objects. As mentioned previously, universes can be filled into cells.
Our geometry consists of the following regions:
pure U-235 cylinder in the center
U235-U238 mixture shell
graphite reflector shell
water outside
No description has been provided for this image
The final structure of the geometry needs to be a single universe.
root=openmc.Universe()
TASK 2a: Create surfaces for the cylinder core with the following properties!
Variable names
Do not modify the provided variable names.
Side
Create a cylinder on the Z-axis (openmc.ZCylinder) with 4 cm diameter.
Top and bottom
Create planes perpendicular to Z-axis (openmc.ZPlane) at 3 cm level for the top, and -3 cm level for the bottom.
The list of surfaces: openmc manual
Some of the surfaces can be seen below.
No description has been provided for this image (Source: openmc.ZPlane; accessed on 16/10/2023)
No description has been provided for this image (Source: openmc.ZCylinder; accessed on 16/10/2023)
Your solution
Cylinder and Planes for Core
cylinder_diameter = 4.0 # cm
core_side = openmc.ZCylinder(r=cylinder_diameter / 2)
core_top = openmc.ZPlane(z0=3.0, boundary_type=‘transmission’) core_bottom = openmc.ZPlane(z0=-3.0, boundary_type=‘transmission’)
DO NOT MODIFY THIS BLOCK
check_core_side(core_side) check_core_top(core_top) check_core_bottom(core_bottom)
Spatial regions (half-spaces) are defined with CSG in the following manner:
Divide space with a surface and select one ‘side’ according to the normal vector of the surface:
The normal vector for x-, y-, z-planes looks in the positive direction on the axis, The normal vector for cylinders and spheres looks in the outwards direction.
Use + to indicate the spatial region along the normal vector
Use - to indicate the spatial region in the opposite direction to the normal vector
Use boolean logic to merge, cut and intersect spatial regions
Use & symbol between half-spaces ond surfaces Use ( and ) to prioritize operations
Use ~ to get the inverse of complex regions
TASK 2b: Create region (half-space) for the cylinder core with the following properties!
Variable name
Do not modify the provided variable name.
Region
Define the region inside the cylinder and between the bottom and top planes.
Tip: You can take a look at the region definitions in other cells after Task 2c for reference.
Your solution
Core Region
core_region = +core_bottom & -core_top & -core_side
DO NOT MODIFY THIS BLOCK
check_core_region(core_side,core_top,core_bottom,core_region)
TASK 2c: Create cell object for the cylinder core with the following properties!
Variable name
Do not modify the provided variable name.
Region
Set the region of the cell object to ‘core_region’.
Identification
Set the ‘cell_id’ parameter of the cell object to 42.
Material
Fill the cell object with MAT_U235 material.
Add to geometry
Add the cell object to the root universe with root.add_cell() method. No description has been provided for this image
(Source: openmc.Cell; accessed on 10/12/2023)
Tip: You can take a look at the definitions in other cells below for reference.
Your solution
Core Cell
cell_core = openmc.Cell(name=‘cylinder_core’, fill=MAT_U235, region=core_region)
DO NOT MODIFY THIS BLOCK
check_core_cell(core_region,MAT_U235,cell_core) check_root(cell_core,root)
U sphere
sphere_1_radius=12
s_sphere_1=openmc.Sphere(r=sphere_1_radius)
cell_fuel=openmc.Cell(region=(-s_sphere_1 & ~core_region),cell_id=1,name=“U_sphere”)
cell_fuel.fill=MAT_fuel
root.add_cell(cell_fuel)
reflector shell
sphere_2_radius=22
s_sphere_2=openmc.Sphere(r=sphere_2_radius)
cell_reflector=openmc.Cell(region=(-s_sphere_2 & +s_sphere_1),cell_id=2,name=“reflector_shell”)
cell_reflector.fill=MAT_graphite
root.add_cell(cell_reflector)
Outer boundary of the problem
s_boundary_x1=openmc.XPlane(x0=-40,boundary_type=‘vacuum’)
s_boundary_x2=openmc.XPlane(x0= 40,boundary_type=‘vacuum’)
s_boundary_y1=openmc.YPlane(y0=-40,boundary_type=‘vacuum’)
s_boundary_y2=openmc.YPlane(y0= 40,boundary_type=‘vacuum’)
s_boundary_z1=openmc.ZPlane(z0=-40,boundary_type=‘vacuum’)
s_boundary_z2=openmc.ZPlane(z0= 40,boundary_type=‘vacuum’)
region_boundary = +s_boundary_x1 & -s_boundary_x2 & +s_boundary_y1 & -s_boundary_y2 & +s_boundary_z1 & -s_boundary_z2
some water
cell_workspace=openmc.Cell(region=(+s_sphere_2 & region_boundary ),cell_id=3,name=“some_water”) cell_workspace.fill=MAT_water
root.add_cell(cell_workspace)
The geometry needs to be exported into a xml file.
geometry = openmc.Geometry(root)
geometry.export_to_xml()
Tip: Examine ‘geometry.xml’ in the file browser.
The geometry can be visualized for example using the openmc.Universe.plot() member, as seen below:
root.plot(
basis = ‘xz’,
width = (100, 100),
pixels = (500, 500),
color_by = ‘material’,
colors = {MAT_water: ‘powderblue’,MAT_fuel: ‘tomato’,MAT_graphite: ‘lightslategray’,MAT_U235: ‘purple’} )
TASK 2 geometry inspection: Compare the plot above to the expected geometry!
Due to that automatic checking is not available for the direction parameter in the surface-to-halfspace conversion, you need to manually make sure that the cylindrical core element in the geometry is properly defined.
Compare to two plots and make sure that all four regions are properly present!
In case of a difference, first try to re-run the notebook, otherwise, check the region definition!
Expected geometry:
No description has been provided for this image
At this point, use the “Kernel > Restart Kernel and Run up to Selected Cell…” option to clean up the workspace of the notebook!
STEP 3 - Defining a simple source
We need a set of source points at the start of the criticality calculation.
In case of a k_eff calculation, we need to have at least one (1) point in fissile medium, but more points may help the source distribution to converge faster.
Isotropic point source
source_n=openmc.Source(space=openmc.stats.Point((0,0,0)))
source_n.particle=‘neutron’ # optional, as ‘neutron’ is the default
STEP 4 - Setting parameters of the simulation
settings = openmc.Settings()
Criticality calculation
settings.batches = 1200
settings.inactive = 200
settings.particles = 5000
Source
settings.source = source_n
Exporting the settings into a xml file.
settings.export_to_xml()
Tip: Examine ‘settings.xml’ in the file browser.
STEP 5 - Defining tallies (“optional”)
Generally speaking, a Monte Carlo code is more than happy to run a simulation without calculating and outputting any quantities.
In case of OpenMC, a series if filters can be defined (energy, space, angle, isotope, etc.), through which it will collect contributions to the desired quantity.
Without getting into the detailes of each object, let’s define a tally to calculate the spatial distribution of the thermal neutron flux along y=0 plane
Defining a spatial filter
x_min=-40
x_max=40
y_min=x_min
y_max=x_max
resolution=160
myMesh_xz=openmc.RegularMesh()
myMesh_xz.dimension=[resolution,1,resolution]
myMesh_xz.lower_left=[x_min,-2,y_min]
myMesh_xz.upper_right=[x_max,2,y_max]
meshfilter_xz=openmc.MeshFilter(myMesh_xz)
…and an energy filter for thermal energies
energybins_th=np.array([0., 0.625])
energy_filter_th = openmc.EnergyFilter(values=energybins_th)
Defining the tally
tally_xz=openmc.Tally(name=‘thermal_flux’) tally_xz.filters = [meshfilter_xz,energy_filter_th] tally_xz.scores = [‘flux’]
Exporting the tally information into a xml file.
tallies=openmc.Tallies()
tallies.append(tally_xz) tallies.export_to_xml()
Tip: Examine ‘tallies.xml’ in the file browser.
STEP 6 - Run the simulation
An OpenMC simulation can be initiated in a number of ways, for example by entering the openmc command in a system terminal, or by using the openmc.run() member inside the API. Note that the program will use the data defined in the ‘.xml’ files, and it is independent from the Jupyter Lab or python environment. It is therefore important to make sure that any desired changes are registered in the ‘.xml’ files.
option=2
if option==1 :
!openmc else:
openmc.run()
STEP 7 - Post processing
Some information is displayed during simulation as well as at the end of the calculation regarding the criticality calculation, but everything else needs to be extracted from the output files. We are going to use the statepoint.*.h5 binary file, which contains all results, but the tally outputs can be found also in the ‘tallies.out’ file.
Tip: Examine the contents ‘tallies.out’ file in the file browser.
Extracting data…
with openmc.StatePoint(‘statepoint.’+str(settings.batches)+‘.h5’) as sp:
tally_flux_xz = sp.get_tally(name=‘thermal_flux’)
keff=sp.keff
print(“total run time :”,int(np.round(sp.runtime[‘total’])),‘s’)
sp.close()
print(‘k_eff =’,np.round(keff.n,5),’ +/- ',np.round(keff.s,5))
Processing data…
dataframe_xz=tally_flux_xz.get_pandas_dataframe()
tally_mean_values_xz=dataframe_xz[‘mean’].values.reshape((resolution,resolution))
Displaying data…
plt.imshow(tally_mean_values_xz, interpolation=‘nearest’, extent=[x_min,x_max,y_max,y_min]) # NOTE: ‘extent’ parameter is “visual”, it doesn’t affect the data range plt.gca().invert_yaxis() # tally data order is different than drawing direction
plt.title(‘thermal flux (y=0 cm)’)
plt.set_cmap(‘gist_ncar’) # prism, gist_ncar plt.xlabel(‘x (cm)’)
plt.ylabel(‘z (cm)’) plt.colorbar()
Reminder: Make sure to complete the associated Matlab Grader assignment.
Tip: After completing the associated Matlab Grader assignment, you are free to experiment with the notebook.