<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 This notebook is intended to introduce the reader to the depletion interface contained in OpenMC. It is recommended that you are moderately familiar with building models using the OpenMC Python API. The earlier examples are excellent starting points, as this notebook will not focus heavily on model building.
If you have a real power reactor, the fuel composition is constantly changing as fission events produce energy, remove some fissile isotopes, and produce fission products. Other reactions, like $(n, \alpha)$ and $(n, \gamma)$ will alter the composition as well. Furthermore, some nuclides undergo spontaneous decay with widely ranging frequencies. Depletion is the process of modeling this behavior.
In this notebook, we will model a simple fuel pin in an infinite lattice using the Python API. We will then build and examine some of the necessary components for performing depletion analysis. Then, we will use the depletion interface in OpenMC to simulate the fuel pin producing power over several months. Lastly, we will wrap up with some helpful tips to improve the fidelity of depletion simulations.
Step1:  Build the Geometry
Much of this section is borrowed from the "Modeling a Pin-Cell" example. If you find yourself not understanding some aspects of this section, feel free to refer to that example, as some details may be glossed over for brevity.
First, we will create our fuel, cladding, and water materials to represent a typical PWR.
Step2:  Here, we are going to use the openmc.model.pin function to build our pin cell. The pin function anticipates concentric cylinders and materials to fill the inner regions. One additional material is needed than the number of cylinders to cover the domain outside the final ring. 
To do this, we define two radii for the outer radius of our fuel pin, and the outer radius of the cladding.
Step3:  Using these radii, we define concentric ZCylinder objects. So long as the cylinders are concentric and increasing in radius, any orientation can be used. We also take advantage of the fact that the openmc.Materials object is a subclass of the list object to assign materials to the regions defined by the surfaces.
Step4:  The first material, in our case fuel, is placed inside the first cylinder in the inner-most region. The second material, clad, fills the space between our cylinders, while water is placed outside the last ring. The pin function returns an openmc.Universe object, and has some additional features we will mention later.
Step5:  Lastly we construct our settings. For the sake of time, a relatively low number of particles will be used.
Step6:  The depletion interface relies on OpenMC to perform the transport simulation and obtain reaction rates and other important information. We then have to create the xml input files that openmc expects, specifically geometry.xml, settings.xml, and materials.xml.
Step7:  Before we write the material file, we must add one bit of information
Step8:  Setting up for depletion
The OpenMC depletion interface can be accessed from the openmc.deplete module, and has a variety of classes that will help us.
Step9:  In order to run the depletion calculation we need the following information
Step10:  The primary entry point for depletion is the openmc.deplete.Operator. It relies on the openmc.deplete.Chain and helper classes to run openmc, retrieve and normalize reaction rates, and other perform other tasks. For a thorough description, please see the full API documentation.
We will create our Operator using the geometry and settings from above, and our simple chain file. The materials are read in automatically using the materials.xml file.
Step11:  We will then simulate our fuel pin operating at linear power of 174 W/cm, or 174 W given a unit height for our problem.
Step12:  For this problem, we will take depletion step sizes of 30 days, and instruct OpenMC to re-run a transport simulation every 30 days until we have modeled the problem over a six month cycle. The depletion interface expects the time to be given in seconds, so we will have to convert. Note that these values are not cumulative.
Step13:  And lastly, we will use the basic predictor, or forward Euler, time integration scheme. Other, more advanced methods are provided to the user through openmc.deplete
Step14:  To perform the simulation, we use the integrate method, and let openmc take care of the rest.
Step15:  Processing the outputs
The depletion simulation produces a few output files. First, the statepoint files from each individual transport simulation are written to openmc_simulation_n&lt;N&gt;.h5, where &lt;N&gt; indicates the current depletion step. Any tallies that we defined in tallies.xml will be included in these files across our simulations. We have 7 such files, one for each our of 6 depletion steps and the initial state.
Step16:  The depletion_results.h5 file contains information that is aggregated over all time steps through depletion. This includes the multiplication factor, as well as concentrations. We can process this file using the openmc.deplete.ResultsList object
Step17:  The first column of k is the value of k-combined at each point in our simulation, while the second column contains the associated uncertainty. We can plot this using matplotlib
Step18:  Due to the low number of particles selected, we have not only a very high uncertainty, but likely a horrendously poor fission source. This pin cell should have $k>1$, but we can still see the decline over time due to fuel consumption.
We can then examine concentrations of atoms in each of our materials. This requires knowing the material ID, which can be obtained from the materials.xml file.
Step19:  We can also examine reaction rates over time using the ResultsList
Step20:  Helpful tips
Depletion is a tricky task to get correct. Use too short of time steps and you may never get your results due to running many transport simulations. Use long of time steps and you may get incorrect answers. Consider the xenon plot from above. Xenon-135 is a fission product with a thermal absorption cross section on the order of millions of barns, but has a half life of ~9 hours. Taking smaller time steps at the beginning of your simulation to build up some equilibrium in your fission products is highly recommended.
When possible, differentiate materials that reappear in multiple places. If we had built an entire core with the single fuel material, every pin would be depleted using the same averaged spectrum and reaction rates which is incorrect. The Operator can differentiate these materials using the diff_burnable_mats argument, but note that the volumes will be copied from the original material.
Using higher-order integrators, like the CECMIntegrator, EPCRK4Integrator with a fourth order Runge-Kutta, or the LEQIIntegrator, can improve the accuracy of a simulation, or at least allow you to take longer depletion steps between transport simulations with similar accuracy.
Fuel pins with integrated burnable absorbers, like gadolinia, experience strong flux gradients until the absorbers are mostly burned away. This means that the spectrum and magnitude of the flux at the edge of the fuel pin can be vastly different than that in the interior. The helper pin function can be used to subdivide regions into equal volume segments, as follows.
Step21:  The innermost region has been divided into 10 equal volume regions. We can pass additional arguments to divide multiple regions, except for the region outside the last cylinder.
Register depletion chain
The depletion chain we created can be registered into the OpenMC cross_sections.xml file, so we don't have to always pass the chain_file argument to the Operator. To do this, we create a DataLibrary using openmc.data. Without any arguments, the from_xml method will look for the file located at OPENMC_CROSS_SECTIONS. For this example, we will just create a bare library.
Step22:  This allows us to make an Operator simply with the geometry and settings arguments, provided we exported our library to OPENMC_CROSS_SECTIONS. For a problem where we built and registered a Chain using all the available nuclear data, we might see something like the following.
Step23:  Choice of depletion step size
A general rule of thumb is to use depletion step sizes around 2 MWd/kgHM, where kgHM is really the initial heavy metal mass in kg. If your problem includes integral burnable absorbers, these typically require shorter time steps at or below 1 MWd/kgHM. These are typically valid for the predictor scheme, as the point of recent schemes is to extend this step size. A good convergence study, where the step size is decreased until some convergence metric is satisfied, is a beneficial exercise.
We can use the Operator to determine our maximum step size using this recommendation. The heavy_metal attribute returns the mass of initial heavy metal in g, which, using our power, can be used to compute this step size. $$\frac{2\,MWd}{kgHM} = \frac{P\times\Delta}{hm_{op}}$$
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import math
import openmc
Explanation: This notebook is intended to introduce the reader to the depletion interface contained in OpenMC. It is recommended that you are moderately familiar with building models using the OpenMC Python API. The earlier examples are excellent starting points, as this notebook will not focus heavily on model building.
If you have a real power reactor, the fuel composition is constantly changing as fission events produce energy, remove some fissile isotopes, and produce fission products. Other reactions, like $(n, \alpha)$ and $(n, \gamma)$ will alter the composition as well. Furthermore, some nuclides undergo spontaneous decay with widely ranging frequencies. Depletion is the process of modeling this behavior.
In this notebook, we will model a simple fuel pin in an infinite lattice using the Python API. We will then build and examine some of the necessary components for performing depletion analysis. Then, we will use the depletion interface in OpenMC to simulate the fuel pin producing power over several months. Lastly, we will wrap up with some helpful tips to improve the fidelity of depletion simulations.
End of explanation
fuel = openmc.Material(name="uo2")
fuel.add_element("U", 1, percent_type="ao", enrichment=4.25)
fuel.add_element("O", 2)
fuel.set_density("g/cc", 10.4)
clad = openmc.Material(name="clad")
clad.add_element("Zr", 1)
clad.set_density("g/cc", 6)
water = openmc.Material(name="water")
water.add_element("O", 1)
water.add_element("H", 2)
water.set_density("g/cc", 1.0)
water.add_s_alpha_beta("c_H_in_H2O")
materials = openmc.Materials([fuel, clad, water])
Explanation: Build the Geometry
Much of this section is borrowed from the "Modeling a Pin-Cell" example. If you find yourself not understanding some aspects of this section, feel free to refer to that example, as some details may be glossed over for brevity.
First, we will create our fuel, cladding, and water materials to represent a typical PWR.
End of explanation
radii = [0.42, 0.45]
Explanation: Here, we are going to use the openmc.model.pin function to build our pin cell. The pin function anticipates concentric cylinders and materials to fill the inner regions. One additional material is needed than the number of cylinders to cover the domain outside the final ring. 
To do this, we define two radii for the outer radius of our fuel pin, and the outer radius of the cladding.
End of explanation
pin_surfaces = [openmc.ZCylinder(r=r) for r in radii]
pin_univ = openmc.model.pin(pin_surfaces, materials)
Explanation: Using these radii, we define concentric ZCylinder objects. So long as the cylinders are concentric and increasing in radius, any orientation can be used. We also take advantage of the fact that the openmc.Materials object is a subclass of the list object to assign materials to the regions defined by the surfaces.
End of explanation
pin_univ.plot()
bound_box = openmc.rectangular_prism(0.62, 0.62, boundary_type="reflective")
root_cell = openmc.Cell(fill=pin_univ, region=bound_box)
root_univ = openmc.Universe(cells=[root_cell])
geometry = openmc.Geometry(root_univ)
Explanation: The first material, in our case fuel, is placed inside the first cylinder in the inner-most region. The second material, clad, fills the space between our cylinders, while water is placed outside the last ring. The pin function returns an openmc.Universe object, and has some additional features we will mention later.
End of explanation
settings = openmc.Settings()
settings.particles = 100
settings.inactive = 10
settings.batches = 50
Explanation: Lastly we construct our settings. For the sake of time, a relatively low number of particles will be used.
End of explanation
geometry.export_to_xml()
settings.export_to_xml()
Explanation: The depletion interface relies on OpenMC to perform the transport simulation and obtain reaction rates and other important information. We then have to create the xml input files that openmc expects, specifically geometry.xml, settings.xml, and materials.xml.
End of explanation
fuel.volume = math.pi * radii[0] ** 2
materials.export_to_xml()
Explanation: Before we write the material file, we must add one bit of information: the volume of our fuel. In order to translate the reaction rates obtained by openmc to meaningful units for depletion, we have to normalize them to a correct power. This requires us to know, or be able to calculate, how much fuel is in our problem. Correctly setting the volumes is a critical step, and can lead to incorrect answers, as the fuel is over- or under-depleted due to poor normalization.
For our problem, we can assign the "volume" to be the cross-sectional area of our fuel. This is identical to modeling our fuel pin inside a box with height of 1 cm.
End of explanation
import openmc.deplete
Explanation: Setting up for depletion
The OpenMC depletion interface can be accessed from the openmc.deplete module, and has a variety of classes that will help us.
End of explanation
chain = openmc.deplete.Chain.from_xml("./chain_simple.xml")
chain.nuclide_dict
Explanation: In order to run the depletion calculation we need the following information:
Nuclide decay, fission yield, and reaction data
Operational power or power density
Desired depletion schedule
Desired time integration scheme
The first item is necessary to determine the paths by which nuclides transmute over the depletion simulation. This includes spontaneous decay, fission product yield distributions, and nuclides produced through neutron-reactions. For example,
* Te129 decays to I129 with a half life of ~70 minutes
* A fission event for U-235 produces fission products like Xe135 according to a distribution
* For thermal problems, Am241 will produce metastable Am242 about 8% of the time during an $(n,\gamma)$ reaction. The other 92% of capture reactions will produce ground state Am242
These data are often distributed with other nuclear data, like incident neutron cross sections with ENDF/B-VII.
OpenMC uses the openmc.deplete.Chain to collect represent the various decay and transmutation pathways in a single object.
While a complete Chain can be created using nuclear data files, users may prefer to download pre-generated XML-representations instead.
Such files can be found at https://openmc.org/depletion-chains/ and include full and compressed chains, with capture branching ratios derived using PWR- or SFR-spectra.
For this problem, we will be using a much smaller depletion chain that contains very few nuclides. In a realistic problem, over 1000 isotopes may be included in the depletion chain.
End of explanation
operator = openmc.deplete.Operator(geometry, settings, "./chain_simple.xml")
Explanation: The primary entry point for depletion is the openmc.deplete.Operator. It relies on the openmc.deplete.Chain and helper classes to run openmc, retrieve and normalize reaction rates, and other perform other tasks. For a thorough description, please see the full API documentation.
We will create our Operator using the geometry and settings from above, and our simple chain file. The materials are read in automatically using the materials.xml file.
End of explanation
power = 174
Explanation: We will then simulate our fuel pin operating at linear power of 174 W/cm, or 174 W given a unit height for our problem.
End of explanation
time_steps = [30 * 24 * 60 * 60] * 6
Explanation: For this problem, we will take depletion step sizes of 30 days, and instruct OpenMC to re-run a transport simulation every 30 days until we have modeled the problem over a six month cycle. The depletion interface expects the time to be given in seconds, so we will have to convert. Note that these values are not cumulative.
End of explanation
integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, power)
Explanation: And lastly, we will use the basic predictor, or forward Euler, time integration scheme. Other, more advanced methods are provided to the user through openmc.deplete
End of explanation
integrator.integrate()
Explanation: To perform the simulation, we use the integrate method, and let openmc take care of the rest.
End of explanation
!ls *.h5
Explanation: Processing the outputs
The depletion simulation produces a few output files. First, the statepoint files from each individual transport simulation are written to openmc_simulation_n&lt;N&gt;.h5, where &lt;N&gt; indicates the current depletion step. Any tallies that we defined in tallies.xml will be included in these files across our simulations. We have 7 such files, one for each our of 6 depletion steps and the initial state.
End of explanation
results = openmc.deplete.ResultsList.from_hdf5("./depletion_results.h5")
time, k = results.get_eigenvalue()
time /= (24 * 60 * 60)  # convert back to days from seconds
k
Explanation: The depletion_results.h5 file contains information that is aggregated over all time steps through depletion. This includes the multiplication factor, as well as concentrations. We can process this file using the openmc.deplete.ResultsList object
End of explanation
from matplotlib import pyplot
pyplot.errorbar(time, k[:, 0], yerr=k[:, 1])
pyplot.xlabel("Time [d]")
pyplot.ylabel("$k_{eff}\pm \sigma$");
Explanation: The first column of k is the value of k-combined at each point in our simulation, while the second column contains the associated uncertainty. We can plot this using matplotlib
End of explanation
_time, u5 = results.get_atoms("1", "U235")
_time, xe135 = results.get_atoms("1", "Xe135")
pyplot.plot(time, u5, label="U235")
pyplot.xlabel("Time [d]")
pyplot.ylabel("Number of atoms - U235");
pyplot.plot(time, xe135, label="Xe135")
pyplot.xlabel("Time [d]")
pyplot.ylabel("Number of atoms - Xe135");
Explanation: Due to the low number of particles selected, we have not only a very high uncertainty, but likely a horrendously poor fission source. This pin cell should have $k>1$, but we can still see the decline over time due to fuel consumption.
We can then examine concentrations of atoms in each of our materials. This requires knowing the material ID, which can be obtained from the materials.xml file.
End of explanation
_time, u5_fission = results.get_reaction_rate("1", "U235", "fission")
pyplot.plot(time, u5_fission)
pyplot.xlabel("Time [d]")
pyplot.ylabel("Fission reactions / s");
Explanation: We can also examine reaction rates over time using the ResultsList
End of explanation
div_surfs_1 = [openmc.ZCylinder(r=1)]
div_1 = openmc.model.pin(div_surfs_1, [fuel, water], subdivisions={0: 10})
div_1.plot(width=(2.0, 2.0))
Explanation: Helpful tips
Depletion is a tricky task to get correct. Use too short of time steps and you may never get your results due to running many transport simulations. Use long of time steps and you may get incorrect answers. Consider the xenon plot from above. Xenon-135 is a fission product with a thermal absorption cross section on the order of millions of barns, but has a half life of ~9 hours. Taking smaller time steps at the beginning of your simulation to build up some equilibrium in your fission products is highly recommended.
When possible, differentiate materials that reappear in multiple places. If we had built an entire core with the single fuel material, every pin would be depleted using the same averaged spectrum and reaction rates which is incorrect. The Operator can differentiate these materials using the diff_burnable_mats argument, but note that the volumes will be copied from the original material.
Using higher-order integrators, like the CECMIntegrator, EPCRK4Integrator with a fourth order Runge-Kutta, or the LEQIIntegrator, can improve the accuracy of a simulation, or at least allow you to take longer depletion steps between transport simulations with similar accuracy.
Fuel pins with integrated burnable absorbers, like gadolinia, experience strong flux gradients until the absorbers are mostly burned away. This means that the spectrum and magnitude of the flux at the edge of the fuel pin can be vastly different than that in the interior. The helper pin function can be used to subdivide regions into equal volume segments, as follows.
End of explanation
data_lib = openmc.data.DataLibrary()
data_lib.register_file("./chain_simple.xml")
data_lib.export_to_xml()
!cat cross_sections.xml
Explanation: The innermost region has been divided into 10 equal volume regions. We can pass additional arguments to divide multiple regions, except for the region outside the last cylinder.
Register depletion chain
The depletion chain we created can be registered into the OpenMC cross_sections.xml file, so we don't have to always pass the chain_file argument to the Operator. To do this, we create a DataLibrary using openmc.data. Without any arguments, the from_xml method will look for the file located at OPENMC_CROSS_SECTIONS. For this example, we will just create a bare library.
End of explanation
new_op = openmc.deplete.Operator(geometry, settings)
len(new_op.chain.nuclide_dict)
[nuc.name for nuc in new_op.chain.nuclides[:10]]
[nuc.name for nuc in new_op.chain.nuclides[-10:]]
Explanation: This allows us to make an Operator simply with the geometry and settings arguments, provided we exported our library to OPENMC_CROSS_SECTIONS. For a problem where we built and registered a Chain using all the available nuclear data, we might see something like the following.
End of explanation
operator.heavy_metal
max_step = 2 * operator.heavy_metal / power * 1E3
print("\"Maximum\" depletion step: {:5.3} [d]".format(max_step))
Explanation: Choice of depletion step size
A general rule of thumb is to use depletion step sizes around 2 MWd/kgHM, where kgHM is really the initial heavy metal mass in kg. If your problem includes integral burnable absorbers, these typically require shorter time steps at or below 1 MWd/kgHM. These are typically valid for the predictor scheme, as the point of recent schemes is to extend this step size. A good convergence study, where the step size is decreased until some convergence metric is satisfied, is a beneficial exercise.
We can use the Operator to determine our maximum step size using this recommendation. The heavy_metal attribute returns the mass of initial heavy metal in g, which, using our power, can be used to compute this step size. $$\frac{2\,MWd}{kgHM} = \frac{P\times\Delta}{hm_{op}}$$
End of explanation
<END_TASK>