# Improve performance of this index creation and accessing code

**URL:** <https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235>\
**Category:** Performance\
**Tags:** indexing, rounding\
**Created:** [November 30, 2025, 3:59pm UTC](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235 "2025-11-30T15:59:23Z")\
**Posts on this page:** 6\
**Page:** 4

<div class="post-metadata">

**Author:** ![raman\_kumar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raman_kumar/32/26782_2.png) [@raman\_kumar](https://discourse.julialang.org/u/raman_kumar)\
**Post date:** [January 4, 2026, 4:18pm UTC](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235/62 "2026-01-04T16:18:32Z")

</div>

> [@eldee](#):
>
> Presumably because you are not interpolating across block boundaries.

Yes, I can’t interpolate them as boundaries are at the starting and ending positions of `re, θe, ϕe` values but `Bp_block` values are given for `r, θ, ϕ` positions. I should do extrapolation.

> [@eldee](#):
>
> I don’t know what the concepts of cell-interface or cell-center positions mean.

**Cell-interfaces** are the corner positions of rooms/cells while **cell-centers** are the _centroid /centre of mass_ positions of rooms/cells.

I interpolate _cell-centred_ `Bp_block` values on _cell-centred_ `r, θ, ϕ` positions but to get plots without gap(see photo in last reply) i need to extrapolate them to `re, θe, ϕe` cell-interface or room corner positions.

> [@eldee](#):
>
> However, if the blocks overlap you can (and probably should?) perform your interpolation in multiple blocks,

Blocks overlap in `re, θe, ϕe` coordinate values. If the blocks are adjacent then atleast one of the `re, θe, ϕe` coordinate’s start/end point will overlap and have same value. But the problem is that `Bp_block` values are given at `r, θ, ϕ` positions and these `r, θ, ϕ` positions don’t overlap so we have to extrapolate to `re, θe, ϕe` positions. Try following code for block iteration and adjacent blocks will show some common values.

```julia-auto
@show first(re),last(re),first(θe),last(θe),first(ϕe),last(ϕe)

```

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [January 5, 2026, 9:18am UTC](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235/63 "2026-01-05T09:18:40Z")

</div>

> [@raman\_kumar](#):
>
> **Cell-interfaces** are the corner positions of rooms/cells while **cell-centers** are the _centroid /centre of mass_ positions of rooms/cells.

As far as I can tell the `r` entries are just the midpoints of the `re` entries, making its inclusion in the dataset completely redundant.

If each cell were store a centre of mass of the physical black hole stuff it contains, these centres would not lie on a nice grid, i.e. you would not be able to represent them using three `Vectors`.

> [@raman\_kumar](#):
>
> Blocks overlap in `re, θe, ϕe` coordinate values. If the blocks are adjacent then atleast one of the `re, θe, ϕe` coordinate’s start/end point will overlap and have same value.

I would not call such blocks ‘overlapping’. I guess this depends on your convention, but to me the `re` define a partition [re\_1, re\_2), [re\_2, re\_3), \ldots. By overlapping I rather meant e.g. grids spanning the same area, but with different resolutions, like [0, 1/2), [1/2, 1) and [0, 1/3), [1/3, 2/3), [2/3, 1). Or simply arbitrarily located measurement data, e.g. one block is in [1/3, 2/3), and another in [\pi/12, e/6).

> [@raman\_kumar](#):
>
> Yes, I can’t interpolate them as boundaries are at the starting and ending positions of `re, θe, ϕe` values but `Bp_block` values are given for `r, θ, ϕ` positions. I should do extrapolation.

The mere fact that you want to remove gaps, shows you need interpolation across different blocks. E.g. if your first block has data at 0 and 1/2 and your second block at 1 and 3/2, then to get a value at 3/4, you need to take into account both the first block’s value at 1/2 and the second block’s at 1.

Apart from being conceptually ill-advised, using extrapolation risks that your per-block extrapolated values don’t meet properly. I.e. you trade gaps for inconsistencies (which might still manifest as gaps).

* * *

> [@raman\_kumar](#):
>
> We have to calculate extrapolation on `re, θe, ϕe = xe_all[:, b], ye_all[:, b], ze_all[:, b]`.

> [@eldee](#):
>
> Why? If you just want to visualise your data, you don’t need this extrapolation.

> [@raman\_kumar](#):
>
> Because `re, θe, ϕe`(cell interface positions) contains points which outside the `r, θ, ϕ`(cell-center positions).

> [@eldee](#):
>
> But the question remains, if all you want to do is visualise the data, why do you apparently need to get values at specific predetermined positions (the `_e` versions)?

> [@raman\_kumar](#):
>
> i need to extrapolate them to `re, θe, ϕe` cell-interface or room corner positions.

This is clearly going nowhere, so I’m going to throw in the towel on this topic.

---

<div class="post-metadata">

**Author:** ![raman\_kumar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raman_kumar/32/26782_2.png) [@raman\_kumar](https://discourse.julialang.org/u/raman_kumar)\
**Post date:** [January 5, 2026, 12:42pm UTC](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235/64 "2026-01-05T12:42:25Z")

</div>

> [@eldee](#):
>
> centre of mass of the physical black hole stuff it contains, these centres would not lie on a nice grid

These cells contain average value of quantity that lies inside the cell. This data is from numerical simulation using Athena++ astrophysical code. This [research paper pdf](https://iopscience.iop.org/article/10.3847/1538-4365/ab929b/pdf) describes how blocks are positioned in mesh using **LogicalLocations** by [**Z-order curve**](https://en.wikipedia.org/wiki/Z-order_curve). **Levels** denotes the relative size of blocks.

```julia-auto
Levels = read(h5["Levels"])
LogicalLocations = read(h5["LogicalLocations"])

```

 ![apjsab929bf3_hr](https://global.discourse-cdn.com/julialang/original/3X/b/a/bacf6e74df7a6ccae25946f3ca922cb6c0bb5a2e.jpeg)

> [@eldee](#):
>
> The mere fact that you want to remove gaps, shows you need interpolation across different blocks.

Yes, if it can be done then we can do interpolation accross different blocks. Then we need not to do extrapolation of **re** values. Please suggest how to do it. It is an important problem I am trying to solve since many months.

> [@eldee](#):
>
> This is clearly going nowhere, so I’m going to throw in the towel on this topic.

Sorry, Now I understood that we can do it without extrapolation. If possible then split this topic as it has diverted from initial question of improving performance.

---

<div class="post-metadata">

**Author:** ![raman\_kumar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raman_kumar/32/26782_2.png) [@raman\_kumar](https://discourse.julialang.org/u/raman_kumar)\
**Post date:** [January 8, 2026, 11:59pm UTC](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235/65 "2026-01-08T23:59:38Z")

</div>

This gives plot with gaps as i am interpolating over cell centres x1v, x2v,x3v etc. If i could interpolate over whole HDF5 data once without interpolating for each block iteratively then these gaps in sphere can be filled.

```julia-auto
using HDF5, Interpolations, GLMakie, BenchmarkTools

h5 = h5open("sane98.prim.01900.athdf", "r")
x1v = h5["x1v"] # (Nx, blocks)
x2v = h5["x2v"] # (Ny, blocks)
x3v = h5["x3v"] # (Nz, blocks)
B = h5["B"] # (Nx, Ny, Nz, blocks, 3)
MeshBlockSize::Vector{Int32} = read((HDF5.attributes(h5))["MeshBlockSize"])
blocks::Int32 = read(HDF5.attributes(h5),"NumMeshBlocks") 
Nblocks = 1:blocks

function compute_bounds(x1v, x2v, x3v, Nblocks) 
    xmin = Inf32; xmax = -Inf32; ymin = Inf32; ymax = -Inf32; zmin = Inf32; zmax = -Inf32
    for b in Nblocks
        @views r, theta, phi = x1v[:, b], x2v[:, b], x3v[:, b]
        rvals = extrema(r)
        sinθs = extrema(sin, theta)
        cosθs = extrema(cos, theta)
        sinϕs = extrema(sin, phi)
        cosϕs = extrema(cos, phi)        
		new_xmin, new_xmax = extrema(r*sinθ*cosϕ for r in rvals for sinθ in sinθs for cosϕ in cosϕs)
        xmax = max(xmax, new_xmax); xmin = min(xmin, new_xmin)
		new_ymin, new_ymax = extrema(r*sinθ*sinϕ for r in rvals for sinθ in sinθs for sinϕ in sinϕs)
        ymax = max(ymax, new_ymax); ymin = min(ymin, new_ymin)
        new_zmin, new_zmax = extrema(r*cosθ for r in rvals for cosθ in cosθs)
        zmax = max(zmax, new_zmax); zmin = min(zmin, new_zmin)
    end
    @show xmin, xmax, ymin, ymax, zmin, zmax
    return xmin, xmax, ymin, ymax, zmin, zmax
end

xmin, xmax, ymin, ymax, zmin, zmax = compute_bounds(x1v, x2v, x3v, Nblocks)

function block_interpolator(x1,x2,x3, Bblock::Array{<:Real,3})
    interpolate((x1, x2, x3), Bblock, Gridded(Linear()))
end

function cartesian_to_spherical(x, y, z)
    r = hypot(x, y, z)
    θ = r == 0 ? 0.0 : acos(z/r)
	ϕ = mod2pi(atan(y, x)) # In [0, 2pi) (like ϕs)
    return r, θ, ϕ
end

interps = Vector{Any}(undef, blocks)

for b in Nblocks
    @views begin
        x1 = x1v[:, b]
        x2 = x2v[:, b]
        x3 = x3v[:, b]
        Br = B[:, :, :, b, 1]
        Bθ = B[:, :, :, b, 2]
        Bϕ = B[:, :, :, b, 3]
        Bmag = @. hypot(Br, Bθ) / Bϕ
		interps[b] = extrapolate(block_interpolator(log.(x1), x2, x3, Bmag),NaN32) 
    end
end

#=for b in Nblocks
	itp = interps[b](log(200), 3, 3)
	itp !== NaN32 ? println("$b	$itp") : continue
end=#
steps = 50
@show rx= length(xmin:steps:xmax); ry= length(ymin:steps:ymax); rz = length(zmin:steps:zmax)
field = fill(NaN32, (rx,ry,rz))

function fill_field(xmin, xmax, ymin, ymax, zmin, zmax, field, step, interps, Nblocks)
	for (i,x) in enumerate(xmin:steps:xmax), (j,y) in enumerate(ymin:steps:ymax), (k,z) in enumerate(zmin:steps:zmax)
		r, θ, ϕ = cartesian_to_spherical(x, y, z)
		for b in Nblocks
			itp = interps[b](log(r), θ, ϕ)
			itp !== NaN32 ? field[i,j,k]= itp : continue
		end
	end
	return field
end

field = fill_field(xmin, xmax, ymin, ymax, zmin, zmax, field, step, interps, Nblocks)

fig = Figure()
ax = LScene(fig[1,1])
volume!(ax, xmin..xmax, ymin..ymax, zmin..zmax, field)
display(fig)

```

---

<div class="post-metadata">

**Author:** ![raman\_kumar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raman_kumar/32/26782_2.png) [@raman\_kumar](https://discourse.julialang.org/u/raman_kumar)\
**Post date:** [January 11, 2026, 11:40pm UTC](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235/66 "2026-01-11T23:40:58Z")

</div>

I found a trick to fill gap by appending arrays at end to match size of `Bp` to cell interface values size. See function `append_repeat_last()`in code given below. You can try code below on [this data](https://szkolaglowna-my.sharepoint.com/:u:/g/personal/ramankumar_student_sgmk_edu_pl/IQBkghVdHym4R4ZADnuvDmJIAfJHrVWlAeuoMiZDw8FxtVI?e=M0c1nc). I think appending arrays in `Bp` will have same effect as of `Flat()` in extrapolation. Am i correct?

> **code appending array**
>
> ```julia-auto
> using HDF5, Interpolations, GLMakie, Base.Threads, Statistics
> 
> @time begin
> function readhdf(h5file::String)
> h5 = h5open(h5file, "r")
> @views x1v::Matrix{Float32}, x2v::Matrix{Float32}, x3v::Matrix{Float32}= Array(h5["x1f"]), Array(h5["x2f"]), Array(h5["x3f"])
> B = h5["B"] # (Nx, Ny, Nz, blocks, 3)
> @views Br::Array{Float32, 4}, Bθ::Array{Float32, 4}, Bϕ::Array{Float32, 4} = B[:,:,:,:,1], B[:,:,:,:,2], B[:,:,:,:,3]
> Bp = @. hypot(Br, Bθ) / Bϕ
> blocks::Int32 = read(HDF5.attributes(h5),"NumMeshBlocks") 
> Nblocks = 1:blocks
> return Bp, x1v, x2v, x3v, blocks, Nblocks
> close(h5)
> end
> 
> Bp, x1v, x2v, x3v, blocks, Nblocks = readhdf("sane98.prim.01900.athdf")
> 
> function compute_bounds(x1v, x2v, x3v, Nblocks) 
> xmin=Inf32; xmax=-Inf32; ymin=Inf32; ymax=-Inf32; zmin=Inf32; zmax=-Inf32
> for b in Nblocks
> @views r, theta, phi = x1v[:, b], x2v[:, b], x3v[:, b]
> rvals = extrema(r)
> sinθs = extrema(sin, theta)
> cosθs = extrema(cos, theta)
> sinϕs = extrema(sin, phi)
> cosϕs = extrema(cos, phi)        
> new_xmin::Float32, new_xmax::Float32 = extrema(r*sinθ*cosϕ for r in rvals for sinθ in sinθs for cosϕ in cosϕs)
> xmax = max(xmax, new_xmax); xmin = min(xmin, new_xmin)
> new_ymin::Float32, new_ymax::Float32 = extrema(r*sinθ*sinϕ for r in rvals for sinθ in sinθs for sinϕ in sinϕs)
> ymax = max(ymax, new_ymax); ymin = min(ymin, new_ymin)
> new_zmin::Float32, new_zmax::Float32 = extrema(r*cosθ for r in rvals for cosθ in cosθs)
> zmax = max(zmax, new_zmax); zmin = min(zmin, new_zmin)
> end
> @show xmin, xmax, ymin, ymax, zmin, zmax
> return xmin, xmax, ymin, ymax, zmin, zmax
> end
> 
> xmin, xmax, ymin, ymax, zmin, zmax = compute_bounds(x1v, x2v, x3v, Nblocks)
> 
> function cartesian_to_spherical(x, y, z)
> r = hypot(x, y, z)
> θ = r == 0 ? 0.0 : acos(z/r)
> ϕ = mod2pi(atan(y, x)) # In [0, 2pi) (like ϕs)
> return r, θ, ϕ
> end
> 
> function append_repeat_last(A::AbstractArray)
> B = similar(A, size(A) .+ 1)
> B[axes(A)...] .= A
> for I in CartesianIndices(A)
> B[I] = A[I]
> end
> for I in CartesianIndices(B)
> idx = ntuple(d -> min(I[d], size(A, d)), ndims(A))
> B[I] = A[idx...]
> end
> return B
> end
> 
> function block_interpolator(x1,x2,x3, Bblock::Array{<:Real,3})
> extrapolate(interpolate((x1, x2, x3), Bblock, Gridded(Linear())), NaN32)
> #extrapolate(interpolate(Bblock, BSpline(Cubic())), NaN32) #,x1, x2, x3)
> end
> 
> function interp(blocks)
> interps = Vector{Any}(undef, blocks)
> @threads for b in Nblocks
> @views begin
> x1 = log.(x1v[:, b]); x2 = x2v[:, b]; x3 = x3v[:, b]
> Bp_block = Bp[:, :, :, b]
> Bp_block = append_repeat_last(Bp_block)
> interps[b] = block_interpolator(x1,x2,x3, Bp_block)
> end
> end
> return interps
> end
> 
> interps = interp(blocks)
> steps::Int32 = 20
> 
> function fill_field(xmin, xmax, ymin, ymax, zmin, zmax, interps, Nblocks, steps)
> @show rx= length(xmin:steps:xmax); ry= length(ymin:steps:ymax); rz = length(zmin:steps:zmax)
> field = fill(NaN32, (rx,ry,rz))
> #=for (i,x) in enumerate(xmin:steps:xmax), (j,y) in enumerate(ymin:steps:ymax), (k,z) in enumerate(zmin:steps:zmax)
> r, θ, ϕ = cartesian_to_spherical(x, y, z)
> for b in Nblocks
> itp = interps[b](log(r), θ, ϕ)
> itp !==NaN32 ? field[i,j,k]= itp : continue
> end
> end=#
> @threads for i in 1:rx
> x = xmin + (i-1)*steps
> for j in 1:ry
> y = ymin + (j-1)*steps
> for k in 1:rz
> z = zmin + (k-1)*steps
> r, θ, ϕ = cartesian_to_spherical(x, y, z)
> lr = log(r)
> @inbounds for b in Nblocks
> v = interps[b](lr, θ, ϕ)
> if !isnan(v)
> field[i,j,k] = v
> break
> end
> end
> end
> end
> end
> 
> return field
> end
> 
> field = fill_field(xmin, xmax, ymin, ymax, zmin, zmax, interps, Nblocks, steps)
> 
> function tdplot(xmin, xmax, ymin, ymax, zmin, zmax, field)
> finite_vals = filter(!isnan, vec(field))
> q_low, q_high = quantile(finite_vals, (0.05, 0.95))
> fig = Figure()
> ax = LScene(fig[1,1], show_axis = false)
> #volume!(ax, xmin..xmax, ymin..ymax, zmin..zmax, field; transparency=true, algorithm = :iso, isorange =0.05, isovalue =median(finite_vals))
> volume!(ax, xmin..xmax, ymin..ymax, zmin..zmax, field; transparency=true,, colorrange=(q_low, q_high))
> display(fig)
> end
> 
> tdplot(xmin, xmax, ymin, ymax, zmin, zmax, field)
> 
> end
> 
> ```

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [January 30, 2026, 3:47pm UTC](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235/67 "2026-01-30T15:47:32Z")

</div>

> [@raman\_kumar](#):
>
> I think appending arrays in `Bp` will have same effect as of `Flat()` in extrapolation. Am i correct?

Appending to the arrays by looking across the boundary would indeed be one way to interpolate across block boundaries.

I see no reason why it would be similar to `Flat()`, nor why you would want that.

* * *

You really should open a new topic, as this has nothing to do with the performance of your `clp`.

* * *

Also, have you tried asking an LLM about your problems? In absence of human responses on this forum, I would expect them to be able to sufficiently help you with your questions on interpolation / extrapolation, and interactions with periodicity, boundaries, etc. Additionally, they’ll ‘gladly’ write tons of code for you. (Of course, it’s still up to you to check and understand what they’re actually writing.)

[Previous page](https://discourse.julialang.org/t/improve-performance-of-this-index-creation-and-accessing-code/134235.md?page=3)