diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 2f758ca2d0395e76634d01e5c0500a5374d14203..920b8413a95b5929255a5f06f4e1afb6a4c6cbf0 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -378,6 +378,52 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
'}\n' +
'delete fs;\n')
return None
+ def add_3D_rFFTW_field(
+ self,
+ name = 'rFFTW_acc'):
+ if self.dtype == np.float32:
+ FFTW = 'fftwf'
+ elif self.dtype == np.float64:
+ FFTW = 'fftw'
+ self.fluid_variables += '{0} *{1};\n'.format(self.C_dtype, name)
+ self.fluid_start += '{0} = {1}_alloc_real(2*fs->cd->local_size);\n'.format(name, FFTW)
+ self.fluid_end += '{0}_free({1});\n'.format(FFTW, name)
+ return None
+ def add_interpolator(
+ self,
+ interp_type = 'spline',
+ kcut = None,
+ neighbours = 1,
+ smoothness = 1,
+ name = 'field_interpolator',
+ field_name = 'fs->rvelocity'):
+ name = quantity + '_' + name
+ self.fluid_includes += '#include "rFFTW_interpolator.hpp"\n'
+ self.fluid_variables += 'rFFTW_interpolator <{0}, {1}> *{2};\n'.format(self.C_dtype, neighbours, name)
+ self.parameters[name + '_type'] = interp_type
+ self.parameters[name + '_neighbours'] = neighbours
+ if interp_type == 'spline':
+ self.parameters[name + '_smoothness'] = smoothness
+ 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});\n'.format(
+ name,
+ self.C_dtype,
+ neighbours,
+ beta_name)
+ self.fluid_end += 'delete {0};\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)
+ if quantity == 'vel':
+ update_fields += ('fs->ift_velocity();\n' +
+ '{0}->read_rFFTW(fs->rvelocity);\n').format(name)
+ elif quantity == 'acc':
+ update_fields += 'fs->compute_Lagrangian_acceleration({0}->field);\n'.format(name)
+ self.fluid_start += update_fields
+ self.fluid_loop += update_fields
+ return None
def add_particle_fields(
self,
interp_type = 'spline',
@@ -385,7 +431,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
neighbours = 1,
smoothness = 1,
name = 'particle_field',
- field_class = 'rFFTW_interpolator'):
+ field_class = 'rFFTW_interpolator',
+ acc_field_name = 'rFFTW_acc'):
self.fluid_includes += '#include "{0}.hpp"\n'.format(field_class)
self.fluid_variables += field_class + '<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(self.C_dtype, neighbours, name)
self.parameters[name + '_type'] = interp_type
@@ -395,19 +442,19 @@ class NavierStokes(bfps.fluid_base.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 += ('vel_{0} = new {1}<{2}, {3}>(fs, {4});\n' +
- 'acc_{0} = new {1}<{2}, {3}>(fs, {4});\n').format(name,
+ self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4}, fs->rvelocity);\n' +
+ 'acc_{0} = new {1}<{2}, {3}>(fs, {4}, {5});\n').format(name,
field_class,
self.C_dtype,
neighbours,
- beta_name)
+ beta_name,
+ acc_field_name)
self.fluid_end += ('delete vel_{0};\n' +
'delete acc_{0};\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' +
- 'vel_{0}->read_rFFTW(fs->rvelocity);\n' +
'fs->compute_Lagrangian_acceleration(acc_{0}->field);\n').format(name)
self.fluid_start += update_fields
self.fluid_loop += update_fields
@@ -419,7 +466,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
kcut = 'fs->kM',
frozen_particles = False,
fields_name = None,
- particle_class = 'rFFTW_particles'):
+ particle_class = 'rFFTW_particles',
+ acc_stats = False):
if integration_method == 'cRK4':
integration_steps = 4
elif integration_method == 'Heun':
diff --git a/bfps/cpp/rFFTW_interpolator.cpp b/bfps/cpp/rFFTW_interpolator.cpp
index a858eb104abf19837f7ccd30ae673caaf841a493..c701fbd64bba8e1c68ab2de953eae0795baead2f 100644
--- a/bfps/cpp/rFFTW_interpolator.cpp
+++ b/bfps/cpp/rFFTW_interpolator.cpp
@@ -32,15 +32,13 @@
template <class rnumber, int interp_neighbours>
rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
fluid_solver_base<rnumber> *fs,
- base_polynomial_values BETA_POLYS)
+ base_polynomial_values BETA_POLYS,
+ rnumber *FIELD)
{
this->descriptor = fs->rd;
this->field_size = 2*fs->cd->local_size;
this->compute_beta = BETA_POLYS;
- if (sizeof(rnumber) == 4)
- this->field = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
- else if (sizeof(rnumber) == 8)
- this->field = (rnumber*)((void*)fftw_alloc_real(this->field_size));
+ this->field = FIELD;
// compute dx, dy, dz;
this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
@@ -59,24 +57,9 @@ rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
template <class rnumber, int interp_neighbours>
rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
{
- if (sizeof(rnumber) == 4)
- fftwf_free((float*)((void*)this->field));
- else if (sizeof(rnumber) == 8)
- fftw_free((double*)((void*)this->field));
delete[] this->compute;
}
-template <class rnumber, int interp_neighbours>
-int rFFTW_interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
-{
- rnumber *src = (rnumber*)void_src;
- /* do big copy of middle stuff */
- std::copy(src,
- src + this->field_size,
- this->field);
- return EXIT_SUCCESS;
-}
-
template <class rnumber, int interp_neighbours>
void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
const int nparticles,
diff --git a/bfps/cpp/rFFTW_interpolator.hpp b/bfps/cpp/rFFTW_interpolator.hpp
index 154296018926202bb53367e4313510079343c4b8..3dccd7150d17bea612128bf6d70b37033d281db8 100644
--- a/bfps/cpp/rFFTW_interpolator.hpp
+++ b/bfps/cpp/rFFTW_interpolator.hpp
@@ -60,7 +60,8 @@ class rFFTW_interpolator
rFFTW_interpolator(
fluid_solver_base<rnumber> *FSOLVER,
- base_polynomial_values BETA_POLYS);
+ base_polynomial_values BETA_POLYS,
+ rnumber *FIELD_DATA);
~rFFTW_interpolator();
/* map real locations to grid coordinates */
@@ -83,7 +84,6 @@ class rFFTW_interpolator
const double *__restrict__ xx,
double *__restrict__ dest,
const int *deriv = NULL);
- int read_rFFTW(void *src);
};
#endif//RFFTW_INTERPOLATOR
diff --git a/tests/base.py b/tests/base.py
index c050245b85910bd35e2c62bf937945ba8d897e52..df1d65614ccedaa68491f961d8e7af7b178d8973 100644
--- a/tests/base.py
+++ b/tests/base.py
@@ -102,6 +102,7 @@ def launch(
c.parameters['famplitude'] = 0.2
c.fill_up_fluid_code()
if c.parameters['nparticles'] > 0:
+ c.add_3D_rFFTW_field()
c.add_particle_fields(
name = 'regular',
neighbours = opt.neighbours,