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

try to use a global particle field

parent f2efedd9
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
bfps/NavierStokes.py
+ 49
− 16
View file @ 4ba8bc24
...@@ -230,6 +230,33 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): ...@@ -230,6 +230,33 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.fluid_output + '\n}\n' + self.fluid_output + '\n}\n' +
'delete fs;\n') 'delete fs;\n')
return None return None
def add_particle_fields(
self,
kcut = None,
neighbours = 1,
name = 'particle_field'):
self.particle_variables += ('buffered_vec_field<{0}> *vel_{1}, *acc_{1};\n' +
'{0} *{1}_tmp;\n').format(self.C_dtype, name)
self.particle_start += ('vel_{0} = new buffered_vec_field<{1}>(fs, {2});\n' +
'acc_{0} = new buffered_vec_field<{1}>(fs, {2});\n' +
'{0}_tmp = new {1}[acc_{0}->src_descriptor->local_size];\n').format(name, self.C_dtype, neighbours+1)
self.particle_end += ('delete vel_{0};\n' +
'delete acc_{0};\n' +
'delete[] {0}_tmp;\n').format(name)
update_fields = 'fs->compute_velocity(fs->cvorticity);\n'
if not type(kcut) == type(None):
update_fields += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
update_fields += ('fs->ift_velocity();\n' +
'clip_zero_padding(fs->rd, fs->rvelocity, 3);\n' +
'vel_{0}->read_rFFTW(fs->rvelocity);\n' +
'fs->compute_velocity(fs->cvorticity);\n' +
'fs->ift_velocity();\n' +
'fs->compute_Lagrangian_acceleration({0}_tmp);\n' +
'clip_zero_padding(fs->rd, {0}_tmp, 3);\n' +
'acc_{0}->read_rFFTW({0}_tmp);\n').format(name)
self.particle_start += update_fields
self.particle_loop += update_fields
return None
def add_particles( def add_particles(
self, self,
neighbours = 1, neighbours = 1,
...@@ -239,7 +266,12 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): ...@@ -239,7 +266,12 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
integration_steps = 2, integration_steps = 2,
kcut = 'fs->kM', kcut = 'fs->kM',
frozen_particles = False, frozen_particles = False,
particle_type = 'old_tracers'): particle_type = 'old_tracers',
fields_name = None):
if integration_method == 'cRK4':
integration_steps = 4
elif integration_method == 'Heun':
integration_steps = 2
self.parameters['integration_method{0}'.format(self.particle_species)] = integration_method self.parameters['integration_method{0}'.format(self.particle_species)] = integration_method
self.parameters['interp_type{0}'.format(self.particle_species)] = interp_type self.parameters['interp_type{0}'.format(self.particle_species)] = interp_type
self.parameters['neighbours{0}'.format(self.particle_species)] = neighbours self.parameters['neighbours{0}'.format(self.particle_species)] = neighbours
...@@ -255,21 +287,25 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): ...@@ -255,21 +287,25 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
H5Gclose(group); H5Gclose(group);
//endcpp //endcpp
""".format(self.particle_species) """.format(self.particle_species)
if particle_type == 'old_tracers':
update_field = 'fs->compute_velocity(fs->cvorticity);\n' update_field = 'fs->compute_velocity(fs->cvorticity);\n'
if kcut != 'fs->kM': if kcut != 'fs->kM':
update_field += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut) update_field += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
update_field += ('fs->ift_velocity();\n' + update_field += ('fs->ift_velocity();\n' +
'clip_zero_padding(fs->rd, fs->rvelocity, 3);\n' + 'clip_zero_padding(fs->rd, fs->rvelocity, 3);\n' +
'ps{0}->rFFTW_to_buffered(fs->rvelocity, ps{0}->data);\n').format(self.particle_species) 'ps{0}->rFFTW_to_buffered(fs->rvelocity, ps{0}->data);\n').format(self.particle_species)
if self.dtype == np.float32:
FFTW = 'fftwf'
elif self.dtype == np.float64:
FFTW = 'fftw'
compute_acc = ('{0} *acc_field_tmp = {1}_alloc_real(fs->rd->local_size);\n' + compute_acc = ('{0} *acc_field_tmp = {1}_alloc_real(fs->rd->local_size);\n' +
'fs->compute_Lagrangian_acceleration(acc_field_tmp);\n' + 'fs->compute_Lagrangian_acceleration(acc_field_tmp);\n' +
'ps{2}->rFFTW_to_buffered(acc_field_tmp, ps{2}->data);\n' + 'ps{2}->rFFTW_to_buffered(acc_field_tmp, ps{2}->data);\n' +
'ps{2}->sample_vec_field(ps{2}->data, acceleration);\n' + 'ps{2}->sample_vec_field(ps{2}->data, acceleration);\n' +
'{1}_free(acc_field_tmp);\n').format(self.C_dtype, FFTW, self.particle_species) '{1}_free(acc_field_tmp);\n').format(self.C_dtype, FFTW, self.particle_species)
else:
update_field = ''
compute_acc = 'ps{0}->sample_vec_field(acc_{1}->f+acc_{1}->buffer_size, acceleration);\n'.format(self.particle_species, fields_name)
if self.dtype == np.float32:
FFTW = 'fftwf'
elif self.dtype == np.float64:
FFTW = 'fftw'
output_vel_acc = """ output_vel_acc = """
//begincpp //begincpp
{{ {{
...@@ -342,16 +378,13 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base): ...@@ -342,16 +378,13 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
multistep, multistep,
neighbours, neighbours,
self.particle_species) self.particle_species)
self.particle_variables += '{0} *particle_field{1};\n'.format(self.C_dtype, self.particle_species)
self.particle_start += ('ps{0} = new particles<VELOCITY_TRACER, {1}, {2}, 3, {3}>(\n' + self.particle_start += ('ps{0} = new particles<VELOCITY_TRACER, {1}, {2}, 3, {3}>(\n' +
'fname, fs,\n' + 'fname, fs,\n' +
'nparticles,\n' + 'nparticles,\n' +
'{4},\n' + '{4},\n' +
'niter_part, integration_steps{0});\n').format( 'niter_part, integration_steps{0});\n').format(
self.particle_species, self.C_dtype, multistep, neighbours, beta_name) self.particle_species, self.C_dtype, multistep, neighbours, beta_name)
self.particle_start += ('particle_field{0} = {1}_alloc_real(ps{0}->buffered_field_descriptor->local_size);\n' + self.particle_start += ('ps{0}->data = vel_{1}->f+vel_{1}->buffer_size;\n').format(self.particle_species, fields_name)
'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)
self.particle_start += ('ps{0}->dt = dt;\n' + self.particle_start += ('ps{0}->dt = dt;\n' +
'ps{0}->iteration = iteration;\n' + 'ps{0}->iteration = iteration;\n' +
update_field + update_field +
......
bfps/cpp/interpolations.hpp
+ 21
− 0
View file @ 4ba8bc24
...@@ -24,6 +24,8 @@ ...@@ -24,6 +24,8 @@
#include "field_descriptor.hpp"
#include "fluid_solver_base.hpp"
#include "spline_n1.hpp" #include "spline_n1.hpp"
#include "spline_n2.hpp" #include "spline_n2.hpp"
#include "spline_n3.hpp" #include "spline_n3.hpp"
...@@ -41,5 +43,24 @@ typedef void (*base_polynomial_values)( ...@@ -41,5 +43,24 @@ typedef void (*base_polynomial_values)(
double fraction, double fraction,
double *destination); double *destination);
template <class rnumber>
class buffered_vec_field
{
public:
int buffer_width;
ptrdiff_t buffer_size;
field_descriptor<rnumber> *descriptor;
field_descriptor<rnumber> *src_descriptor;
rnumber *f;
buffered_vec_field(
fluid_solver_base<rnumber> *FSOLVER,
const int BUFFER_WIDTH);
~buffered_vec_field();
/* destroys input */
int read_rFFTW(void *src);
};
#endif//INTERPOLATIONS #endif//INTERPOLATIONS
bfps/cpp/particles.cpp
+ 3
− 65
View file @ 4ba8bc24
...@@ -147,7 +147,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours ...@@ -147,7 +147,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
/* get grid coordinates */ /* get grid coordinates */
int *xg = new int[this->array_size]; int *xg = new int[this->array_size];
double *xx = new double[this->array_size]; double *xx = new double[this->array_size];
rnumber *vel = this->data + this->buffer_size;
this->get_grid_coordinates(x, xg, xx); this->get_grid_coordinates(x, xg, xx);
for (int p=0; p<this->nparticles; p++) for (int p=0; p<this->nparticles; p++)
{ {
...@@ -156,7 +155,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours ...@@ -156,7 +155,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
int crank = this->get_rank(x[p*3 + 2]); int crank = this->get_rank(x[p*3 + 2]);
if (this->fs->rd->myrank == crank) if (this->fs->rd->myrank == crank)
{ {
this->interpolation_formula(vel, xg + p*3, xx + p*3, y + p*3, deriv); this->interpolation_formula(this->data, xg + p*3, xx + p*3, y + p*3, deriv);
DEBUG_MSG( DEBUG_MSG(
"position is %g %g %g %d %d %d %g %g %g, result is %g %g %g\n", "position is %g %g %g %d %d %d %g %g %g, result is %g %g %g\n",
x[p*3], x[p*3+1], x[p*3+2], x[p*3], x[p*3+1], x[p*3+2],
...@@ -192,7 +191,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours ...@@ -192,7 +191,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
int deriv[] = {0, 0, 0}; int deriv[] = {0, 0, 0};
int *xg = new int[this->array_size]; int *xg = new int[this->array_size];
double *xx = new double[this->array_size]; double *xx = new double[this->array_size];
rnumber *vel = this->data + this->buffer_size;
double tmp[3]; double tmp[3];
/* get grid coordinates */ /* get grid coordinates */
this->get_grid_coordinates(this->state, xg, xx); this->get_grid_coordinates(this->state, xg, xx);
...@@ -200,7 +198,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours ...@@ -200,7 +198,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
/* perform interpolation */ /* perform interpolation */
for (int p=0; p<this->nparticles; p++) if (this->fs->rd->myrank == this->computing[p]) for (int p=0; p<this->nparticles; p++) if (this->fs->rd->myrank == this->computing[p])
{ {
this->interpolation_formula(vel, xg + p*3, xx + p*3, tmp, deriv); this->interpolation_formula(this->data, xg + p*3, xx + p*3, tmp, deriv);
jump[p] = fabs(3*this->dt * tmp[2]); jump[p] = fabs(3*this->dt * tmp[2]);
if (jump[p] < this->dz*1.01) if (jump[p] < this->dz*1.01)
jump[p] = this->dz*1.01; jump[p] = this->dz*1.01;
...@@ -693,7 +691,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours ...@@ -693,7 +691,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
template <int particle_type, class rnumber, bool multistep, int ncomponents, int 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) 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; rnumber *tmp_field = ((rnumber*)vec_field);
double *vec_local = new double[3*this->nparticles]; double *vec_local = new double[3*this->nparticles];
std::fill_n(vec_local, 3*this->nparticles, 0.0); std::fill_n(vec_local, 3*this->nparticles, 0.0);
int deriv[] = {0, 0, 0}; int deriv[] = {0, 0, 0};
...@@ -722,66 +720,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours ...@@ -722,66 +720,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
delete[] vec_local; delete[] vec_local;
} }
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);
}
}
/*****************************************************************************/ /*****************************************************************************/
template class particles<VELOCITY_TRACER, float, true, 3, 1>; template class particles<VELOCITY_TRACER, float, true, 3, 1>;
......
setup.py
+ 2
− 1
View file @ 4ba8bc24
...@@ -46,6 +46,7 @@ except: ...@@ -46,6 +46,7 @@ except:
VERSION = date_name VERSION = date_name
src_file_list = ['field_descriptor', src_file_list = ['field_descriptor',
'interpolations',
'fftw_tools', 'fftw_tools',
'fluid_solver_base', 'fluid_solver_base',
'fluid_solver', 'fluid_solver',
...@@ -60,7 +61,7 @@ src_file_list = ['field_descriptor', ...@@ -60,7 +61,7 @@ src_file_list = ['field_descriptor',
'spline_n6', 'spline_n6',
'Lagrange_polys'] 'Lagrange_polys']
header_list = ['cpp/base.hpp', 'cpp/interpolations.hpp'] + ['cpp/' + fname + '.hpp' for fname in src_file_list] header_list = ['cpp/base.hpp'] + ['cpp/' + fname + '.hpp' for fname in src_file_list]
# not sure we need the MANIFEST.in file, but I might as well # not sure we need the MANIFEST.in file, but I might as well
#with open('MANIFEST.in', 'w') as manifest_in_file: #with open('MANIFEST.in', 'w') as manifest_in_file:
......
tests/test_base.py
+ 23
− 19
View file @ 4ba8bc24
...@@ -92,38 +92,42 @@ def launch( ...@@ -92,38 +92,42 @@ def launch(
c.parameters['niter_part'] = 1 c.parameters['niter_part'] = 1
c.parameters['famplitude'] = 0.2 c.parameters['famplitude'] = 0.2
if c.parameters['nparticles'] > 0: if c.parameters['nparticles'] > 0:
c.add_particles(
kcut = 'fs->kM/2',
integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness,
particle_type = 'new_tracers')
#c.add_particles( #c.add_particles(
# kcut = 'fs->kM/2', # kcut = 'fs->kM/2',
# integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness, # integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness,
# particle_type = 'old_tracers') # particle_type = 'old_tracers')
c.add_particle_fields(name = 'regular', neighbours = 6)
c.add_particle_fields(kcut = 'fs->kM/2', name = 'filtered', neighbours = opt.neighbours)
c.add_particles(
kcut = 'fs->kM/2',
integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness,
particle_type = 'new_tracers',
fields_name = 'filtered')
for integr_steps in range(1, 7): for integr_steps in range(1, 7):
c.add_particles( c.add_particles(
integration_steps = 1, integration_steps = 1,
neighbours = opt.neighbours, neighbours = opt.neighbours,
smoothness = opt.smoothness, smoothness = opt.smoothness,
particle_type = 'new_tracers') particle_type = 'new_tracers',
for info in [(2, 1, 0), (2, 1, 1), (2, 6, 1), (2, 6, 2)]: fields_name = 'regular')
for info in [(2, 1, 0, 'spline', 'AdamsBashforth'),
(2, 1, 1, 'spline', 'AdamsBashforth'),
(2, 6, 1, 'spline', 'AdamsBashforth'),
(2, 6, 2, 'spline', 'AdamsBashforth'),
(2, 6, 1, 'spline', 'Heun'),
(4, 6, 1, 'spline', 'cRK4'),
(2, 1, 1, 'Lagrange', 'AdamsBashforth'),
(2, 6, 1, 'Lagrange', 'AdamsBashforth'),
(2, 6, 1, 'Lagrange', 'Heun'),
(4, 6, 1, 'Lagrange', 'cRK4')]:
c.add_particles( c.add_particles(
integration_steps = info[0], integration_steps = info[0],
neighbours = info[1], neighbours = info[1],
smoothness = info[2], smoothness = info[2],
particle_type = 'new_tracers') interp_type = info[3],
c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 1, integration_method = 'Heun', integration_method = info[4],
particle_type = 'new_tracers') particle_type = 'new_tracers',
c.add_particles(integration_steps = 4, neighbours = 6, smoothness = 1, integration_method = 'cRK4', fields_name = 'regular')
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.fill_up_fluid_code()
c.finalize_code() c.finalize_code()
c.write_src() c.write_src()
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment