-
Cristian Lalescu authoredCristian Lalescu authored
p2p_computer.hpp 5.42 KiB
/******************************************************************************
* *
* Copyright 2019 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 P2P_COMPUTER_HPP
#define P2P_COMPUTER_HPP
#include <cstring>
#include <cassert>
template <class real_number, class partsize_t>
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:
p2p_computer() : isActive(true){}
template <int size_particle_rhs>
void init_result_array(real_number rhs[], const partsize_t nbParticles) const{
memset(rhs, 0, sizeof(real_number)*nbParticles*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{
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){
// 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];
}
}
}
template <int size_particle_positions, int size_particle_rhs>
void compute_interaction(const partsize_t /*idx_part1*/,
const real_number pos_part1[], real_number rhs_part1[],
const partsize_t /*idx_part2*/,
const real_number pos_part2[], real_number rhs_part2[],
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{
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 6 values per particle's rhs");
// 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.
const double ww = dumb_distance_weight(dist_pow2, cutoff);
///
/// term in equation is:
///
/// \f[
/// (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j
/// \f]
///
const double dot_product = (pos_part1[3+IDXC_X]*pos_part2[3+IDXC_X] +
pos_part1[3+IDXC_Y]*pos_part2[3+IDXC_Y] +
pos_part1[3+IDXC_Z]*pos_part2[3+IDXC_Z]);
rhs_part1[3+IDXC_X] += pos_part2[3+IDXC_X] * 4 * ww * dot_product;
rhs_part1[3+IDXC_Y] += pos_part2[3+IDXC_Y] * 4 * ww * dot_product;
rhs_part1[3+IDXC_Z] += pos_part2[3+IDXC_Z] * 4 * ww * dot_product;
rhs_part2[3+IDXC_X] += pos_part1[3+IDXC_X] * 4 * ww * dot_product;
rhs_part2[3+IDXC_Y] += pos_part1[3+IDXC_Y] * 4 * ww * dot_product;
rhs_part2[3+IDXC_Z] += pos_part1[3+IDXC_Z] * 4 * ww * dot_product;
}
void merge(const p2p_computer&){}
bool isEnable() const {
return isActive;
}
void setEnable(const bool inIsActive) {
isActive = inIsActive;
}
};
#endif