diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 29d9a8fc1bf4ca9cc2285d7e3851f8ace1ba583a..a25a6e78acd05ff81c8a708cd016e627705024e8 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -234,14 +234,23 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
return None
def add_particle_fields(
self,
+ interp_type = 'spline',
kcut = None,
neighbours = 1,
+ smoothness = 1,
name = 'particle_field'):
self.particle_variables += ('interpolator<{0}, {1}> *vel_{2}, *acc_{2};\n' +
'{0} *{2}_tmp;\n').format(self.C_dtype, neighbours, name)
- self.particle_start += ('vel_{0} = new interpolator<{1}, {2}>(fs);\n' +
- 'acc_{0} = new interpolator<{1}, {2}>(fs);\n' +
- '{0}_tmp = new {1}[acc_{0}->unbuffered_descriptor->local_size];\n').format(name, self.C_dtype, neighbours)
+ if interp_type == 'spline':
+ beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
+ elif interp_type == 'Lagrange':
+ beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
+ self.particle_start += ('vel_{0} = new interpolator<{1}, {2}>(fs, {3});\n' +
+ 'acc_{0} = new interpolator<{1}, {2}>(fs, {3});\n' +
+ '{0}_tmp = new {1}[acc_{0}->unbuffered_descriptor->local_size];\n').format(name,
+ self.C_dtype,
+ neighbours,
+ beta_name)
self.particle_end += ('delete vel_{0};\n' +
'delete acc_{0};\n' +
'delete[] {0}_tmp;\n').format(name)
diff --git a/bfps/cpp/interpolator.cpp b/bfps/cpp/interpolator.cpp
index 7609387b43fadeb34503379dfaae659322ef7ae7..ed8e16ac8ce981a5a04394edc629da832674f12f 100644
--- a/bfps/cpp/interpolator.cpp
+++ b/bfps/cpp/interpolator.cpp
@@ -28,11 +28,13 @@
template <class rnumber, int interp_neighbours>
interpolator<rnumber, interp_neighbours>::interpolator(
- fluid_solver_base<rnumber> *fs)
+ fluid_solver_base<rnumber> *fs,
+ base_polynomial_values 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];
@@ -116,6 +118,39 @@ int interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
return EXIT_SUCCESS;
}
+template <class rnumber, int interp_neighbours>
+void interpolator<rnumber, interp_neighbours>::operator()(
+ int *xg,
+ double *xx,
+ double *dest,
+ int *deriv)
+{
+ double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
+ if (deriv == NULL)
+ {
+ this->compute_beta(0, xx[0], bx);
+ this->compute_beta(0, xx[1], by);
+ this->compute_beta(0, xx[2], bz);
+ }
+ else
+ {
+ this->compute_beta(deriv[0], xx[0], bx);
+ this->compute_beta(deriv[1], xx[1], by);
+ this->compute_beta(deriv[2], xx[2], bz);
+ }
+ std::fill_n(dest, 3, 0);
+ rnumber *field = this->f + this->buffer_size;
+ for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
+ for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
+ for (int ix = -interp_neighbours; ix <= interp_neighbours+1; ix++)
+ for (int c=0; c<3; c++)
+ dest[c] += field[((ptrdiff_t( xg[2]+iz ) *this->descriptor->sizes[1] +
+ ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1])))*this->descriptor->sizes[2] +
+ ptrdiff_t(MOD(xg[0]+ix, this->descriptor->sizes[2])))*3+c]*(bz[iz+interp_neighbours]*
+ by[iy+interp_neighbours]*
+ bx[ix+interp_neighbours]);
+}
+
template class interpolator<float, 1>;
template class interpolator<float, 2>;
template class interpolator<float, 3>;
diff --git a/bfps/cpp/interpolator.hpp b/bfps/cpp/interpolator.hpp
index e6576453500b69354aa617b0c11a7ec4d2ba024e..8d1bc8c71c518e6cc2d066562967e0150c3e1a56 100644
--- a/bfps/cpp/interpolator.hpp
+++ b/bfps/cpp/interpolator.hpp
@@ -48,14 +48,17 @@ class interpolator
{
public:
ptrdiff_t buffer_size;
+ base_polynomial_values compute_beta;
field_descriptor<rnumber> *descriptor;
field_descriptor<rnumber> *unbuffered_descriptor;
rnumber *f;
interpolator(
- fluid_solver_base<rnumber> *FSOLVER);
+ fluid_solver_base<rnumber> *FSOLVER,
+ base_polynomial_values BETA_POLYS);
~interpolator();
+ void operator()(int *xg, double *xx, double *dest, int *deriv = NULL);
/* destroys input */
int read_rFFTW(void *src);
};
diff --git a/tests/test_base.py b/tests/test_base.py
index 25be39b49505cdb56ac0ece9371d8caf8e71a085..6ec5aa63a2109585940c941b0667d975acb8b180 100755
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -93,34 +93,28 @@ def launch(
c.parameters['niter_part'] = 1
c.parameters['famplitude'] = 0.2
if c.parameters['nparticles'] > 0:
- c.add_particle_fields(name = 'regular', neighbours = 6)
+ c.add_particle_fields(name = 'regular', neighbours = opt.neighbours, smoothness = opt.smoothness)
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,
fields_name = 'filtered')
- for integr_steps in range(1, 7):
- c.add_particles(
- integration_steps = 1,
- neighbours = opt.neighbours,
- smoothness = opt.smoothness,
- 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')]:
+ #for integr_steps in range(1, 7):
+ # c.add_particles(
+ # integration_steps = integr_steps,
+ # neighbours = opt.neighbours,
+ # smoothness = opt.smoothness,
+ # fields_name = 'regular')
+ for info in [(2, 1, 'spline', 'Heun'),
+ (4, 1, 'spline', 'cRK4'),
+ (2, 1, 'Lagrange', 'Heun'),
+ (4, 1, 'Lagrange', 'cRK4')]:
c.add_particles(
integration_steps = info[0],
- neighbours = info[1],
- smoothness = info[2],
- interp_type = info[3],
- integration_method = info[4],
+ neighbours = opt.neighbours,
+ smoothness = info[1],
+ interp_type = info[2],
+ integration_method = info[3],
fields_name = 'regular')
c.fill_up_fluid_code()
c.finalize_code()