diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py index a9a940a55b1e3d5b38da62c32d8dae18e99a278a..9fd928033910fa43327a1064316dc0c680276037 100644 --- a/bfps/NavierStokes.py +++ b/bfps/NavierStokes.py @@ -440,7 +440,7 @@ class NavierStokes(_fluid_particle_base): interpolator = 'field_interpolator', frozen_particles = False, acc_name = None, - class_name = 'rFFTW_particles'): + class_name = 'particles'): """Adds code for tracking a series of particle species, each consisting of `nparticles` particles. diff --git a/bfps/cpp/rFFTW_particles.cpp b/bfps/cpp/rFFTW_particles.cpp deleted file mode 100644 index 70db8a86cd79dfd67ffeafee50fbbe1dd0088787..0000000000000000000000000000000000000000 --- a/bfps/cpp/rFFTW_particles.cpp +++ /dev/null @@ -1,341 +0,0 @@ -/********************************************************************** -* * -* Copyright 2015 Max Planck Institute * -* for Dynamics and Self-Organization * -* * -* This file is part of bfps. * -* * -* bfps is free software: you can redistribute it and/or modify * -* it under the terms of the GNU General Public License as published * -* by the Free Software Foundation, either version 3 of the License, * -* or (at your option) any later version. * -* * -* bfps is distributed in the hope that it will be useful, * -* but WITHOUT ANY WARRANTY; without even the implied warranty of * -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * -* GNU General Public License for more details. * -* * -* You should have received a copy of the GNU General Public License * -* along with bfps. If not, see <http://www.gnu.org/licenses/> * -* * -* Contact: Cristian.Lalescu@ds.mpg.de * -* * -**********************************************************************/ - - - -#define NDEBUG - -#include <cmath> -#include <cassert> -#include <cstring> -#include <string> -#include <sstream> - -#include "base.hpp" -#include "rFFTW_particles.hpp" -#include "fftw_tools.hpp" - - -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, - const hid_t data_file_id, - rFFTW_interpolator<rnumber, interp_neighbours> *FIELD, - const int NPARTICLES, - const int TRAJ_SKIP, - const int INTEGRATION_STEPS) -{ - switch(particle_type) - { - case VELOCITY_TRACER: - this->ncomponents = 3; - break; - default: - this->ncomponents = 3; - break; - } - assert((INTEGRATION_STEPS <= 6) && - (INTEGRATION_STEPS >= 1)); - strncpy(this->name, NAME, 256); - this->nparticles = NPARTICLES; - this->vel = FIELD; - this->integration_steps = INTEGRATION_STEPS; - this->traj_skip = TRAJ_SKIP; - this->myrank = this->vel->descriptor->myrank; - this->nprocs = this->vel->descriptor->nprocs; - this->comm = this->vel->descriptor->comm; - this->array_size = this->nparticles * this->ncomponents; - this->state = new double[this->array_size]; - std::fill_n(this->state, this->array_size, 0.0); - for (int i=0; i < this->integration_steps; i++) - { - this->rhs[i] = new double[this->array_size]; - std::fill_n(this->rhs[i], this->array_size, 0.0); - } -} - -template <int particle_type, class rnumber, int interp_neighbours> -rFFTW_particles<particle_type, rnumber, interp_neighbours>::~rFFTW_particles() -{ - delete[] this->state; - for (int i=0; i < this->integration_steps; i++) - { - delete[] this->rhs[i]; - } -} - -template <int particle_type, class rnumber, int interp_neighbours> -void rFFTW_particles<particle_type, rnumber, interp_neighbours>::get_rhs(double *x, double *y) -{ - switch(particle_type) - { - case VELOCITY_TRACER: - this->vel->sample(this->nparticles, this->ncomponents, x, y); - break; - } -} - -template <int particle_type, class rnumber, int interp_neighbours> -void rFFTW_particles<particle_type, rnumber, interp_neighbours>::roll_rhs() -{ - for (int i=this->integration_steps-2; i>=0; i--) - std::copy(this->rhs[i], - this->rhs[i] + this->array_size, - this->rhs[i+1]); -} - - - -template <int particle_type, class rnumber, int interp_neighbours> -void rFFTW_particles<particle_type, rnumber, interp_neighbours>::AdamsBashforth( - const int nsteps) -{ - ptrdiff_t ii; - this->get_rhs(this->state, this->rhs[0]); - switch(nsteps) - { - case 1: - for (int p=0; p<this->nparticles; p++) - for (int i=0; i<this->ncomponents; i++) - { - ii = p*this->ncomponents+i; - this->state[ii] += this->dt*this->rhs[0][ii]; - } - break; - case 2: - for (int p=0; p<this->nparticles; p++) - for (int i=0; i<this->ncomponents; i++) - { - ii = p*this->ncomponents+i; - this->state[ii] += this->dt*(3*this->rhs[0][ii] - - this->rhs[1][ii])/2; - } - break; - case 3: - for (int p=0; p<this->nparticles; p++) - for (int i=0; i<this->ncomponents; i++) - { - ii = p*this->ncomponents+i; - this->state[ii] += this->dt*(23*this->rhs[0][ii] - - 16*this->rhs[1][ii] - + 5*this->rhs[2][ii])/12; - } - break; - case 4: - for (int p=0; p<this->nparticles; p++) - for (int i=0; i<this->ncomponents; i++) - { - ii = p*this->ncomponents+i; - this->state[ii] += this->dt*(55*this->rhs[0][ii] - - 59*this->rhs[1][ii] - + 37*this->rhs[2][ii] - - 9*this->rhs[3][ii])/24; - } - break; - case 5: - for (int p=0; p<this->nparticles; p++) - for (int i=0; i<this->ncomponents; i++) - { - ii = p*this->ncomponents+i; - this->state[ii] += this->dt*(1901*this->rhs[0][ii] - - 2774*this->rhs[1][ii] - + 2616*this->rhs[2][ii] - - 1274*this->rhs[3][ii] - + 251*this->rhs[4][ii])/720; - } - break; - case 6: - for (int p=0; p<this->nparticles; p++) - for (int i=0; i<this->ncomponents; i++) - { - ii = p*this->ncomponents+i; - this->state[ii] += this->dt*(4277*this->rhs[0][ii] - - 7923*this->rhs[1][ii] - + 9982*this->rhs[2][ii] - - 7298*this->rhs[3][ii] - + 2877*this->rhs[4][ii] - - 475*this->rhs[5][ii])/1440; - } - break; - } - this->roll_rhs(); -} - - -template <int particle_type, class rnumber, int interp_neighbours> -void rFFTW_particles<particle_type, rnumber, interp_neighbours>::step() -{ - this->AdamsBashforth((this->iteration < this->integration_steps) ? - this->iteration+1 : - this->integration_steps); - this->iteration++; -} - - -template <int particle_type, class rnumber, int interp_neighbours> -void rFFTW_particles<particle_type, rnumber, interp_neighbours>::read( - const hid_t data_file_id) -{ - if (this->myrank == 0) - { - std::string temp_string = (std::string("/") + - std::string(this->name) + - std::string("/state")); - hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT); - hid_t mspace, rspace; - hsize_t count[4], offset[4]; - rspace = H5Dget_space(dset); - H5Sget_simple_extent_dims(rspace, count, NULL); - count[0] = 1; - offset[0] = this->iteration / this->traj_skip; - offset[1] = 0; - offset[2] = 0; - mspace = H5Screate_simple(3, count, NULL); - H5Sselect_hyperslab(rspace, H5S_SELECT_SET, offset, NULL, count, NULL); - H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, this->state); - H5Sclose(mspace); - H5Sclose(rspace); - H5Dclose(dset); - if (this->iteration > 0) - { - temp_string = (std::string("/") + - std::string(this->name) + - std::string("/rhs")); - dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT); - rspace = H5Dget_space(dset); - H5Sget_simple_extent_dims(rspace, count, NULL); - //reading from last available position - offset[0] = count[0] - 1; - offset[3] = 0; - count[0] = 1; - count[1] = 1; - mspace = H5Screate_simple(4, count, NULL); - for (int i=0; i<this->integration_steps; i++) - { - offset[1] = i; - H5Sselect_hyperslab(rspace, H5S_SELECT_SET, offset, NULL, count, NULL); - H5Dread(dset, H5T_NATIVE_DOUBLE, mspace, rspace, H5P_DEFAULT, this->rhs[i]); - } - H5Sclose(mspace); - H5Sclose(rspace); - H5Dclose(dset); - } - } - MPI_Bcast( - this->state, - this->array_size, - MPI_DOUBLE, - 0, - this->comm); - for (int i = 0; i<this->integration_steps; i++) - MPI_Bcast( - this->rhs[i], - this->array_size, - MPI_DOUBLE, - 0, - this->comm); -} - -template <int particle_type, class rnumber, int interp_neighbours> -void rFFTW_particles<particle_type, rnumber, interp_neighbours>::write( - const hid_t data_file_id, - const char *dset_name, - const double *data) -{ - std::string temp_string = (std::string(this->name) + - std::string("/") + - std::string(dset_name)); - hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT); - hid_t mspace, wspace; - hsize_t count[3], offset[3]; - wspace = H5Dget_space(dset); - H5Sget_simple_extent_dims(wspace, count, NULL); - count[0] = 1; - offset[0] = this->iteration / this->traj_skip; - offset[1] = 0; - offset[2] = 0; - mspace = H5Screate_simple(3, count, NULL); - H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL); - H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, data); - H5Sclose(mspace); - H5Sclose(wspace); - H5Dclose(dset); -} - -template <int particle_type, class rnumber, int interp_neighbours> -void rFFTW_particles<particle_type, rnumber, interp_neighbours>::write( - const hid_t data_file_id, - const bool write_rhs) -{ - if (this->myrank == 0) - { - this->write(data_file_id, "state", this->state); - if (write_rhs) - { - std::string temp_string = ( - std::string("/") + - std::string(this->name) + - std::string("/rhs")); - hid_t dset = H5Dopen(data_file_id, temp_string.c_str(), H5P_DEFAULT); - hid_t wspace = H5Dget_space(dset); - hsize_t count[4], offset[4]; - H5Sget_simple_extent_dims(wspace, count, NULL); - //writing to last available position - offset[0] = count[0] - 1; - offset[1] = 0; - offset[2] = 0; - offset[3] = 0; - count[0] = 1; - count[1] = 1; - hid_t mspace = H5Screate_simple(4, count, NULL); - for (int i=0; i<this->integration_steps; i++) - { - offset[1] = i; - H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL); - H5Dwrite(dset, H5T_NATIVE_DOUBLE, mspace, wspace, H5P_DEFAULT, this->rhs[i]); - } - H5Sclose(mspace); - H5Sclose(wspace); - H5Dclose(dset); - } - } -} - - -/*****************************************************************************/ -template class rFFTW_particles<VELOCITY_TRACER, float, 1>; -template class rFFTW_particles<VELOCITY_TRACER, float, 2>; -template class rFFTW_particles<VELOCITY_TRACER, float, 3>; -template class rFFTW_particles<VELOCITY_TRACER, float, 4>; -template class rFFTW_particles<VELOCITY_TRACER, float, 5>; -template class rFFTW_particles<VELOCITY_TRACER, float, 6>; -template class rFFTW_particles<VELOCITY_TRACER, double, 1>; -template class rFFTW_particles<VELOCITY_TRACER, double, 2>; -template class rFFTW_particles<VELOCITY_TRACER, double, 3>; -template class rFFTW_particles<VELOCITY_TRACER, double, 4>; -template class rFFTW_particles<VELOCITY_TRACER, double, 5>; -template class rFFTW_particles<VELOCITY_TRACER, double, 6>; -/*****************************************************************************/ diff --git a/bfps/cpp/rFFTW_particles.hpp b/bfps/cpp/rFFTW_particles.hpp deleted file mode 100644 index 3ef20c5c3c4f7ce45e34005ef8eaf263017c4d32..0000000000000000000000000000000000000000 --- a/bfps/cpp/rFFTW_particles.hpp +++ /dev/null @@ -1,119 +0,0 @@ -/********************************************************************** -* * -* Copyright 2015 Max Planck Institute * -* for Dynamics and Self-Organization * -* * -* This file is part of bfps. * -* * -* bfps is free software: you can redistribute it and/or modify * -* it under the terms of the GNU General Public License as published * -* by the Free Software Foundation, either version 3 of the License, * -* or (at your option) any later version. * -* * -* bfps is distributed in the hope that it will be useful, * -* but WITHOUT ANY WARRANTY; without even the implied warranty of * -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * -* GNU General Public License for more details. * -* * -* You should have received a copy of the GNU General Public License * -* along with bfps. If not, see <http://www.gnu.org/licenses/> * -* * -* Contact: Cristian.Lalescu@ds.mpg.de * -* * -**********************************************************************/ - - - -#include <stdio.h> -#include <stdlib.h> -#include <iostream> -#include <hdf5.h> -#include "base.hpp" -#include "particles_base.hpp" -#include "fluid_solver_base.hpp" -#include "rFFTW_interpolator.hpp" - -#ifndef RFFTW_PARTICLES - -#define RFFTW_PARTICLES - -template <int particle_type, class rnumber, int interp_neighbours> -class rFFTW_particles -{ - public: - int myrank, nprocs; - MPI_Comm comm; - rnumber *data; - - /* 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. - * */ - double *state; - double *rhs[6]; - int nparticles; - int ncomponents; - int array_size; - int integration_steps; - int traj_skip; - rFFTW_interpolator<rnumber, interp_neighbours> *vel; - - /* simulation parameters */ - char name[256]; - int iteration; - double dt; - - /* methods */ - - /* constructor and destructor. - * allocate and deallocate: - * this->state - * this->rhs - * */ - rFFTW_particles( - const char *NAME, - const hid_t data_file_id, - rFFTW_interpolator<rnumber, interp_neighbours> *FIELD, - const int NPARTICLES, - const int TRAJ_SKIP, - const int INTEGRATION_STEPS = 2); - ~rFFTW_particles(); - - inline const char *get_name() - { - return this->name; - } - - inline void sample(interpolator_base<rnumber, interp_neighbours> *field, double *y) - { - field->sample(this->nparticles, this->ncomponents, this->state, y); - } - - inline void sample( - interpolator_base<rnumber, interp_neighbours> *field, - const hid_t data_file_id, - const char *dset_name) - { - double *y = new double[this->nparticles*this->ncomponents]; - field->sample(this->nparticles, this->ncomponents, this->state, y); - this->write(data_file_id, dset_name, y); - delete[] y; - } - - void get_rhs( - double *__restrict__ x, - double *__restrict__ rhs); - - /* input/output */ - void read(const hid_t data_file_id); - void write(const hid_t data_file_id, const char *dset_name, const double *data); - void write(const hid_t data_file_id, const bool write_rhs = true); - - /* solvers */ - void step(); - void roll_rhs(); - void AdamsBashforth(const int nsteps); -}; - -#endif//RFFTW_PARTICLES - diff --git a/setup.py b/setup.py index 3bc291a2cf0fec3959f1c5b888c6dc007c4be9e0..8f982144e9188dbc74eea1bb4d6ef1f32508b45c 100644 --- a/setup.py +++ b/setup.py @@ -97,7 +97,6 @@ src_file_list = ['field_descriptor', 'interpolator', 'particles', 'rFFTW_interpolator', - 'rFFTW_particles', 'fluid_solver_base', 'fluid_solver', 'fftw_tools', diff --git a/tests/base.py b/tests/base.py index b859eb5fc8521dfd4b8d924c0b10032212013be8..2c616aeb1bbfbdc1b346ae365836429daafa93b0 100644 --- a/tests/base.py +++ b/tests/base.py @@ -116,7 +116,7 @@ def launch( tracer_state_file = None, vorticity_field = None, code_class = bfps.NavierStokes, - particle_class = 'rFFTW_particles', + particle_class = 'particles', interpolator_class = 'rFFTW_interpolator'): c = code_class( work_dir = opt.work_dir, diff --git a/tests/test_plain.py b/tests/test_plain.py index e400b185a81957b8ee476fd4c2ec58274f6ffe4f..e57d6a3b9c82aa197b5f272475a3501fb9474383 100644 --- a/tests/test_plain.py +++ b/tests/test_plain.py @@ -48,8 +48,6 @@ parser.add_argument( def plain(opt): wd = opt.work_dir opt.work_dir = wd + '/N{0:0>3x}_1'.format(opt.n) - if 'rFFTW' in opt.particle_class: - opt.interpolator_class = 'rFFTW_interpolator' c0 = launch(opt, dt = 0.2/opt.n, particle_class = opt.particle_class, interpolator_class = opt.interpolator_class)