Commit 78b0ebbc authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'develop'

parents db2b5427 a88bb74c
......@@ -6,3 +6,7 @@ meta/.ipynb_checkpoints/*
dist
build
bfps.egg-info
MANIFEST.in
bfps/install_info.pickle
bfps/lib*.a
obj
......@@ -408,8 +408,7 @@ 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}->temp);\n' +
'acc_{0}->read_rFFTW(acc_{0}->temp);\n').format(name)
'fs->compute_Lagrangian_acceleration(acc_{0}->field);\n').format(name)
self.fluid_start += update_fields
self.fluid_loop += update_fields
return None
......@@ -450,7 +449,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
double *velocity = new double[ps{0}->array_size];
ps{0}->sample_vec_field(vel_{1}, velocity);
ps{0}->sample_vec_field(acc_{1}, acceleration);
if (ps{0}->fs->rd->myrank == 0)
if (ps{0}->myrank == 0)
{{
//VELOCITY begin
std::string temp_string = (std::string("/particles/") +
......@@ -512,7 +511,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
neighbours,
self.particle_species)
self.particle_start += ('ps{0} = new rFFTW_particles<VELOCITY_TRACER, {1}, {2}>(\n' +
'fname, fs, vel_{3},\n' +
'fname, vel_{3},\n' +
'nparticles,\n' +
'niter_part, tracers{0}_integration_steps);\n').format(
self.particle_species, self.C_dtype, neighbours, fields_name)
......@@ -536,10 +535,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.particle_loop += 'ps{0}->synchronize();\n'.format(self.particle_species)
elif particle_class == 'rFFTW_particles':
self.particle_loop += 'ps{0}->step();\n'.format(self.particle_species)
self.particle_loop += (('if (ps{0}->iteration % niter_part == 0)\n' +
'{{\n' +
'ps{0}->write(stat_file, false);\n').format(self.particle_species) +
output_vel_acc + '}\n')
self.particle_stat_src += (
('if (ps{0}->iteration % niter_part == 0)\n' +
'{{\n' +
'ps{0}->write(stat_file, false);\n').format(self.particle_species) +
output_vel_acc + '}\n')
self.particle_species += 1
return None
def get_data_file(self):
......
......@@ -26,6 +26,7 @@
#define NDEBUG
#include <cmath>
#include "rFFTW_interpolator.hpp"
template <class rnumber, int interp_neighbours>
......@@ -37,58 +38,104 @@ rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
this->field_size = 2*fs->cd->local_size;
this->compute_beta = BETA_POLYS;
if (sizeof(rnumber) == 4)
{
this->f0 = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
this->f1 = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
}
this->field = (rnumber*)((void*)fftwf_alloc_real(this->field_size));
else if (sizeof(rnumber) == 8)
{
this->f0 = (rnumber*)((void*)fftw_alloc_real(this->field_size));
this->f1 = (rnumber*)((void*)fftw_alloc_real(this->field_size));
}
this->temp = this->f1;
this->field = (rnumber*)((void*)fftw_alloc_real(this->field_size));
// compute dx, dy, dz;
this->dx = 4*acos(0) / (fs->dkx*this->descriptor->sizes[2]);
this->dy = 4*acos(0) / (fs->dky*this->descriptor->sizes[1]);
this->dz = 4*acos(0) / (fs->dkz*this->descriptor->sizes[0]);
// generate compute array
this->compute = new bool[this->descriptor->sizes[0]];
std::fill_n(this->compute, this->descriptor->sizes[0], false);
for (int iz = this->descriptor->starts[0]-interp_neighbours-1;
iz <= this->descriptor->starts[0]+this->descriptor->subsizes[0]+interp_neighbours;
iz++)
this->compute[((iz + this->descriptor->sizes[0]) % this->descriptor->sizes[0])] = true;
}
template <class rnumber, int interp_neighbours>
rFFTW_interpolator<rnumber, interp_neighbours>::~rFFTW_interpolator()
{
if (sizeof(rnumber) == 4)
{
fftwf_free((float*)((void*)this->f0));
fftwf_free((float*)((void*)this->f1));
}
fftwf_free((float*)((void*)this->field));
else if (sizeof(rnumber) == 8)
{
fftw_free((double*)((void*)this->f0));
fftw_free((double*)((void*)this->f1));
}
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)
{
/* first, roll fields */
rnumber *tmp = this->f0;
this->f0 = this->f1;
this->f1 = tmp;
this->temp = this->f0;
/* now do regular things */
rnumber *src = (rnumber*)void_src;
rnumber *dst = this->f1;
/* do big copy of middle stuff */
std::copy(src,
src + this->field_size,
dst);
this->field);
return EXIT_SUCCESS;
}
template <class rnumber, int interp_neighbours>
void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
double t,
void rFFTW_interpolator<rnumber, interp_neighbours>::get_grid_coordinates(
const int nparticles,
const int pdimension,
const double *x,
int *xg,
double *xx,
double *xx)
{
static double grid_size[] = {this->dx, this->dy, this->dz};
double tval;
std::fill_n(xg, nparticles*3, 0);
std::fill_n(xx, nparticles*3, 0.0);
for (int p=0; p<nparticles; p++)
{
for (int c=0; c<3; c++)
{
tval = floor(x[p*pdimension+c]/grid_size[c]);
xg[p*3+c] = MOD(int(tval), this->descriptor->sizes[2-c]);
xx[p*3+c] = (x[p*pdimension+c] - tval*grid_size[c]) / grid_size[c];
}
}
}
template <class rnumber, int interp_neighbours>
void rFFTW_interpolator<rnumber, interp_neighbours>::sample(
const int nparticles,
const int pdimension,
const double *__restrict__ x,
double *__restrict__ y,
const int *deriv)
{
/* get grid coordinates */
int *xg = new int[3*nparticles];
double *xx = new double[3*nparticles];
double *yy = new double[3*nparticles];
std::fill_n(yy, 3*nparticles, 0.0);
this->get_grid_coordinates(nparticles, pdimension, x, xg, xx);
/* perform interpolation */
for (int p=0; p<nparticles; p++)
if (this->compute[xg[p*3+2]])
this->operator()(xg + p*3, xx + p*3, yy + p*3, deriv);
MPI_Allreduce(
yy,
y,
3*nparticles,
MPI_DOUBLE,
MPI_SUM,
this->descriptor->comm);
delete[] yy;
delete[] xg;
delete[] xx;
}
template <class rnumber, int interp_neighbours>
void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
const int *xg,
const double *xx,
double *dest,
int *deriv)
const int *deriv)
{
double bx[interp_neighbours*2+2], by[interp_neighbours*2+2], bz[interp_neighbours*2+2];
if (deriv == NULL)
......@@ -105,13 +152,11 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
}
std::fill_n(dest, 3, 0);
ptrdiff_t bigiz, bigiy, bigix;
double tval[3];
for (int iz = -interp_neighbours; iz <= interp_neighbours+1; iz++)
{
bigiz = ptrdiff_t(((xg[2]+iz) + this->descriptor->sizes[0]) % this->descriptor->sizes[0]);
if (this->descriptor->myrank == this->descriptor->rank[bigiz])
{
std::fill_n(tval, 3, 0);
for (int iy = -interp_neighbours; iy <= interp_neighbours+1; iy++)
{
bigiy = ptrdiff_t(MOD(xg[1]+iy, this->descriptor->sizes[1]));
......@@ -122,17 +167,11 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
bigiy)*(this->descriptor->sizes[2]+2) +
bigix)*3;
for (int c=0; c<3; c++)
{
dest[c] += (this->f0[tindex+c]*(1-t) + t*this->f1[tindex+c])*(bz[iz+interp_neighbours]*
by[iy+interp_neighbours]*
bx[ix+interp_neighbours]);
tval[c] += (this->f0[tindex+c]*(1-t) + t*this->f1[tindex+c])*(bz[iz+interp_neighbours]*
by[iy+interp_neighbours]*
bx[ix+interp_neighbours]);
}
dest[c] += this->field[tindex+c]*(bz[iz+interp_neighbours]*
by[iy+interp_neighbours]*
bx[ix+interp_neighbours]);
}
}
DEBUG_MSG("%ld %d %d %g %g %g\n", bigiz, xg[1], xg[0], tval[0], tval[1], tval[2]);
}
}
}
......
......@@ -44,17 +44,52 @@ template <class rnumber, int interp_neighbours>
class rFFTW_interpolator
{
public:
/* size of field to interpolate */
ptrdiff_t field_size;
/* pointer to polynomial function */
base_polynomial_values compute_beta;
/* descriptor of field to interpolate */
field_descriptor<rnumber> *descriptor;
rnumber *f0, *f1, *temp;
/* pointers to fields that are to be interpolated
* */
rnumber *field;
/* physical parameters of field */
double dx, dy, dz;
/* compute[iz] is true if .
* local_zstart - neighbours <= iz <= local_zend + 1 + neighbours
* */
bool *compute;
rFFTW_interpolator(
fluid_solver_base<rnumber> *FSOLVER,
base_polynomial_values BETA_POLYS);
~rFFTW_interpolator();
void operator()(double t, int *__restrict__ xg, double *__restrict__ xx, double *__restrict__ dest, int *deriv = NULL);
/* map real locations to grid coordinates */
void get_grid_coordinates(
const int nparticles,
const int pdimension,
const double *__restrict__ x,
int *__restrict__ xg,
double *__restrict__ xx);
/* interpolate field at an array of locations */
void sample(
const int nparticles,
const int pdimension,
const double *__restrict__ x,
double *__restrict__ y,
const int *deriv = NULL);
/* interpolate 1 point */
void operator()(
const int *__restrict__ xg,
const double *__restrict__ xx,
double *__restrict__ dest,
const int *deriv = NULL);
int read_rFFTW(void *src);
};
......
......@@ -24,7 +24,7 @@
//#define NDEBUG
#define NDEBUG
#include <cmath>
#include <cassert>
......@@ -42,7 +42,6 @@ extern int myrank, nprocs;
template <int particle_type, class rnumber, int interp_neighbours>
rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
const char *NAME,
fluid_solver_base<rnumber> *FSOLVER,
rFFTW_interpolator<rnumber, interp_neighbours> *FIELD,
const int NPARTICLES,
const int TRAJ_SKIP,
......@@ -60,7 +59,6 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
assert((INTEGRATION_STEPS <= 6) &&
(INTEGRATION_STEPS >= 1));
strncpy(this->name, NAME, 256);
this->fs = FSOLVER;
this->nparticles = NPARTICLES;
this->vel = FIELD;
this->integration_steps = INTEGRATION_STEPS;
......@@ -76,41 +74,6 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
this->rhs[i] = new double[this->array_size];
std::fill_n(this->rhs[i], this->array_size, 0.0);
}
// compute dx, dy, dz;
this->dx = 4*acos(0) / (this->fs->dkx*this->fs->rd->sizes[2]);
this->dy = 4*acos(0) / (this->fs->dky*this->fs->rd->sizes[1]);
this->dz = 4*acos(0) / (this->fs->dkz*this->fs->rd->sizes[0]);
// compute lower and upper bounds
this->lbound = new double[nprocs];
this->ubound = new double[nprocs];
double *tbound = new double[nprocs];
std::fill_n(tbound, nprocs, 0.0);
tbound[this->myrank] = this->fs->rd->starts[0]*this->dz;
MPI_Allreduce(
tbound,
this->lbound,
nprocs,
MPI_DOUBLE,
MPI_SUM,
this->comm);
std::fill_n(tbound, nprocs, 0.0);
tbound[this->myrank] = (this->fs->rd->starts[0] + this->fs->rd->subsizes[0])*this->dz;
MPI_Allreduce(
tbound,
this->ubound,
nprocs,
MPI_DOUBLE,
MPI_SUM,
this->comm);
delete[] tbound;
//for (int r = 0; r<nprocs; r++)
// DEBUG_MSG(
// "lbound[%d] = %lg, ubound[%d] = %lg\n",
// r, this->lbound[r],
// r, this->ubound[r]
// );
}
template <int particle_type, class rnumber, int interp_neighbours>
......@@ -121,59 +84,19 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::~rFFTW_particles()
{
delete[] this->rhs[i];
}
delete[] this->lbound;
delete[] this->ubound;
}
template <int particle_type, class rnumber, int interp_neighbours>
void rFFTW_particles<particle_type, rnumber, interp_neighbours>::sample_vec_field(
rFFTW_interpolator<rnumber, interp_neighbours> *vec,
double t,
double *x,
double *y,
int *deriv)
{
/* get grid coordinates */
int *xg = new int[3*this->nparticles];
double *xx = new double[3*this->nparticles];
double *yy = new double[3*this->nparticles];
std::fill_n(yy, 3*this->nparticles, 0.0);
this->get_grid_coordinates(x, xg, xx);
/* perform interpolation */
for (int p=0; p<this->nparticles; p++)
(*vec)(t, xg + p*3, xx + p*3, yy + p*3, deriv);
MPI_Allreduce(
yy,
y,
3*this->nparticles,
MPI_DOUBLE,
MPI_SUM,
this->comm);
delete[] yy;
delete[] xg;
delete[] xx;
}
template <int particle_type, class rnumber, int interp_neighbours>
void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rhs(double t, double *x, double *y)
void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y)
{
switch(particle_type)
{
case VELOCITY_TRACER:
this->sample_vec_field(this->vel, t, x, y);
this->vel->sample(this->nparticles, this->ncomponents, x, y);
break;
}
}
template <int particle_type, class rnumber, int interp_neighbours>
int rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rank(double z)
{
int tmp = this->fs->rd->rank[MOD(int(floor(z/this->dz)), this->fs->rd->sizes[0])];
assert(tmp >= 0 && tmp < this->fs->rd->nprocs);
return tmp;
}
template <int particle_type, class rnumber, int interp_neighbours>
void rFFTW_particles<particle_type, rnumber, interp_neighbours>::roll_rhs()
{
......@@ -189,7 +112,7 @@ template <int particle_type, class rnumber, int interp_neighbours>
void rFFTW_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth(int nsteps)
{
ptrdiff_t ii;
this->get_rhs(0, this->state, this->rhs[0]);
this->get_rhs(this->state, this->rhs[0]);
switch(nsteps)
{
case 1:
......@@ -270,29 +193,6 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::step()
}
template <int particle_type, class rnumber, int interp_neighbours>
void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_grid_coordinates(double *x, int *xg, double *xx)
{
static double grid_size[] = {this->dx, this->dy, this->dz};
double tval;
std::fill_n(xg, this->nparticles*3, 0);
std::fill_n(xx, this->nparticles*3, 0.0);
for (int p=0; p<this->nparticles; p++)
{
for (int c=0; c<3; c++)
{
tval = floor(x[p*this->ncomponents+c]/grid_size[c]);
xg[p*3+c] = MOD(int(tval), this->fs->rd->sizes[2-c]);
xx[p*3+c] = (x[p*this->ncomponents+c] - tval*grid_size[c]) / grid_size[c];
}
/*xg[p*3+2] -= this->fs->rd->starts[0];
if (this->myrank == this->fs->rd->rank[0] &&
xg[p*3+2] > this->fs->rd->subsizes[0])
xg[p*3+2] -= this->fs->rd->sizes[0];*/
}
}
template <int particle_type, class rnumber, int interp_neighbours>
void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read(hid_t data_file_id)
{
......
......@@ -43,20 +43,8 @@ class rFFTW_particles
public:
int myrank, nprocs;
MPI_Comm comm;
fluid_solver_base<rnumber> *fs;
rnumber *data;
/* watching is an array of shape [nparticles], with
* watching[p] being true if particle p is in the domain of myrank
* or in the buffer regions.
* only used if multistep is false.
* */
bool *watching;
/* computing is an array of shape [nparticles], with
* computing[p] being the rank that is currently working on particle p
* */
int *computing;
/* state will generally hold all the information about the particles.
* in the beginning, we will only need to solve 3D ODEs, but I figured
* a general ncomponents is better, since we may change our minds.
......@@ -68,8 +56,6 @@ class rFFTW_particles
int array_size;
int integration_steps;
int traj_skip;
double *lbound;
double *ubound;
rFFTW_interpolator<rnumber, interp_neighbours> *vel;
/* simulation parameters */
......@@ -77,23 +63,15 @@ class rFFTW_particles
int iteration;
double dt;
/* physical parameters of field */
double dx, dy, dz;
/* methods */
/* constructor and destructor.
* allocate and deallocate:
* this->state
* this->rhs
* this->lbound
* this->ubound
* this->computing
* this->watching
* */
rFFTW_particles(
const char *NAME,
fluid_solver_base<rnumber> *FSOLVER,
rFFTW_interpolator<rnumber, interp_neighbours> *FIELD,
const int NPARTICLES,
const int TRAJ_SKIP,
......@@ -101,19 +79,10 @@ class rFFTW_particles
~rFFTW_particles();
void get_rhs(double *__restrict__ x, double *__restrict__ rhs);
void get_rhs(double t, double *__restrict__ x, double *__restrict__ rhs);
int get_rank(double z); // get rank for given value of z
void get_grid_coordinates(double *__restrict__ x, int *__restrict__ xg, double *__restrict__ xx);
void sample_vec_field(
rFFTW_interpolator<rnumber, interp_neighbours> *vec,
double t,
double *__restrict__ x,
double *__restrict__ y,
int *deriv = NULL);
inline void sample_vec_field(rFFTW_interpolator<rnumber, interp_neighbours> *field, double *vec_values)
{
this->sample_vec_field(field, 1.0, this->state, vec_values, NULL);
field->sample(this->nparticles, this->ncomponents, this->state, vec_values, NULL);
}
/* input/output */
......
......@@ -90,13 +90,14 @@ class fluid_particle_base(bfps.code):
self.fluid_loop = ''
self.fluid_end = ''
self.fluid_output = ''
self.stat_src = ''
self.particle_includes = ''
self.particle_variables = ''
self.particle_definitions = ''
self.particle_start = ''
self.particle_loop = ''
self.particle_end = ''
self.stat_src = ''
self.particle_stat_src = ''
self.file_datasets_grow = ''
return None
def finalize_code(self):
......@@ -144,6 +145,7 @@ class fluid_particle_base(bfps.code):
'return file_problems;\n'
'}\n')
self.definitions += 'void do_stats()\n{\n' + self.stat_src + '}\n'
self.definitions += 'void do_particle_stats()\n{\n' + self.particle_stat_src + '}\n'
# take care of wisdom
if self.use_fftw_wisdom:
if self.dtype == np.float32:
......@@ -191,6 +193,7 @@ class fluid_particle_base(bfps.code):
return EXIT_SUCCESS;
}
do_stats();
do_particle_stats();
//endcpp
"""
output_time_difference = ('time1 = clock();\n' +
......@@ -202,16 +205,21 @@ class fluid_particle_base(bfps.code):
'<< iteration << " took " ' +
'<< time_difference/nprocs << " seconds" << std::endl;\n' +
'time0 = time1;\n')
self.main += 'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)\n{\n'
self.main += 'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)\n'
self.main += '{\n'
self.main += output_time_difference
self.main += self.fluid_loop
if self.particle_species > 0:
self.main += self.particle_loop
self.main += 'if (iteration % niter_stat == 0) do_stats();\n}\n'
self.main += self.fluid_loop
self.main += 'if (iteration % niter_stat == 0) do_stats();\n'
if self.particle_species > 0:
self.main += 'if (iteration % niter_part == 0) do_particle_stats();\n'
self.main += '}\n'
self.main += output_time_difference
self.main += 'do_stats();\n'
self.main += 'do_particle_stats();\n'
if self.particle_species > 0:
self.main += self.particle_end
self.main += 'do_stats();\n'
self.main += self.fluid_end
return None
def read_rfield(
......
x 2015-12-04 make code py3 compatible @python3
x 2015-12-23 decide on versioning system +merge0
x 2015-12-24 move get_grid coords to interpolator @optimization +v1.0
x 2015-12-25 get rid of temporal interpolation @optimization +v1.0
x 2015-12-26 call interpolation only when needed @optimization +v1.0
......@@ -25,9 +25,13 @@
import os
import sys
import subprocess
import pickle
import numpy as np
import matplotlib.pyplot as plt
import bfps
from bfps import fluid_resize
......@@ -211,6 +215,93 @@ def acceleration_test(c, m = 3, species = 0):