[bug] add test for p2p data shuffling

TurTLE/test/particle_set/ADNS.py

+import os
+import sys
+import argparse
+import getpass
+
+import numpy as np
+import h5py
+
+import TurTLE
+import TurTLE._base
+import TurTLE.DNS
+from TurTLE._code import _code
+
+cpp_location = os.path.join(os.path.join(
+        TurTLE.lib_dir, 'test'), 'particle_set')
+
+
+class ADNS(TurTLE.DNS):
+    def __init__(
+            self,
+            **kwargs):
+        TurTLE.DNS.__init__(self, **kwargs)
+        return None
+    def write_src(
+            self):
+        self.version_message = (
+                '/***********************************************************************\n' +
+                '* this code automatically generated by TurTLE\n' +
+                '* version {0}\n'.format(TurTLE.__version__) +
+                '***********************************************************************/\n\n\n')
+        self.include_list = [
+                '"base.hpp"',
+                '"scope_timer.hpp"',
+                '"fftw_interface.hpp"',
+                '"full_code/main_code.hpp"',
+                '"full_code/NSVEparticles.hpp"',
+                '<cmath>',
+                '<iostream>',
+                '<hdf5.h>',
+                '<string>',
+                '<cstring>',
+                '<fftw3-mpi.h>',
+                '<omp.h>',
+                '<cfenv>',
+                '<cstdlib>']
+        self.main = """
+            int main(int argc, char *argv[])
+            {{
+                bool fpe = (
+                    (getenv("BFPS_FPE_OFF") == nullptr) ||
+                    (getenv("BFPS_FPE_OFF") != std::string("TRUE")));
+                return main_code< {0} >(argc, argv, fpe);
+            }}
+            """.format(self.dns_type + '<{0}>'.format(self.C_field_dtype))
+        self.includes = '\n'.join(
+                ['#include ' + hh
+                 for hh in self.include_list])
+        self.name = 'NSVEparticle_set'
+        self.dns_type = 'NSVEparticle_set'
+        self.definitions += open(
+            os.path.join(
+                cpp_location, 'NSVEparticle_set.hpp'), 'r').read()
+        self.definitions += open(
+            os.path.join(
+                cpp_location, 'NSVEparticle_set.cpp'), 'r').read()
+        with open(self.name + '.cpp', 'w') as outfile:
+            outfile.write(self.version_message + '\n\n')
+            outfile.write(self.includes + '\n\n')
+            outfile.write(self.definitions + '\n\n')
+            outfile.write(self.main + '\n')
+        self.check_current_vorticity_exists = True
+        return None
+    def generate_default_parameters(self):
+        # these parameters are relevant for all DNS classes
+        self.parameters['fftw_plan_rigor'] = 'FFTW_ESTIMATE'
+        self.parameters['dealias_type'] = int(1)
+        self.parameters['dkx'] = float(1.0)
+        self.parameters['dky'] = float(1.0)
+        self.parameters['dkz'] = float(1.0)
+        self.parameters['niter_todo'] = int(200)
+        self.parameters['niter_stat'] = int(20)
+        self.parameters['niter_out'] = int(200)
+        self.parameters['checkpoints_per_file'] = int(1)
+        self.parameters['dt'] = float(0.001)
+        self.parameters['histogram_bins'] = int(64)
+        self.parameters['max_velocity_estimate'] = float(1)
+        self.parameters['max_vorticity_estimate'] = float(1)
+        self.parameters['niter_part'] = int(1)
+        self.parameters['niter_part_fine_period'] = int(2)
+        self.parameters['niter_part_fine_duration'] = int(0)
+        self.parameters['nparticles'] = int(1000)
+        self.parameters['tracers0_integration_steps'] = int(0)
+        self.parameters['tracers0_neighbours'] = int(3)
+        self.parameters['tracers0_smoothness'] = int(2)
+        self.parameters['tracers0_enable_p2p'] = int(0)
+        self.parameters['tracers0_enable_inner'] = int(0)
+        self.parameters['tracers0_enable_vorticity_omega'] = int(0)
+        self.parameters['tracers0_cutoff'] = float(0.314)
+        self.parameters['tracers0_inner_v0'] = float(1)
+        self.parameters['tracers0_lambda'] = float(1)
+
+        self.parameters['nu'] = float(0.1)
+        self.parameters['fmode'] = int(1)
+        self.parameters['famplitude'] = float(0.5)
+        self.parameters['friction_coefficient'] = float(0.5)
+        self.parameters['injection_rate'] = float(0.4)
+        self.parameters['fk0'] = float(2.0)
+        self.parameters['fk1'] = float(4.0)
+        self.parameters['energy'] = float(0.5)
+        self.parameters['forcing_type'] = 'fixed_energy_injection_rate'
+        return None
+    def generate_tracer_state(
+            self,
+            rseed = None,
+            species = 0,
+            integration_steps = None,
+            ncomponents = 3):
+        try:
+            if type(integration_steps) == type(None):
+                integration_steps = self.NSVEp_extra_parameters['tracers0_integration_steps']
+            if 'tracers{0}_integration_steps'.format(species) in self.parameters.keys():
+                integration_steps = self.parameters['tracers{0}_integration_steps'.format(species)]
+            with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file:
+                nn = self.parameters['nparticles']
+                if not 'tracers{0}'.format(species) in data_file.keys():
+                    data_file.create_group('tracers{0}'.format(species))
+                    data_file.create_group('tracers{0}/rhs'.format(species))
+                    data_file.create_group('tracers{0}/state'.format(species))
+                data_file['tracers{0}/rhs'.format(species)].create_dataset(
+                        '0',
+                        shape = (integration_steps, nn, ncomponents,),
+                        dtype = np.float)
+                dset = data_file['tracers{0}/state'.format(species)].create_dataset(
+                        '0',
+                        shape = (nn, ncomponents,),
+                        dtype = np.float)
+                if not type(rseed) == type(None):
+                    np.random.seed(rseed)
+                cc = int(0)
+                batch_size = int(1e6)
+                def get_random_phases(npoints):
+                    return np.random.random(
+                                (npoints, 3))*2*np.pi
+                def get_random_versors(npoints):
+                    bla = np.random.normal(
+                            size = (npoints, 3))
+                    bla  /= np.sum(bla**2, axis = 1)[:, None]**.5
+                    return bla
+                while nn > 0:
+                    #TODO introduce initial vel parameter as in established Stokes class, now set to zero
+                    if nn > batch_size:
+                        dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size)
+                        if dset.shape[1] == 6:
+                            dset[cc*batch_size:(cc+1)*batch_size, 3:] = 0.0* get_random_versors(batch_size)
+                        nn -= batch_size
+                    else:
+                        dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
+                        if dset.shape[1] == 6:
+                            dset[cc*batch_size:cc*batch_size+nn, 3:] = 0.0* get_random_versors(nn)
+                        nn = 0
+                    cc += 1
+        except Exception as e:
+            print(e)
+        return None
+    def prepare_launch(
+            self,
+            args = [],
+            **kwargs):
+        parser = argparse.ArgumentParser('turtle ' + type(self).__name__)
+        self.job_parser_arguments(parser)
+        self.simulation_parser_arguments(parser)
+        self.particle_parser_arguments(parser)
+        self.parameters_to_parser_arguments(parser)
+        opt = parser.parse_args(args)
+        self.simname=opt.simname
+        self.set_precision(opt.precision)
+        self.dns_type = 'NSVEparticle_set'
+        self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
+        if ((self.parameters['niter_todo'] % self.parameters['niter_out']) != 0):
+            self.parameters['niter_out'] = self.parameters['niter_todo']
+        if len(opt.src_work_dir) == 0:
+            opt.src_work_dir = os.path.realpath(opt.work_dir)
+        if type(opt.dkx) == type(None):
+            opt.dkx = 2. / opt.Lx
+        if type(opt.dky) == type(None):
+            opt.dky = 2. / opt.Ly
+        if type(opt.dkz) == type(None):
+            opt.dkz = 2. / opt.Lz
+        if type(opt.nx) == type(None):
+            opt.nx = opt.n
+        if type(opt.ny) == type(None):
+            opt.ny = opt.n
+        if type(opt.nz) == type(None):
+            opt.nz = opt.n
+        if type(opt.fk0) == type(None):
+            opt.fk0 = self.parameters['fk0']
+        if type(opt.fk1) == type(None):
+            opt.fk1 = self.parameters['fk1']
+        if type(opt.injection_rate) == type(None):
+            opt.injection_rate = self.parameters['injection_rate']
+        if type(opt.dealias_type) == type(None):
+            opt.dealias_type = self.parameters['dealias_type']
+        if (opt.nx > opt.n or
+            opt.ny > opt.n or
+            opt.nz > opt.n):
+            opt.n = min(opt.nx, opt.ny, opt.nz)
+            print("Warning: '-n' parameter changed to minimum of nx, ny, nz. This affects the computation of nu.")
+        self.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
+        # check value of kMax
+        kM = opt.n * 0.5
+        if opt.dealias_type == 1:
+            kM *= 0.8
+        # tweak forcing/viscosity based on forcint type
+        if opt.forcing_type == 'linear':
+            # custom famplitude for 288 and 576
+            if opt.n == 288:
+                self.parameters['famplitude'] = 0.45
+            elif opt.n == 576:
+                self.parameters['famplitude'] = 0.47
+        elif opt.forcing_type == 'fixed_energy_injection_rate':
+            # use the fact that mean dissipation rate is equal to injection rate
+            self.parameters['nu'] = (
+                    opt.injection_rate *
+                    (opt.kMeta / kM)**4)**(1./3)
+        elif opt.forcing_type == 'fixed_energy':
+            kf = 1. / (1./opt.fk0 +
+                       1./opt.fk1)
+            self.parameters['nu'] = (
+                    (opt.kMeta / kM)**(4./3) *
+                    (np.pi / kf)**(1./3) *
+                    (2*self.parameters['energy'] / 3)**0.5)
+        if type(opt.checkpoints_per_file) == type(None):
+            # hardcoded FFTW complex representation size
+            field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize
+            checkpoint_size = field_size
+        self.set_host_info(TurTLE.host_info)
+        if type(opt.environment) != type(None):
+            self.host_info['environment'] = opt.environment
+        self.pars_from_namespace(opt)
+        return opt
+
+    def write_par(
+            self,
+            iter0 = 0):
+        assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0)
+        assert (self.parameters['niter_todo'] % self.parameters['niter_out']  == 0)
+        assert (self.parameters['niter_out']  % self.parameters['niter_stat'] == 0)
+        assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
+        assert (self.parameters['niter_out']  % self.parameters['niter_part'] == 0)
+        _code.write_par(self, iter0 = iter0)
+        with h5py.File(self.get_data_file_name(), 'r+') as ofile:
+            ofile['code_info/exec_name'] = self.name
+            kspace = self.get_kspace()
+            for k in kspace.keys():
+                ofile['kspace/' + k] = kspace[k]
+            nshells = kspace['nshell'].shape[0]
+            kspace = self.get_kspace()
+            nshells = kspace['nshell'].shape[0]
+            vec_stat_datasets = ['velocity', 'vorticity']
+            scal_stat_datasets = []
+            for k in vec_stat_datasets:
+                time_chunk = 2**20//(8*3*3*nshells)
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/spectra/' + k + '_' + k,
+                                     (1, nshells, 3, 3),
+                                     chunks = (time_chunk, nshells, 3, 3),
+                                     maxshape = (None, nshells, 3, 3),
+                                     dtype = np.float64)
+                time_chunk = 2**20//(8*4*10)
+                time_chunk = max(time_chunk, 1)
+                a = ofile.create_dataset('statistics/moments/' + k,
+                                     (1, 10, 4),
+                                     chunks = (time_chunk, 10, 4),
+                                     maxshape = (None, 10, 4),
+                                     dtype = np.float64)
+                time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/histograms/' + k,
+                                     (1,
+                                      self.parameters['histogram_bins'],
+                                      4),
+                                     chunks = (time_chunk,
+                                               self.parameters['histogram_bins'],
+                                               4),
+                                     maxshape = (None,
+                                                 self.parameters['histogram_bins'],
+                                                 4),
+                                     dtype = np.int64)
+            #TODO move collision stats to particle stats
+            for i in range(self.parameters['niter_out']):
+                ofile.create_dataset('statistics/collisions/'+str(i),
+                                     shape=(1,),
+                                     dtype = np.int64)
+            ofile['checkpoint'] = int(0)
+        if (self.dns_type in ['NSVE', 'NSVE_no_output']):
+            return None
+        return None
+
+    def launch(
+            self,
+            args = [],
+            **kwargs):
+        opt = self.prepare_launch(args = args)
+        if not os.path.exists(self.get_data_file_name()):
+            self.write_par()
+            self.generate_initial_condition(opt = opt)
+            with h5py.File(self.get_particle_file_name(), 'w') as particle_file:
+                particle_file.create_group('tracers0/position')
+                particle_file.create_group('tracers0/velocity')
+                particle_file.create_group('collisions0')
+            self.generate_tracer_state(integration_steps = self.parameters['tracers0_integration_steps'],
+                                       rseed = opt.particle_rand_seed,
+                                       ncomponents = 3)
+        self.run(
+                nb_processes = opt.nb_processes,
+                nb_threads_per_process = opt.nb_threads_per_process,
+                njobs = opt.njobs,
+                hours = opt.minutes // 60,
+                minutes = opt.minutes % 60,
+                no_submit = opt.no_submit,
+                no_debug = opt.no_debug)
+        return None
+
+def main():
+    bla = ADNS()
+    bla.launch(sys.argv[1:])
+    df = bla.get_particle_file()
+    tt = 0
+    x = []
+    for ii in range(200):
+        x.append(df['tracers0/position/{0}'.format(ii)][tt])
+    x = np.array(x)
+    df.close()
+    try:
+        import matplotlib.pyplot as plt
+        f = plt.figure()
+        a = f.add_subplot(111)
+        a.plot(x)
+        f.tight_layout()
+        f.savefig('traj.pdf')
+    except ImportError:
+        print('couldn\'t import matplotlib, no figure')
+    return None
+
+if __name__ == '__main__':
+    main()

TurTLE/test/particle_set/NSVEparticle_set.cpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2019 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 <string>
+#include <cmath>
+//#include "NSVEparticle_set.hpp"
+#include "scope_timer.hpp"
+
+template <typename rnumber>
+int NSVEparticle_set<rnumber>::initialize(void)
+{
+    TIMEZONE("NSVEparticle_set::intialize");
+    this->NSVE<rnumber>::initialize();
+
+
+    // allocate previous velocity
+    this->previous_velocity = new field<rnumber, FFTW, THREE>(
+            this->fs->cvelocity->rlayout->sizes[2],
+            this->fs->cvelocity->rlayout->sizes[1],
+            this->fs->cvelocity->rlayout->sizes[0],
+            this->fs->cvelocity->rlayout->comm,
+            this->fs->cvelocity->fftw_plan_rigor);
+    *this->previous_velocity = *this->fs->cvelocity;
+    this->previous_velocity->ift();
+
+    this->fti = new field_tinterpolator<rnumber, FFTW, THREE, LINEAR>();
+    this->trhs = new tracer_with_collision_counter_rhs<rnumber, FFTW, LINEAR>();
+
+    // We're not using Adams-Bashforth in this code.
+    // This assert is there to ensure
+    // user knows what's happening.
+    assert(tracers0_integration_steps == 0);
+    // neighbours and smoothness are second and third template parameters of particle_set
+    assert(tracers0_neighbours == 3);
+    assert(tracers0_smoothness == 2);
+    this->pset = new particle_set<3, 3, 2>(
+            this->fs->cvelocity->rlayout,
+            this->dkx,
+            this->dky,
+            this->dkz,
+            0.5);
+
+    // set two fields for the temporal interpolation
+    this->fti->set_field(this->previous_velocity, 0);
+    this->fti->set_field(this->fs->cvelocity, 1);
+    this->trhs->setVelocity(this->fti);
+
+    particles_input_hdf5<long long int, double, 3, 3> particle_reader(
+            this->comm,
+            this->fs->get_current_fname(),
+            std::string("/tracers0/state/") + std::to_string(this->fs->iteration), // dataset name for initial input
+            std::string("/tracers0/rhs/")  + std::to_string(this->fs->iteration),  // dataset name for initial input
+            this->pset->getSpatialLowLimitZ(),
+            this->pset->getSpatialUpLimitZ());
+    this->pset->init(particle_reader);
+
+    this->psolver = new particle_solver(*(this->pset), 0);
+
+    this->particles_output_writer_mpi = new particles_output_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                "tracers0",
+                nparticles,
+                tracers0_integration_steps);
+    this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                this->pset->getTotalNumberOfParticles(),
+                (this->simname + "_particles.h5"),
+                "tracers0",
+                "position/0");
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEparticle_set<rnumber>::step(void)
+{
+    TIMEZONE("NSVEparticle_set::step");
+    // compute next step Navier-Stokes
+    this->NSVE<rnumber>::step();
+    // update velocity 1
+    this->fs->compute_velocity(this->fs->cvorticity);
+    this->fs->cvelocity->ift();
+    // compute particle step using velocity 0 and velocity 1
+    this->psolver->Heun(this->dt, *(this->trhs));
+    // update velocity 0
+    *this->previous_velocity = *this->fs->cvelocity;
+    this->psolver->setIteration(this->fs->iteration);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEparticle_set<rnumber>::write_checkpoint(void)
+{
+    TIMEZONE("NSVEparticle_set::write_checkpoint");
+    this->NSVE<rnumber>::write_checkpoint();
+    this->psolver->setIteration(this->fs->iteration);
+    this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
+    this->psolver->template writeCheckpoint<3>(this->particles_output_writer_mpi);
+    this->particles_output_writer_mpi->close_file();
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEparticle_set<rnumber>::finalize(void)
+{
+    TIMEZONE("NSVEparticle_set::finalize");
+    delete this->particles_output_writer_mpi;
+    delete this->particles_sample_writer_mpi;
+    delete this->psolver;
+    delete this->pset;
+    delete this->trhs;
+    delete this->fti;
+    this->NSVE<rnumber>::finalize();
+    return EXIT_SUCCESS;
+}
+
+/** \brief Compute fluid stats and sample fields at particle locations.
+ */
+
+template <typename rnumber>
+int NSVEparticle_set<rnumber>::do_stats()
+{
+    TIMEZONE("NSVEparticle_set::do_stats");
+    /// fluid stats go here
+    this->NSVE<rnumber>::do_stats();
+
+
+    /// either one of two conditions suffices to compute statistics:
+    /// 1) current iteration is a multiple of niter_part
+    /// 2) we are within niter_part_fine_duration/2 of a multiple of niter_part_fine_period
+    if (!(this->iteration % this->niter_part == 0 ||
+          ((this->iteration + this->niter_part_fine_duration/2) % this->niter_part_fine_period <=
+           this->niter_part_fine_duration)))
+        return EXIT_SUCCESS;
+
+    // sample position
+    this->pset->writeStateTriplet(
+            0,
+            this->particles_sample_writer_mpi,
+            "tracers0",
+            "position",
+            psolver->getIteration());
+
+    // sample velocity
+    // first ensure velocity field is computed, and is in real space
+    if (!(this->iteration % this->niter_stat == 0))
+    {
+        // we need to compute velocity field manually, because it didn't happen in NSVE::do_stats()
+        this->fs->compute_velocity(this->fs->cvorticity);
+        *this->tmp_vec_field = this->fs->cvelocity->get_cdata();
+        this->tmp_vec_field->ift();
+    }
+    this->pset->writeSample(
+            this->tmp_vec_field,
+            this->particles_sample_writer_mpi,
+            "tracers0",
+            "velocity",
+            psolver->getIteration());
+
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVEparticle_set<rnumber>::read_parameters(void)
+{
+    TIMEZONE("NSVEparticle_set::read_parameters");
+    this->NSVE<rnumber>::read_parameters();
+    hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+    this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+    this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
+    this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
+    this->nparticles = hdf5_tools::read_value<int>(parameter_file, "parameters/nparticles");
+    this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
+    this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
+    this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
+    H5Fclose(parameter_file);
+    return EXIT_SUCCESS;
+}
+
+template class NSVEparticle_set<float>;
+template class NSVEparticle_set<double>;
+

TurTLE/test/particle_set/NSVEparticle_set.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2017 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                                 *
+*                                                                     *
+**********************************************************************/
+
+
+
+#ifndef NSVEPARTICLE_SET_HPP
+#define NSVEPARTICLE_SET_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/interpolation/field_tinterpolator.hpp"
+#include "particles/interpolation/particle_set.hpp"
+#include "particles/particle_solver.hpp"
+#include "particles/rhs/tracer_with_collision_counter_rhs.hpp"
+
+/** \brief Test of particle set class hierarchy.
+ *
+ *  Child of Navier Stokes vorticity equation solver, this class calls all the
+ *  methods from `NSVE`, and in addition integrates simple Lagrangian tracers
+ *  in the resulting velocity field.
+ *  It also executes a p2p kernel that doesn't affect the particle trajectories
+ *  in any way.
+ *
+ *  The solver is meant to test that the particle set class hierarchy is working
+ *  as desired.
+ */
+
+template <typename rnumber>
+class NSVEparticle_set: public NSVE<rnumber>
+{
+    public:
+
+        /* parameters that are read in read_parameters */
+        int niter_part;
+        int niter_part_fine_period;
+        int niter_part_fine_duration;
+        int nparticles;
+        int tracers0_integration_steps;
+        int tracers0_neighbours;
+        int tracers0_smoothness;
+
+        /* other stuff */
+        field_tinterpolator<rnumber, FFTW, THREE, LINEAR> *fti;
+        tracer_with_collision_counter_rhs<rnumber, FFTW, LINEAR> *trhs;
+        particle_set<3, 3, 2> *pset;
+        particle_solver *psolver;
+
+        field<rnumber, FFTW, THREE> *previous_velocity;
+
+        particles_output_hdf5<long long int, double, 3> *particles_output_writer_mpi;
+        particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+
+        NSVEparticle_set(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            NSVE<rnumber>(
+                    COMMUNICATOR,
+                    simulation_name){}
+        ~NSVEparticle_set(){}
+
+        int initialize(void);
+        int step(void);
+        int finalize(void);
+
+        int read_parameters(void);
+        int write_checkpoint(void);
+        int do_stats(void);
+};
+
+#endif//NSVEPARTICLE_SET_HPP
+

cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp

@@ -55,9 +55,21 @@ class tracer_with_collision_counter_rhs: public tracer_rhs<rnumber, be, tt>
             std::vector<std::unique_ptr<abstract_particle_set::particle_rnumber[]>> bla;
             bla.resize(1);
             bla[0].reset(new abstract_particle_set::particle_rnumber[pset.getLocalNumberOfParticles()*pset.getStateSize()]);
+            // copy rhs values to temporary array
+            pset.LOOP(
+                    [&](const abstract_particle_set::partsize_t idx)
+                    {
+                        bla[0][idx] = result[idx];
+                    });
             pset.template applyP2PKernel<3, p2p_ghost_collisions<abstract_particle_set::particle_rnumber, abstract_particle_set::partsize_t>>(
                     this->p2p_gc,
                     bla);
+            // copy shuffled rhs values
+            pset.LOOP(
+                    [&](const abstract_particle_set::partsize_t idx)
+                    {
+                        result[idx] = bla[0][idx];
+                    });
             return interpolation_result;
         }
 

setup.py

@@ -77,7 +77,11 @@ setup(
         name = 'TurTLE',
         packages = ['TurTLE', 'TurTLE/test'],
         install_requires = ['numpy>=1.8', 'h5py>=2.2.1'],
-        package_data = {'TurTLE': ['test/B32p1e4_checkpoint_0.h5']},
+        package_data = {'TurTLE': ['test/B32p1e4_checkpoint_0.h5',
+                                   'test/particle_set/NSVEparticle_set.hpp',
+                                   'test/particle_set/NSVEparticle_set.cpp',
+                                   'test/particle_set/ADNS.py',
+                                   ]},
         entry_points = {
             'console_scripts': [
                 'turtle = TurTLE.__main__:main',