diff --git a/CMakeLists.txt b/CMakeLists.txt index 0119e047edc6d5a2f3740217b819f35614b1802a..5b9abaf32ccf87fcf0f87d761340709059a9f96d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -301,6 +301,7 @@ set(cpp_for_lib ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/field_tinterpolator.cpp ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/particle_set.cpp ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_rhs.cpp + ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_with_collision_counter_rhs.cpp ${PROJECT_SOURCE_DIR}/cpp/particles/particle_solver.cpp ${PROJECT_SOURCE_DIR}/cpp/scope_timer.cpp ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp @@ -359,6 +360,7 @@ set(hpp_for_lib ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/field_tinterpolator.hpp ${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particle_rhs.hpp ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_rhs.hpp + ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp ${PROJECT_SOURCE_DIR}/cpp/particles/particle_solver.hpp ${PROJECT_SOURCE_DIR}/cpp/scope_timer.hpp ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.hpp diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp index 74112e870fcfdaaac173f23aafae28b17fe8e8bb..1012652a0be8890945885e11518fcf867ec7d70d 100644 --- a/cpp/full_code/test_tracer_set.cpp +++ b/cpp/full_code/test_tracer_set.cpp @@ -28,6 +28,7 @@ #include <random> #include "test_tracer_set.hpp" #include "particles/rhs/tracer_rhs.hpp" +#include "particles/rhs/tracer_with_collision_counter_rhs.hpp" #include "particles/interpolation/particle_set.hpp" #include "particles/particle_solver.hpp" #include "particles/particles_input_random.hpp" @@ -71,7 +72,7 @@ int test_tracer_set<rnumber>::do_work(void) vec_field->real_space_representation = true; field_tinterpolator<rnumber, FFTW, THREE, NONE> *fti = new field_tinterpolator<rnumber, FFTW, THREE, NONE>(); - tracer_rhs<rnumber, FFTW, NONE> *trhs = new tracer_rhs<rnumber, FFTW, NONE>(); + tracer_with_collision_counter_rhs<rnumber, FFTW, NONE> *trhs = new tracer_with_collision_counter_rhs<rnumber, FFTW, NONE>(); particle_set<3, 2, 1> pset( vec_field->rlayout, diff --git a/cpp/particles/abstract_particle_rhs.hpp b/cpp/particles/abstract_particle_rhs.hpp index 0a5fddc6e684ebcbddb25fc4dd8f5aab12b32978..3687997809c041c16ca2243f2d9f624c3b79400c 100644 --- a/cpp/particles/abstract_particle_rhs.hpp +++ b/cpp/particles/abstract_particle_rhs.hpp @@ -42,6 +42,11 @@ class abstract_particle_rhs double t, abstract_particle_set &pset, particle_rnumber *result) const = 0; + + // for post-timestep actions such as + // normalization of orientation vector + virtual int imposeModelConstraints( + abstract_particle_set &pset) const = 0; }; #endif//ABSTRACT_PARTICLE_RHS_HPP diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp index c452a50a2dd173fc886cd72a980112017114aaac..232ff2c5c92e4a04cd954c57f4849bb9a26b1bf9 100644 --- a/cpp/particles/interpolation/abstract_particle_set.hpp +++ b/cpp/particles/interpolation/abstract_particle_set.hpp @@ -31,6 +31,7 @@ #include <vector> #include <memory> #include "field.hpp" +#include "particles/p2p/p2p_distr_mpi.hpp" @@ -56,8 +57,8 @@ class abstract_particle_set virtual ~abstract_particle_set(){} // extract particle information - virtual const particle_rnumber* getParticleState() const = 0; - virtual const partsize_t* getParticleIndex() const = 0; + virtual particle_rnumber* getParticleState() const = 0; + virtual partsize_t* getParticleIndex() const = 0; virtual int setParticleState(particle_rnumber *) = 0; virtual std::unique_ptr<particle_rnumber[]> extractFromParticleState( @@ -66,8 +67,12 @@ class abstract_particle_set virtual partsize_t getLocalNumberOfParticles() const = 0; virtual partsize_t getTotalNumberOfParticles() const = 0; + virtual partsize_t* getParticlesPerPartition() const = 0; virtual const int getStateSize() const = 0; + // get p2p computer + virtual p2p_distr_mpi<partsize_t, particle_rnumber>* getP2PComputer() = 0; + // generic loop over data template <class func_type> int LOOP(func_type expression) @@ -125,6 +130,28 @@ class abstract_particle_set FFTW, THREExTHREE> &field_to_sample, particle_rnumber *result) = 0; + + // apply p2p kernel + template <int state_size, + class p2p_kernel_class> + int applyP2PKernel( + p2p_kernel_class p2p_kernel, + std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data) + { + // there must be at least one array with additional data, + // since the p2p kernel expects an array where to store the result of the computation + assert(int(additional_data.size()) > int(0)); + this->getP2PComputer()->template compute_distr<p2p_kernel_class, + state_size, // I'd prefer to use this->getStateSize() here. not possible because virtual function is not constexpr... + state_size>( + p2p_kernel, + this->getParticlesPerPartition(), + this->getParticleState(), + additional_data.data(), + int(additional_data.size()), + this->getParticleIndex()); + return EXIT_SUCCESS; + } }; #endif//ABSTRACT_PARTICLE_SET_HPP diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp index 6f474008b5f4bf542cc16d3be70260ad5882a7e4..57fd2cf5300e53870d1b9b1616a779c00fcc163d 100644 --- a/cpp/particles/interpolation/particle_set.hpp +++ b/cpp/particles/interpolation/particle_set.hpp @@ -38,6 +38,7 @@ #include "particles/interpolation/particles_field_computer.hpp" #include "particles/interpolation/particles_generic_interp.hpp" #include "particles/interpolation/abstract_particle_set.hpp" +#include "particles/p2p/p2p_distr_mpi.hpp" @@ -80,22 +81,24 @@ class particle_set: public abstract_particle_set std::vector<hsize_t> particle_file_layout; // raw objects for distributing data in MPI, and for interpolating fields - distributor_class pdistributor; + distributor_class pDistributor; particles_generic_interp<particle_rnumber, neighbours, smoothness> interpolation_formula; - interpolator_class pinterpolator; + interpolator_class pInterpolator; + p2p_distr_mpi<partsize_t, particle_rnumber> p2pComputer; public: particle_set( const abstract_field_layout *fl, const double dkx, const double dky, - const double dkz): + const double dkz, + const double p2p_cutoff = 1.0): mpi_comm(fl->getMPIComm()), current_partition_interval( {fl->getStart(IDXV_Z), fl->getStart(IDXV_Z) + fl->getSubSize(IDXV_Z)}), partition_interval_size(current_partition_interval.second - current_partition_interval.first), - pdistributor( + pDistributor( mpi_comm, current_partition_interval, std::array<size_t, 3>({ //field_grid_dim from particles_system_builder @@ -104,7 +107,7 @@ class particle_set: public abstract_particle_set fl->getSize(IDXV_Z), })), interpolation_formula(), - pinterpolator( + pInterpolator( std::array<size_t, 3>({ //field_grid_dim from particles_system_builder fl->getSize(IDXV_X), fl->getSize(IDXV_Y), @@ -126,7 +129,27 @@ class particle_set: public abstract_particle_set (4*acos(0) / dkx) / fl->getSize(IDXV_X), (4*acos(0) / dky) / fl->getSize(IDXV_Y), (4*acos(0) / dkz) / fl->getSize(IDXV_Z) - })) + })), + p2pComputer( + mpi_comm, + current_partition_interval, + std::array<size_t, 3>({ //field_grid_dim from particles_system_builder + fl->getSize(IDXV_X), + fl->getSize(IDXV_Y), + fl->getSize(IDXV_Z), + }), + std::array<particle_rnumber, 3>({ // spatial_box_width from particles_system_builder + 4*acos(0) / dkx, + 4*acos(0) / dky, + 4*acos(0) / dkz, + }), + std::array<particle_rnumber, 3>({ // spatial_box_offset from particles_system_builder + 0, + 0, + 0, + }), + p2p_cutoff) + { // if these assertions fail, // it means the constructor prototype is broken @@ -145,7 +168,7 @@ class particle_set: public abstract_particle_set ~particle_set(){} - const particle_rnumber* getParticleState() const + particle_rnumber* getParticleState() const { return this->local_state.get(); } @@ -159,7 +182,7 @@ class particle_set: public abstract_particle_set return EXIT_SUCCESS; } - const partsize_t* getParticleIndex() const + partsize_t* getParticleIndex() const { return this->local_index.get(); } @@ -174,6 +197,11 @@ class particle_set: public abstract_particle_set return this->total_number_of_particles; } + partsize_t* getParticlesPerPartition() const + { + return this->number_particles_per_partition.get(); + } + const int getStateSize() const { return state_size; @@ -217,11 +245,11 @@ class particle_set: public abstract_particle_set { // attention: compute_distr adds result on top of existing values. // please clean up result as appropriate before call. - this->pdistributor.template compute_distr<interpolator_class, + this->pDistributor.template compute_distr<interpolator_class, field<field_rnumber, be, fc>, state_size, ncomp(fc)>( - this->pinterpolator, + this->pInterpolator, field_to_sample, this->number_particles_per_partition.get(), this->local_state.get(), @@ -263,11 +291,11 @@ class particle_set: public abstract_particle_set int redistribute(std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data) { - this->pdistributor.template redistribute<interpolator_class, + this->pDistributor.template redistribute<interpolator_class, state_size, state_size, 1>( - this->pinterpolator, + this->pInterpolator, this->number_particles_per_partition.get(), &this->local_number_of_particles, &this->local_state, @@ -292,7 +320,7 @@ class particle_set: public abstract_particle_set this->number_particles_per_partition.get(), this->offset_particles_for_partition.get(), [&](const particle_rnumber& z_pos){ - const int partition_level = this->pinterpolator.pbc_field_layer(z_pos, IDXC_Z); + const int partition_level = this->pInterpolator.pbc_field_layer(z_pos, IDXC_Z); assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second); return partition_level - current_partition_interval.first; }, @@ -305,12 +333,17 @@ class particle_set: public abstract_particle_set particle_rnumber getSpatialLowLimitZ() { - return this->pinterpolator.getSpatialLowLimitZ(); + return this->pInterpolator.getSpatialLowLimitZ(); } particle_rnumber getSpatialUpLimitZ() { - return this->pinterpolator.getSpatialUpLimitZ(); + return this->pInterpolator.getSpatialUpLimitZ(); + } + + p2p_distr_mpi<partsize_t, particle_rnumber> *getP2PComputer() + { + return &(this->p2pComputer); } }; diff --git a/cpp/particles/particle_solver.cpp b/cpp/particles/particle_solver.cpp index e90e16731ce1006ab619dcc5fc61025486f76ff2..a67691a1d877d1225c22e0133cdcc1ce1c1da6fd 100644 --- a/cpp/particles/particle_solver.cpp +++ b/cpp/particles/particle_solver.cpp @@ -51,6 +51,9 @@ int particle_solver::Euler( // deallocate temporary arrays delete[] rhs_val; delete[] xx; + + // ensure new state is consistent with physical model + rhs.imposeModelConstraints(*(this->pset)); return EXIT_SUCCESS; } diff --git a/cpp/particles/rhs/tracer_rhs.hpp b/cpp/particles/rhs/tracer_rhs.hpp index c4f0342870d783911ab6ebd899ccff1fa8786e15..676789ae1ad5ffdcce91d13d2e456fadcce591e0 100644 --- a/cpp/particles/rhs/tracer_rhs.hpp +++ b/cpp/particles/rhs/tracer_rhs.hpp @@ -58,6 +58,12 @@ class tracer_rhs: public abstract_particle_rhs std::fill_n(result, pset.getLocalNumberOfParticles()*3, 0); return (*(this->velocity))(t, pset, result); } + + int imposeModelConstraints( + abstract_particle_set &pset) const + { + return EXIT_SUCCESS; + } }; #endif//TRACER_RHS_HPP diff --git a/cpp/particles/rhs/tracer_with_collision_counter_rhs.cpp b/cpp/particles/rhs/tracer_with_collision_counter_rhs.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3f9efd293aae9da1048c856418a31cf5b2e1f82c --- /dev/null +++ b/cpp/particles/rhs/tracer_with_collision_counter_rhs.cpp @@ -0,0 +1,31 @@ +/****************************************************************************** +* * +* Copyright 2020 Max Planck Institute for Dynamics and Self-Organization * +* * +* This file is part of TurTLE. * +* * +* TurTLE 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. * +* * +* TurTLE 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 TurTLE. If not, see <http://www.gnu.org/licenses/> * +* * +* Contact: Cristian.Lalescu@ds.mpg.de * +* * +******************************************************************************/ + + + +#include "particles/rhs/tracer_with_collision_counter_rhs.hpp" + + +template class tracer_with_collision_counter_rhs<float, FFTW, NONE>; +template class tracer_with_collision_counter_rhs<double, FFTW, NONE>; + diff --git a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp new file mode 100644 index 0000000000000000000000000000000000000000..9e33b3a2a95e25c1219f84a7db581295b21b2fed --- /dev/null +++ b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp @@ -0,0 +1,66 @@ +/****************************************************************************** +* * +* Copyright 2020 Max Planck Institute for Dynamics and Self-Organization * +* * +* This file is part of TurTLE. * +* * +* TurTLE 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. * +* * +* TurTLE 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 TurTLE. If not, see <http://www.gnu.org/licenses/> * +* * +* Contact: Cristian.Lalescu@ds.mpg.de * +* * +******************************************************************************/ + + + +#ifndef TRACER_WITH_COLLISION_COUNTER_RHS_HPP +#define TRACER_WITH_COLLISION_COUNTER_RHS_HPP + +#include "particles/interpolation/field_tinterpolator.hpp" +#include "particles/abstract_particle_rhs.hpp" +#include "particles/rhs/tracer_rhs.hpp" +#include "particles/p2p/p2p_ghost_collisions.hpp" + + +template <typename rnumber, + field_backend be, + temporal_interpolation_type tt> +class tracer_with_collision_counter_rhs: public tracer_rhs<rnumber, be, tt> +{ + protected: + p2p_ghost_collisions<abstract_particle_set::particle_rnumber, abstract_particle_set::partsize_t> p2p_gc; + + public: + tracer_with_collision_counter_rhs(){} + ~tracer_with_collision_counter_rhs(){} + + int operator()( + double t, + abstract_particle_set &pset, + abstract_particle_set::particle_rnumber *result) const + { + // interpolation adds on top of existing values, so result must be cleared. + std::fill_n(result, pset.getLocalNumberOfParticles()*3, 0); + int interpolation_result = (*(this->velocity))(t, pset, result); + 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.template applyP2PKernel<3, p2p_ghost_collisions<abstract_particle_set::particle_rnumber, abstract_particle_set::partsize_t>>( + this->p2p_gc, + bla); + return interpolation_result; + } +}; + +#endif//TRACER_WITH_COLLISION_COUNTER_RHS_HPP +