Skip to content
Snippets Groups Projects
Commit fda98c1a authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

clean up particle code

parent c0ec5a97
Branches
Tags
No related merge requests found
...@@ -536,10 +536,9 @@ void slab_field_particles<rnumber>::get_grid_coordinates(double *x, int *xg, dou ...@@ -536,10 +536,9 @@ void slab_field_particles<rnumber>::get_grid_coordinates(double *x, int *xg, dou
} }
template <class rnumber> 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]; 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[0], xx[0], bx);
this->compute_beta(deriv[1], xx[1], by); this->compute_beta(deriv[1], xx[1], by);
this->compute_beta(deriv[2], xx[2], bz); this->compute_beta(deriv[2], xx[2], bz);
...@@ -564,58 +563,6 @@ void slab_field_particles<rnumber>::spline_formula(rnumber *field, int *xg, doub ...@@ -564,58 +563,6 @@ void slab_field_particles<rnumber>::spline_formula(rnumber *field, int *xg, doub
by[iy+this->interp_neighbours]* by[iy+this->interp_neighbours]*
bx[ix+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> template <class rnumber>
......
...@@ -138,10 +138,7 @@ class slab_field_particles ...@@ -138,10 +138,7 @@ class slab_field_particles
void synchronize_single_particle_state(int p, double *x, int source_id = -1); void synchronize_single_particle_state(int p, double *x, int source_id = -1);
void get_grid_coordinates(double *x, int *xg, double *xx); void get_grid_coordinates(double *x, int *xg, double *xx);
void linear_interpolation(rnumber *field, int *xg, double *xx, double *dest, int *deriv); 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 interpolation_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 rFFTW_to_buffered(rnumber *src, rnumber *dst); void rFFTW_to_buffered(rnumber *src, rnumber *dst);
......
...@@ -47,7 +47,7 @@ void tracers<rnumber>::jump_estimate(double *jump) ...@@ -47,7 +47,7 @@ void tracers<rnumber>::jump_estimate(double *jump)
/* perform interpolation */ /* perform interpolation */
for (int p=0; p<this->nparticles; p++) if (this->fs->rd->myrank == this->computing[p]) 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]); jump[p] = fabs(3*this->dt * tmp[2]);
if (jump[p] < this->dz*1.01) if (jump[p] < this->dz*1.01)
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) ...@@ -79,7 +79,7 @@ void tracers<rnumber>::get_rhs(double *x, double *y)
int crank = this->get_rank(x[p*3 + 2]); int crank = this->get_rank(x[p*3 + 2]);
if (this->fs->rd->myrank == crank) 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( DEBUG_MSG(
"position is %g %g %g %d %d %d %g %g %g, result is %g %g %g\n", "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], 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) \ ...@@ -159,11 +159,12 @@ void tracers<R>::sample_vec_field(R *vec_field, double *vec_values) \
/* perform interpolation */ \ /* perform interpolation */ \
for (int p=0; p<this->nparticles; p++) \ for (int p=0; p<this->nparticles; p++) \
if (this->fs->rd->myrank == this->computing[p]) \ if (this->fs->rd->myrank == this->computing[p]) \
this->spline_formula(vec_field, \ this->interpolation_formula( \
xg + p*3, \ vec_field, \
xx + p*3, \ xg + p*3, \
vec_local + p*3, \ xx + p*3, \
deriv); \ vec_local + p*3, \
deriv); \
MPI_Allreduce( \ MPI_Allreduce( \
vec_local, \ vec_local, \
vec_values, \ vec_values, \
......
...@@ -92,28 +92,28 @@ def launch( ...@@ -92,28 +92,28 @@ def launch(
c.parameters['niter_part'] = 1 c.parameters['niter_part'] = 1
c.parameters['famplitude'] = 0.2 c.parameters['famplitude'] = 0.2
if c.parameters['nparticles'] > 0: if c.parameters['nparticles'] > 0:
#c.add_particles(kcut = 'fs->kM/2', c.add_particles(kcut = 'fs->kM/2',
# integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness) integration_steps = 1, neighbours = opt.neighbours, smoothness = opt.smoothness)
#c.add_particles( c.add_particles(
# integration_steps = 1, integration_steps = 1,
# neighbours = opt.neighbours, neighbours = opt.neighbours,
# smoothness = opt.smoothness, smoothness = opt.smoothness,
# force_vel_reset = True) force_vel_reset = True)
#c.add_particles(integration_steps = 2, neighbours = opt.neighbours, smoothness = opt.smoothness) 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 = 3, neighbours = opt.neighbours, smoothness = opt.smoothness)
#c.add_particles(integration_steps = 4, 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 = 5, neighbours = opt.neighbours, smoothness = opt.smoothness)
#c.add_particles(integration_steps = 6, 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 = 0)
#c.add_particles(integration_steps = 2, neighbours = 1, smoothness = 1) 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 = 1)
#c.add_particles(integration_steps = 2, neighbours = 6, smoothness = 2) 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 = 1, interp_type = 'Lagrange')
#c.add_particles(integration_steps = 2, neighbours = 6, 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, smoothness = 1, integration_method = 'Heun')
#c.add_particles(integration_steps = 2, neighbours = 6, interp_type = 'Lagrange', 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, smoothness = 1, integration_method = 'cRK4')
#c.add_particles(integration_steps = 4, neighbours = 6, interp_type = 'Lagrange', integration_method = 'cRK4') c.add_particles(integration_steps = 4, neighbours = 6, interp_type = 'Lagrange', integration_method = 'cRK4')
c.fill_up_fluid_code() c.fill_up_fluid_code()
c.finalize_code() c.finalize_code()
c.write_src() c.write_src()
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment