diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index c77cb3098f204a39ab5cc47de9337b2e66feaea5..8b3a30ca6881a8051c039c1b39cb07b8d59579a5 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -230,6 +230,33 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.fluid_output + '\n}\n' +
'delete fs;\n')
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(
self,
neighbours = 1,
@@ -239,7 +266,12 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
integration_steps = 2,
kcut = 'fs->kM',
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['interp_type{0}'.format(self.particle_species)] = interp_type
self.parameters['neighbours{0}'.format(self.particle_species)] = neighbours
@@ -255,21 +287,25 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
H5Gclose(group);
//endcpp
""".format(self.particle_species)
- update_field = 'fs->compute_velocity(fs->cvorticity);\n'
- if kcut != 'fs->kM':
- update_field += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
- update_field += ('fs->ift_velocity();\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)
+ if particle_type == 'old_tracers':
+ update_field = 'fs->compute_velocity(fs->cvorticity);\n'
+ if kcut != 'fs->kM':
+ update_field += 'fs->low_pass_Fourier(fs->cvelocity, 3, {0});\n'.format(kcut)
+ update_field += ('fs->ift_velocity();\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)
+ compute_acc = ('{0} *acc_field_tmp = {1}_alloc_real(fs->rd->local_size);\n' +
+ 'fs->compute_Lagrangian_acceleration(acc_field_tmp);\n' +
+ 'ps{2}->rFFTW_to_buffered(acc_field_tmp, ps{2}->data);\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)
+ 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'
- compute_acc = ('{0} *acc_field_tmp = {1}_alloc_real(fs->rd->local_size);\n' +
- 'fs->compute_Lagrangian_acceleration(acc_field_tmp);\n' +
- 'ps{2}->rFFTW_to_buffered(acc_field_tmp, ps{2}->data);\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)
output_vel_acc = """
//begincpp
{{
@@ -342,16 +378,13 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
multistep,
neighbours,
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' +
'fname, fs,\n' +
'nparticles,\n' +
'{4},\n' +
'niter_part, integration_steps{0});\n').format(
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' +
- '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}->data = vel_{1}->f+vel_{1}->buffer_size;\n').format(self.particle_species, fields_name)
self.particle_start += ('ps{0}->dt = dt;\n' +
'ps{0}->iteration = iteration;\n' +
update_field +
diff --git a/bfps/cpp/interpolations.hpp b/bfps/cpp/interpolations.hpp
index 73ae78b9a48b1d0dd931c7d6e0583cdfc6fb3cfd..426ac20be6e3cc82752489ea072b5872dbb084e7 100644
--- a/bfps/cpp/interpolations.hpp
+++ b/bfps/cpp/interpolations.hpp
@@ -24,6 +24,8 @@
+#include "field_descriptor.hpp"
+#include "fluid_solver_base.hpp"
#include "spline_n1.hpp"
#include "spline_n2.hpp"
#include "spline_n3.hpp"
@@ -41,5 +43,24 @@ typedef void (*base_polynomial_values)(
double fraction,
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
diff --git a/bfps/cpp/particles.cpp b/bfps/cpp/particles.cpp
index 6e412300a221a39c3075b440bc9e61e57c158464..0ffbdc0bb36b39e9101aa02f077ffd93bbc83c28 100644
--- a/bfps/cpp/particles.cpp
+++ b/bfps/cpp/particles.cpp
@@ -147,7 +147,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
/* get grid coordinates */
int *xg = new int[this->array_size];
double *xx = new double[this->array_size];
- rnumber *vel = this->data + this->buffer_size;
this->get_grid_coordinates(x, xg, xx);
for (int p=0; p<this->nparticles; p++)
{
@@ -156,7 +155,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
int crank = this->get_rank(x[p*3 + 2]);
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(
"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],
@@ -192,7 +191,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
int deriv[] = {0, 0, 0};
int *xg = new int[this->array_size];
double *xx = new double[this->array_size];
- rnumber *vel = this->data + this->buffer_size;
double tmp[3];
/* get grid coordinates */
this->get_grid_coordinates(this->state, xg, xx);
@@ -200,7 +198,7 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
/* perform interpolation */
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]);
if (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
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;
+ rnumber *tmp_field = ((rnumber*)vec_field);
double *vec_local = new double[3*this->nparticles];
std::fill_n(vec_local, 3*this->nparticles, 0.0);
int deriv[] = {0, 0, 0};
@@ -722,66 +720,6 @@ void particles<particle_type, rnumber, multistep, ncomponents, interp_neighbours
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>;
diff --git a/setup.py b/setup.py
index 3910cec7a87642ebaf6573f12379981180b09a99..9b387cf9af9b8f1a345e395d1ca3eeb1c5e80d9d 100644
--- a/setup.py
+++ b/setup.py
@@ -46,6 +46,7 @@ except:
VERSION = date_name
src_file_list = ['field_descriptor',
+ 'interpolations',
'fftw_tools',
'fluid_solver_base',
'fluid_solver',
@@ -60,7 +61,7 @@ src_file_list = ['field_descriptor',
'spline_n6',
'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
#with open('MANIFEST.in', 'w') as manifest_in_file:
diff --git a/tests/test_base.py b/tests/test_base.py
index 4e7ae85050305e190a9dbc85cb5b580f83af3b13..bf85671d8189aadfc47b74bbd2a42bef25c66f52 100755
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -92,38 +92,42 @@ 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,
- particle_type = 'new_tracers')
#c.add_particles(
# kcut = 'fs->kM/2',
# integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness,
# 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):
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)]:
+ particle_type = 'new_tracers',
+ 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(
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')
+ interp_type = info[3],
+ integration_method = info[4],
+ particle_type = 'new_tracers',
+ fields_name = 'regular')
c.fill_up_fluid_code()
c.finalize_code()
c.write_src()