check.py 4.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
""" @package ./examples/cosmo_box_gravity_only_3d/check.py
Code that checks results of gravity only structure formation simulation

created by Rainer Weinberger, last modified 25.02.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
Rainer Weinberger's avatar
Rainer Weinberger committed
11
12
13
import os      # file specific calls
import matplotlib.pyplot as plt    # needs to be active for plotting!
plt.rcParams['text.usetex'] = True
14
15

simulation_directory = str(sys.argv[1])
Rainer Weinberger's avatar
Rainer Weinberger committed
16
print("cosmo_box_gravity_only_3d: checking simulation output in directory " + simulation_directory) 
17
18
19
20
21
22
23
24
25
26
27
28

FloatType = np.float64  # double precision: np.float64, for single use np.float32

Boxsize = 50 ##Mpc/h
HubbleParam = 0.6774
UnitMass = 1.0e10
Volume = Boxsize * Boxsize * Boxsize / HubbleParam / HubbleParam / HubbleParam

Redshifts = [1, 0]
status = 0

CompareAgainstReferenceRun = True ## comparison for small L50m32 box; deactivate this when comparing against self-created ICs
Rainer Weinberger's avatar
Rainer Weinberger committed
29
30
31
32
33
34
makeplots = True
if len(sys.argv) > 2:
  if sys.argv[2] == "True":
    makeplots = True
  else:
    makeplots = False
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

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 """
    ## simulation data
    GrpPos = np.array(data["Group"]["GroupPos"], dtype=FloatType) / HubbleParam / 1000.
    GrpR200c = np.array(data["Group"]["Group_R_Crit200"], dtype=FloatType) / HubbleParam / 1000.
    M200c = np.array(data["Group"]["GroupMass"], dtype=FloatType) * UnitMass
    M200c = np.sort(M200c)[::-1]
    CumMassFunction = np.cumsum(np.ones(M200c.shape) ) / Volume
    
    if CompareAgainstReferenceRun:
        ## comparison to reference run (sorted list of M200)
Rainer Weinberger's avatar
Rainer Weinberger committed
56
        M200c_ref = np.loadtxt(simulation_directory+"/Masses_L50n32_z%.1d.txt"% z)
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
        
        minLen = np.min([len(M200c), len(M200c_ref)])
        i_select = np.arange(minLen)
        
        delta = (M200c[i_select]-M200c_ref[i_select]) / M200c_ref[i_select]
        
        ## empirically based tolerances
        tolerance_average = 0.01
        tolerance_std = 0.05
        if np.abs(np.average(delta)) > tolerance_average or np.abs(np.std(delta)) > tolerance_std:
            status = 1
            print("ERROR: z=%g difference in halo masses exceeding limits!" % z)
            print("relative mass error (=delta)")
            print(delta)
            print("average delta (tolerance: %g)" % tolerance_average)
            print(np.average(delta))
            print("stddev delta (tolerance: %g)" % tolerance_std)
            print(np.std(delta))
    
    """ optional figure """
Rainer Weinberger's avatar
Rainer Weinberger committed
77
    if makeplots:
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
        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.
    
        fig = plt.figure(figsize=(6.9,6.9))
        ax = plt.axes([0.1,0.1,0.87,0.87])
        
        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])
        ax.scatter(pos[i_select, 0], pos[i_select, 1], marker='.', s=0.05, alpha=0.5, rasterized=True)
        for i in np.arange(GrpR200c.shape[0]):
            ax.add_artist(plt.Circle((GrpPos[i,0], GrpPos[i,1]), GrpR200c[i], color='k', fill=False))
        
        ax.set_xlim([0,Boxsize/HubbleParam])
        ax.set_ylim([0,Boxsize/HubbleParam])
        ax.set_xlabel('[Mpc]')
        ax.set_ylabel('[Mpc]')
        
        bx = plt.axes([0.70,0.74,0.26,0.22])
        bx.plot(M200c, CumMassFunction)
        bx.set_xscale('log')
        bx.set_yscale('log')
        bx.set_xlim([9e12,5e14])
        bx.set_ylim([9e-7,2e-4])
Rainer Weinberger's avatar
Rainer Weinberger committed
108
109
110
111
112
113
        bx.set_xlabel(r"$M_{200,c}\ \mathrm{[M_\odot]}$")
        bx.set_ylabel(r"$n(>M)\ \mathrm{[Mpc^{-3}]}$")
        
        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)
114
115
116

## if everything is ok: 0 else: 1
sys.exit(status)