plot_gravity_checks.py 12.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#!/usr/bin/env python

import sys
import glob
import re
import numpy as np
import matplotlib.pyplot as plt

params = {'axes.labelsize': 14,
'axes.titlesize': 18,
'font.size': 12,
'legend.fontsize': 12,
'xtick.labelsize': 14,
'ytick.labelsize': 14,
'text.usetex': True,
16
'figure.figsize': (12, 10),
17
18
19
'figure.subplot.left'    : 0.06,
'figure.subplot.right'   : 0.99  ,
'figure.subplot.bottom'  : 0.06  ,
20
'figure.subplot.top'     : 0.99  ,
21
22
23
24
25
26
27
28
29
'figure.subplot.wspace'  : 0.14  ,
'figure.subplot.hspace'  : 0.14  ,
'lines.markersize' : 6,
'lines.linewidth' : 3.,
'text.latex.unicode': True
}
plt.rcParams.update(params)
plt.rc('font',**{'family':'sans-serif','sans-serif':['Times']})

30
31
32
min_error = 1e-7
max_error = 3e-1
num_bins = 64
33
34
35
36
37
38
39
40
41

# Construct the bins
bin_edges = np.linspace(np.log10(min_error), np.log10(max_error), num_bins + 1)
bin_size = (np.log10(max_error) - np.log10(min_error)) / num_bins
bins = 0.5*(bin_edges[1:] + bin_edges[:-1])
bin_edges = 10**bin_edges
bins = 10**bins

# Colours
42
cols = ['#332288', '#88CCEE', '#117733', '#DDCC77', '#CC6677']
43
44
45

# Time-step to plot
step = int(sys.argv[1])
46
periodic = int(sys.argv[2])
47
48

# Find the files for the different expansion orders
49
order_list = glob.glob("gravity_checks_swift_step%.4d_order*.dat"%step)
50
51
52
53
54
num_order = len(order_list)

# Get the multipole orders
order = np.zeros(num_order)
for i in range(num_order):
55
    order[i] = int(order_list[i][35])
56
57
order = sorted(order)
order_list = sorted(order_list)
58
59

# Read the exact accelerations first
60
if periodic:
61
    data = np.loadtxt('gravity_checks_exact_periodic_step%.4d.dat'%step)
62
else:
63
    data = np.loadtxt('gravity_checks_exact_step%.4d.dat'%step)
64
65
66
exact_ids = data[:,0]
exact_pos = data[:,1:4]
exact_a = data[:,4:7]
67
exact_pot = data[:,7]
68
69
70
71
72
# Sort stuff
sort_index = np.argsort(exact_ids)
exact_ids = exact_ids[sort_index]
exact_pos = exact_pos[sort_index, :]
exact_a = exact_a[sort_index, :]        
73
exact_pot = exact_pot[sort_index]
74
exact_a_norm = np.sqrt(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2)
75
76

print "Number of particles tested:", np.size(exact_ids)
77
78
79
80
    
# Start the plot
plt.figure()

81
82
count = 0

83
# Get the Gadget-2 data if existing
84
85
86
87
if periodic:
    gadget2_file_list = glob.glob("forcetest_gadget2_periodic.txt")
else:
    gadget2_file_list = glob.glob("forcetest_gadget2.txt")
88
89
90
91
92
93
94
95
96
97
98
99
100
if len(gadget2_file_list) != 0:

    gadget2_data = np.loadtxt(gadget2_file_list[0])
    gadget2_ids = gadget2_data[:,0]
    gadget2_pos = gadget2_data[:,1:4]
    gadget2_a_exact = gadget2_data[:,4:7]
    gadget2_a_grav = gadget2_data[:, 7:10]

    # Sort stuff
    sort_index = np.argsort(gadget2_ids)
    gadget2_ids = gadget2_ids[sort_index]
    gadget2_pos = gadget2_pos[sort_index, :]
    gadget2_a_exact = gadget2_a_exact[sort_index, :]
101
    gadget2_exact_a_norm = np.sqrt(gadget2_a_exact[:,0]**2 + gadget2_a_exact[:,1]**2 + gadget2_a_exact[:,2]**2)
102
    gadget2_a_grav = gadget2_a_grav[sort_index, :]
103
104
105
106
107
108
109
110
111
112

    # Cross-checks
    if not np.array_equal(exact_ids, gadget2_ids):
        print "Comparing different IDs !"

    if np.max(np.abs(exact_pos - gadget2_pos)/np.abs(gadget2_pos)) > 1e-6:
        print "Comparing different positions ! max difference:"
        index = np.argmax(exact_pos[:,0]**2 + exact_pos[:,1]**2 + exact_pos[:,2]**2 - gadget2_pos[:,0]**2 - gadget2_pos[:,1]**2 - gadget2_pos[:,2]**2)
        print "Gadget2 (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "exact (id=%d):"%exact_ids[index], exact_pos[index,:], "\n"

113
114
115
116
117
118
119
120
    diff = np.abs(exact_a_norm - gadget2_exact_a_norm) / np.abs(gadget2_exact_a_norm)
    max_diff = np.max(diff)
    if max_diff > 2e-6:
        print "Comparing different exact accelerations !"
        print "Median=", np.median(diff), "Mean=", np.mean(diff), "99%=", np.percentile(diff, 99)
        print "max difference ( relative diff =", max_diff, "):"
        #index = np.argmax(exact_a[:,0]**2 + exact_a[:,1]**2 + exact_a[:,2]**2 - gadget2_a_exact[:,0]**2 - gadget2_a_exact[:,1]**2 - gadget2_a_exact[:,2]**2)
        index = np.argmax(diff)
121
        print "a_exact --- Gadget2:", gadget2_a_exact[index,:], "exact:", exact_a[index,:]
122
        print "pos ---     Gadget2: (id=%d):"%gadget2_ids[index], gadget2_pos[index,:], "exact (id=%d):"%gadget2_ids[index], gadget2_pos[index,:],"\n"
123

124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
    
    # Compute the error norm
    diff = gadget2_a_exact - gadget2_a_grav

    norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)
    norm_a = np.sqrt(gadget2_a_exact[:,0]**2 + gadget2_a_exact[:,1]**2 + gadget2_a_exact[:,2]**2)

    norm_error = norm_diff / norm_a
    error_x = abs(diff[:,0]) / norm_a
    error_y = abs(diff[:,1]) / norm_a
    error_z = abs(diff[:,2]) / norm_a
    
    # Bin the error
    norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
    error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
    error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
    error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
    
    norm_median = np.median(norm_error)
    median_x = np.median(error_x)
    median_y = np.median(error_y)
    median_z = np.median(error_z)

147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    norm_per99 = np.percentile(norm_error,99)
    per99_x = np.percentile(error_x,99)
    per99_y = np.percentile(error_y,99)
    per99_z = np.percentile(error_z,99)

    norm_max = np.max(norm_error)
    max_x = np.max(error_x)
    max_y = np.max(error_y)
    max_z = np.max(error_z)

    print "Gadget-2 ---- "
    print "Norm: median= %f 99%%= %f max= %f"%(norm_median, norm_per99, norm_max)
    print "X   : median= %f 99%%= %f max= %f"%(median_x, per99_x, max_x)
    print "Y   : median= %f 99%%= %f max= %f"%(median_y, per99_y, max_y)
    print "Z   : median= %f 99%%= %f max= %f"%(median_z, per99_z, max_z)
    print ""
163

164
    plt.subplot(231)    
165
166
    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", alpha=0.8)
    plt.semilogx(bins, norm_error_hist, 'k--', label="Gadget-2", alpha=0.8)
167
    plt.subplot(232)
168
169
    plt.semilogx(bins, error_x_hist, 'k--', label="Gadget-2", alpha=0.8)
    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_x, per99_x), ha="left", va="top", alpha=0.8)
170
    plt.subplot(233)    
171
172
    plt.semilogx(bins, error_y_hist, 'k--', label="Gadget-2", alpha=0.8)
    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_y, per99_y), ha="left", va="top", alpha=0.8)
173
    plt.subplot(234)    
174
175
176
177
    plt.semilogx(bins, error_z_hist, 'k--', label="Gadget-2", alpha=0.8)
    plt.text(min_error * 1.5, 1.55, "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_z, per99_z), ha="left", va="top", alpha=0.8)
    
    count += 1
178
179
180


# Plot the different histograms
181
for i in range(num_order):
182
183
184
    data = np.loadtxt(order_list[i])
    ids = data[:,0]
    pos = data[:,1:4]
185
    a_grav = data[:, 4:7]
186
    pot = data[:, 7]
187
188
189
190
191

    # Sort stuff
    sort_index = np.argsort(ids)
    ids = ids[sort_index]
    pos = pos[sort_index, :]
192
    a_grav = a_grav[sort_index, :]        
193
    pot = pot[sort_index]
194
195

    # Cross-checks
196
    if not np.array_equal(exact_ids, ids):
197
198
        print "Comparing different IDs !"

199
    if np.max(np.abs(exact_pos - pos)/np.abs(pos)) > 1e-6:
200
        print "Comparing different positions ! max difference:"
201
202
203
        index = np.argmax(exact_pos[:,0]**2 + exact_pos[:,1]**2 + exact_pos[:,2]**2 - pos[:,0]**2 - pos[:,1]**2 - pos[:,2]**2)
        print "SWIFT (id=%d):"%ids[index], pos[index,:], "exact (id=%d):"%exact_ids[index], exact_pos[index,:], "\n"
    
204
    # Compute the error norm
205
    diff = exact_a - a_grav
206
    diff_pot = exact_pot - pot
207

208
209
210
211
212
213
214
    # Correct for different normalization of potential
    print "Difference in normalization of potential:", np.mean(diff_pot),
    print "std_dev=", np.std(diff_pot), "99-percentile:", np.percentile(diff_pot, 99)-np.median(diff_pot), "1-percentile:", np.median(diff_pot) - np.percentile(diff_pot, 1)

    exact_pot -= np.mean(diff_pot)
    diff_pot = exact_pot - pot

215
216
    norm_diff = np.sqrt(diff[:,0]**2 + diff[:,1]**2 + diff[:,2]**2)

217
218
219
220
    norm_error = norm_diff / exact_a_norm
    error_x = abs(diff[:,0]) / exact_a_norm
    error_y = abs(diff[:,1]) / exact_a_norm
    error_z = abs(diff[:,2]) / exact_a_norm
221
    error_pot = abs(diff_pot) / abs(exact_pot)
222
223
224
225
226
227
    
    # Bin the error
    norm_error_hist,_ = np.histogram(norm_error, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
    error_x_hist,_ = np.histogram(error_x, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
    error_y_hist,_ = np.histogram(error_y, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
    error_z_hist,_ = np.histogram(error_z, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
228
    error_pot_hist,_ = np.histogram(error_pot, bins=bin_edges, density=False) / (np.size(norm_error) * bin_size)
229
230
231
232
233

    norm_median = np.median(norm_error)
    median_x = np.median(error_x)
    median_y = np.median(error_y)
    median_z = np.median(error_z)
234
    median_pot = np.median(error_pot)
235

236
237
238
239
    norm_per99 = np.percentile(norm_error,99)
    per99_x = np.percentile(error_x,99)
    per99_y = np.percentile(error_y,99)
    per99_z = np.percentile(error_z,99)
240
    per99_pot = np.percentile(error_pot, 99)
241
242
243
244
245

    norm_max = np.max(norm_error)
    max_x = np.max(error_x)
    max_y = np.max(error_y)
    max_z = np.max(error_z)
246
    max_pot = np.max(error_pot)
247
248
249
250
251
252

    print "Order %d ---- "%order[i]
    print "Norm: median= %f 99%%= %f max= %f"%(norm_median, norm_per99, norm_max)
    print "X   : median= %f 99%%= %f max= %f"%(median_x, per99_x, max_x)
    print "Y   : median= %f 99%%= %f max= %f"%(median_y, per99_y, max_y)
    print "Z   : median= %f 99%%= %f max= %f"%(median_z, per99_z, max_z)
253
    print "Pot : median= %f 99%%= %f max= %f"%(median_pot, per99_pot, max_pot)
254
    print ""
255
    
256
    plt.subplot(231)    
257
    plt.semilogx(bins, error_x_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
258
    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_x, per99_x), ha="left", va="top", color=cols[i])
259
    plt.subplot(232)    
260
    plt.semilogx(bins, error_y_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
261
    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_y, per99_y), ha="left", va="top", color=cols[i])
262
    plt.subplot(233)    
263
    plt.semilogx(bins, error_z_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
264
    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_z, per99_z), ha="left", va="top", color=cols[i])
265
266
267
268
269
270
    plt.subplot(234)
    plt.semilogx(bins, norm_error_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(norm_median, norm_per99), ha="left", va="top", color=cols[i])
    plt.subplot(235)    
    plt.semilogx(bins, error_pot_hist, color=cols[i],label="SWIFT m-poles order %d"%order[i])
    plt.text(min_error * 1.5, 1.5 - count/10., "$50\\%%\\rightarrow%.4f~~ 99\\%%\\rightarrow%.4f$"%(median_pot, per99_pot), ha="left", va="top", color=cols[i])
271
272

    count += 1
273

274
plt.subplot(231)    
275
plt.xlabel("$\delta a_x/|\overrightarrow{a}_{exact}|$")
276
277
278
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
279
#plt.legend(loc="center left")
280
plt.subplot(232)    
281
plt.xlabel("$\delta a_y/|\overrightarrow{a}_{exact}|$")
282
283
284
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
285
#plt.legend(loc="center left")
286
plt.subplot(233)    
287
plt.xlabel("$\delta a_z/|\overrightarrow{a}_{exact}|$")
288
289
290
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
291
292
293
294
295
296
297
298
299
300
301
plt.subplot(234)
plt.xlabel("$|\delta \overrightarrow{a}|/|\overrightarrow{a}_{exact}|$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,2.5)
plt.legend(loc="upper left")
plt.subplot(235)    
plt.xlabel("$\delta \phi/\phi_{exact}$")
#plt.ylabel("Density")
plt.xlim(min_error, max_error)
plt.ylim(0,1.75)
302
303
304
305
#plt.legend(loc="center left")



306
307
plt.savefig("gravity_checks_step%.4d.png"%step, dpi=200)
plt.savefig("gravity_checks_step%.4d.pdf"%step, dpi=200)