{
"metadata": {
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Neutron Diffusion in Python \n",
"\n",
"This notebook is an entirely self-contained solution to a basic [neutron diffision](http://mragheb.com/NPRE%20402%20ME%20405%20Nuclear%20Power%20Engineering/One%20Group%20Reactor%20Theory.pdf) equation for a reactor *rx* made up of a single fuel rod. The one-group diffusion equation that we will be stepping through time and space is, \n",
"\n",
"$\\frac{1}{v}\\frac{\\partial \\phi}{\\partial t} = D \\nabla^2 \\phi + (k - 1) \\Sigma_a \\phi + S$\n",
"\n",
"where \n",
"\n",
"* $\\phi$ is the neutron flux [n/cm$^2$/s],\n",
"* $D$ is the diffusion coefficient [cm],\n",
"* $k$ is the multiplication factor of the material [unitless],\n",
"* $S$ is a static source term [n/cm$^2$/s], and\n",
"* $v$ is the neutron velocity, which for [thermal neutrons](http://en.wikipedia.org/wiki/Neutron_temperature) is 2.2e5 [cm/s]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## How Does This Work?!\n",
"\n",
"Here we use [ITAPS / MOAB](https://trac.mcs.anl.gov/projects/ITAPS/wiki/MOAB), via [PyTAPS](http://pythonhosted.org/PyTAPS/), as our underlying mesh representation. This allows us to keep a reperesentation of the mesh in-memory or dump it out to disk. The reason for using MOAB, other than that it saves us a ton of time and effort since we don't have to write our own mesh package, is that MOAB is now understood by both [PyNE](http://pynesim.org/) and [yt](http://yt-project.org/)!\n",
"\n",
"PyNE is used to compute all of the nuclear data needs here and for a structured / semi-structured represntation of MOAB Hex8 meshes. yt has been supped up to understand both plain old MOAB Hex8 meshes as well as the PyNE variant. The simulation, analysis, and visulaization here takes place entirely within memory.\n",
"\n",
"The software stack here is gloriously deep. Refresh yourself on this well of knowledge."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from itertools import product \n",
"from pyne.mesh import Mesh\n",
"from pyne.xs.cache import XSCache\n",
"from pyne.xs.data_source import CinderDataSource, SimpleDataSource, NullDataSource\n",
"from pyne.xs.channels import sigma_a, sigma_s\n",
"from pyne.material import Material, from_atom_frac\n",
"import numpy as np\n",
"from yt.config import ytcfg; ytcfg[\"yt\",\"suppressStreamLogging\"] = \"True\"\n",
"from yt.mods import *\n",
"from itaps import iBase, iMesh\n",
"from matplotlib import animation\n",
"from JSAnimation import IPython_display\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.backends.backend_agg import FigureCanvasAgg\n",
"from IPython.display import HTML"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"xsc = XSCache([0.026e-6, 0.024e-6], (SimpleDataSource, NullDataSource))"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Laplacian\n",
"\n",
"The functions in the following cell solve for the laplacian ($\\nabla^2$) for any index in in the mesh using a [3 point stencil](http://en.wikipedia.org/wiki/Five-point_stencil) along each axis. This implements relfecting boundary conditions along the edges of the domain."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def lpoint(idx, n, coords, shape, m):\n",
" lidx = list(idx)\n",
" lidx[n] += 1 if idx[n] == 0 else -1\n",
" left = m.structured_get_hex(*lidx)\n",
" l = m.mesh.getVtxCoords(left)[n]\n",
" if idx[n] == 0:\n",
" l = 2*coords[n] - l \n",
" return left, l\n",
"\n",
"def rpoint(idx, n, coords, shape, m):\n",
" ridx = list(idx)\n",
" ridx[n] += -1 if idx[n] == shape[n]-2 else 1\n",
" right = m.structured_get_hex(*ridx)\n",
" r = m.mesh.getVtxCoords(right)[n]\n",
" if idx[n] == shape[n]-2:\n",
" r = 2*coords[n] - r\n",
" return right, r\n",
"\n",
"def laplace(tag, idx, m, shape):\n",
" ent = m.structured_get_hex(*idx)\n",
" coords = m.mesh.getVtxCoords(ent)\n",
" lptag = 0.0\n",
" for n in range(3):\n",
" left, l = lpoint(idx, n, coords, shape, m)\n",
" right, r = rpoint(idx, n, coords, shape, m)\n",
" c = coords[n]\n",
" lptag += (((tag[right] - tag[ent])/(r-c)) - ((tag[ent] - tag[left])/(c-l))) / ((r-l)/2)\n",
" return lptag"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solve in space\n",
"\n",
"The ``timestep()`` fucntion sweeps through the entire mesh and computes the new flux everywhere. This opperation takes place enritely on the mesh object."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def timestep(m, dt):\n",
" nx = len(m.structured_get_divisions(\"x\"))\n",
" ny = len(m.structured_get_divisions(\"y\"))\n",
" nz = len(m.structured_get_divisions(\"z\"))\n",
" shape = (nx, ny, nz)\n",
" D = m.mesh.getTagHandle(\"D\")\n",
" k = m.mesh.getTagHandle(\"k\")\n",
" S = m.mesh.getTagHandle(\"S\")\n",
" Sigma_a = m.mesh.getTagHandle(\"Sigma_a\")\n",
" phi = m.mesh.getTagHandle(\"phi\")\n",
" phi_next = m.mesh.getTagHandle(\"phi_next\")\n",
" for idx in product(*[range(xyz-1) for xyz in shape]):\n",
" ent = m.structured_get_hex(*idx)\n",
" phi_next[ent] = (max(D[ent] * laplace(phi, idx, m, shape) + \n",
" (k[ent] - 1.0) * Sigma_a[ent] * phi[ent], 0.0) + S[ent])*dt*2.2e5 + phi[ent]\n",
" ents = m.mesh.getEntities(iBase.Type.region)\n",
" phi[ents] = phi_next[ents]"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 4
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solve in time\n",
"\n",
"The ``render()`` function steps through time calling the ``timestep()`` function and then creating an image. The images that are generated are then dumped into a movie."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def render(m, dt, axis=\"z\", field=\"phi\", frames=100):\n",
" timestep(m, dt)\n",
" pf = PyneMoabHex8StaticOutput(m)\n",
" s = SlicePlot(pf, axis, field, origin='native')\n",
" fig = s.plots['gas', field].figure\n",
" fig.canvas = FigureCanvasAgg(fig)\n",
" axim = fig.axes[0].images[0]\n",
"\n",
" def init():\n",
" axim = s.plots['gas', 'phi'].image\n",
" return axim\n",
"\n",
" def animate(i):\n",
" s = SlicePlot(pf, axis, field, origin='native')\n",
" axim.set_data(s._frb['gas', field])\n",
" timestep(m, dt)\n",
" return axim\n",
"\n",
" return animation.FuncAnimation(fig, animate, init_func=init, frames=frames, interval=100, blit=False)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Reactor\n",
"\n",
"This setups up a simple light water reactor fuel pin in a water cell. Note that our cells are allowed to have varing aspect ratios. This alows us to be corsely refined inside of the pin, finely refined around the edge of the pin, and then have a differernt coarse refinement out in the coolant."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def create_reactor(multfact=1.0, radius=0.4):\n",
" fuel = from_atom_frac({'U235': 0.045, 'U238': 0.955, 'O16': 2.0}, density=10.7)\n",
" cool = from_atom_frac({'H1': 2.0, 'O16': 1.0}, density=1.0)\n",
" xpoints = np.linspace(0.0, 1.0, 21)\n",
" ypoints = np.linspace(0.0, 1.0, 21)\n",
" zpoints = np.linspace(0.0, 1.0, 11)\n",
" rx = Mesh(structured_coords=[xpoints, ypoints, zpoints], structured=True)\n",
" D = rx.mesh.createTag(\"D\", 1, float)\n",
" k = rx.mesh.createTag(\"k\", 1, float)\n",
" S = rx.mesh.createTag('S', 1, float)\n",
" Sigma_a = rx.mesh.createTag(\"Sigma_a\", 1, float)\n",
" phi = rx.mesh.createTag(\"phi\", 1, float)\n",
" phi_next = rx.mesh.createTag(\"phi_next\", 1, float)\n",
" for ent in rx.structured_iterate_hex(\"xyz\"):\n",
" D[ent] = 1.0 / (3.0 * cool.density * 1e-24 * sigma_s(cool, xs_cache=xsc))\n",
" Sigma_a[ent] = cool.density * 1e-24 * sigma_a(cool, xs_cache=xsc)\n",
" coords = rx.mesh.getVtxCoords(ent)\n",
" phi[ent] = 0.0\n",
" k[ent] = 0.0\n",
" S[ent] = 0.0\n",
" phi_next[ent] = 0.0\n",
" ent = rx.structured_get_hex(len(xpoints)/2, len(ypoints)/2, len(zpoints)/2)\n",
" S[ent] = 4.0e14\n",
" return rx"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"rx = create_reactor()\n",
"render(rx, dt=2.5e-31, frames=40)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"\n",
"\n",
"\n",
"