1 /*============================================================================
2 * Extract information from lagrangian particles.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <stdarg.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *----------------------------------------------------------------------------*/
42
43 #include "bft_mem.h"
44 #include "bft_printf.h"
45
46 #include "cs_base.h"
47 #include "cs_field.h"
48 #include "cs_log.h"
49 #include "cs_mesh.h"
50 #include "cs_parall.h"
51 #include "cs_lagr_tracking.h"
52 #include "cs_selector.h"
53 #include "cs_timer.h"
54 #include "cs_time_step.h"
55
56 /*----------------------------------------------------------------------------
57 * Header for the current file
58 *----------------------------------------------------------------------------*/
59
60 #include "cs_lagr_extract.h"
61
62 /*----------------------------------------------------------------------------*/
63
64 BEGIN_C_DECLS
65
66 /*=============================================================================
67 * Additional Doxygen documentation
68 *============================================================================*/
69
70 /*!
71 \file cs_lagr_extract.c
72
73 \brief Extract information from lagrangian particles.
74 */
75
76 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
77
78 /*============================================================================
79 * Type definitions
80 *============================================================================*/
81
82 /*============================================================================
83 * Static global variables
84 *============================================================================*/
85
86 /*============================================================================
87 * Private function definitions
88 *============================================================================*/
89
90 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
91
92 /*============================================================================
93 * Public function definitions
94 *============================================================================*/
95
96 /*----------------------------------------------------------------------------*/
97 /*!
98 * \brief Get the local number of particles.
99 *
100 * \return current number of particles.
101 */
102 /*----------------------------------------------------------------------------*/
103
104 cs_lnum_t
cs_lagr_get_n_particles(void)105 cs_lagr_get_n_particles(void)
106 {
107 cs_lnum_t retval = 0;
108
109 const cs_lagr_particle_set_t *p_set = cs_glob_lagr_particle_set;
110 if (p_set != NULL)
111 retval = p_set->n_particles;
112
113 return retval;
114 }
115
116 /*----------------------------------------------------------------------------*/
117 /*!
118 * \brief Extract a list of particles using an optional cell filter and
119 * statistical density filter.
120 *
121 * The output array must have been allocated by the caller and be of
122 * sufficient size.
123 *
124 * \param[in] n_cells number of cells in filter
125 * \param[in] cell_list optional list of containing cells filter
126 * (1 to n numbering)
127 * \param[in] density if < 1, fraction of particles to select
128 * \param[out] n_particles number of selected particles, or NULL
129 * \param[out] particle_list particle_list (1 to n numbering), or NULL
130 */
131 /*----------------------------------------------------------------------------*/
132
133 void
cs_lagr_get_particle_list(cs_lnum_t n_cells,const cs_lnum_t cell_list[],double density,cs_lnum_t * n_particles,cs_lnum_t * particle_list)134 cs_lagr_get_particle_list(cs_lnum_t n_cells,
135 const cs_lnum_t cell_list[],
136 double density,
137 cs_lnum_t *n_particles,
138 cs_lnum_t *particle_list)
139 {
140 cs_lnum_t i;
141
142 ptrdiff_t displ = 0;
143
144 cs_lnum_t p_count = 0;
145
146 bool *cell_flag = NULL;
147
148 const cs_mesh_t *mesh = cs_glob_mesh;
149
150 const cs_lagr_particle_set_t *p_set = cs_glob_lagr_particle_set;
151
152 assert(p_set != NULL);
153
154 size_t extents = p_set->p_am->extents;
155
156 if (density < 1) {
157
158 size_t _extents, size;
159 cs_datatype_t datatype;
160 int count;
161
162 cs_lagr_get_attr_info(p_set,
163 0,
164 CS_LAGR_RANDOM_VALUE,
165 &_extents, &size, &displ,
166 &datatype, &count);
167
168 assert( (displ > 0 && count == 1 && datatype == CS_REAL_TYPE)
169 || (displ < 0 && count == 0 && datatype == CS_DATATYPE_NULL));
170
171 }
172
173 /* Case where we have a filter */
174
175 if (n_cells < mesh->n_cells) {
176
177 /* Build cel filter */
178
179 BFT_MALLOC(cell_flag, mesh->n_cells, bool);
180
181 for (i = 0; i < mesh->n_cells; i++)
182 cell_flag[i] = false;
183
184 if (cell_list != NULL) {
185 for (i = 0; i < n_cells; i++)
186 cell_flag[cell_list[i] - 1] = true;
187 }
188 else {
189 for (i = 0; i < n_cells; i++)
190 cell_flag[i] = true;
191 }
192
193 }
194
195 /* Loop on particles */
196
197 for (i = 0; i < p_set->n_particles; i++) {
198
199 /* If density < 1, randomly select which particles are added;
200 normally, particles maintain a random_value for this purpose,
201 but we also plan for the case where this could be optional. */
202
203 if (density < 1) {
204 double r;
205 if (displ < 0)
206 r = (double)rand() / RAND_MAX;
207 else {
208 const unsigned char
209 *p = (const unsigned char *)(p_set->p_buffer + i*extents) + displ;
210 r = *(const cs_real_t *)p;
211 }
212 if (r > density)
213 continue;
214 }
215
216 /* Check for filter cell */
217
218 if (cell_flag != NULL) {
219 cs_lnum_t cur_cell_id
220 = cs_lagr_particles_get_lnum(p_set, i, CS_LAGR_CELL_ID);
221 if (cell_flag[cur_cell_id] == false)
222 continue;
223 }
224
225 if (particle_list != NULL)
226 particle_list[p_count] = i+1;
227
228 p_count += 1;
229
230 }
231
232 if (cell_flag != NULL)
233 BFT_FREE(cell_flag);
234
235 *n_particles = p_count;
236 }
237
238 /*----------------------------------------------------------------------------*/
239 /*!
240 * \brief Extract values for a set of particles.
241 *
242 * The output array must have been allocated by the caller and be of
243 * sufficient size.
244 *
245 * \param[in] particles associated particle set
246 * \param[in] attr attribute whose values are required
247 * \param[in] datatype associated value type
248 * \param[in] stride number of values per particle
249 * \param[in] component_id if -1 : extract the whole attribute
250 * if >0 : id of the component to extract
251 * \param[in] n_particles number of particles in filter
252 * \param[in] particle_list particle_list (1 to n numbering), or NULL
253 * \param[out] values particle values for given attribute
254 *
255 * \return 0 in case of success, 1 if attribute is not present
256 */
257 /*----------------------------------------------------------------------------*/
258
259 int
cs_lagr_get_particle_values(const cs_lagr_particle_set_t * particles,cs_lagr_attribute_t attr,cs_datatype_t datatype,int stride,int component_id,cs_lnum_t n_particles,const cs_lnum_t particle_list[],void * values)260 cs_lagr_get_particle_values(const cs_lagr_particle_set_t *particles,
261 cs_lagr_attribute_t attr,
262 cs_datatype_t datatype,
263 int stride,
264 int component_id,
265 cs_lnum_t n_particles,
266 const cs_lnum_t particle_list[],
267 void *values)
268 {
269 size_t j;
270 cs_lnum_t i;
271
272 size_t extents, size, _length;
273 ptrdiff_t displ;
274 cs_datatype_t _datatype;
275 int _count;
276 unsigned char *_values = values;
277
278 assert(particles != NULL);
279
280 cs_lagr_get_attr_info(particles, 0, attr,
281 &extents, &size, &displ, &_datatype, &_count);
282
283 if (_count == 0)
284 return 1;
285 else {
286 if (component_id == -1)
287 _length = size;
288 else
289 _length = size/_count;
290 }
291
292 /* Check consistency */
293
294 if (cs_lagr_check_attr_query(particles,
295 attr,
296 datatype,
297 stride,
298 component_id) != 0)
299 return 1;
300
301 /* No offset if export of the whole attribute */
302 if (component_id == -1)
303 component_id = 0;
304
305 /* Case where we have no filter */
306
307 if (particle_list == NULL) {
308 for (i = 0; i < n_particles; i++) {
309 unsigned char *dest = _values + i*_length;
310 const unsigned char
311 *src = (const unsigned char *)(particles->p_buffer + i*extents)
312 + displ
313 + component_id * _length;
314 for (j = 0; j < _length; j++)
315 dest[j] = src[j];
316 }
317 }
318
319 /* Case where we have a filter list */
320 else {
321 for (i = 0; i < n_particles; i++) {
322 cs_lnum_t p_id = particle_list[i] - 1;
323 unsigned char *dest = _values + i*_length;
324 const unsigned char
325 *src = (const unsigned char *)(particles->p_buffer + p_id*extents)
326 + displ
327 + component_id * _length;
328 for (j = 0; j < _length; j++)
329 dest[j] = src[j];
330 }
331 }
332
333 return 0;
334 }
335
336 /*----------------------------------------------------------------------------*/
337 /*!
338 * \brief Extract trajectory values for a set of particles.
339 *
340 * Trajectories are defined as a mesh of segments, whose start and end
341 * points are copied in an interleaved manner in the segment_values array
342 * (p1_old, p1_new, p2_old, p2_new, ... pn_old, pn_new).
343 *
344 * The output array must have been allocated by the caller and be of
345 * sufficient size.
346 *
347 * \param[in] particles associated particle set
348 * \param[in] attr attribute whose values are required
349 * \param[in] datatype associated value type
350 * \param[in] stride number of values per particle
351 * \param[in] component_id if -1 : extract the whole attribute
352 * if >0 : id of the component to extract
353 * \param[in] n_particles number of particles in filter
354 * \param[in] particle_list particle_list (1 to n numbering), or NULL
355 * \param[out] segment_values particle segment values
356 *
357 * \return 0 in case of success, 1 if attribute is not present
358 */
359 /*----------------------------------------------------------------------------*/
360
361 int
cs_lagr_get_trajectory_values(const cs_lagr_particle_set_t * particles,cs_lagr_attribute_t attr,cs_datatype_t datatype,int stride,int component_id,cs_lnum_t n_particles,const cs_lnum_t particle_list[],void * segment_values)362 cs_lagr_get_trajectory_values(const cs_lagr_particle_set_t *particles,
363 cs_lagr_attribute_t attr,
364 cs_datatype_t datatype,
365 int stride,
366 int component_id,
367 cs_lnum_t n_particles,
368 const cs_lnum_t particle_list[],
369 void *segment_values)
370 {
371 size_t j;
372 cs_lnum_t i;
373
374 size_t extents, size, _length;
375 ptrdiff_t displ, displ_p;
376 cs_datatype_t _datatype;
377 int _count;
378 unsigned char *_values = segment_values;
379
380 assert(particles != NULL);
381
382 const unsigned char *p_buffer = particles->p_buffer;
383
384 cs_lagr_get_attr_info(particles, 0, attr,
385 &extents, &size, &displ, &_datatype, &_count);
386
387 if (_count == 0)
388 return 1;
389 else {
390 if (component_id == -1)
391 _length = size;
392 else
393 _length = size/_count;
394 }
395
396 if (particles->p_am->count[1][attr] > 0)
397 cs_lagr_get_attr_info(particles, 1, attr,
398 &extents, NULL, &displ_p, NULL, NULL);
399
400 /* Check consistency */
401
402 if (cs_lagr_check_attr_query(particles,
403 attr,
404 datatype,
405 stride,
406 component_id) != 0)
407 return 1;
408
409 /* No offset if export of the whole attribute */
410 if (component_id == -1)
411 component_id = 0;
412
413 /* Case where we have no filter */
414
415 if (particle_list == NULL) {
416
417 if (particles->p_am->count[1][attr] > 0) {
418
419 for (i = 0; i < n_particles; i++) {
420 unsigned char *dest = _values + i*_length*2;
421 const unsigned char
422 *src = (const unsigned char *)(p_buffer + i*extents)
423 + displ
424 + component_id * _length;
425 const unsigned char
426 *srcp = (const unsigned char *)(p_buffer + i*extents)
427 + displ_p
428 + component_id * _length;
429 for (j = 0; j < _length; j++) {
430 dest[j] = src[j];
431 dest[j + _length] = srcp[j];
432 }
433
434 }
435
436 }
437 else { /* With no previous value available; copy current value */
438
439 for (i = 0; i < n_particles; i++) {
440 unsigned char *dest = _values + i*_length*2;
441 const unsigned char
442 *src = (const unsigned char *)(p_buffer + i*extents)
443 + displ
444 + component_id * _length;
445 for (j = 0; j < _length; j++) {
446 dest[j] = src[j];
447 dest[j + _length] = src[j];
448 }
449
450 }
451
452 }
453 }
454
455 /* Case where we have a filter list */
456 else {
457
458 if (particles->p_am->count[1][attr] > 0) {
459
460 for (i = 0; i < n_particles; i++) {
461 cs_lnum_t p_id = particle_list[i] - 1;
462 unsigned char *dest = _values + i*_length*2;
463 const unsigned char
464 *src = (const unsigned char *)(p_buffer + p_id*extents)
465 + displ
466 + component_id * _length;
467 const unsigned char
468 *srcp = (const unsigned char *)(p_buffer + p_id*extents)
469 + displ_p
470 + component_id * _length;
471 for (j = 0; j < _length; j++) {
472 dest[j] = src[j];
473 dest[j + _length] = srcp[j];
474 }
475 }
476
477 }
478
479 else { /* With no previous value available; copy current value */
480
481 for (i = 0; i < n_particles; i++) {
482 cs_lnum_t p_id = particle_list[i] - 1;
483 unsigned char *dest = _values + i*_length*2;
484 const unsigned char
485 *src = (const unsigned char *)(p_buffer + p_id*extents)
486 + displ
487 + component_id * _length;
488 for (j = 0; j < _length; j++) {
489 dest[j] = src[j];
490 dest[j + _length] = src[j];
491 }
492 }
493
494 }
495
496 }
497
498 return 0;
499 }
500
501 /*----------------------------------------------------------------------------*/
502
503 END_C_DECLS
504