Skip to content

GitLab

Commit 31b0ec63 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/separate-particle-file' into develop

parents 5456d7b9 b880d863
Hide whitespace changes
Inline Side-by-side
bfps/NavierStokes.py
View file @ 31b0ec63
......@@ -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,
......
bfps/_code.py
View file @ 31b0ec63
......@@ -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();
......
bfps/_fluid_base.py
View file @ 31b0ec63
......@@ -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(
......
bfps/cpp/rFFTW_particles.cpp
View file @ 31b0ec63
......@@ -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;
......
bfps/cpp/rFFTW_particles.hpp
View file @ 31b0ec63
......@@ -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
......
bfps/tools.py
View file @ 31b0ec63
......@@ -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
......
documentation/_static/development.rst
View file @ 31b0ec63
......@@ -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.
tests/base.py
View file @ 31b0ec63
......@@ -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))
......
tests/test_plain.py
View file @ 31b0ec63
......@@ -70,7 +70,6 @@ if __name__ == '__main__':
'--nparticles', '1000',
'--niter_todo', '48',
'--precision', 'single',
'--multiplejob',
'--wd', 'data/single'] +
sys.argv[1:])
plain(opt)
......