VTK support

[henry2004y]

1. The VTK files does not have timestep information. To allow for further time series processing in Paraview, a script create_pvd.jl is provided for generating the pvd container that stores the filenames and timestamps calculated from the filenames. This then becomes the new function create_pvd in vtk.jl. After dicussing with the author of WriteVTK.jl, the timestamp information is now added to each file through a scalar that is not related to grid in VTK, i.e. "global data" in the VTK world.

2. Currently, multi-block (VTM), rectilinear (VTR), and unstructured (VTU) conversions are supported. By default the file size will be reduced with compression level 6, but the actual compression ratio depends on the original data.

3. Currently, BATSRUS output contains only cell center coordinates and cell center values (referring as tcp in PARAM.in), or node coordinates and interpolated nodal values (referring as tec). It is recommended to save in the tcp format to avoid interpolation. In principle VTK also supports the combination of node coordinates and cell center values, but it is not necessary here.

4. The simple multiblock VTK format conversion can only deal with data within a block, but it cannot figure out the connectivities between neighboring blocks. To fill the gap between blocks, we have to retrieve the tree data stored in .tree files. This requires a in-depth understanding of the grid tree structure in BATSRUS, and it took me a whole week to think, write and debug the code!

Information contained in iTree_IA:

• the status of the node (used, unused, to be refined, to be coarsened, etc.);
• the current, the maximum allowed and minimum allowed AMR levels for this node;
• the three integer coordinates with respect to the whole grid;
• the index of the parent node (if any);
• the indexes of the children nodes (if any);
• the processor index where the block is stored for active nodes;
• the local block index for active nodes.

Basically, it requires several steps:

1. Obtain the neighbor block indexes.
2. Obtain the relative AMR level between neighbor blocks.
3. Calculate the global cell indexes.
4. Fill in the connectivity list.

Several issues worth noticing:

• The blocks are load balanced in Morton curve ordering.
• Rules are needed to decide which block is in charge of writing the connectivities. Following the ModWriteTecplot implementation,
  • connectivities between blocks on the same AMR level are written by the left neighbors.
  • connectivities between blocks on different AMR levels are written by the finer blocks. This choice may simplify the logic in writing.
• 2D formatted outputs (with dxPlot⩾0) will generate duplicate points for 2D runs after postprocessing. I don't fully understand why.
• 2D unformatted outputs (with dxPlot<0) can generate correct point coordinates, but the output points in postprocessing are sorted based on x+e*y+e^2*z to allow removing duplicate points. This is mostly useful for 2D cuts of 3D simulations where duplicate points are possible to exist, and should be merged/averaged. This also means that it is literally impossible to figure out the original cell ordering, and impossible to convert to 2D VTU format (unless you construct some new connectivities based on triangulation)! However, it guarantees that the cut outputs from parallel runs are identical.
• 3D unformatted output as stored as it is.
• The simulation is typically done with multiple MPI processes. In the tree structure, each node only records the MPI process local block index. What makes it more complicated is that the local block indexes on a given MPI process may have skipped numbers because they are not in the physical simulation domain. This means that in order to get the true global block indexes, one needs to count how many blocks are used on all MPI processes with lower ranks, and also the order of the current node based on local block index on this MPI process.
• It would be pretty difficult to count the number of elements in the connectivity list. Appending to an array is a very expensive operation: counting the elements in the first loop, preallocating the array and then repeat the same computation to fill in the list results in over 1000 times speedup.
• The implementation in v0.2.3 inherits the bug from BATSRUS ModWriteTecplot that will generate duplicate voxel grid. For example, imagine a coarse block which is surrounded by refined blocks all around. On edge 8, the two refined face neighbor will both think they are in charge of writing the connectivity for the cells around the edge. This does not cause any serious issues as of now, so I don't fix it. Speaking of fixing that, there are two ways:
  • fillCellNeighbors! not only deals with coarsened neighbors, but also the refined neighbors. This means that literally we need to implement the message_pass_cell completely. We still do a two-round process, write everything without checking, and in the end before IO remove the duplicate rows. This may be even faster then implementing some complicated writing rules to avoid duplicating.
  • Polish the writing rules even further. This requires me to consider all possible cases. As a first step, for the above example case, neither of the two face neighbor should write the connectivity, but instead the edge neighbor should write.
• Many function names are inherited from BATL. I may consider rename them in the future if more general task is required.

[henry2004y]

I will close this for now. Reopen it if some new features are required in the future.