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 optionsThe 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.
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.
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.
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.