p2p_computer.hpp 3.22 KB
 Berenger Bramas committed Sep 20, 2017 1 2 3 4 5 #ifndef P2P_COMPUTER_HPP #define P2P_COMPUTER_HPP #include  Cristian Lalescu committed Oct 10, 2017 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; }  Berenger Bramas committed Sep 20, 2017 22 23 24 template class p2p_computer{ public:  Berenger Bramas committed Sep 22, 2017 25 26  constexpr static int size_data = 3;  Berenger Bramas committed Sep 20, 2017 27 28  template void init_result_array(real_number rhs[], const partsize_t nbParticles) const{  Berenger Bramas committed Sep 21, 2017 29  memset(rhs, 0, sizeof(real_number)*nbParticles*size_particle_rhs);  Berenger Bramas committed Sep 20, 2017 30 31 32 33 34 35 36 37 38 39 40  } template 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]; } } }  Berenger Bramas committed Sep 22, 2017 41 42 43  template 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 committed Sep 21, 2017 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 committed Oct 10, 2017 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");  Cristian Lalescu committed Oct 10, 2017 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 committed Oct 10, 2017 52  const double ww = dumb_distance_weight(distance, max_distance);  Cristian Lalescu committed Oct 10, 2017 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 committed Oct 10, 2017 60  const double dot_product = (data_part1[0]*data_part2[0] +  Cristian Lalescu committed Oct 10, 2017 61  data_part1[1]*data_part2[1] +  Berenger Bramas committed Oct 10, 2017 62  data_part1[2]*data_part2[2]);  Cristian Lalescu committed Oct 10, 2017 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;  Berenger Bramas committed Sep 20, 2017 69  }  Berenger Bramas committed Sep 22, 2017 70   Cristian Lalescu committed Oct 10, 2017 71   Berenger Bramas committed Sep 22, 2017 72 73 74  constexpr static bool isEmpty() { return false; }  Berenger Bramas committed Sep 20, 2017 75 76 77 }; #endif`