abstract_particle_set.hpp 14.9 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
26
27
28
29
30
31
/******************************************************************************
*                                                                             *
*  Copyright 2020 Max Planck Institute for Dynamics and Self-Organization     *
*                                                                             *
*  This file is part of TurTLE.                                               *
*                                                                             *
*  TurTLE 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.                                     *
*                                                                             *
*  TurTLE 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 TurTLE.  If not, see <http://www.gnu.org/licenses/>             *
*                                                                             *
* Contact: Cristian.Lalescu@ds.mpg.de                                         *
*                                                                             *
******************************************************************************/



#ifndef ABSTRACT_PARTICLE_SET_HPP
#define ABSTRACT_PARTICLE_SET_HPP

#include <array>
#include <cassert>
#include <vector>
32
#include <memory>
33
#include "field.hpp"
34
#include "particles/p2p/p2p_distr_mpi.hpp"
35
#include "particles/particles_sampling.hpp"
36

37
38
39
40
41
42
43
44
45
46
47


/** Brings together particle information with interpolation functionality.
 *
 * This is an abstract class that defines the functionality required by the
 * particle solver code.
 * We define methods for:
 *  - accessing particle set data and miscelaneous information.
 *  - updating particle information and redistributing data among MPI processes
 *    accordingly.
 *  - approximating fields at particle locations ("sampling").
Cristian Lalescu's avatar
Cristian Lalescu committed
48
 *  - applying a generic particle-to-particle interaction object (P2PKernel)
49
50
51
 *
 */

52
53
class abstract_particle_set
{
54
    public:
55
56
57
        using partsize_t = long long int;
        using particle_rnumber = double;

58
        // virtual destructor
59
60
61
        virtual ~abstract_particle_set(){}

        // extract particle information
62
63
        virtual particle_rnumber* getParticleState() const = 0;
        virtual partsize_t* getParticleIndex() const = 0;
64
        virtual int setParticleState(particle_rnumber *) = 0;
Cristian Lalescu's avatar
Cristian Lalescu committed
65
        virtual int getParticleState(particle_rnumber *) = 0;
66
67
68
69
70
71
72

        virtual std::unique_ptr<particle_rnumber[]> extractFromParticleState(
                const int firstComponent,
                const int lastComponent) const = 0;

        virtual partsize_t getLocalNumberOfParticles() const = 0;
        virtual partsize_t getTotalNumberOfParticles() const = 0;
73
        virtual partsize_t* getParticlesPerPartition() const = 0;
74
        virtual int getStateSize() const = 0;
75

76
77
        virtual std::vector<hsize_t> getParticleFileLayout() = 0;

78
79
80
        // get p2p computer
        virtual p2p_distr_mpi<partsize_t, particle_rnumber>* getP2PComputer() = 0;

81
82
83
84
        // generic loop over data
        template <class func_type>
        int LOOP(func_type expression)
        {
Cristian Lalescu's avatar
Cristian Lalescu committed
85
            TIMEZONE("abstract_particle_set::LOOP");
86
87
88
89
90
91
92
            for (partsize_t idx = 0; idx < this->getLocalNumberOfParticles()*this->getStateSize(); idx++)
                expression(idx);
            return EXIT_SUCCESS;
        }
        template <class func_type>
        int LOOP_state(func_type expression)
        {
Cristian Lalescu's avatar
Cristian Lalescu committed
93
            TIMEZONE("abstract_particle_set::LOOP_state");
94
95
96
97
98
            for (partsize_t idx_part = 0; idx_part < this->getLocalNumberOfParticles(); idx_part++)
            for (unsigned int cc = 0; cc < this->getStateSize(); cc++)
                expression(idx_part, cc);
            return EXIT_SUCCESS;
        }
99

100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
        int copy_state_tofrom(
                particle_rnumber *dst,
                const particle_rnumber *src)
        {
            return this->LOOP(
                    [&](const partsize_t idx)
                    {
                        dst[idx] = src[idx];
                    });
        }

        int copy_state_tofrom(
                std::unique_ptr<particle_rnumber[]> &dst,
                const std::unique_ptr<particle_rnumber[]> &src)
        {
            return this->LOOP(
                    [&](const partsize_t idx)
                    {
                        dst[idx] = src[idx];
                    });
        }

122
123
124
125
126
127
128
        /** Reorganize particles within MPI domain.
         *
         * Based on current particle locations, redistribute the full state
         * data among the MPI processes.
         * Optional: also redistribute a list of arrays of the same shape.
         */
        virtual int redistribute(
129
                std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data) = 0;
130
131
132
133
134
135

        // sample field values at particle location
        virtual int sample(
                const field<float,
                            FFTW,
                            ONE>         &field_to_sample,
136
                particle_rnumber *result) = 0;
137
138
139
140
        virtual int sample(
                const field<float,
                            FFTW,
                            THREE>       &field_to_sample,
141
                particle_rnumber *result) = 0;
142
143
144
145
        virtual int sample(
                const field<float,
                            FFTW,
                            THREExTHREE> &field_to_sample,
146
                particle_rnumber *result) = 0;
147
148
149
150
        virtual int sample(
                const field<double,
                            FFTW,
                            ONE>         &field_to_sample,
151
                particle_rnumber *result) = 0;
152
153
154
155
        virtual int sample(
                const field<double,
                            FFTW,
                            THREE>       &field_to_sample,
156
                particle_rnumber *result) = 0;
157
158
159
160
        virtual int sample(
                const field<double,
                            FFTW,
                            THREExTHREE> &field_to_sample,
161
                particle_rnumber *result) = 0;
162
163

        // apply p2p kernel
Cristian Lalescu's avatar
Cristian Lalescu committed
164
165
166
167
        // IMPORTANT: this method shuffles particle arrays.
        // If you need some data to be shuffled according to the same permutation,
        // please place it in `additional_data`.
        // This applies for instance to rhs data during time steps.
168
169
170
        template <int state_size,
                  class p2p_kernel_class>
        int applyP2PKernel(
171
                p2p_kernel_class &p2p_kernel,
172
173
                std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data)
        {
Cristian Lalescu's avatar
Cristian Lalescu committed
174
            TIMEZONE("abstract_particle_set::applyP2PKernel");
175
176
177
178
179
180
181
182
183
184
185
186
187
188
            // there must be at least one array with additional data,
            // since the p2p kernel expects an array where to store the result of the computation
            assert(int(additional_data.size()) > int(0));
            this->getP2PComputer()->template compute_distr<p2p_kernel_class,
                                                           state_size, // I'd prefer to use this->getStateSize() here. not possible because virtual function is not constexpr...
                                                           state_size>(
                    p2p_kernel,
                    this->getParticlesPerPartition(),
                    this->getParticleState(),
                    additional_data.data(),
                    int(additional_data.size()),
                    this->getParticleIndex());
            return EXIT_SUCCESS;
        }
189
190
191
192

        template <typename field_rnumber,
                  field_backend be,
                  field_components fc>
193
        int writeSample(
194
195
196
197
198
199
                field<field_rnumber, be, fc> *field_to_sample,
                particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
                const std::string species_name,
                const std::string field_name,
                const int iteration)
        {
Cristian Lalescu's avatar
Cristian Lalescu committed
200
            TIMEZONE("abstract_particle_set::writeSample");
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
            // set file layout
            particle_sample_writer->setParticleFileLayout(this->getParticleFileLayout());
            // allocate position array
            std::unique_ptr<particle_rnumber[]> xx = this->extractFromParticleState(0, 3);
            // allocate temporary array
            std::unique_ptr<particle_rnumber[]> pdata(new particle_rnumber[ncomp(fc)*this->getLocalNumberOfParticles()]);
            // clean up temporary array
            std::fill_n(pdata.get(), ncomp(fc)*this->getLocalNumberOfParticles(), 0);
            this->sample(*field_to_sample, pdata.get());
            particle_sample_writer->template save_dataset<ncomp(fc)>(
                    species_name,
                    field_name,
                    xx.get(),
                    &pdata,
                    this->getParticleIndex(),
                    this->getLocalNumberOfParticles(),
                    iteration);
            // deallocate temporary array
            delete[] pdata.release();
            // deallocate position array
            delete[] xx.release();
            return EXIT_SUCCESS;
        }

225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
        template <typename field_rnumber,
                  field_backend be,
                  field_components fc>
        int writeSample(
                field<field_rnumber, be, fc> *field_to_sample,
                const std::string file_name,
                const std::string species_name,
                const std::string field_name,
                const int iteration)
        {
            TIMEZONE("abstract_particle_set::writeSample_direct_to_file");
            // allocate temporary array
            std::unique_ptr<particle_rnumber[]> pdata(new particle_rnumber[ncomp(fc)*this->getLocalNumberOfParticles()]);
            // clean up temporary array
            set_particle_data_to_zero<partsize_t, particle_rnumber, ncomp(fc)>(
                    pdata.get(), this->getLocalNumberOfParticles());
            this->sample(*field_to_sample, pdata.get());
            hid_t file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
            hdf5_tools::gather_and_write_with_single_rank(
                    field_to_sample->myrank,
                    0,
                    field_to_sample->comm,
                    pdata.get(),
                    this->getLocalNumberOfParticles(),
                    file_id,
                    species_name + std::string("/") + field_name + std::string("/") + std::to_string(iteration));
            H5Fclose(file_id);
            // deallocate temporary array
            delete[] pdata.release();
            return EXIT_SUCCESS;
        }

257
        int writeStateChunk(
258
                const int i0,
259
                const int contiguous_state_chunk,
260
261
262
263
264
                particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
                const std::string species_name,
                const std::string field_name,
                const int iteration)
        {
265
            TIMEZONE("abstract_particle_set::writeStateChunk");
266
            assert(i0 >= 0);
267
            assert(i0 <= this->getStateSize()-contiguous_state_chunk);
268
269
270
271
272
            // set file layout
            particle_sample_writer->setParticleFileLayout(this->getParticleFileLayout());
            // allocate position array
            std::unique_ptr<particle_rnumber[]> xx = this->extractFromParticleState(0, 3);
            // allocate temporary array
273
            std::unique_ptr<particle_rnumber[]> yy = this->extractFromParticleState(i0, i0+contiguous_state_chunk);
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
            switch(contiguous_state_chunk)
            {
                case 1:
                    particle_sample_writer->template save_dataset<1>(
                            species_name,
                            field_name,
                            xx.get(),
                            &yy,
                            this->getParticleIndex(),
                            this->getLocalNumberOfParticles(),
                            iteration);
                    break;
                case 3:
                    particle_sample_writer->template save_dataset<3>(
                            species_name,
                            field_name,
                            xx.get(),
                            &yy,
                            this->getParticleIndex(),
                            this->getLocalNumberOfParticles(),
                            iteration);
                    break;
                case 9:
                    particle_sample_writer->template save_dataset<9>(
                            species_name,
                            field_name,
                            xx.get(),
                            &yy,
                            this->getParticleIndex(),
                            this->getLocalNumberOfParticles(),
                            iteration);
                    break;
                default:
                    DEBUG_MSG("code not specialized for this value of contiguous_state_chunk in abstract_particle_set::writeStateChunk\n");
                    std::cerr << "error in abstract_particle_set::writeStateChunk. please contact maintainer.\n" << std::endl;
            }
310
311
312
313
314
315
            // deallocate temporary array
            delete[] yy.release();
            // deallocate position array
            delete[] xx.release();
            return EXIT_SUCCESS;
        }
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347

        int writeStateTriplet(
                const int i0,
                particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
                const std::string species_name,
                const std::string field_name,
                const int iteration)
        {
            TIMEZONE("abstract_particle_set::writeStateTriplet");
            return this->writeStateChunk(
                    i0, 3,
                    particle_sample_writer,
                    species_name,
                    field_name,
                    iteration);
        }

        int writeStateComponent(
                const int i0,
                particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
                const std::string species_name,
                const std::string field_name,
                const int iteration)
        {
            TIMEZONE("abstract_particle_set::writeStateComponent");
            return this->writeStateChunk(
                    i0, 1,
                    particle_sample_writer,
                    species_name,
                    field_name,
                    iteration);
        }
348
349
350
351
};

#endif//ABSTRACT_PARTICLE_SET_HPP