From 2d38e4042f90d9a37626fea017ad281794df3c2b Mon Sep 17 00:00:00 2001 From: Lukas Bentkamp <bentkamp@lashgarak.lmp.ds.mpg.de> Date: Tue, 5 Nov 2019 13:25:52 +0100 Subject: [PATCH] More Gaussian field test stuff --- TurTLE/test/test_Gaussian_field.py | 41 ++++++++++++++++++++++-------- cpp/full_code/Gauss_field_test.cpp | 7 ++--- 2 files changed, 34 insertions(+), 14 deletions(-) diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py index f28abf74..7f4c53b1 100644 --- a/TurTLE/test/test_Gaussian_field.py +++ b/TurTLE/test/test_Gaussian_field.py @@ -1,9 +1,11 @@ #! /usr/bin/env python import numpy as np -from scipy.special import gamma +from scipy import trapz +from scipy.stats import norm import h5py import sys +import time import TurTLE from TurTLE import TEST @@ -16,11 +18,12 @@ except: def main(): c = TEST() # size of grid - n = 64 - slope = 1. - k_cutoff = 20. - coeff = np.sqrt(1./0.93058) - bin_no = 33 + n = 128 + slope = -5./3. + k_cutoff = 16. + coeff = 1. + bin_no = 129 + rseed = int(time.time()) c.launch( ['Gauss_field_test', @@ -36,20 +39,23 @@ def main(): '--spectrum_slope', str(slope), '--spectrum_k_cutoff', str(k_cutoff), '--spectrum_coefficient', str(coeff), - '--field_random_seed', '33'] + + '--field_random_seed', str(rseed)] + sys.argv[1:]) df = h5py.File(c.simname + '.h5', 'r') f = plt.figure() + # EVERYTHING SHOULD BE READ FROM FILE BECAUSE WE CANT BE SURE THE CODE RAN AGAIN + # test spectrum a = f.add_subplot(121) kk = df['kspace/kshell'][...] - phi_ij = df['statistics/spectra/velocity_velocity'][0] + print('dk: {}'.format(df['kspace/dk'][()])) + phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()] energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2 - a.plot(kk[1:-2], energy[1:-2]) + a.scatter(kk[1:-2], energy[1:-2]) a.plot(kk[1:-2], coeff*kk[1:-2]**slope*np.exp(-kk[1:-2]/k_cutoff), ls='--', c='C0') a.set_xscale('log') a.set_yscale('log') @@ -58,8 +64,21 @@ def main(): max_vel_estimate = df['parameters/max_velocity_estimate'][()] velbinsize = 2*max_vel_estimate/bin_no - np.linspace(-max_vel_estimate+velbinsize/2, max_vel_estimate-velbinsize/2, bin_no) - a.plot(df['statistics/histograms/velocity'][0, :, :3]) + 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 + + print('Energy integral: {}'.format(trapz(energy[1:-2], kk[1:-2]))) + print('Energy sum: {}'.format(np.sum(energy[1:-2]))) + 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) f.tight_layout() f.savefig('spectrum_isotropy_test.pdf') df.close() diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp index 7e26954d..2b6926da 100644 --- a/cpp/full_code/Gauss_field_test.cpp +++ b/cpp/full_code/Gauss_field_test.cpp @@ -82,9 +82,10 @@ int make_gaussian_random_field( } case THREE: for (int cc = 0; cc<3; cc++) - { - output_field->cval(cindex,cc,0) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4. - 1) * exp(- sqrt(k2)/k_cutoff/2.); - output_field->cval(cindex,cc,1) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4. - 1) * exp(- sqrt(k2)/k_cutoff/2.); + { // In the slope, a factor of 1/2 compensates the k^2, a factor of 1/2 goes from the energy to the Fourier coefficient + // and the 2*pi*k^2 divides out the surface of the shell + output_field->cval(cindex,cc,0) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, (slope - 2)/4.) * exp(- sqrt(k2)/k_cutoff/2.); + output_field->cval(cindex,cc,1) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, (slope - 2)/4.) * exp(- sqrt(k2)/k_cutoff/2.); } // TODO: ADJUST OTHER FIELD TYPES break; case THREExTHREE: -- GitLab