Commit 287bb9af authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

rFFTW no longer has own field

parent 32301028
......@@ -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':
......
......@@ -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,
......
......@@ -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
......
......@@ -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,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment