p4est 2020 HCM Summer School: I/O

We are organizing a p4est summer school sponsored by the Hausdorff Center for Mathematics at the University of Bonn, Germany. Please see also the school's home page and application forms.

Graphics output and mesh load and save

This tutorial is about file system I/O. The first topic is writing specially formatted text files in the VTK format for interactive visualization. The second topic is about writing the current state of the mesh to disk, and about reading it back in, with the same or a different number of MPI ranks. To this end, we use the MPI I/O interface standard.

Dependencies
Build; roughly understand example/simple/; have a running p4est main program to fiddle with.
Required skills
Interactive visualization of VTK files using e.g. Paraview or Visit; reading up on the MPI I/O interface.
Skills acquired
You will have written cell- and corner-based numerical data to graphics files and looked at them in third-party visualization software. We also cover saving and loading the p4est mesh in a parallel partition-independent format, which is great for flexible checkpointing and overall reproducibility.
Next steps
Extend your hello-world program by adding a nonlinear, smooth geometry transformation that is passed to the VTK output routine. Extend the VTK output routine to write one big MPI file in the binary appended VTK format (instead of many VTU files and one PVTU file as we do currently).

VTK graphics output

We are creating and manipulating meshes, which are complex objects embedded in 2- or 3-dimensional space, and want to look at them. To this end, p4est provides the files p{4,8}est_vtk.{c,h} with routines that write files to disk in the VTK graphics format. The easiest call is 

p4est_vtk_write_file (p4est, NULL, "filename");

We take a valid p4est object, query its communicator for metadata, and write the local quadrants using a filename that is postfixed with _rank.vtu, where the rank is zero-padded to four digits. For each quadrant, we write its MPI rank, tree number and level. These are available as input fields within the visualization program. The NULL argument indicates to use the vertex coordinates from the connectivity. You may alternatively specify a p4est_geometry object to provide your own coordinate map.

While the call is not strictly MPI-collective, the output is only complete if every process calls it with the same filename. Rank zero in addition writes a meta file ending in .pvtu that contains a directory of all .vtu files written. You may now use a visualization program to read these files and display them.

Exercise I1: to visualize a connectivity structure, create a minimal p4est at level zero and visualize that. Experiment with varying MPI sizes and play with the wrap_rank parameter.

Exercise I2: Add some calls to p4est_vtk_write_file to your favorite example program, using a different file name on each occasion, and examine the results. Verify the tree number, MPI rank, and quadrant level interactively.

A more powerful interface to writing VTK files works as follows: Create an opaque p4est VTK context, set parameters to the context, write the file header, zero or more data arrays, and the file footer. Writing the footer automatically deallocates the context and closes all files. Please see the inline documentation in the p4est_vtk header file for the prototypes.

This interface allows to add custom data to the output. The data may be provided per corner or per cell (quadrant). Writing per corner requires a consistent corner numbering, which will be discussed in a later tutorial. The number of local quadrants, on the other hand, is easily available from the p4est object, and the input data array must match in length.

Exercise I3: Create an sc_array (see the files sc/src/sc_containers.{c,h} for details) of correct length. Populate some per-quadrant data: Write a loop over the process-local trees, a nested loop over the tree-local quadrants, access each quadrant and define the data written for it as some function of its coordinates, level, tree or sequence number or anything else. Look at the visualization to make sure everything is right. Don't forget to deallocate the array after writing.

Partition-independent file I/O

The p4est object stores the topology of the parallel adapted mesh and optionally some fixed amount of data per quadrant, the size of which is specified with p4est_new or p4est_copy. This data is preserved on repartitioning. It is also visible to the replace callbacks in p4est_refine_ext, p4est_coarsen_ext and p4est_balance_ext, which allows the user to process it appropriately on any local mesh change. It is convenient and generally not harmful to use this data for numerical state, however, for anything more than a couple variables we recommend to store numerical data in application memory. In that case, use your own interpolation/projection logic on adaptation and the p4est_transfer functions on repartition.

Both the mesh and the quadrant-local data can be written to one large MPI file using the collective call 

p4est_save_ext ("filename", p4est, save_data, save_partition);

The boolean parameter save_data enables or disables storing per-quadrant data in the file. The boolean parameter save_partition enables or disables storing the exact partition. If we do not store the partition, the code on loading will equi-partition the quadrants in the file among the processes used for reading. This is fast; there is no need to worry about partitioning in p4est. It will then not be possible to recreate the parallel ownership of each quadrant prior to saving, but then such is rarely required. If, on the other hand, the partition is saved, the header of the file is prepended with this information. The volume of this information is small, but there is the caveat that the file size written, for the same global mesh, will differ between different MPI sizes. Such is not acceptable for strictly partition-independent storage.

On loading the mesh, you may use the call 

p4est = p4est_load_ext ("filename", mpicomm, data_size, load_data,
                        autopartition, 1, NULL, NULL);

This function allocates a new p4est object and reads it from disk. If loading data is enabled, the data size provided must be equal to the data size in the file. If loading data is false, the new p4est object will have zero data size. If autopartition is true, we discard any saved partition and read the quadrants from the file equi-partitioned. If it is false, the call will crash if no saved partition is present.

Exercise I4: Write a function p4est_can_load that determines without crashing whether the file on disk exists and will be acceptable to load for a given set of parameters. Such a function shall only read the header on rank zero and broadcast the result to all other processes. It may optionally populate a structure passed by reference with some statistics and information from the file.

Exercise I5: Extend your favorite program example to save and load a p4est. Create a partition-independent mesh refinement and write the same mesh using different MPI sizes (and file names). Check under which conditions the files are identical, and if this matches with the documentation. Verify that the function indeed crashes if used out of specification.