pfvtksave - A new pftool with advanced graphics options

by Nick Engdahl

It's no secret that the outputs of ParFlow have historically been...well, let's be polite and say they have left much to be desired in the area of graphical presentation. Silo is fine, as long as you have it installed, but we've never had a portable (i.e. doesn't require a specific library) graphical output that allows direct visualization. There's also the point that terrain following grids (TFGs) are an internal coordinate transformation of the problem (see this post), so the resulting models look like a flat, constant thickness orthogonal grid (OG), even thought they might not be. Finally, there are also those mysterious single-file CLM outputs (see this other post) and they're just a giant stack of outputs that are equally difficult to visualize.

Starting with ParFlow revision 713, there is now a pftool that is designed specifically to help you see your results faster, in a portable format, that can handle OGs, TFGs, variable layer thicknesses, and single-file CLM outputs. The design of the tool mimics pfsave but with optional flags that specify what the tool will do and this blog post walks you through the tool, called pfvtksave. We'll also cover an example that is now included within the "test/clm" directory to help you get started using it.

Two quick preliminary notes: 1) The visualization toolkit legacy binary format is used as the standard for writing the files, hence the vtk in the name. This is a very widely used format and readily works with VisIt. As long as your specified filename(s) (see below) have a vtk extension, these should be very easy to use. 2) You need to be at or above ParFlow release 713 to use the pfvtksave tool. If you aren't, update and recompile first. There are many blog posts to help guide you through that.

pfvtksave - syntax and options

The basic syntax for the command is:

pfvtksave dataset filetype filename [options]

The order of the required arguments has to follow this format, but the options can be in any order (see below). Without options, this will create a VTK file exactly the same as a PFB. dataset and filename are pretty straight forward and are exactly the same as pfsave.

The filetype must be -vtk OR -clmvtk with the latter writing a file that contains a single layer and each of the output CLM variables are stored as properties by name.

The options:

Any combination of these can be used and they can be specified in any order as long as the required elements immediately follow each option.

-var specifies what the variable written to the dataset will be called. This is followed by a text string, like "Pressure" or "Saturation" to define the name of the data that will be written to the VTK. If this isn't specified, you'll get a property written to the file creatively called "Variable". This option is ignored if you are using -clmvtk since all its variables are predefined.

-dem specifies that a DEM is to be used. The argument following -dem MUST be the handle of the dataset containing the elevations. If it cannot be found, the tool ignores it and reverts to non-dem mode. If the nx and ny dimensions of the grids don’t match, the tool will error out. This option shifts the layers so that the top of the domain coincides with the land surface defined by the DEM. Regardless of the actual number of layers in the DEM file, the tool only uses the elevations in the top layer of this dataset, meaning a 1-layer PFB can be used.

-flt tells the tool to write the data as type float instead of double. Since the VTKs are really only used for visualization, this reduces the file size and speeds up plotting.

-tfg causes the tool to override the specified dz in the dataset PFB and uses a user specified list of layer thicknesses instead. This is designed for terrain following grids and can only be used in conjunction with a DEM. The argument following the flag is a text string containing the number of layers and the dz list of actual layer thicknesses (not dz multipliers) for each layer from the bottom up such as: -tfg "5 200.0 1.0 0.7 0.2 0.1"

Note that the quotation marks around the list are necessary. 

Caution about common errors:

Many programs will insert a slightly different character for an opening vs. closing parentheses and will auto-display a double dash -- as a slightly elongated dash. Either of these cases may cause the tool to not recognize an option as specified, making it look like the tool isn't working (this happens with other pftools too). Basic text editors like vim and mvim do not create any of these issues because they use the normal ascii charters the tool looks for by default. Many GUI text editors are fine as long as you create the tcl within them, but copying and pasting from other places can inadvertently insert long dashes and funky quotes that are hard to see at a glance. If an option doesn't seem to be working, go back, check your text, and you'll probably fix the problem. 

Example script

The following is the bottom part of the clm_vtk.tcl script, now included in "test/clm"

file copy -force CLM_dem.cpfb CLM_dem.pfb

set CLMdat [pfload -pfb clm.out.clm_output.00005.C.pfb]

set Pdat [pfload -pfb clm.out.press.00005.pfb]

set Perm [pfload -pfb clm.out.perm_x.pfb]

set DEMdat [pfload -pfb CLM_dem.pfb]

set dzlist "10 6.0 5.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5"

pfvtksave $Pdat -vtk "CLM.out.Press.00005a.vtk" -var "Press" 

pfvtksave $Pdat -vtk "CLM.out.Press.00005b.vtk" -var "Press" -flt 

pfvtksave $Pdat -vtk "CLM.out.Press.00005c.vtk" -var "Press" -dem $DEMdat 

pfvtksave $Pdat -vtk "CLM.out.Press.00005d.vtk" -var "Press" -dem $DEMdat -flt

pfvtksave $Pdat -vtk "CLM.out.Press.00005e.vtk" -var "Press" -dem $DEMdat -flt -tfg $dzlist

pfvtksave $Perm -vtk "CLM.out.Perm.00005.vtk" -var "Perm" -flt -dem $DEMdat -tfg $dzlist

pfvtksave $CLMdat -clmvtk "CLM.out.CLM.00005.vtk" -flt 

pfvtksave $CLMdat -clmvtk "CLM.out.CLM.00005.vtk" -flt -dem $DEMdat

pfvtksave $DEMdat -vtk "CLM.out.Elev.00000.vtk" -flt -var "Elevation" -dem $DEMdat

The DEM for this example is entirely fictitious and is only used for illustration. It comes packaged with a .cpfb extension to prevent the make clean script fro deleting it, so the first line just copies the DEM over to a .pfb. The next three lines load the various data files and store their handles. 

The 6th line illustrates the preferred way for specifying the dzlist as a tcl variable. Note that these are user specified and do not have to coincide with the dz used for any particular model. For example, the actual clm_vtk.tcl test case uses a constant layer thickness of 0.5m, but we can override this for plotting. The land surface will always be defined by the DEM.

The nine pfvtksave commands illustrate the functional modes of the tool

1. Write the loaded pressure dataset as a flat OG and name the output variable "Press"

2. The same as 1, but compress the output b writing it as float instead of double

3. The same as 1 but use the loaded DEM dataset to define the surface

4. The same as 3, but reduce the size of the output using the float option

5. Do the same as 4, but now use variable thickness layers

6. This is the same as 5, but the permeability dataset is being written

7. Write and compress the CLM outputs as a flat OG

8. This is the same as 7 but uses a DEM to specify the elevations

9. Write the DEM itself as the dataset, compress it, and call it "Elevation"

That's all there is to it. Hopefully, you find this tool useful!

Lastly, a brief technical note most of you can skip:

To reduce the size of the points that have to written for the terrain following (structured points) output, the x, y, z coordinates of the points are always written as float; this only applies when a valid DEM is used. Prior to this tool, there was a bug in the way floats were written, which should be fixed at this point. If there are still errors and the thing won’t work, there is a provision in the source to write the points as double but this requires a comment/uncomment in printdatabox.c in the PrintTFG_VTK() function and pftools will need to be recompiled. The operation should be pretty straightforward once you find that part of the code. 

Using Indicator Files as input to ParFlow

A common question is how to get spatially-varying input, such as soil types, into ParFlow. Though there are a number of ways to do this, a most common approach is to use an indicator pfb file with integer values for each land coverage. This is what we've done for the coupled model where we have different land surface types (such as from IGBP) distributed over the ground surface. An example of this type of coverage is in this paper.

What we do is write the indicator value to the corresponding spatial grid location in a pfb file, then load this file in. An indicator field is a .pfb file with integer numbers for every cell in the computational domain. An example in two-dimensions might look as follows.


These are mapped to geometry names within the parflow .tcl input file and the Tcl/TK input script would look something like:

pfset GeomInput.indinput.InputType IndicatorField
pfset GeomInput.indinput.GeomNames "field channel subsurf"
pfset Geom.indinput.FileName "my_indicator_input.pfb"

pfset GeomInput.field.Value 1
pfset GeomInput.channel.Value 2
pfset GeomInput.subsurf.Value 3
You then can use the names to assign properties, say permeability

pfset Geom.Perm.Names "field channel subsurf"

pfset Geom.field.Perm.Type Constant
pfset Geom.field.Perm.Value 0.2

pfset Geom.channel.Perm.Type Constant
pfset Geom.channel.Perm.Value 0.001

pfset Geom.subsurf.Perm.Type Constant
pfset Geom.subsurf.Perm.Value 1.0
This will assign these parameter values anywhere these indicators are located in the file. You need to distribute the file before reading it in. To write the file you can modify the FORTRAN code pf_write.f90 located in the /helper_codes directory or use another program to write an ascii file and then use pftools to convert it to a .pfb file.

Subsurface Setup Options

by Nick Engdahl & Lindsay Bearup

Keywords: Solid Files, Indicator Files, heterogeneity. 

Manual Location

Manual version: July, 2014 v.693

Defining the Problem – Chapter 3.1

Little Washita Example – Chapter 3.6.2

TFG formulation – Chapter 5.3

ParFlow Solid Files (.pfsol) – Chapter 6.5

So, you’ve decided to use ParFlow and got tired of only being able to make heterogeneities that are rectangular boxes…and now you have no idea what’s going on. Domain setup is tricky, with a lot of options. Usually, we rely on two approaches and we hope this post will help to at least conceptually clarify how those two options differ.

To start, we’re going to assume that you have a processed DEM and that your domain is already defined (slopes, limits, etc…) already as either a terrain following grid (TFG) or an orthogonal (or rectilinear) grid (OG). We’ll also assume that you’ve gotten a homogeneous version of the domain to run, which is a highly recommended first-step in the modeling process. With all that done, the task we’re going to look at is how to add some complexity within those domains with spatial heterogeneity. There are two ways to do this: Indicator Files and Solid Files, and either can be used with both grid types (TFG or OG)

Common Features

The concept here is to define different zones within your model and that’s the common framework that indicators and solid files share. The model will know that certain places in the domain belong to Category 1, for example, and you can use that reference to define any of the properties ParFlow uses. You can think of this as defining hydrofacies or hydrogeologic units. Both kinds of files will need to be created externally then brought in to ParFlow and that’s where the similarities stop.

Indicator Files

An indicator file is a .pfb file that assigns an integer to each cell. In the main ParFlow .tcl script, each integer needs to be linked to a name that is then used throughout the script to assign parameters (permeability, porosity, Van Genuchten parameters, etc). In the Little Washita example, this looks something like this:

pfset GeomInput.indi_input.InputType   IndicatorField

pfset GeomInput.indi_input.GeomNames   "s1 s2 s3 s4”

pfset Geom.indi_input.FileName         "IndicatorFile_Gleeson.50z.pfb"

pfset GeomInput.s1.Value                1

pfset GeomInput.s2.Value                2

pfset GeomInput.s3.Value                3

pfset GeomInput.s4.Value                4

…

pfset Geom.Perm.Names         "domain s1 s2 s3 s4 "

pfset Geom.domain.Perm.Type   Constant

pfset Geom.domain.Perm.Value  0.2

pfset Geom.s1.Perm.Type       Constant

pfset Geom.s1.Perm.Value      0.269022595

pfset Geom.s2.Perm.Type       Constant

pfset Geom.s2.Perm.Value      0.043630356

pfset Geom.s3.Perm.Type       Constant

pfset Geom.s3.Perm.Value      0.015841225

pfset Geom.s4.Perm.Type       Constant

pfset Geom.s4.Perm.Value      0.007582087

Anytime you have specifications like this, there must be complete entries for every entry in the “Names” list or your model won’t run. If you wanted to fill-in the space represented by any of the zones, like “s3” with a distribution of permeability instead of a single value, you can swap “Constant” for turning bands or another option (see the manual for those).

Like other input files, the indicator file needs to be distributed and undistributed to run on multiple processors:

pfdist IndicatorFile_Gleeson.50z.pfb

set runname "LW"

puts $runname

pfrun    $runname

pfundist IndicatorFile_Gleeson.50z.pfb

For more information on indicator files, check out an earlier blog post, here.

Indicator files can be used for OG and TFG domains but it is important to remember that all positions in a TFG are relative to the model top. Numerically, the grid is still just a rectangle but the internal geometric transform using the slopes makes the model think it is tilted. The variable thickness is the same way; it’s a multiplier within the solver, which is unknown at the time the program does its domain setup tasks.

Many people would like a step-by-step guide describing how to make 3-D shapes in an indicator file but there are far too many ways for us to cover. The simplest way is to use some the common indicator geostatistical methods (like sequential indicator simulation, or transition probability based indicators) that produce a 3-D array of indicators. If you convert those outputs in to a PFB using pftools, you’ll be ready to go.

Solid Files

Solid files are a very different approach and are basically exactly what you would expect: solid objects. Instead of defining zones at every cell in the domain, we define the volumetric extents of solids and let ParFlow decide what portion of the domain is inside each. This has the useful property that the solids are independent of any particular grid as long as the bounds of the domain completely contain the solids. For example, if your initial grid uses a dx=5.0 with nx=10, the same solid file can be used for dx=1.0 and nx=50, where dx and nx are the cell spacing and number of cells, respectively.

The most common use of solid files is to define the extent of the domain (see the crater2d.tcl example in the /test directory and this blog post), but they can also be used to identify zones within the domain, just like indicators. In our obviously biased opinion, a great example of this is The East Inlet model described in Engdahl & Maxwell (2015). The different zones corresponding to different geologic environments are clearly visible in Figure 2 of that article and each has been filled in with a Gaussian random permeability field.

Solid files are powerful tools for representing complex shapes but they are not easy to create. This isn’t a ParFlow specific issue; no matter what your application, it can be tricky defining 3-D objects that are complete, enclosed volumes. On the plus side, the pfsol file is plain text, so you can read it. Its contents are a collection of vertices, connections and object id’s that make up a collection of small, planar objects and the outward normal of these triangular planes define the geometry. An earlier blog post explains how a utility in the pftools directory can help with the task. The simplest approach is to use objects with vertical edges along every side but the top and bottom surfaces can vary significantly. Again, the East Inlet model is a good example of the complexity the surfaces can take.

In prepping this post, we realized that none of the ParFlow test cases to date (up to release 710) contain examples of assigning properties using solid files for zones; only defining domain limits. We’ll address this in a future post and add an example that walks you through the creation of the East Inlet domain using solid files for both OG and TFG equivalents. For now, we’ll just look at the necessary keys in crater2d.tcl:

pfset GeomInput.Names        "solidinput $Zones background"

pfset GeomInput.solidinput.InputType  SolidFile

pfset GeomInput.solidinput.GeomNames  domain

pfset GeomInput.solidinput.FileName   crater2D.pfsol

This is used to define the active domain, from which the mask is computed. If you plot the mask, you’ll see only the portions of the domain that reside within the pfsol object. If another solid file was read in and given different GeomName, it could be referred to just like “domain” for any property in the domain. The same procedure shown here is merely repeated to create other zones as geometric objects.

An important point is that when you are defining zones within the domain, and not just the domain limits, the order that the solid files are read in specifies their hierarchical relationships. In other words, if two overlap, the last one that is read in wins. Another point is that strange things will happen if the solid file that defines your domain (if you are using one) doesn’t fit within the computational grid.

Again, we realize this post might not answer all your questions just yet, but stay tuned for a post with new examples soon!

References & Links

Engdahl, N.B. and Maxwell, R.M. (2015) Quantifying changes in age distributions and the hydrologic balance of a high-mountain watershed from climate induced variations in recharge, Journal of Hydrology, 522, p 152-162, doi: http://dx.doi.org/10.1016/j.jhydrol.2014.12.032