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 gives correct scatter plot but i want 3D volume plot. My grid is non-uniform so i have trouble getting volume plot. Please guide me how to get volume plot?

using HDF5, GLMakie, Colors, LinearAlgebra
GLMakie.activate!(ssao=true)
GLMakie.closeall() # close any open screen

function plothdf(h5file::String)

    h5 = h5open(h5file, "r")
    num_blocks = read(HDF5.attributes(h5), "NumMeshBlocks")
    ρ = read(h5["prim"])[:,:,:,:,1]  # shape: (Nx, Ny, Nz, Nvars, Nblocks)
    x_all = read(h5["x1v"])          # shape: (Nx, Nblocks)
    y_all = read(h5["x2v"])
    z_all = read(h5["x3v"])
    close(h5)

    all_xs, all_ys, all_zs, all_density = Float32[], Float32[], Float32[], Float32[]

    for b in 1:num_blocks
        x = x_all[:, b]
        y = y_all[:, b]
        z = z_all[:, b]
        density = ρ[:,:,:,b]  # 3D array

        for i in 1:length(x), j in 1:length(y), k in 1:length(z)
            r, θ, ϕ = x[i], y[j], z[k]
            X = r * sin(θ) * cos(ϕ)
            Y = r * sin(θ) * sin(ϕ)
            Z = r * cos(θ)

            push!(all_xs, X)
            push!(all_ys, Y)
            push!(all_zs, Z)
            push!(all_density, density[i, j, k])
        end
    end

    # Normalize density
    norm_density = log10.(abs.(normalize!(all_density)))
    cmap = cgrad(:plasma)
    colors = get.(Ref(cmap), norm_density)
    fig = Figure(size = (1000, 800))
    ssao = Makie.SSAO(radius = 5.0, blur = 3)
    ax = LScene(fig[1, 1], scenekw = (ssao=ssao,))
    ax.scene.ssao.bias[] = 1.95
    plt = scatter!(ax, all_xs, all_ys, all_zs, color = norm_density, markersize = 5)
    Colorbar(fig[1, 2], plt, label = "log₁₀(norm(ρ))")
    display(fig)
end
plothdf("sane00.prim.01800.athdf")
12

In Makie.jl, the volume function does not natively support non-uniform grids via direct x, y, z vectors — it assumes a uniform grid internally. My x, y and z vectors contain non-uniform data. How can i get uniform grid using interpolation? My x, y and z vectors are points of Mesh with many level of refinement.

image
image1035×646 122 KB

13

Yes, that is what I remember from the last time I looked into it. If Makie.jl doesn’t offer a function for this curently, we can’t use it to optimize the vizualization of 3D irregular grids.

1 Like
14

The way to go in these scenarios is usually to interpolate onto a regular grid and use that

1 Like
15

That is a good point. You can use the InterpolateNeighbors transform to interpolate your RectilinearGrid values into a RegularGrid and then call viz. The transform is documented in GeoStats.jl and the GDSJL book.