Improve performance of this index creation and accessing code

[abraemer]

@raman_kumar

function bound(xe_all, ye_all, ze_all, MeshBlockSize, blocks)
	XX,YY,ZZ = ntuple(_ -> zeros(Float32, blocks* prod(MeshBlockSize.+1)), 3)
	icx=1
	for b in 1:blocks
		@views r, θ, ϕ = xe_all[:, b], ye_all[:, b], ze_all[:, b]
		for rr in r, th in θ, ph in ϕ,  
			XX[icx] = rr * sin(th) * cos(ph)
			YY[icx] = rr * sin(th) * sin(ph)
			ZZ[icx] = rr * cos(th)
			icx += 1
		end
	end
	xmin,xmax = extrema(XX); ymin,ymax = extrema(YY); zmin,zmax = extrema(ZZ)
	println("  x ∈ [$xmin, $xmax], y ∈ [$ymin, $ymax], z ∈ [$zmin, $zmax]")
	return xmin, xmax, ymin, ymax, zmin, zmax
end

This function is very inefficient: you allocate a huge array only to get the extrema. Also since you take all combinations anyhow, you can extremize each axis by itself. Here is my suggestion which should take a fraction of the original runtime:

function bound(xe_all, ye_all, ze_all)
    xmin = Inf; xmax = -Inf
    ymin = Inf; ymax = -Inf
    zmin = Inf; zmax = -Inf
    for block in axes(xe_all, 2)
        @views r, theta, phi = xe_all[:, b], ye_all[:, b], ze_all[:, b]
        rvals = extrema(r)
        theta_sin_vals = extrema(sin, theta)
        theta_cos_vals = extrema(cos, theta)
        phi_sin_vals = extrema(sin, phi)
        phi_cos_vals = extrema(cos, phi)
        
        new_zmin, new_zmax = extrema(r*cos_theta for r in rvals for cos_theta in theta_cos_vals)
        zmax = max(zmax, new_zmax)
        zmin = min(zmin, new_zmin)
        # do the others yourself
    end
    return xmin, xmax, ymin, ymax, zmin, zmax
end

@raman_kumar

		@inbounds for ii in 1:nr, jj in 1:nθ, kk in 1:nϕ
			cart(iex,ii,jj,kk,r_eval,θ_eval,ϕ_eval,ext_itp,X,Y,Z,C)
			clp(X[iex],Y[iex],Z[iex],C[iex],xmin,ymin,zmin,sx,sy,sz,nx,ny,nz,field)
			iex += 1
		end

Note that the xmin, ymin, zmin parameters are globals again. In this inner loop that will cost a lot of performance.