From fda98c1ad3587e31990ebeb112e48d751ff84a38 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Wed, 4 Nov 2015 13:51:34 +0100
Subject: [PATCH] clean up particle code
---
bfps/cpp/slab_field_particles.cpp | 55 +------------------------------
bfps/cpp/slab_field_particles.hpp | 5 +--
bfps/cpp/tracers.cpp | 15 +++++----
tests/test_base.py | 42 +++++++++++------------
4 files changed, 31 insertions(+), 86 deletions(-)
diff --git a/bfps/cpp/slab_field_particles.cpp b/bfps/cpp/slab_field_particles.cpp
index c8e95a98..e2ff9e1b 100644
--- a/bfps/cpp/slab_field_particles.cpp
+++ b/bfps/cpp/slab_field_particles.cpp
@@ -536,10 +536,9 @@ void slab_field_particles<rnumber>::get_grid_coordinates(double *x, int *xg, dou
}
template <class rnumber>
-void slab_field_particles<rnumber>::spline_formula(rnumber *field, int *xg, double *xx, double *dest, int *deriv)
+void slab_field_particles<rnumber>::interpolation_formula(rnumber *field, int *xg, double *xx, double *dest, int *deriv)
{
double bx[this->interp_neighbours*2+2], by[this->interp_neighbours*2+2], bz[this->interp_neighbours*2+2];
- //DEBUG_MSG("entered spline_formula\n");
this->compute_beta(deriv[0], xx[0], bx);
this->compute_beta(deriv[1], xx[1], by);
this->compute_beta(deriv[2], xx[2], bz);
@@ -564,58 +563,6 @@ void slab_field_particles<rnumber>::spline_formula(rnumber *field, int *xg, doub
by[iy+this->interp_neighbours]*
bx[ix+this->interp_neighbours]);
}
- //DEBUG_MSG("finishing spline_formula\n");
-}
-
-template <class rnumber>
-void slab_field_particles<rnumber>::spline_n1_formula(rnumber *field, int *xg, double *xx, double *dest, int *deriv)
-{
- double bx[4], by[4], bz[4];
- 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);
- for (int iz = -1; iz <= 2; iz++)
- for (int iy = -1; iy <= 2; iy++)
- for (int ix = -1; ix <= 2; ix++)
- for (int c=0; c<3; c++)
- dest[c] += field[((ptrdiff_t( xg[2]+iz ) *this->fs->rd->subsizes[1] +
- ptrdiff_t(MOD(xg[1]+iy, this->fs->rd->sizes[1])))*this->fs->rd->subsizes[2] +
- ptrdiff_t(MOD(xg[0]+ix, this->fs->rd->sizes[2])))*3+c]*bz[iz+1]*by[iy+1]*bx[ix+1];
-}
-
-template <class rnumber>
-void slab_field_particles<rnumber>::spline_n2_formula(rnumber *field, int *xg, double *xx, double *dest, int *deriv)
-{
- double bx[6], by[6], bz[6];
- 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);
- for (int iz = -2; iz <= 3; iz++)
- for (int iy = -2; iy <= 3; iy++)
- for (int ix = -2; ix <= 3; ix++)
- for (int c=0; c<3; c++)
- dest[c] += field[((ptrdiff_t(xg[2]+iz) *this->fs->rd->subsizes[1] +
- ptrdiff_t(xg[1]+iy))*this->fs->rd->subsizes[2] +
- ptrdiff_t(xg[0]+ix))*3+c]*bz[iz+2]*by[iy+2]*bx[ix+2];
-}
-
-template <class rnumber>
-void slab_field_particles<rnumber>::spline_n3_formula(rnumber *field, int *xg, double *xx, double *dest, int *deriv)
-{
- double bx[8], by[8], bz[8];
- 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);
- for (int iz = -3; iz <= 4; iz++)
- for (int iy = -3; iy <= 4; iy++)
- for (int ix = -3; ix <= 4; ix++)
- for (int c=0; c<3; c++)
- dest[c] += field[((ptrdiff_t(xg[2]+iz) *this->fs->rd->subsizes[1] +
- ptrdiff_t(xg[1]+iy))*this->fs->rd->subsizes[2] +
- ptrdiff_t(xg[0]+ix))*3+c]*bz[iz+3]*by[iy+3]*bx[ix+3];
}
template <class rnumber>
diff --git a/bfps/cpp/slab_field_particles.hpp b/bfps/cpp/slab_field_particles.hpp
index d75e778c..834a3608 100644
--- a/bfps/cpp/slab_field_particles.hpp
+++ b/bfps/cpp/slab_field_particles.hpp
@@ -138,10 +138,7 @@ class slab_field_particles
void synchronize_single_particle_state(int p, double *x, int source_id = -1);
void get_grid_coordinates(double *x, int *xg, double *xx);
void linear_interpolation(rnumber *field, int *xg, double *xx, double *dest, int *deriv);
- void spline_formula( rnumber *field, int *xg, double *xx, double *dest, int *deriv);
- void spline_n1_formula( rnumber *field, int *xg, double *xx, double *dest, int *deriv);
- void spline_n2_formula( rnumber *field, int *xg, double *xx, double *dest, int *deriv);
- void spline_n3_formula( rnumber *field, int *xg, double *xx, double *dest, int *deriv);
+ void interpolation_formula(rnumber *field, int *xg, double *xx, double *dest, int *deriv);
void rFFTW_to_buffered(rnumber *src, rnumber *dst);
diff --git a/bfps/cpp/tracers.cpp b/bfps/cpp/tracers.cpp
index e3c53ccc..d0fd2c27 100644
--- a/bfps/cpp/tracers.cpp
+++ b/bfps/cpp/tracers.cpp
@@ -47,7 +47,7 @@ void tracers<rnumber>::jump_estimate(double *jump)
/* perform interpolation */
for (int p=0; p<this->nparticles; p++) if (this->fs->rd->myrank == this->computing[p])
{
- this->spline_formula(vel, xg + p*3, xx + p*3, tmp, deriv);
+ this->interpolation_formula(vel, 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;
@@ -79,7 +79,7 @@ void tracers<rnumber>::get_rhs(double *x, double *y)
int crank = this->get_rank(x[p*3 + 2]);
if (this->fs->rd->myrank == crank)
{
- this->spline_formula(vel, xg + p*3, xx + p*3, y + p*3, deriv);
+ this->interpolation_formula(vel, 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],
@@ -159,11 +159,12 @@ void tracers<R>::sample_vec_field(R *vec_field, double *vec_values) \
/* perform interpolation */ \
for (int p=0; p<this->nparticles; p++) \
if (this->fs->rd->myrank == this->computing[p]) \
- this->spline_formula(vec_field, \
- xg + p*3, \
- xx + p*3, \
- vec_local + p*3, \
- deriv); \
+ this->interpolation_formula( \
+ vec_field, \
+ xg + p*3, \
+ xx + p*3, \
+ vec_local + p*3, \
+ deriv); \
MPI_Allreduce( \
vec_local, \
vec_values, \
diff --git a/tests/test_base.py b/tests/test_base.py
index 3ceca37b..7fd4df46 100755
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -92,28 +92,28 @@ 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)
- #c.add_particles(
- # integration_steps = 1,
- # neighbours = opt.neighbours,
- # smoothness = opt.smoothness,
- # force_vel_reset = True)
- #c.add_particles(integration_steps = 2, neighbours = opt.neighbours, smoothness = opt.smoothness)
- #c.add_particles(integration_steps = 3, neighbours = opt.neighbours, smoothness = opt.smoothness)
- #c.add_particles(integration_steps = 4, neighbours = opt.neighbours, smoothness = opt.smoothness)
- #c.add_particles(integration_steps = 5, neighbours = opt.neighbours, smoothness = opt.smoothness)
- #c.add_particles(integration_steps = 6, neighbours = opt.neighbours, smoothness = opt.smoothness)
- #c.add_particles(integration_steps = 2, neighbours = 1, smoothness = 0)
- #c.add_particles(integration_steps = 2, neighbours = 1, smoothness = 1)
- #c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 1)
- #c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 2)
- #c.add_particles(integration_steps = 2, neighbours = 1, interp_type = 'Lagrange')
- #c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange')
+ c.add_particles(kcut = 'fs->kM/2',
+ integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness)
+ c.add_particles(
+ integration_steps = 1,
+ neighbours = opt.neighbours,
+ smoothness = opt.smoothness,
+ force_vel_reset = True)
+ c.add_particles(integration_steps = 2, neighbours = opt.neighbours, smoothness = opt.smoothness)
+ c.add_particles(integration_steps = 3, neighbours = opt.neighbours, smoothness = opt.smoothness)
+ c.add_particles(integration_steps = 4, neighbours = opt.neighbours, smoothness = opt.smoothness)
+ c.add_particles(integration_steps = 5, neighbours = opt.neighbours, smoothness = opt.smoothness)
+ c.add_particles(integration_steps = 6, neighbours = opt.neighbours, smoothness = opt.smoothness)
+ c.add_particles(integration_steps = 2, neighbours = 1, smoothness = 0)
+ c.add_particles(integration_steps = 2, neighbours = 1, smoothness = 1)
+ c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 1)
+ c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 2)
+ c.add_particles(integration_steps = 2, neighbours = 1, interp_type = 'Lagrange')
+ c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange')
c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 1, integration_method = 'Heun')
- #c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange', integration_method = 'Heun')
- #c.add_particles(integration_steps = 4, neighbours = 6, smoothness = 1, integration_method = 'cRK4')
- #c.add_particles(integration_steps = 4, neighbours = 6, interp_type = 'Lagrange', integration_method = 'cRK4')
+ c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange', integration_method = 'Heun')
+ c.add_particles(integration_steps = 4, neighbours = 6, smoothness = 1, integration_method = 'cRK4')
+ c.add_particles(integration_steps = 4, neighbours = 6, interp_type = 'Lagrange', integration_method = 'cRK4')
c.fill_up_fluid_code()
c.finalize_code()
c.write_src()
--
GitLab