Testing the premultiplied spectrum

TurTLE/TEST.py

@@ -188,6 +188,7 @@ class TEST(_code):
             scal_rspace_stats = []
             if self.dns_type in ['Gauss_field_test']:
                 vec_spectra_stats.append('velocity')
+                vec_spectra_stats.append('k*velocity')
                 vec4_rspace_stats.append('velocity')
                 tens_rspace_stats.append('velocity_gradient')
                 scal_rspace_stats.append('velocity_divergence')

TurTLE/test/test_Gaussian_field.py

@@ -19,31 +19,32 @@ except:
 def main():
     c = TEST()
     # size of grid
-    n = 256
+    n = 1024
     slope = -5./3.
-    k_cutoff = 30.
+    k_cutoff = 64.
     func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
-    total_energy = quad(func, 1, k_cutoff*4)[0]
+    total_energy = quad(func, 0.6, k_cutoff*8)[0]
     coeff = 1./total_energy
     bin_no = 100
     rseed = int(time.time())
+    simname = 'Gaussianity_test'
 
-    c.launch(
-            ['Gauss_field_test',
-             '--nx', str(n),
-             '--ny', str(n),
-             '--nz', str(n),
-             '--simname', 'Gaussianity_test',
-             '--np', '4',
-             '--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:])
+    opt = c.launch(
+        ['Gauss_field_test',
+         '--nx', str(n),
+         '--ny', str(n),
+         '--nz', str(n),
+         '--simname', simname,
+         '--np', '4',
+         '--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(c.simname, total_energy = total_energy)
     return None
 
@@ -82,6 +83,7 @@ def plot_stuff(simname, total_energy = 1.):
     f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize
 
     print('Energy analytically: {}'.format(total_energy))
+    print(np.sum(energy*np.arange(len(energy))**2))
     print('Energy sum: {}'.format(np.sum(energy*df['kspace/dk'][()])))
     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)))
@@ -101,6 +103,17 @@ def plot_stuff(simname, total_energy = 1.):
             df['statistics/moments/velocity_divergence'][0, 2]))
     print('Gradient second moment is: {0}'.format(
             df['statistics/moments/velocity_gradient'][0, 2].mean()))
+
+    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]
+    print('Analytically: {}'.format(k2sum_analytic))
+    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])
+    print('Energy sum: {}'.format(np.sum(k2spec_trace*df['kspace/dk'][()])))
+
     df.close()
     return None
 

cpp/full_code/Gauss_field_test.cpp

@@ -144,6 +144,14 @@ int Gauss_field_test<rnumber>::do_work(void)
             this->vector_field,
             this->scalar_field);
 
+    /// compute integral of premultiplied energy spectrum
+    this->kk->cospectrum<rnumber, THREE>(
+        this->vector_field->get_cdata(),
+        stat_group,
+        "k*velocity_k*velocity",
+        0,
+        2.);
+
     /// compute basic statistics of Gaussian field
     this->vector_field->compute_stats(
             this->kk,

cpp/kspace.cpp

@@ -783,6 +783,7 @@ void kspace<be, dt>::cospectrum(
             {
                 double* spec_local = spec_local_thread.getMine();
                 int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
+                
                 for (hsize_t i=0; i<ncomp(fc); i++)
                 for (hsize_t j=0; j<ncomp(fc); j++){
                     spec_local[tmp_int + i*ncomp(fc)+j] += pow(k2, wavenumber_exp/2.)*