Let's Explore PyNE Meshes!

In [2]:
from yt.config import ytcfg; ytcfg["yt","suppressStreamLogging"] = "True"
from yt.mods import *
from pyne.xs.channels import sigma_t
from pyne.material import Material, from_atom_frac
from pyne.mesh import Mesh, IMeshTag, MetadataTag, ComputedTag
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-2-c9a37f39ff47> in <module>()
      1 from yt.config import ytcfg; ytcfg["yt","suppressStreamLogging"] = "True"
----> 2 from yt.mods import *
      3 from pyne.xs.channels import sigma_t
      4 from pyne.material import Material, from_atom_frac
      5 from pyne.mesh import Mesh, IMeshTag, MetadataTag, ComputedTag

/home/scopatz/.local/lib/python2.7/site-packages/yt-3.0dev-py2.7-linux-x86_64.egg/yt/mods.py in <module>()
     99     RAMSESStaticOutput, RAMSESFieldInfo
    100 
--> 101 from yt.frontends.halo_catalogs.api import \
    102     HaloCatalogStaticOutput, HaloCatalogFieldInfo, \
    103     RockstarStaticOutput, RockstarFieldInfo

/home/scopatz/.local/lib/python2.7/site-packages/yt-3.0dev-py2.7-linux-x86_64.egg/yt/frontends/halo_catalogs/api.py in <module>()
     15 #-----------------------------------------------------------------------------
     16 
---> 17 from halo_catalog.api import \
     18      HaloCatalogStaticOutput, \
     19      IOHandlerHaloCatalogHDF5, \

ImportError: No module named halo_catalog.api
In [3]:
import yt
In [4]:
yt.__file__
Out[4]:
'/home/scopatz/.local/lib/python2.7/site-packages/yt-3.0dev-py2.7-linux-x86_64.egg/yt/__init__.pyc'

First, let's create a 2D structured mesh, whose grid points in x & y are Cantor dust. Most of the mesh will have water as a coolant. Some volume elements will have fuel, though.

In [2]:
def cantor(n):
    return [0.] + cant(0., 1., n) + [1.]

def cant(x, y, n):
    if n == 0:
        return []
    new_pts = [2.*x/3. + y/3., x/3. + 2.*y/3.]
    return cant(x, new_pts[0], n-1) + new_pts + cant(new_pts[1], y, n-1)
In [3]:
c5 = cantor(5)
coords = [c5, c5, [0.0, 1.0]]

We access volume elements (ve) on the mesh through a volume element index (ve_idx). This index is a unique integer i that gives a sort ordering to an otherwise unorder mesh. The volume element index is defined on the range from 0 (inclusive) to the number of volumes in the mesh (exclusive). You may access a volume element's material through the mats attribute and indexing with i.

In [4]:
m = Mesh(structured_coords=coords, structured=True)
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)
for i in range(len(m)):
    m.mats[i] = cool
m.mats[len(m)/2] = fuel
m.mats[len(m)/4] = fuel
for i, c in enumerate(c5[:-1]):
    m.mats[i*len(c5)] = fuel

Tags

Tags - sometimes known as fields - are generic way of storing data on a mesh. Tags can be accessed as attributes on the mesh class itself. For example, the density tag may be grabbed via:

In [5]:
m.density 
Out[5]:
MaterialPropertyTag(name='density', doc='the density [g/cc]')

To get the value of the density, you have to provide the location by its index.

In [6]:
m.density[42]
Out[6]:
1.0

You can also slice, mask, or fancy index tags:

In [7]:
print m.density[::100]  # slice
[ 10.7   1.    1.    1.    1.    1.    1.    1.    1.    1.    1.    1.
   1.    1.    1.    1.   10.7   1.    1.    1.    1.    1.    1.    1.
   1.    1.    1.    1.    1.    1.    1.    1.   10.7   1.    1.    1.
   1.    1.    1.    1. ]
In [8]:
print m.density[m.density[:] >= 10]  # mask
[ 10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7
  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7
  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7
  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7
  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7  10.7
  10.7  10.7  10.7  10.7]
In [9]:
print m.density[[10, 0, 11, 100]]  # fancy index is fancy
[  1.   10.7   1.    1. ]

To get a list of all current tag names, uses the tags dictionary on the mesh:

In [10]:
m.tags.keys()
Out[10]:
['sub_act',
 'sub_lan',
 'number_density',
 'density',
 'comp',
 'expand_elements',
 'to_atom_frac',
 'molecular_weight',
 'mass_density',
 'atoms_per_mol',
 'sub_tru',
 'mass',
 'attrs',
 'mult_by_mass',
 'sub_ma',
 've_idx',
 'sub_fp']

Analysis & Visulaization

PyNE meshes are supported by the yt project. Use yt's plotting infrastructure to display tags on your mesh.

In [11]:
pf = PyneMoabHex8StaticOutput(m)
s = SlicePlot(pf, 'z', 'density', origin='native')
s.display()

Computed Tags

Computed tags - also known as 'derived fields' - are a way of having a lazily evaluated 'virtual tag.' A computed tag is defined with a function, lambda, or any other callable object. The function must have the following interface:

def f(mesh, i):
    """mesh is a pyne.mesh.Mesh() object and i is the volume element index
    to compute.
    """
    # ... do some work ...
    return anything_you_want

Here is a somewhat silly example which squares the density.

In [12]:
m.density2 = ComputedTag(lambda mesh, i: mesh.density[i]**2)

Note that the bounds on the color bar have changed.

In [13]:
pf = PyneMoabHex8StaticOutput(m)
s = SlicePlot(pf, 'z', 'density2', origin='native')
s.display()

Here is a more serious example which uses PyNE's cross section tools to compute the one-group total cross section $\sigma_t$ everywhere on the mesh.

In [14]:
m.sigma_t = ComputedTag(lambda mesh, i: sigma_t(mesh.mats[i], group_struct=[10.0, 1e-6], phi_g=[1.0])[0])
In [15]:
pf = PyneMoabHex8StaticOutput(m)
s = SlicePlot(pf, 'z', 'sigma_t', origin='native')
s.display()

Yes, that was only one line of code.