Skip to content

GitLab

Commit d4227289 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/distributed_particles' into develop

parents c0894829 26a7a8cc
Hide whitespace changes
Inline Side-by-side
bfps/NavierStokes.py
View file @ d4227289
......@@ -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),
......
bfps/_fluid_base.py
View file @ d4227289
......@@ -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' +
......
bfps/cpp/distributed_particles.cpp 0 → 100644
View file @ d4227289
/**********************************************************************
* *
* 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>