p2p_computer.hpp 3.22 KB
Newer Older
1
2
3
4
5
#ifndef P2P_COMPUTER_HPP
#define P2P_COMPUTER_HPP

#include <cstring>

6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
/** \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;
}

22
23
24
template <class real_number, class partsize_t>
class p2p_computer{
public:
25
26
    constexpr static int size_data = 3;

27
28
    template <int size_particle_rhs>
    void init_result_array(real_number rhs[], const partsize_t nbParticles) const{
Berenger Bramas's avatar
Berenger Bramas committed
29
        memset(rhs, 0, sizeof(real_number)*nbParticles*size_particle_rhs);
30
31
32
33
34
35
36
37
38
39
40
    }

    template <int size_particle_rhs>
    void reduce_particles_rhs(real_number rhs_dst[], const real_number rhs_src[], const partsize_t nbParticles) const{
        for(int idx_part = 0 ; idx_part < nbParticles ; ++idx_part){
            for(int idx_rhs = 0 ; 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];
            }
        }
    }

41
42
43
    template <int size_particle_positions, int size_particle_data, int size_particle_rhs>
    void compute_interaction(const real_number pos_part1[], const real_number data_part1[], real_number rhs_part1[],
                             const real_number pos_part2[], const real_number data_part2[], real_number rhs_part2[],
Berenger Bramas's avatar
Berenger Bramas committed
44
45
                             const real_number dist_pow2,
                             const real_number xshift_coef, const real_number yshift_coef, const real_number zshift_coef) const{
Berenger Bramas's avatar
Berenger Bramas committed
46
47
48
49
        static_assert(size_particle_positions == 3, "This kernel works only with 3 values for one particle's position");
        static_assert(size_particle_rhs >= 6, "This kernel works only with more than 6 values per particle's rhs");
        static_assert(size_particle_data == 3, "This kernel works only with 3 values per particle's' data");

50
51
        // 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.
Berenger Bramas's avatar
Berenger Bramas committed
52
        const double ww = dumb_distance_weight(distance, max_distance);
53
54
55
56
57
58
59
        ///
        /// term in equation is:
        ///
        /// \f[
        ///     (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j
        /// \f]
        ///
Berenger Bramas's avatar
Berenger Bramas committed
60
        const double dot_product = (data_part1[0]*data_part2[0] +
61
                              data_part1[1]*data_part2[1] +
Berenger Bramas's avatar
Berenger Bramas committed
62
                              data_part1[2]*data_part2[2]);
63
64
65
66
67
68
        rhs_part1[3] += data_part2[0] * 4 * ww * dot_product;
        rhs_part1[4] += data_part2[1] * 4 * ww * dot_product;
        rhs_part1[5] += data_part2[2] * 4 * ww * dot_product;
        rhs_part2[3] += data_part1[0] * 4 * ww * dot_product;
        rhs_part2[4] += data_part1[1] * 4 * ww * dot_product;
        rhs_part2[5] += data_part1[2] * 4 * ww * dot_product;
69
    }
70

71

72
73
74
    constexpr static bool isEmpty() {
        return false;
    }
75
76
77
};

#endif