diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 33a820330795e907b6a6fff574d04c3a0c5acdd7..7f82de445391d8628e69bd4177f2ca201cda5779 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -24,6 +24,7 @@
+import sys
import os
import numpy as np
import h5py
@@ -74,7 +75,6 @@ class NavierStokes(_fluid_particle_base):
self.parameters['max_R_estimate'] = 1.0
self.file_datasets_grow = """
//begincpp
- std::string temp_string;
hsize_t dims[4];
hid_t group;
hid_t Cspace, Cdset;
@@ -338,6 +338,7 @@ class NavierStokes(_fluid_particle_base):
self.fluid_includes += '#include <cstring>\n'
self.fluid_variables += ('fluid_solver<{0}> *fs;\n'.format(self.C_dtype) +
'int *kindices;\n' +
+ 'hid_t particle_file;\n' +
'hid_t H5T_field_complex;\n')
self.fluid_definitions += """
typedef struct {{
@@ -373,6 +374,15 @@ class NavierStokes(_fluid_particle_base):
H5Tinsert(H5T_field_complex, "r", HOFFSET(tmp_complex_type, re), {2});
H5Tinsert(H5T_field_complex, "i", HOFFSET(tmp_complex_type, im), {2});
}}
+ if (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 cache for particles I got %d\\n", cache_err);
+ sprintf(fname, "%s_particles.h5", simname);
+ particle_file = H5Fopen(fname, H5F_ACC_RDWR, fapl);
+ }}
//endcpp
""".format(self.C_dtype, self.fftw_plan_rigor, field_H5T)
if not self.frozen_fields:
@@ -480,10 +490,8 @@ class NavierStokes(_fluid_particle_base):
self.parameters['tracers{0}_integration_steps'.format(s0 + s)] = integration_steps[s]
self.file_datasets_grow += """
//begincpp
- temp_string = (std::string("/particles/") +
- std::string(ps{0}->name));
- group = H5Gopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
- grow_particle_datasets(group, temp_string.c_str(), NULL, NULL);
+ group = H5Gopen(particle_file, ps{0}->name, H5P_DEFAULT);
+ grow_particle_datasets(group, "", NULL, NULL);
H5Gclose(group);
//endcpp
""".format(s0 + s)
@@ -515,40 +523,13 @@ class NavierStokes(_fluid_particle_base):
{0}->field = {1};
ps{2}->sample_vec_field({0}, acceleration);
""".format(interpolator[s], acc_name, s0 + s)
- output_vel_acc += """
- if (myrank == 0)
- {{
- //VELOCITY begin
- std::string temp_string = (std::string("/particles/") +
- std::string(ps{0}->name) +
- std::string("/velocity"));
- hid_t Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
- hid_t mspace, wspace;
- int ndims;
- hsize_t count[3], offset[3];
- wspace = H5Dget_space(Cdset);
- ndims = H5Sget_simple_extent_dims(wspace, count, NULL);
- count[0] = 1;
- offset[0] = ps{0}->iteration / ps{0}->traj_skip;
- offset[1] = 0;
- offset[2] = 0;
- mspace = H5Screate_simple(ndims, count, NULL);
- H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, velocity);
- H5Dclose(Cdset);
- //VELOCITY end\n""".format(s0 + s)
+ output_vel_acc += (
+ 'if (myrank == 0)\n' +
+ '{\n' +
+ 'ps{0}->write(particle_file, "velocity", velocity);\n'.format(s0 + s))
if not type(acc_name) == type(None):
- output_vel_acc += """
- //ACCELERATION begin
- temp_string = (std::string("/particles/") +
- std::string(ps{0}->name) +
- std::string("/acceleration"));
- Cdset = H5Dopen(stat_file, temp_string.c_str(), H5P_DEFAULT);
- H5Dwrite(Cdset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, acceleration);
- H5Sclose(mspace);
- H5Sclose(wspace);
- H5Dclose(Cdset);
- //ACCELERATION end\n""".format(s0 + s)
+ output_vel_acc += (
+ 'ps{0}->write(particle_file, "acceleration", acceleration);\n'.format(s0 + s))
output_vel_acc += '}\n'
output_vel_acc += 'delete[] velocity;\n'
if not type(acc_name) == type(None):
@@ -570,7 +551,7 @@ class NavierStokes(_fluid_particle_base):
for s in range(nspecies):
neighbours = self.parameters[interpolator[s] + '_neighbours']
self.particle_start += 'sprintf(fname, "tracers{0}");\n'.format(s0 + s)
- self.particle_end += ('ps{0}->write(stat_file);\n' +
+ self.particle_end += ('ps{0}->write(particle_file);\n' +
'delete ps{0};\n').format(s0 + s)
self.particle_variables += 'rFFTW_particles<VELOCITY_TRACER, {0}, {1}> *ps{2};\n'.format(
self.C_dtype,
@@ -586,7 +567,7 @@ class NavierStokes(_fluid_particle_base):
interpolator[s])
self.particle_start += ('ps{0}->dt = dt;\n' +
'ps{0}->iteration = iteration;\n' +
- 'ps{0}->read(stat_file);\n').format(s0 + s)
+ 'ps{0}->read(particle_file);\n').format(s0 + s)
if not frozen_particles:
if type(kcut) == list:
update_field = ('fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut[s]) +
@@ -594,8 +575,7 @@ class NavierStokes(_fluid_particle_base):
self.particle_loop += update_field
self.particle_loop += '{0}->field = fs->rvelocity;\n'.format(interpolator[s])
self.particle_loop += 'ps{0}->step();\n'.format(s0 + s)
- self.particle_stat_src += 'ps{0}->write(stat_file, false);\n'.format(s0 + s)
- self.particle_start += output_vel_acc
+ self.particle_stat_src += 'ps{0}->write(particle_file, false);\n'.format(s0 + s)
self.particle_stat_src += output_vel_acc
self.particle_stat_src += '}\n'
self.particle_species += nspecies
@@ -604,6 +584,10 @@ class NavierStokes(_fluid_particle_base):
return os.path.join(self.work_dir, self.simname + '.h5')
def get_data_file(self):
return h5py.File(self.get_data_file_name(), 'r')
+ def get_particle_file_name(self):
+ return os.path.join(self.work_dir, self.simname + '_particles.h5')
+ def get_particle_file(self):
+ return h5py.File(self.get_particle_file_name(), 'r')
def get_postprocess_file_name(self):
return os.path.join(self.work_dir, self.simname + '_postprocess.h5')
def get_postprocess_file(self):
@@ -627,11 +611,6 @@ class NavierStokes(_fluid_particle_base):
self.statistics['kshell'] = data_file['kspace/kshell'].value
self.statistics['kM'] = data_file['kspace/kM'].value
self.statistics['dk'] = data_file['kspace/dk'].value
- if self.particle_species > 0:
- self.trajectories = [data_file['particles/' + key + '/state'][
- iter0//self.parameters['niter_part'] :
- iter1//self.parameters['niter_part']+1]
- for key in data_file['particles'].keys()]
computation_needed = True
pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
if 'ii0' in pp_file.keys():
@@ -746,7 +725,7 @@ class NavierStokes(_fluid_particle_base):
3))
def write_par(self, iter0 = 0):
_fluid_particle_base.write_par(self, iter0 = iter0)
- with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r+') as ofile:
+ with h5py.File(self.get_data_file_name(), 'r+') as ofile:
kspace = self.get_kspace()
nshells = kspace['nshell'].shape[0]
for k in ['velocity', 'vorticity']:
@@ -788,44 +767,6 @@ class NavierStokes(_fluid_particle_base):
4),
dtype = np.int64,
compression = 'gzip')
- for s in range(self.particle_species):
- time_chunk = 2**20 // (8*3*
- self.parameters['nparticles']*
- self.parameters['tracers{0}_integration_steps'.format(s)])
- time_chunk = max(time_chunk, 1)
- ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
- (1,
- self.parameters['tracers{0}_integration_steps'.format(s)],
- self.parameters['nparticles'],
- 3),
- maxshape = (None,
- self.parameters['tracers{0}_integration_steps'.format(s)],
- self.parameters['nparticles'],
- 3),
- chunks = (time_chunk,
- self.parameters['tracers{0}_integration_steps'.format(s)],
- self.parameters['nparticles'],
- 3),
- dtype = np.float64)
- time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
- time_chunk = max(time_chunk, 1)
- ofile.create_dataset(
- '/particles/tracers{0}/velocity'.format(s),
- (1,
- self.parameters['nparticles'],
- 3),
- chunks = (time_chunk, self.parameters['nparticles'], 3),
- maxshape = (None, self.parameters['nparticles'], 3),
- dtype = np.float64)
- if self.parameters['tracers{0}_acc_on'.format(s)]:
- ofile.create_dataset(
- '/particles/tracers{0}/acceleration'.format(s),
- (1,
- self.parameters['nparticles'],
- 3),
- chunks = (time_chunk, self.parameters['nparticles'], 3),
- maxshape = (None, self.parameters['nparticles'], 3),
- dtype = np.float64)
if self.QR_stats_on:
time_chunk = 2**20//(8*3*self.parameters['histogram_bins'])
time_chunk = max(time_chunk, 1)
@@ -888,6 +829,80 @@ class NavierStokes(_fluid_particle_base):
self.parameters['QR2D_histogram_bins']),
dtype = np.int64,
compression = 'gzip')
+ if self.particle_species == 0:
+ return None
+ def create_particle_dataset(
+ data_file,
+ dset_name,
+ dset_shape,
+ dset_maxshape,
+ dset_chunks,
+ # maybe something more general can be used here
+ dset_dtype = h5py.h5t.IEEE_F64LE):
+ # create the dataspace.
+ space_id = h5py.h5s.create_simple(
+ dset_shape,
+ dset_maxshape)
+ # create the dataset creation property list.
+ dcpl = h5py.h5p.create(h5py.h5p.DATASET_CREATE)
+ # set the allocation time to "early".
+ dcpl.set_alloc_time(h5py.h5d.ALLOC_TIME_EARLY)
+ dcpl.set_chunk(dset_chunks)
+ # and now create dataset
+ if sys.version_info[0] == 3:
+ dset_name = dset_name.encode()
+ return h5py.h5d.create(
+ data_file.id,
+ dset_name,
+ dset_dtype,
+ space_id,
+ dcpl,
+ h5py.h5p.DEFAULT)
+
+ with h5py.File(self.get_particle_file_name(), 'a') as ofile:
+ for s in range(self.particle_species):
+ ofile.create_group('tracers{0}'.format(s))
+ time_chunk = 2**20 // (8*3*
+ self.parameters['nparticles']*
+ self.parameters['tracers{0}_integration_steps'.format(s)])
+ time_chunk = max(time_chunk, 1)
+ dims = (1,
+ self.parameters['tracers{0}_integration_steps'.format(s)],
+ self.parameters['nparticles'],
+ 3)
+ maxshape = (h5py.h5s.UNLIMITED,
+ self.parameters['tracers{0}_integration_steps'.format(s)],
+ self.parameters['nparticles'],
+ 3)
+ chunks = (time_chunk,
+ self.parameters['tracers{0}_integration_steps'.format(s)],
+ self.parameters['nparticles'],
+ 3)
+ create_particle_dataset(
+ ofile,
+ '/tracers{0}/rhs'.format(s),
+ dims, maxshape, chunks)
+ time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
+ time_chunk = max(time_chunk, 1)
+ create_particle_dataset(
+ ofile,
+ '/tracers{0}/state'.format(s),
+ (1, self.parameters['nparticles'], 3),
+ (h5py.h5s.UNLIMITED, self.parameters['nparticles'], 3),
+ (time_chunk, self.parameters['nparticles'], 3))
+ create_particle_dataset(
+ ofile,
+ '/tracers{0}/velocity'.format(s),
+ (1, self.parameters['nparticles'], 3),
+ (h5py.h5s.UNLIMITED, self.parameters['nparticles'], 3),
+ (time_chunk, self.parameters['nparticles'], 3))
+ if self.parameters['tracers{0}_acc_on'.format(s)]:
+ create_particle_dataset(
+ ofile,
+ '/tracers{0}/acceleration'.format(s),
+ (1, self.parameters['nparticles'], 3),
+ (h5py.h5s.UNLIMITED, self.parameters['nparticles'], 3),
+ (time_chunk, self.parameters['nparticles'], 3))
return None
def add_particle_fields(
self,
diff --git a/bfps/_code.py b/bfps/_code.py
index 3589689cebd20104dd5bc034ecd6bcba63a2df97..80f1fdf4c5f5b3b7e78d4b4d107b096b0050da86 100644
--- a/bfps/_code.py
+++ b/bfps/_code.py
@@ -93,7 +93,13 @@ class _code(_base):
read_parameters(parameter_file);
H5Fclose(parameter_file);
if (myrank == 0)
- stat_file = H5Fopen(fname, H5F_ACC_RDWR, H5P_DEFAULT);
+ {
+ // 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(fname, H5F_ACC_RDWR, fapl);
+ }
//endcpp
"""
for ostream in ['cout', 'cerr']:
@@ -107,6 +113,7 @@ class _code(_base):
H5Dwrite(Cdset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &iteration);
H5Dclose(Cdset);
H5Fclose(stat_file);
+ H5Fclose(particle_file);
}
fftwf_mpi_cleanup();
fftw_mpi_cleanup();
diff --git a/bfps/_fluid_base.py b/bfps/_fluid_base.py
index 0e407098d746ae964f2e6138468f2e3774859169..aece581924900cf9b01548c8305eb852c2529815 100644
--- a/bfps/_fluid_base.py
+++ b/bfps/_fluid_base.py
@@ -117,17 +117,15 @@ class _fluid_particle_base(_code):
'H5Dclose(dset);\n}\n' +
'return 0;\n}\n')
self.definitions += ('herr_t grow_particle_datasets(hid_t g_id, const char *name, const H5L_info_t *info, void *op_data)\n{\n' +
- 'std::string full_name;\n' +
'hsize_t dset;\n')
for key in ['state', 'velocity', 'acceleration']:
- self.definitions += ('full_name = (std::string(name) + std::string("/{0}"));\n'.format(key) +
- 'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
- 'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
+ self.definitions += ('if (H5Lexists(g_id, "{0}", H5P_DEFAULT))\n'.format(key) +
+ '{\n' +
+ 'dset = H5Dopen(g_id, "{0}", H5P_DEFAULT);\n'.format(key) +
'grow_single_dataset(dset, niter_todo/niter_part);\n' +
'H5Dclose(dset);\n}\n')
- self.definitions += ('full_name = (std::string(name) + std::string("/rhs"));\n' +
- 'if (H5Lexists(g_id, full_name.c_str(), H5P_DEFAULT))\n{\n' +
- 'dset = H5Dopen(g_id, full_name.c_str(), H5P_DEFAULT);\n' +
+ self.definitions += ('if (H5Lexists(g_id, "rhs", H5P_DEFAULT))\n{\n' +
+ 'dset = H5Dopen(g_id, "rhs", H5P_DEFAULT);\n' +
'grow_single_dataset(dset, 1);\n' +
'H5Dclose(dset);\n}\n' +
'return 0;\n}\n')
@@ -206,6 +204,7 @@ class _fluid_particle_base(_code):
self.main += '}\n'
self.main += 'do_stats();\n'
self.main += 'do_particle_stats();\n'
+ self.main += output_time_difference
if self.particle_species > 0:
self.main += self.particle_end
self.main += self.fluid_end
@@ -328,19 +327,8 @@ class _fluid_particle_base(_code):
if testing:
#data[0] = np.array([3.26434, 4.24418, 3.12157])
data[0] = np.array([ 0.72086101, 2.59043666, 6.27501953])
- with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r+') as data_file:
- time_chunk = 2**20 // (8*ncomponents*
- self.parameters['nparticles'])
- time_chunk = max(time_chunk, 1)
- dset = data_file.create_dataset(
- '/particles/tracers{0}/state'.format(species),
- (1,
- self.parameters['nparticles'],
- ncomponents),
- chunks = (time_chunk, self.parameters['nparticles'], ncomponents),
- maxshape = (None, self.parameters['nparticles'], ncomponents),
- dtype = np.float64)
- dset[0] = data
+ with h5py.File(self.get_particle_file_name(), 'r+') as data_file:
+ data_file['tracers{0}/state'.format(species)][0] = data
if write_to_file:
data.tofile(
os.path.join(
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
index 4cd79c1259940861c1b0a8f07d147bc5a1356250..a71c0574326340dcdd80ee9ca4a29278b62e4a12 100644
--- a/bfps/cpp/rFFTW_particles.cpp
+++ b/bfps/cpp/rFFTW_particles.cpp
@@ -109,7 +109,8 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(int nsteps)
+void rFFTW_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
+ const int nsteps)
{
ptrdiff_t ii;
this->get_rhs(this->state, this->rhs[0]);
@@ -194,11 +195,12 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::step()
template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read(hid_t data_file_id)
+void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read(
+ const hid_t data_file_id)
{
if (this->myrank == 0)
{
- std::string temp_string = (std::string("/particles/") +
+ std::string temp_string = (std::string("/") +
std::string(this->name) +
std::string("/state"));
hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
@@ -218,7 +220,7 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read(hid_t data
H5Dclose(dset);
if (this->iteration > 0)
{
- temp_string = (std::string("/particles/") +
+ temp_string = (std::string("/") +
std::string(this->name) +
std::string("/rhs"));
dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
@@ -257,42 +259,57 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read(hid_t data
}
template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::write(hid_t data_file_id, bool write_rhs)
+void rFFTW_particles<particle_type, rnumber, interp_neighbours>::write(
+ const hid_t data_file_id,
+ const char *dset_name,
+ const double *data)
+{
+ std::string temp_string = (std::string(this->name) +
+ std::string("/") +
+ std::string(dset_name));
+ hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
+ hid_t mspace, wspace;
+ hsize_t count[3], offset[3];
+ wspace = H5Dget_space(dset);
+ H5Sget_simple_extent_dims(wspace, count, NULL);
+ count[0] = 1;
+ offset[0] = this->iteration / this->traj_skip;
+ offset[1] = 0;
+ offset[2] = 0;
+ mspace = H5Screate_simple(3, count, NULL);
+ H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+ H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, data);
+ H5Sclose(mspace);
+ H5Sclose(wspace);
+ H5Dclose(dset);
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void rFFTW_particles<particle_type, rnumber, interp_neighbours>::write(
+ const hid_t data_file_id,
+ const bool write_rhs)
{
if (this->myrank == 0)
{
- std::string temp_string = (std::string("/particles/") +
- std::string(this->name) +
- std::string("/state"));
- hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
- hid_t mspace, wspace;
- hsize_t count[4], offset[4];
- wspace = H5Dget_space(dset);
- H5Sget_simple_extent_dims(wspace, count, NULL);
- count[0] = 1;
- offset[0] = this->iteration / this->traj_skip;
- offset[1] = 0;
- offset[2] = 0;
- mspace = H5Screate_simple(3, count, NULL);
- H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, this->state);
- H5Sclose(mspace);
- H5Sclose(wspace);
- H5Dclose(dset);
+ this->write(data_file_id, "state", this->state);
if (write_rhs)
{
- temp_string = (std::string("/particles/") +
- std::string(this->name) +
- std::string("/rhs"));
- dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
- wspace = H5Dget_space(dset);
+ std::string temp_string = (
+ std::string("/") +
+ std::string(this->name) +
+ std::string("/rhs"));
+ hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
+ hid_t wspace = H5Dget_space(dset);
+ hsize_t count[4], offset[4];
H5Sget_simple_extent_dims(wspace, count, NULL);
//writing to last available position
offset[0] = count[0] - 1;
+ offset[1] = 0;
+ offset[2] = 0;
+ offset[3] = 0;
count[0] = 1;
count[1] = 1;
- offset[3] = 0;
- mspace = H5Screate_simple(4, count, NULL);
+ hid_t mspace = H5Screate_simple(4, count, NULL);
for (int i=0; i<this->integration_steps; i++)
{
offset[1] = i;
diff --git a/bfps/cpp/rFFTW_particles.hpp b/bfps/cpp/rFFTW_particles.hpp
index b577da055b77f58703cbf175c8215c6178cd7600..488e6612fbb16f954efcca6b85e30bf74f8375fe 100644
--- a/bfps/cpp/rFFTW_particles.hpp
+++ b/bfps/cpp/rFFTW_particles.hpp
@@ -78,21 +78,26 @@ class rFFTW_particles
const int INTEGRATION_STEPS = 2);
~rFFTW_particles();
- void get_rhs(double *__restrict__ x, double *__restrict__ rhs);
+ void get_rhs(
+ double *__restrict__ x,
+ double *__restrict__ rhs);
- inline void sample_vec_field(rFFTW_interpolator<rnumber, interp_neighbours> *field, double *vec_values)
+ inline void sample_vec_field(
+ rFFTW_interpolator<rnumber, interp_neighbours> *field,
+ double *vec_values)
{
field->sample(this->nparticles, this->ncomponents, this->state, vec_values, NULL);
}
/* input/output */
- void read(hid_t data_file_id);
- void write(hid_t data_file_id, bool write_rhs = true);
+ void read(const hid_t data_file_id);
+ void write(const hid_t data_file_id, const char *dset_name, const double *data);
+ void write(const hid_t data_file_id, const bool write_rhs = true);
/* solvers */
void step();
void roll_rhs();
- void AdamsBashforth(int nsteps);
+ void AdamsBashforth(const int nsteps);
};
#endif//RFFTW_PARTICLES
diff --git a/bfps/tools.py b/bfps/tools.py
index 0ecfc768f6a972cfb5520b2998bc6eced110e46e..7502468810fd48f7b5f32fa3a51f66f57906ffd8 100644
--- a/bfps/tools.py
+++ b/bfps/tools.py
@@ -198,8 +198,12 @@ def particle_finite_diff_test(
.. seealso:: :func:`get_fornberg_coeffs`
"""
- d = c.get_data_file()
- group = d['particles/tracers{0}'.format(species)]
+ df = c.get_data_file()
+ if 'particles' in df.keys():
+ group = df['particles/tracers{0}'.format(species)]
+ else:
+ pf = c.get_particle_file()
+ group = pf['tracers{0}'.format(species)]
acc_on = 'acceleration' in group.keys()
pos = group['state'].value
vel = group['velocity'].value
@@ -207,7 +211,7 @@ def particle_finite_diff_test(
acc = group['acceleration'].value
n = m
fc = get_fornberg_coeffs(0, range(-n, n+1))
- dt = d['parameters/dt'].value*d['parameters/niter_part'].value
+ dt = c.parameters['dt']*c.parameters['niter_part']
num_vel1 = sum(fc[1, n-i]*pos[1+n-i:pos.shape[0]-i-n-1] for i in range(-n, n+1)) / dt
if acc_on:
@@ -220,7 +224,7 @@ def particle_finite_diff_test(
pid = np.argmin(SNR(num_acc1, acc[n+1:-n-1]))
else:
pid = np.argmin(SNR(num_vel1, vel[n+1:-n-1]))
- pars = d['parameters']
+ pars = df['parameters']
to_print = (
'steps={0}, interp={1}, neighbours={2}, '.format(
pars['tracers{0}_integration_steps'.format(species)].value,
@@ -246,7 +250,6 @@ def particle_finite_diff_test(
for n in range(1, m):
fc = get_fornberg_coeffs(0, range(-n, n+1))
- dt = d['parameters/dt'].value*d['parameters/niter_part'].value
num_acc1 = sum(fc[1, n-i]*vel[n-i:vel.shape[0]-i-n] for i in range(-n, n+1)) / dt
num_acc2 = sum(fc[2, n-i]*pos[n-i:pos.shape[0]-i-n] for i in range(-n, n+1)) / dt**2
diff --git a/documentation/_static/development.rst b/documentation/_static/development.rst
index 446bc1ab765a281f05410939c74e8c08dabaf3e6..423a448946041d4ec0aca77f29a2311e3f0eb177 100644
--- a/documentation/_static/development.rst
+++ b/documentation/_static/development.rst
@@ -57,3 +57,30 @@ Code testing
------------
Testing for :mod:`bfps` is done with ``tox``.
+
+----------------
+Work in progress
+----------------
+
+HDF5 field I/O
+--------------
+
+As you can tell from the ``todo.txt`` file, in the future the code will
+use HDF5 for input/output of fields.
+For now, field I/O is done with binary files, and the field data type is
+stored in the parameter file.
+
+HDF5 particle I/O
+-----------------
+
+Particle I/O seems to be very slow, no idea why.
+Relevant links:
+
+ 1. http://api.h5py.org/index.html
+ 2. https://www.hdfgroup.org/ftp/HDF5/examples/python/hdf5examples-py/low_level/h5ex_d_alloc.py
+
+Code flexibility
+----------------
+
+Version 2.0 will be designed so that new fluid equations can be coded in
+from python classes, rather than as new C++ files.
diff --git a/tests/base.py b/tests/base.py
index 712b5a24b88eb1bb43a3806d75018df52f237790..7b75080a44e12ffbd11243749fbbcd394671731d 100644
--- a/tests/base.py
+++ b/tests/base.py
@@ -191,7 +191,8 @@ def compare_stats(
for i in range(c0.particle_species):
print('maximum traj difference species {0} is {1}'.format(
i,
- np.max(np.abs(c0.trajectories[i] - c1.trajectories[i]))))
+ np.max(np.abs(c0.get_particle_file()['tracers{0}/state'.format(i)][:] -
+ c1.get_particle_file()['tracers{0}/state'.format(i)][:]))))
if plots_on:
# plot energy and enstrophy
fig = plt.figure(figsize = (12, 12))
diff --git a/tests/test_plain.py b/tests/test_plain.py
index 5d526cfd3412b097068f9ab1b93b2ca86523b34d..ab745d2e3cd4027010e3aab527be6bd477baf4ad 100644
--- a/tests/test_plain.py
+++ b/tests/test_plain.py
@@ -70,7 +70,6 @@ if __name__ == '__main__':
'--nparticles', '1000',
'--niter_todo', '48',
'--precision', 'single',
- '--multiplejob',
'--wd', 'data/single'] +
sys.argv[1:])
plain(opt)