Commit 7d3bb095 authored by Berenger Bramas's avatar Berenger Bramas
Browse files

Update the kernels and make it possible to update parameters at runtime

parent 359ef9ab
Pipeline #19548 failed with stage
in 2 minutes and 25 seconds
...@@ -58,6 +58,10 @@ class field ...@@ -58,6 +58,10 @@ class field
private: private:
rnumber *__restrict__ data; /**< data array */ rnumber *__restrict__ data; /**< data array */
public: public:
constexpr int nb_components() {
return ncomp(fc);
}
hsize_t npoints; /**< total number of grid points. Useful for normalization. */ hsize_t npoints; /**< total number of grid points. Useful for normalization. */
bool real_space_representation; /**< `true` if field is in real space representation. */ bool real_space_representation; /**< `true` if field is in real space representation. */
......
...@@ -4,13 +4,20 @@ ...@@ -4,13 +4,20 @@
#include "scope_timer.hpp" #include "scope_timer.hpp"
#include "particles/particles_sampling.hpp" #include "particles/particles_sampling.hpp"
#include "particles/p2p_computer.hpp" #include "particles/p2p_computer.hpp"
#include "particles/particles_inner_computer.hpp"
template <typename rnumber> template <typename rnumber>
int NSVEparticlesP2P<rnumber>::initialize(void) int NSVEparticlesP2P<rnumber>::initialize(void)
{ {
this->NSVE<rnumber>::initialize(); this->NSVE<rnumber>::initialize();
this->ps = particles_system_builder( p2p_computer<rnumber, long long int> current_p2p_computer;
current_p2p_computer.setEnable(enable_p2p);
particles_inner_computer<rnumber, long long int> current_particles_inner_computer(inner_v0);
current_particles_inner_computer.setEnable(enable_inner);
this->ps = particles_system_builder_with_p2p(
this->fs->cvelocity, // (field object) this->fs->cvelocity, // (field object)
this->fs->kk, // (kspace object, contains dkx, dky, dkz) this->fs->kk, // (kspace object, contains dkx, dky, dkz)
tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs) tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
...@@ -21,10 +28,13 @@ int NSVEparticlesP2P<rnumber>::initialize(void) ...@@ -21,10 +28,13 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
tracers0_neighbours, // parameter (interpolation no neighbours) tracers0_neighbours, // parameter (interpolation no neighbours)
tracers0_smoothness, // parameter tracers0_smoothness, // parameter
this->comm, this->comm,
this->fs->iteration+1); this->fs->iteration+1,
// TODO P2P write particle data too std::move(current_p2p_computer),
std::move(current_particles_inner_computer),
cutoff);
this->particles_output_writer_mpi = new particles_output_hdf5< this->particles_output_writer_mpi = new particles_output_hdf5<
long long int, double, 3, 3>( long long int, double, 6, 6>(
MPI_COMM_WORLD, MPI_COMM_WORLD,
"tracers0", "tracers0",
nparticles, nparticles,
...@@ -37,7 +47,12 @@ int NSVEparticlesP2P<rnumber>::step(void) ...@@ -37,7 +47,12 @@ int NSVEparticlesP2P<rnumber>::step(void)
{ {
this->fs->compute_velocity(this->fs->cvorticity); this->fs->compute_velocity(this->fs->cvorticity);
this->fs->cvelocity->ift(); this->fs->cvelocity->ift();
this->ps->completeLoop(this->dt); if(enable_vorticity_omega){
this->ps->completeLoopWithVorticity(this->dt, *this->fs->cvorticity);
}
else{
this->ps->completeLoop(this->dt);
}
this->NSVE<rnumber>::step(); this->NSVE<rnumber>::step();
return EXIT_SUCCESS; return EXIT_SUCCESS;
} }
......
...@@ -55,10 +55,15 @@ class NSVEparticlesP2P: public NSVE<rnumber> ...@@ -55,10 +55,15 @@ class NSVEparticlesP2P: public NSVE<rnumber>
int tracers0_neighbours; int tracers0_neighbours;
int tracers0_smoothness; int tracers0_smoothness;
double cutoff;
bool enable_p2p;
bool enable_inner;
bool enable_vorticity_omega;
/* other stuff */ /* other stuff */
std::unique_ptr<abstract_particles_system<long long int, double>> ps; std::unique_ptr<abstract_particles_system<long long int, double>> ps;
// TODO P2P use a reader with particle data // TODO P2P use a reader with particle data
particles_output_hdf5<long long int, double,3,3> *particles_output_writer_mpi; particles_output_hdf5<long long int, double,6,6> *particles_output_writer_mpi;
NSVEparticlesP2P( NSVEparticlesP2P(
...@@ -66,7 +71,8 @@ class NSVEparticlesP2P: public NSVE<rnumber> ...@@ -66,7 +71,8 @@ class NSVEparticlesP2P: public NSVE<rnumber>
const std::string &simulation_name): const std::string &simulation_name):
NSVE<rnumber>( NSVE<rnumber>(
COMMUNICATOR, COMMUNICATOR,
simulation_name){} simulation_name),
cutoff(std::numeric_limits<double>::max()){}
~NSVEparticlesP2P(){} ~NSVEparticlesP2P(){}
int initialize(void); int initialize(void);
......
...@@ -16,6 +16,10 @@ public: ...@@ -16,6 +16,10 @@ public:
virtual void compute_p2p() = 0; virtual void compute_p2p() = 0;
virtual void compute_particles_inner() = 0;
virtual void compute_particles_inner(const real_number particle_extra_rhs[]) = 0;
virtual void move(const real_number dt) = 0; virtual void move(const real_number dt) = 0;
virtual void redistribute() = 0; virtual void redistribute() = 0;
...@@ -26,6 +30,9 @@ public: ...@@ -26,6 +30,9 @@ public:
virtual void completeLoop(const real_number dt) = 0; virtual void completeLoop(const real_number dt) = 0;
virtual void completeLoopWithVorticity(const real_number dt,
const real_number particle_extra_rhs[]) = 0;
virtual const real_number* getParticlesPositions() const = 0; virtual const real_number* getParticlesPositions() const = 0;
virtual const std::unique_ptr<real_number[]>* getParticlesRhs() const = 0; virtual const std::unique_ptr<real_number[]>* getParticlesRhs() const = 0;
...@@ -54,6 +61,15 @@ public: ...@@ -54,6 +61,15 @@ public:
virtual void sample_compute_field(const field<double, FFTW, THREExTHREE>& sample_field, virtual void sample_compute_field(const field<double, FFTW, THREExTHREE>& sample_field,
real_number sample_rhs[]) = 0; real_number sample_rhs[]) = 0;
//- Not generic to enable sampling end //- Not generic to enable sampling end
template <typename rnumber, field_backend be, field_components fc>
void completeLoopWithVorticity(const real_number dt,
const field<rnumber, be, fc>& in_field) {
static_assert(fc == THREE, "only THREE is supported for now");
std::unique_ptr<real_number[]> extra_rhs(new real_number[getLocalNbParticles()*3]());
sample_compute_field(in_field, extra_rhs.get());
completeLoopWithVorticity(dt, extra_rhs.get());
}
}; };
#endif #endif
...@@ -2,27 +2,31 @@ ...@@ -2,27 +2,31 @@
#define P2P_COMPUTER_HPP #define P2P_COMPUTER_HPP
#include <cstring> #include <cstring>
#include <cassert>
/** \brief A simple distance weighting function.
*
* This function returns 1 if a distance is smaller than a cut-off length,
* i.e. particle 1 interacts with particle 2 if particle 2 is inside a
* sphere of radius `cutoff' centered on particle 1.
*/
double dumb_distance_weight(
double distance,
double cutoff)
{
// this function should only be called for interacting particles,
// and particles interact if they are closer than cutoff.
assert(distance < cutoff);
return 1.0;
}
template <class real_number, class partsize_t> template <class real_number, class partsize_t>
class p2p_computer{ class p2p_computer{
bool isActive;
/** \brief A simple distance weighting function.
*
* This function returns 1 if a distance is smaller than a cut-off length,
* i.e. particle 1 interacts with particle 2 if particle 2 is inside a
* sphere of radius `cutoff' centered on particle 1.
*/
static double dumb_distance_weight(
const double dist_pow2,
const double cutoff){
// this function should only be called for interacting particles,
// and particles interact if they are closer than cutoff.
assert(dist_pow2 < cutoff*cutoff);
return 1.0;
}
public: public:
constexpr static int size_data = 3; p2p_computer() : isActive(true){}
template <int size_particle_rhs> template <int size_particle_rhs>
void init_result_array(real_number rhs[], const partsize_t nbParticles) const{ void init_result_array(real_number rhs[], const partsize_t nbParticles) const{
...@@ -31,25 +35,26 @@ public: ...@@ -31,25 +35,26 @@ public:
template <int size_particle_rhs> template <int size_particle_rhs>
void reduce_particles_rhs(real_number rhs_dst[], const real_number rhs_src[], const partsize_t nbParticles) const{ void reduce_particles_rhs(real_number rhs_dst[], const real_number rhs_src[], const partsize_t nbParticles) const{
static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
for(int idx_part = 0 ; idx_part < nbParticles ; ++idx_part){ for(int idx_part = 0 ; idx_part < nbParticles ; ++idx_part){
for(int idx_rhs = 0 ; idx_rhs < size_particle_rhs ; ++idx_rhs){ // We merge only the values modified by the current kernel (3-5)
for(int idx_rhs = 3 ; idx_rhs < size_particle_rhs ; ++idx_rhs){
rhs_dst[idx_part*size_particle_rhs+idx_rhs] += rhs_src[idx_part*size_particle_rhs+idx_rhs]; rhs_dst[idx_part*size_particle_rhs+idx_rhs] += rhs_src[idx_part*size_particle_rhs+idx_rhs];
} }
} }
} }
template <int size_particle_positions, int size_particle_data, int size_particle_rhs> template <int size_particle_positions, int size_particle_rhs>
void compute_interaction(const real_number pos_part1[], const real_number data_part1[], real_number rhs_part1[], void compute_interaction(const real_number pos_part1[], real_number rhs_part1[],
const real_number pos_part2[], const real_number data_part2[], real_number rhs_part2[], const real_number pos_part2[], real_number rhs_part2[],
const real_number dist_pow2, const real_number dist_pow2, const real_number cutoff,
const real_number xshift_coef, const real_number yshift_coef, const real_number zshift_coef) const{ const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{
static_assert(size_particle_positions == 3, "This kernel works only with 3 values for one particle's position"); static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position+orientation");
static_assert(size_particle_rhs >= 6, "This kernel works only with more than 6 values per particle's rhs"); static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
static_assert(size_particle_data == 3, "This kernel works only with 3 values per particle's' data");
// TODO: a reasonable way of choosing between different distance_weight functions should be thought of. // TODO: a reasonable way of choosing between different distance_weight functions should be thought of.
// We need to ask Michael about how flexible this distance_weight needs to be. // We need to ask Michael about how flexible this distance_weight needs to be.
const double ww = dumb_distance_weight(distance, max_distance); const double ww = dumb_distance_weight(dist_pow2, cutoff);
/// ///
/// term in equation is: /// term in equation is:
/// ///
...@@ -57,20 +62,23 @@ public: ...@@ -57,20 +62,23 @@ public:
/// (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j /// (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j
/// \f] /// \f]
/// ///
const double dot_product = (data_part1[0]*data_part2[0] + const double dot_product = (pos_part1[3+IDX_X]*pos_part2[3+IDX_X] +
data_part1[1]*data_part2[1] + pos_part1[3+IDX_Y]*pos_part2[3+IDX_Y] +
data_part1[2]*data_part2[2]); pos_part1[3+IDX_Z]*pos_part2[3+IDX_Z]);
rhs_part1[3] += data_part2[0] * 4 * ww * dot_product; rhs_part1[3+IDX_X] += pos_part2[3+IDX_X] * 4 * ww * dot_product;
rhs_part1[4] += data_part2[1] * 4 * ww * dot_product; rhs_part1[3+IDX_Y] += pos_part2[3+IDX_Y] * 4 * ww * dot_product;
rhs_part1[5] += data_part2[2] * 4 * ww * dot_product; rhs_part1[3+IDX_Z] += pos_part2[3+IDX_Z] * 4 * ww * dot_product;
rhs_part2[3] += data_part1[0] * 4 * ww * dot_product; rhs_part2[3+IDX_X] += pos_part1[3+IDX_X] * 4 * ww * dot_product;
rhs_part2[4] += data_part1[1] * 4 * ww * dot_product; rhs_part2[3+IDX_Y] += pos_part1[3+IDX_Y] * 4 * ww * dot_product;
rhs_part2[5] += data_part1[2] * 4 * ww * dot_product; rhs_part2[3+IDX_Z] += pos_part1[3+IDX_Z] * 4 * ww * dot_product;
} }
bool isEnable() const {
return isActive;
}
constexpr static bool isEmpty() { void setEnable(const bool inIsActive) {
return false; isActive = inIsActive;
} }
}; };
......
...@@ -6,8 +6,6 @@ ...@@ -6,8 +6,6 @@
template <class real_number, class partsize_t> template <class real_number, class partsize_t>
class p2p_computer_empty{ class p2p_computer_empty{
public: public:
constexpr int static size_data = 0;
template <int size_particle_rhs> template <int size_particle_rhs>
void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{ void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
} }
...@@ -16,15 +14,15 @@ public: ...@@ -16,15 +14,15 @@ public:
void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{ void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
} }
template <int size_particle_positions, int size_particle_data, int size_particle_rhs> template <int size_particle_positions, int size_particle_rhs>
void compute_interaction(const real_number /*pos_part1*/[], const real_number /*data_part1*/[], real_number /*rhs_part1*/[], void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
const real_number /*pos_part2*/[], const real_number /*data_part2*/[], real_number /*rhs_part2*/[], const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
const real_number /*dist_pow2*/, const real_number /*dist_pow2*/, const real_number /*cutoff*/,
const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{ const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{
} }
constexpr static bool isEmpty() { constexpr static bool isEnable() {
return true; return false;
} }
}; };
......
...@@ -36,7 +36,6 @@ protected: ...@@ -36,7 +36,6 @@ protected:
std::unique_ptr<real_number[]> toRecvAndMerge; std::unique_ptr<real_number[]> toRecvAndMerge;
std::unique_ptr<real_number[]> toCompute; std::unique_ptr<real_number[]> toCompute;
std::unique_ptr<real_number[]> toData;
std::unique_ptr<real_number[]> results; std::unique_ptr<real_number[]> results;
}; };
...@@ -248,11 +247,10 @@ public: ...@@ -248,11 +247,10 @@ public:
return (diff_x*diff_x) + (diff_y*diff_y) + (diff_z*diff_z); return (diff_x*diff_x) + (diff_y*diff_y) + (diff_z*diff_z);
} }
template <class computer_class, int size_particle_positions, int size_particle_data, int size_particle_rhs> template <class computer_class, int size_particle_positions, int size_particle_rhs>
void compute_distr(computer_class& in_computer, void compute_distr(computer_class& in_computer,
const partsize_t current_my_nb_particles_per_partition[], const partsize_t current_my_nb_particles_per_partition[],
real_number particles_positions[], real_number particles_positions[],
real_number particles_data[],
real_number particles_current_rhs[], real_number particles_current_rhs[],
partsize_t inout_index_particles[]){ partsize_t inout_index_particles[]){
TIMEZONE("compute_distr"); TIMEZONE("compute_distr");
...@@ -313,9 +311,6 @@ public: ...@@ -313,9 +311,6 @@ public:
permute_copy<real_number, size_particle_positions>(current_offset_particles_for_partition[idxPartition], permute_copy<real_number, size_particle_positions>(current_offset_particles_for_partition[idxPartition],
current_my_nb_particles_per_partition[idxPartition], current_my_nb_particles_per_partition[idxPartition],
part_to_sort.data(), particles_positions, &buffer); part_to_sort.data(), particles_positions, &buffer);
permute_copy<real_number, size_particle_data>(current_offset_particles_for_partition[idxPartition],
current_my_nb_particles_per_partition[idxPartition],
part_to_sort.data(), particles_data, &buffer);
permute_copy<real_number, size_particle_rhs>(current_offset_particles_for_partition[idxPartition], permute_copy<real_number, size_particle_rhs>(current_offset_particles_for_partition[idxPartition],
current_my_nb_particles_per_partition[idxPartition], current_my_nb_particles_per_partition[idxPartition],
part_to_sort.data(), particles_current_rhs, &buffer); part_to_sort.data(), particles_current_rhs, &buffer);
...@@ -457,14 +452,6 @@ public: ...@@ -457,14 +452,6 @@ public:
descriptor.destProc, TAG_POSITION_PARTICLES, descriptor.destProc, TAG_POSITION_PARTICLES,
current_com, &mpiRequests.back())); current_com, &mpiRequests.back()));
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(descriptor.nbParticlesToExchange*size_particle_data < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_data[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_data]),
int(descriptor.nbParticlesToExchange*size_particle_data), particles_utils::GetMpiType(real_number()),
descriptor.destProc, TAG_POSITION_PARTICLES,
current_com, &mpiRequests.back()));
assert(descriptor.toRecvAndMerge == nullptr); assert(descriptor.toRecvAndMerge == nullptr);
descriptor.toRecvAndMerge.reset(new real_number[descriptor.nbParticlesToExchange*size_particle_rhs]); descriptor.toRecvAndMerge.reset(new real_number[descriptor.nbParticlesToExchange*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr}); whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
...@@ -526,15 +513,6 @@ public: ...@@ -526,15 +513,6 @@ public:
AssertMpi(MPI_Irecv(descriptor.toCompute.get(), int(NbParticlesToReceive*size_particle_positions), AssertMpi(MPI_Irecv(descriptor.toCompute.get(), int(NbParticlesToReceive*size_particle_positions),
particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES, particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES,
current_com, &mpiRequests.back())); current_com, &mpiRequests.back()));
descriptor.toData.reset(new real_number[NbParticlesToReceive*size_particle_data]);
whatNext.emplace_back(std::pair<Action,int>{COMPUTE_PARTICLES, releasedAction.second});
mpiRequests.emplace_back();
assert(NbParticlesToReceive*size_particle_data < std::numeric_limits<int>::max());
AssertMpi(MPI_Irecv(descriptor.toData.get(), int(NbParticlesToReceive*size_particle_data),
particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES,
current_com, &mpiRequests.back()));
} }
} }
...@@ -548,7 +526,6 @@ public: ...@@ -548,7 +526,6 @@ public:
const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange; const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
assert(descriptor.toCompute != nullptr); assert(descriptor.toCompute != nullptr);
assert(descriptor.toData != nullptr);
assert(descriptor.positionsReceived == true); assert(descriptor.positionsReceived == true);
descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]); descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
in_computer.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive); in_computer.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
...@@ -589,14 +566,12 @@ public: ...@@ -589,14 +566,12 @@ public:
particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z], particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]); shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){ if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
in_computer.template compute_interaction<size_particle_positions,size_particle_data, size_particle_rhs>( in_computer.template compute_interaction<size_particle_positions, size_particle_rhs>(
&descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions], &descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions],
&descriptor.toData[(idxPart+idx_p1)*size_particle_data],
&descriptor.results[(idxPart+idx_p1)*size_particle_rhs], &descriptor.results[(idxPart+idx_p1)*size_particle_rhs],
&particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions], &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
&particles_data[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_data],
&particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs], &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
dist_r2, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]); dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
} }
} }
} }
...@@ -681,14 +656,12 @@ public: ...@@ -681,14 +656,12 @@ public:
particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Z], particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Z],
0, 0, 0); 0, 0, 0);
if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){ if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
in_computer.template compute_interaction<size_particle_positions,size_particle_data,size_particle_rhs>( in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
&particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions], &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
&particles_data[(intervals[idx_1].first+idx_p1)*size_particle_data],
&particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs], &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
&particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions], &particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions],
&particles_data[(intervals[idx_1].first+idx_p2)*size_particle_data],
&particles_current_rhs[(intervals[idx_1].first+idx_p2)*size_particle_rhs], &particles_current_rhs[(intervals[idx_1].first+idx_p2)*size_particle_rhs],
dist_r2, 0, 0, 0); dist_r2, cutoff_radius_compute, 0, 0, 0);
} }
} }
} }
...@@ -705,14 +678,12 @@ public: ...@@ -705,14 +678,12 @@ public:
particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Z], particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
0, 0, 0); 0, 0, 0);
if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){ if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
in_computer.template compute_interaction<size_particle_positions,size_particle_data,size_particle_rhs>( in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
&particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions], &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
&particles_data[(intervals[idx_1].first+idx_p1)*size_particle_data],
&particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs], &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
&particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions], &particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions],
&particles_data[(intervals[idx_2].first+idx_p2)*size_particle_data],
&particles_current_rhs[(intervals[idx_2].first+idx_p2)*size_particle_rhs], &particles_current_rhs[(intervals[idx_2].first+idx_p2)*size_particle_rhs],
dist_r2, 0, 0, 0); dist_r2, cutoff_radius_compute, 0, 0, 0);
} }
} }
} }
...@@ -741,14 +712,12 @@ public: ...@@ -741,14 +712,12 @@ public:
particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z], particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]); shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){ if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
in_computer.template compute_interaction<size_particle_positions,size_particle_data,size_particle_rhs>( in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
&particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions], &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
&particles_data[(intervals[idx_1].first+idx_p1)*size_particle_data],
&particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs], &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
&particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions], &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
&particles_data[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_data],
&particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs], &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
dist_r2, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]); dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
} }
} }
} }
......
...@@ -7,11 +7,30 @@ ...@@ -7,11 +7,30 @@
#include "scope_timer.hpp" #include "scope_timer.hpp"
#include "particles_utils.hpp" #include "particles_utils.hpp"
template <class partsize_t, class real_number, int size_particle_positions = 3, int size_particle_rhs = 3> template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
class particles_adams_bashforth { class particles_adams_bashforth;
static_assert(size_particle_positions == size_particle_rhs,