test whether project_divfree actually works

TurTLE/TEST.py

@@ -51,6 +51,7 @@ class TEST(_code):
                 work_dir = work_dir,
                 simname = simname)
         self.generate_default_parameters()
+        self.check_current_vorticity_exists = False
         return None
     def set_precision(
             self,
@@ -181,11 +182,16 @@ class TEST(_code):
             kspace = self.get_kspace()
             nshells = kspace['nshell'].shape[0]
             ofile['checkpoint'] = int(0)
-            vec_stat_datasets = []
-            scal_stat_datasets = []
+            vec_spectra_stats = []
+            tens_rspace_stats = []
+            vec4_rspace_stats = []
+            scal_rspace_stats = []
             if self.dns_type in ['Gauss_field_test']:
-                vec_stat_datasets.append('velocity')
-            for k in vec_stat_datasets:
+                vec_spectra_stats.append('velocity')
+                vec4_rspace_stats.append('velocity')
+                tens_rspace_stats.append('velocity_gradient')
+                scal_rspace_stats.append('velocity_divergence')
+            for k in vec_spectra_stats:
                 time_chunk = 2**20//(8*3*3*nshells)
                 time_chunk = max(time_chunk, 1)
                 ofile.create_dataset('statistics/spectra/' + k + '_' + k,
@@ -193,6 +199,25 @@ class TEST(_code):
                                      chunks = (time_chunk, nshells, 3, 3),
                                      maxshape = (None, nshells, 3, 3),
                                      dtype = np.float64)
+            for k in scal_rspace_stats:
+                time_chunk = 2**20//(8*10)
+                time_chunk = max(time_chunk, 1)
+                a = ofile.create_dataset('statistics/moments/' + k,
+                                     (1, 10),
+                                     chunks = (time_chunk, 10),
+                                     maxshape = (None, 10),
+                                     dtype = np.float64)
+                time_chunk = 2**20//(8*self.parameters['histogram_bins'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/histograms/' + k,
+                                     (1,
+                                      self.parameters['histogram_bins']),
+                                     chunks = (time_chunk,
+                                               self.parameters['histogram_bins']),
+                                     maxshape = (None,
+                                                 self.parameters['histogram_bins']),
+                                     dtype = np.int64)
+            for k in vec4_rspace_stats:
                 time_chunk = 2**20//(8*4*10)
                 time_chunk = max(time_chunk, 1)
                 a = ofile.create_dataset('statistics/moments/' + k,
@@ -213,6 +238,27 @@ class TEST(_code):
                                                  self.parameters['histogram_bins'],
                                                  4),
                                      dtype = np.int64)
+            for k in tens_rspace_stats:
+                time_chunk = 2**20//(8*9*10)
+                time_chunk = max(time_chunk, 1)
+                a = ofile.create_dataset('statistics/moments/' + k,
+                                     (1, 10, 3, 3),
+                                     chunks = (time_chunk, 10, 3, 3),
+                                     maxshape = (None, 10, 3, 3),
+                                     dtype = np.float64)
+                time_chunk = 2**20//(8*9*self.parameters['histogram_bins'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/histograms/' + k,
+                                     (1,
+                                      self.parameters['histogram_bins'],
+                                      3, 3),
+                                     chunks = (time_chunk,
+                                               self.parameters['histogram_bins'],
+                                               3, 3),
+                                     maxshape = (None,
+                                                 self.parameters['histogram_bins'],
+                                                 3, 3),
+                                     dtype = np.int64)
         return None
     def job_parser_arguments(
             self,

TurTLE/test/test_Gaussian_field.py

@@ -94,6 +94,13 @@ def plot_stuff(simname, total_energy = 1.):
     a.set_ylim(1e-5,1)
     f.tight_layout()
     f.savefig('spectrum_isotropy_test.pdf')
+    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()))
     df.close()
     return None
 

cpp/full_code/Gauss_field_test.cpp

@@ -130,13 +130,43 @@ int Gauss_field_test<rnumber>::do_work(void)
         stat_group = 0;
     }
 
+    /// compute gradient and divergence of field
+    field<rnumber, FFTW, THREExTHREE> *tensor_field = new field<rnumber, FFTW, THREExTHREE>(
+            nx, ny, nz,
+            this->comm,
+            FFTW_ESTIMATE);
+    compute_gradient(
+            this->kk,
+            this->vector_field,
+            tensor_field);
+    compute_divergence(
+            this->kk,
+            this->vector_field,
+            this->scalar_field);
+
     /// compute basic statistics of Gaussian field
     this->vector_field->compute_stats(
             this->kk,
             stat_group,
             "velocity",
             0,
-            8.);
+            this->max_velocity_estimate);
+
+    /// verify divergence free method
+    tensor_field->ift();
+    /// for the stats max value doesn't really matter, I just want the moments
+    tensor_field->compute_rspace_stats(
+            stat_group,
+            "velocity_gradient",
+            0,
+            std::vector<double>({1, 1, 1, 1, 1, 1, 1, 1, 1}));
+    this->scalar_field->ift();
+    this->scalar_field->compute_rspace_stats(
+            stat_group,
+            "velocity_divergence",
+            0,
+            std::vector<double>({1}));
+    delete tensor_field;
 
     /// close stat file
     if (this->myrank == 0)