This notebook is an entirely self-contained solution to a basic neutron diffision 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,
$\frac{1}{v}\frac{\partial \phi}{\partial t} = D \nabla^2 \phi + (k - 1) \Sigma_a \phi + S$
where
In [1]:
from itertools import product
from pyne.mesh import Mesh, IMeshTag
from pyne.xs.cache import XSCache
from pyne.xs.data_source import CinderDataSource, SimpleDataSource, NullDataSource
from pyne.xs.channels import sigma_a, sigma_s
from pyne.material import Material, from_atom_frac
import numpy as np
from yt.config import ytcfg; ytcfg["yt","suppressStreamLogging"] = "True"
from yt.mods import *
from itaps import iBase, iMesh
from matplotlib import animation
from JSAnimation import IPython_display
import matplotlib.pyplot as plt
from matplotlib.backends.backend_agg import FigureCanvasAgg
from IPython.display import HTML
In [2]:
xsc = XSCache([0.026e-6, 0.024e-6], (SimpleDataSource, NullDataSource))
The functions in the following cell solve for the laplacian ($\nabla^2$) for any index in in the mesh using a 3 point stencil along each axis. This implements reflecting boundary conditions along the edges of the domain.
In [3]:
def lpoint(idx, n, coords, shape, m):
lidx = list(idx)
lidx[n] += 1 if idx[n] == 0 else -1
left = m.structured_get_hex(*lidx)
l = m.mesh.getVtxCoords(left)[n]
if idx[n] == 0:
l = 2*coords[n] - l
return left, l
def rpoint(idx, n, coords, shape, m):
ridx = list(idx)
ridx[n] += -1 if idx[n] == shape[n]-2 else 1
right = m.structured_get_hex(*ridx)
r = m.mesh.getVtxCoords(right)[n]
if idx[n] == shape[n]-2:
r = 2*coords[n] - r
return right, r
def laplace(tag, idx, m, shape):
ent = m.structured_get_hex(*idx)
coords = m.mesh.getVtxCoords(ent)
lptag = 0.0
for n in range(3):
left, l = lpoint(idx, n, coords, shape, m)
right, r = rpoint(idx, n, coords, shape, m)
c = coords[n]
lptag += (((tag[right] - tag[ent])/(r-c)) - ((tag[ent] - tag[left])/(c-l))) / ((r-l)/2)
return lptag
In [36]:
def timestep(m, dt):
nx = len(m.structured_get_divisions("x"))
ny = len(m.structured_get_divisions("y"))
nz = len(m.structured_get_divisions("z"))
shape = (nx, ny, nz)
D = m.mesh.getTagHandle("D")
k = m.mesh.getTagHandle("k")
S = m.mesh.getTagHandle("S")
Sigma_a = m.mesh.getTagHandle("Sigma_a")
phi = m.mesh.getTagHandle("phi")
phi_next = m.mesh.getTagHandle("phi_next")
for idx in product(*[range(xyz-1) for xyz in shape]):
ent = m.structured_get_hex(*idx)
phi_next[ent] = (max(D[ent] * laplace(phi, idx, m, shape) +
(k[ent] - 1.0) * Sigma_a[ent] * phi[ent], 0.0) + S[ent])*dt*2.2e5 + phi[ent]
ents = m.mesh.getEntities(iBase.Type.region)
phi[ents] = phi_next[ents]
In [37]:
def render(m, dt, axis="z", field="phi", frames=100):
pf = PyneMoabHex8StaticOutput(m)
s = SlicePlot(pf, axis, field, origin='native')
fig = s.plots['gas', field].figure
fig.canvas = FigureCanvasAgg(fig)
axim = fig.axes[0].images[0]
def init():
axim = s.plots['gas', 'phi'].image
return axim
def animate(i):
s = SlicePlot(pf, axis, field, origin='native')
axim.set_data(s._frb['gas', field])
timestep(m, dt)
return axim
return animation.FuncAnimation(fig, animate, init_func=init, frames=frames, interval=100, blit=False)
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 allows us to be coarsely refined inside of the pin, finely refined around the edge of the pin, and then have a different coarse refinement out in the coolant.
In [38]:
def isinrod(ent, rx, radius=0.4):
"""returns whether an entity is in a control rod"""
coord = rx.mesh.getVtxCoords(ent)
return (coord[0]**2 + coord[1]**2) <= radius**2
def create_reactor(multfact=1.0, radius=0.4):
fuel = from_atom_frac({'U235': 0.045, 'U238': 0.955, 'O16': 2.0}, density=10.7)
cool = from_atom_frac({'H1': 2.0, 'O16': 1.0}, density=1.0)
xpoints = [0.0, 0.075, 0.15, 0.225] + list(np.arange(0.3, 0.7, 0.025)) + list(np.arange(0.7, 1.201, 0.05))
ypoints = xpoints
zpoints = np.linspace(0.0, 1.0, 8)
# Make Mesh
rx = Mesh(structured_coords=[xpoints, ypoints, zpoints], structured=True)
# Add Tags
rx.D = IMeshTag(size=1, dtype=float)
rx.k = IMeshTag(size=1, dtype=float)
rx.S = IMeshTag(size=1, dtype=float)
rx.Sigma_a = IMeshTag(size=1, dtype=float)
rx.phi = IMeshTag(size=1, dtype=float)
rx.phi_next = IMeshTag(size=1, dtype=float)
# Assign initial conditions
Ds = []; Sigma_as = []; phis = []; ks = [];
for i, mat, ent in rx:
if isinrod(ent, rx, radius=radius):
Ds.append(1.0 / (3.0 * fuel.density * 1e-24 * sigma_s(fuel, xs_cache=xsc)))
Sigma_as.append(fuel.density * 1e-24 * sigma_a(fuel, xs_cache=xsc))
phis.append(4e14)
ks.append(multfact)
else:
Ds.append(1.0 / (3.0 * cool.density * 1e-24 * sigma_s(cool, xs_cache=xsc)))
Sigma_as.append(cool.density * 1e-24 * sigma_a(cool, xs_cache=xsc))
r2 = (rx.mesh.getVtxCoords(ent)[:2]**2).sum()
phis.append(4e14 * radius**2 / r2 if r2 < 0.7**2 else 0.0)
ks.append(0.0)
rx.D[:] = Ds
rx.Sigma_a[:] = Sigma_as
rx.phi[:] = phis
rx.k[:] = ks
rx.S[:] = 0.0
rx.phi_next[:] = 0.0
return rx
In [39]:
rx = create_reactor()
In [43]:
render(rx, dt=2.5e-31, frames=10)
Out[43]: