diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 2dc2283c27d21bd981abc0cf1df39881aa59d995..3da0c8095b505521c66b6c66d24e92486d960f74 100644
--- a/TurTLE/TEST.py
+++ b/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,
diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index 5fdc3506f87f5fc7e964bb779007d303cbfd0ef7..149073e7c1abe931368df1a026ef8f02ac4bc070 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/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
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index 6251d558efec714385c134ebde681d940f7512eb..6822b6d8c02f87a532ffb0d0f63ed9bc767a8564 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/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)