diff --git a/bfps/PP.py b/bfps/PP.py
index ec9994bcabf43f30f6c4b2e4f901481812ba678f..6e02f2aefd5db2e9790f3a16cbc2bfa3c85ab37b 100644
--- a/bfps/PP.py
+++ b/bfps/PP.py
@@ -117,7 +117,7 @@ class PP(_code):
outfile.write(self.main + '\n')
return None
def generate_default_parameters(self):
- # these parameters are relevant for all DNS classes
+ # these parameters are relevant for all PP classes
self.parameters['dealias_type'] = int(1)
self.parameters['dkx'] = float(1.0)
self.parameters['dky'] = float(1.0)
@@ -131,6 +131,15 @@ class PP(_code):
self.pp_parameters = {}
self.pp_parameters['iteration_list'] = np.zeros(1).astype(np.int)
return None
+ def extra_postprocessing_parameters(
+ self,
+ dns_type = 'joint_acc_vel_stats'):
+ pars = {}
+ if dns_type == 'joint_acc_vel_stats':
+ pars['max_acceleration_estimate'] = float(10)
+ pars['max_velocity_estimate'] = float(1)
+ pars['histogram_bins'] = int(129)
+ return pars
def get_data_file_name(self):
return os.path.join(self.work_dir, self.simname + '.h5')
def get_data_file(self):
@@ -426,6 +435,15 @@ class PP(_code):
self.simulation_parser_arguments(parser_get_rfields)
self.job_parser_arguments(parser_get_rfields)
self.parameters_to_parser_arguments(parser_get_rfields)
+ parser_joint_acc_vel_stats = subparsers.add_parser(
+ 'joint_acc_vel_stats',
+ help = 'get joint acceleration and velocity statistics')
+ self.simulation_parser_arguments(parser_joint_acc_vel_stats)
+ self.job_parser_arguments(parser_joint_acc_vel_stats)
+ self.parameters_to_parser_arguments(parser_joint_acc_vel_stats)
+ self.parameters_to_parser_arguments(
+ parser_joint_acc_vel_stats,
+ parameters = self.extra_postprocessing_parameters('joint_acc_vel_stats'))
return None
def prepare_launch(
self,
@@ -463,7 +481,8 @@ class PP(_code):
for k in self.pp_parameters.keys():
self.parameters[k] = self.pp_parameters[k]
self.pars_from_namespace(opt)
- niter_out = self.get_data_file()['parameters/niter_todo'].value
+ niter_out = self.get_data_file()['parameters/niter_out'].value
+ assert(opt.iter0 % niter_out == 0)
self.pp_parameters['iteration_list'] = np.arange(
opt.iter0, opt.iter1+niter_out, niter_out, dtype = np.int)
return opt
@@ -622,13 +641,170 @@ class PP(_code):
dst_file[dst_dset_name][kz,:min_shape[1], :min_shape[2]] = \
src_file[src_dset_name][kz, :min_shape[1], :min_shape[2]]
return None
+ def prepare_post_file(
+ self,
+ opt = None):
+ self.pp_parameters.update(
+ self.extra_postprocessing_parameters(self.dns_type))
+ self.pars_from_namespace(
+ opt,
+ parameters = self.pp_parameters,
+ get_sim_info = False)
+ for kk in ['nx', 'ny', 'nz']:
+ self.parameters[kk] = self.get_data_file()['parameters/' + kk].value
+ n = self.parameters['nx']
+ if self.dns_type in ['filtered_slices',
+ 'filtered_acceleration']:
+ if type(opt.klist_kmax) == type(None):
+ opt.klist_kmax = n / 3.
+ if type(opt.klist_kmin) == type(None):
+ opt.klist_kmin = 6.
+ kvals = bfps_addons.tools.power_space_array(
+ power = opt.klist_power,
+ size = opt.klist_size,
+ vmin = opt.klist_kmin,
+ vmax = opt.klist_kmax)
+ if opt.test_klist:
+ for i in range(opt.klist_size):
+ print('kcut{0} = {1}, ell{0} = {2:.3e}'.format(
+ i, kvals[i], 2*np.pi / kvals[i]))
+ opt.no_submit = True
+ self.pp_parameters['kcut'] = kvals
+ self.rewrite_par(
+ group = self.dns_type + '/parameters',
+ parameters = self.pp_parameters,
+ file_name = os.path.join(self.work_dir, self.simname + '_post.h5'))
+ histogram_bins = opt.histogram_bins
+ if (type(histogram_bins) == type(None) and
+ 'histogram_bins' in self.pp_parameters.keys()):
+ histogram_bins = self.pp_parameters['histogram_bins']
+ with h5py.File(os.path.join(self.work_dir, self.simname + '_post.h5'), 'r+') as ofile:
+ group = ofile[self.dns_type]
+ group.require_group('histograms')
+ group.require_group('moments')
+ group.require_group('spectra')
+ vec_spectra_stats = []
+ vec4_rspace_stats = []
+ scal_rspace_stats = []
+ if self.dns_type == 'joint_acc_vel_stats':
+ vec_spectra_stats.append('velocity')
+ vec4_rspace_stats.append('velocity')
+ vec_spectra_stats.append('acceleration')
+ vec4_rspace_stats.append('acceleration')
+ for quantity in scal_rspace_stats:
+ if quantity not in group['histograms'].keys():
+ time_chunk = 2**20 // (8*histogram_bins)
+ time_chunk = max(time_chunk, 1)
+ group['histograms'].create_dataset(
+ quantity,
+ (1, histogram_bins),
+ chunks = (time_chunk, histogram_bins),
+ maxshape = (None, histogram_bins),
+ dtype = np.int64)
+ else:
+ assert(histogram_bins ==
+ group['histograms/' + quantity].shape[1])
+ if quantity not in group['moments'].keys():
+ time_chunk = 2**20 // (8*10)
+ time_chunk = max(time_chunk, 1)
+ group['moments'].create_dataset(
+ quantity,
+ (1, 10),
+ chunks = (time_chunk, 10),
+ maxshape = (None, 10),
+ dtype = np.float64)
+ if self.dns_type == 'joint_acc_vel_stats':
+ quantity = 'acceleration_and_velocity_components'
+ if quantity not in group['histograms'].keys():
+ time_chunk = 2**20 // (8*9*histogram_bins**2)
+ time_chunk = max(time_chunk, 1)
+ group['histograms'].create_dataset(
+ quantity,
+ (1, histogram_bins, histogram_bins, 3, 3),
+ chunks = (time_chunk, histogram_bins, histogram_bins, 3, 3),
+ maxshape = (None, histogram_bins, histogram_bins, 3, 3),
+ dtype = np.int64)
+ quantity = 'acceleration_and_velocity_magnitudes'
+ if quantity not in group['histograms'].keys():
+ time_chunk = 2**20 // (8*histogram_bins**2)
+ time_chunk = max(time_chunk, 1)
+ group['histograms'].create_dataset(
+ quantity,
+ (1, histogram_bins, histogram_bins),
+ chunks = (time_chunk, histogram_bins, histogram_bins),
+ maxshape = (None, histogram_bins, histogram_bins),
+ dtype = np.int64)
+ ncomps = 4
+ for quantity in vec4_rspace_stats:
+ if quantity not in group['histograms'].keys():
+ time_chunk = 2**20 // (8*histogram_bins*ncomps)
+ time_chunk = max(time_chunk, 1)
+ group['histograms'].create_dataset(
+ quantity,
+ (1, histogram_bins, ncomps),
+ chunks = (time_chunk, histogram_bins, ncomps),
+ maxshape = (None, histogram_bins, ncomps),
+ dtype = np.int64)
+ if quantity not in group['moments'].keys():
+ time_chunk = 2**20 // (8*10*ncomps)
+ time_chunk = max(time_chunk, 1)
+ group['moments'].create_dataset(
+ quantity,
+ (1, 10, ncomps),
+ chunks = (time_chunk, 10, ncomps),
+ maxshape = (None, 10, ncomps),
+ dtype = np.float64)
+ time_chunk = 2**20 // (
+ 4*3*
+ self.parameters['nx']*self.parameters['ny'])
+ time_chunk = max(time_chunk, 1)
+ for quantity in vec_spectra_stats:
+ df = self.get_data_file()
+ if quantity + '_' + quantity not in group['spectra'].keys():
+ spec_chunks = df['statistics/spectra/velocity_velocity'].chunks
+ spec_shape = df['statistics/spectra/velocity_velocity'].shape
+ spec_maxshape = df['statistics/spectra/velocity_velocity'].maxshape
+ group['spectra'].create_dataset(
+ quantity + '_' + quantity,
+ spec_shape,
+ chunks = spec_chunks,
+ maxshape = spec_maxshape,
+ dtype = np.float64)
+ df.close()
+ return None
+ def prepare_field_file(self):
+ df = self.get_data_file()
+ if 'field_dtype' in df.keys():
+ # we don't need to do anything, raw binary files are used
+ return None
+ last_iteration = df['iteration'].value
+ cppf = df['parameters/checkpoints_per_file'].value
+ niter_out = df['parameters/niter_out'].value
+ with h5py.File(os.path.join(self.work_dir, self.simname + '_fields.h5'), 'a') as ff:
+ ff.require_group('vorticity')
+ ff.require_group('vorticity/complex')
+ checkpoint = 0
+ while True:
+ cpf_name = os.path.join(
+ self.work_dir,
+ self.simname + '_checkpoint_{0}.h5'.format(checkpoint))
+ if os.path.exists(cpf_name):
+ cpf = h5py.File(cpf_name, 'r')
+ for iter_name in cpf['vorticity/complex'].keys():
+ if iter_name not in ff['vorticity/complex'].keys():
+ ff['vorticity/complex/' + iter_name] = h5py.ExternalLink(
+ cpf_name,
+ 'vorticity/complex/' + iter_name)
+ checkpoint += 1
+ else:
+ break
+ return None
def launch_jobs(
self,
opt = None,
particle_initial_condition = None):
- self.rewrite_par(
- group = self.dns_type,
- parameters = self.pp_parameters)
+ self.prepare_post_file(opt)
+ self.prepare_field_file()
self.run(
nb_processes = opt.nb_processes,
nb_threads_per_process = opt.nb_threads_per_process,
diff --git a/bfps/cpp/field.cpp b/bfps/cpp/field.cpp
index cc28ced9c829d82132cf061affe5f4e5262f1d2e..0e05983875d3288987c8af153b961e163cac711f 100644
--- a/bfps/cpp/field.cpp
+++ b/bfps/cpp/field.cpp
@@ -1135,6 +1135,183 @@ int invert_curl(
return EXIT_SUCCESS;
}
+template <typename rnumber,
+ field_backend be>
+int joint_rspace_PDF(
+ field<rnumber, be, THREE> *f1,
+ field<rnumber, be, THREE> *f2,
+ const hid_t group,
+ const std::string dset_name,
+ const hsize_t toffset,
+ const std::vector<double> max_f1_estimate,
+ const std::vector<double> max_f2_estimate)
+{
+ TIMEZONE("joint_rspace_PDF");
+ assert(f1->real_space_representation);
+ assert(f2->real_space_representation);
+ assert(max_f1_estimate.size() == 4);
+ assert(max_f2_estimate.size() == 4);
+ int nbins;
+ if (f1->myrank == 0)
+ {
+ hid_t dset, wspace;
+ hsize_t dims[5];
+ int ndims;
+ dset = H5Dopen(
+ group,
+ ("histograms/" + dset_name + "_components").c_str(),
+ H5P_DEFAULT);
+ wspace = H5Dget_space(dset);
+ ndims = H5Sget_simple_extent_dims(wspace, dims, NULL);
+ DEBUG_MSG("number of dimensions is %d\n", ndims);
+ assert(ndims == 5);
+ assert(dims[3] == 3);
+ assert(dims[4] == 3);
+ H5Sclose(wspace);
+ H5Dclose(dset);
+ dset = H5Dopen(
+ group,
+ ("histograms/" + dset_name + "_magnitudes").c_str(),
+ H5P_DEFAULT);
+ wspace = H5Dget_space(dset);
+ ndims = H5Sget_simple_extent_dims(wspace, dims, NULL);
+ assert(ndims == 3);
+ H5Sclose(wspace);
+ H5Dclose(dset);
+ nbins = dims[1];
+ }
+ {
+ TIMEZONE("MPI_Bcast");
+ MPI_Bcast(&nbins, 1, MPI_INT, 0, f1->comm);
+ }
+
+ /// histogram components
+ shared_array<ptrdiff_t> local_histc_threaded(
+ nbins*nbins*9,
+ [&](ptrdiff_t* local_hist){
+ std::fill_n(local_hist, nbins*nbins*9, 0);
+ });
+
+ /// histogram magnitudes
+ shared_array<ptrdiff_t> local_histm_threaded(
+ nbins*nbins,
+ [&](ptrdiff_t* local_hist){
+ std::fill_n(local_hist, nbins*nbins, 0);
+ });
+
+ /// set up bin sizes
+ std::vector<double> bin1size, bin2size;
+ bin1size.resize(4);
+ bin2size.resize(4);
+ for (unsigned int i=0; i<3; i++)
+ {
+ bin1size[i] = 2*max_f1_estimate[i] / nbins;
+ bin2size[i] = 2*max_f2_estimate[i] / nbins;
+ }
+ bin1size[3] = max_f1_estimate[3] / nbins;
+ bin2size[3] = max_f2_estimate[3] / nbins;
+
+
+ {
+ TIMEZONE("field::RLOOP");
+ f1->RLOOP(
+ [&](ptrdiff_t rindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex){
+ ptrdiff_t *local_histc = local_histc_threaded.getMine();
+ ptrdiff_t *local_histm = local_histm_threaded.getMine();
+
+ double mag1, mag2;
+ mag1 = 0;
+ mag2 = 0;
+
+ for (unsigned int i=0; i<3; i++)
+ {
+ double val1 = f1->rval(rindex, i);
+ mag1 += val1*val1;
+ int bin1 = int(floor((val1 + max_f1_estimate[i])/bin1size[i]));
+ for (unsigned int j=0; j<3; j++)
+ {
+ double val2 = f2->rval(rindex, i);
+ mag2 += val2*val2;
+ int bin2 = int(floor((val2 + max_f2_estimate[i])/bin2size[i]));
+ if ((bin1 >= 0 && bin1 < nbins) &&
+ (bin2 >= 0 && bin2 < nbins))
+ local_histc[(bin1*nbins + bin2)*9 + i*3 + j]++;
+ }
+ }
+ int bin1 = int(floor(sqrt(mag1)/bin1size[3]));
+ int bin2 = int(floor(sqrt(mag2)/bin2size[3]));
+ if ((bin1 >= 0 && bin1 < nbins) &&
+ (bin2 >= 0 && bin2 < nbins))
+ local_histm[bin1*nbins + bin2]++;
+ });
+ }
+ local_histm_threaded.mergeParallel();
+ local_histc_threaded.mergeParallel();
+ ptrdiff_t *histm = new ptrdiff_t[nbins*nbins];
+ ptrdiff_t *histc = new ptrdiff_t[nbins*nbins*9];
+ {
+ MPI_Allreduce(
+ (void*)local_histm_threaded.getMasterData(),
+ (void*)histm,
+ nbins*nbins,
+ MPI_INT64_T, MPI_SUM, f1->comm);
+ MPI_Allreduce(
+ (void*)local_histc_threaded.getMasterData(),
+ (void*)histc,
+ nbins*nbins*9,
+ MPI_INT64_T, MPI_SUM, f1->comm);
+ }
+
+ if (f1->myrank == 0)
+ {
+ TIMEZONE("root-work");
+ hid_t dset, wspace, mspace;
+ hsize_t count[5], offset[5];
+ dset = H5Dopen(group, ("histograms/" + dset_name + "_components").c_str(), H5P_DEFAULT);
+ assert(dset > 0);
+ wspace = H5Dget_space(dset);
+ offset[0] = toffset;
+ offset[1] = 0;
+ offset[2] = 0;
+ offset[3] = 0;
+ offset[4] = 0;
+ count[0] = 1;
+ count[1] = nbins;
+ count[2] = nbins;
+ count[3] = 3;
+ count[4] = 3;
+ mspace = H5Screate_simple(5, count, NULL);
+ H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+ H5Dwrite(dset, H5T_NATIVE_INT64, mspace, wspace, H5P_DEFAULT, histc);
+ H5Sclose(wspace);
+ H5Sclose(mspace);
+ H5Dclose(dset);
+ dset = H5Dopen(group, ("histograms/" + dset_name + "_magnitudes").c_str(), H5P_DEFAULT);
+ assert(dset > 0);
+ offset[0] = toffset;
+ offset[1] = 0;
+ offset[2] = 0;
+ count[0] = 1;
+ count[1] = nbins;
+ count[2] = nbins;
+ mspace = H5Screate_simple(3, count, NULL);
+ wspace = H5Dget_space(dset);
+ H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+ H5Dwrite(dset, H5T_NATIVE_INT64, mspace, wspace, H5P_DEFAULT, histm);
+ H5Sclose(wspace);
+ H5Sclose(mspace);
+ H5Dclose(dset);
+ }
+
+ delete[] histm;
+ delete[] histc;
+
+ return EXIT_SUCCESS;
+}
+
template class field<float, FFTW, ONE>;
template class field<float, FFTW, THREE>;
template class field<float, FFTW, THREExTHREE>;
@@ -1209,3 +1386,20 @@ template int invert_curl<double, FFTW, SMOOTH>(
field<double, FFTW, THREE> *,
field<double, FFTW, THREE> *);
+template int joint_rspace_PDF<float, FFTW>(
+ field<float, FFTW, THREE> *,
+ field<float, FFTW, THREE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const std::vector<double>);
+template int joint_rspace_PDF<double, FFTW>(
+ field<double, FFTW, THREE> *,
+ field<double, FFTW, THREE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const std::vector<double>);
+
diff --git a/bfps/cpp/field.hpp b/bfps/cpp/field.hpp
index 82ce7afaf2fb7ba0c3d3f13ed54f4759b45425e9..519bd576fde8ff1f4a0d0f9ac0152c549bd03415 100644
--- a/bfps/cpp/field.hpp
+++ b/bfps/cpp/field.hpp
@@ -306,5 +306,16 @@ int invert_curl(
field<rnumber, be, THREE> *source,
field<rnumber, be, THREE> *destination);
+template <typename rnumber,
+ field_backend be>
+int joint_rspace_PDF(
+ field<rnumber, be, THREE> *f1,
+ field<rnumber, be, THREE> *f2,
+ const hid_t group,
+ const std::string dset_name,
+ const hsize_t toffset,
+ const std::vector<double> max_f1_estimate,
+ const std::vector<double> max_f2_estimate);
+
#endif//FIELD_HPP
diff --git a/bfps/cpp/full_code/joint_acc_vel_stats.cpp b/bfps/cpp/full_code/joint_acc_vel_stats.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d4bd4502e11844b7c9fd19be73f29029d2befa1d
--- /dev/null
+++ b/bfps/cpp/full_code/joint_acc_vel_stats.cpp
@@ -0,0 +1,167 @@
+#include <string>
+#include <cmath>
+#include "joint_acc_vel_stats.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int joint_acc_vel_stats<rnumber>::initialize(void)
+{
+ this->NSVE_field_stats<rnumber>::initialize();
+ this->kk = new kspace<FFTW, SMOOTH>(
+ this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+ this->ve = new vorticity_equation<rnumber, FFTW>(
+ this->simname.c_str(),
+ this->nx,
+ this->ny,
+ this->nz,
+ this->dkx,
+ this->dky,
+ this->dkz,
+ this->vorticity->fftw_plan_rigor);
+ hid_t parameter_file = H5Fopen(
+ (this->simname + std::string(".h5")).c_str(),
+ H5F_ACC_RDONLY,
+ H5P_DEFAULT);
+ hid_t dset = H5Dopen(parameter_file, "/parameters/niter_out", H5P_DEFAULT);
+ H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->niter_out);
+ H5Dclose(dset);
+ if (H5Lexists(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT))
+ {
+ dset = H5Dopen(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT);
+ H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->checkpoints_per_file);
+ H5Dclose(dset);
+ }
+ else
+ this->checkpoints_per_file = 1;
+ H5Fclose(parameter_file);
+ parameter_file = H5Fopen(
+ (this->simname + std::string("_post.h5")).c_str(),
+ H5F_ACC_RDONLY,
+ H5P_DEFAULT);
+ this->iteration_list = hdf5_tools::read_vector<int>(
+ parameter_file,
+ "/joint_acc_vel_stats/parameters/iteration_list");
+ dset = H5Dopen(parameter_file, "joint_acc_vel_stats/parameters/max_acceleration_estimate", H5P_DEFAULT);
+ H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->max_acceleration_estimate);
+ H5Dclose(dset);
+ dset = H5Dopen(parameter_file, "joint_acc_vel_stats/parameters/max_velocity_estimate", H5P_DEFAULT);
+ H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->max_velocity_estimate);
+ H5Dclose(dset);
+ H5Fclose(parameter_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);
+ this->stat_file = H5Fopen(
+ (this->simname + "_post.h5").c_str(),
+ H5F_ACC_RDWR,
+ fapl);
+ }
+ else
+ {
+ this->stat_file = 0;
+ }
+ int data_file_problem;
+ if (this->myrank == 0)
+ data_file_problem = hdf5_tools::require_size_file_datasets(
+ this->stat_file,
+ "joint_acc_vel_stats",
+ (this->iteration_list.back() / this->niter_out) + 1);
+ MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
+ if (data_file_problem > 0)
+ {
+ std::cerr <<
+ data_file_problem <<
+ " problems setting sizes of file datasets.\ntrying to exit now." <<
+ std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int joint_acc_vel_stats<rnumber>::work_on_current_iteration(void)
+{
+ DEBUG_MSG("entered joint_acc_vel_stats::work_on_current_iteration\n");
+ /// read current vorticity, place it in this->ve->cvorticity
+ this->read_current_cvorticity();
+ *this->ve->cvorticity = this->vorticity->get_cdata();
+ /// after the previous instruction, we are free to use this->vorticity
+ /// for any other purpose
+
+ /// initialize `stat_group`.
+ hid_t stat_group;
+ if (this->myrank == 0)
+ stat_group = H5Gopen(
+ this->stat_file,
+ "joint_acc_vel_stats",
+ H5P_DEFAULT);
+ else
+ stat_group = 0;
+
+ field<rnumber, FFTW, THREE> *vel;
+ field<rnumber, FFTW, THREE> *acc;
+
+ /// compute velocity
+ vel = new field<rnumber, FFTW, THREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ DEFAULT_FFTW_FLAG);
+ invert_curl(kk, this->ve->cvorticity, vel);
+ vel->ift();
+
+ /// compute Lagrangian acceleration
+ /// use this->vorticity as temporary field
+ acc = this->vorticity;
+ this->ve->compute_Lagrangian_acceleration(acc);
+ acc->ift();
+
+ /// compute single field real space statistics
+ std::vector<double> max_acc_estimate;
+ std::vector<double> max_vel_estimate;
+
+ max_acc_estimate.resize(4, max_acceleration_estimate);
+ max_vel_estimate.resize(4, max_velocity_estimate);
+
+ acc->compute_rspace_stats(
+ stat_group,
+ "acceleration",
+ this->iteration / this->niter_out,
+ max_acc_estimate);
+ vel->compute_rspace_stats(
+ stat_group,
+ "velocity",
+ this->iteration / this->niter_out,
+ max_vel_estimate);
+
+ /// compute joint PDF
+ joint_rspace_PDF(
+ acc, vel,
+ stat_group,
+ "acceleration_and_velocity",
+ this->iteration / this->niter_out,
+ max_acc_estimate,
+ max_vel_estimate);
+
+ delete vel;
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int joint_acc_vel_stats<rnumber>::finalize(void)
+{
+ delete this->ve;
+ delete this->kk;
+ if (this->myrank == 0)
+ H5Fclose(this->stat_file);
+ this->NSVE_field_stats<rnumber>::finalize();
+ return EXIT_SUCCESS;
+}
+
+template class joint_acc_vel_stats<float>;
+template class joint_acc_vel_stats<double>;
+
diff --git a/bfps/cpp/full_code/joint_acc_vel_stats.hpp b/bfps/cpp/full_code/joint_acc_vel_stats.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..721698b7620dbe0ae282da28bd38835b475d50ca
--- /dev/null
+++ b/bfps/cpp/full_code/joint_acc_vel_stats.hpp
@@ -0,0 +1,67 @@
+/**********************************************************************
+* *
+* Copyright 2017 Max Planck Institute *
+* for Dynamics and Self-Organization *
+* *
+* This file is part of bfps. *
+* *
+* bfps is free software: you can redistribute it and/or modify *
+* it under the terms of the GNU General Public License as published *
+* by the Free Software Foundation, either version 3 of the License, *
+* or (at your option) any later version. *
+* *
+* bfps is distributed in the hope that it will be useful, *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+* GNU General Public License for more details. *
+* *
+* You should have received a copy of the GNU General Public License *
+* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@ds.mpg.de *
+* *
+**********************************************************************/
+
+
+
+#ifndef JOINT_ACC_VEL_STATS_HPP
+#define JOINT_ACC_VEL_STATS_HPP
+
+#include <cstdlib>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <vector>
+#include "base.hpp"
+#include "field.hpp"
+#include "field_binary_IO.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE_field_stats.hpp"
+
+template <typename rnumber>
+class joint_acc_vel_stats: public NSVE_field_stats<rnumber>
+{
+ public:
+ hid_t stat_file;
+ kspace<FFTW, SMOOTH> *kk;
+ vorticity_equation<rnumber, FFTW> *ve;
+
+ int checkpoints_per_file;
+ int niter_out;
+ double max_acceleration_estimate;
+ double max_velocity_estimate;
+
+ joint_acc_vel_stats(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ NSVE_field_stats<rnumber>(
+ COMMUNICATOR,
+ simulation_name){}
+ virtual ~joint_acc_vel_stats(){}
+
+ int initialize(void);
+ int work_on_current_iteration(void);
+ int finalize(void);
+};
+
+#endif//JOINT_ACC_VEL_STATS_HPP
+
diff --git a/bfps/cpp/vorticity_equation.cpp b/bfps/cpp/vorticity_equation.cpp
index 82f3c23dd7f445a11d9a680493271c9e9a5081e8..300e9a24f8aa06ffcd34c751c6761c316759cea0 100644
--- a/bfps/cpp/vorticity_equation.cpp
+++ b/bfps/cpp/vorticity_equation.cpp
@@ -140,8 +140,6 @@ vorticity_equation<rnumber, be>::vorticity_equation(
this->cvelocity = new field<rnumber, be, THREE>(
nx, ny, nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
- this->rvelocity = new field<rnumber, be, THREE>(
- nx, ny, nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
this->u = this->cvelocity;
/* initialize kspace */
@@ -168,7 +166,6 @@ vorticity_equation<rnumber, be>::~vorticity_equation()
delete this->v[1];
delete this->v[2];
delete this->cvelocity;
- delete this->rvelocity;
}
template <class rnumber,
diff --git a/bfps/cpp/vorticity_equation.hpp b/bfps/cpp/vorticity_equation.hpp
index 880f05d2d93ada14390f9ed0a003c5cc299633e2..e8bd1d843f730d39439bc99703956dc623ca4e42 100644
--- a/bfps/cpp/vorticity_equation.hpp
+++ b/bfps/cpp/vorticity_equation.hpp
@@ -58,7 +58,7 @@ class vorticity_equation
/* fields */
field<rnumber, be, THREE> *cvorticity, *cvelocity;
- field<rnumber, be, THREE> *rvorticity, *rvelocity;
+ field<rnumber, be, THREE> *rvorticity;
kspace<be, SMOOTH> *kk;
diff --git a/setup.py b/setup.py
index 78db9dcf3d66aa954f0a56173d28f6c10af50cc4..9bba17014aabf36c685395843b806f650604face 100644
--- a/setup.py
+++ b/setup.py
@@ -88,7 +88,8 @@ print('This is bfps version ' + VERSION)
### lists of files and MANIFEST.in
-src_file_list = ['full_code/test',
+src_file_list = ['full_code/joint_acc_vel_stats',
+ 'full_code/test',
'full_code/filter_test',
'hdf5_tools',
'full_code/get_rfields',