MAJ terminée. Nous sommes passés en version 14.6.2 . Pour consulter les "releases notes" associées c'est ici :

https://about.gitlab.com/releases/2022/01/11/security-release-gitlab-14-6-2-released/
https://about.gitlab.com/releases/2022/01/04/gitlab-14-6-1-released/

FAdaptiveSequential.hpp 12.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
17
#include "Kernels/FKernelConcepts.hpp"

template<class _Tree, class _Kernel,
18
19
20
21
22
23
24
25
26
         typename std::enable_if<scalfmm::sfinae::has_P2M<_Tree, _Kernel>::value>::type* = nullptr,
         typename std::enable_if<scalfmm::sfinae::has_M2M<_Tree, _Kernel>::value>::type* = nullptr,
         typename std::enable_if<scalfmm::sfinae::has_M2L<_Tree, _Kernel>::value>::type* = nullptr,
         typename std::enable_if<scalfmm::sfinae::has_L2L<_Tree, _Kernel>::value>::type* = nullptr,
         typename std::enable_if<scalfmm::sfinae::has_L2P<_Tree, _Kernel>::value>::type* = nullptr,
         typename std::enable_if<scalfmm::sfinae::has_P2P<_Tree, _Kernel>::value>::type* = nullptr,
         typename std::enable_if<scalfmm::sfinae::has_M2P<_Tree, _Kernel>::value>::type* = nullptr,
         typename std::enable_if<scalfmm::sfinae::has_P2L<_Tree, _Kernel>::value>::type* = nullptr
         >
27
class FAdaptiveSequential : public FAlgorithmInterface, public FAlgorithmTimers {
28
29
30
public:
    using tree_t = _Tree;
    using kernel_t = _Kernel;
31
    using FReal = typename tree_t::FReal;
32
33
34
35
36
37
38
39
40
41
42
43
44

    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*;


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
59
60
61
62
63
    void setup_kernel(...) {}

public:

    FAdaptiveSequential(tree_t& tree, kernel_t& kernel) :
        _tree(tree),
        _kernel(kernel) {
    }

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
156
157
            if(leaf->getParticleContainer()->size() != 0) {
                _kernel.P2M(leaf->getData(), leaf->getParticleContainer());
            }
158
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()) {

                // Fill array of children data
167
                typename node_t::data_t* child_data[node_t::child_count] = {};
168
169
170
171
172
173
174
175
176
177
178
                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() {
179
180
181
        std::vector<decltype(std::declval<const node_t>().getData())> v_item_data_list;
        std::vector<int> v_item_indices;

182
        for(node_t& node : _tree.in_order_walk()) {
183
184
185
186
            if(node.is_leaf() && node.getParticleContainer()->size() == 0) {
                continue;
            }

187
188
189
            v_item_data_list.clear();
            v_item_indices.clear();

190
191
            // Needed to compute offset between boxes
            for(node_t* v_item : node.V) {
192
193
194
195
                if(v_item->is_leaf()
                   && v_item->getParticleContainer()->size() == 0) {
                    continue;
                }
196
197
                v_item_indices.push_back(compute_box_offset_index(&node, v_item, 3));
                v_item_data_list.push_back(v_item->getData());
198
            }
199
200
201

            // Call kernel M2L operator
            _kernel.M2L(node.getData(), v_item_data_list.data(), v_item_indices.data(), static_cast<int>(v_item_data_list.size()), static_cast<int>(node.getDepth()));
202
203
204
205
206
        }
    }

    // P2L
    void x_list_step() {
207
208
209
210
211
212
213
214
        /* 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.
         */
        for(node_t* node : _tree.leaves()) {
            if(node->getParticleContainer()->size() > 0) {
                for(node_t* w_item : node->W) {
215
216
217
                    if(w_item->is_leaf() && w_item->getParticleContainer()->size() == 0) {
                        continue;
                    }
218
219
                    _kernel.P2L(w_item->getData(), node->getParticleContainer());
                }
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
            }
        }
    }

    // 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() {
239
240
241
        for(node_t* node : _tree.leaves()) {
            if(node->getParticleContainer()->size() > 0) {
                for(node_t* w_item : node->W) {
242
243
244
                    if(w_item->is_leaf() && w_item->getParticleContainer()->size() == 0) {
                        continue;
                    }
245
246
                    _kernel.M2P(w_item->getData(), node->getParticleContainer());
                }
247
248
249
250
251
252
253
            }
        }
    }

    // L2P
    void down_to_target() {
        for(node_t* leaf : _tree.leaves()) {
254
255
256
            if(leaf->getParticleContainer()->size() != 0) {
                _kernel.L2P(leaf->getData(), leaf->getParticleContainer());
            }
257
258
259
        }
    }

260

261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    /** \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.
     *
     */
287
288
    void u_list_step() {
        using container_t = typename node_t::particle_container_t;
289

290
        // Vectors to pass to kernel P2P
291
292
        std::vector<container_t*> u_item_source_particle_containers;
        std::vector<int> u_item_indices;
293

294
        for(node_t* leaf : _tree.leaves()) {
295

296
297
            container_t* const leaf_source_particle_container =
                leaf->getParticleContainer();
298
            container_t* const leaf_target_particle_container =
299
300
                leaf->getParticleContainer();

301
            // Skip empty leaves
302
303
            if( leaf_source_particle_container->size() == 0
                && leaf_target_particle_container->size() == 0) {
304
                continue;
305
306
            }

307
308
309
            u_item_source_particle_containers.clear();
            u_item_indices.clear();

310
            for(node_t* u_item : leaf->U) {
311
312
313
314
315
                // The kernels do not consider leaf to be adjacent to
                // itself. Skip empty u_items
                if(u_item == leaf
                   || u_item->getParticleContainer()->size() == 0)
                {
316
317
                    continue;
                }
318
                u_item_source_particle_containers.push_back(u_item->getParticleContainer());
319
                u_item_indices.push_back(compute_box_offset_index(leaf, u_item, 1));
320
            }
321

322
            // Call P2P on vectors data
323
324
325
            _kernel.P2P(
                FTreeCoordinate(MortonIndex(leaf->getIndex())),
                leaf_target_particle_container,
326
                leaf_source_particle_container,
327
328
                u_item_source_particle_containers.data(),
                u_item_indices.data(),
329
                static_cast<int>(u_item_source_particle_containers.size())
330
                );
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
        }
    }


    /** 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
352
353
354
     *
     * \warning If the nodes are not at the same level, the lowest one's
     * ancestor at the highest one's level is used.
355
356
     */
    int compute_box_offset_index(node_t* node, node_t* other_node, const std::size_t n) {
357
358
359
360
361
362
363
        while(node->getDepth() > other_node->getDepth()) {
            node = node->getParent();
        }
        while(node->getDepth() < other_node->getDepth()) {
            other_node = other_node->getParent();
        }

364
365
366
367
368
        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) {
369
            other_node_index = other_node_index * (2 * n + 1)
370
371
372
373
374
375
376
377
378
379
380
                + static_cast<unsigned long>(std::lround(center_offset[i] / box_width))
                + n;
        }

        return static_cast<int>(other_node_index);
    }


};

#endif