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.
