The cosmology simulation used to make the beautiful header image contains many Adaptive Mesh Refinement (AMR) levels, many variables, and took a team of professionals many months to make look that good. Let's start a little simpler.

The problem with AMR is that Houdini likes nice, simple, uniform grids - but astrophysicists oftentimes don't. So, how do we render AMR data with Houdini? Here is the formula for Enzo, Athena, and FLASH (FITS coming soon). There are two main parts:

  1. Part One: External Scripting
  2. Part Two: In Houdini

Before starting this tutorial make sure you've installed pyopenvdb as discussed in the Install PyOpenVDB Tutorial and make sure all of your yt and pyopenvdb paths are in your PATH variable as discussed in both the Install PyOpenVDB Tutorial and the Getting Started page.





1.1 The Data

Begin by downloading the Enzo Tiny Cosmology sample dataset from here. Take note of the directory where this is being saved. Go there, and unzip the folder.





1.2 The Code

Download the writeAMRVDB.py Python script from our Bitbucket repository. Take note of the directory where this is being saved.

Open the file in a text editor. Search for the line that starts with datafilename =. Write in the path to the data file you downloaded.

Search for the line that starts with outfilepath =. Write in the path to the directory where you want to write the output VDB files.


First, source your yt installation (which should include pyopenvdb if you completed the pre-requisite tutorials).

To run the code, open a terminal window. Navigate to the directory where you saved the writeAMRVDB.py file. Then write: python writeAMRVDB.py

This should only take a few seconds to run. It will write out three .vdb files to the location you specified with outfilepath.

If you're impatient, or having trouble, you can download the resulting files directly, from here. Download the three "Enzo Tiny Cosmology VDB" files.


#! /usr/bin/env python
import sys
import pyopenvdb as vdb
import numpy as np
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 ###
#########################################################################

datafilename = '/home/kalina/Downloads/enzo_tiny_cosmology/DD0010/DD0010'
outfilepath = '/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 data
ds = yt.load(datafilename)
		      
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
min_level = 0
max_level = int( ds.all_data()['grid_level'].max() )
		      
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.


#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(min_level, max_level+1):

   #select specific level(#level) of grids set from dataset(ds)
   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 FloatGrid.
   data_cube = vdb.FloatGrid()
   mask_cube = 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)):

       sub_grid_all = gs[index]
		      
Each level may contain multiple grids. Here we cycle through all of the grids.


       #extract grid(w/o ghost zone) with specific varible
       sub_grid = sub_grid_all[variable]

       #extract grid(w/ ghost zone) with specific varible
       sub_grid_GZ_raw = sub_grid_all.retrieve_ghost_zones(n_zones=gz_size, fields=variable)
       sub_grid_GZ = sub_grid_GZ_raw[variable]
		      
Here we extract just the variable that we are interested in (e.g. "Density"), with ghost zones.


       #extra mask grid (mask) (eg. {[1 0 0 1],[0 1 0 1]...})
       mask = sub_grid_all.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 = sub_grid_all.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
        data_cube.copyFromArray(sub_grid_GZ,ijk=(ijkout[0],ijkout[1],ijkout[2])) 
        mask_cube.copyFromArray(mask,ijk=(ijkout[0],ijkout[1],ijkout[2]))   
		      
We cope the data and the mask into the VDB containers we have made.


    #Calculate voxel size to make data fit into a unit box (size 1 1 1)
    voxel_numbers = ds.domain_dimensions*ds.refine_by**level
    vsize = 1 / float(voxel_numbers[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.


   #Transformation matrices: scale and translate
   data_matrix = [[vsize,0,0,0],[0,vsize,0,0],[0,0,vsize,0],[-vsize/2,-vsize/2,-vsize/2,1]]
   mask_matrix = [[vsize,0,0,0],[0,vsize,0,0],[0,0,vsize,0],[ vsize/2, vsize/2, vsize/2,1]]
   data_cube.transform = vdb.createLinearTransform(data_matrix)
   mask_cube.transform = vdb.createLinearTransform(mask_matrix)
		      
We scale the data and mask by the same amount with the calculated voxel size.
We translate the data by 1/2 voxel in the negative direction. This is to account for AMR data being cell-centered, whereas in Houdini it is vertex-centered. 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], whereas typcially the axes would align with the volume edges.
We translate the mask by 1/2 voxel in the opposite direction. This is because the data has an extra ghost zone voxel all around its boundary, and the mask doesn't. So we need the data and the mask to be off by 1 voxel.


   #Write out the generated vdb
   output = []
   data_cube.name = 'density'
   mask_cube.name = 'mask'
   output.append(mask_cube)
   output.append(data_cube)
   outvdbname = outfilepath + '/' + variable + '_level' + str(level) + '.vdb'
   vdb.write(outvdbname, grids=output)
		      
Finally, we write out the data and mask to a vdb file.




2.1 The Geometry

Now that you have created your VDB files, open up Houdini.

  1. Use the TAB menu in the network view window to create a Geometry node. (This is an advanced tutorial, so I assume you already know how to do this. If you don't, take a look at an introductory tutorial
  2. to learn the basic steps.)
  3. Click on the "geo1" text next to the geo1 node and rename the node to "geo0", since this will contain our level 0 VDB. This is an important step - if you do this, the rest of the default naming will fall into place correctly.


Double-click to step inside the Geometry node, and point the File node to the path of your level 0 VDB.

Step back out to the object level, and repeat the steps above to create a Geometry and File node for each AMR level.




2.2 The Shader

Download the AMR shader Digital Asset from here.

Load the shader into your project:

  1. Navigate to File -> Import -> Houdini Digital Asset...
  2. Select the downloaded amr_shader.hdanc file, and clicking "Install and Create".
  3. You should now have a node called "amr_shader1" in your shop network. Click on the "amr_shader1" text next to the node and rename it to "amr_shader0".

Hover over image to zoom

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.


Now, we set the data parameters:

VDB Path = The path to Density_level0.vdb, the same file that you pointed to in your Geometry File node.

Min Data Value = 0

Max Data Value = 7

If you are working with a different dataset, you can figure out a good min and max to use by using a Slice node, as described in this tutorial.


Step out of the shop network and back to the obj network.

Select your geo0 node, and go to the Render tab.

Under Material, assign the shader we just created at /shop/amr_shader0

Repeat this for geo1 and geo2 as well, typing in /shop/amr_shader1 and /shop/amr_shader2, respectively. We haven't created these yet, but we will, soon.


Return to the shop view.

Now that we have the beginnings of a shader assigned to an object, we should start being able to see something. Click on the "Render View" tab. Then click the "Render" button. You should see a black and white box of data.

Tip: If you don't see anything, place your mouse over the render view window and press the G key on your keyboard. This will frame the image so that you can see it.


Select your amr_shader0 node. Go to the Color tab and set the Color Field to density. We are using the same data field, density, to drive both the opacity and the color.

Since this data is self-luminescent, we want to give it a higher emission. Set Emission to 10.

The "Emission Color Ramp" area is where you color the data. Double-click on it to set a key and select a color for it. (Or, you could single-click to create a key first, and then double-click to change its color). Create a color map that you like. You should see that your render is updating to show these colors.


We have finished the shader for for a single level, level 0. Now we want to create the other two shaders, but we want to link the parameters to these values, so that we only have to make changes in one place. Right-click on the amr_shader0 node and navigate to Actions -> Create Reference Copy. This will create an referenced copy called "amr_shader1", which links each parameter to amr_shader0.


Select the amr_shader1 node. Right-click on the "VDB Path" field and select "Delete Channels". This will unlink this field from amr_shader0, and will allow you to edit the text. Update the path to point to your "Density_level1.vdb" path.

Make a copy of your new amr_shader1 node (CTRL+C and CTRL+V) to create a amr_shader2 node. Again, update the "VDB Path" to point to the correct vdb file, Density_level2.vdb.


Congratulations, you have rendered AMR data! To tweak the look of it, you can go back and modify the color, opacity, emission, or other shader variables. Apply these changes to amr_shader0, and they will be reflected across all three shaders.




More about shaders

Next Tutorial