Commit ca624182 authored by Rainer Weinberger's avatar Rainer Weinberger
Browse files

added example cosmo_zoom_gravity_only_3d (outcommented in test.sh, since it...

added example cosmo_zoom_gravity_only_3d (outcommented in test.sh, since it requires external IC code in create.py which might not run on all systems)
parent 68a6f07c
#!/bin/bash # this line only there to enable syntax highlighting in this file
###################################################
## examples/cosmo_zoom_gravity_only_3d/Config.sh ##
###################################################
#--------------------------------------- Gravity treatment
SELFGRAVITY # gravitational intraction between simulation particles/cells
HIERARCHICAL_GRAVITY # use hierarchical splitting of the time integration of the gravity
ALLOW_DIRECT_SUMMATION # Performed direct summation instead of tree-based gravity if number of active particles < DIRECT_SUMMATION_THRESHOLD (= 3000 unless specified differently here)
DIRECT_SUMMATION_THRESHOLD=1024 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation
EVALPOTENTIAL # computes gravitational potential
#--------------------------------------- TreePM Options
PMGRID=256 # Enables particle mesh; number of cells used for grid in each dimension
RCUT=5.5 # This can be used to override the maximum radius in which the short-range tree-force is evaluated (in case the TreePM algorithm is used). The default value is 4.5, given in mesh-cells.
PLACEHIGHRESREGION=2 # Places a second, high-resolution PM grid; number encodes high-res particle types from which the extent of this second PM grid is calculated.
ENLARGEREGION=1.1 # Factor by which high res region is increased with respect to max-min position of high-res particles
GRIDBOOST=2 # Factor by which PMGRID is increased in non-periodic (or high res) PM calculation (if not set, code uses 2)
PM_ZOOM_OPTIMIZED # Particle-mesh calculation that is optimized for cosmological zoom simulations; disable for cosmological volume simulations.
#--------------------------------------- Gravity softening
MULTIPLE_NODE_SOFTENING # If a tree node is to be used which is softened, this is done with the softenings of its different mass components
INDIVIDUAL_GRAVITY_SOFTENING=4+8+16+32 # bitmask with particle types where the softenig type should be chosen with that of parttype 1 as a reference type
#--------------------------------------- Time integration options
TREE_BASED_TIMESTEPS # non-local timestep criterion (take 'signal speed' into account)
#--------------------------------------- Single/Double Precision
DOUBLEPRECISION=1 # Mode of double precision: not defined: single; 1: full double precision 2: mixed, 3: mixed, fewer single precisions; unless short of memory, use 1.
DOUBLEPRECISION_FFTW # FFTW calculation in double precision
#--------------------------------------- On the fly FOF groupfinder
FOF # enable FoF output
FOF_PRIMARY_LINK_TYPES=2 # 2^type for the primary dark matter type
FOF_SECONDARY_LINK_TYPES=1+16+32 # 2^type for the types linked to nearest primaries
#--------------------------------------- Subfind
SUBFIND # enables substructure finder
#--------------------------------------- output fields in snapshots--Default output filds are: position, velocity, ID, mass, spec. internal energy (gas), density(gas)
OUTPUTPOTENTIAL # output potential at particle position
#--------------------------------------- output options
PROCESS_TIMES_OF_OUTPUTLIST # goes through times of output list prior to starting the simulaiton to ensure that outputs are written as close to the desired time as possible (as opposed to at next possible time if this flag is not active)
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps
#!/bin/bash # this line only there to enable syntax highlighting in this file
##########################################################
## examples/cosmo_zoom_gravity_only_3d/Config_parent.sh ##
##########################################################
#--------------------------------------- Gravity treatment
SELFGRAVITY # gravitational intraction between simulation particles/cells
HIERARCHICAL_GRAVITY # use hierarchical splitting of the time integration of the gravity
ALLOW_DIRECT_SUMMATION # Performed direct summation instead of tree-based gravity if number of active particles < DIRECT_SUMMATION_THRESHOLD (= 3000 unless specified differently here)
DIRECT_SUMMATION_THRESHOLD=1024 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation
#--------------------------------------- TreePM Options
PMGRID=256 # Enables particle mesh; number of cells used for grid in each dimension
#--------------------------------------- Gravity softening
MULTIPLE_NODE_SOFTENING # If a tree node is to be used which is softened, this is done with the softenings of its different mass components
#--------------------------------------- Time integration options
TREE_BASED_TIMESTEPS # non-local timestep criterion (take 'signal speed' into account)
#--------------------------------------- Single/Double Precision
DOUBLEPRECISION=1 # Mode of double precision: not defined: single; 1: full double precision 2: mixed, 3: mixed, fewer single precisions; unless short of memory, use 1.
DOUBLEPRECISION_FFTW # FFTW calculation in double precision
#--------------------------------------- On the fly FOF groupfinder
FOF # enable FoF output
FOF_PRIMARY_LINK_TYPES=2 # 2^type for the primary dark matter type
#--------------------------------------- Subfind
SUBFIND # enables substructure finder
#--------------------------------------- output options
PROCESS_TIMES_OF_OUTPUTLIST # goes through times of output list prior to starting the simulaiton to ensure that outputs are written as close to the desired time as possible (as opposed to at next possible time if this flag is not active)
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps
7.051804088933790625e+13
1.851817588029021729e+12
1.123913657112673340e+12
5.412134689182386475e+11
4.950123266868264771e+11
3.837524067169450684e+11
3.733807268978959961e+11
3.535802293253184204e+11
3.234080519336917725e+11
3.196365192008997192e+11
2.762640194798557739e+11
2.687209821711749268e+11
2.206340594949150696e+11
1.980049334904209290e+11
1.895190059593162537e+11
1.876332395929201965e+11
1.866903634489479980e+11
1.631183472220299988e+11
1.536895435469531250e+11
1.404892071390841675e+11
1.329461698304033508e+11
1.310604034640073090e+11
1.272888848096669006e+11
1.159743147681939850e+11
1.037168685827489014e+11
9.805958356201243591e+10
9.145942239730380249e+10
8.863077988693557739e+10
7.071604868075401306e+10
6.883028231435797119e+10
6.411588047681954193e+10
6.128723796645131683e+10
6.034436182247913361e+10
5.657283612891287994e+10
5.562995646532778168e+10
5.468707680174267578e+10
5.468707680174267578e+10
4.902979178100621796e+10
4.620114927063799286e+10
4.525826960705288696e+10
4.431538994346778870e+10
4.431538994346778870e+10
4.431538994346778870e+10
4.242962709668466187e+10
4.148674743309955597e+10
4.054386776951445770e+10
4.054386776951445770e+10
4.054386776951445770e+10
4.054386776951445770e+10
3.771522525914622498e+10
3.677234559556112671e+10
3.677234559556112671e+10
3.488658274877799988e+10
3.394370308519289780e+10
3.300081990199487305e+10
3.300081990199487305e+10
3.111506057482466888e+10
3.111506057482466888e+10
3.111506057482466888e+10
3.017218091123956680e+10
3.017218091123956680e+10
2.922929948784800339e+10
2.922929948784800339e+10
2.828641806445643997e+10
2.828641806445643997e+10
2.734353840087133789e+10
2.734353840087133789e+10
2.734353840087133789e+10
2.640065697747977448e+10
2.640065697747977448e+10
2.640065697747977448e+10
2.545777731389467239e+10
2.545777731389467239e+10
2.451489589050310898e+10
2.451489589050310898e+10
2.357201622691800690e+10
2.357201622691800690e+10
2.357201622691800690e+10
2.357201622691800690e+10
2.357201622691800690e+10
2.262913480352644348e+10
2.168625338013488388e+10
1.980049229315821838e+10
1.980049229315821838e+10
1.885761262957311249e+10
1.064267082203092656e+14
2.644780184056226074e+12
1.144657073064578369e+12
7.891910914513154297e+11
6.977316690540117188e+11
3.865810492273132935e+11
3.507515868149501953e+11
3.450943017942137451e+11
2.847499470109604492e+11
2.423203093554370422e+11
2.055479707991017456e+11
1.923476484696844788e+11
1.734900270410790710e+11
1.725471508971068726e+11
1.536895435469531250e+11
1.338890459743755341e+11
1.254031184432708588e+11
1.159743147681939850e+11
1.037168685827489014e+11
9.805958356201243591e+10
9.711670037881442261e+10
7.825909302866067505e+10
7.825909302866067505e+10
7.731620984546266174e+10
7.260180800792422485e+10
7.165892482472621155e+10
7.165892482472621155e+10
7.071604868075401306e+10
6.600163980398974609e+10
6.411588047681954193e+10
6.317300433284735870e+10
6.223012114964933777e+10
5.657283612891287994e+10
5.562995646532778168e+10
5.374419361854465485e+10
5.280131395495954895e+10
5.280131395495954895e+10
5.280131395495954895e+10
4.997267144459132385e+10
4.997267144459132385e+10
4.902979178100621796e+10
4.620114927063799286e+10
4.337250676026976776e+10
4.242962709668466187e+10
4.148674743309955597e+10
4.148674743309955597e+10
3.960098458631643677e+10
3.771522525914622498e+10
3.771522525914622498e+10
3.677234559556112671e+10
3.677234559556112671e+10
3.582946241236310577e+10
3.488658274877799988e+10
3.394370308519289780e+10
3.394370308519289780e+10
3.300081990199487305e+10
3.300081990199487305e+10
3.205794023840977097e+10
2.922929948784800339e+10
2.828641806445643997e+10
2.828641806445643997e+10
2.734353840087133789e+10
2.734353840087133789e+10
2.734353840087133789e+10
2.734353840087133789e+10
2.640065697747977448e+10
2.545777731389467239e+10
2.545777731389467239e+10
2.451489589050310898e+10
2.451489589050310898e+10
2.451489589050310898e+10
2.451489589050310898e+10
2.357201622691800690e+10
2.262913480352644348e+10
2.168625338013488388e+10
2.168625338013488388e+10
2.168625338013488388e+10
2.168625338013488388e+10
2.168625338013488388e+10
2.168625338013488388e+10
2.074337371654977798e+10
1.980049229315821838e+10
1.980049229315821838e+10
1.980049229315821838e+10
1.885761262957311249e+10
1.885761262957311249e+10
4.187615763166150000e+13
1.683041955786255371e+12
1.614211617158090576e+12
7.005603115643798828e+11
6.515305268225994873e+11
5.242416138560292358e+11
3.828095446514245605e+11
2.724924867470636902e+11
2.611779448624941101e+11
2.366630243347005920e+11
2.281771108820475769e+11
2.140338983302064209e+11
1.404892071390841675e+11
1.188029572785622253e+11
1.140885554410237885e+11
1.027739853995508728e+11
9.145942239730380249e+10
9.051653921410577393e+10
7.260180800792422485e+10
6.883028231435797119e+10
6.883028231435797119e+10
6.883028231435797119e+10
6.317300433284735870e+10
5.940147863928111267e+10
5.845859897569600677e+10
5.374419361854465485e+10
5.091555462778934479e+10
4.902979178100621796e+10
4.902979178100621796e+10
4.714403245383601379e+10
4.431538994346778870e+10
4.337250676026976776e+10
4.337250676026976776e+10
4.242962709668466187e+10
4.148674743309955597e+10
4.054386776951445770e+10
3.960098458631643677e+10
3.865810492273133087e+10
3.771522525914622498e+10
3.582946241236310577e+10
3.488658274877799988e+10
3.300081990199487305e+10
3.205794023840977097e+10
3.111506057482466888e+10
3.017218091123956680e+10
2.922929948784800339e+10
2.828641806445643997e+10
2.734353840087133789e+10
2.734353840087133789e+10
2.640065697747977448e+10
2.640065697747977448e+10
2.545777731389467239e+10
2.545777731389467239e+10
2.357201622691800690e+10
2.357201622691800690e+10
2.262913480352644348e+10
2.262913480352644348e+10
2.262913480352644348e+10
2.262913480352644348e+10
1.885761262957311249e+10
1.499840212952326953e+13
2.762640194798557739e+11
1.527466533245292664e+11
1.414320973615080261e+11
1.301175273200351257e+11
1.301175273200351257e+11
1.216315997889304504e+11
8.863077988693557739e+10
8.768789670373754883e+10
8.674501352053953552e+10
7.920196917263287354e+10
6.505876366001756287e+10
5.845859897569600677e+10
5.657283612891287994e+10
4.431538994346778870e+10
3.771522525914622498e+10
3.582946241236310577e+10
3.582946241236310577e+10
3.300081990199487305e+10
3.017218091123956680e+10
2.828641806445643997e+10
2.640065697747977448e+10
2.262913480352644348e+10
1.885761262957311249e+10
7.730678291422166016e+12
3.252938042216361694e+11
8.957365603090776062e+10
8.108773553902891541e+10
8.014485235583088684e+10
6.600163980398974609e+10
4.242962709668466187e+10
3.111506057482466888e+10
2.757925884470955078e+12
2.489204845985973511e+11
3.677234559556112671e+10
2.545777731389467239e+10
1.885761262957311249e+10
""" @package ./examples/cosmo_zoom_3d/check.py
Code that checks results of gravity only cosmological zoom simulation
created by Rainer Weinberger, last modified 05.03.2019
"""
""" load libraries """
import sys ## load sys; needed for exit codes
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
from matplotlib import pyplot as plt ## optional plots
createReferenceSolution = False
createFigures = False
simulation_directory = str(sys.argv[1])
print("examples/cosmo_zoom_3d/check.py: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32
Redshifts = np.array([50,4,3,2,1,0.5,0])
HubbleParam = 0.6774
UnitMass = 1.0e10
status = 0
for i_file, z in enumerate(Redshifts):
""" try to read in snapshot """
directory = simulation_directory+"/output/"
filename = "fof_subhalo_tab_%03d.hdf5" % (i_file)
try:
data = h5py.File(directory+filename, "r")
except:
print("could not open "+directory+filename)
sys.exit(1)
""" get simulation data """
zpos = 62.
if i_file != 0:
SubhaloMass = np.array(data["Subhalo"]["SubhaloMass"], dtype=FloatType) * UnitMass / HubbleParam
SubhaloPos = np.array(data["Subhalo"]["SubhaloPos"], dtype=FloatType)/ HubbleParam / 1000.
SubhaloRad = np.array(data["Subhalo"]["SubhaloHalfmassRad"], dtype=FloatType)/ HubbleParam / 1000.
GrpPos = np.array(data["Group"]["GroupPos"], dtype=FloatType)/ HubbleParam / 1000.
GrpR200c = np.array(data["Group"]["Group_R_Crit200"], dtype=FloatType)/ HubbleParam / 1000.
zpos = GrpPos[0,2]
## select subhalos within main halo only (2*R200c)
dist = (SubhaloPos[:,0] - GrpPos[0,0])**2
dist += (SubhaloPos[:,1] - GrpPos[0,1])**2
dist += (SubhaloPos[:,2] - GrpPos[0,2])**2
dist = np.sqrt(dist)
i_select, = np.where(dist < 2.0 * GrpR200c[0])
n_subs = len(SubhaloMass[i_select])
SubhaloMass = np.sort(SubhaloMass[i_select])[::-1]
if createReferenceSolution:
# save reference masses
np.savetxt("./examples/cosmo_zoom_gravity_only_3d/Masses_z%.1g.txt"% z, SubhaloMass)
""" optional figure """
if createFigures:
filename = "snap_%03d.hdf5" % (i_file)
try:
data = h5py.File(directory+filename, "r")
except:
print("could not open "+directory+filename)
sys.exit(1)
pos = np.array(data["PartType1"]["Coordinates"], dtype=FloatType) / HubbleParam / 1000.
pos2 = np.array(data["PartType2"]["Coordinates"], dtype=FloatType) / HubbleParam / 1000.
vel2 = np.array(data["PartType2"]["Velocities"], dtype=FloatType)
mass2 = np.array(data["PartType2"]["Masses"], dtype=FloatType) / HubbleParam * UnitMass
id2 = np.array(data["PartType2"]["ParticleIDs"], dtype=IntType)
if i_file != 0:
cont_dist = (pos2[:,0] - GrpPos[0,0])**2
cont_dist += (pos2[:,1] - GrpPos[0,1])**2
cont_dist += (pos2[:,2] - GrpPos[0,2])**2
cont_dist = np.sqrt(cont_dist)
print('minimum distance of contaminating (low-res) particle: %g'%np.min(cont_dist) )
i_issue = np.where(cont_dist < GrpR200c[0])[0]
if len(i_issue) > 0:
filename = "snap_%03d.hdf5" % (0)
data = h5py.File(directory+filename, "r")
id2_ic = np.array(data["PartType2"]["ParticleIDs"], dtype=IntType)
pos2_ic = np.array(data["PartType2"]["Coordinates"], dtype=FloatType)
ids, i_select1, dummy = np.intersect1d(id2_ic, id2[i_issue], return_indices=True)
print('problem at positions (code units!)')
print(pos2_ic[i_select1,:])
print(id2_ic[i_select1])
exit(1)
fig = plt.figure(figsize=(6.9,6.9))
ax = plt.axes([0.1,0.1,0.87,0.87])
particle_select = np.where( (pos[:,2]>zpos-2.5) & (pos[:,2]<zpos+2.5) )[0]
particle2_select = np.where( (pos2[:,2]>zpos-2.5) & (pos2[:,2]<zpos+2.5) )[0]
ax.scatter(pos[particle_select, 0], pos[particle_select, 1], marker='.', s=0.05, alpha=0.5, rasterized=True)
ax.scatter(pos2[particle2_select, 0],pos2[particle2_select,1], marker='.', s=0.05, c='r', alpha=0.5, rasterized=True)
if i_file != 0:
ax.add_artist(plt.Circle((GrpPos[0,0], GrpPos[0,1]),GrpR200c[0],linestyle='--', color='k', fill=False))
for i in np.arange(SubhaloRad.shape[0]):
if (dist[i] < 2.0 * GrpR200c[0]):
ax.add_artist(plt.Circle((SubhaloPos[i,0], SubhaloPos[i,1]), SubhaloRad[i], color='k', fill=False))
ax.set_xlim([0,100./HubbleParam])
ax.set_ylim([0,100./HubbleParam])
if i_file != 0:
ax.set_xlim([ GrpPos[0,0]-1.75, GrpPos[0,0]+1.25])
ax.set_ylim([ GrpPos[0,1]-1.25, GrpPos[0,1]+1.75])
ax.set_xlabel('[Mpc]')
ax.set_ylabel('[Mpc]')
if i_file != 0:
CumMassFunction = np.cumsum( np.ones(SubhaloMass.shape[0]))
bx = plt.axes([0.20,0.74,0.26,0.22])
bx.plot(SubhaloMass, CumMassFunction)
bx.set_xscale('log')
bx.set_yscale('log')
bx.set_xlim([9e9,3e12])
bx.set_ylim([0.9,150])
bx.set_xlabel(r"$M_{h}$ [M$_\odot$]")
bx.set_ylabel(r"$N$(>$M$)")
fig.savefig(simulation_directory+'/haloStructure_z%.1d.pdf'%z, dpi=300)
## comparison to reference run (sorted list of SubhaloMass)
if i_file != 0:
SubhaloMass_ref = np.loadtxt("./examples/cosmo_zoom_gravity_only_3d/Masses_z%.1g.txt"% z)
minLen = np.min([len(SubhaloMass), len(SubhaloMass_ref)])
i_select2 = np.arange(minLen)
delta = np.array(SubhaloMass[i_select2]-SubhaloMass_ref[i_select2], dtype=np.float64)
tolerance_average = 1e11
tolerance_std = 1e12
if np.abs(np.average(delta)) > tolerance_average or np.abs(np.std(delta)) > tolerance_std:
status = 1
print("ERROR: z=%g difference in subhalo masses exceeding limits!" % z)
print("mass error (=delta)")
print(delta)
print("average (tolerance: %g)" % tolerance_average)
print(np.average(delta))
print("stddev (tolerance: %g)" % tolerance_std)
print(np.std(delta))
## if everything is ok: 0 else: 1
sys.exit(status)
""" @package examples/cosmo_zoom_3d/create.py
Code that creates ics and outuput list for a cosmological zoom
simulaiton;
created by Rainer Weinberger, last modified 05.03.2019
"""
""" load libraries """
import sys # system calls
import numpy as np # scientific computing package
import h5py # hdf5 format
import os # operating system interface
from subprocess import call # execute bash commands
"""
The usual procedure of producing zoom ICs is to run a uniform
box fist, select the halo in there and trace back the particles
to the initlal condition to define a Lagrangian region there
which will be refined in the zoom ICs. This is done in this
script if the flag 'runParent' is active. If it is inactive,
the zoom simulation is produced right away from a previously
selected meaningful region.
Note that in Arepo, the high-resolution region must not
include a (periodic) boundary at any time of the simulation.
"""
runParent=False
simulation_directory = str(sys.argv[1])
print("examples/cosmo_zoom_gravity_only_3d/create.py " + simulation_directory)
""" output times """
outputTimes = np.array([0.0197,0.2,0.25,0.33,0.5,0.66,1], dtype=np.float64)
ones = np.ones(outputTimes.shape, dtype=np.int)
data = np.array([outputTimes, ones]).T
np.savetxt(simulation_directory+"/output_list.txt",data, fmt="%g %1.f" )
""" IC generating code MUSIC """
# get code from repository and compile IC generating code
status = call(['hg', 'clone', 'https://bitbucket.org/ohahn/music', simulation_directory+'/music'])
if status != 0:
print('CREATE: ERROR: hg clone failed!')
sys.exit(status)
cwd = os.getcwd()
os.chdir(simulation_directory+'/music/')
status = call(['make'])
if status != 0:
print('CREATE: ERROR: make failed!')
sys.exit(status)
os.chdir(cwd)
if runParent:
""" parent simulation """
# execute ic generating code for parent box ICs (change working directory to do this)
os.chdir(simulation_directory+'/music/')
status = call(['./MUSIC', cwd+'/examples/cosmo_zoom_gravity_only_3d/param_music_parent.txt'])
if status != 0:
print('CREATE: ERROR: IC creation failed!')
sys.exit(status)
os.chdir(cwd)