adds stat computation for Gaussian field test

TurTLE/TEST.py

@@ -129,6 +129,12 @@ class TEST(_code):
             pars['tracers0_integration_steps'] = int(4)
             pars['tracers0_neighbours'] = int(1)
             pars['tracers0_smoothness'] = int(1)
+        if dns_type == 'Gauss_field_test':
+            pars['histogram_bins'] = 129
+            pars['max_velocity_estimate'] = 8.
+            pars['spectrum_slope'] = 1.
+            pars['spectrum_k_cutoff'] = 16.
+            pars['spectrum_coefficient'] = 1.
         return pars
     def get_kspace(self):
         kspace = {}
@@ -175,6 +181,38 @@ class TEST(_code):
             kspace = self.get_kspace()
             nshells = kspace['nshell'].shape[0]
             ofile['checkpoint'] = int(0)
+            vec_stat_datasets = []
+            scal_stat_datasets = []
+            if self.dns_type in ['Gauss_field_test']:
+                vec_stat_datasets.append('velocity')
+            for k in vec_stat_datasets:
+                time_chunk = 2**20//(8*3*3*nshells)
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/spectra/' + k + '_' + k,
+                                     (1, nshells, 3, 3),
+                                     chunks = (time_chunk, nshells, 3, 3),
+                                     maxshape = (None, nshells, 3, 3),
+                                     dtype = np.float64)
+                time_chunk = 2**20//(8*4*10)
+                time_chunk = max(time_chunk, 1)
+                a = ofile.create_dataset('statistics/moments/' + k,
+                                     (1, 10, 4),
+                                     chunks = (time_chunk, 10, 4),
+                                     maxshape = (None, 10, 4),
+                                     dtype = np.float64)
+                time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/histograms/' + k,
+                                     (1,
+                                      self.parameters['histogram_bins'],
+                                      4),
+                                     chunks = (time_chunk,
+                                               self.parameters['histogram_bins'],
+                                               4),
+                                     maxshape = (None,
+                                                 self.parameters['histogram_bins'],
+                                                 4),
+                                     dtype = np.int64)
         return None
     def job_parser_arguments(
             self,

TurTLE/test/test_Gaussian_field.py

+#! /usr/bin/env python
+
+import numpy as np
+import h5py
+import sys
+
+import TurTLE
+from TurTLE import TEST
+
+try:
+    import matplotlib.pyplot as plt
+except:
+    plt = None
+
+def main():
+    c = TEST()
+    # size of grid
+    n = 64
+    c.launch(
+            ['Gauss_field_test',
+             '--nx', str(n),
+             '--ny', str(n),
+             '--nz', str(n),
+             '--simname', 'Gaussianity_test',
+             '--np', '4',
+             '--ntpp', '1',
+             '--wd', './'] +
+             sys.argv[1:])
+    df = h5py.File(c.simname + '.h5', 'r')
+    f = plt.figure()
+
+    # test spectrum
+    a = f.add_subplot(121)
+    kk = df['kspace/kshell'][...]
+    phi_ij = df['statistics/spectra/velocity_velocity'][0]
+    energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
+    a.plot(kk, energy)
+    a.set_xscale('log')
+    a.set_yscale('log')
+    # test isotropy
+    a = f.add_subplot(122)
+    a.plot(df['statistics/histograms/velocity'][0, :, :3])
+    f.tight_layout()
+    f.savefig('spectrum_isotropy_test.pdf')
+    df.close()
+    return None
+
+if __name__ == '__main__':
+    main()
+

cpp/full_code/Gauss_field_test.cpp

@@ -160,17 +160,50 @@ template <typename rnumber>
 int Gauss_field_test<rnumber>::do_work(void)
 {
     TIMEZONE("Gauss_field_test::do_work");
-    make_gaussian_random_field(this->kk, this->scalar_field);
+
+    /// initialize vector Gaussian random field
     make_gaussian_random_field(this->kk, this->vector_field, 3);
-    DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
+
+    /// impose divergence free condition while maintaining the energy of the field
     this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata(), true);
-    DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
-    this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
-    DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
-    make_gaussian_random_field(this->kk, this->vector_field, 3);
-    DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
-    this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
-    DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
+
+    /// initialize statistics file.
+    hid_t stat_group, stat_file;
+    if (this->myrank == 0)
+    {
+        // set caching parameters
+        hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
+        herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
+        DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
+        stat_file = H5Fopen(
+                (this->simname + ".h5").c_str(),
+                H5F_ACC_RDWR,
+                fapl);
+        stat_group = H5Gopen(
+                stat_file,
+                "statistics",
+                H5P_DEFAULT);
+    }
+    else
+    {
+        stat_file = 0;
+        stat_group = 0;
+    }
+
+    /// compute basic statistics of Gaussian field
+    this->vector_field->compute_stats(
+            this->kk,
+            stat_group,
+            "velocity",
+            0,
+            8.);
+
+    /// close stat file
+    if (this->myrank == 0)
+    {
+        H5Gclose(stat_group);
+        H5Fclose(stat_file);
+    }
     return EXIT_SUCCESS;
 }