Out of memory error in `viz` of 3D `RectilinearGrid` from Meshes.jl

Specific Domains Visualization
, , ,
1

I want to plot a grid using Meshes.jl and GLMakie.jl. I am reading the HDF5 data file output of Athena++ astrophysical code using HDF5.jl library. In Athena++ mesh grid is divided into smaller parts called Blocks to use block level parallelization. In code below i create x ,y ad z direction coordinates points for 872 blocks using rand in matrices. I got out of memory error. How can i plot the grid containing x, y and z coordinate points for every block?

using HDF5, Meshes, GLMakie
x = rand(Float32, 18,872)
y = rand(Float32, 4,872)
z = rand(Float32, 16,872)

grid = RectilinearGrid(vec(x), vec(y), vec(z))
viz(grid, alpha = 0.3)
2

The RectilinearGrid expects x, y and z vectors. The resulting grid is the Cartesian product of these vectors. In your example, you are defining matrices instead. Are you sure you don’t want StructuredGrid?

Please check the documentation:

3

I am getting such matrices when i read Athena++ output HDF5 file. In those matrices 872 represents the number of blocks that have been used to create a grid during simulation.

image
image811×573 52.3 KB

image
image529×250 58 KB

4

I have no idea what Athena++ is or does. You have to manipulate the data to feed the expected inputs to Meshes.jl

5

I am using vec() to convert matrix to vector.

I want to create a grid of those points for all blocks. I don’t mind whether it is ‘StructuredGrid’ or not.

6

The grid is successfully created. The out-of-memory error is thrown in the viz call.

Visualization of RectilinearGrid is optimized in the 2D case. If you can plot the grid with Makie.jl directly, please share the code used, and we can add the optimization.

Currently, the visualization of 3D grid with arbitrary spacing is using the fallback visualization recipe for grids, which is memory-hungry.

7

This code creates plot on using only some data points in GLMakie.jl.

using HDF5, Meshes, GLMakie

# Load Athena++ grid from file
filename = "athena_output.00000.h5"
h5 = h5open("sane00.prim.01800.athdf", "r")

@timev x1v = read(h5["x1v"])  # cell centers in x
x2v = read(h5["x2v"])  # cell centers in y
x3v = read(h5["x3v"])  # cell centers in z

close(h5)
@show size(x2v)
# Convert to RectilinearGrid (if the grid is structured and aligned)
#grid = RectilinearGrid(vec(x1v), vec(x2v), vec(x3v))
#return viz(grid, alpha = 0.3)

# Use stride to subsample (reduce memory usage)
stride = 100
x_sample = x1v[1:stride:end]
y_sample = x2v[1:stride:end]
z_sample = x3v[1:stride:end]

# Create Point3f0 objects
points = [GLMakie.Point(x, y, z) for x in x_sample, y in y_sample, z in z_sample]
flat_points = vec(points)  # convert 3D array to 1D vector

# Plotting
fig = Figure(size = (800, 600))
ax = Axis3(fig[1, 1], title = "Grid Points", xlabel = "x", ylabel = "y", zlabel = "z")
scatter!(ax, flat_points, color = :blue, markersize = 9.5)
fig

If i use stride = 1 then also i get out of memory error. My Athena++ data is at this link. Please explain me solution in detail. :speech_balloon: Over this grid later i will plot density or Magnetic fields etc.

8

The problem is that you are instantiating the Cartesian product of x y z. You need to find the function in Makie.jl that plots the corresponding grid directiy without manual construction of points.

9

I am able to plot grid of single block using Meshes.jl. Can i merge all these block’s grids to get final grid?

10

Plotting individual blocks one by one will take an enormous amount of time with multiple calls to the GPU. I repeat: find the Makie.jl function that can plot the grid directly with the x, y and z vectors. If this function exists, we can call it in our viz recipes.

11

This code creates meshes for blocks and merge them. But it hangs a bit. :innocent:

using HDF5, GLMakie

function plot_blockwise_amr(filename::String)
    h5 = h5open(filename, "r")
    num_blocks = read(HDF5.attributes(h5), "NumMeshBlocks")

    fig = Figure(size = (1000, 800))
    ax = Axis3(fig[1, 1], xlabel = "x", ylabel = "y", zlabel = "z")

    for b in 1:num_blocks
        # Load 1D coordinate arrays for block `b`
        x = read(h5["x1v"])[:, b]
        y = read(h5["x2v"])[:, b]
        z = read(h5["x3v"])[:, b]

        # Plot grid lines in each direction (avoid Cartesian product)
        for xi in x
            lines!(ax, [xi, xi], [minimum(y), maximum(y)], [minimum(z), minimum(z)], color = :gray)
        end
        for yi in y
            lines!(ax, [minimum(x), maximum(x)], [yi, yi], [minimum(z), minimum(z)], color = :gray)
        end
        for zi in z
            lines!(ax, [minimum(x), minimum(x)], [minimum(y), maximum(y)], [zi, zi], color = :gray)
        end
    end

    close(h5)
    fig
end

fig = plot_blockwise_amr("sane00.prim.01800.athdf")