Improve performance of this index creation and accessing code

[raman_kumar]

@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.

@show first(re),last(re),first(θe),last(θe),first(ϕe),last(ϕe)

[eldee]

@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.

[raman_kumar]

@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 describes how blocks are positioned in mesh using LogicalLocations by Z-order curve. Levels denotes the relative size of blocks.

Levels = read(h5["Levels"])
LogicalLocations = read(h5["LogicalLocations"])

@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.

[raman_kumar]

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.

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)

[raman_kumar]

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. I think appending arrays in Bp will have same effect as of Flat() in extrapolation. Am i correct?

code appending array

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

[eldee]

@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.)