particles_system_builder.hpp 18.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
/******************************************************************************
*                                                                             *
*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
*                                                                             *
*  This file is part of bfps.                                                 *
*                                                                             *
*  bfps is free software: you can redistribute it and/or modify               *
*  it under the terms of the GNU General Public License as published          *
*  by the Free Software Foundation, either version 3 of the License,          *
*  or (at your option) any later version.                                     *
*                                                                             *
*  bfps is distributed in the hope that it will be useful,                    *
*  but WITHOUT ANY WARRANTY; without even the implied warranty of             *
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
*  GNU General Public License for more details.                               *
*                                                                             *
*  You should have received a copy of the GNU General Public License          *
*  along with bfps.  If not, see <http://www.gnu.org/licenses/>               *
*                                                                             *
* Contact: Cristian.Lalescu@ds.mpg.de                                         *
*                                                                             *
******************************************************************************/



26
27
28
29
30
#ifndef PARTICLES_SYSTEM_BUILDER_HPP
#define PARTICLES_SYSTEM_BUILDER_HPP

#include <string>

Cristian Lalescu's avatar
Cristian Lalescu committed
31
#include <cmath>
32
33
34
#include "abstract_particles_system.hpp"
#include "particles_system.hpp"
#include "particles_input_hdf5.hpp"
35
#include "particles_generic_interp.hpp"
36
#include "p2p_computer_empty.hpp"
37
#include "particles_inner_computer_empty.hpp"
38
39
40
41
42
43
44
45
46
47

#include "field.hpp"
#include "kspace.hpp"



//////////////////////////////////////////////////////////////////////////////
///
/// Double template "for"
///
Cristian Lalescu's avatar
Cristian Lalescu committed
48
49
50
51
52
53
54
/// "Template_double_for_if" is used to systematically generate specialized
/// templates for two ranges of template parameters.
/// For the interpolation we have `n` and `m` that designate "number of
/// neighbours" and "smoothness of interpolant".
/// Calling Template_double_for_if is essentially equivalent to having a couple
/// of nested "switch" statements, but without all the repetitive code lines.
///
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
//////////////////////////////////////////////////////////////////////////////

namespace Template_double_for_if{

template <class RetType,
          class IterType1, IterType1 CurrentIter1,
          class IterType2, const IterType2 CurrentIter2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, bool IsNotOver, typename... Args>
struct For2{
    static RetType evaluate(IterType2 value2, Args... args){
        if(CurrentIter2 == value2){
            return std::move(Func::template instanciate<CurrentIter1, CurrentIter2>(args...));
        }
        else{
            return std::move(For2<RetType,
                                        IterType1, CurrentIter1,
                                        IterType2, CurrentIter2+IterStep2, iterTo2, IterStep2,
                                        Func, (CurrentIter2+IterStep2 < iterTo2), Args...>::evaluate(value2, args...));
        }
    }
};

template <class RetType,
          class IterType1, IterType1 CurrentIter1,
          class IterType2, const IterType2 CurrentIter2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, typename... Args>
struct For2<RetType,
                  IterType1, CurrentIter1,
                  IterType2, CurrentIter2, iterTo2, IterStep2,
                  Func, false, Args...>{
    static RetType evaluate(IterType2 value2, Args... args){
86
        std::cout << __FUNCTION__ << "[ERROR] template values for loop 2 " << value2 << " does not exist\n";
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
        return RetType();
    }
};

template <class RetType,
          class IterType1, const IterType1 CurrentIter1, const IterType1 iterTo1, const IterType1 IterStep1,
          class IterType2, const IterType2 IterFrom2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, bool IsNotOver, typename... Args>
struct For1{
    static RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
        if(CurrentIter1 == value1){
            return std::move(For2<RetType,
                                        IterType1, CurrentIter1,
                                        IterType2, IterFrom2, iterTo2, IterStep2,
                                        Func, (IterFrom2<iterTo2), Args...>::evaluate(value2, args...));
        }
        else{
            return std::move(For1<RetType,
                              IterType1, CurrentIter1+IterStep1, iterTo1, IterStep1,
                              IterType2, IterFrom2, iterTo2, IterStep2,
                              Func, (CurrentIter1+IterStep1 < iterTo1), Args...>::evaluate(value1, value2, args...));
        }
    }
};

template <class RetType,
          class IterType1, const IterType1 IterFrom1, const IterType1 iterTo1, const IterType1 IterStep1,
          class IterType2, const IterType2 IterFrom2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, typename... Args>
struct For1<RetType,
                IterType1, IterFrom1, iterTo1, IterStep1,
                IterType2, IterFrom2, iterTo2, IterStep2,
                Func, false, Args...>{
    static RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
121
        std::cout << __FUNCTION__ << "[ERROR] template values for loop 1 " << value1 << " does not exist\n";
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
        return RetType();
    }
};

template <class RetType,
          class IterType1, const IterType1 IterFrom1, const IterType1 iterTo1, const IterType1 IterStep1,
          class IterType2, const IterType2 IterFrom2, const IterType2 iterTo2, const IterType2 IterStep2,
          class Func, typename... Args>
inline RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
    return std::move(For1<RetType,
            IterType1, IterFrom1, iterTo1, IterStep1,
            IterType2, IterFrom2, iterTo2, IterStep2,
            Func, (IterFrom1<iterTo1), Args...>::evaluate(value1, value2, args...));
}

}


//////////////////////////////////////////////////////////////////////////////
///
/// Builder Functions
///
//////////////////////////////////////////////////////////////////////////////

146
147
template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber, class p2p_computer_class,
          class particles_inner_computer_class, int size_particle_positions, int size_particle_rhs>
148
149
struct particles_system_build_container {
    template <const int interpolation_size, const int spline_mode>
150
    static std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>> instanciate(
151
             const field<field_rnumber, be, fc>* fs_field, // (field object)
152
153
             const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
             const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
154
             const partsize_t nparticles, // to check coherency between parameters and hdf input file
155
             const std::string& fname_input, // particles input filename
156
            const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
157
             MPI_Comm mpi_comm,
158
159
160
161
            const int in_current_iteration,
            p2p_computer_class p2p_computer,
            particles_inner_computer_class inner_computer,
            const particles_rnumber cutoff = std::numeric_limits<particles_rnumber>::max()){
162

163
        // The size of the field grid (global size) all_size seems
164
        std::array<size_t,3> field_grid_dim;
165
166
167
        field_grid_dim[IDXC_X] = fs_field->rlayout->sizes[IDXV_X];// nx
        field_grid_dim[IDXC_Y] = fs_field->rlayout->sizes[IDXV_Y];// nx
        field_grid_dim[IDXC_Z] = fs_field->rlayout->sizes[IDXV_Z];// nz
168
169
170

        // The size of the local field grid (the field nodes that belong to current process)
        std::array<size_t,3> local_field_dims;
171
172
173
        local_field_dims[IDXC_X] = fs_field->rlayout->subsizes[IDXV_X];
        local_field_dims[IDXC_Y] = fs_field->rlayout->subsizes[IDXV_Y];
        local_field_dims[IDXC_Z] = fs_field->rlayout->subsizes[IDXV_Z];
174
175
176

        // The offset of the local field grid
        std::array<size_t,3> local_field_offset;
177
178
179
        local_field_offset[IDXC_X] = fs_field->rlayout->starts[IDXV_X];
        local_field_offset[IDXC_Y] = fs_field->rlayout->starts[IDXV_Y];
        local_field_offset[IDXC_Z] = fs_field->rlayout->starts[IDXV_Z];
180
181
182
183
184
185
186


        // Retreive split from fftw to know processes that have no work
        int my_rank, nb_processes;
        AssertMpi(MPI_Comm_rank(mpi_comm, &my_rank));
        AssertMpi(MPI_Comm_size(mpi_comm, &nb_processes));

187
188
        const int split_step = (int(field_grid_dim[IDXC_Z])+nb_processes-1)/nb_processes;
        const int nb_processes_involved = (int(field_grid_dim[IDXC_Z])+split_step-1)/split_step;
189

190
191
192
        assert((my_rank < nb_processes_involved && local_field_dims[IDXC_Z] != 0)
               || (nb_processes_involved <= my_rank && local_field_dims[IDXC_Z] == 0));
        assert(nb_processes_involved <= int(field_grid_dim[IDXC_Z]));
193
194
195

        // Make the idle processes starting from the limit (and not 0 as set by fftw)
        if(nb_processes_involved <= my_rank){
196
            local_field_offset[IDXC_Z] = field_grid_dim[IDXC_Z];
197
198
        }

199
200
        // Ensure that 1D partitioning is used
        {
201
202
203
204
205
206
207
208
209
            assert(local_field_offset[IDXC_X] == 0);
            assert(local_field_offset[IDXC_Y] == 0);
            assert(local_field_dims[IDXC_X] == field_grid_dim[IDXC_X]);
            assert(local_field_dims[IDXC_Y] == field_grid_dim[IDXC_Y]);

            assert(my_rank >= nb_processes_involved || ((my_rank == 0 && local_field_offset[IDXC_Z] == 0)
                   || (my_rank != 0 && local_field_offset[IDXC_Z] != 0)));
            assert(my_rank >= nb_processes_involved || ((my_rank == nb_processes_involved-1 && local_field_offset[IDXC_Z]+local_field_dims[IDXC_Z] == field_grid_dim[IDXC_Z])
                   || (my_rank != nb_processes_involved-1 && local_field_offset[IDXC_Z]+local_field_dims[IDXC_Z] != field_grid_dim[IDXC_Z])));
210
        }
211

212
        // The spatial box size (all particles should be included inside)
213
        std::array<particles_rnumber,3> spatial_box_width;
214
215
216
        spatial_box_width[IDXC_X] = 4 * acos(0) / (fs_kk->dkx);
        spatial_box_width[IDXC_Y] = 4 * acos(0) / (fs_kk->dky);
        spatial_box_width[IDXC_Z] = 4 * acos(0) / (fs_kk->dkz);
217

218
219
        // Box is in the corner
        std::array<particles_rnumber,3> spatial_box_offset;
220
221
222
        spatial_box_offset[IDXC_X] = 0;
        spatial_box_offset[IDXC_Y] = 0;
        spatial_box_offset[IDXC_Z] = 0;
223

224
        // The distance between two field nodes in z
225
        std::array<particles_rnumber,3> spatial_partition_width;
226
227
228
        spatial_partition_width[IDXC_X] = spatial_box_width[IDXC_X]/particles_rnumber(field_grid_dim[IDXC_X]);
        spatial_partition_width[IDXC_Y] = spatial_box_width[IDXC_Y]/particles_rnumber(field_grid_dim[IDXC_Y]);
        spatial_partition_width[IDXC_Z] = spatial_box_width[IDXC_Z]/particles_rnumber(field_grid_dim[IDXC_Z]);
229
        // The spatial interval of the current process
230
231
        const particles_rnumber my_spatial_low_limit_z = particles_rnumber(local_field_offset[IDXC_Z])*spatial_partition_width[IDXC_Z];
        const particles_rnumber my_spatial_up_limit_z = particles_rnumber(local_field_offset[IDXC_Z]+local_field_dims[IDXC_Z])*spatial_partition_width[IDXC_Z];
232
233

        // Create the particles system
234
235
236
        using particles_system_type = particles_system<partsize_t, particles_rnumber, field_rnumber,
                                                       field<field_rnumber, be, fc>,
                                                       particles_generic_interp<particles_rnumber, interpolation_size,spline_mode>,
237
238
239
240
                                                       interpolation_size,
                                                       size_particle_positions, size_particle_rhs,
                                                       p2p_computer_class,
                                                       particles_inner_computer_class>;
241
242
243
244
245
246
247
248
249
250
        particles_system_type* part_sys = new particles_system_type(field_grid_dim,
                                               spatial_box_width,
                                               spatial_box_offset,
                                               spatial_partition_width,
                                               my_spatial_low_limit_z,
                                               my_spatial_up_limit_z,
                                               local_field_dims,
                                               local_field_offset,
                                               (*fs_field),
                                               mpi_comm,
251
                                               nparticles,
252
253
254
                                               cutoff,
                                               p2p_computer,
                                               inner_computer,
255
                                               in_current_iteration);
256

257
        // TODO P2P load particle data too
258
        // Load particles from hdf5
259
        particles_input_hdf5<partsize_t, particles_rnumber, size_particle_positions, size_particle_rhs> generator(mpi_comm, fname_input,
260
                                            inDatanameState, inDatanameRhs, my_spatial_low_limit_z, my_spatial_up_limit_z);
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275

        // Ensure parameters match the input file
        if(generator.getNbRhs() != nsteps){
            std::runtime_error(std::string("Nb steps is ") + std::to_string(nsteps)
                               + " in the parameters but " + std::to_string(generator.getNbRhs()) + " in the particles file.");
        }
        // Ensure parameters match the input file
        if(generator.getTotalNbParticles() != nparticles){
            std::runtime_error(std::string("Nb particles is ") + std::to_string(nparticles)
                               + " in the parameters but " + std::to_string(generator.getTotalNbParticles()) + " in the particles file.");
        }

        // Load the particles and move them to the particles system
        part_sys->init(generator);

Berenger Bramas's avatar
Berenger Bramas committed
276
277
        assert(part_sys->getNbRhs() == nsteps);

278
279
280
        // store particle file layout
        part_sys->setParticleFileLayout(generator.getParticleFileLayout());

281
        // Return the created particles system
282
        return std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>>(part_sys);
283
284
285
286
    }
};


287
template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber = double>
288
inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>> particles_system_builder(
289
        const field<field_rnumber, be, fc>* fs_field, // (field object)
290
291
        const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
        const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
292
        const partsize_t nparticles, // to check coherency between parameters and hdf input file
293
        const std::string& fname_input, // particles input filename
294
        const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
295
296
        const int interpolation_size,
        const int spline_mode,
297
298
        MPI_Comm mpi_comm,
        const int in_current_iteration){
299
    return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
300
                       int, 1, 11, 1, // interpolation_size
301
                       int, 0, 5, 1, // spline_mode
302
303
304
305
                       particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
                                                        p2p_computer_empty<particles_rnumber,partsize_t>,
                                                        particles_inner_computer_empty<particles_rnumber,partsize_t>,
                                                        3,3>>(
306
307
                           interpolation_size, // template iterator 1
                           spline_mode, // template iterator 2
308
309
                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration,
                           p2p_computer_empty<particles_rnumber,partsize_t>(), particles_inner_computer_empty<particles_rnumber,partsize_t>());
310
311
}

312
313
314
template <class partsize_t, class field_rnumber, field_backend be, field_components fc,
          class p2p_computer_class, class particles_inner_computer_class,
          class particles_rnumber = double>
315
inline std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>> particles_system_builder_with_p2p(
316
317
318
319
320
321
322
323
324
        const field<field_rnumber, be, fc>* fs_field, // (field object)
        const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
        const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
        const partsize_t nparticles, // to check coherency between parameters and hdf input file
        const std::string& fname_input, // particles input filename
        const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
        const int interpolation_size,
        const int spline_mode,
        MPI_Comm mpi_comm,
325
326
327
328
        const int in_current_iteration,
        p2p_computer_class p2p_computer,
        particles_inner_computer_class inner_computer,
        const particles_rnumber cutoff){
329
    return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>>,
330
                       int, 1, 11, 1, // interpolation_size
331
                       int, 0, 5, 1, // spline_mode
332
333
334
335
                       particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
                                                        p2p_computer_class,
                                                        particles_inner_computer_class,
                                                        6,6>>(
336
337
                           interpolation_size, // template iterator 1
                           spline_mode, // template iterator 2
338
339
                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration,
                           std::move(p2p_computer), std::move(inner_computer), cutoff);
340
341
342
343
}


#endif