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

add linear_interpolation method

parent 50e4e414
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
bfps/test.py
+ 0
− 1
View file @ d0bda39b
...@@ -116,7 +116,6 @@ class convergence_test(bfps.code): ...@@ -116,7 +116,6 @@ class convergence_test(bfps.code):
ps->iteration++; ps->iteration++;
ps->synchronize(); ps->synchronize();
do_stats(fs, ps); do_stats(fs, ps);
DEBUG_MSG("iteration %d hello\\n", fs->iteration);
} }
if (myrank == 0) if (myrank == 0)
{ {
......
src/tracers.cpp
+ 21
− 21
View file @ d0bda39b
...@@ -24,6 +24,20 @@ ...@@ -24,6 +24,20 @@
#include "fftw_tools.hpp" #include "fftw_tools.hpp"
#include "tracers.hpp" #include "tracers.hpp"
template <class rnumber>
void tracers<rnumber>::linear_interpolation(float *field, int *xg, double *xx, double *dest)
{
for (int c=0; c<3; c++)
dest[c] = (field[((ptrdiff_t(xg[2] )*this->fs->rd->subsizes[1]+xg[1] )*this->fs->rd->subsizes[0]+xg[0] )*3+c]*((1-xx[0])*(1-xx[1])*(1-xx[2])) +
field[((ptrdiff_t(xg[2] )*this->fs->rd->subsizes[1]+xg[1] )*this->fs->rd->subsizes[0]+xg[0]+1)*3+c]*(( xx[0])*(1-xx[1])*(1-xx[2])) +
field[((ptrdiff_t(xg[2] )*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[0]+xg[0] )*3+c]*((1-xx[0])*( xx[1])*(1-xx[2])) +
field[((ptrdiff_t(xg[2] )*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[0]+xg[0]+1)*3+c]*(( xx[0])*( xx[1])*(1-xx[2])) +
field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1] )*this->fs->rd->subsizes[0]+xg[0] )*3+c]*((1-xx[0])*(1-xx[1])*( xx[2])) +
field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1] )*this->fs->rd->subsizes[0]+xg[0]+1)*3+c]*(( xx[0])*(1-xx[1])*( xx[2])) +
field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[0]+xg[0] )*3+c]*((1-xx[0])*( xx[1])*( xx[2])) +
field[((ptrdiff_t(xg[2]+1)*this->fs->rd->subsizes[1]+xg[1]+1)*this->fs->rd->subsizes[0]+xg[0]+1)*3+c]*(( xx[0])*( xx[1])*( xx[2])));
}
template <class rnumber> template <class rnumber>
void tracers<rnumber>::jump_estimate(double *jump) void tracers<rnumber>::jump_estimate(double *jump)
{ {
...@@ -32,20 +46,14 @@ void tracers<rnumber>::jump_estimate(double *jump) ...@@ -32,20 +46,14 @@ void tracers<rnumber>::jump_estimate(double *jump)
int *xg = new int[this->array_size]; int *xg = new int[this->array_size];
double *xx = new double[this->array_size]; double *xx = new double[this->array_size];
float *vel = this->data + this->buffer_size*this->fs->rd->slice_size; float *vel = this->data + this->buffer_size*this->fs->rd->slice_size;
double tmp[3];
this->get_grid_coordinates(this->state, xg, xx); this->get_grid_coordinates(this->state, xg, xx);
; ;
/* perform interpolation */ /* perform interpolation */
for (int p=0; p<this->nparticles; p++) if (this->is_active[this->fs->rd->myrank][p]) for (int p=0; p<this->nparticles; p++) if (this->is_active[this->fs->rd->myrank][p])
{ {
jump[p] = fabs(2*this->dt * \ this->linear_interpolation(vel, xg + p*3, xx + p*3, tmp);
(vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3] )*3+2]*((1-xx[p*3 ])*(1-xx[p*3+1])*(1-xx[p*3+2])) + jump[p] = fabs(2*this->dt * tmp[2]);
vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+2]*(( xx[p*3 ])*(1-xx[p*3+1])*(1-xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3] )*3+2]*((1-xx[p*3 ])*( xx[p*3+1])*(1-xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+2]*(( xx[p*3 ])*( xx[p*3+1])*(1-xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3] )*3+2]*((1-xx[p*3 ])*(1-xx[p*3+1])*( xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+2]*(( xx[p*3 ])*(1-xx[p*3+1])*( xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3] )*3+2]*((1-xx[p*3 ])*( xx[p*3+1])*( xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+2]*(( xx[p*3 ])*( xx[p*3+1])*( xx[p*3+2]))));
} }
delete[] xg; delete[] xg;
delete[] xx; delete[] xx;
...@@ -62,17 +70,9 @@ void tracers<rnumber>::get_rhs(double *x, double *y) ...@@ -62,17 +70,9 @@ void tracers<rnumber>::get_rhs(double *x, double *y)
this->get_grid_coordinates(x, xg, xx); this->get_grid_coordinates(x, xg, xx);
/* perform interpolation */ /* perform interpolation */
for (int p=0; p<this->nparticles; p++) if (this->is_active[this->fs->rd->myrank][p]) for (int p=0; p<this->nparticles; p++) if (this->is_active[this->fs->rd->myrank][p])
for (int c=0; c<3; c++)
{ {
y[p*3+c] = (vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3] )*3+c]*((1-xx[p*3 ])*(1-xx[p*3+1])*(1-xx[p*3+2])) + this->linear_interpolation(vel, xg + p*3, xx + p*3, y + p*3);
vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+c]*(( xx[p*3 ])*(1-xx[p*3+1])*(1-xx[p*3+2])) + DEBUG_MSG("particle %d found y %lg %lg %lg\n", p, y[p*3], y[p*3+1], y[p*3+2]);
vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3] )*3+c]*((1-xx[p*3 ])*( xx[p*3+1])*(1-xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2] )*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+c]*(( xx[p*3 ])*( xx[p*3+1])*(1-xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3] )*3+c]*((1-xx[p*3 ])*(1-xx[p*3+1])*( xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1] )*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+c]*(( xx[p*3 ])*(1-xx[p*3+1])*( xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3] )*3+c]*((1-xx[p*3 ])*( xx[p*3+1])*( xx[p*3+2])) +
vel[((ptrdiff_t(xg[p*3+2]+1)*this->fs->rd->subsizes[1]+xg[p*3+1]+1)*this->fs->rd->subsizes[0]+xg[p*3]+1)*3+c]*(( xx[p*3 ])*( xx[p*3+1])*( xx[p*3+2])));
DEBUG_MSG("particle %d found y %lg\n", p, y[p*3+c]);
} }
delete[] xg; delete[] xg;
delete[] xx; delete[] xx;
......
src/tracers.hpp
+ 1
− 0
View file @ d0bda39b
...@@ -47,6 +47,7 @@ class tracers:public slab_field_particles<rnumber> ...@@ -47,6 +47,7 @@ class tracers:public slab_field_particles<rnumber>
void update_field(bool clip_on = true); void update_field(bool clip_on = true);
virtual void get_rhs(double *x, double *rhs); virtual void get_rhs(double *x, double *rhs);
virtual void jump_estimate(double *jump_length); virtual void jump_estimate(double *jump_length);
void linear_interpolation(float *field, int *xg, double *xx, double *dest);
}; };
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment