diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 535897672e9c243adc501eb767c3e904bd2a2abc..32ef1ef79fef5b7365310b5b26f02bf6b2baaa43 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -78,7 +78,6 @@ class NavierStokes(_fluid_particle_base):
hsize_t dims[4];
hid_t group;
hid_t Cspace, Cdset;
- int ndims;
// store kspace information
Cdset = H5Dopen(stat_file, "/kspace/kshell", H5P_DEFAULT);
Cspace = H5Dget_space(Cdset);
@@ -337,8 +336,7 @@ class NavierStokes(_fluid_particle_base):
def fill_up_fluid_code(self):
self.fluid_includes += '#include <cstring>\n'
self.fluid_variables += ('fluid_solver<{0}> *fs;\n'.format(self.C_dtype) +
- 'hid_t particle_file;\n' +
- 'hid_t H5T_field_complex;\n')
+ 'hid_t particle_file;\n')
self.fluid_definitions += """
typedef struct {{
{0} re;
@@ -367,12 +365,6 @@ class NavierStokes(_fluid_particle_base):
strncpy(fs->forcing_type, forcing_type, 128);
fs->iteration = iteration;
fs->read('v', 'c');
- if (fs->cd->myrank == 0)
- {{
- H5T_field_complex = H5Tcreate(H5T_COMPOUND, sizeof(tmp_complex_type));
- H5Tinsert(H5T_field_complex, "r", HOFFSET(tmp_complex_type, re), {2});
- H5Tinsert(H5T_field_complex, "i", HOFFSET(tmp_complex_type, im), {2});
- }}
//endcpp
""".format(self.C_dtype, self.fftw_plan_rigor, field_H5T)
if self.parameters['nparticles'] > 0:
@@ -395,13 +387,12 @@ class NavierStokes(_fluid_particle_base):
self.fluid_output + '\n}\n')
self.fluid_end = ('if (fs->iteration % niter_out != 0)\n{\n' +
self.fluid_output + '\n}\n' +
- 'if (fs->cd->myrank == 0)\n' +
- '{\n' +
- 'H5Tclose(H5T_field_complex);\n' +
- '}\n' +
'delete fs;\n')
if self.parameters['nparticles'] > 0:
- self.fluid_end += 'H5Fclose(particle_file);\n'
+ self.fluid_end += ('if (myrank == 0)\n' +
+ '{\n' +
+ 'H5Fclose(particle_file);\n' +
+ '}\n')
return None
def add_3D_rFFTW_field(
self,
@@ -420,10 +411,11 @@ class NavierStokes(_fluid_particle_base):
neighbours = 1,
smoothness = 1,
name = 'field_interpolator',
- field_name = 'fs->rvelocity'):
- self.fluid_includes += '#include "rFFTW_interpolator.hpp"\n'
- self.fluid_variables += 'rFFTW_interpolator <{0}, {1}> *{2};\n'.format(
- self.C_dtype, neighbours, name)
+ field_name = 'fs->rvelocity',
+ class_name = 'rFFTW_interpolator'):
+ self.fluid_includes += '#include "{0}.hpp"\n'.format(class_name)
+ self.fluid_variables += '{0} <{1}, {2}> *{3};\n'.format(
+ class_name, self.C_dtype, neighbours, name)
self.parameters[name + '_type'] = interp_type
self.parameters[name + '_neighbours'] = neighbours
if interp_type == 'spline':
@@ -431,8 +423,9 @@ class NavierStokes(_fluid_particle_base):
beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
elif interp_type == 'Lagrange':
beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
- self.fluid_start += '{0} = new rFFTW_interpolator<{1}, {2}>(fs, {3}, {4});\n'.format(
+ self.fluid_start += '{0} = new {1}<{2}, {3}>(fs, {4}, {5});\n'.format(
name,
+ class_name,
self.C_dtype,
neighbours,
beta_name,
@@ -445,7 +438,8 @@ class NavierStokes(_fluid_particle_base):
kcut = None,
interpolator = 'field_interpolator',
frozen_particles = False,
- acc_name = None):
+ acc_name = None,
+ class_name = 'particles'):
"""Adds code for tracking a series of particle species, each
consisting of `nparticles` particles.
@@ -493,7 +487,7 @@ class NavierStokes(_fluid_particle_base):
self.parameters['tracers{0}_integration_steps'.format(s0 + s)] = integration_steps[s]
self.file_datasets_grow += """
//begincpp
- group = H5Gopen(particle_file, ps{0}->name, H5P_DEFAULT);
+ group = H5Gopen(particle_file, ps{0}->get_name(), H5P_DEFAULT);
grow_particle_datasets(group, "", NULL, NULL);
H5Gclose(group);
//endcpp
@@ -504,39 +498,26 @@ class NavierStokes(_fluid_particle_base):
# array for putting sampled velocity in
# must compute velocity, just in case it was messed up by some
# other particle species before the stats
- output_vel_acc += ('double *velocity = new double[3*nparticles];\n' +
- 'fs->compute_velocity(fs->cvorticity);\n')
+ output_vel_acc += 'fs->compute_velocity(fs->cvorticity);\n'
if not type(kcut) == list:
output_vel_acc += 'fs->ift_velocity();\n'
if not type(acc_name) == type(None):
# array for putting sampled acceleration in
# must compute acceleration
- output_vel_acc += 'double *acceleration = new double[3*nparticles];\n'
output_vel_acc += 'fs->compute_Lagrangian_acceleration({0});\n'.format(acc_name)
for s in range(nspecies):
if type(kcut) == list:
output_vel_acc += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut[s])
output_vel_acc += 'fs->ift_velocity();\n'
output_vel_acc += """
- {0}->field = fs->rvelocity;
- ps{1}->sample_vec_field({0}, velocity);
+ {0}->read_rFFTW(fs->rvelocity);
+ ps{1}->sample({0}, "velocity");
""".format(interpolator[s], s0 + s)
if not type(acc_name) == type(None):
output_vel_acc += """
- {0}->field = {1};
- ps{2}->sample_vec_field({0}, acceleration);
+ {0}->read_rFFTW({1});
+ ps{2}->sample({0}, "acceleration");
""".format(interpolator[s], acc_name, 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 += (
- '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):
- output_vel_acc += 'delete[] acceleration;\n'
output_vel_acc += '}\n'
#### initialize, stepping and finalize code
@@ -547,38 +528,39 @@ class NavierStokes(_fluid_particle_base):
self.particle_loop += update_fields
else:
self.particle_loop += 'fs->compute_velocity(fs->cvorticity);\n'
- self.particle_includes += '#include "rFFTW_particles.hpp"\n'
+ self.particle_includes += '#include "{0}.hpp"\n'.format(class_name)
self.particle_stat_src += (
'if (ps0->iteration % niter_part == 0)\n' +
'{\n')
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(particle_file);\n' +
+ self.particle_end += ('ps{0}->write();\n' +
'delete ps{0};\n').format(s0 + s)
- self.particle_variables += 'rFFTW_particles<VELOCITY_TRACER, {0}, {1}> *ps{2};\n'.format(
+ self.particle_variables += '{0}<VELOCITY_TRACER, {1}, {2}> *ps{3};\n'.format(
+ class_name,
self.C_dtype,
neighbours,
s0 + s)
- self.particle_start += ('ps{0} = new rFFTW_particles<VELOCITY_TRACER, {1}, {2}>(\n' +
- 'fname, {3},\n' +
- 'nparticles,\n' +
+ self.particle_start += ('ps{0} = new {1}<VELOCITY_TRACER, {2}, {3}>(\n' +
+ 'fname, particle_file, {4},\n' +
'niter_part, tracers{0}_integration_steps);\n').format(
s0 + s,
+ class_name,
self.C_dtype,
neighbours,
interpolator[s])
self.particle_start += ('ps{0}->dt = dt;\n' +
'ps{0}->iteration = iteration;\n' +
- 'ps{0}->read(particle_file);\n').format(s0 + s)
+ 'ps{0}->read();\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]) +
'fs->ift_velocity();\n')
self.particle_loop += update_field
- self.particle_loop += '{0}->field = fs->rvelocity;\n'.format(interpolator[s])
+ self.particle_loop += '{0}->read_rFFTW(fs->rvelocity);\n'.format(interpolator[s])
self.particle_loop += 'ps{0}->step();\n'.format(s0 + s)
- self.particle_stat_src += 'ps{0}->write(particle_file, false);\n'.format(s0 + s)
+ self.particle_stat_src += 'ps{0}->write(false);\n'.format(s0 + s)
self.particle_stat_src += output_vel_acc
self.particle_stat_src += '}\n'
self.particle_species += nspecies
@@ -898,9 +880,7 @@ class NavierStokes(_fluid_particle_base):
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 = 2**20 // (8*3*self.parameters['nparticles'])
time_chunk = max(time_chunk, 1)
dims = (1,
self.parameters['tracers{0}_integration_steps'.format(s)],
@@ -911,15 +891,13 @@ class NavierStokes(_fluid_particle_base):
self.parameters['nparticles'],
3)
chunks = (time_chunk,
- self.parameters['tracers{0}_integration_steps'.format(s)],
+ 1,
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),
diff --git a/bfps/_fluid_base.py b/bfps/_fluid_base.py
index aece581924900cf9b01548c8305eb852c2529815..5bb539f795228ac3c4562fb69e198816b3fd090e 100644
--- a/bfps/_fluid_base.py
+++ b/bfps/_fluid_base.py
@@ -102,13 +102,15 @@ class _fluid_particle_base(_code):
self.definitions += self.particle_definitions
self.definitions += ('int grow_single_dataset(hid_t dset, int tincrement)\n{\n' +
'int ndims;\n' +
- 'hsize_t dims[5];\n' +
'hsize_t space;\n' +
'space = H5Dget_space(dset);\n' +
- 'ndims = H5Sget_simple_extent_dims(space, dims, NULL);\n' +
+ 'ndims = H5Sget_simple_extent_ndims(space);\n' +
+ 'hsize_t *dims = new hsize_t[ndims];\n' +
+ 'H5Sget_simple_extent_dims(space, dims, NULL);\n' +
'dims[0] += tincrement;\n' +
'H5Dset_extent(dset, dims);\n' +
'H5Sclose(space);\n' +
+ 'delete[] dims;\n' +
'return EXIT_SUCCESS;\n}\n')
self.definitions += ('herr_t grow_statistics_dataset(hid_t o_id, const char *name, const H5O_info_t *info, void *op_data)\n{\n' +
'if (info->type == H5O_TYPE_DATASET)\n{\n' +
diff --git a/bfps/cpp/distributed_particles.cpp b/bfps/cpp/distributed_particles.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b0a445af03ab51d79dd1aded6ba230fc4b80adbf
--- /dev/null
+++ b/bfps/cpp/distributed_particles.cpp
@@ -0,0 +1,461 @@
+/**********************************************************************
+* *
+* Copyright 2015 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 *
+* *
+**********************************************************************/
+
+
+
+#define NDEBUG
+
+#include <cmath>
+#include <cassert>
+#include <cstring>
+#include <string>
+#include <sstream>
+
+#include "base.hpp"
+#include "distributed_particles.hpp"
+#include "fftw_tools.hpp"
+
+
+extern int myrank, nprocs;
+
+template <int particle_type, class rnumber, int interp_neighbours>
+distributed_particles<particle_type, rnumber, interp_neighbours>::distributed_particles(
+ const char *NAME,
+ const hid_t data_file_id,
+ interpolator<rnumber, interp_neighbours> *FIELD,
+ const int TRAJ_SKIP,
+ const int INTEGRATION_STEPS) : particles_io_base<particle_type>(
+ NAME,
+ TRAJ_SKIP,
+ data_file_id,
+ FIELD->descriptor->comm)
+{
+ assert((INTEGRATION_STEPS <= 6) &&
+ (INTEGRATION_STEPS >= 1));
+ this->vel = FIELD;
+ this->rhs.resize(INTEGRATION_STEPS);
+ this->integration_steps = INTEGRATION_STEPS;
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+distributed_particles<particle_type, rnumber, interp_neighbours>::~distributed_particles()
+{
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
+ interpolator<rnumber, interp_neighbours> *field,
+ std::unordered_map<int, single_particle_state<POINT3D>> &y)
+{
+ double *yy = new double[3];
+ y.clear();
+ for (auto &pp: this->state)
+ {
+ (*field)(pp.second.data, yy);
+ y[pp.first] = yy;
+ }
+ delete[] yy;
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::get_rhs(
+ const std::unordered_map<int, single_particle_state<particle_type>> &x,
+ std::unordered_map<int, single_particle_state<particle_type>> &y)
+{
+ double *yy = new double[this->ncomponents];
+ y.clear();
+ for (auto &pp: x)
+ {
+ (*this->vel)(pp.second.data, yy);
+ y[pp.first] = yy;
+ }
+ delete[] yy;
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::sample(
+ interpolator<rnumber, interp_neighbours> *field,
+ const char *dset_name)
+{
+ std::unordered_map<int, single_particle_state<POINT3D>> y;
+ this->sample(field, y);
+ this->write(dset_name, y);
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
+{
+ for (int i=this->integration_steps-2; i>=0; i--)
+ rhs[i+1] = rhs[i];
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::redistribute(
+ std::unordered_map<int, single_particle_state<particle_type>> &x,
+ std::vector<std::unordered_map<int, single_particle_state<particle_type>>> &vals)
+{
+ //DEBUG_MSG("entered redistribute\n");
+ /* neighbouring rank offsets */
+ int ro[2];
+ ro[0] = -1;
+ ro[1] = 1;
+ /* neighbouring ranks */
+ int nr[2];
+ nr[0] = MOD(this->myrank+ro[0], this->nprocs);
+ nr[1] = MOD(this->myrank+ro[1], this->nprocs);
+ /* particles to send, particles to receive */
+ std::vector<int> ps[2], pr[2];
+ /* number of particles to send, number of particles to receive */
+ int nps[2], npr[2];
+ int rsrc, rdst;
+ /* get list of id-s to send */
+ for (auto &pp: x)
+ for (int i=0; i<2; i++)
+ if (this->vel->get_rank(pp.second.data[2]) == nr[i])
+ ps[i].push_back(pp.first);
+ /* prepare data for send recv */
+ for (int i=0; i<2; i++)
+ nps[i] = ps[i].size();
+ for (rsrc = 0; rsrc<this->nprocs; rsrc++)
+ for (int i=0; i<2; i++)
+ {
+ rdst = MOD(rsrc+ro[i], this->nprocs);
+ if (this->myrank == rsrc)
+ MPI_Send(
+ nps+i,
+ 1,
+ MPI_INTEGER,
+ rdst,
+ 2*(rsrc*this->nprocs + rdst)+i,
+ this->comm);
+ if (this->myrank == rdst)
+ MPI_Recv(
+ npr+1-i,
+ 1,
+ MPI_INTEGER,
+ rsrc,
+ 2*(rsrc*this->nprocs + rdst)+i,
+ this->comm,
+ MPI_STATUS_IGNORE);
+ }
+ //DEBUG_MSG("I have to send %d %d particles\n", nps[0], nps[1]);
+ //DEBUG_MSG("I have to recv %d %d particles\n", npr[0], npr[1]);
+ for (int i=0; i<2; i++)
+ pr[i].resize(npr[i]);
+
+ int buffer_size = (nps[0] > nps[1]) ? nps[0] : nps[1];
+ buffer_size = (buffer_size > npr[0])? buffer_size : npr[0];
+ buffer_size = (buffer_size > npr[1])? buffer_size : npr[1];
+ //DEBUG_MSG("buffer size is %d\n", buffer_size);
+ double *buffer = new double[buffer_size*this->ncomponents*(1+vals.size())];
+ for (rsrc = 0; rsrc<this->nprocs; rsrc++)
+ for (int i=0; i<2; i++)
+ {
+ rdst = MOD(rsrc+ro[i], this->nprocs);
+ if (this->myrank == rsrc && nps[i] > 0)
+ {
+ MPI_Send(
+ &ps[i].front(),
+ nps[i],
+ MPI_INTEGER,
+ rdst,
+ 2*(rsrc*this->nprocs + rdst),
+ this->comm);
+ int pcounter = 0;
+ for (int p: ps[i])
+ {
+ std::copy(x[p].data,
+ x[p].data + this->ncomponents,
+ buffer + pcounter*(1+vals.size())*this->ncomponents);
+ x.erase(p);
+ for (int tindex=0; tindex<vals.size(); tindex++)
+ {
+ std::copy(vals[tindex][p].data,
+ vals[tindex][p].data + this->ncomponents,
+ buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents);
+ vals[tindex].erase(p);
+ }
+ pcounter++;
+ }
+ MPI_Send(
+ buffer,
+ nps[i]*(1+vals.size())*this->ncomponents,
+ MPI_DOUBLE,
+ rdst,
+ 2*(rsrc*this->nprocs + rdst)+1,
+ this->comm);
+ }
+ if (this->myrank == rdst && npr[1-i] > 0)
+ {
+ MPI_Recv(
+ &pr[1-i].front(),
+ npr[1-i],
+ MPI_INTEGER,
+ rsrc,
+ 2*(rsrc*this->nprocs + rdst),
+ this->comm,
+ MPI_STATUS_IGNORE);
+ MPI_Recv(
+ buffer,
+ npr[1-i]*(1+vals.size())*this->ncomponents,
+ MPI_DOUBLE,
+ rsrc,
+ 2*(rsrc*this->nprocs + rdst)+1,
+ this->comm,
+ MPI_STATUS_IGNORE);
+ int pcounter = 0;
+ for (int p: pr[1-i])
+ {
+ x[p] = buffer + (pcounter*(1+vals.size()))*this->ncomponents;
+ for (int tindex=0; tindex<vals.size(); tindex++)
+ {
+ vals[tindex][p] = buffer + (pcounter*(1+vals.size()) + tindex+1)*this->ncomponents;
+ }
+ pcounter++;
+ }
+ }
+ }
+ delete[] buffer;
+
+
+#ifndef NDEBUG
+ /* check that all particles at x are local */
+ for (auto &pp: x)
+ if (this->vel->get_rank(pp.second.data[2]) != this->myrank)
+ {
+ DEBUG_MSG("found particle %d with rank %d\n",
+ pp.first,
+ this->vel->get_rank(pp.second.data[2]));
+ assert(false);
+ }
+#endif
+ //DEBUG_MSG("exiting redistribute\n");
+}
+
+
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
+ const int nsteps)
+{
+ this->get_rhs(this->state, this->rhs[0]);
+ for (auto &pp: this->state)
+ for (int i=0; i<this->ncomponents; i++)
+ switch(nsteps)
+ {
+ case 1:
+ pp.second[i] += this->dt*this->rhs[0][pp.first][i];
+ break;
+ case 2:
+ pp.second[i] += this->dt*(3*this->rhs[0][pp.first][i]
+ - this->rhs[1][pp.first][i])/2;
+ break;
+ case 3:
+ pp.second[i] += this->dt*(23*this->rhs[0][pp.first][i]
+ - 16*this->rhs[1][pp.first][i]
+ + 5*this->rhs[2][pp.first][i])/12;
+ break;
+ case 4:
+ pp.second[i] += this->dt*(55*this->rhs[0][pp.first][i]
+ - 59*this->rhs[1][pp.first][i]
+ + 37*this->rhs[2][pp.first][i]
+ - 9*this->rhs[3][pp.first][i])/24;
+ break;
+ case 5:
+ pp.second[i] += this->dt*(1901*this->rhs[0][pp.first][i]
+ - 2774*this->rhs[1][pp.first][i]
+ + 2616*this->rhs[2][pp.first][i]
+ - 1274*this->rhs[3][pp.first][i]
+ + 251*this->rhs[4][pp.first][i])/720;
+ break;
+ case 6:
+ pp.second[i] += this->dt*(4277*this->rhs[0][pp.first][i]
+ - 7923*this->rhs[1][pp.first][i]
+ + 9982*this->rhs[2][pp.first][i]
+ - 7298*this->rhs[3][pp.first][i]
+ + 2877*this->rhs[4][pp.first][i]
+ - 475*this->rhs[5][pp.first][i])/1440;
+ break;
+ }
+ this->redistribute(this->state, this->rhs);
+ this->roll_rhs();
+}
+
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::step()
+{
+ this->AdamsBashforth((this->iteration < this->integration_steps) ?
+ this->iteration+1 :
+ this->integration_steps);
+ this->iteration++;
+}
+
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::read()
+{
+ double *temp = new double[this->chunk_size*this->ncomponents];
+ for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
+ {
+ //read state
+ if (this->myrank == 0)
+ this->read_state_chunk(cindex, temp);
+ MPI_Bcast(
+ temp,
+ this->chunk_size*this->ncomponents,
+ MPI_DOUBLE,
+ 0,
+ this->comm);
+ for (int p=0; p<this->chunk_size; p++)
+ {
+ if (this->vel->get_rank(temp[this->ncomponents*p+2]) == this->myrank)
+ this->state[p+cindex*this->chunk_size] = temp + this->ncomponents*p;
+ }
+ //read rhs
+ if (this->iteration > 0)
+ for (int i=0; i<this->integration_steps; i++)
+ {
+ if (this->myrank == 0)
+ this->read_rhs_chunk(cindex, i, temp);
+ MPI_Bcast(
+ temp,
+ this->chunk_size*this->ncomponents,
+ MPI_DOUBLE,
+ 0,
+ this->comm);
+ for (int p=0; p<this->chunk_size; p++)
+ {
+ auto pp = this->state.find(p+cindex*this->chunk_size);
+ if (pp != this->state.end())
+ this->rhs[i][p+cindex*this->chunk_size] = temp + this->ncomponents*p;
+ }
+ }
+ }
+ DEBUG_MSG("%s->state.size = %ld\n", this->name.c_str(), this->state.size());
+ delete[] temp;
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
+ const char *dset_name,
+ std::unordered_map<int, single_particle_state<POINT3D>> &y)
+{
+ double *data = new double[this->nparticles*3];
+ double *yy = new double[this->nparticles*3];
+ for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
+ {
+ std::fill_n(yy, this->chunk_size*3, 0);
+ for (int p=0; p<this->chunk_size; p++)
+ {
+ auto pp = y.find(p+cindex*this->chunk_size);
+ if (pp != y.end())
+ std::copy(pp->second.data,
+ pp->second.data + 3,
+ yy + pp->first*3);
+ }
+ MPI_Allreduce(
+ yy,
+ data,
+ 3*this->nparticles,
+ MPI_DOUBLE,
+ MPI_SUM,
+ this->comm);
+ if (this->myrank == 0)
+ this->write_point3D_chunk(dset_name, cindex, data);
+ }
+ delete[] yy;
+ delete[] data;
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void distributed_particles<particle_type, rnumber, interp_neighbours>::write(
+ const bool write_rhs)
+{
+ double *temp0 = new double[this->chunk_size*this->ncomponents];
+ double *temp1 = new double[this->chunk_size*this->ncomponents];
+ for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
+ {
+ //write state
+ std::fill_n(temp0, this->ncomponents*this->chunk_size, 0);
+ for (int p=0; p<this->chunk_size; p++)
+ {
+ auto pp = this->state.find(p + cindex*this->chunk_size);
+ if (pp != this->state.end())
+ std::copy(pp->second.data,
+ pp->second.data + this->ncomponents,
+ temp0 + pp->first*this->ncomponents);
+ }
+ MPI_Allreduce(
+ temp0,
+ temp1,
+ this->ncomponents*this->chunk_size,
+ MPI_DOUBLE,
+ MPI_SUM,
+ this->comm);
+ if (this->myrank == 0)
+ this->write_state_chunk(cindex, temp1);
+ //write rhs
+ if (write_rhs)
+ for (int i=0; i<this->integration_steps; i++)
+ {
+ std::fill_n(temp0, this->ncomponents*this->chunk_size, 0);
+ for (int p=0; p<this->chunk_size; p++)
+ {
+ auto pp = this->rhs[i].find(p + cindex*this->chunk_size);
+ if (pp != this->rhs[i].end())
+ std::copy(pp->second.data,
+ pp->second.data + this->ncomponents,
+ temp0 + pp->first*this->ncomponents);
+ }
+ MPI_Allreduce(
+ temp0,
+ temp1,
+ this->ncomponents*this->chunk_size,
+ MPI_DOUBLE,
+ MPI_SUM,
+ this->comm);
+ if (this->myrank == 0)
+ this->write_rhs_chunk(cindex, i, temp1);
+ }
+ }
+ delete[] temp0;
+ delete[] temp1;
+}
+
+
+/*****************************************************************************/
+template class distributed_particles<VELOCITY_TRACER, float, 1>;
+template class distributed_particles<VELOCITY_TRACER, float, 2>;
+template class distributed_particles<VELOCITY_TRACER, float, 3>;
+template class distributed_particles<VELOCITY_TRACER, float, 4>;
+template class distributed_particles<VELOCITY_TRACER, float, 5>;
+template class distributed_particles<VELOCITY_TRACER, float, 6>;
+template class distributed_particles<VELOCITY_TRACER, double, 1>;
+template class distributed_particles<VELOCITY_TRACER, double, 2>;
+template class distributed_particles<VELOCITY_TRACER, double, 3>;
+template class distributed_particles<VELOCITY_TRACER, double, 4>;
+template class distributed_particles<VELOCITY_TRACER, double, 5>;
+template class distributed_particles<VELOCITY_TRACER, double, 6>;
+/*****************************************************************************/
diff --git a/bfps/cpp/rFFTW_particles.hpp b/bfps/cpp/distributed_particles.hpp
similarity index 58%
rename from bfps/cpp/rFFTW_particles.hpp
rename to bfps/cpp/distributed_particles.hpp
index 488e6612fbb16f954efcca6b85e30bf74f8375fe..e972f69943814823e8b5527422b9727ab7331a8d 100644
--- a/bfps/cpp/rFFTW_particles.hpp
+++ b/bfps/cpp/distributed_particles.hpp
@@ -27,40 +27,31 @@
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
+#include <unordered_map>
+#include <vector>
#include <hdf5.h>
#include "base.hpp"
#include "particles_base.hpp"
#include "fluid_solver_base.hpp"
-#include "rFFTW_interpolator.hpp"
+#include "interpolator.hpp"
-#ifndef RFFTW_PARTICLES
+#ifndef DISTRIBUTED_PARTICLES
-#define RFFTW_PARTICLES
+#define DISTRIBUTED_PARTICLES
template <int particle_type, class rnumber, int interp_neighbours>
-class rFFTW_particles
+class distributed_particles: public particles_io_base<particle_type>
{
- public:
- int myrank, nprocs;
- MPI_Comm comm;
- rnumber *data;
+ private:
+ std::unordered_map<int, single_particle_state<particle_type> > state;
+ std::vector<std::unordered_map<int, single_particle_state<particle_type>>> rhs;
- /* state will generally hold all the information about the particles.
- * in the beginning, we will only need to solve 3D ODEs, but I figured
- * a general ncomponents is better, since we may change our minds.
- * */
- double *state;
- double *rhs[6];
- int nparticles;
- int ncomponents;
- int array_size;
+ public:
int integration_steps;
- int traj_skip;
- rFFTW_interpolator<rnumber, interp_neighbours> *vel;
+ // this class only works with buffered interpolator
+ interpolator<rnumber, interp_neighbours> *vel;
/* simulation parameters */
- char name[256];
- int iteration;
double dt;
/* methods */
@@ -70,29 +61,38 @@ class rFFTW_particles
* this->state
* this->rhs
* */
- rFFTW_particles(
+ distributed_particles(
const char *NAME,
- rFFTW_interpolator<rnumber, interp_neighbours> *FIELD,
- const int NPARTICLES,
+ const hid_t data_file_id,
+ interpolator<rnumber, interp_neighbours> *FIELD,
const int TRAJ_SKIP,
const int INTEGRATION_STEPS = 2);
- ~rFFTW_particles();
-
+ ~distributed_particles();
+
+ void sample(
+ interpolator<rnumber, interp_neighbours> *field,
+ const char *dset_name);
+ void sample(
+ interpolator<rnumber, interp_neighbours> *field,
+ std::unordered_map<int, single_particle_state<POINT3D>> &y);
void get_rhs(
- double *__restrict__ x,
- double *__restrict__ rhs);
+ const std::unordered_map<int, single_particle_state<particle_type>> &x,
+ std::unordered_map<int, single_particle_state<particle_type>> &y);
+
+ void redistribute(
+ std::unordered_map<int, single_particle_state<particle_type>> &x,
+ std::vector<std::unordered_map<int, single_particle_state<particle_type>>> &vals);
- 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(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);
+ void read();
+ void write(
+ const char *dset_name,
+ std::unordered_map<int, single_particle_state<POINT3D>> &y);
+ void write(
+ const char *dset_name,
+ std::unordered_map<int, single_particle_state<particle_type>> &y);
+ void write(const bool write_rhs = true);
/* solvers */
void step();
@@ -100,5 +100,5 @@ class rFFTW_particles
void AdamsBashforth(const int nsteps);
};
-#endif//RFFTW_PARTICLES
+#endif//DISTRIBUTED_PARTICLES
diff --git a/bfps/cpp/fluid_solver_base.cpp b/bfps/cpp/fluid_solver_base.cpp
index db1cd59771c243d5fff007dfe8f39d66c22a1952..6d826d675c0addd605cf1027ae9871b5c81059ef 100644
--- a/bfps/cpp/fluid_solver_base.cpp
+++ b/bfps/cpp/fluid_solver_base.cpp
@@ -346,7 +346,7 @@ fluid_solver_base<R>::fluid_solver_base( \
(void*)(this->kshell), \
this->nshells, \
MPI_DOUBLE, MPI_SUM, this->cd->comm); \
- for (int n=0; n<this->nshells; n++) \
+ for (unsigned int n=0; n<this->nshells; n++) \
{ \
this->kshell[n] /= this->nshell[n]; \
} \
diff --git a/bfps/cpp/fluid_solver_base.hpp b/bfps/cpp/fluid_solver_base.hpp
index 26ae5ded6399e1d3ac36c760f2ea53b2e10f011c..fc7f67f2da68e7805c9079edc173152d3f234360 100644
--- a/bfps/cpp/fluid_solver_base.hpp
+++ b/bfps/cpp/fluid_solver_base.hpp
@@ -67,7 +67,7 @@ class fluid_solver_base
std::unordered_map<int, double> Fourier_filter;
double *kshell;
int64_t *nshell;
- int nshells;
+ unsigned int nshells;
/* methods */
diff --git a/bfps/cpp/interpolator.cpp b/bfps/cpp/interpolator.cpp
index cf04af668c0cee55b1e4181ebc2654d1d4612965..ef53742a4fdeb2545f02954f10c47d2bcb3f6538 100644
--- a/bfps/cpp/interpolator.cpp
+++ b/bfps/cpp/interpolator.cpp
@@ -24,129 +24,134 @@
+#define NDEBUG
+
#include "interpolator.hpp"
template <class rnumber, int interp_neighbours>
interpolator<rnumber, interp_neighbours>::interpolator(
fluid_solver_base<rnumber> *fs,
- base_polynomial_values BETA_POLYS)
+ base_polynomial_values BETA_POLYS,
+ ...) : interpolator_base<rnumber, interp_neighbours>(fs, BETA_POLYS)
{
int tdims[4];
- this->unbuffered_descriptor = fs->rd;
- this->buffer_size = (interp_neighbours+1)*this->unbuffered_descriptor->slice_size;
this->compute_beta = BETA_POLYS;
- tdims[0] = (interp_neighbours+1)*2*this->unbuffered_descriptor->nprocs + this->unbuffered_descriptor->sizes[0];
- tdims[1] = this->unbuffered_descriptor->sizes[1];
- tdims[2] = this->unbuffered_descriptor->sizes[2];
- tdims[3] = this->unbuffered_descriptor->sizes[3];
- this->descriptor = new field_descriptor<rnumber>(
+ tdims[0] = (interp_neighbours+1)*2*this->descriptor->nprocs + this->descriptor->sizes[0];
+ tdims[1] = this->descriptor->sizes[1];
+ tdims[2] = this->descriptor->sizes[2]+2;
+ tdims[3] = this->descriptor->sizes[3];
+ this->buffered_descriptor = new field_descriptor<rnumber>(
4, tdims,
- this->unbuffered_descriptor->mpi_dtype,
- this->unbuffered_descriptor->comm);
- if (sizeof(rnumber) == 4)
- {
- this->f0 = (rnumber*)((void*)fftwf_alloc_real(this->descriptor->local_size));
- this->f1 = (rnumber*)((void*)fftwf_alloc_real(this->descriptor->local_size));
- }
- else if (sizeof(rnumber) == 8)
- {
- this->f0 = (rnumber*)((void*)fftw_alloc_real(this->descriptor->local_size));
- this->f1 = (rnumber*)((void*)fftw_alloc_real(this->descriptor->local_size));
- }
- this->temp = this->f1 + this->buffer_size;
+ this->descriptor->mpi_dtype,
+ this->descriptor->comm);
+ this->buffer_size = (interp_neighbours+1)*this->buffered_descriptor->slice_size;
+ this->field = new rnumber[this->buffered_descriptor->local_size];
}
template <class rnumber, int interp_neighbours>
interpolator<rnumber, interp_neighbours>::~interpolator()
{
- if (sizeof(rnumber) == 4)
- {
- fftwf_free((float*)((void*)this->f0));
- fftwf_free((float*)((void*)this->f1));
- }
- else if (sizeof(rnumber) == 8)
- {
- fftw_free((double*)((void*)this->f0));
- fftw_free((double*)((void*)this->f1));
- }
- delete this->descriptor;
+ delete[] this->field;
+ delete this->buffered_descriptor;
}
template <class rnumber, int interp_neighbours>
-int interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
+int interpolator<rnumber, interp_neighbours>::read_rFFTW(const void *void_src)
{
- /* first, roll fields */
- rnumber *tmp = this->f0;
- this->f0 = this->f1;
- this->f1 = tmp;
- this->temp = this->f0 + this->buffer_size;
- /* now do regular things */
rnumber *src = (rnumber*)void_src;
- rnumber *dst = this->f1;
+ rnumber *dst = this->field;
/* do big copy of middle stuff */
- clip_zero_padding<rnumber>(this->unbuffered_descriptor, src, 3);
std::copy(src,
- src + this->unbuffered_descriptor->local_size,
+ src + this->buffered_descriptor->slice_size*this->descriptor->subsizes[0],
dst + this->buffer_size);
MPI_Datatype MPI_RNUM = (sizeof(rnumber) == 4) ? MPI_FLOAT : MPI_DOUBLE;
int rsrc;
/* get upper slices */
- for (int rdst = 0; rdst < this->unbuffered_descriptor->nprocs; rdst++)
+ for (int rdst = 0; rdst < this->descriptor->nprocs; rdst++)
{
- rsrc = this->unbuffered_descriptor->rank[(this->unbuffered_descriptor->all_start0[rdst] +
- this->unbuffered_descriptor->all_size0[rdst]) %
- this->unbuffered_descriptor->sizes[0]];
- if (this->unbuffered_descriptor->myrank == rsrc)
+ rsrc = this->descriptor->rank[(this->descriptor->all_start0[rdst] +
+ this->descriptor->all_size0[rdst]) %
+ this->descriptor->sizes[0]];
+ if (this->descriptor->myrank == rsrc)
MPI_Send(
src,
this->buffer_size,
MPI_RNUM,
rdst,
- 2*(rsrc*this->unbuffered_descriptor->nprocs + rdst),
- this->descriptor->comm);
- if (this->unbuffered_descriptor->myrank == rdst)
+ 2*(rsrc*this->descriptor->nprocs + rdst),
+ this->buffered_descriptor->comm);
+ if (this->descriptor->myrank == rdst)
MPI_Recv(
- dst + this->buffer_size + this->unbuffered_descriptor->local_size,
+ dst + this->buffer_size + this->buffered_descriptor->slice_size*this->descriptor->subsizes[0],
this->buffer_size,
MPI_RNUM,
rsrc,
- 2*(rsrc*this->unbuffered_descriptor->nprocs + rdst),
- this->descriptor->comm,
+ 2*(rsrc*this->descriptor->nprocs + rdst),
+ this->buffered_descriptor->comm,
MPI_STATUS_IGNORE);
}
/* get lower slices */
- for (int rdst = 0; rdst < this->unbuffered_descriptor->nprocs; rdst++)
+ for (int rdst = 0; rdst < this->descriptor->nprocs; rdst++)
{
- rsrc = this->unbuffered_descriptor->rank[MOD(this->unbuffered_descriptor->all_start0[rdst] - 1,
- this->unbuffered_descriptor->sizes[0])];
- if (this->unbuffered_descriptor->myrank == rsrc)
+ rsrc = this->descriptor->rank[MOD(this->descriptor->all_start0[rdst] - 1,
+ this->descriptor->sizes[0])];
+ if (this->descriptor->myrank == rsrc)
MPI_Send(
- src + this->unbuffered_descriptor->local_size - this->buffer_size,
+ src + this->buffered_descriptor->slice_size*this->descriptor->subsizes[0] - this->buffer_size,
this->buffer_size,
MPI_RNUM,
rdst,
- 2*(rsrc*this->unbuffered_descriptor->nprocs + rdst)+1,
- this->unbuffered_descriptor->comm);
- if (this->unbuffered_descriptor->myrank == rdst)
+ 2*(rsrc*this->descriptor->nprocs + rdst)+1,
+ this->descriptor->comm);
+ if (this->descriptor->myrank == rdst)
MPI_Recv(
dst,
this->buffer_size,
MPI_RNUM,
rsrc,
- 2*(rsrc*this->unbuffered_descriptor->nprocs + rdst)+1,
- this->unbuffered_descriptor->comm,
+ 2*(rsrc*this->descriptor->nprocs + rdst)+1,
+ this->descriptor->comm,
MPI_STATUS_IGNORE);
}
return EXIT_SUCCESS;
}
+template <class rnumber, int interp_neighbours>
+void interpolator<rnumber, interp_neighbours>::sample(
+ const int nparticles,
+ const int pdimension,
+ const double *__restrict__ x,
+ double *__restrict__ y,
+ const int *deriv)
+{
+ /* get grid coordinates */
+ int *xg = new int[3*nparticles];
+ double *xx = new double[3*nparticles];
+ double *yy = new double[3*nparticles];
+ std::fill_n(yy, 3*nparticles, 0.0);
+ this->get_grid_coordinates(nparticles, pdimension, x, xg, xx);
+ /* perform interpolation */
+ for (int p=0; p<nparticles; p++)
+ if (this->descriptor->rank[MOD(xg[p*3+2], this->descriptor->sizes[0])] == this->descriptor->myrank)
+ this->operator()(xg + p*3, xx + p*3, yy + p*3, deriv);
+ MPI_Allreduce(
+ yy,
+ y,
+ 3*nparticles,
+ MPI_DOUBLE,
+ MPI_SUM,
+ this->descriptor->comm);
+ delete[] yy;
+ delete[] xg;
+ delete[] xx;
+}
+
template <class rnumber, int interp_neighbours>
void interpolator<rnumber, interp_neighbours>::operator()(
- double t,
- int *xg,
- double *xx,
+ const int *xg,
+ const double *xx,
double *dest,
- int *deriv)
+ const int *deriv)
{
double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
if (deriv == NULL)
@@ -163,32 +168,26 @@ void interpolator<rnumber, interp_neighbours>::operator()(
}
std::fill_n(dest, 3, 0);
ptrdiff_t bigiz, bigiy, bigix;
- double tval[3];
for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
{
- bigiz = ptrdiff_t(xg[2]+iz);
- std::fill_n(tval, 3, 0);
+ bigiz = ptrdiff_t(xg[2]+iz)-this->descriptor->starts[0];
for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
{
bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
for (int ix = -interp_neighbours; ix <= interp_neighbours+1; ix++)
{
bigix = ptrdiff_t(MOD(xg[0]+ix, this->descriptor->sizes[2]));
+ ptrdiff_t tindex = ((bigiz *this->buffered_descriptor->sizes[1] +
+ bigiy)*this->buffered_descriptor->sizes[2] +
+ bigix)*3 + this->buffer_size;
for (int c=0; c<3; c++)
{
- ptrdiff_t tindex = ((bigiz *this->descriptor->sizes[1] +
- bigiy)*this->descriptor->sizes[2] +
- bigix)*3+c + this->buffer_size;
- dest[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
- by[iy+interp_neighbours]*
- bx[ix+interp_neighbours]);
- tval[c] += (this->f0[tindex]*(1-t) + t*this->f1[tindex])*(bz[iz+interp_neighbours]*
- by[iy+interp_neighbours]*
- bx[ix+interp_neighbours]);
+ dest[c] += this->field[tindex+c]*(bz[iz+interp_neighbours]*
+ by[iy+interp_neighbours]*
+ bx[ix+interp_neighbours]);
}
}
}
- DEBUG_MSG("%ld %d %d %g %g %g\n", bigiz, xg[1], xg[0], tval[0], tval[1], tval[2]);
}
}
diff --git a/bfps/cpp/interpolator.hpp b/bfps/cpp/interpolator.hpp
index e083bd56c6c27c1cc048d0ecfdfb7ee416082651..d20ee6ee1b5023d6a5e177d15e8ddc92a339d7ca 100644
--- a/bfps/cpp/interpolator.hpp
+++ b/bfps/cpp/interpolator.hpp
@@ -24,6 +24,7 @@
+#include <cmath>
#include "field_descriptor.hpp"
#include "fftw_tools.hpp"
#include "fluid_solver_base.hpp"
@@ -34,23 +35,54 @@
#define INTERPOLATOR
template <class rnumber, int interp_neighbours>
-class interpolator
+class interpolator:public interpolator_base<rnumber, interp_neighbours>
{
+ private:
+ /* pointer to buffered field */
+ rnumber *field;
+
public:
ptrdiff_t buffer_size;
- base_polynomial_values compute_beta;
- field_descriptor<rnumber> *descriptor;
- field_descriptor<rnumber> *unbuffered_descriptor;
- rnumber *f0, *f1, *temp;
+
+ /* descriptor for buffered field */
+ field_descriptor<rnumber> *buffered_descriptor;
interpolator(
fluid_solver_base<rnumber> *FSOLVER,
- base_polynomial_values BETA_POLYS);
+ base_polynomial_values BETA_POLYS,
+ ...);
~interpolator();
- void operator()(double t, int *__restrict__ xg, double *__restrict__ xx, double *__restrict__ dest, int *deriv = NULL);
- /* destroys input */
- int read_rFFTW(void *src);
+ int read_rFFTW(const void *src);
+
+ inline int get_rank(double z)
+ {
+ return this->descriptor->rank[MOD(int(floor(z/this->dz)), this->descriptor->sizes[0])];
+ }
+
+ /* interpolate field at an array of locations */
+ void sample(
+ const int nparticles,
+ const int pdimension,
+ const double *__restrict__ x,
+ double *__restrict__ y,
+ const int *deriv = NULL);
+ /* interpolate 1 point */
+ inline void operator()(
+ const double *__restrict__ x,
+ double *__restrict__ dest,
+ const int *deriv = NULL)
+ {
+ int xg[3];
+ double xx[3];
+ this->get_grid_coordinates(x, xg, xx);
+ (*this)(xg, xx, dest, deriv);
+ }
+ void operator()(
+ const int *__restrict__ xg,
+ const double *__restrict__ xx,
+ double *__restrict__ dest,
+ const int *deriv = NULL);
};
#endif//INTERPOLATOR
diff --git a/bfps/cpp/interpolator_base.cpp b/bfps/cpp/interpolator_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..58bf57cf13382f0704da4537dae9d21bb4a841da
--- /dev/null
+++ b/bfps/cpp/interpolator_base.cpp
@@ -0,0 +1,91 @@
+/**********************************************************************
+* *
+* Copyright 2015 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 *
+* *
+**********************************************************************/
+
+
+
+#define NDEBUG
+
+#include <cmath>
+#include "interpolator_base.hpp"
+
+template <class rnumber, int interp_neighbours>
+interpolator_base<rnumber, interp_neighbours>::interpolator_base(
+ fluid_solver_base<rnumber> *fs,
+ base_polynomial_values BETA_POLYS)
+{
+ this->descriptor = fs->rd;
+ this->compute_beta = BETA_POLYS;
+
+ // compute dx, dy, dz;
+ this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
+ this->dy = 4*acos(0) / (fs->dky*this->descriptor->sizes[1]);
+ this->dz = 4*acos(0) / (fs->dkz*this->descriptor->sizes[0]);
+}
+
+template <class rnumber, int interp_neighbours>
+void interpolator_base<rnumber, interp_neighbours>::get_grid_coordinates(
+ const int nparticles,
+ const int pdimension,
+ const double *x,
+ int *xg,
+ double *xx)
+{
+ for (int p=0; p<nparticles; p++)
+ this->get_grid_coordinates(
+ x + p*pdimension,
+ xg + p*3,
+ xx + p*3);
+}
+
+template <class rnumber, int interp_neighbours>
+void interpolator_base<rnumber, interp_neighbours>::get_grid_coordinates(
+ const double *x,
+ int *xg,
+ double *xx)
+{
+ static double grid_size[] = {this->dx, this->dy, this->dz};
+ double tval;
+ for (int c=0; c<3; c++)
+ {
+ tval = floor(x[c]/grid_size[c]);
+ xg[c] = MOD(int(tval), this->descriptor->sizes[2-c]);
+ xx[c] = (x[c] - tval*grid_size[c]) / grid_size[c];
+ }
+}
+
+
+
+template class interpolator_base<float, 1>;
+template class interpolator_base<float, 2>;
+template class interpolator_base<float, 3>;
+template class interpolator_base<float, 4>;
+template class interpolator_base<float, 5>;
+template class interpolator_base<float, 6>;
+template class interpolator_base<double, 1>;
+template class interpolator_base<double, 2>;
+template class interpolator_base<double, 3>;
+template class interpolator_base<double, 4>;
+template class interpolator_base<double, 5>;
+template class interpolator_base<double, 6>;
+
diff --git a/bfps/cpp/interpolator_base.hpp b/bfps/cpp/interpolator_base.hpp
index 69c2d4f27bb3b029fc5fd4c4e85d6160b23fe32c..1f29a5c05e2afa2dda45f45e5f87a71a63d2b25d 100644
--- a/bfps/cpp/interpolator_base.hpp
+++ b/bfps/cpp/interpolator_base.hpp
@@ -24,6 +24,7 @@
+#include "fluid_solver_base.hpp"
#include "spline_n1.hpp"
#include "spline_n2.hpp"
#include "spline_n3.hpp"
@@ -41,5 +42,52 @@ typedef void (*base_polynomial_values)(
const double fraction,
double *__restrict__ destination);
+template <class rnumber, int interp_neighbours>
+class interpolator_base
+{
+ public:
+ /* pointer to polynomial function */
+ base_polynomial_values compute_beta;
+
+ /* descriptor of field to interpolate */
+ field_descriptor<rnumber> *descriptor;
+
+ /* physical parameters of field */
+ double dx, dy, dz;
+
+ interpolator_base(
+ fluid_solver_base<rnumber> *FSOLVER,
+ base_polynomial_values BETA_POLYS);
+ virtual ~interpolator_base(){}
+
+ /* may not destroy input */
+ virtual int read_rFFTW(const void *src) = 0;
+
+ /* map real locations to grid coordinates */
+ void get_grid_coordinates(
+ const int nparticles,
+ const int pdimension,
+ const double *__restrict__ x,
+ int *__restrict__ xg,
+ double *__restrict__ xx);
+ void get_grid_coordinates(
+ const double *__restrict__ x,
+ int *__restrict__ xg,
+ double *__restrict__ xx);
+ /* interpolate field at an array of locations */
+ virtual void sample(
+ const int nparticles,
+ const int pdimension,
+ const double *__restrict__ x,
+ double *__restrict__ y,
+ const int *deriv = NULL) = 0;
+ /* interpolate 1 point */
+ virtual void operator()(
+ const int *__restrict__ xg,
+ const double *__restrict__ xx,
+ double *__restrict__ dest,
+ const int *deriv = NULL) = 0;
+};
+
#endif//INTERPOLATOR_BASE
diff --git a/bfps/cpp/particles.cpp b/bfps/cpp/particles.cpp
index 6d7e14290f1bd69decff1ffa5c38d0e015001d90..2182e769191d587aea2f5b174993ccae407bc109 100644
--- a/bfps/cpp/particles.cpp
+++ b/bfps/cpp/particles.cpp
@@ -24,8 +24,7 @@
-#define NDEBUG
-
+//#define NDEBUG
#include <cmath>
#include <cassert>
@@ -40,311 +39,55 @@
extern int myrank, nprocs;
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-particles<particle_type, rnumber, multistep, interp_neighbours>::particles(
+template <int particle_type, class rnumber, int interp_neighbours>
+particles<particle_type, rnumber, interp_neighbours>::particles(
const char *NAME,
- fluid_solver_base<rnumber> *FSOLVER,
- interpolator<rnumber, interp_neighbours> *FIELD,
- const int NPARTICLES,
+ const hid_t data_file_id,
+ interpolator_base<rnumber, interp_neighbours> *FIELD,
const int TRAJ_SKIP,
- const int INTEGRATION_STEPS)
+ const int INTEGRATION_STEPS) : particles_io_base<particle_type>(
+ NAME,
+ TRAJ_SKIP,
+ data_file_id,
+ FIELD->descriptor->comm)
{
- switch(particle_type)
- {
- case VELOCITY_TRACER:
- this->ncomponents = 3;
- break;
- default:
- this->ncomponents = 3;
- break;
- }
assert((INTEGRATION_STEPS <= 6) &&
(INTEGRATION_STEPS >= 1));
- strncpy(this->name, NAME, 256);
- this->fs = FSOLVER;
- this->nparticles = NPARTICLES;
this->vel = FIELD;
this->integration_steps = INTEGRATION_STEPS;
- this->traj_skip = TRAJ_SKIP;
- this->myrank = this->vel->descriptor->myrank;
- this->nprocs = this->vel->descriptor->nprocs;
- this->comm = this->vel->descriptor->comm;
this->array_size = this->nparticles * this->ncomponents;
- this->state = fftw_alloc_real(this->array_size);
+ this->state = new double[this->array_size];
std::fill_n(this->state, this->array_size, 0.0);
for (int i=0; i < this->integration_steps; i++)
{
- this->rhs[i] = fftw_alloc_real(this->array_size);
+ this->rhs[i] = new double[this->array_size];
std::fill_n(this->rhs[i], this->array_size, 0.0);
}
- this->watching = new bool[this->fs->rd->nprocs*nparticles];
- std::fill_n(this->watching, this->fs->rd->nprocs*this->nparticles, false);
- this->computing = new int[nparticles];
-
- // compute dx, dy, dz;
- this->dx = 4*acos(0) / (this->fs->dkx*this->fs->rd->sizes[2]);
- this->dy = 4*acos(0) / (this->fs->dky*this->fs->rd->sizes[1]);
- this->dz = 4*acos(0) / (this->fs->dkz*this->fs->rd->sizes[0]);
-
- // compute lower and upper bounds
- this->lbound = new double[nprocs];
- this->ubound = new double[nprocs];
- double *tbound = new double[nprocs];
- std::fill_n(tbound, nprocs, 0.0);
- tbound[this->myrank] = this->fs->rd->starts[0]*this->dz;
- MPI_Allreduce(
- tbound,
- this->lbound,
- nprocs,
- MPI_DOUBLE,
- MPI_SUM,
- this->comm);
- std::fill_n(tbound, nprocs, 0.0);
- tbound[this->myrank] = (this->fs->rd->starts[0] + this->fs->rd->subsizes[0])*this->dz;
- MPI_Allreduce(
- tbound,
- this->ubound,
- nprocs,
- MPI_DOUBLE,
- MPI_SUM,
- this->comm);
- delete[] tbound;
- //for (int r = 0; r<nprocs; r++)
- // DEBUG_MSG(
- // "lbound[%d] = %lg, ubound[%d] = %lg\n",
- // r, this->lbound[r],
- // r, this->ubound[r]
- // );
}
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-particles<particle_type, rnumber, multistep, interp_neighbours>::~particles()
+template <int particle_type, class rnumber, int interp_neighbours>
+particles<particle_type, rnumber, interp_neighbours>::~particles()
{
- delete[] this->computing;
- delete[] this->watching;
- fftw_free(this->state);
+ delete[] this->state;
for (int i=0; i < this->integration_steps; i++)
{
- fftw_free(this->rhs[i]);
- }
- delete[] this->lbound;
- delete[] this->ubound;
-}
-
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::sample_vec_field(
- interpolator<rnumber, interp_neighbours> *vec,
- double t,
- double *x,
- double *y,
- const bool synch,
- int *deriv)
-{
- /* get grid coordinates */
- int *xg = new int[3*this->nparticles];
- double *xx = new double[3*this->nparticles];
- this->get_grid_coordinates(x, xg, xx);
- /* perform interpolation */
- std::fill_n(y, 3*this->nparticles, 0.0);
- if (multistep)
- {
- for (int p=0; p<this->nparticles; p++)
- {
- if (this->myrank == this->computing[p])
- (*vec)(t, xg + p*3, xx + p*3, y + p*3, deriv);
- }
- }
- else
- {
- for (int p=0; p<this->nparticles; p++)
- {
- if (this->watching[this->myrank*this->nparticles+p])
- {
- int crank = this->get_rank(x[p*3 + 2]);
- if (this->myrank == crank)
- (*vec)(t, xg + p*3, xx + p*3, y + p*3, deriv);
- if (crank != this->computing[p])
- this->synchronize_single_particle_state(p, y, crank);
- }
- }
- }
- if (synch)
- {
- double *yy = new double[3*this->nparticles];
- std::fill_n(yy, 3*this->nparticles, 0.0);
- for (int p=0; p<this->nparticles; p++)
- {
- if (this->myrank == this->computing[p])
- std::copy(y+3*p, y+3*(p+1), yy+3*p);
- }
- MPI_Allreduce(
- yy,
- y,
- 3*this->nparticles,
- MPI_DOUBLE,
- MPI_SUM,
- this->comm);
- delete[] yy;
+ delete[] this->rhs[i];
}
- delete[] xg;
- delete[] xx;
}
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::get_rhs(double t, double *x, double *y)
+template <int particle_type, class rnumber, int interp_neighbours>
+void particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y)
{
switch(particle_type)
{
case VELOCITY_TRACER:
- this->sample_vec_field(this->vel, t, x, y, false);
+ this->vel->sample(this->nparticles, this->ncomponents, x, y);
break;
}
}
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::jump_estimate(double *jump)
-{
- std::fill_n(jump, this->nparticles, 0.0);
- switch(particle_type)
- {
- case VELOCITY_TRACER:
- double *y = new double[3*this->nparticles];
- this->sample_vec_field(this->vel, 1.0, this->state, y, false);
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
- {
- jump[p] = fabs(2*this->dt * y[3*p+2]);
- if (jump[p] < this->dz*1.01)
- jump[p] = this->dz*1.01;
- }
- delete[] y;
- break;
- }
-}
-
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-int particles<particle_type, rnumber, multistep, interp_neighbours>::get_rank(double z)
-{
- int tmp = this->fs->rd->rank[MOD(int(floor(z/this->dz)), this->fs->rd->sizes[0])];
- assert(tmp >= 0 && tmp < this->fs->rd->nprocs);
- return tmp;
-}
-
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::synchronize_single_particle_state(int p, double *x, int source)
-{
- if (!multistep)
- {
- if (source == -1) source = this->computing[p];
- if (this->watching[this->myrank*this->nparticles+p]) for (int r=0; r<this->fs->rd->nprocs; r++)
- if (r != source &&
- this->watching[r*this->nparticles+p])
- {
- //DEBUG_MSG("synchronizing state %d from %d to %d\n", p, this->computing[p], r);
- if (this->myrank == source)
- MPI_Send(
- x+p*this->ncomponents,
- this->ncomponents,
- MPI_DOUBLE,
- r,
- p+this->computing[p]*this->nparticles,
- this->comm);
- if (this->myrank == r)
- MPI_Recv(
- x+p*this->ncomponents,
- this->ncomponents,
- MPI_DOUBLE,
- source,
- p+this->computing[p]*this->nparticles,
- this->comm,
- MPI_STATUS_IGNORE);
- }
- }
-}
-
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::synchronize()
-{
- double *tstate = fftw_alloc_real(this->array_size);
- // first, synchronize state and jump across CPUs
- std::fill_n(tstate, this->array_size, 0.0);
- for (int p=0; p<this->nparticles; p++)
- {
- //if (this->watching[this->myrank*this->nparticles + p])
- //DEBUG_MSG(
- // "in synchronize, position for particle %d is %g %g %g\n",
- // p,
- // this->state[p*this->ncomponents],
- // this->state[p*this->ncomponents+1],
- // this->state[p*this->ncomponents+2]);
- if (this->myrank == this->computing[p])
- std::copy(this->state + p*this->ncomponents,
- this->state + (p+1)*this->ncomponents,
- tstate + p*this->ncomponents);
- }
- MPI_Allreduce(
- tstate,
- this->state,
- this->array_size,
- MPI_DOUBLE,
- MPI_SUM,
- this->comm);
- if (this->integration_steps >= 1 && multistep)
- {
- for (int i=0; i<this->integration_steps; i++)
- {
- std::fill_n(tstate, this->array_size, 0.0);
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
- std::copy(this->rhs[i] + p*this->ncomponents,
- this->rhs[i] + (p+1)*this->ncomponents,
- tstate + p*this->ncomponents);
- std::fill_n(this->rhs[i], this->array_size, 0.0);
- MPI_Allreduce(
- tstate,
- this->rhs[i],
- this->array_size,
- MPI_DOUBLE,
- MPI_SUM,
- this->comm);
- }
- }
- fftw_free(tstate);
- // assignment of particles
- for (int p=0; p<this->nparticles; p++)
- {
- this->computing[p] = this->get_rank(this->state[p*this->ncomponents + 2]);
- for (int r=0; r<this->nprocs; r++)
- this->watching[r*this->nparticles+p] = (r == this->computing[p]);
- //DEBUG_MSG("synchronizing particles, particle %d computing is %d\n", p, this->computing[p]);
- }
- if (!multistep)
- {
- double *jump = fftw_alloc_real(this->nparticles);
- this->jump_estimate(jump);
- // now, see who needs to watch
- bool *local_watching = new bool[this->fs->rd->nprocs*this->nparticles];
- std::fill_n(local_watching, this->fs->rd->nprocs*this->nparticles, false);
- for (int p=0; p<this->nparticles; p++)
- if (this->myrank == this->computing[p])
- {
- local_watching[this->get_rank(this->state[this->ncomponents*p+2] )*this->nparticles+p] = true;
- local_watching[this->get_rank(this->state[this->ncomponents*p+2]-jump[p])*this->nparticles+p] = true;
- local_watching[this->get_rank(this->state[this->ncomponents*p+2]+jump[p])*this->nparticles+p] = true;
- }
- fftw_free(jump);
- MPI_Allreduce(
- local_watching,
- this->watching,
- this->nparticles*this->fs->rd->nprocs,
- MPI_C_BOOL,
- MPI_LOR,
- this->comm);
- delete[] local_watching;
- }
-}
-
-
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::roll_rhs()
+template <int particle_type, class rnumber, int interp_neighbours>
+void particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
{
for (int i=this->integration_steps-2; i>=0; i--)
std::copy(this->rhs[i],
@@ -354,15 +97,16 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::roll_rhs()
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashforth(int nsteps)
+template <int particle_type, class rnumber, int interp_neighbours>
+void particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
+ const int nsteps)
{
ptrdiff_t ii;
- this->get_rhs(0, this->state, this->rhs[0]);
+ this->get_rhs(this->state, this->rhs[0]);
switch(nsteps)
{
case 1:
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
+ for (int p=0; p<this->nparticles; p++)
for (int i=0; i<this->ncomponents; i++)
{
ii = p*this->ncomponents+i;
@@ -370,7 +114,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashf
}
break;
case 2:
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
+ for (int p=0; p<this->nparticles; p++)
for (int i=0; i<this->ncomponents; i++)
{
ii = p*this->ncomponents+i;
@@ -379,7 +123,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashf
}
break;
case 3:
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
+ for (int p=0; p<this->nparticles; p++)
for (int i=0; i<this->ncomponents; i++)
{
ii = p*this->ncomponents+i;
@@ -389,7 +133,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashf
}
break;
case 4:
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
+ for (int p=0; p<this->nparticles; p++)
for (int i=0; i<this->ncomponents; i++)
{
ii = p*this->ncomponents+i;
@@ -400,7 +144,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashf
}
break;
case 5:
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
+ for (int p=0; p<this->nparticles; p++)
for (int i=0; i<this->ncomponents; i++)
{
ii = p*this->ncomponents+i;
@@ -412,7 +156,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashf
}
break;
case 6:
- for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
+ for (int p=0; p<this->nparticles; p++)
for (int i=0; i<this->ncomponents; i++)
{
ii = p*this->ncomponents+i;
@@ -429,284 +173,82 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashf
}
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::step()
+template <int particle_type, class rnumber, int interp_neighbours>
+void particles<particle_type, rnumber, interp_neighbours>::step()
{
- if (multistep)
- this->AdamsBashforth((this->iteration < this->integration_steps) ?
- this->iteration+1 :
- this->integration_steps);
- else
- this->Heun();
+ this->AdamsBashforth((this->iteration < this->integration_steps) ?
+ this->iteration+1 :
+ this->integration_steps);
this->iteration++;
- this->synchronize();
-}
-
-
-
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::Heun()
-{
- if (!multistep)
- {
- double *y = new double[this->array_size];
- double dtfactor[2];
- double factor[] = {0.0, 1.0};
- dtfactor[0] = factor[0]*this->dt;
- dtfactor[1] = factor[1]*this->dt;
- this->get_rhs(0, this->state, this->rhs[0]);
- for (int p=0; p<this->nparticles; p++)
- this->synchronize_single_particle_state(p, this->rhs[0]);
- for (int kindex = 1; kindex < 2; kindex++)
- {
- for (int p=0; p<this->nparticles; p++)
- {
- if (this->watching[this->myrank*this->nparticles+p])
- for (int i=0; i<this->ncomponents; i++)
- {
- ptrdiff_t tindex = ptrdiff_t(p)*this->ncomponents + i;
- y[tindex] = this->state[tindex] + dtfactor[kindex]*this->rhs[kindex-1][tindex];
- }
- }
- for (int p=0; p<this->nparticles; p++)
- this->synchronize_single_particle_state(p, y);
- this->get_rhs(factor[kindex], y, this->rhs[kindex]);
- for (int p=0; p<this->nparticles; p++)
- this->synchronize_single_particle_state(p, this->rhs[kindex]);
- }
- for (int p=0; p<this->nparticles; p++)
- {
- if (this->watching[this->myrank*this->nparticles+p])
- {
- for (int i=0; i<this->ncomponents; i++)
- {
- ptrdiff_t tindex = ptrdiff_t(p)*this->ncomponents + i;
- this->state[tindex] += this->dt*(this->rhs[0][tindex] + this->rhs[1][tindex])/2;
- }
- }
- }
- delete[] y;
- }
-}
-
-
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::cRK4()
-{
- if (!multistep)
- {
- static double factor[] = {0.0, 0.5, 0.5, 1.0};
- static double dtfactor[] = {0.0, this->dt/2, this->dt/2, this->dt};
- double *y = new double[this->array_size];
- this->get_rhs(0, this->state, this->rhs[0]);
- for (int p=0; p<this->nparticles; p++)
- this->synchronize_single_particle_state(p, this->rhs[0]);
- for (int kindex = 1; kindex < 4; kindex++)
- {
- for (int p=0; p<this->nparticles; p++)
- {
- if (this->watching[this->myrank*this->nparticles+p])
- for (int i=0; i<this->ncomponents; i++)
- {
- ptrdiff_t tindex = ptrdiff_t(p)*this->ncomponents + i;
- y[tindex] = this->state[tindex] + dtfactor[kindex]*this->rhs[kindex-1][tindex];
- }
- }
- for (int p=0; p<this->nparticles; p++)
- this->synchronize_single_particle_state(p, y);
- this->get_rhs(factor[kindex], y, this->rhs[kindex]);
- for (int p=0; p<this->nparticles; p++)
- this->synchronize_single_particle_state(p, this->rhs[kindex]);
- }
- for (int p=0; p<this->nparticles; p++)
- {
- if (this->watching[this->myrank*this->nparticles+p])
- for (int i=0; i<this->ncomponents; i++)
- {
- ptrdiff_t tindex = ptrdiff_t(p)*this->ncomponents + i;
- this->state[tindex] += this->dt*(this->rhs[0][tindex] +
- 2*(this->rhs[1][tindex] + this->rhs[2][tindex]) +
- this->rhs[3][tindex])/6;
- }
- }
- delete[] y;
- }
}
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::get_grid_coordinates(double *x, int *xg, double *xx)
-{
- static double grid_size[] = {this->dx, this->dy, this->dz};
- double tval;
- std::fill_n(xg, this->nparticles*3, 0);
- std::fill_n(xx, this->nparticles*3, 0.0);
- for (int p=0; p<this->nparticles; p++) if (this->watching[this->myrank*this->nparticles+p])
- {
- for (int c=0; c<3; c++)
- {
- tval = floor(x[p*this->ncomponents+c]/grid_size[c]);
- xg[p*3+c] = MOD(int(tval), this->fs->rd->sizes[2-c]);
- xx[p*3+c] = (x[p*this->ncomponents+c] - tval*grid_size[c]) / grid_size[c];
- }
- xg[p*3+2] -= this->fs->rd->starts[0];
- if (this->myrank == this->fs->rd->rank[0] &&
- xg[p*3+2] > this->fs->rd->subsizes[0])
- xg[p*3+2] -= this->fs->rd->sizes[0];
- //DEBUG_MSG(
- // "particle %d x is %lg %lg %lg xx is %lg %lg %lg xg is %d %d %d\n",
- // p,
- // x[p*3], x[p*3+1], x[p*3+2],
- // xx[p*3], xx[p*3+1], xx[p*3+2],
- // xg[p*3], xg[p*3+1], xg[p*3+2]);
- }
-}
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::read(hid_t data_file_id)
+template <int particle_type, class rnumber, int interp_neighbours>
+void particles<particle_type, rnumber, interp_neighbours>::read()
{
- //DEBUG_MSG("aloha\n");
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, rspace;
- hsize_t count[4], offset[4];
- rspace = H5Dget_space(dset);
- H5Sget_simple_extent_dims(rspace, 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(rspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, this->state);
- H5Sclose(mspace);
- H5Sclose(rspace);
- H5Dclose(dset);
- if (this->iteration > 0 && multistep)
+ for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
{
- temp_string = (std::string("/particles/") +
- std::string(this->name) +
- std::string("/rhs"));
- dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
- rspace = H5Dget_space(dset);
- H5Sget_simple_extent_dims(rspace, count, NULL);
- //reading from last available position
- offset[0] = count[0] - 1;
- offset[3] = 0;
- count[0] = 1;
- count[1] = 1;
- mspace = H5Screate_simple(4, count, NULL);
- for (int i=0; i<this->integration_steps; i++)
- {
- offset[1] = i;
- H5Sselect_hyperslab(rspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, this->rhs[i]);
- }
- H5Sclose(mspace);
- H5Sclose(rspace);
- H5Dclose(dset);
+ this->read_state_chunk(cindex, this->state+cindex*this->chunk_size*this->ncomponents);
+ if (this->iteration > 0)
+ for (int i=0; i<this->integration_steps; i++)
+ this->read_rhs_chunk(cindex, i, this->rhs[i]+cindex*this->chunk_size*this->ncomponents);
}
- }
MPI_Bcast(
this->state,
this->array_size,
MPI_DOUBLE,
0,
this->comm);
- if (multistep) for (int i = 0; i<this->integration_steps; i++)
- MPI_Bcast(
- this->rhs[i],
- this->array_size,
- MPI_DOUBLE,
- 0,
- this->comm);
- // initial assignment of particles
- for (int p=0; p<this->nparticles; p++)
- {
- this->computing[p] = this->get_rank(this->state[p*this->ncomponents + 2]);
- //DEBUG_MSG("reading particles, particle %d computing is %d\n", p, this->computing[p]);
- }
- // now actual synchronization
- this->synchronize();
+ if (this->iteration > 0)
+ for (int i = 0; i<this->integration_steps; i++)
+ MPI_Bcast(
+ this->rhs[i],
+ this->array_size,
+ MPI_DOUBLE,
+ 0,
+ this->comm);
}
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::write(hid_t data_file_id, bool write_rhs)
+template <int particle_type, class rnumber, int interp_neighbours>
+void particles<particle_type, rnumber, interp_neighbours>::write(
+ 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);
- if (write_rhs && multistep)
+ for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
{
- 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);
- H5Sget_simple_extent_dims(wspace, count, NULL);
- //writing to last available position
- offset[0] = count[0] - 1;
- count[0] = 1;
- count[1] = 1;
- offset[3] = 0;
- mspace = H5Screate_simple(4, count, NULL);
- for (int i=0; i<this->integration_steps; i++)
- {
- offset[1] = i;
- H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, this->rhs[i]);
- }
- H5Sclose(mspace);
- H5Sclose(wspace);
- H5Dclose(dset);
+ this->write_state_chunk(cindex, this->state+cindex*this->chunk_size*this->ncomponents);
+ if (write_rhs)
+ for (int i=0; i<this->integration_steps; i++)
+ this->write_rhs_chunk(cindex, i, this->rhs[i]+cindex*this->chunk_size*this->ncomponents);
}
- }
+}
+
+template <int particle_type, class rnumber, int interp_neighbours>
+void particles<particle_type, rnumber, interp_neighbours>::sample(
+ interpolator_base<rnumber, interp_neighbours> *field,
+ const char *dset_name)
+{
+ double *y = new double[this->nparticles*3];
+ field->sample(this->nparticles, this->ncomponents, this->state, y);
+ if (this->myrank == 0)
+ for (int cindex=0; cindex<this->get_number_of_chunks(); cindex++)
+ this->write_point3D_chunk(dset_name, cindex, y+cindex*this->chunk_size*3);
+ delete[] y;
}
/*****************************************************************************/
-template class particles<VELOCITY_TRACER, float, true, 1>;
-template class particles<VELOCITY_TRACER, float, true, 2>;
-template class particles<VELOCITY_TRACER, float, true, 3>;
-template class particles<VELOCITY_TRACER, float, true, 4>;
-template class particles<VELOCITY_TRACER, float, true, 5>;
-template class particles<VELOCITY_TRACER, float, true, 6>;
-template class particles<VELOCITY_TRACER, float, false, 1>;
-template class particles<VELOCITY_TRACER, float, false, 2>;
-template class particles<VELOCITY_TRACER, float, false, 3>;
-template class particles<VELOCITY_TRACER, float, false, 4>;
-template class particles<VELOCITY_TRACER, float, false, 5>;
-template class particles<VELOCITY_TRACER, float, false, 6>;
-template class particles<VELOCITY_TRACER, double, true, 1>;
-template class particles<VELOCITY_TRACER, double, true, 2>;
-template class particles<VELOCITY_TRACER, double, true, 3>;
-template class particles<VELOCITY_TRACER, double, true, 4>;
-template class particles<VELOCITY_TRACER, double, true, 5>;
-template class particles<VELOCITY_TRACER, double, true, 6>;
-template class particles<VELOCITY_TRACER, double, false, 1>;
-template class particles<VELOCITY_TRACER, double, false, 2>;
-template class particles<VELOCITY_TRACER, double, false, 3>;
-template class particles<VELOCITY_TRACER, double, false, 4>;
-template class particles<VELOCITY_TRACER, double, false, 5>;
-template class particles<VELOCITY_TRACER, double, false, 6>;
+template class particles<VELOCITY_TRACER, float, 1>;
+template class particles<VELOCITY_TRACER, float, 2>;
+template class particles<VELOCITY_TRACER, float, 3>;
+template class particles<VELOCITY_TRACER, float, 4>;
+template class particles<VELOCITY_TRACER, float, 5>;
+template class particles<VELOCITY_TRACER, float, 6>;
+template class particles<VELOCITY_TRACER, double, 1>;
+template class particles<VELOCITY_TRACER, double, 2>;
+template class particles<VELOCITY_TRACER, double, 3>;
+template class particles<VELOCITY_TRACER, double, 4>;
+template class particles<VELOCITY_TRACER, double, 5>;
+template class particles<VELOCITY_TRACER, double, 6>;
/*****************************************************************************/
diff --git a/bfps/cpp/particles.hpp b/bfps/cpp/particles.hpp
index 6fa379c458f38d364765efa4d7e58bb25e9daff5..cf67b1fc32509b16d7a0774c0fc40a74759e39e0 100644
--- a/bfps/cpp/particles.hpp
+++ b/bfps/cpp/particles.hpp
@@ -31,110 +31,68 @@
#include "base.hpp"
#include "particles_base.hpp"
#include "fluid_solver_base.hpp"
-#include "interpolator.hpp"
+#include "interpolator_base.hpp"
#ifndef PARTICLES
#define PARTICLES
-template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-class particles
+template <int particle_type, class rnumber, int interp_neighbours>
+class particles: public particles_io_base<particle_type>
{
- public:
- int myrank, nprocs;
- MPI_Comm comm;
- fluid_solver_base<rnumber> *fs;
- rnumber *data;
-
- /* watching is an array of shape [nparticles], with
- * watching[p] being true if particle p is in the domain of myrank
- * or in the buffer regions.
- * only used if multistep is false.
- * */
- bool *watching;
- /* computing is an array of shape [nparticles], with
- * computing[p] being the rank that is currently working on particle p
- * */
- int *computing;
-
- /* state will generally hold all the information about the particles.
- * in the beginning, we will only need to solve 3D ODEs, but I figured
- * a general ncomponents is better, since we may change our minds.
- * */
+ private:
double *state;
double *rhs[6];
- int nparticles;
- int ncomponents;
+
+ public:
int array_size;
int integration_steps;
- int traj_skip;
- int buffer_width;
- ptrdiff_t buffer_size;
- double *lbound;
- double *ubound;
- interpolator<rnumber, interp_neighbours> *vel;
+ interpolator_base<rnumber, interp_neighbours> *vel;
/* simulation parameters */
- char name[256];
- int iteration;
double dt;
- /* physical parameters of field */
- double dx, dy, dz;
-
/* methods */
/* constructor and destructor.
* allocate and deallocate:
* this->state
* this->rhs
- * this->lbound
- * this->ubound
- * this->computing
- * this->watching
* */
particles(
const char *NAME,
- fluid_solver_base<rnumber> *FSOLVER,
- interpolator<rnumber, interp_neighbours> *FIELD,
- const int NPARTICLES,
+ const hid_t data_file_id,
+ interpolator_base<rnumber, interp_neighbours> *FIELD,
const int TRAJ_SKIP,
const int INTEGRATION_STEPS = 2);
~particles();
- /* an Euler step is needed to compute an estimate of future positions,
- * which is needed for synchronization.
- * */
- void jump_estimate(double *__restrict__ jump_length);
- void get_rhs(double *__restrict__ x, double *__restrict__ rhs);
- void get_rhs(double t, double *__restrict__ x, double *__restrict__ rhs);
-
- int get_rank(double z); // get rank for given value of z
- void synchronize();
- void synchronize_single_particle_state(int p, double *__restrict__ x, int source_id = -1);
- void get_grid_coordinates(double *__restrict__ x, int *__restrict__ xg, double *__restrict__ xx);
- void sample_vec_field(
- interpolator<rnumber, interp_neighbours> *vec,
- double t,
- double *__restrict__ x,
- double *__restrict__ y,
- const bool synch = false,
- int *deriv = NULL);
- inline void sample_vec_field(interpolator<rnumber, interp_neighbours> *field, double *vec_values)
+ void sample(
+ interpolator_base<rnumber, interp_neighbours> *field,
+ const char *dset_name);
+
+ inline void sample(
+ interpolator_base<rnumber, interp_neighbours> *field,
+ double *y)
{
- this->sample_vec_field(field, 1.0, this->state, vec_values, true, NULL);
+ field->sample(this->nparticles, this->ncomponents, this->state, y);
}
+ void get_rhs(
+ double *__restrict__ x,
+ double *__restrict__ rhs);
+
/* input/output */
- void read(hid_t data_file_id);
- void write(hid_t data_file_id, bool write_rhs = true);
+ void read();
+ void write(
+ const char *dset_name,
+ const double *data);
+ void write(const bool write_rhs = true);
/* solvers */
void step();
void roll_rhs();
- void AdamsBashforth(int nsteps);
- void Heun();
- void cRK4();
+ void AdamsBashforth(const int nsteps);
};
#endif//PARTICLES
diff --git a/bfps/cpp/particles_base.cpp b/bfps/cpp/particles_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c0fc631799cf55e1fbd9b8c93db945f47aa9e980
--- /dev/null
+++ b/bfps/cpp/particles_base.cpp
@@ -0,0 +1,430 @@
+/**********************************************************************
+* *
+* Copyright 2015 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 *
+* *
+**********************************************************************/
+
+
+
+#define NDEBUG
+
+#include <algorithm>
+#include <cassert>
+#include "particles_base.hpp"
+
+template <int particle_type>
+single_particle_state<particle_type>::single_particle_state()
+{
+ switch(particle_type)
+ {
+ default:
+ this->data = new double[3];
+ std::fill_n(this->data, 3, 0);
+ break;
+ }
+}
+
+template <int particle_type>
+single_particle_state<particle_type>::single_particle_state(
+ const single_particle_state<particle_type> &src)
+{
+ switch(particle_type)
+ {
+ default:
+ this->data = new double[3];
+ std::copy(src.data, src.data + 3, this->data);
+ break;
+ }
+}
+
+template <int particle_type>
+single_particle_state<particle_type>::single_particle_state(
+ const double *src)
+{
+ switch(particle_type)
+ {
+ default:
+ this->data = new double[3];
+ std::copy(src, src + 3, this->data);
+ break;
+ }
+}
+
+template <int particle_type>
+single_particle_state<particle_type>::~single_particle_state()
+{
+ switch(particle_type)
+ {
+ default:
+ delete[] this->data;
+ break;
+ }
+}
+
+template <int particle_type>
+single_particle_state<particle_type> &single_particle_state<particle_type>::operator=(
+ const single_particle_state &src)
+{
+ switch(particle_type)
+ {
+ default:
+ std::copy(src.data, src.data + 3, this->data);
+ break;
+ }
+ return *this;
+}
+
+template <int particle_type>
+single_particle_state<particle_type> &single_particle_state<particle_type>::operator=(
+ const double *src)
+{
+ switch(particle_type)
+ {
+ default:
+ std::copy(src, src + 3, this->data);
+ break;
+ }
+ return *this;
+}
+
+int get_chunk_offsets(
+ std::vector<hsize_t> data_dims,
+ std::vector<hsize_t> chnk_dims,
+ std::vector<std::vector<hsize_t>> &co)
+{
+ std::vector<hsize_t> nchunks(data_dims);
+ int total_number_of_chunks = 1;
+ for (unsigned i=0; i<nchunks.size(); i++)
+ {
+ nchunks[i] = data_dims[i] / chnk_dims[i];
+ total_number_of_chunks *= nchunks[i];
+ }
+ co.resize(total_number_of_chunks);
+ DEBUG_MSG("total number of chunks is %d\n", total_number_of_chunks);
+ for (int cindex=0; cindex < total_number_of_chunks; cindex++)
+ {
+ int cc = cindex;
+ for (unsigned i=0; i<nchunks.size(); i++)
+ {
+ int ii = nchunks.size()-1-i;
+ co[cindex].resize(nchunks.size());
+ co[cindex][ii] = cc % nchunks[ii];
+ cc = (cc - co[cindex][ii]) / nchunks[ii];
+ co[cindex][ii] *= chnk_dims[ii];
+ }
+ }
+ return EXIT_SUCCESS;
+}
+
+template <int particle_type>
+particles_io_base<particle_type>::particles_io_base(
+ const char *NAME,
+ const int TRAJ_SKIP,
+ const hid_t data_file_id,
+ MPI_Comm COMM)
+{
+ switch(particle_type)
+ {
+ default:
+ this->ncomponents = 3;
+ break;
+ }
+ this->name = std::string(NAME);
+ this->traj_skip = TRAJ_SKIP;
+ this->comm = COMM;
+ MPI_Comm_rank(COMM, &this->myrank);
+ MPI_Comm_size(COMM, &this->nprocs);
+
+ if (this->myrank == 0)
+ {
+ hid_t dset, prop_list, dspace;
+ this->hdf5_group_id = H5Gopen(data_file_id, this->name.c_str(), H5P_DEFAULT);
+ dset = H5Dopen(this->hdf5_group_id, "state", H5P_DEFAULT);
+ dspace = H5Dget_space(dset);
+ this->hdf5_state_dims.resize(H5Sget_simple_extent_ndims(dspace));
+ H5Sget_simple_extent_dims(dspace, &this->hdf5_state_dims.front(), NULL);
+ assert(this->hdf5_state_dims[this->hdf5_state_dims.size()-1] == this->ncomponents);
+ this->nparticles = 1;
+ for (int i=1; i<this->hdf5_state_dims.size()-1; i++)
+ this->nparticles *= this->hdf5_state_dims[i];
+ prop_list = H5Dget_create_plist(dset);
+ this->hdf5_state_chunks.resize(this->hdf5_state_dims.size());
+ H5Pget_chunk(prop_list, this->hdf5_state_dims.size(), &this->hdf5_state_chunks.front());
+ H5Pclose(prop_list);
+ H5Sclose(dspace);
+ H5Dclose(dset);
+ this->chunk_size = 1;
+ for (auto i=1; i<this->hdf5_state_dims.size()-1; i++)
+ this->chunk_size *= this->hdf5_state_chunks[i];
+ dset = H5Dopen(this->hdf5_group_id, "rhs", H5P_DEFAULT);
+ dspace = H5Dget_space(dset);
+ this->hdf5_rhs_dims.resize(H5Sget_simple_extent_ndims(dspace));
+ H5Sget_simple_extent_dims(dspace, &this->hdf5_rhs_dims.front(), NULL);
+ prop_list = H5Dget_create_plist(dset);
+ this->hdf5_rhs_chunks.resize(this->hdf5_rhs_dims.size());
+ H5Pget_chunk(prop_list, this->hdf5_rhs_dims.size(), &this->hdf5_rhs_chunks.front());
+ H5Pclose(prop_list);
+ H5Sclose(dspace);
+ H5Dclose(dset);
+ }
+
+ int tmp;
+ tmp = this->hdf5_state_dims.size();
+ MPI_Bcast(
+ &tmp,
+ 1,
+ MPI_INTEGER,
+ 0,
+ this->comm);
+ if (this->myrank != 0)
+ {
+ this->hdf5_state_dims.resize(tmp);
+ this->hdf5_state_chunks.resize(tmp);
+ }
+ MPI_Bcast(
+ &this->hdf5_state_dims.front(),
+ this->hdf5_state_dims.size(),
+ MPI_INTEGER,
+ 0,
+ this->comm);
+ MPI_Bcast(
+ &this->hdf5_state_chunks.front(),
+ this->hdf5_state_chunks.size(),
+ MPI_INTEGER,
+ 0,
+ this->comm);
+ std::vector<hsize_t> tdims(this->hdf5_state_dims), tchnk(this->hdf5_state_chunks);
+ tdims.erase(tdims.begin()+0);
+ tchnk.erase(tchnk.begin()+0);
+ tdims.erase(tdims.end()-1);
+ tchnk.erase(tchnk.end()-1);
+ get_chunk_offsets(tdims, tchnk, this->chunk_offsets);
+ MPI_Bcast(
+ &this->chunk_size,
+ 1,
+ MPI_INTEGER,
+ 0,
+ this->comm);
+ MPI_Bcast(
+ &this->nparticles,
+ 1,
+ MPI_INTEGER,
+ 0,
+ this->comm);
+ DEBUG_MSG("nparticles = %d, chunk_size = %d\n",
+ this->nparticles,
+ this->chunk_size);
+ DEBUG_MSG("exiting particles_io_base constructor\n");
+}
+
+template <int particle_type>
+particles_io_base<particle_type>::~particles_io_base()
+{
+ if(this->myrank == 0)
+ H5Gclose(this->hdf5_group_id);
+}
+
+template <int particle_type>
+void particles_io_base<particle_type>::read_state_chunk(
+ const int cindex,
+ double *data)
+{
+ DEBUG_MSG("entered read_state_chunk\n");
+ hid_t dset = H5Dopen(this->hdf5_group_id, "state", H5P_DEFAULT);
+ hid_t rspace = H5Dget_space(dset);
+ std::vector<hsize_t> mem_dims(this->hdf5_state_chunks);
+ mem_dims[0] = 1;
+ hid_t mspace = H5Screate_simple(
+ this->hdf5_state_dims.size(),
+ &mem_dims.front(),
+ NULL);
+ hsize_t *offset = new hsize_t[this->hdf5_state_dims.size()];
+ offset[0] = this->iteration / this->traj_skip;
+ for (int i=1; i<this->hdf5_state_dims.size()-1; i++)
+ offset[i] = this->chunk_offsets[cindex][i-1];
+ offset[this->hdf5_state_dims.size()-1] = 0;
+ H5Sselect_hyperslab(
+ rspace,
+ H5S_SELECT_SET,
+ offset,
+ NULL,
+ &mem_dims.front(),
+ NULL);
+ H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, data);
+ H5Sclose(mspace);
+ H5Sclose(rspace);
+ H5Dclose(dset);
+ delete[] offset;
+ DEBUG_MSG("exiting read_state_chunk\n");
+}
+
+template <int particle_type>
+void particles_io_base<particle_type>::write_state_chunk(
+ const int cindex,
+ const double *data)
+{
+ hid_t dset = H5Dopen(this->hdf5_group_id, "state", H5P_DEFAULT);
+ hid_t rspace = H5Dget_space(dset);
+ std::vector<hsize_t> mem_dims(this->hdf5_state_chunks);
+ mem_dims[0] = 1;
+ hid_t mspace = H5Screate_simple(
+ this->hdf5_state_dims.size(),
+ &mem_dims.front(),
+ NULL);
+ hsize_t *offset = new hsize_t[this->hdf5_state_dims.size()];
+ offset[0] = this->iteration / this->traj_skip;
+ for (int i=1; i<this->hdf5_state_dims.size()-1; i++)
+ offset[i] = this->chunk_offsets[cindex][i-1];
+ offset[this->hdf5_state_dims.size()-1] = 0;
+ H5Sselect_hyperslab(
+ rspace,
+ H5S_SELECT_SET,
+ offset,
+ NULL,
+ &mem_dims.front(),
+ NULL);
+ H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, data);
+ H5Sclose(mspace);
+ H5Sclose(rspace);
+ H5Dclose(dset);
+ delete[] offset;
+}
+
+template <int particle_type>
+void particles_io_base<particle_type>::read_rhs_chunk(
+ const int cindex,
+ const int rhsindex,
+ double *data)
+{
+ //DEBUG_MSG("entered read_rhs_chunk\n");
+ hid_t dset = H5Dopen(this->hdf5_group_id, "rhs", H5P_DEFAULT);
+ hid_t rspace = H5Dget_space(dset);
+ std::vector<hsize_t> mem_dims(this->hdf5_rhs_chunks);
+ mem_dims[0] = 1;
+ mem_dims[1] = 1;
+ hid_t mspace = H5Screate_simple(
+ this->hdf5_rhs_dims.size(),
+ &mem_dims.front(),
+ NULL);
+ hsize_t *offset = new hsize_t[this->hdf5_rhs_dims.size()];
+ offset[0] = this->hdf5_rhs_dims[0]-1;
+ offset[1] = rhsindex;
+ for (int i=2; i<this->hdf5_rhs_dims.size()-1; i++)
+ offset[i] = this->chunk_offsets[cindex][i-2];
+ offset[this->hdf5_rhs_dims.size()-1] = 0;
+ //for (int i=0; i<this->hdf5_rhs_dims.size(); i++)
+ // DEBUG_MSG("rhs dim %d: size=%d chunk=%d offset=%d\n",
+ // i, this->hdf5_rhs_dims[i], this->hdf5_rhs_chunks[i], offset[i]);
+ H5Sselect_hyperslab(
+ rspace,
+ H5S_SELECT_SET,
+ offset,
+ NULL,
+ &mem_dims.front(),
+ NULL);
+ //DEBUG_MSG("selected hyperslab\n");
+ H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, data);
+ //DEBUG_MSG("data has been read\n");
+ H5Sclose(mspace);
+ H5Sclose(rspace);
+ H5Dclose(dset);
+ delete[] offset;
+ //DEBUG_MSG("exiting read_rhs_chunk\n");
+}
+
+template <int particle_type>
+void particles_io_base<particle_type>::write_rhs_chunk(
+ const int cindex,
+ const int rhsindex,
+ const double *data)
+{
+ hid_t dset = H5Dopen(this->hdf5_group_id, "rhs", H5P_DEFAULT);
+ hid_t rspace = H5Dget_space(dset);
+ std::vector<hsize_t> mem_dims(this->hdf5_rhs_chunks);
+ mem_dims[0] = 1;
+ mem_dims[1] = 1;
+ hid_t mspace = H5Screate_simple(
+ this->hdf5_rhs_dims.size(),
+ &mem_dims.front(),
+ NULL);
+ hsize_t *offset = new hsize_t[this->hdf5_rhs_dims.size()];
+ offset[0] = this->hdf5_rhs_dims[0];
+ offset[1] = rhsindex;
+ for (int i=2; i<this->hdf5_rhs_dims.size()-1; i++)
+ offset[i] = this->chunk_offsets[cindex][i-2];
+ offset[this->hdf5_rhs_dims.size()-1] = 0;
+ DEBUG_MSG("rhs write offsets are %d %d %d %d\n",
+ offset[0], offset[1], offset[2], offset[3]);
+ H5Sselect_hyperslab(
+ rspace,
+ H5S_SELECT_SET,
+ offset,
+ NULL,
+ &mem_dims.front(),
+ NULL);
+ H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, data);
+ H5Sclose(mspace);
+ H5Sclose(rspace);
+ H5Dclose(dset);
+ delete[] offset;
+}
+
+template <int particle_type>
+void particles_io_base<particle_type>::write_point3D_chunk(
+ const std::string dset_name,
+ const int cindex,
+ const double *data)
+{
+ hid_t dset = H5Dopen(this->hdf5_group_id, dset_name.c_str(), H5P_DEFAULT);
+ hid_t rspace = H5Dget_space(dset);
+ std::vector<hsize_t> mem_dims(this->hdf5_state_chunks);
+ mem_dims[0] = 1;
+ mem_dims[mem_dims.size()-1] = 3;
+ hid_t mspace = H5Screate_simple(
+ this->hdf5_state_dims.size(),
+ &mem_dims.front(),
+ NULL);
+ hsize_t *offset = new hsize_t[this->hdf5_state_dims.size()];
+ offset[0] = this->iteration / this->traj_skip;
+ for (int i=1; i<this->hdf5_state_dims.size()-1; i++)
+ offset[i] = this->chunk_offsets[cindex][i-1];
+ offset[this->hdf5_state_dims.size()-1] = 0;
+ H5Sselect_hyperslab(
+ rspace,
+ H5S_SELECT_SET,
+ offset,
+ NULL,
+ &mem_dims.front(),
+ NULL);
+ H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, data);
+ H5Sclose(mspace);
+ H5Sclose(rspace);
+ H5Dclose(dset);
+ delete[] offset;
+}
+
+/*****************************************************************************/
+template class single_particle_state<POINT3D>;
+template class single_particle_state<VELOCITY_TRACER>;
+
+template class particles_io_base<VELOCITY_TRACER>;
+/*****************************************************************************/
+
diff --git a/bfps/cpp/particles_base.hpp b/bfps/cpp/particles_base.hpp
index cf61f38efadf0030a67a10ff94c21a92e79bdd64..923c29253ba1ddabff1791ca40d03f79109eaf20 100644
--- a/bfps/cpp/particles_base.hpp
+++ b/bfps/cpp/particles_base.hpp
@@ -24,6 +24,8 @@
+#include <vector>
+#include <hdf5.h>
#include "interpolator_base.hpp"
#ifndef PARTICLES_BASE
@@ -31,7 +33,96 @@
#define PARTICLES_BASE
/* particle types */
-enum particle_types {VELOCITY_TRACER};
+enum particle_types {POINT3D, VELOCITY_TRACER};
+
+/* 1 particle state type */
+
+template <int particle_type>
+class single_particle_state
+{
+ public:
+ double *data;
+
+ single_particle_state();
+ single_particle_state(const single_particle_state &src);
+ single_particle_state(const double *src);
+ ~single_particle_state();
+
+ single_particle_state<particle_type> &operator=(const single_particle_state &src);
+ single_particle_state<particle_type> &operator=(const double *src);
+
+ inline double &operator[](const int i)
+ {
+ return this->data[i];
+ }
+};
+
+std::vector<std::vector<hsize_t>> get_chunk_offsets(
+ std::vector<hsize_t> data_dims,
+ std::vector<hsize_t> chnk_dims);
+
+template <int particle_type>
+class particles_io_base
+{
+ protected:
+ int myrank, nprocs;
+ MPI_Comm comm;
+
+ int nparticles;
+ int ncomponents;
+
+ std::string name;
+ int chunk_size;
+ int traj_skip;
+
+ hid_t hdf5_group_id;
+ std::vector<hsize_t> hdf5_state_dims, hdf5_state_chunks;
+ std::vector<hsize_t> hdf5_rhs_dims, hdf5_rhs_chunks;
+
+ std::vector<std::vector<hsize_t>> chunk_offsets;
+
+ particles_io_base(
+ const char *NAME,
+ const int TRAJ_SKIP,
+ const hid_t data_file_id,
+ MPI_Comm COMM);
+ virtual ~particles_io_base();
+
+ void read_state_chunk(
+ const int cindex,
+ double *__restrict__ data);
+ void write_state_chunk(
+ const int cindex,
+ const double *data);
+ void read_rhs_chunk(
+ const int cindex,
+ const int rhsindex,
+ double *__restrict__ data);
+ void write_rhs_chunk(
+ const int cindex,
+ const int rhsindex,
+ const double *data);
+
+ void write_point3D_chunk(
+ const std::string dset_name,
+ const int cindex,
+ const double *data);
+
+ public:
+ int iteration;
+
+ inline const char *get_name()
+ {
+ return this->name.c_str();
+ }
+ inline const unsigned int get_number_of_chunks()
+ {
+ return this->chunk_offsets.size();
+ }
+ inline const unsigned int get_number_of_rhs_chunks();
+ virtual void read() = 0;
+ virtual void write(const bool write_rhs = true) = 0;
+};
#endif//PARTICLES_BASE
diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index c701fbd64bba8e1c68ab2de953eae0795baead2f..facc09c42164d44478ffa01b47354a05b0decce0 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -33,17 +33,11 @@ template <class rnumber, int interp_neighbours>
rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
fluid_solver_base<rnumber> *fs,
base_polynomial_values BETA_POLYS,
- rnumber *FIELD)
+ rnumber *FIELD) : interpolator_base<rnumber, interp_neighbours>(fs, BETA_POLYS)
{
- this->descriptor = fs->rd;
this->field_size = 2*fs->cd->local_size;
- this->compute_beta = BETA_POLYS;
this->field = FIELD;
- // compute dx, dy, dz;
- this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
- this->dy = 4*acos(0) / (fs->dky*this->descriptor->sizes[1]);
- this->dz = 4*acos(0) / (fs->dkz*this->descriptor->sizes[0]);
// generate compute array
this->compute = new bool[this->descriptor->sizes[0]];
@@ -60,29 +54,6 @@ rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
delete[] this->compute;
}
-template <class rnumber, int interp_neighbours>
-void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
- const int nparticles,
- const int pdimension,
- const double *x,
- int *xg,
- double *xx)
-{
- static double grid_size[] = {this->dx, this->dy, this->dz};
- double tval;
- std::fill_n(xg, nparticles*3, 0);
- std::fill_n(xx, nparticles*3, 0.0);
- for (int p=0; p<nparticles; p++)
- {
- for (int c=0; c<3; c++)
- {
- tval = floor(x[p*pdimension+c]/grid_size[c]);
- xg[p*3+c] = MOD(int(tval), this->descriptor->sizes[2-c]);
- xx[p*3+c] = (x[p*pdimension+c] - tval*grid_size[c]) / grid_size[c];
- }
- }
-}
-
template <class rnumber, int interp_neighbours>
void rFFTW_interpolator<rnumber, interp_neighbours>::sample(
const int nparticles,
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
index 3dccd7150d17bea612128bf6d70b37033d281db8..fb96963eee8dd1361a59edf00b38b6460a27d9e9 100644
--- a/bfps/cpp/rFFTW_interpolator.hpp
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -34,25 +34,16 @@
#define RFFTW_INTERPOLATOR
template <class rnumber, int interp_neighbours>
-class rFFTW_interpolator
+class rFFTW_interpolator:public interpolator_base<rnumber, interp_neighbours>
{
public:
/* size of field to interpolate */
ptrdiff_t field_size;
- /* pointer to polynomial function */
- base_polynomial_values compute_beta;
-
- /* descriptor of field to interpolate */
- field_descriptor<rnumber> *descriptor;
-
/* pointers to fields that are to be interpolated
* */
rnumber *field;
- /* physical parameters of field */
- double dx, dy, dz;
-
/* compute[iz] is true if .
* local_zstart - neighbours <= iz <= local_zend + 1 + neighbours
* */
@@ -64,13 +55,13 @@ class rFFTW_interpolator
rnumber *FIELD_DATA);
~rFFTW_interpolator();
- /* map real locations to grid coordinates */
- void get_grid_coordinates(
- const int nparticles,
- const int pdimension,
- const double *__restrict__ x,
- int *__restrict__ xg,
- double *__restrict__ xx);
+ /* does not destroy input */
+ inline int read_rFFTW(const void *src)
+ {
+ this->field = (rnumber*)src;
+ return EXIT_SUCCESS;
+ }
+
/* interpolate field at an array of locations */
void sample(
const int nparticles,
diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp
deleted file mode 100644
index a71c0574326340dcdd80ee9ca4a29278b62e4a12..0000000000000000000000000000000000000000
--- a/bfps/cpp/rFFTW_particles.cpp
+++ /dev/null
@@ -1,340 +0,0 @@
-/**********************************************************************
-* *
-* Copyright 2015 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 *
-* *
-**********************************************************************/
-
-
-
-#define NDEBUG
-
-#include <cmath>
-#include <cassert>
-#include <cstring>
-#include <string>
-#include <sstream>
-
-#include "base.hpp"
-#include "rFFTW_particles.hpp"
-#include "fftw_tools.hpp"
-
-
-extern int myrank, nprocs;
-
-template <int particle_type, class rnumber, int interp_neighbours>
-rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
- const char *NAME,
- rFFTW_interpolator<rnumber, interp_neighbours> *FIELD,
- const int NPARTICLES,
- const int TRAJ_SKIP,
- const int INTEGRATION_STEPS)
-{
- switch(particle_type)
- {
- case VELOCITY_TRACER:
- this->ncomponents = 3;
- break;
- default:
- this->ncomponents = 3;
- break;
- }
- assert((INTEGRATION_STEPS <= 6) &&
- (INTEGRATION_STEPS >= 1));
- strncpy(this->name, NAME, 256);
- this->nparticles = NPARTICLES;
- this->vel = FIELD;
- this->integration_steps = INTEGRATION_STEPS;
- this->traj_skip = TRAJ_SKIP;
- this->myrank = this->vel->descriptor->myrank;
- this->nprocs = this->vel->descriptor->nprocs;
- this->comm = this->vel->descriptor->comm;
- this->array_size = this->nparticles * this->ncomponents;
- this->state = new double[this->array_size];
- std::fill_n(this->state, this->array_size, 0.0);
- for (int i=0; i < this->integration_steps; i++)
- {
- this->rhs[i] = new double[this->array_size];
- std::fill_n(this->rhs[i], this->array_size, 0.0);
- }
-}
-
-template <int particle_type, class rnumber, int interp_neighbours>
-rFFTW_particles<particle_type, rnumber, interp_neighbours>::~rFFTW_particles()
-{
- delete[] this->state;
- for (int i=0; i < this->integration_steps; i++)
- {
- delete[] this->rhs[i];
- }
-}
-
-template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y)
-{
- switch(particle_type)
- {
- case VELOCITY_TRACER:
- this->vel->sample(this->nparticles, this->ncomponents, x, y);
- break;
- }
-}
-
-template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
-{
- for (int i=this->integration_steps-2; i>=0; i--)
- std::copy(this->rhs[i],
- this->rhs[i] + this->array_size,
- this->rhs[i+1]);
-}
-
-
-
-template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(
- const int nsteps)
-{
- ptrdiff_t ii;
- this->get_rhs(this->state, this->rhs[0]);
- switch(nsteps)
- {
- case 1:
- for (int p=0; p<this->nparticles; p++)
- for (int i=0; i<this->ncomponents; i++)
- {
- ii = p*this->ncomponents+i;
- this->state[ii] += this->dt*this->rhs[0][ii];
- }
- break;
- case 2:
- for (int p=0; p<this->nparticles; p++)
- for (int i=0; i<this->ncomponents; i++)
- {
- ii = p*this->ncomponents+i;
- this->state[ii] += this->dt*(3*this->rhs[0][ii]
- - this->rhs[1][ii])/2;
- }
- break;
- case 3:
- for (int p=0; p<this->nparticles; p++)
- for (int i=0; i<this->ncomponents; i++)
- {
- ii = p*this->ncomponents+i;
- this->state[ii] += this->dt*(23*this->rhs[0][ii]
- - 16*this->rhs[1][ii]
- + 5*this->rhs[2][ii])/12;
- }
- break;
- case 4:
- for (int p=0; p<this->nparticles; p++)
- for (int i=0; i<this->ncomponents; i++)
- {
- ii = p*this->ncomponents+i;
- this->state[ii] += this->dt*(55*this->rhs[0][ii]
- - 59*this->rhs[1][ii]
- + 37*this->rhs[2][ii]
- - 9*this->rhs[3][ii])/24;
- }
- break;
- case 5:
- for (int p=0; p<this->nparticles; p++)
- for (int i=0; i<this->ncomponents; i++)
- {
- ii = p*this->ncomponents+i;
- this->state[ii] += this->dt*(1901*this->rhs[0][ii]
- - 2774*this->rhs[1][ii]
- + 2616*this->rhs[2][ii]
- - 1274*this->rhs[3][ii]
- + 251*this->rhs[4][ii])/720;
- }
- break;
- case 6:
- for (int p=0; p<this->nparticles; p++)
- for (int i=0; i<this->ncomponents; i++)
- {
- ii = p*this->ncomponents+i;
- this->state[ii] += this->dt*(4277*this->rhs[0][ii]
- - 7923*this->rhs[1][ii]
- + 9982*this->rhs[2][ii]
- - 7298*this->rhs[3][ii]
- + 2877*this->rhs[4][ii]
- - 475*this->rhs[5][ii])/1440;
- }
- break;
- }
- this->roll_rhs();
-}
-
-
-template <int particle_type, class rnumber, int interp_neighbours>
-void rFFTW_particles<particle_type, rnumber, interp_neighbours>::step()
-{
- this->AdamsBashforth((this->iteration < this->integration_steps) ?
- this->iteration+1 :
- this->integration_steps);
- this->iteration++;
-}
-
-
-template <int particle_type, class rnumber, int interp_neighbours>
-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("/") +
- std::string(this->name) +
- std::string("/state"));
- hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
- hid_t mspace, rspace;
- hsize_t count[4], offset[4];
- rspace = H5Dget_space(dset);
- H5Sget_simple_extent_dims(rspace, 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(rspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, this->state);
- H5Sclose(mspace);
- H5Sclose(rspace);
- H5Dclose(dset);
- if (this->iteration > 0)
- {
- temp_string = (std::string("/") +
- std::string(this->name) +
- std::string("/rhs"));
- dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT);
- rspace = H5Dget_space(dset);
- H5Sget_simple_extent_dims(rspace, count, NULL);
- //reading from last available position
- offset[0] = count[0] - 1;
- offset[3] = 0;
- count[0] = 1;
- count[1] = 1;
- mspace = H5Screate_simple(4, count, NULL);
- for (int i=0; i<this->integration_steps; i++)
- {
- offset[1] = i;
- H5Sselect_hyperslab(rspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, this->rhs[i]);
- }
- H5Sclose(mspace);
- H5Sclose(rspace);
- H5Dclose(dset);
- }
- }
- MPI_Bcast(
- this->state,
- this->array_size,
- MPI_DOUBLE,
- 0,
- this->comm);
- for (int i = 0; i<this->integration_steps; i++)
- MPI_Bcast(
- this->rhs[i],
- this->array_size,
- MPI_DOUBLE,
- 0,
- this->comm);
-}
-
-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 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)
- {
- this->write(data_file_id, "state", this->state);
- if (write_rhs)
- {
- 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;
- hid_t mspace = H5Screate_simple(4, count, NULL);
- for (int i=0; i<this->integration_steps; i++)
- {
- offset[1] = i;
- H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
- H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, this->rhs[i]);
- }
- H5Sclose(mspace);
- H5Sclose(wspace);
- H5Dclose(dset);
- }
- }
-}
-
-
-/*****************************************************************************/
-template class rFFTW_particles<VELOCITY_TRACER, float, 1>;
-template class rFFTW_particles<VELOCITY_TRACER, float, 2>;
-template class rFFTW_particles<VELOCITY_TRACER, float, 3>;
-template class rFFTW_particles<VELOCITY_TRACER, float, 4>;
-template class rFFTW_particles<VELOCITY_TRACER, float, 5>;
-template class rFFTW_particles<VELOCITY_TRACER, float, 6>;
-template class rFFTW_particles<VELOCITY_TRACER, double, 1>;
-template class rFFTW_particles<VELOCITY_TRACER, double, 2>;
-template class rFFTW_particles<VELOCITY_TRACER, double, 3>;
-template class rFFTW_particles<VELOCITY_TRACER, double, 4>;
-template class rFFTW_particles<VELOCITY_TRACER, double, 5>;
-template class rFFTW_particles<VELOCITY_TRACER, double, 6>;
-/*****************************************************************************/
diff --git a/documentation/_static/cpp.rst b/documentation/_static/cpp.rst
new file mode 100644
index 0000000000000000000000000000000000000000..b54e600897920696dbf67fd880c5a270a4286707
--- /dev/null
+++ b/documentation/_static/cpp.rst
@@ -0,0 +1,64 @@
+===========
+C++ classes
+===========
+
+-------------
+Interpolators
+-------------
+
+The point of these classes is to take care of interpolation.
+Given the number of neighbours used, and the type of interpolation, they
+compute the polynomials, and they compute the big sum.
+They need to have at least the following public methods:
+
+.. code:: cpp
+
+ /* map real locations to grid coordinates */
+ void get_grid_coordinates(
+ const int nparticles,
+ const int pdimension,
+ const double *__restrict__ x,
+ int *__restrict__ xg,
+ double *__restrict__ xx);
+ /* interpolate field at an array of locations */
+ void sample(
+ const int nparticles,
+ const int pdimension,
+ const double *__restrict__ x,
+ double *__restrict__ y,
+ const int *deriv = NULL);
+ /* interpolate 1 point */
+ void operator()(
+ const int *__restrict__ xg,
+ const double *__restrict__ xx,
+ double *__restrict__ dest,
+ const int *deriv = NULL);
+
+interpolator
+------------
+
+fields need to be padded, so that interpolation is a 1cpu job.
+
+rFFTW_interpolator
+------------------
+
+fields are not padded, computation is synchronized across all CPUs.
+
+-----------------
+Particle trackers
+-----------------
+
+The point of these classes is to solve ODEs.
+They need to be coupled to an interpolator.
+
+particles
+---------
+
+*work in progress* distributed solver.
+
+rFFTW_particles
+---------------
+
+*NOT* distributed.
+Particle data is syncrhonized across all CPUs.
+
diff --git a/documentation/index.rst b/documentation/index.rst
index 4f175c7f95191d414ba040daba075fb44fea38bd..9b0abc23648a08a9bd6692d2bc7651323c1b09d8 100644
--- a/documentation/index.rst
+++ b/documentation/index.rst
@@ -14,6 +14,7 @@ Welcome to bfps's documentation!
/_static/overview
/_static/development
/_static/api
+ /_static/cpp
Indices and tables
diff --git a/setup.py b/setup.py
index cb4ab8c3c2a456d073ba846803113b05e6e95264..8f982144e9188dbc74eea1bb4d6ef1f32508b45c 100644
--- a/setup.py
+++ b/setup.py
@@ -81,7 +81,8 @@ else:
if (('develop' in git_branch) or
('feature' in git_branch) or
('bugfix' in git_branch)):
- VERSION = subprocess.check_output(['git', 'describe', '--tags']).strip().decode()
+ VERSION = subprocess.check_output(
+ ['git', 'describe', '--tags', '--dirty']).strip().decode().replace('-g', '+g').replace('-dirty', '.dirty').replace('-', '.post')
else:
VERSION = subprocess.check_output(['git', 'describe', '--tags']).strip().decode().split('-')[0]
print('This is bfps version ' + VERSION)
@@ -90,10 +91,14 @@ print('This is bfps version ' + VERSION)
### lists of files and MANIFEST.in
src_file_list = ['field_descriptor',
+ 'interpolator_base',
+ 'distributed_particles',
+ 'particles_base',
+ 'interpolator',
+ 'particles',
+ 'rFFTW_interpolator',
'fluid_solver_base',
'fluid_solver',
- 'rFFTW_interpolator',
- 'rFFTW_particles',
'fftw_tools',
'spline_n1',
'spline_n2',
@@ -103,9 +108,7 @@ src_file_list = ['field_descriptor',
'spline_n6',
'Lagrange_polys']
-header_list = (['cpp/base.hpp',
- 'cpp/particles_base.hpp',
- 'cpp/interpolator_base.hpp'] +
+header_list = (['cpp/base.hpp'] +
['cpp/' + fname + '.hpp'
for fname in src_file_list])
diff --git a/tests/base.py b/tests/base.py
index 7b75080a44e12ffbd11243749fbbcd394671731d..2c616aeb1bbfbdc1b346ae365836429daafa93b0 100644
--- a/tests/base.py
+++ b/tests/base.py
@@ -115,7 +115,9 @@ def launch(
dt = None,
tracer_state_file = None,
vorticity_field = None,
- code_class = bfps.NavierStokes):
+ code_class = bfps.NavierStokes,
+ particle_class = 'particles',
+ interpolator_class = 'rFFTW_interpolator'):
c = code_class(
work_dir = opt.work_dir,
fluid_precision = opt.precision,
@@ -136,19 +138,23 @@ def launch(
c.parameters['famplitude'] = 0.2
c.fill_up_fluid_code()
if c.parameters['nparticles'] > 0:
+ c.name += '-' + particle_class
c.add_3D_rFFTW_field(name = 'rFFTW_acc')
c.add_interpolator(
name = 'spline',
neighbours = opt.neighbours,
- smoothness = opt.smoothness)
+ smoothness = opt.smoothness,
+ class_name = interpolator_class)
c.add_particles(
kcut = ['fs->kM/2', 'fs->kM/3'],
integration_steps = 3,
- interpolator = 'spline')
+ interpolator = 'spline',
+ class_name = particle_class)
c.add_particles(
integration_steps = [2, 3, 4, 6],
interpolator = 'spline',
- acc_name = 'rFFTW_acc')
+ acc_name = 'rFFTW_acc',
+ class_name = particle_class)
c.finalize_code()
c.write_src()
c.write_par()
@@ -189,10 +195,12 @@ def compare_stats(
key,
np.max(np.abs(c0.statistics[key + '(t)'] - c0.statistics[key + '(t)']))))
for i in range(c0.particle_species):
- print('maximum traj difference species {0} is {1}'.format(
+ print('maximum traj difference species {0} is {1} {2}'.format(
i,
np.max(np.abs(c0.get_particle_file()['tracers{0}/state'.format(i)][:] -
- c1.get_particle_file()['tracers{0}/state'.format(i)][:]))))
+ c1.get_particle_file()['tracers{0}/state'.format(i)][:])),
+ np.nanmax(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_interpolation.py b/tests/test_interpolation.py
index 29755bd6a5207c322c7e442df89be28b7a996543..924b3508ae3f82a275c14ae3df27a417c95b5bb8 100644
--- a/tests/test_interpolation.py
+++ b/tests/test_interpolation.py
@@ -49,20 +49,24 @@ if __name__ == '__main__':
frozen_fields = True)
c.pars_from_namespace(opt)
c.fill_up_fluid_code()
- c.add_particle_fields(
+ c.add_interpolator(
name = 'n{0}m{1}'.format(opt.neighbours, opt.smoothness),
neighbours = opt.neighbours,
- smoothness = opt.smoothness)
- c.add_particle_fields(
+ smoothness = opt.smoothness,
+ class_name = 'interpolator')
+ c.add_interpolator(
name = 'n{0}'.format(opt.neighbours),
neighbours = opt.neighbours,
- interp_type = 'Lagrange')
+ interp_type = 'Lagrange',
+ class_name = 'interpolator')
c.add_particles(
integration_steps = 1,
- fields_name = 'n{0}m{1}'.format(opt.neighbours, opt.smoothness))
+ interpolator = 'n{0}m{1}'.format(opt.neighbours, opt.smoothness),
+ class_name = 'particles')
c.add_particles(
integration_steps = 1,
- fields_name = 'n{0}'.format(opt.neighbours))
+ interpolator = 'n{0}'.format(opt.neighbours),
+ class_name = 'particles')
c.finalize_code()
c.write_src()
c.write_par()
@@ -84,17 +88,17 @@ if __name__ == '__main__':
c.set_host_info({'type' : 'pc'})
if opt.run:
c.run(ncpu = opt.ncpu)
- df = c.get_data_file()
+ df = c.get_particle_file()
fig = plt.figure(figsize=(10,5))
a = fig.add_subplot(121)
a.contour(
- df['particles/tracers0/state'][0, :, 0].reshape(nn, nn),
- df['particles/tracers0/state'][0, :, 2].reshape(nn, nn),
- df['particles/tracers0/velocity'][1, :, 0].reshape(nn, nn))
+ df['tracers0/state'][0, :, 0].reshape(nn, nn),
+ df['tracers0/state'][0, :, 2].reshape(nn, nn),
+ df['tracers0/velocity'][1, :, 0].reshape(nn, nn))
a = fig.add_subplot(122)
a.contour(
- df['particles/tracers1/state'][0, :, 0].reshape(nn, nn),
- df['particles/tracers1/state'][0, :, 2].reshape(nn, nn),
- df['particles/tracers1/velocity'][1, :, 0].reshape(nn, nn))
+ df['tracers1/state'][0, :, 0].reshape(nn, nn),
+ df['tracers1/state'][0, :, 2].reshape(nn, nn),
+ df['tracers1/velocity'][1, :, 0].reshape(nn, nn))
fig.savefig('interp.pdf')
diff --git a/tests/test_interpolation2.py b/tests/test_interpolation2.py
new file mode 100644
index 0000000000000000000000000000000000000000..b9c366c1bbb4b063c51a14642eac3d7b048955a7
--- /dev/null
+++ b/tests/test_interpolation2.py
@@ -0,0 +1,104 @@
+#######################################################################
+# #
+# Copyright 2015 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 #
+# #
+#######################################################################
+
+
+
+import numpy as np
+from base import *
+import matplotlib.pyplot as plt
+
+if __name__ == '__main__':
+ opt = parser.parse_args(
+ ['-n', '32',
+ '--run',
+ '--ncpu', '2',
+ '--initialize',
+ '--neighbours', '4',
+ '--smoothness', '1',
+ '--nparticles', '4225',
+ '--niter_todo', '1',
+ '--niter_stat', '1',
+ '--niter_part', '1',
+ '--niter_out', '1',
+ '--wd', 'interp'] +
+ sys.argv[1:])
+ c = bfps.NavierStokes(
+ name = 'test_interpolation',
+ use_fftw_wisdom = False,
+ frozen_fields = True)
+ c.pars_from_namespace(opt)
+ c.fill_up_fluid_code()
+ c.add_interpolator(
+ name = 'buffered',
+ neighbours = opt.neighbours,
+ smoothness = opt.smoothness,
+ class_name = 'interpolator')
+ c.add_interpolator(
+ name = 'rFFTW',
+ neighbours = opt.neighbours,
+ smoothness = opt.smoothness,
+ class_name = 'rFFTW_interpolator')
+ c.add_particles(
+ integration_steps = 1,
+ interpolator = 'buffered',
+ class_name = 'particles')
+ c.add_particles(
+ integration_steps = 1,
+ interpolator = 'rFFTW',
+ class_name = 'particles')
+ c.finalize_code()
+ c.write_src()
+ c.write_par()
+ nn = int(c.parameters['nparticles']**.5)
+ pos = np.zeros((nn, nn, 3), dtype = np.float64)
+ pos[:, :, 0] = np.linspace(0, 2*np.pi, nn)[None, :]
+ pos[:, :, 2] = np.linspace(0, 2*np.pi, nn)[:, None]
+ pos[:, :, 1] = np.random.random()*2*np.pi
+ c.generate_vector_field(write_to_file = True,
+ spectra_slope = 1.)
+ c.generate_tracer_state(
+ species = 0,
+ write_to_file = False,
+ data = pos.reshape(-1, 3))
+ c.generate_tracer_state(
+ species = 1,
+ write_to_file = False,
+ data = pos.reshape(-1, 3))
+ c.set_host_info({'type' : 'pc'})
+ if opt.run:
+ c.run(ncpu = opt.ncpu)
+ df = c.get_particle_file()
+ fig = plt.figure(figsize=(10,5))
+ a = fig.add_subplot(121)
+ a.contour(
+ df['tracers0/state'][0, :, 0].reshape(nn, nn),
+ df['tracers0/state'][0, :, 2].reshape(nn, nn),
+ df['tracers0/velocity'][1, :, 0].reshape(nn, nn))
+ a = fig.add_subplot(122)
+ a.contour(
+ df['tracers1/state'][0, :, 0].reshape(nn, nn),
+ df['tracers1/state'][0, :, 2].reshape(nn, nn),
+ df['tracers1/velocity'][1, :, 0].reshape(nn, nn))
+ fig.savefig('interp.pdf')
+
diff --git a/tests/test_particle_classes.py b/tests/test_particle_classes.py
new file mode 100644
index 0000000000000000000000000000000000000000..6f38859a22c4f2e5a0ae71627899cd580af740d4
--- /dev/null
+++ b/tests/test_particle_classes.py
@@ -0,0 +1,68 @@
+#! /usr/bin/env python
+#######################################################################
+# #
+# Copyright 2015 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 #
+# #
+#######################################################################
+
+
+
+from base import *
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def scaling(opt):
+ wd = opt.work_dir
+ opt.work_dir = wd + '/N{0:0>3x}_1'.format(opt.n)
+ c0 = launch(opt, dt = 0.2/opt.n, particle_class = 'distributed_particles', interpolator_class = 'interpolator')
+ c0.compute_statistics()
+ print ('Re = {0:.0f}'.format(c0.statistics['Re']))
+ print ('Rlambda = {0:.0f}'.format(c0.statistics['Rlambda']))
+ print ('Lint = {0:.4e}, etaK = {1:.4e}'.format(c0.statistics['Lint'], c0.statistics['etaK']))
+ print ('Tint = {0:.4e}, tauK = {1:.4e}'.format(c0.statistics['Tint'], c0.statistics['tauK']))
+ print ('kMetaK = {0:.4e}'.format(c0.statistics['kMeta']))
+ for s in range(c0.particle_species):
+ acceleration_test(c0, species = s, m = 1)
+ opt.work_dir = wd + '/N{0:0>3x}_2'.format(opt.n)
+ c1 = launch(opt, dt = c0.parameters['dt'], particle_class = 'particles')
+ c1.compute_statistics()
+ compare_stats(opt, c0, c1)
+ opt.work_dir = wd + '/N{0:0>3x}_3'.format(opt.n)
+ c2 = launch(opt, dt = c0.parameters['dt'], particle_class = 'particles', interpolator_class = 'interpolator')
+ c2.compute_statistics()
+ compare_stats(opt, c0, c2)
+ compare_stats(opt, c1, c2)
+ return None
+
+if __name__ == '__main__':
+ opt = parser.parse_args(
+ ['-n', '32',
+ '--run',
+ '--initialize',
+ '--ncpu', '4',
+ '--nparticles', '10000',
+ '--niter_todo', '8',
+ '--precision', 'single',
+ '--wd', 'data/single'] +
+ sys.argv[1:])
+ scaling(opt)
+
diff --git a/tests/test_plain.py b/tests/test_plain.py
index ab745d2e3cd4027010e3aab527be6bd477baf4ad..e57d6a3b9c82aa197b5f272475a3501fb9474383 100644
--- a/tests/test_plain.py
+++ b/tests/test_plain.py
@@ -33,10 +33,24 @@ import matplotlib.pyplot as plt
parser.add_argument('--multiplejob',
dest = 'multiplejob', action = 'store_true')
+parser.add_argument(
+ '--particle-class',
+ default = 'particles',
+ dest = 'particle_class',
+ type = str)
+
+parser.add_argument(
+ '--interpolator-class',
+ default = 'interpolator',
+ dest = 'interpolator_class',
+ type = str)
+
def plain(opt):
wd = opt.work_dir
opt.work_dir = wd + '/N{0:0>3x}_1'.format(opt.n)
- c0 = launch(opt, dt = 0.2/opt.n)
+ c0 = launch(opt, dt = 0.2/opt.n,
+ particle_class = opt.particle_class,
+ interpolator_class = opt.interpolator_class)
c0.compute_statistics()
print ('Re = {0:.0f}'.format(c0.statistics['Re']))
print ('Rlambda = {0:.0f}'.format(c0.statistics['Rlambda']))
@@ -51,14 +65,19 @@ def plain(opt):
opt.work_dir = wd + '/N{0:0>3x}_2'.format(opt.n)
opt.njobs *= 2
opt.niter_todo = opt.niter_todo//2
- c1 = launch(opt, dt = c0.parameters['dt'])
+ c1 = launch(opt, dt = c0.parameters['dt'],
+ particle_class = opt.particle_class,
+ interpolator_class = opt.interpolator_class)
c1.compute_statistics()
opt.work_dir = wd + '/N{0:0>3x}_3'.format(opt.n)
opt.njobs = 3*opt.njobs//2
opt.niter_todo = 2*opt.niter_todo//3
- c2 = launch(opt, dt = c0.parameters['dt'])
+ c2 = launch(opt, dt = c0.parameters['dt'],
+ particle_class = opt.particle_class,
+ interpolator_class = opt.interpolator_class)
c2.compute_statistics()
compare_stats(opt, c0, c1)
+ compare_stats(opt, c0, c2)
return None
if __name__ == '__main__':
@@ -66,7 +85,7 @@ if __name__ == '__main__':
['-n', '32',
'--run',
'--initialize',
- '--ncpu', '2',
+ '--ncpu', '4',
'--nparticles', '1000',
'--niter_todo', '48',
'--precision', 'single',