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

updated 3d examples

parent f9dfcf8c
...@@ -41,7 +41,7 @@ FOF_SECONDARY_LINK_TYPES=1+16+32 # 2^type for the types linked to neares ...@@ -41,7 +41,7 @@ FOF_SECONDARY_LINK_TYPES=1+16+32 # 2^type for the types linked to neares
SUBFIND # enables substructure finder SUBFIND # enables substructure finder
SAVE_HSML_IN_SNAPSHOT # stores hsml, density, and velocity dispersion values in the snapshot files SAVE_HSML_IN_SNAPSHOT # stores hsml, density, and velocity dispersion values in the snapshot files
SUBFIND_CALC_MORE # calculates also the velocity dispersion in the local density estimate (this is automatically enabled by several other options, e.g. SAVE_HSML_IN_SNAPSHOT) SUBFIND_CALC_MORE # calculates also the velocity dispersion in the local density estimate (this is automatically enabled by several other options, e.g. SAVE_HSML_IN_SNAPSHOT)
SUBFIND_EXTENDED_PROPERTIES # adds calculation of further quantities related to angular momentum in different components #SUBFIND_EXTENDED_PROPERTIES # adds calculation of further quantities related to angular momentum in different components
#-------------------------------------------- Things for special behaviour #-------------------------------------------- Things for special behaviour
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) 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)
......
...@@ -8,10 +8,12 @@ created by Rainer Weinberger, last modified 25.02.2019 ...@@ -8,10 +8,12 @@ created by Rainer Weinberger, last modified 25.02.2019
import sys ## load sys; needed for exit codes import sys ## load sys; needed for exit codes
import numpy as np ## load numpy import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt ## plot stuff import os # file specific calls
import matplotlib.pyplot as plt # needs to be active for plotting!
plt.rcParams['text.usetex'] = True
simulation_directory = str(sys.argv[1]) simulation_directory = str(sys.argv[1])
print("examples/cosmo_box_gravity_only_3d/check.py: checking simulation output in directory " + simulation_directory) print("cosmo_box_gravity_only_3d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32 FloatType = np.float64 # double precision: np.float64, for single use np.float32
...@@ -24,7 +26,12 @@ Redshifts = [1, 0] ...@@ -24,7 +26,12 @@ Redshifts = [1, 0]
status = 0 status = 0
CompareAgainstReferenceRun = True ## comparison for small L50m32 box; deactivate this when comparing against self-created ICs CompareAgainstReferenceRun = True ## comparison for small L50m32 box; deactivate this when comparing against self-created ICs
createFigures = False makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
for i_file, z in enumerate(Redshifts): for i_file, z in enumerate(Redshifts):
""" try to read in snapshot """ """ try to read in snapshot """
...@@ -46,7 +53,7 @@ for i_file, z in enumerate(Redshifts): ...@@ -46,7 +53,7 @@ for i_file, z in enumerate(Redshifts):
if CompareAgainstReferenceRun: if CompareAgainstReferenceRun:
## comparison to reference run (sorted list of M200) ## comparison to reference run (sorted list of M200)
M200c_ref = np.loadtxt("./examples/cosmo_box_gravity_only_3d/Masses_L50n32_z%.1d.txt"% z) M200c_ref = np.loadtxt(simulation_directory+"/Masses_L50n32_z%.1d.txt"% z)
minLen = np.min([len(M200c), len(M200c_ref)]) minLen = np.min([len(M200c), len(M200c_ref)])
i_select = np.arange(minLen) i_select = np.arange(minLen)
...@@ -67,7 +74,7 @@ for i_file, z in enumerate(Redshifts): ...@@ -67,7 +74,7 @@ for i_file, z in enumerate(Redshifts):
print(np.std(delta)) print(np.std(delta))
""" optional figure """ """ optional figure """
if createFigures: if makeplots:
filename = "snap_%03d.hdf5" % (i_file) filename = "snap_%03d.hdf5" % (i_file)
try: try:
data = h5py.File(directory+filename, "r") data = h5py.File(directory+filename, "r")
...@@ -98,9 +105,12 @@ for i_file, z in enumerate(Redshifts): ...@@ -98,9 +105,12 @@ for i_file, z in enumerate(Redshifts):
bx.set_yscale('log') bx.set_yscale('log')
bx.set_xlim([9e12,5e14]) bx.set_xlim([9e12,5e14])
bx.set_ylim([9e-7,2e-4]) bx.set_ylim([9e-7,2e-4])
bx.set_xlabel(r"M$_{200,c}$ [M$_\odot$]") bx.set_xlabel(r"$M_{200,c}\ \mathrm{[M_\odot]}$")
bx.set_ylabel(r"$n$(>M) [Mpc$^{-3}$]") bx.set_ylabel(r"$n(>M)\ \mathrm{[Mpc^{-3}]}$")
fig.savefig(simulation_directory+'/largeScaleStructure_z%.1d.pdf'%z, dpi=300)
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
fig.savefig(simulation_directory+'/plots/largeScaleStructure_z%.1d.pdf'%z, dpi=300)
## if everything is ok: 0 else: 1 ## if everything is ok: 0 else: 1
sys.exit(status) sys.exit(status)
...@@ -20,13 +20,13 @@ ic_creation='copy' ...@@ -20,13 +20,13 @@ ic_creation='copy'
""" input """ """ input """
simulation_directory = str(sys.argv[1]) simulation_directory = str(sys.argv[1])
print("examples/cosmo_box_gravity_only_3d/create.py " + simulation_directory) print("create.py " + simulation_directory)
""" initial conditions: either copy or create with code """ """ initial conditions: either copy or create with code """
if ic_creation == 'copy': if ic_creation == 'copy':
## copy and use provided initial conditions ## copy and use provided initial conditions
call(['cp', './examples/cosmo_box_gravity_only_3d/L50n32/ics', simulation_directory+'/ics']) call(['cp', simulation_directory+'/L50n32/ics', simulation_directory+'/ics'])
elif ic_creation == 'music': elif ic_creation == 'music':
## create new initial conditions with the MUSIC code ## create new initial conditions with the MUSIC code
...@@ -40,7 +40,7 @@ elif ic_creation == 'music': ...@@ -40,7 +40,7 @@ elif ic_creation == 'music':
if status != 0: if status != 0:
print('CREATE: ERROR: make failed!') print('CREATE: ERROR: make failed!')
sys.exit(status) sys.exit(status)
status = call(['./MUSIC',cwd+'/examples/cosmo_box_gravity_only_3d/param_music.txt']) status = call(['./MUSIC',cwd+'/param_music.txt'])
if status != 0: if status != 0:
print('CREATE: ERROR: execution failed!') print('CREATE: ERROR: execution failed!')
sys.exit(status) sys.exit(status)
...@@ -58,7 +58,7 @@ elif ic_creation == 'ngenic': ...@@ -58,7 +58,7 @@ elif ic_creation == 'ngenic':
if status != 0: if status != 0:
print('CREATE: ERROR: make failed!') print('CREATE: ERROR: make failed!')
sys.exit(status) sys.exit(status)
status = call(['mpiexec','-np','1','./N-GenIC',cwd+'/examples/cosmo_box_gravity_only_3d/param_ngenic.txt']) status = call(['mpiexec','-np','1','./N-GenIC',cwd+'/param_ngenic.txt'])
if status != 0: if status != 0:
print('CREATE: ERROR: execution failed!') print('CREATE: ERROR: execution failed!')
sys.exit(status) sys.exit(status)
......
%---- Relevant files %---- Relevant files
InitCondFile ics InitCondFile ./ics
OutputDir output OutputDir ./output
SnapshotFileBase snap SnapshotFileBase snap
OutputListFilename output_list.txt OutputListFilename ./output_list.txt
%---- File formats %---- File formats
ICFormat 1 ICFormat 1
......
...@@ -50,7 +50,7 @@ FOF_SECONDARY_LINK_TYPES=1+16+32 # 2^type for the types linked t ...@@ -50,7 +50,7 @@ FOF_SECONDARY_LINK_TYPES=1+16+32 # 2^type for the types linked t
SUBFIND # enables substructure finder SUBFIND # enables substructure finder
SAVE_HSML_IN_SNAPSHOT # stores hsml, density, and velocity dispersion values in the snapshot files SAVE_HSML_IN_SNAPSHOT # stores hsml, density, and velocity dispersion values in the snapshot files
SUBFIND_CALC_MORE # calculates also the velocity dispersion in the local density estimate (this is automatically enabled by several other options, e.g. SAVE_HSML_IN_SNAPSHOT) SUBFIND_CALC_MORE # calculates also the velocity dispersion in the local density estimate (this is automatically enabled by several other options, e.g. SAVE_HSML_IN_SNAPSHOT)
SUBFIND_EXTENDED_PROPERTIES # adds calculation of further quantities related to angular momentum in different components #SUBFIND_EXTENDED_PROPERTIES # adds calculation of further quantities related to angular momentum in different components
#------------------------------------------------ Things for special behaviour #------------------------------------------------ Things for special behaviour
......
...@@ -8,16 +8,24 @@ created by Rainer Weinberger, last modified 28.02.2019 ...@@ -8,16 +8,24 @@ created by Rainer Weinberger, last modified 28.02.2019
import sys ## system calls import sys ## system calls
import numpy as np ## load numpy import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt ## plot stuff import os # file specific calls
import matplotlib.pyplot as plt # needs to be active for plotting!
plt.rcParams['text.usetex'] = True
simulation_directory = str(sys.argv[1]) simulation_directory = str(sys.argv[1])
print("examples/cosmo_box_star_formation_3d/check.py: checking simulation output in directory " + simulation_directory) print("cosmo_box_star_formation_3d: checking simulation output in directory " + simulation_directory)
""" script settings """ """ script settings """
createReferenceSolution = False createReferenceSolution = False
createFigures = False
CompareToReferenceRun = True CompareToReferenceRun = True
makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
FloatType = np.float64 # double precision: np.float64, for single use np.float32 FloatType = np.float64 # double precision: np.float64, for single use np.float32
Boxsize = 7.5 ##Mpc/h Boxsize = 7.5 ##Mpc/h
...@@ -35,7 +43,7 @@ filename = "sfr.txt" ...@@ -35,7 +43,7 @@ filename = "sfr.txt"
data = np.loadtxt(directory+filename) data = np.loadtxt(directory+filename)
if createReferenceSolution: if createReferenceSolution:
np.savetxt("./examples/cosmo_box_star_formation_3d/sfr_ref_L8n32_reduced.txt", data[::50,:]) np.savetxt(simulation_directory+"/sfr_ref_L8n32_reduced.txt", data[::50,:])
scalefactor = data[:,0] scalefactor = data[:,0]
redshift_plus_one = 1.0 / data[:,0] redshift_plus_one = 1.0 / data[:,0]
...@@ -43,7 +51,7 @@ sfrDensity = data[:,3] / Volume ...@@ -43,7 +51,7 @@ sfrDensity = data[:,3] / Volume
cum_mass_stars = data[:,5] * UnitMass / Volume cum_mass_stars = data[:,5] * UnitMass / Volume
""" compare to reference output """ """ compare to reference output """
data_ref = np.loadtxt("./examples/cosmo_box_star_formation_3d/sfr_ref_L8n32_reduced.txt") data_ref = np.loadtxt(simulation_directory+"/sfr_ref_L8n32_reduced.txt")
scalefactor_ref = data_ref[:,0] scalefactor_ref = data_ref[:,0]
redshift_plus_one_ref = 1.0 / data_ref[:,0] redshift_plus_one_ref = 1.0 / data_ref[:,0]
cum_mass_stars_ref = data_ref[:,5]* UnitMass / Volume cum_mass_stars_ref = data_ref[:,5]* UnitMass / Volume
...@@ -52,7 +60,10 @@ scalefactor_to_probe = np.array([0.25,0.3333,0.5,0.66667,1]) ...@@ -52,7 +60,10 @@ scalefactor_to_probe = np.array([0.25,0.3333,0.5,0.66667,1])
delta_mass_stars = np.zeros(scalefactor_to_probe.shape) delta_mass_stars = np.zeros(scalefactor_to_probe.shape)
avg_mass_stars = np.zeros(scalefactor_to_probe.shape) avg_mass_stars = np.zeros(scalefactor_to_probe.shape)
if createFigures: if makeplots:
if not os.path.exists( simulation_directory+"/plots" ):
os.mkdir( simulation_directory+"/plots" )
## figure -- dark matter and stellar positions ## figure -- dark matter and stellar positions
fig, ax = plt.subplots( 2, 2, figsize=np.array([6.9,6.9]), sharex=True, sharey=True ) fig, ax = plt.subplots( 2, 2, figsize=np.array([6.9,6.9]), sharex=True, sharey=True )
fig.subplots_adjust(left = 0.09, bottom = 0.09,right = 0.98, top = 0.98, hspace=0.0, wspace=0.0) fig.subplots_adjust(left = 0.09, bottom = 0.09,right = 0.98, top = 0.98, hspace=0.0, wspace=0.0)
...@@ -94,7 +105,71 @@ if createFigures: ...@@ -94,7 +105,71 @@ if createFigures:
ax[1][0].set_xlabel('[Mpc]') ax[1][0].set_xlabel('[Mpc]')
ax[1][0].set_ylabel('[Mpc]') ax[1][0].set_ylabel('[Mpc]')
ax[1][1].set_xlabel('[Mpc]') ax[1][1].set_xlabel('[Mpc]')
fig.savefig(simulation_directory+'/dark_matter_and_stars_evolution.pdf', dpi=300) fig.savefig(simulation_directory+'/plots/dark_matter_and_stars_evolution.pdf', dpi=300)
## figure -- gas density
fig, ax = plt.subplots( 2, 2, figsize=np.array([6.9,6.9]), sharex=True, sharey=True )
fig.subplots_adjust(left = 0.09, bottom = 0.09,right = 0.98, top = 0.98, hspace=0.0, wspace=0.0)
for i, a in enumerate(scalefactor_to_probe):
directory = simulation_directory+"/output/"
filename = "fof_subhalo_tab_%03d.hdf5" % (i)
try:
data = h5py.File(directory+filename, "r")
except:
print("could not open "+directory+filename)
sys.exit(1)
""" get simulation data """
## simulation data
GrpPos = np.array(data["Group"]["GroupPos"], dtype=FloatType) / HubbleParam / 1000.
delta_mass_stars[i] = np.interp(a, scalefactor, cum_mass_stars)
delta_mass_stars[i] -= np.interp(a, scalefactor_ref, cum_mass_stars_ref)
avg_mass_stars[i] = 0.5 * np.interp(a, scalefactor, cum_mass_stars)
avg_mass_stars[i] += 0.5 * np.interp(a, scalefactor_ref, cum_mass_stars_ref)
filename = "snap_%03d.hdf5" % (i)
try:
data = h5py.File(directory+filename, "r")
except:
print("could not open "+directory+filename)
sys.exit(1)
VoronoiPos = np.array(data["PartType0"]["Coordinates"], dtype = FloatType) / HubbleParam / 1000.
Density = np.array(data["PartType0"]["Density"], dtype = FloatType)
if(pos.shape[0] > 32**3):
i_select = np.random.uniform(low=0.0, high=pos.shape[0], size=32**3).astype(np.int)
else:
i_select = np.arange(pos.shape[0])
## Try to show the different gas phases and their spatial distribution
z = 1./ a - 1
if i == 0:
continue
ax_col = np.int((i-1)%2)
ax_row = np.int((i-1)/2)
Nplot = 512
from scipy import spatial # needed for KDTree that we use for nearest neighbour search and Voronoi mesh
Edges1d = np.linspace(0., Boxsize/HubbleParam, Nplot+1, endpoint=True, dtype=FloatType)
Grid1d = 0.5 * (Edges1d[1:] + Edges1d[:-1])
xx, yy = np.meshgrid(Grid1d, Grid1d)
Grid2D = np.array( [xx.reshape(Nplot**2), yy.reshape(Nplot**2), np.ones(Nplot**2)*GrpPos[0,2]] ).T
dist, cells = spatial.KDTree( VoronoiPos ).query( Grid2D, k=1 )
gas_density_slice = Density[cells].reshape((Nplot,Nplot))
ax[ax_row][ax_col].pcolormesh( Edges1d, Edges1d, np.log10(gas_density_slice), rasterized=True, cmap=plt.get_cmap('plasma') )
ax[ax_row][ax_col].text(5.5,10.0,'redshift %.1f'%z, color='w')
ax[ax_row][ax_col].set_xlim([0,Boxsize/HubbleParam])
ax[ax_row][ax_col].set_ylim([0,Boxsize/HubbleParam])
ax[0][0].set_ylabel('[Mpc]')
ax[1][0].set_xlabel('[Mpc]')
ax[1][0].set_ylabel('[Mpc]')
ax[1][1].set_xlabel('[Mpc]')
fig.savefig(simulation_directory+'/plots/gas_density_evolution.pdf', dpi=300)
""" figure of stellar mass density """ """ figure of stellar mass density """
fig=plt.figure() fig=plt.figure()
...@@ -108,20 +183,27 @@ if createFigures: ...@@ -108,20 +183,27 @@ if createFigures:
ax.set_xticklabels(["0", "1", "2", "4", "8"]) ax.set_xticklabels(["0", "1", "2", "4", "8"])
ax.set_xlabel(r"redshift") ax.set_xlabel(r"redshift")
ax.set_ylabel(r"Stellar mass denstiy [M$_\odot$ Mpc$^{-3}$]") ax.set_ylabel(r"Stellar mass denstiy [M$_\odot$ Mpc$^{-3}$]")
fig.savefig(simulation_directory+"/Mstar_redshift.pdf") fig.savefig(simulation_directory+"/plots/Mstar_redshift.pdf")
if CompareToReferenceRun: if CompareToReferenceRun:
""" check if deviations within tolerance """ """ check if deviations within tolerance """
abs_tolerance = 5e7 ## depends on box size... abs_tolerance = 5e7 ## depends on box size...
if np.max( np.abs(delta_mass_stars ) ) > abs_tolerance:
print("Error: stellar mass deviaties from reference more than %g (Msun/Mpc^3): "% abs_tolerance ) for i, a in enumerate(scalefactor_to_probe):
print("redshifts:") delta_mass_stars[i] = np.interp(a, scalefactor, cum_mass_stars)
print(1./scalefactor_to_probe - 1) delta_mass_stars[i] -= np.interp(a, scalefactor_ref, cum_mass_stars_ref)
print("relative deviation:") avg_mass_stars[i] = 0.5 * np.interp(a, scalefactor, cum_mass_stars)
print(delta_mass_stars / avg_mass_stars) avg_mass_stars[i] += 0.5 * np.interp(a, scalefactor_ref, cum_mass_stars_ref)
print("absolute deviation (1e6 Msun/Mpc^3):")
print(delta_mass_stars/1.0e6) if np.max( np.abs(delta_mass_stars ) ) > abs_tolerance:
sys.exit(1) print("Error: stellar mass deviaties from reference more than %g (Msun/Mpc^3): "% abs_tolerance )
print("redshifts:")
print(1./scalefactor_to_probe - 1)
print("relative deviation:")
print(delta_mass_stars / avg_mass_stars)
print("absolute deviation (1e6 Msun/Mpc^3):")
print(delta_mass_stars/1.0e6)
sys.exit(1)
""" if everything is ok """ """ if everything is ok """
sys.exit(0) sys.exit(0)
...@@ -19,12 +19,12 @@ ic_creation='copy' ...@@ -19,12 +19,12 @@ ic_creation='copy'
""" input """ """ input """
simulation_directory = str(sys.argv[1]) simulation_directory = str(sys.argv[1])
print("examples/cosmo_box_star_formation_3d/create.py " + simulation_directory) print("create.py " + simulation_directory)
""" initial conditions: either copy or create with code """ """ initial conditions: either copy or create with code """
if ic_creation == 'copy': if ic_creation == 'copy':
call(['cp', './examples/cosmo_box_star_formation_3d/L8n32/ics', simulation_directory+'/ics']) call(['cp', simulation_directory+'/L8n32/ics', simulation_directory+'/ics'])
elif ic_creation == 'music': elif ic_creation == 'music':
status = call(['hg', 'clone', 'https://bitbucket.org/ohahn/music', simulation_directory+'/music']) status = call(['hg', 'clone', 'https://bitbucket.org/ohahn/music', simulation_directory+'/music'])
if status != 0: if status != 0:
...@@ -36,7 +36,7 @@ elif ic_creation == 'music': ...@@ -36,7 +36,7 @@ elif ic_creation == 'music':
if status != 0: if status != 0:
print('CREATE: ERROR: make failed!') print('CREATE: ERROR: make failed!')
sys.exit(status) sys.exit(status)
status = call(['./MUSIC',cwd+'/examples/cosmo_box_star_formation_3d/param_music.txt']) status = call(['./MUSIC',cwd+'/param_music.txt'])
if status != 0: if status != 0:
print('CREATE: ERROR: execution failed!') print('CREATE: ERROR: execution failed!')
sys.exit(status) sys.exit(status)
...@@ -52,7 +52,7 @@ elif ic_creation == 'ngenic': ...@@ -52,7 +52,7 @@ elif ic_creation == 'ngenic':
if status != 0: if status != 0:
print('CREATE: ERROR: make failed!') print('CREATE: ERROR: make failed!')
sys.exit(status) sys.exit(status)
status = call(['mpiexec','-np','1','./N-GenIC',cwd+'/examples/cosmo_box_star_formation_3d/param_ngenic.txt']) status = call(['mpiexec','-np','1','./N-GenIC',cwd+'/param_ngenic.txt'])
if status != 0: if status != 0:
print('CREATE: ERROR: execution failed!') print('CREATE: ERROR: execution failed!')
sys.exit(status) sys.exit(status)
...@@ -66,15 +66,17 @@ else: ...@@ -66,15 +66,17 @@ else:
outputTimes = np.array([0.2,0.25,0.33,0.5,0.66,1], dtype=np.float64) outputTimes = np.array([0.2,0.25,0.33,0.5,0.66,1], dtype=np.float64)
ones = np.ones(outputTimes.shape, dtype=np.int) ones = np.ones(outputTimes.shape, dtype=np.int)
""" copy treecool file to run directory """
if len(sys.argv) > 2:
arepopath = sys.argv[2]
else:
arepopath = '.'
call(['cp', arepopath + '/data/TREECOOL_ep', simulation_directory+'/TREECOOL_ep'])
""" write output list file """ """ write output list file """
data = np.array([outputTimes, ones]).T data = np.array([outputTimes, ones]).T
np.savetxt(simulation_directory+"/output_list.txt",data, fmt="%g %1.f" ) np.savetxt(simulation_directory+"/output_list.txt",data, fmt="%g %1.f" )
""" copy treecool file to run directory """
call(['cp', './data/TREECOOL_ep', simulation_directory+'/TREECOOL_ep'])
""" normal exit """ """ normal exit """
sys.exit(0) sys.exit(0)
%% parameters for cosmo_box_star_formation_3d %% parameters for cosmo_box_star_formation_3d
%---- Relevant files %---- Relevant files
InitCondFile ics InitCondFile ./ics
OutputDir output OutputDir ./output
SnapshotFileBase snap SnapshotFileBase snap
OutputListFilename output_list.txt OutputListFilename ./output_list.txt
%---- File formats %---- File formats
ICFormat 1 ICFormat 1
...@@ -120,4 +120,4 @@ FactorEVP 573.0 ...@@ -120,4 +120,4 @@ FactorEVP 573.0
TemperatureThresh 1e+06 TemperatureThresh 1e+06
FactorSN 0.1 FactorSN 0.1
TreecoolFile TREECOOL_ep TreecoolFile ./TREECOOL_ep
...@@ -9,13 +9,20 @@ created by Rainer Weinberger, last modified 05.03.2019 ...@@ -9,13 +9,20 @@ created by Rainer Weinberger, last modified 05.03.2019
import sys ## load sys; needed for exit codes import sys ## load sys; needed for exit codes
import numpy as np ## load numpy import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots import h5py ## load h5py; needed to read snapshots
from matplotlib import pyplot as plt ## optional plots import os # file specific calls
import matplotlib.pyplot as plt # needs to be active for plotting!
plt.rcParams['text.usetex'] = True
createReferenceSolution = False createReferenceSolution = False
createFigures = False makeplots = True
if len(sys.argv) > 2:
if sys.argv[2] == "True":
makeplots = True
else:
makeplots = False
simulation_directory = str(sys.argv[1]) simulation_directory = str(sys.argv[1])
print("examples/cosmo_zoom_3d/check.py: checking simulation output in directory " + simulation_directory) print("cosmo_zoom_3d: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32 FloatType = np.float64 # double precision: np.float64, for single use np.float32
IntType = np.int32 IntType = np.int32
...@@ -60,10 +67,10 @@ for i_file, z in enumerate(Redshifts): ...@@ -60,10 +67,10 @@ for i_file, z in enumerate(Redshifts):
if createReferenceSolution: if createReferenceSolution:
# save reference masses # save reference masses
np.savetxt("./examples/cosmo_zoom_gravity_only_3d/Masses_z%.1g.txt"% z, SubhaloMass) np.savetxt(simulation_directory+"/Masses_z%.1g.txt"% z, SubhaloMass)
""" optional figure """ """ optional figure """
if createFigures: if makeplots:
filename = "snap_%03d.hdf5" % (i_file) filename = "snap_%03d.hdf5" % (i_file)
try: try:
data = h5py.File(directory+filename, "r") data = h5py.File(directory+filename, "r")
...@@ -101,22 +108,22 @@ for i_file, z in enumerate(Redshifts): ...@@ -101,22 +108,22 @@ for i_file, z in enumerate(Redshifts):
particle_select = np.where( (pos[:,2]>zpos-2.5) & (pos[:,2]<zpos+2.5) )[0] 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] 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(pos[particle_select, 0], pos[particle_select, 1], marker='.', s=0.05, alpha=0.5)
ax.scatter(pos2[particle2_select, 0],pos2[particle2_select,1], marker='.', s=0.05, c='r', alpha=0.5, rasterized=True) ax.scatter(pos2[particle2_select, 0],pos2[particle2_select,1], marker='.', s=0.05, c='r', alpha=0.5)
if i_file != 0: if i_file != 0:
ax.add_artist(plt.Circle((GrpPos[0,0], GrpPos[0,1]),GrpR200c[0],linestyle='--', color='k', fill=False)) 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]): for i in np.arange(SubhaloRad.shape[0]):
if (dist[i] < 2.0 * GrpR200c[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.add_artist(plt.Circle((SubhaloPos[i,0], SubhaloPos[i,1]), SubhaloRad[i], color='k', fill=False))
ax.set_xlim([0,100./HubbleParam]) ax.set_xlim([0,100./HubbleParam])
ax.set_ylim([0,100./HubbleParam]) ax.set_ylim([0,100./HubbleParam])
if i_file != 0: if i_file != 0:
ax.set_xlim([ GrpPos[0,0]-1.75, GrpPos[0,0]+1.25]) 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_ylim( GrpPos[0,1]-1.25, GrpPos[0,1]+1.75 )
ax.set_xlabel('[Mpc]') ax.set_xlabel('[Mpc]')
ax.set_ylabel('[Mpc]') ax.set_ylabel('[Mpc]')
...@@ -128,19 +135,21 @@ for i_file, z in enumerate(Redshifts): ...@@ -128,19 +135,21 @@ for i_file, z in enumerate(Redshifts):
bx.set_yscale('log') bx.set_yscale('log')
bx.set_xlim([9e9,3e12]) bx.set_xlim([9e9,3e12])