Commit ccb4319f authored by Quentin Khan's avatar Quentin Khan

FAdaptiveSequential: creation, adaptive FMM implementation

parent f37e2926
#ifndef SCALFMM_SEQUENTIAL_ALGO_HPP_
#define SCALFMM_SEQUENTIAL_ALGO_HPP_
#include <cmath> // Used to round box differences
#include "Containers/FTreeCoordinate.hpp"
template<class _Tree, class _Kernel>
class FAdaptiveSequential {
public:
using tree_t = _Tree;
using kernel_t = _Kernel;
using node_t = typename tree_t::node_t;
private:
tree_t& _tree;
kernel_t& _kernel;
template<typename K, typename Ret, typename... Args>
using has_setup = typename std::enable_if<
std::is_same<Ret, decltype(std::declval<K>().setup(std::declval<Args>()...))>::value
>::type*;
template<typename K, has_setup<K, void> = nullptr>
void setup_kernel() {
_kernel.setup(_tree);
}
void setup_kernel(...) {}
public:
FAdaptiveSequential(tree_t& tree, kernel_t& kernel) :
_tree(tree),
_kernel(kernel) {
}
FAdaptiveSequential(tree_t& tree, kernel_t* kernel) :
FAdaptiveSequential(tree, *kernel)
{}
FAdaptiveSequential(tree_t* tree, kernel_t& kernel) :
FAdaptiveSequential(*tree, kernel)
{}
FAdaptiveSequential(tree_t* tree, kernel_t* kernel) :
FAdaptiveSequential(*tree, *kernel)
{}
void execute() {
this->run();
}
void run() {
this->setup_kernel();
// 1. source to up, P2M
source_to_up();
// 2. up to up, M2M
up_to_up();
// 3a V-list, M2L
v_list_step();
// 3b X-list, P2L
x_list_step();
// 4. down to down, L2L
down_to_down();
// 5a W-list, M2P
w_list_step();
// 5b down to target, L2P
down_to_target();
// A. U-list, P2P
u_list_step();
}
// P2M
void source_to_up() {
for(node_t* leaf : _tree.leaves()) {
_kernel.P2M(leaf->getData(), leaf->getParticleContainer());
}
}
// M2M
void up_to_up() {
for(node_t& node : _tree.post_order_walk()) {
if(! node.is_leaf()) {
// Fill array of children data
typename node_t::data_t* child_data[8] = {};
for(node_t* child : node.getChildren()) {
child_data[child->getIndex() & (node_t::child_count-1)] = child->getData();
}
// Call kernel module
_kernel.M2M(node.getData(), child_data, static_cast<int>(node.getDepth()));
}
}
}
// M2L
void v_list_step() {
for(node_t& node : _tree.in_order_walk()) {
// Needed to compute offset between boxes
for(node_t* v_item : node.V) {
int v_item_index[1] = {compute_box_offset_index(&node, v_item, 3)};
const typename node_t::data_t* v_item_data[1] = {v_item->getData()};
// Call kernel M2L operator
_kernel.M2L(node.getData(), v_item_data, v_item_index, 1, static_cast<int>(node.getDepth()));
}
}
}
// P2L
void x_list_step() {
for(node_t& node : _tree.in_order_walk()) {
for(node_t* x_item : node.X) {
_kernel.P2L(node.getData(), x_item->getParticleContainer());
}
}
}
// L2L
void down_to_down() {
for(node_t& node : _tree.pre_order_walk()) {
if(! node.is_leaf()) {
typename node_t::data_t* child_data[node_t::child_count];
for(std::size_t i = 0; i < node_t::child_count; ++i) {
child_data[i] = node.getChild(i)->getData();
}
_kernel.L2L(node.getData(), child_data, static_cast<int>(node.getDepth()));
}
}
}
// M2P
void w_list_step() {
for(node_t& node : _tree.in_order_walk()) {
for(node_t* w_item : node.W) {
_kernel.M2P(w_item->getData(), node.getParticleContainer());
}
}
}
// L2P
void down_to_target() {
for(node_t* leaf : _tree.leaves()) {
_kernel.L2P(leaf->getData(), leaf->getParticleContainer());
}
}
// P2P
void u_list_step() {
using container_t = typename node_t::particle_container_t;
for(node_t* leaf : _tree.leaves()) {
for(node_t* u_item : leaf->U) {
int u_item_index = compute_box_offset_index(leaf, u_item, 1);
container_t* const leaf_source_particle_container =
leaf->getParticleContainer();
const container_t* const leaf_target_particle_container =
leaf->getParticleContainer();
container_t* const u_item_source_particle_container[1] =
{u_item->getParticleContainer()};
_kernel.P2P(FTreeCoordinate(MortonIndex(leaf->getIndex()), static_cast<int>(leaf->getDepth())),
leaf_source_particle_container,
leaf_target_particle_container,
u_item_source_particle_container,
&u_item_index,
1);
}
}
}
/** Compute resulting index in pre-computation array.
*
* Kernels precompute some operators according to node offsets. This method
* returns the `other_node` index in the precomputation array.
*
* The index has the following form, with d the space dimension:
*
* x+n * 7^d + y+n * 7^(d-1) + ... + z+n
*
* Because an `other_node` is at most `n` boxes away in every direction from a
* `node`, in order to get a positive index, we add `n` to the offset in every
* direction.
*
* \param node The target node
* \param other_node The interacting node
* \param n The longest possible distance between node and other_node in
* terms of boxes on an axis
*/
int compute_box_offset_index(node_t* node, node_t* other_node, const std::size_t n) {
typename node_t::FReal box_width = (node->getBox().c2 - node->getBox().c1)[0];
auto center_offset = other_node->getBox().center() - node->getBox().center();
std::size_t other_node_index = 0;
for(std::size_t i = 0; i < node_t::Dim; ++i) {
other_node_index = other_node_index * 7
+ static_cast<unsigned long>(std::lround(center_offset[i] / box_width))
+ n;
}
return static_cast<int>(other_node_index);
}
};
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment