test_Gaussian_field.py 4.28 KB
Newer Older
1
2
3
#! /usr/bin/env python

import numpy as np
4
5
from scipy import trapz
from scipy.stats import norm
Lukas Bentkamp's avatar
Lukas Bentkamp committed
6
from scipy.integrate import quad
7
import h5py
8
import sys, os
9
import time
10
11
12
13
14
15
16
17
18
19
20

import TurTLE
from TurTLE import TEST

try:
    import matplotlib.pyplot as plt
except:
    plt = None

def main():
    # size of grid
21
    n = 1024
22
    slope = -5./3.
23
    k_cutoff = 64.
24
    func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
25
    total_energy = quad(func, 0.6, k_cutoff*8)[0]
26
    coeff = 1./total_energy
Lukas Bentkamp's avatar
Lukas Bentkamp committed
27
    bin_no = 100
28
    rseed = int(time.time())
29
    simname = 'Gaussianity_test'
30

31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
    if not os.path.exists(simname+'.h5'):
        c = TEST()
        opt = c.launch(
            ['Gauss_field_test',
             '--nx', str(n),
             '--ny', str(n),
             '--nz', str(n),
             '--simname', simname,
             '--np', '4',
             '--environment', 'short',
             '--minutes', '60',
             '--ntpp', '1',
             '--wd', './',
             '--histogram_bins', str(bin_no),
             '--max_velocity_estimate', '8.',
             '--spectrum_slope', str(slope),
             '--spectrum_k_cutoff', str(k_cutoff),
             '--spectrum_coefficient', str(coeff),
             '--field_random_seed', str(rseed)] +
             sys.argv[1:])
    plot_stuff(simname, total_energy = total_energy)
52
    return None
53

54
def plot_stuff(simname, total_energy = 1.):
Lukas Bentkamp's avatar
Lukas Bentkamp committed
55
    df = h5py.File(simname + '.h5', 'r')
Cristian Lalescu's avatar
Cristian Lalescu committed
56
57
58
59
60
61
    for kk in ['spectrum_slope',
               'spectrum_k_cutoff',
               'spectrum_coefficient',
               'field_random_seed',
               'histogram_bins']:
        print(kk, df['parameters/' + kk][...])
Lukas Bentkamp's avatar
Lukas Bentkamp committed
62
63
64
65
    coeff = df['parameters/spectrum_coefficient'][()]
    k_cutoff = df['parameters/spectrum_k_cutoff'][()]
    slope = df['parameters/spectrum_slope'][()]
    bin_no = df['parameters/histogram_bins'][()]
66

67
68
69
70
    f = plt.figure()
    # test spectrum
    a = f.add_subplot(121)
    kk = df['kspace/kshell'][...]
71
72
    print('dk: {}'.format(df['kspace/dk'][()]))
    phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()]
73
    energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
74

75
    a.scatter(kk[1:-2], energy[1:-2])
Lukas Bentkamp's avatar
Lukas Bentkamp committed
76
    a.plot(kk[1:-2], coeff*kk[1:-2]**slope*np.exp(-kk[1:-2]/k_cutoff), ls='--', c='C0')
77
78
79
80
    a.set_xscale('log')
    a.set_yscale('log')
    # test isotropy
    a = f.add_subplot(122)
81
82
83

    max_vel_estimate = df['parameters/max_velocity_estimate'][()]
    velbinsize = 2*max_vel_estimate/bin_no
84
85
86
87
    vel = np.linspace(-max_vel_estimate+velbinsize/2, max_vel_estimate-velbinsize/2, bin_no)
    hist_vel = df['statistics/histograms/velocity'][0, :, :3]
    f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize

Lukas Bentkamp's avatar
Lukas Bentkamp committed
88
    print('Energy analytically: {}'.format(total_energy))
89
    print(np.sum(energy*np.arange(len(energy))**2))
Lukas Bentkamp's avatar
Lukas Bentkamp committed
90
    print('Energy sum: {}'.format(np.sum(energy*df['kspace/dk'][()])))
91
92
93
94
95
96
97
98
99
    print('Moment sum: {}'.format(df['statistics/moments/velocity'][0,2,3]/2))
    print('Velocity variances: {}'.format(trapz(vel[:,None]**2*f_vel, vel[:,None], axis=0)))

    vel_variance = df['statistics/moments/velocity'][0,2,3]/3.
    a.plot(vel[:,None]/np.sqrt(vel_variance), f_vel*np.sqrt(vel_variance))
    a.plot(vel/np.sqrt(vel_variance), norm.pdf(vel/np.sqrt(vel_variance)), ls='--')
    a.set_yscale('log')
    a.set_xlim(-5,5)
    a.set_ylim(1e-5,1)
100
101
    f.tight_layout()
    f.savefig('spectrum_isotropy_test.pdf')
102
103
104
105
106
107
108
    plt.close(f)

    ### check divergence
    print('Divergence second moment is: {0}'.format(
            df['statistics/moments/velocity_divergence'][0, 2]))
    print('Gradient second moment is: {0}'.format(
            df['statistics/moments/velocity_gradient'][0, 2].mean()))
109
110
111
112

    print('----------- k2-premultiplied spectrum -----------')
    k2func = lambda k, k_c=k_cutoff, s=slope : k**(2+s)*np.exp(-k/k_c)
    k2sum_analytic = quad(k2func, 0, k_cutoff*20)[0]
113
    print('Analytically: {}'.format(k2sum_analytic))
114
115
116
117
    k2spec_trace = (
            df['statistics/spectra/k*velocity_k*velocity'][..., 0, 0]
            + df['statistics/spectra/k*velocity_k*velocity'][..., 1, 1]
            + df['statistics/spectra/k*velocity_k*velocity'][..., 2, 2])
118
    print('Energy sum: {}'.format(np.sum(k2spec_trace*df['kspace/dk'][()])/2./coeff))
119

120
121
122
123
124
125
    df.close()
    return None

if __name__ == '__main__':
    main()