FAdaptiveSequential.hpp 14.7 KB
Newer Older
1 2 3
#ifndef SCALFMM_SEQUENTIAL_ALGO_HPP_
#define SCALFMM_SEQUENTIAL_ALGO_HPP_

4
#include <algorithm>
5
#include <cmath> // Used to round box differences
6
#include <functional>
7
#include <map>
8
#include <vector>
9

10

11
#include "Core/FCoreCommon.hpp"
12
#include "Containers/FTreeCoordinate.hpp"
13 14
#include "Utils/FAlgorithmTimers.hpp"

15 16
#include "Kernels/FKernelConcepts.hpp"

17 18 19 20 21 22
template<
    class _Tree, class _Kernel,
    class = inria::require<
        scalfmm::meta::adaptive_compatible<_Tree,_Kernel>
        >
    >
23
class FAdaptiveSequential : public FAlgorithmInterface, public FAlgorithmTimers {
24 25 26
public:
    using tree_t = _Tree;
    using kernel_t = _Kernel;
27
    using FReal = typename tree_t::FReal;
28 29 30

    using node_t = typename tree_t::node_t;

31 32 33 34
    using multipole_t = typename node_t::data_t::multipole_t;
    using local_expansion_t = typename node_t::data_t::local_expansion_t;
    using symbolic_data_t = typename node_t::symbolic_data_t;

35 36 37 38 39 40 41 42 43 44
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*;


45
    template<typename K, has_setup<K, void, tree_t> = nullptr>
46
    void setup_kernel(K*) {
47 48 49
        _kernel.setup(_tree);
    }

50
    template<typename K, has_setup<K, void> = nullptr>
51
    void setup_kernel(K*) {
52 53 54
        _kernel.setup();
    }

55 56 57 58
    void setup_kernel(...) {}

public:

59 60 61
    FAdaptiveSequential(tree_t* tree, kernel_t* kernel) :
        _tree(*tree),
        _kernel(*kernel) {
62 63
    }

64 65 66 67 68 69 70 71
    std::string name() const override {
        return "Sequential adaptive algorithm";
    }

    std::string description() const override {
        return "";
    }

72 73
    using FAlgorithmInterface::execute;

74 75 76 77 78 79
    /** \brief Run specific steps of the algorithm
     *
     * \param operations Specifies the algorithm operations to run, see
     * FFmmOperations.
     */
    void execute(const unsigned int operations) override {
80
        this->run(operations);
81 82
    }

83
    void run(int operations = FFmmNearAndFarFields) {
84

85
        this->setup_kernel(&(this->_kernel));
86

87 88
        if(operations & FFmmP2M) {
            // 1. source to up, P2M
89 90 91 92
            Timers["P2M"].tic();
            this->source_to_up();;
            Timers["P2M"].tac();
            std::cout << "    P2M: " << Timers["P2M"].elapsed() << std::endl;
93
        }
94

95 96
        if(operations & FFmmM2M) {
            // 2. up to up, M2M
97 98 99 100
            Timers["M2M"].tic();
            this->up_to_up();
            Timers["M2M"].tac();
            std::cout << "    M2M: " << Timers["M2M"].elapsed() << std::endl;
101
        }
102

103 104
        if(operations & FFmmM2L) {
            // 3a V-list, M2L
105 106 107 108
            Timers["M2L"].tic();
            this->v_list_step();
            Timers["M2L"].tac();
            std::cout << "    M2L: " << Timers["M2L"].elapsed() << std::endl;
109
        }
110

111 112
        if(operations & FFmmP2L) {
            // 3b X-list, P2L
113 114 115 116
            Timers["P2L"].tic();
            this->x_list_step();
            Timers["P2L"].tac();
            std::cout << "    P2L: " << Timers["P2L"].elapsed() << std::endl;
117
        }
118

119 120
        if(operations & FFmmL2L) {
            // 4. down to down, L2L
121 122 123 124
            Timers["L2L"].tic();
            this->down_to_down();
            Timers["L2L"].tac();
            std::cout << "    L2L: " << Timers["L2L"].elapsed() << std::endl;
125
        }
126

127
        if(operations & FFmmM2P) {
128
        // 5a W-list, M2P
129 130 131 132
            Timers["M2P"].tic();
            this->w_list_step();
            Timers["M2P"].tac();
            std::cout << "    M2P: " << Timers["M2P"].elapsed() << std::endl;
133
        }
134

135
        if(operations & FFmmL2P) {
136
            // 5b down to target, L2P
137 138 139 140
            Timers["L2P"].tic();
            this->down_to_target();
            Timers["L2P"].tac();
            std::cout << "    L2P: " << Timers["L2P"].elapsed() << std::endl;
141
        }
142

143 144
        if(operations & FFmmP2P) {
            // A. U-list, P2P
145 146 147 148
            Timers["P2P"].tic();
            this->u_list_step();
            Timers["P2P"].tac();
            std::cout << "    P2P: " << Timers["P2P"].elapsed() << std::endl;
149
        }
150 151 152 153 154
    }

    // P2M
    void source_to_up() {
        for(node_t* leaf : _tree.leaves()) {
155
            if(leaf->getParticleContainer()->size() != 0) {
156 157 158
                _kernel.P2M(&(leaf->getData()->getMultipoleData()),
                            &(leaf->getSymbolicData()),
                            leaf->getParticleContainer());
159
            }
160 161 162 163 164 165 166
        }
    }

    // M2M
    void up_to_up() {
        for(node_t& node : _tree.post_order_walk()) {
            if(! node.is_leaf()) {
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182
                // Fill array of child data
                std::array<multipole_t*, node_t::child_count> child_multipoles {};
                std::transform(std::begin(node.getChildren()), std::end(node.getChildren()),
                               child_multipoles.begin(),
                               [](node_t* n) {
                                   return n == nullptr ? nullptr
                                       : &(n->getData()->getMultipoleData());
                               });

                std::array<symbolic_data_t*, node_t::child_count> child_symbolics {};
                std::transform(std::begin(node.getChildren()), std::end(node.getChildren()),
                               child_symbolics.begin(),
                               [](node_t* n) {
                                   return n == nullptr ? nullptr
                                       : &(n->getSymbolicData());
                               });
183 184

                // Call kernel module
185
                _kernel.M2M(&(node.getData()->getMultipoleData()),
186 187 188 189
                            &(node.getSymbolicData()),
                            child_multipoles.data(),
                            child_symbolics.data()
                    );
190 191 192 193 194 195
            }
        }
    }

    // M2L
    void v_list_step() {
196 197
        std::vector<multipole_t*> v_item_multipoles;
        std::vector<symbolic_data_t*> v_item_symbolics;
198 199
        std::vector<int> v_item_indices;

200
        for(node_t& node : _tree.in_order_walk()) {
201 202 203 204
            if(node.is_leaf() && node.getParticleContainer()->size() == 0) {
                continue;
            }

205 206
            v_item_multipoles.clear();
            v_item_symbolics.clear();
207 208
            v_item_indices.clear();

209 210
            // Needed to compute offset between boxes
            for(node_t* v_item : node.V) {
211 212 213 214
                if(v_item->is_leaf()
                   && v_item->getParticleContainer()->size() == 0) {
                    continue;
                }
215
                v_item_indices.push_back(compute_box_offset_index(&node, v_item, 3));
216 217
                v_item_multipoles.push_back(&(v_item->getData()->getMultipoleData()));
                v_item_symbolics.push_back(&(v_item->getSymbolicData()));
218
            }
219 220

            // Call kernel M2L operator
221
            _kernel.M2L(&(node.getData()->getLocalExpansionData()),
222 223 224
                        &(node.getSymbolicData()),
                        v_item_multipoles.data(),
                        v_item_symbolics.data(),
225
                        v_item_indices.data(),
226 227
                        static_cast<int>(v_item_multipoles.size())
                );
228 229 230 231 232
        }
    }

    // P2L
    void x_list_step() {
233 234 235 236 237
        /* NOTE: the X list and W list are complementary: if A is in X(B) then B
         * is in W(A).
         *
         * We loop over the leaves first to detect the empty ones early on.
         */
238 239 240
        for(node_t* leaf : _tree.leaves()) {
            if(leaf->getParticleContainer()->size() > 0) {
                for(node_t* w_item : leaf->W) {
241 242 243
                    if(w_item->is_leaf() && w_item->getParticleContainer()->size() == 0) {
                        continue;
                    }
244 245
                    _kernel.P2L(&(w_item->getData()->getLocalExpansionData()),
                                &(w_item->getSymbolicData()),
246 247 248
                                leaf->getParticleContainer(),
                                &(leaf->getSymbolicData())
                        );
249
                }
250 251 252 253 254 255 256 257
            }
        }
    }

    // L2L
    void down_to_down() {
        for(node_t& node : _tree.pre_order_walk()) {
            if(! node.is_leaf()) {
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274

                std::array<local_expansion_t*, node_t::child_count> child_local_expansions {};
                std::transform(std::begin(node.getChildren()), std::end(node.getChildren()),
                               child_local_expansions.begin(),
                               [](node_t* n) {
                                   return n == nullptr ? nullptr
                                       : &(n->getData()->getLocalExpansionData());
                               });

                std::array<symbolic_data_t*, node_t::child_count> child_symbolics {};
                std::transform(std::begin(node.getChildren()), std::end(node.getChildren()),
                               child_symbolics.begin(),
                               [](node_t* n) {
                                   return n == nullptr ? nullptr
                                       : &(n->getSymbolicData());
                               });

275
                _kernel.L2L(&(node.getData()->getLocalExpansionData()),
276 277 278 279
                            &(node.getSymbolicData()),
                            child_local_expansions.data(),
                            child_symbolics.data()
                    );
280 281 282 283 284 285
            }
        }
    }

    // M2P
    void w_list_step() {
286 287 288
        for(node_t* leaf : _tree.leaves()) {
            if(leaf->getParticleContainer()->size() > 0) {
                for(node_t* w_item : leaf->W) {
289 290 291
                    if(w_item->is_leaf() && w_item->getParticleContainer()->size() == 0) {
                        continue;
                    }
292 293
                    _kernel.M2P(&(w_item->getData()->getMultipoleData()),
                                &(w_item->getSymbolicData()),
294 295 296
                                leaf->getParticleContainer(),
                                &(leaf->getSymbolicData())
                        );
297
                }
298 299 300 301 302 303 304
            }
        }
    }

    // L2P
    void down_to_target() {
        for(node_t* leaf : _tree.leaves()) {
305
            if(leaf->getParticleContainer()->size() != 0) {
306 307 308
                _kernel.L2P(&(leaf->getData()->getLocalExpansionData()),
                            &(leaf->getSymbolicData()),
                            leaf->getParticleContainer());
309
            }
310 311 312
        }
    }

313

314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339
    /** \brief Direct computation step (P2P)
     *
     * For each tree leaf, a direct computation is done with the adjacent
     * leaves.
     *
     * The symmetric computation kernels expect to receive the list of adjacent
     * leaves, stored in a leaf's U list, sorted by offset index.
     *
     * The offset index is of a node compared to another follows the
     * numerotation of the adjacent nodes, on each axis. In 3D, this woud give,
     * with L the reference leaf:
     *
     *          x          x            x
     *        0 1 2  |   9 10 11  |  18 19 20
     *     y  3 4 5  |  12  L 14  |  21 22 23
     *        6 7 8  |  15 16 17  |  24 25 26
     *         z=0        z=1          z=2
     *
     * When two node that are to be compared are not on the same level, the
     * lowest one's ancestor that is at the highest one's level is used to
     * compute the index.
     *
     * To statisfy this condition, the leaves have to be sorted before being
     * given to the kernel P2P method.
     *
     */
340 341
    void u_list_step() {
        using container_t = typename node_t::particle_container_t;
342

343
        // Vectors to pass to kernel P2P
344 345
        std::vector<container_t*> u_item_source_particle_containers;
        std::vector<int> u_item_indices;
346

347
        for(node_t* leaf : _tree.leaves()) {
348

349 350
            container_t* const leaf_source_particle_container =
                leaf->getParticleContainer();
351
            container_t* const leaf_target_particle_container =
352 353
                leaf->getParticleContainer();

354
            // Skip empty leaves
355 356
            if( leaf_source_particle_container->size() == 0
                && leaf_target_particle_container->size() == 0) {
357
                continue;
358 359
            }

360 361 362
            u_item_source_particle_containers.clear();
            u_item_indices.clear();

363
            for(node_t* u_item : leaf->U) {
364 365 366 367 368
                // The kernels do not consider leaf to be adjacent to
                // itself. Skip empty u_items
                if(u_item == leaf
                   || u_item->getParticleContainer()->size() == 0)
                {
369 370
                    continue;
                }
371
                u_item_source_particle_containers.push_back(u_item->getParticleContainer());
372
                u_item_indices.push_back(compute_box_offset_index(leaf, u_item, 1));
373
            }
374

375
            // Call P2P on vectors data
376 377 378
            _kernel.P2P(
                FTreeCoordinate(MortonIndex(leaf->getIndex())),
                leaf_target_particle_container,
379
                leaf_source_particle_container,
380 381
                u_item_source_particle_containers.data(),
                u_item_indices.data(),
382
                static_cast<int>(u_item_source_particle_containers.size())
383
                );
384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404
        }
    }


    /** 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
405 406 407
     *
     * \warning If the nodes are not at the same level, the lowest one's
     * ancestor at the highest one's level is used.
408 409
     */
    int compute_box_offset_index(node_t* node, node_t* other_node, const std::size_t n) {
410 411 412 413 414 415 416
        while(node->getDepth() > other_node->getDepth()) {
            node = node->getParent();
        }
        while(node->getDepth() < other_node->getDepth()) {
            other_node = other_node->getParent();
        }

417
        typename node_t::FReal box_width = node->getBox().width(0);
418 419 420 421

        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) {
422
            other_node_index = other_node_index * (2 * n + 1)
423 424 425 426 427 428 429 430 431 432 433
                + static_cast<unsigned long>(std::lround(center_offset[i] / box_width))
                + n;
        }

        return static_cast<int>(other_node_index);
    }


};

#endif