Building Minimal Cosmological Zoomin ICs

Comments

Source

Cubes, Ellipsoids, and Convex Hulls are non-optimal

When building a cosmological zoom-in simulation, the strategy is to run a low resolution, dark matter only simulation, select halos of interest, trace the particles that form that halo back to the initial conditions, and then build a new set of ICs with the regions that form the halo given higher resolution. In order to minimize the number of high-resolution elements needed in the zoom-in IC, the volume that is refined should hug the particles of interest as tightly as possible. Unfortunately, because the cosmic web is composed of sheets and filaments, the regions in an IC that need refinement can have complex shapes (they often look like prawn crackers). This means that simple shapes that enclose the region (cubes, ellipsoids and convex hulls are often used) will frequently contain many times the volume just traced by the particles, giving much larger (and computationally expensive). THERE HAS TO BE ANOTHER WAY!

A Grid Solution

A dead simple, efficient solution would be to simply grid the data (with half the spacing of the low resolution simulation, to prevent oversampling and generating holes inside the high-resolution), and find grid cells that contain particles. In order to ensure that resolution changes smoothly, once the high-resolution region is found, you can sweep down to the original resolution, tagging neighbour cells of high-resolution regions until you finally reach the original resolution. That way, every cell is neighbouring a cell that has no more than 1 level (a factor of 2) change in resolution.

In [159]:
cd ~/ssd/IC_Slice/

/home/kellerbw/ssd/IC_Slice
In [160]:
res = 128 #grid resolution
reflevel = 6 #How much extra refinement to use (results in an effective grid resolution of res*2^reflevel)
#Load particle data that defines zoom-in region
particles = genfromtxt('3Rvir.dat')

#Build grid
grid = zeros((res,res,res))
In [161]:
#Place the particles on the grid
for i,j,k in particles:
    x = int((i+0.5)*res)
    y = int((j+0.5)*res)
    z = int((k+0.5)*res)
    grid[x,y,z] += 1
In [162]:
idx = where(grid > 0)
In [163]:
grid[idx] = reflevel
In [164]:
#Go through the levels of refinement all the way back to level 0
import itertools
nn = (0,1,-1)
for curlevel in range(reflevel, 1, -1):
    refidx = where(grid == curlevel)
    for i,j,k in transpose(array(refidx)):
        offsets = itertools.product(nn,nn,nn)
        for x,y,z in offsets:
            if grid[i+x,j+y,k+z] == 0:
                grid[i+x,j+y,k+z] = curlevel-1
In [177]:
#visualize it
from tempfile import NamedTemporaryFile
from mpl_toolkits.mplot3d import Axes3D
from IPython.display import HTML
import matplotlib.animation as animation

VIDEO_TAG = """<video controls>
 <source src="data:video/x-webm;base64,{0}" type="video/webm">
 Your browser does not support the video tag.
</video>"""

def anim_to_html(anim):
    if not hasattr(anim, '_encoded_video'):
        with NamedTemporaryFile(suffix='.webm') as f:
            anim.save(f.name, fps=30, bitrate=1000, extra_args=['-vcodec', 'libvpx'])
            video = open(f.name, "rb").read()
        anim._encoded_video = video.encode("base64")
    
    return VIDEO_TAG.format(anim._encoded_video)


def display_animation(anim):
    plt.close(anim._fig)
    return HTML(anim_to_html(anim))

plotidx = nonzero(grid)
fig = figure()
ax = fig.add_subplot(111, projection='3d')
def initplot():
    ax.scatter(plotidx[0]/float(res)-0.5, plotidx[1]/float(res)-0.5, plotidx[2]/float(res)-0.5, c=grid[plotidx], marker=',', alpha=0.1)
    fig.set_size_inches(10,10)
def animate(i):
    ax.view_init(azim=i)
    ax.set_axis_off()
anim = animation.FuncAnimation(fig, animate, init_func=initplot,
                               frames=360, interval=10, blit=True)
display_animation(anim)
Out[177]:

Voila, a region that contains everything that will fall into $3R_{vir}$ of a halo, transitioning in steps of 2 from an effective resolution of $128^3$ to $4096^3$.

Comments

Comments powered by Disqus