{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Cubes, Ellipsoids, and Convex Hulls are non-optimal"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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]( https://www.deliverance.co.uk/Images/MenuItemProduct/211/4066)). 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!](https://youtu.be/bR3S690EF2U?t=10s) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# A Grid Solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 159,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/kellerbw/ssd/IC_Slice\n"
]
}
],
"source": [
"cd ~/ssd/IC_Slice/"
]
},
{
"cell_type": "code",
"execution_count": 160,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"res = 128 #grid resolution\n",
"reflevel = 6 #How much extra refinement to use (results in an effective grid resolution of res*2^reflevel)\n",
"#Load particle data that defines zoom-in region\n",
"particles = genfromtxt('3Rvir.dat')\n",
"\n",
"#Build grid\n",
"grid = zeros((res,res,res))"
]
},
{
"cell_type": "code",
"execution_count": 161,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Place the particles on the grid\n",
"for i,j,k in particles:\n",
" x = int((i+0.5)*res)\n",
" y = int((j+0.5)*res)\n",
" z = int((k+0.5)*res)\n",
" grid[x,y,z] += 1"
]
},
{
"cell_type": "code",
"execution_count": 162,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"idx = where(grid > 0)"
]
},
{
"cell_type": "code",
"execution_count": 163,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"grid[idx] = reflevel"
]
},
{
"cell_type": "code",
"execution_count": 164,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Go through the levels of refinement all the way back to level 0\n",
"import itertools\n",
"nn = (0,1,-1)\n",
"for curlevel in range(reflevel, 1, -1):\n",
" refidx = where(grid == curlevel)\n",
" for i,j,k in transpose(array(refidx)):\n",
" offsets = itertools.product(nn,nn,nn)\n",
" for x,y,z in offsets:\n",
" if grid[i+x,j+y,k+z] == 0:\n",
" grid[i+x,j+y,k+z] = curlevel-1"
]
},
{
"cell_type": "code",
"execution_count": 177,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 177,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#visualize it\n",
"from tempfile import NamedTemporaryFile\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from IPython.display import HTML\n",
"import matplotlib.animation as animation\n",
"\n",
"VIDEO_TAG = \"\"\"\"\"\"\n",
"\n",
"def anim_to_html(anim):\n",
" if not hasattr(anim, '_encoded_video'):\n",
" with NamedTemporaryFile(suffix='.webm') as f:\n",
" anim.save(f.name, fps=30, bitrate=1000, extra_args=['-vcodec', 'libvpx'])\n",
" video = open(f.name, \"rb\").read()\n",
" anim._encoded_video = video.encode(\"base64\")\n",
" \n",
" return VIDEO_TAG.format(anim._encoded_video)\n",
"\n",
"\n",
"def display_animation(anim):\n",
" plt.close(anim._fig)\n",
" return HTML(anim_to_html(anim))\n",
"\n",
"plotidx = nonzero(grid)\n",
"fig = figure()\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"def initplot():\n",
" 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)\n",
" fig.set_size_inches(10,10)\n",
"def animate(i):\n",
" ax.view_init(azim=i)\n",
" ax.set_axis_off()\n",
"anim = animation.FuncAnimation(fig, animate, init_func=initplot,\n",
" frames=360, interval=10, blit=True)\n",
"display_animation(anim)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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$."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"nikola": {
"category": "research",
"date": "2016-02-29 16:40:49 UTC-05:00",
"description": "",
"link": "",
"slug": "building-minimal-cosmological-zoomin-ics",
"tags": "guides,linux,simulations,mathjax",
"title": "Building Minimal Cosmological Zoomin ICs",
"type": "text"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.10+"
}
},
"nbformat": 4,
"nbformat_minor": 0
}