AMR Data Data Within Data Within Data

#! /usr/bin/env python
import sys
import pyopenvdb as vdb
import yt
		      

We start by importing the libraries that we will be using. Pretty straightforward, nothing terribly exciting going on here.

#########################################################################
### Modify these values to point to your own data on your own machine ###
#########################################################################

dataFilePath = '/home/kalina/Downloads/enzo_tiny_cosmology/DD0010/DD0010'
outFileDir = '/fe2/amr/data/tutorial'
variable = 'Density'
isFlash = False

#########################################################################
#########################################################################
		      

These are the editable parameters. Modify these values to work with other datasets and write out other variables. FLASH data works a tiny bit differently than Enzo and Athena, so you need to set "isFlash = True" if you are working with FLASH data.

Psst: If you don't know what variables are available in a dataset, after calling ds = yt.load() in the next step, you can call ds.field_list and ds.derived_field_list to see all of your options.

# Load dataset
ds = yt.load(dataFilePath)
		      

We use yt to read the data. Here is where you can add a print(ds.field_list) and/or print(ds.derived_field_list) if you want a print out of what variables are available to you.

# More advanced parameters here
minLevel = 0
maxLevel = ds.index.max_level
vsize0 = None
		      

By default, we want to write out all of the levels is a data set. Here, we start at level 0, and go to the maximum level that is available. You can change these values if you do not want to write all of the levels - though note that this will still cut out "holes", even when you don't write out anything to fill them. A future tutorial will talk about how to fix this.
We will describe the purpose of vsize0 later on.

# Error checking: is this variable in the data?
if not [item for item in ds.field_list if item[1] == variable]:
   print >> sys.stderr, 'ERROR: Invalid field name: ' + variable
   exit()
		      

Here we are doing a check to make sure that the specified variable is actually in this dataset. If it isn't, we exit the script.

# This is required to be able to write out ghost zones for FLASH data
if isFlash:
    ds.periodicity = (True, True, True)
		      

This is a little hack. Normally, you can't get ghost zones around the very outside edges of FLASH data. But, you can by saying that the data is periodic (even though it isn't).

# Iterate through all levels
for level in range(minLevel, maxLevel+1):

   # Select specific level of grids set from dataset
   gs = ds.index.select_grids(level)
		      

We iterate over each level of the data separately, because Houdini needs us to write out each as a separate file.

   # Initiate OpenVDB FloatGrids
   maskCube = vdb.FloatGrid()
   dataCube = vdb.FloatGrid()
		      

We need to write out the data itself, as well as a mask, to tell Houdini where to cut "holes" in the data, so the AMR levels fit nicely within each other. Here we initialize the mask and data grids.

   # Go over all grids in current level
   for index in range(len(gs)):

       subGrid = gs[index]
		      

Each level may contain multiple grids. Here we cycle through all of the grids.

       # Extract grid (without ghost zone) with specific varible
       subGridVar = subGrid[variable]

       # Extract grid (with ghost zone) with specific variable
       subGridVarWithGZ = sub_grid_all.retrieve_ghost_zones(n_zones=1, fields=variable)[variable]
		      

Here we extract just the variable that we are interested in (e.g. "Density"), with ghost zones.

       # Extract mask grid (eg. {[1 0 0 1],[0 1 0 1]...})
       mask = subGrid.child_mask
		      

Here we extract the child mask, so we know where to make "holes" in this level, for the next level to be able to fit into it.

        # ijkout is the global x,y,z index in OpenVDB FloatGrid
        ijkout = subGrid.get_global_startindex()
		      

We need to know the position of each level.

        # Copy data from grid to OpenVDB FloatGrid starting from global x,y,z index in OpenVDB FloatGrid
        maskCube.copyFromArray(mask,ijk=(ijkout[0],ijkout[1],ijkout[2])) 
        dataCube.copyFromArray(subGridVarWithGZ,ijk=(ijkout[0],ijkout[1],ijkout[2]))   
		      

We cope the data and the mask into the VDB containers we have made.

    # Calculate a reasonable voxel size
    resolution = ds.domain_dimensions*ds.refine_by**level
    vSize = 1 / float(resolution[0])
		      

Houdini doesn't like very large objects, and loses precision as positions get farther away from 0. So, we compute a voxel size that will make the data fit into 1 Houdini unit.

   # Keep track of level 0 voxel size              
   if level==minLevel:
       vSize0 = vSize
		      

We want to keep track of the voxel size of level 0. This is so that we can translate all of our grids by this amount, to account for there being an extra ghost zone voxel. If we didn't shift, that would make a ghost zone voxel be at the origin, instead of having our first data value at the origin.

                      
   # Scale and translate
   dataMatrix = [[vSize,0,0,0],[0,vSize,0,0],[0,0,vSize,0],[-vSize/2-vSize0,-vSize/2-vSize0,-vSize/2-vSize0,1]]
   maskMatrix = [[vSize,0,0,0],[0,vSize,0,0],[0,0,vSize,0],[vSize/2-vSize0, vSize/2-vSize0, vSize/2-vSize0,1]]
   dataCube.transform = vdb.createLinearTransform(dataMatrix)
   maskCube.transform = vdb.createLinearTransform(maskMatrix)
		      

We scale the data and mask by the same amount, vsize, with the calculated voxel size.

First we translate the data by -vSize/2. This is to account for AMR data being cell-centered, whereas in Houdini it is vertex-centered. Then we subtract vSize0, the size of a level0 voxel, to adjust for the extra ghost zone voxel that isn't actually in the data. I.e., This makes it so that the origin (0,0,0) in Houdini goes through the center of data voxel at [0][0][0].

We translate the mask by +vSize/2, which is in the opposite direction. This is because the data has an extra ghost zone voxel, but the mask doesn't. So we need the data and the mask to be off by 1 voxel. Then we also shift this by vSize0, so the data and mask line up.

   # Write out the generated VDB
   output = []
   dataCube.name = 'density'
   maskCube.name = 'mask'
   output.append(maskCube)
   output.append(dataCube)
   outFilePath = '%s/%s_level%d.vdb' % (outFileDir, variable, level)
   vdb.write(outFilePath, grids=output)
		      

Finally, we write out the data and mask to a vdb file.

This shader is based on the "Billowy Smoke" material, which you learned about in the More About Shaders tutorial. If you double-click on the node to step inside, the nodes that are in red are the ones relevant to AMR data.

Let's go over what they are doing.

"vdbpath" is the one parameter that we need to import and use in the shader.

"global2" and "transform2": We are taking the position and transforming it from current to world coordinates.

"volumepostoindexfile1" and "volumeindexfile1": We are using the position to get the index in the volume, and then using the index to look up the value of the mask (0 or 1). This may seem obtuse, but there is a reason we are not simply using the position to look up the value -- Mantra/Houdini automatically interpolates volumes to create a smooth falloff at the edges. In this case, we do not want that interpolation; we want a clean cut so that the AMR boxes fit perfectly within each other. So, by looking up the mask value by index rather than by position, we ensure that edges and corners do not fall off to zero, but remain constant through the entire voxel.

"multiply5": We are multiplying the data by 0 or 1 from the mask, to decide whether or not to display anything at that position. Then we pass this value to the rest of the shader.