From 63e10688e4ec39df97b5b64412cfde7c0af1a50e Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Fri, 6 Nov 2015 14:17:09 +0100
Subject: [PATCH] introduce linear interpolation in time
I think. still checking...
---
bfps/NavierStokes.py | 4 ++--
bfps/cpp/interpolator.cpp | 34 +++++++++++++++++++++-------------
bfps/cpp/interpolator.hpp | 4 ++--
bfps/cpp/particles.cpp | 29 +++++++++++++++++------------
bfps/cpp/particles.hpp | 4 +++-
5 files changed, 45 insertions(+), 30 deletions(-)
diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index bf49357b..6f345644 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -256,8 +256,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
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}->f+acc_{0}->buffer_size);\n' +
- 'acc_{0}->read_rFFTW(acc_{0}->f+acc_{0}->buffer_size);\n').format(name)
+ 'fs->compute_Lagrangian_acceleration(acc_{0}->temp);\n' +
+ 'acc_{0}->read_rFFTW(acc_{0}->temp);\n').format(name)
self.particle_start += update_fields
self.particle_loop += update_fields
return None
diff --git a/bfps/cpp/interpolator.cpp b/bfps/cpp/interpolator.cpp
index a0940038..254213e2 100644
--- a/bfps/cpp/interpolator.cpp
+++ b/bfps/cpp/interpolator.cpp
@@ -43,25 +43,30 @@ interpolator<rnumber, interp_neighbours>::interpolator(
4, tdims,
this->unbuffered_descriptor->mpi_dtype,
this->unbuffered_descriptor->comm);
- this->f = new rnumber[this->descriptor->local_size];
- //if (sizeof(rnumber) == 4)
- // this->f = fftwf_alloc_real(this->descriptor->local_size);
- //else if (sizeof(rnumber) == 8)
- // this->f = fftw_alloc_real(this->descriptor->local_size);
+ this->f0 = new rnumber[this->descriptor->local_size];
+ this->f1 = new rnumber[this->descriptor->local_size];
+ this->temp = this->f1 + this->buffer_size;
}
template <class rnumber, int interp_neighbours>
interpolator<rnumber, interp_neighbours>::~interpolator()
{
- delete[] this->f;
+ delete[] this->f0;
+ delete[] this->f1;
delete this->descriptor;
}
template <class rnumber, int interp_neighbours>
int interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
{
+ /* first, roll fields */
+ rnumber *tmp = this->f0;
+ this->f0 = this->f1;
+ this->f1 = tmp;
+ this->temp = this->f0 + this->buffer_size;
+ /* now do regular things */
rnumber *src = (rnumber*)void_src;
- rnumber *dst = this->f;
+ rnumber *dst = this->f1;
/* do big copy of middle stuff */
clip_zero_padding<rnumber>(this->unbuffered_descriptor, src, 3);
std::copy(src,
@@ -121,6 +126,7 @@ int interpolator<rnumber, interp_neighbours>::read_rFFTW(void *void_src)
template <class rnumber, int interp_neighbours>
void interpolator<rnumber, interp_neighbours>::operator()(
+ double t,
int *xg,
double *xx,
double *dest,
@@ -140,16 +146,18 @@ void interpolator<rnumber, interp_neighbours>::operator()(
this->compute_beta(deriv[2], xx[2], bz);
}
std::fill_n(dest, 3, 0);
- rnumber *field = this->f + this->buffer_size;
+ rnumber *ff0 = this->f0 + this->buffer_size;
+ rnumber *ff1 = this->f1 + 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]);
+ {
+ ptrdiff_t tindex = ((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;
+ dest[c] += (ff0[tindex]*(1-t) + t*ff1[tindex])*(bz[iz+interp_neighbours]*by[iy+interp_neighbours]*bx[ix+interp_neighbours]);
+ }
}
template class interpolator<float, 1>;
diff --git a/bfps/cpp/interpolator.hpp b/bfps/cpp/interpolator.hpp
index 3e996f4a..4f2716a2 100644
--- a/bfps/cpp/interpolator.hpp
+++ b/bfps/cpp/interpolator.hpp
@@ -52,14 +52,14 @@ class interpolator
base_polynomial_values compute_beta;
field_descriptor<rnumber> *descriptor;
field_descriptor<rnumber> *unbuffered_descriptor;
- rnumber *f;
+ rnumber *f0, *f1, *temp;
interpolator(
fluid_solver_base<rnumber> *FSOLVER,
base_polynomial_values BETA_POLYS);
~interpolator();
- void operator()(int *xg, double *xx, double *dest, int *deriv = NULL);
+ void operator()(double t, int *xg, double *xx, double *dest, int *deriv = NULL);
/* destroys input */
int read_rFFTW(void *src);
};
diff --git a/bfps/cpp/particles.cpp b/bfps/cpp/particles.cpp
index 2afc73c9..7109f409 100644
--- a/bfps/cpp/particles.cpp
+++ b/bfps/cpp/particles.cpp
@@ -135,6 +135,7 @@ particles<particle_type, rnumber, multistep, interp_neighbours>::~particles()
template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
void particles<particle_type, rnumber, multistep, interp_neighbours>::sample_vec_field(
interpolator<rnumber, interp_neighbours> *vec,
+ double t,
double *x,
double *y,
const bool synch,
@@ -151,7 +152,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::sample_vec
for (int p=0; p<this->nparticles; p++)
{
if (this->myrank == this->computing[p])
- (*vec)(xg + p*3, xx + p*3, y + p*3, deriv);
+ (*vec)(t, xg + p*3, xx + p*3, y + p*3, deriv);
}
}
else
@@ -162,7 +163,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::sample_vec
{
int crank = this->get_rank(x[p*3 + 2]);
if (this->myrank == crank)
- (*vec)(xg + p*3, xx + p*3, y + p*3, deriv);
+ (*vec)(t, xg + p*3, xx + p*3, y + p*3, deriv);
if (crank != this->computing[p])
this->synchronize_single_particle_state(p, y, crank);
}
@@ -191,12 +192,12 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::sample_vec
}
template <int particle_type, class rnumber, bool multistep, int interp_neighbours>
-void particles<particle_type, rnumber, multistep, interp_neighbours>::get_rhs(double *x, double *y)
+void particles<particle_type, rnumber, multistep, interp_neighbours>::get_rhs(double t, double *x, double *y)
{
switch(particle_type)
{
case VELOCITY_TRACER:
- this->sample_vec_field(this->vel, x, y, false);
+ this->sample_vec_field(this->vel, t, x, y, false);
break;
}
}
@@ -209,7 +210,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::jump_estim
{
case VELOCITY_TRACER:
double *y = new double[3*this->nparticles];
- this->sample_vec_field(this->vel, this->state, y, false);
+ this->sample_vec_field(this->vel, 1.0, this->state, y, false);
for (int p=0; p<this->nparticles; p++) if (this->myrank == this->computing[p])
{
jump[p] = fabs(2*this->dt * y[3*p+2]);
@@ -358,7 +359,7 @@ template <int particle_type, class rnumber, bool multistep, int interp_neighbour
void particles<particle_type, rnumber, multistep, interp_neighbours>::AdamsBashforth(int nsteps)
{
ptrdiff_t ii;
- this->get_rhs(this->state, this->rhs[0]);
+ this->get_rhs(0, this->state, this->rhs[0]);
switch(nsteps)
{
case 1:
@@ -450,8 +451,11 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::Heun()
if (!multistep)
{
double *y = new double[this->array_size];
- double dtfactor[] = {0.0, this->dt};
- this->get_rhs(this->state, this->rhs[0]);
+ double dtfactor[2];
+ double factor[] = {0.0, 1.0};
+ dtfactor[0] = factor[0]*this->dt;
+ dtfactor[1] = factor[1]*this->dt;
+ this->get_rhs(0, this->state, this->rhs[0]);
for (int p=0; p<this->nparticles; p++)
this->synchronize_single_particle_state(p, this->rhs[0]);
for (int kindex = 1; kindex < 2; kindex++)
@@ -467,7 +471,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::Heun()
}
for (int p=0; p<this->nparticles; p++)
this->synchronize_single_particle_state(p, y);
- this->get_rhs(y, this->rhs[kindex]);
+ this->get_rhs(factor[kindex], y, this->rhs[kindex]);
for (int p=0; p<this->nparticles; p++)
this->synchronize_single_particle_state(p, this->rhs[kindex]);
}
@@ -492,9 +496,10 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::cRK4()
{
if (!multistep)
{
+ static double factor[] = {0.0, 0.5, 0.5, 1.0};
+ static double dtfactor[] = {0.0, this->dt/2, this->dt/2, this->dt};
double *y = new double[this->array_size];
- double dtfactor[] = {0.0, this->dt/2, this->dt/2, this->dt};
- this->get_rhs(this->state, this->rhs[0]);
+ this->get_rhs(0, this->state, this->rhs[0]);
for (int p=0; p<this->nparticles; p++)
this->synchronize_single_particle_state(p, this->rhs[0]);
for (int kindex = 1; kindex < 4; kindex++)
@@ -510,7 +515,7 @@ void particles<particle_type, rnumber, multistep, interp_neighbours>::cRK4()
}
for (int p=0; p<this->nparticles; p++)
this->synchronize_single_particle_state(p, y);
- this->get_rhs(y, this->rhs[kindex]);
+ this->get_rhs(factor[kindex], y, this->rhs[kindex]);
for (int p=0; p<this->nparticles; p++)
this->synchronize_single_particle_state(p, this->rhs[kindex]);
}
diff --git a/bfps/cpp/particles.hpp b/bfps/cpp/particles.hpp
index b42a4f16..49537de1 100644
--- a/bfps/cpp/particles.hpp
+++ b/bfps/cpp/particles.hpp
@@ -109,6 +109,7 @@ class particles
* */
void jump_estimate(double *jump_length);
void get_rhs(double *x, double *rhs);
+ void get_rhs(double t, double *x, double *rhs);
int get_rank(double z); // get rank for given value of z
void synchronize();
@@ -116,13 +117,14 @@ class particles
void get_grid_coordinates(double *x, int *xg, double *xx);
void sample_vec_field(
interpolator<rnumber, interp_neighbours> *vec,
+ double t,
double *x,
double *y,
const bool synch = false,
int *deriv = NULL);
inline void sample_vec_field(interpolator<rnumber, interp_neighbours> *field, double *vec_values)
{
- this->sample_vec_field(field, this->state, vec_values, true, NULL);
+ this->sample_vec_field(field, 1.0, this->state, vec_values, true, NULL);
}
/* input/output */
--
GitLab