Skip to content
Snippets Groups Projects
Commit 593f68fb authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

code compiles with "particles" class

parent 7404d39a
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
bfps/NavierStokes.py
+ 2
− 3
View file @ 593f68fb
......@@ -238,7 +238,6 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
integration_method = 'AdamsBashforth',
integration_steps = 2,
kcut = 'fs->kM',
force_vel_reset = True,
frozen_particles = False,
particle_type = 'old_tracers'):
self.parameters['integration_method{0}'.format(self.particle_species)] = integration_method
......@@ -342,8 +341,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.particle_start += ('ps{0} = new particles<VELOCITY_TRACER, {1}, true, 3, {2}>(\n' +
'fname, fs,\n' +
'nparticles,\n' +
'{2},\n' +
'niter_part, integration_steps{1});\n').format(self.particle_species, self.C_dtype, neighbours, beta_name)
'{3},\n' +
'niter_part, integration_steps{0});\n').format(self.particle_species, self.C_dtype, neighbours, beta_name)
self.particle_start += ('particle_field{0} = {1}_alloc_real(fs->rd->local_size);\n' +
'ps{0}->data = particle_field{0};\n').format(self.particle_species, FFTW)
self.particle_end += '{0}_free(particle_field{1});\n'.format(FFTW, self.particle_species)
......
bfps/cpp/particles.cpp
+ 90
− 75
View file @ 593f68fb
......@@ -687,82 +687,97 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
}
}
template <int particle_type, class rnumber, bool multistep, int ncomponents, int interp_neighbours>
void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours>::sample_vec_field(void *vec_field, double *vec_values)
{
rnumber *tmp_field = ((rnumber*)vec_field) + this->buffer_size;
double *vec_local = new double[3*this->nparticles];
std::fill_n(vec_local, 3*this->nparticles, 0.0);
int deriv[] = {0, 0, 0};
/* get grid coordinates */
int *xg = new int[3*this->nparticles];
double *xx = new double[3*this->nparticles];
this->get_grid_coordinates(this->state, xg, xx);
/* perform interpolation */
for (int p=0; p<this->nparticles; p++)
if (this->fs->rd->myrank == this->computing[p])
this->interpolation_formula(
tmp_field,
xg + p*3,
xx + p*3,
vec_local + p*3,
deriv);
MPI_Allreduce(
vec_local,
vec_values,
3*this->nparticles,
MPI_DOUBLE,
MPI_SUM,
this->fs->rd->comm);
delete[] xg;
delete[] xx;
delete[] vec_local;
}
/*****************************************************************************/
/* macro for specializations to numeric types compatible with FFTW */
#define PARTICLES_DEFINITIONS(R, MPI_RNUM) \
\
template <int particle_type, class rnumber, bool multistep, int ncomponents, int interp_neighbours> \
void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours>::rFFTW_to_buffered(R *src, R *dst) \
{ \
/* do big copy of middle stuff */ \
std::copy(src, \
src + this->fs->rd->local_size, \
dst + this->buffer_size); \
int rsrc; \
/* get upper slices */ \
for (int rdst = 0; rdst < this->fs->rd->nprocs; rdst++) \
{ \
rsrc = this->fs->rd->rank[(this->fs->rd->all_start0[rdst] + \
this->fs->rd->all_size0[rdst]) % \
this->fs->rd->sizes[0]]; \
if (this->fs->rd->myrank == rsrc) \
MPI_Send( \
(void*)(src), \
this->buffer_size, \
MPI_RNUM, \
rdst, \
2*(rsrc*this->fs->rd->nprocs + rdst), \
this->fs->rd->comm); \
if (this->fs->rd->myrank == rdst) \
MPI_Recv( \
(void*)(dst + this->buffer_size + this->fs->rd->local_size), \
this->buffer_size, \
MPI_RNUM, \
rsrc, \
2*(rsrc*this->fs->rd->nprocs + rdst), \
this->fs->rd->comm, \
MPI_STATUS_IGNORE); \
} \
/* get lower slices */ \
for (int rdst = 0; rdst < this->fs->rd->nprocs; rdst++) \
{ \
rsrc = this->fs->rd->rank[MOD(this->fs->rd->all_start0[rdst] - 1, \
this->fs->rd->sizes[0])]; \
if (this->fs->rd->myrank == rsrc) \
MPI_Send( \
(void*)(src + this->fs->rd->local_size - this->buffer_size), \
this->buffer_size, \
MPI_RNUM, \
rdst, \
2*(rsrc*this->fs->rd->nprocs + rdst)+1, \
this->fs->rd->comm); \
if (this->fs->rd->myrank == rdst) \
MPI_Recv( \
(void*)(dst), \
this->buffer_size, \
MPI_RNUM, \
rsrc, \
2*(rsrc*this->fs->rd->nprocs + rdst)+1, \
this->fs->rd->comm, \
MPI_STATUS_IGNORE); \
} \
} \
/*****************************************************************************/
/*****************************************************************************/
/* now actually use the macro defined above */
PARTICLES_DEFINITIONS(
float,
MPI_FLOAT)
PARTICLES_DEFINITIONS(
double,
MPI_DOUBLE)
/*****************************************************************************/
template <int particle_type, class rnumber, bool multistep, int ncomponents, int interp_neighbours>
void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours>::rFFTW_to_buffered(void *void_src, void *void_dst)
{
rnumber *src = (rnumber*)void_src;
rnumber *dst = (rnumber*)void_dst;
/* do big copy of middle stuff */
std::copy(src,
src + this->fs->rd->local_size,
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->fs->rd->nprocs; rdst++)
{
rsrc = this->fs->rd->rank[(this->fs->rd->all_start0[rdst] +
this->fs->rd->all_size0[rdst]) %
this->fs->rd->sizes[0]];
if (this->fs->rd->myrank == rsrc)
MPI_Send(
src,
this->buffer_size,
MPI_RNUM,
rdst,
2*(rsrc*this->fs->rd->nprocs + rdst),
this->fs->rd->comm);
if (this->fs->rd->myrank == rdst)
MPI_Recv(
dst + this->buffer_size + this->fs->rd->local_size,
this->buffer_size,
MPI_RNUM,
rsrc,
2*(rsrc*this->fs->rd->nprocs + rdst),
this->fs->rd->comm,
MPI_STATUS_IGNORE);
}
/* get lower slices */
for (int rdst = 0; rdst < this->fs->rd->nprocs; rdst++)
{
rsrc = this->fs->rd->rank[MOD(this->fs->rd->all_start0[rdst] - 1,
this->fs->rd->sizes[0])];
if (this->fs->rd->myrank == rsrc)
MPI_Send(
src + this->fs->rd->local_size - this->buffer_size,
this->buffer_size,
MPI_RNUM,
rdst,
2*(rsrc*this->fs->rd->nprocs + rdst)+1,
this->fs->rd->comm);
if (this->fs->rd->myrank == rdst)
MPI_Recv(
dst,
this->buffer_size,
MPI_RNUM,
rsrc,
2*(rsrc*this->fs->rd->nprocs + rdst)+1,
this->fs->rd->comm,
MPI_STATUS_IGNORE);
}
}
/*****************************************************************************/
......
bfps/cpp/particles.hpp
+ 2
− 3
View file @ 593f68fb
......@@ -114,8 +114,7 @@ class particles
const int TRAJ_SKIP,
const int INTEGRATION_STEPS = 2);
~particles();
void rFFTW_to_buffered(float *src, float *dst);
void rFFTW_to_buffered(double *src, double *dst);
void rFFTW_to_buffered(void *src, void *dst);
/* an Euler step is needed to compute an estimate of future positions,
* which is needed for synchronization.
......@@ -128,7 +127,7 @@ class particles
void synchronize_single_particle_state(int p, double *x, int source_id = -1);
void get_grid_coordinates(double *x, int *xg, double *xx);
void interpolation_formula(rnumber *field, int *xg, double *xx, double *dest, int *deriv);
void sample_vec_field(void *vec_field, double *vec_values);
/* input/output */
void read(hid_t data_file_id);
......
tests/test_base.py
+ 27
− 21
View file @ 593f68fb
......@@ -92,28 +92,34 @@ def launch(
c.parameters['niter_part'] = 1
c.parameters['famplitude'] = 0.2
if c.parameters['nparticles'] > 0:
c.add_particles(kcut = 'fs->kM/2',
integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness)
c.add_particles(
integration_steps = 1,
neighbours = opt.neighbours,
smoothness = opt.smoothness,
force_vel_reset = True)
c.add_particles(integration_steps = 2, neighbours = opt.neighbours, smoothness = opt.smoothness)
c.add_particles(integration_steps = 3, neighbours = opt.neighbours, smoothness = opt.smoothness)
c.add_particles(integration_steps = 4, neighbours = opt.neighbours, smoothness = opt.smoothness)
c.add_particles(integration_steps = 5, neighbours = opt.neighbours, smoothness = opt.smoothness)
c.add_particles(integration_steps = 6, neighbours = opt.neighbours, smoothness = opt.smoothness)
c.add_particles(integration_steps = 2, neighbours = 1, smoothness = 0)
c.add_particles(integration_steps = 2, neighbours = 1, smoothness = 1)
c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 1)
c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 2)
c.add_particles(integration_steps = 2, neighbours = 1, interp_type = 'Lagrange')
c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange')
c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 1, integration_method = 'Heun')
c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange', integration_method = 'Heun')
c.add_particles(integration_steps = 4, neighbours = 6, smoothness = 1, integration_method = 'cRK4')
c.add_particles(integration_steps = 4, neighbours = 6, interp_type = 'Lagrange', integration_method = 'cRK4')
kcut = 'fs->kM/2',
integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness,
particle_type = 'new_tracers')
for integr_steps in range(1, 7):
c.add_particles(
integration_steps = 1,
neighbours = opt.neighbours,
smoothness = opt.smoothness,
particle_type = 'new_tracers')
for info in [(2, 1, 0), (2, 1, 1), (2, 6, 1), (2, 6, 2)]:
c.add_particles(
integration_steps = info[0],
neighbours = info[1],
smoothness = info[2],
particle_type = 'new_tracers')
c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 1, integration_method = 'Heun',
particle_type = 'new_tracers')
c.add_particles(integration_steps = 4, neighbours = 6, smoothness = 1, integration_method = 'cRK4',
particle_type = 'new_tracers')
c.add_particles(integration_steps = 2, neighbours = 1, interp_type = 'Lagrange',
particle_type = 'new_tracers')
c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange',
particle_type = 'new_tracers')
c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange', integration_method = 'Heun',
particle_type = 'new_tracers')
c.add_particles(integration_steps = 4, neighbours = 6, interp_type = 'Lagrange', integration_method = 'cRK4',
particle_type = 'new_tracers')
c.fill_up_fluid_code()
c.finalize_code()
c.write_src()
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment