WIP: Dg multi mesh hdf5 output

[Davknapp]

Develop HDF5-output for DGMultiMeshes. Write Mesh informations such as element_type, vertices, connection between elements and vertices in a seperate file. Write informations about the solution each time the save_solution_callback is called. Both informations will be combined via Trixi2Vtk in the future.

[sloede]

My current developments can be found in draft https://github.com/trixi-framework/Trixi.jl/pull/1201 . I got some difficulties to

• extract the vertices in a format suitable for the hdf5 format
• transform u such that it can be stored in 'file["variables_$v"]

I think I need some advanced "Julia-magic" here, I would be happy about some hints ;)

Regarding the first point I think the best person to get in touch with is @jlchan. Regarding the second point: If you manage to extract the solution of a single variable in a (multi-dimensional) array, you should be able to just dump that array into the HDF5 file.

In the current callback we use a lot of Julia magic to optimize away copies, but you can just start with creating a copy ("make it work") and the optimize it afterwards ("make it nice").

[jlchan]

@Davknapp - to clarify:

extract the vertices in a format suitable for the hdf5 format

Are you intending to extract the vertex values of the high order solution and use Paraview's low order visualization capabilities?

transform u such that it can be stored in 'file["variables_$v"]

What format does u need to be in order to be stored in file["variables_$v"]? Is it a flat array? Vector?

[sloede]

What format does u need to be in order to be stored in file["variables_$v"]? Is it a flat array? Vector?

IIRC, for the other files we use flat arrays. However, this is imho not strictly necessary, since you cannot do anything with this data anyways, unless you know how to interpret it. Since we are also in control of interpreting the data (via Trixi2Vtk), I think a multi-dimensional array would be just as well.

[jlchan]

Thanks for the clarification @sloede. What do you mean by "we are also in control of interpreting the data (via Trixi2Vtk)"? Do we specify elsewhere how to plot from the saved data?

@Davknapp - for a basic implementation, we can either

1. extract vertex nodes from each solution component for low order visualization. This would look like

u_component = StructArrays.component(u, index) 
# `StructArrays.component` can also be replaced by `getindex.(u, index)` for 
# more general behavior. StructArrays.component does not work for `DGMulti`'s 
# `SBP()` solver type due to a different data layout, but `getindex` allocates.

file["variables_$v"] = u_component[vertex_indices, :] 
# rows correspond to nodes, columns correspond to elements. extracting vertex 
# dofs is then equivalent to selecting the appropriate indices. I can add something 
# to do that in `StartUpDG.jl` if that helps.

2. interpolate the solution to a set of high order nodes for high order visualization. This is a little more involved, but would look something like

u_component = StructArrays.component(u, index)
u_interpolated = dg.basis.Vp * u_component

Here, dg.basis.Vp is an matrix which interpolates nodal values to plotting nodes. We could then extract low order elements from u_interpolated (I believe we do this for the Plots.jl and Makie.jl recipes currently).

If you have a preference please let me know!

[ranocha]

What do you mean by "we are also in control of interpreting the data (via Trixi2Vtk)"? Do we specify elsewhere how to plot from the saved data?

The save solution callback just spits out the data into an HDF5 file (at the moment - maybe something else such as ADIOS2 later?) in a format that is not specified further by third-party software. The sole purpose is to be able to read it using Trixi2Vtk to process it for VTK-based visualization. Thus, the save solution callback just needs to store enough data such that we are able to construct VTK files from it later.

[Davknapp]

@Davknapp - to clarify:

extract the vertices in a format suitable for the hdf5 format

Are you intending to extract the vertex values of the high order solution and use Paraview's low order visualization capabilities?

For a first version I aim to extract the vertex values of the low-order solution (flat triangles). Therefore, it is sufficient to know the coordinates of the triangle-corners. To write the mesh in a vtk-format eventually it is most handy to just have an array storing the coordinates of all vertices and another one storing the indices (in the vertex-array) of the corners of each element. Having a look at MeshData these informations should be stored in VXYZ and EToV

[jlchan]

For clarification, would this approach force a continuous solution (since there is a global vertex array)?

VXYZ does indeed contain vertex coordinates, but to extract unique vertex coordinates for the solution would also require some sort of averaging.

[sloede]

I think any postprocessing should be done in Trixi2Vtk.jl. IMHO, the goal for the callback is to dump all information with minimum effort into files, which are then to be read and interpreted by the postprocessing tool(s).

[jlchan]

Hi @Davknapp - sorry for the delay. If you wish to extract the vertex node values for the solution on a 2D triangular mesh, you can access them via

vertex_ids = [1, dg.basis.N+1, dg.basis.Np]
# `u[vertex_ids, :]` returns a size `(3, num_elements)` array whose three rows are 

This assumes a specific ordering of the nodes on triangles, but it should hold for all DGMulti solvers.

If you'd like to do something similar for Quad elements, vertex_ids = [1, rd.N+1, (rd.N+1)^2 - rd.N, rd.Np] will work.