1 /*============================================================================
2 * Particle injection for lagrangian module.
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 <limits.h>
34 #include <stdio.h>
35 #include <stddef.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <math.h>
39 #include <ctype.h>
40 #include <float.h>
41 #include <assert.h>
42
43 /*----------------------------------------------------------------------------
44 * Local headers
45 *----------------------------------------------------------------------------*/
46
47 #include "bft_printf.h"
48 #include "bft_error.h"
49 #include "bft_mem.h"
50
51 #include "cs_base.h"
52 #include "cs_math.h"
53
54 #include "cs_boundary_zone.h"
55 #include "cs_volume_zone.h"
56
57 #include "cs_ht_convert.h"
58 #include "cs_mesh.h"
59 #include "cs_mesh_quantities.h"
60 #include "cs_parall.h"
61 #include "cs_thermal_model.h"
62 #include "cs_parameters.h"
63 #include "cs_physical_model.h"
64 #include "cs_physical_constants.h"
65 #include "cs_prototypes.h"
66 #include "cs_time_step.h"
67
68 #include "cs_field.h"
69 #include "cs_field_pointer.h"
70 #include "cs_random.h"
71
72 #include "cs_lagr.h"
73 #include "cs_lagr_gradients.h"
74 #include "cs_lagr_tracking.h"
75 #include "cs_lagr_new.h"
76 #include "cs_lagr_precipitation_model.h"
77 #include "cs_lagr_prototypes.h"
78
79 /*----------------------------------------------------------------------------
80 * Header for the current file
81 *----------------------------------------------------------------------------*/
82
83 #include "cs_lagr_injection.h"
84
85 /*----------------------------------------------------------------------------*/
86
87 BEGIN_C_DECLS
88
89 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
90
91 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
92
93 /*============================================================================
94 * Private function definitions
95 *============================================================================*/
96
97 /*----------------------------------------------------------------------------*/
98 /*!
99 * \brief Binary search for a given local id in a given array of
100 * ordered values.
101 *
102 * We assume the id is present in the array.
103 *
104 * \param[in] n number of values
105 * \param[in] x value to locate
106 * \param[in] a array of ordered values (size n)
107 *
108 * \return index of x in array (smallest i such that a[i] >= x)
109 */
110 /*----------------------------------------------------------------------------*/
111
112 static inline cs_lnum_t
_segment_binary_search(cs_lnum_t n,double x,const double a[])113 _segment_binary_search(cs_lnum_t n,
114 double x,
115 const double a[])
116 {
117 cs_lnum_t start_id = 0;
118 cs_lnum_t end_id = n-1;
119 cs_lnum_t mid_id = (end_id - start_id) / 2;
120
121 x = CS_MIN(x, a[end_id]); /* precaution: force in range */
122
123 while (start_id < end_id) {
124 if (a[mid_id] < x)
125 start_id = mid_id + 1;
126 else /* if (a[mid_id] >= x) */
127 end_id = mid_id;
128 mid_id = start_id + ((end_id - start_id) / 2);
129 }
130
131 assert(mid_id >= 0 && mid_id < n);
132
133 return mid_id;
134 }
135
136 /*----------------------------------------------------------------------------*/
137 /*!
138 * \brief Distribute new particles in a given region.
139 *
140 * \param[in] n_g_particles global number of particles to inject
141 * \param[in] n_elts number of elements in region
142 * \param[in] elt_id element ids (or NULL)
143 * \param[in] elt_weight parent element weights
144 * (i.e. all local surfaces or volumes)
145 * \param[in] elt_profile optional profile values for elements (or NULL)
146 * \param[out] elt_particle_idx start index of added particles for each
147 * element (size: n_elts + 1)
148 *
149 * \return number of particles added on local rank
150 */
151 /*----------------------------------------------------------------------------*/
152
153 static cs_lnum_t
_distribute_particles(cs_gnum_t n_g_particles,cs_lnum_t n_elts,const cs_lnum_t elt_id[],const cs_real_t elt_weight[],const cs_real_t * elt_profile,cs_lnum_t elt_particle_idx[])154 _distribute_particles(cs_gnum_t n_g_particles,
155 cs_lnum_t n_elts,
156 const cs_lnum_t elt_id[],
157 const cs_real_t elt_weight[],
158 const cs_real_t *elt_profile,
159 cs_lnum_t elt_particle_idx[])
160 {
161 cs_lnum_t n_particles = (cs_glob_n_ranks > 1) ? 0 : n_g_particles;
162
163 /* Compute local element weight */
164
165 cs_real_t *elt_cm_weight = NULL;
166
167 BFT_MALLOC(elt_cm_weight, n_elts, cs_real_t);
168
169 if (elt_id != NULL) {
170 if (elt_profile != NULL) {
171 for (cs_lnum_t i = 0; i < n_elts; i++)
172 elt_cm_weight[i] = elt_weight[elt_id[i]]*elt_profile[i];
173 }
174 else {
175 for (cs_lnum_t i = 0; i < n_elts; i++)
176 elt_cm_weight[i] = elt_weight[elt_id[i]];
177 }
178 }
179 else {
180 if (elt_profile != NULL) {
181 for (cs_lnum_t i = 0; i < n_elts; i++)
182 elt_cm_weight[i] = elt_weight[i]*elt_profile[i];
183 }
184 else {
185 for (cs_lnum_t i = 0; i < n_elts; i++)
186 elt_cm_weight[i] = elt_weight[i];
187 }
188 }
189
190 /* Transform to cumulative weight using Kahan summation */
191
192 double l_weight = 0;
193 {
194 double d = 0., c = 0.;
195 for (cs_lnum_t i = 0; i < n_elts; i++) {
196 double z = elt_cm_weight[i] - c;
197 double t = d + z;
198 c = (t - d) - z;
199 d = t;
200 elt_cm_weight[i] = d;
201 }
202 l_weight = d;
203 }
204
205 #if defined(HAVE_MPI)
206
207 /* Pre_distribution to various ranks; we assume that the number of
208 injected particles at a given time is not huge, so it is cheaper
209 to precompute the distribution on a single rank and broadcast it.
210 For a higher number of particles, computing by blocs and then
211 redistributing (with "all to all" operations) could be more efficient. */
212
213 if (cs_glob_n_ranks > 1) {
214
215 int n_ranks = cs_glob_n_ranks;
216 int l_rank = cs_glob_rank_id;
217 int r_rank = 0; /* Root rank for serialized operations */
218
219 cs_lnum_t *n_rank_particles = NULL;
220 double *cm_weight = NULL;
221
222 if (l_rank == r_rank) {
223
224 BFT_MALLOC(n_rank_particles, n_ranks, cs_lnum_t);
225 BFT_MALLOC(cm_weight, n_ranks, double);
226
227 for (int i = 0; i < n_ranks; i++)
228 n_rank_particles[i] = 0;
229
230 }
231
232 MPI_Gather(&l_weight, 1, MPI_DOUBLE, cm_weight, 1, MPI_DOUBLE,
233 r_rank, cs_glob_mpi_comm);
234
235 if (l_rank == r_rank) {
236
237 /* Scan (cumulative sum) operation */
238 for (int i = 1; i < n_ranks; i++)
239 cm_weight[i] += cm_weight[i-1];
240
241 /* Scale to [0, 1] */
242 double tot_weight = cm_weight[n_ranks-1];
243
244 if (tot_weight > 0.) {
245
246 for (int i = 0; i < n_ranks; i++)
247 cm_weight[i] /= tot_weight;
248
249 /* Compute distribution */
250
251 for (cs_gnum_t i = 0; i < n_g_particles; i++) {
252 cs_real_t r;
253 cs_random_uniform(1, &r);
254 int r_id = _segment_binary_search(n_ranks, r, cm_weight);
255 n_rank_particles[r_id] += 1;
256 }
257
258 }
259
260 BFT_FREE(cm_weight);
261 }
262
263 MPI_Scatter(n_rank_particles, 1, CS_MPI_LNUM,
264 &n_particles, 1, CS_MPI_LNUM,
265 r_rank, cs_glob_mpi_comm);
266
267 BFT_FREE(n_rank_particles);
268 }
269
270 #endif /* defined(HAVE_MPI) */
271
272 /* Check for empty zones */
273
274 if (n_particles > 0 && n_elts < 1)
275 n_particles = 0;
276
277 /* Now distribute locally */
278
279 for (cs_lnum_t i = 0; i < n_elts; i++)
280 elt_particle_idx[i] = 0;
281 elt_particle_idx[n_elts] = 0;
282
283 for (cs_lnum_t i = 0; i < n_elts; i++)
284 elt_cm_weight[i] /= l_weight;
285
286 /* Compute distribution */
287
288 for (cs_lnum_t i = 0; i < n_particles; i++) {
289 cs_real_t r;
290 cs_random_uniform(1, &r);
291 cs_lnum_t e_id = _segment_binary_search(n_elts, r, elt_cm_weight);
292 elt_particle_idx[e_id+1] += 1;
293 }
294
295 BFT_FREE(elt_cm_weight);
296
297 /* transform count to index */
298
299 for (cs_lnum_t i = 0; i < n_elts; i++)
300 elt_particle_idx[i+1] += elt_particle_idx[i];
301
302 assert(elt_particle_idx[n_elts] == n_particles);
303
304 return n_particles;
305 }
306
307 /*----------------------------------------------------------------------------*/
308 /*!
309 * \brief Check injection parameters are valid.
310 *
311 * \param[in] zis pointer to injection data for a given zone and set
312 */
313 /*----------------------------------------------------------------------------*/
314
315 static void
_injection_check(const cs_lagr_injection_set_t * zis)316 _injection_check(const cs_lagr_injection_set_t *zis)
317 {
318 const char _profile_err_fmt_i[]
319 = N_("Lagrangian %s zone %d, set %d\n"
320 " %s profile value (%d) is invalid.");
321 const char _profile_err_fmt_d[]
322 = N_("Lagrangian %s zone %d, set %d\n"
323 " %s profile value (%g) is invalid.");
324
325 char z_type_name[32] = "unknown";
326 if (zis->location_id == CS_MESH_LOCATION_BOUNDARY_FACES)
327 strncpy(z_type_name, _("boundary"), 31);
328 else if (zis->location_id == CS_MESH_LOCATION_CELLS)
329 strncpy(z_type_name, _("volume"), 31);
330 z_type_name[31] = '\0';
331
332 int z_id = zis->zone_id;
333 int set_id = zis->set_id;
334
335 cs_lagr_extra_module_t *extra = cs_get_lagr_extra_module();
336
337 /* Verification of particle classes */
338
339 if (cs_glob_lagr_model->n_stat_classes > 0) {
340 if ( zis->cluster < 0
341 || zis->cluster > cs_glob_lagr_model->n_stat_classes)
342 bft_error(__FILE__, __LINE__, 0,
343 _("Lagrangian module: \n"
344 " number of clusters = %d is either not defined (negative)\n"
345 " or > to the number of statistical classes %d\n"
346 " for zone %d and set %d."),
347 (int)zis->cluster,
348 (int)cs_glob_lagr_model->n_stat_classes,
349 z_id,
350 set_id);
351 }
352
353 if (cs_glob_lagr_model->agglomeration == 1) {
354 if (zis->aggregat_class_id < 1)
355 bft_error(__FILE__, __LINE__, 0,
356 _("Lagrangian module: \n"
357 " id of the class of aggregates = %d is \n"
358 " either not defined (negative)\n"
359 " or smaller than 1 \n"
360 " for zone %d and set %d."),
361 (int)zis->aggregat_class_id,
362 z_id,
363 set_id);
364 if (zis->aggregat_fractal_dim > 3.)
365 bft_error(__FILE__, __LINE__, 0,
366 _("Lagrangian module: \n"
367 " value of fractal dimension = %d for aggregates \n"
368 " is either not defined (negative) \n"
369 " or greater than 3 \n"
370 " for zone %d and set %d."),
371 (int)zis->aggregat_fractal_dim,
372 z_id,
373 set_id);
374
375 }
376
377 /* temperature */
378 if ( cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_HEAT
379 && ( cs_glob_lagr_specific_physics->itpvar == 1
380 || cs_glob_lagr_specific_physics->idpvar == 1
381 || cs_glob_lagr_specific_physics->impvar == 1)) {
382 if (zis->temperature_profile < 0 || zis->temperature_profile > 1)
383 bft_error(__FILE__, __LINE__, 0, _profile_err_fmt_i,
384 z_type_name, z_id, set_id,
385 _("temperature"), (int)zis->temperature_profile);
386 }
387
388 /* velocity */
389 if ( zis->location_id != CS_MESH_LOCATION_BOUNDARY_FACES
390 && zis->velocity_profile == CS_LAGR_IN_IMPOSED_NORM)
391 bft_error(__FILE__, __LINE__, 0,
392 _("Lagrangian %s zone %d, set %d:\n"
393 " velocity profile type CS_LAGR_IN_IMPOSED_NORM may not be used\n"
394 " for volume zones, as it requires surface normals."),
395 z_type_name, z_id, set_id);
396 else if (zis->velocity_profile < CS_LAGR_IN_IMPOSED_FLUID_VALUE
397 || zis->velocity_profile > CS_LAGR_IN_IMPOSED_COMPONENTS)
398 bft_error(__FILE__, __LINE__, 0, _profile_err_fmt_i,
399 z_type_name, z_id, set_id,
400 _("velocity"), (int)zis->velocity_profile);
401
402 /* statistical weight */
403 if (zis->stat_weight <= 0.0 && zis->flow_rate <= 0.0)
404 bft_error(__FILE__, __LINE__, 0, _profile_err_fmt_d,
405 z_type_name, z_id, set_id,
406 _("statistical weight"), (double)zis->stat_weight);
407
408 /* mass flow rate */
409 if (zis->flow_rate > 0.0 && zis->n_inject == 0)
410 bft_error(__FILE__, __LINE__, 0,
411 _("Lagrangian %s zone %d, set %d:\n"
412 " flow rate is positive (%g)\n"
413 " while number injected particles is 0."),
414 z_type_name, z_id, set_id,
415 (double)zis->flow_rate);
416
417 /* particle properties: diameter, variance, and rho */
418 if (cs_glob_lagr_model->physical_model != CS_LAGR_PHYS_COAL) {
419 if ( zis->density < 0.0
420 || zis->diameter < 0.0
421 || zis->diameter_variance < 0.0)
422 bft_error
423 (__FILE__, __LINE__, 0,
424 _("Lagrangian %s zone %d, set %d:\n"
425 " error on particle properties definition:\n"
426 " rho = %g, diameter = %g,\n"
427 " diameter standard deviation = %g\n"
428 "This may lead to injection of particles with negative diameters."),
429 z_type_name, z_id, set_id,
430 (double)zis->density,
431 (double)zis->diameter,
432 (double)zis->diameter_variance);
433 }
434
435 if (zis->diameter < 3.0 * zis->diameter_variance)
436 bft_error(__FILE__, __LINE__, 0,
437 _("Lagrangian %s zone %d, set %d:\n"
438 " diameter (%g) is smaller than 3 times\n"
439 " its standard deviation (%g)."),
440 z_type_name, z_id, set_id,
441 (double)zis->diameter,
442 (double)zis->diameter_variance);
443
444 /* Ellipsoidal particle properties: radii */
445 const int shape = cs_glob_lagr_model->shape;
446 if (shape == CS_LAGR_SHAPE_SPHEROID_JEFFERY_MODEL) {
447 if ( zis->radii[0] < 0.0
448 || zis->radii[1] < 0.0
449 || zis->radii[2] < 0.0)
450 bft_error
451 (__FILE__, __LINE__, 0,
452 _("Lagrangian %s zone %d, set %d:\n"
453 " error on particle properties definition:\n"
454 " Ellispoid radii = %g, %g, %g\n"
455 "This may lead to injection of particles with negative radii."),
456 z_type_name, z_id, set_id,
457 (double)zis->radii[0],
458 (double)zis->radii[1],
459 (double)zis->radii[2]);
460 }
461
462 /* temperature and Cp */
463 if ( cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_HEAT
464 && cs_glob_lagr_specific_physics->itpvar == 1) {
465 cs_real_t tkelvn = -cs_physical_constants_celsius_to_kelvin;
466 if (zis->cp < 0.0)
467 bft_error(__FILE__, __LINE__, 0,
468 _("Lagrangian %s zone %d, set %d:\n"
469 " specific heat capacity (%g) is negative."),
470 z_type_name, z_id, set_id, (double)zis->cp);
471 if (zis->temperature_profile > 0 && zis->temperature < tkelvn)
472 bft_error(__FILE__, __LINE__, 0,
473 _("Lagrangian %s zone %d, set %d:\n"
474 " temperature (%g) is lower than %g."),
475 z_type_name, z_id, set_id,
476 (double)zis->temperature,
477 (double)tkelvn);
478 }
479
480 /* emissivity */
481 if ( cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_HEAT
482 && cs_glob_lagr_specific_physics->itpvar == 1
483 && extra->radiative_model > 0) {
484
485 if (zis->emissivity < 0.0 || zis->emissivity > 1.0)
486 bft_error(__FILE__, __LINE__, 0,
487 _("Lagrangian %s zone %d, set %d:\n"
488 " particle emissivity (%g) is not properly set."),
489 z_type_name, z_id, set_id,
490 (double)zis->emissivity);
491
492 }
493
494 /* Coal */
495 if (cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_COAL) {
496
497 cs_real_t tkelvi = cs_physical_constants_celsius_to_kelvin;
498
499 if (zis->coal_number < 1 && zis->coal_number > extra->ncharb)
500 bft_error
501 (__FILE__, __LINE__, 0,
502 _("Lagrangian %s zone %d, set %d:\n"
503 " the coal number %d for the injected particle is either negative\n"
504 " or greater than the maximum number of coals defined (%d)."),
505 z_type_name, z_id, set_id,
506 (int)zis->coal_number, (int)extra->ncharb);
507
508 int coal_id = zis->coal_number - 1;
509
510 /* properties of coal particles */
511 if (zis->temperature < tkelvi)
512 bft_error(__FILE__, __LINE__, 0,
513 _("Lagrangian %s zone %d, set %d:\n"
514 " temperature is not properly set: %g."),
515 z_type_name, z_id, set_id,
516 (double)zis->temperature);
517
518 /* Properties of coal particles */
519
520 /* Composition of coal defined in XML file (DP_FCP) */
521
522 cs_real_t *xashch = cs_glob_lagr_coal_comb->xashch;
523 cs_real_t *cp2ch = cs_glob_lagr_coal_comb->cp2ch;
524 cs_real_t *xwatch = cs_glob_lagr_coal_comb->xwatch;
525 cs_real_t *rho0ch = cs_glob_lagr_coal_comb->rho0ch;
526
527 if ( rho0ch[coal_id] < 0.0
528 || cp2ch[coal_id] < 0.0
529 || xwatch[coal_id] < 0.0
530 || xwatch[coal_id] > 1.0
531 || xashch[coal_id] < 0.0
532 || xashch[coal_id] > 1.0)
533 bft_error(__FILE__, __LINE__, 0,
534 _("Lagrangian %s zone %d, set %d:\n"
535 " wrong conditions for coal number %d.\n"
536 " coal density = %g\n"
537 " Cp CP2CH = %g\n"
538 " water mass fraction = %g\n"
539 " ashes mass fraction = %g."),
540 z_type_name, z_id, set_id,
541 (int)coal_id,
542 (double)rho0ch[coal_id],
543 (double)cp2ch[coal_id],
544 (double)xwatch[coal_id],
545 (double)xashch[coal_id]);
546
547 if (xwatch[coal_id] + xashch[coal_id] > 1.0)
548 bft_error(__FILE__, __LINE__, 0,
549 _("Lagrangian %s zone %d, set %d:\n"
550 " wrong conditions for coal number %d.\n"
551 " water mass fraction = %g\n"
552 " ashes mass fraction = %g\n"
553 " mass fraction is larger than 1: %g."),
554 z_type_name, z_id, set_id,
555 (int)zis->coal_number,
556 (double)xwatch[coal_id],
557 (double)xashch[coal_id],
558 (double)(xwatch[coal_id] + xashch[coal_id]));
559
560 }
561 }
562
563 /*----------------------------------------------------------------------------*/
564 /*!
565 * \brief Build particle injection face ids array for a given boundary
566 * zone and set.
567 *
568 * The caller is responsible for freeing the returned array.
569 *
570 * \param[in] n_faces number of elements in zone
571 * \param[in] face_ids matching face ids
572 * \param[in] face_particle_idx starting id of new particles for a given
573 * face (size: n_faces+1)
574 *
575 * \return array of ids of faces for injected particles
576 * (size: face_particle_idx[n_faces])
577 */
578 /*----------------------------------------------------------------------------*/
579
580 static cs_lnum_t *
_get_particle_face_ids(cs_lnum_t n_faces,const cs_lnum_t face_ids[],const cs_lnum_t face_particle_idx[])581 _get_particle_face_ids(cs_lnum_t n_faces,
582 const cs_lnum_t face_ids[],
583 const cs_lnum_t face_particle_idx[])
584 {
585 cs_lnum_t *particle_face_id = NULL;
586
587 cs_lnum_t n_p_new = face_particle_idx[n_faces];
588
589 BFT_MALLOC(particle_face_id, n_p_new, cs_lnum_t);
590
591 /* Loop on zone elements where particles are injected */
592
593 n_p_new = 0;
594
595 for (cs_lnum_t i = 0; i < n_faces; i++) {
596
597 /* Loop on particles added for this face */
598
599 for (cs_lnum_t j = face_particle_idx[i]; j < face_particle_idx[i+1]; j++)
600 particle_face_id[j] = face_ids[i];
601
602 }
603
604 return(particle_face_id);
605 }
606
607 /*----------------------------------------------------------------------------*/
608 /*!
609 * \brief Initialize particle values
610 *
611 * \param[in,out] p_set particle set
612 * \param[in] zis injection data for this zone and set
613 * \param[in] time_id time step indicator for fields
614 * 0: use fields at current time step
615 * 1: use fields at previous time step
616 * \param[in] n_elts number of elements in zone
617 * \param[in] face_ids matching face ids if zone is a boundary
618 * \param[in] elt_particle_idx starting id of new particles for a given
619 * element (size: n_elts+1)
620 */
621 /*----------------------------------------------------------------------------*/
622
623 static void
_init_particles(cs_lagr_particle_set_t * p_set,const cs_lagr_injection_set_t * zis,int time_id,cs_lnum_t n_elts,const cs_lnum_t * face_ids,const cs_lnum_t elt_particle_idx[])624 _init_particles(cs_lagr_particle_set_t *p_set,
625 const cs_lagr_injection_set_t *zis,
626 int time_id,
627 cs_lnum_t n_elts,
628 const cs_lnum_t *face_ids,
629 const cs_lnum_t elt_particle_idx[])
630 {
631 const cs_lagr_attribute_map_t *p_am = p_set->p_am;
632
633 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
634
635 cs_lagr_extra_module_t *extra = cs_get_lagr_extra_module();
636
637 /* Non-lagrangian fields */
638
639 const cs_real_t *xashch = cs_glob_lagr_coal_comb->xashch;
640 const cs_real_t *cp2ch = cs_glob_lagr_coal_comb->cp2ch;
641 const cs_real_t *xwatch = cs_glob_lagr_coal_comb->xwatch;
642 const cs_real_t *rho0ch = cs_glob_lagr_coal_comb->rho0ch;
643
644 const cs_real_t *vela = extra->vel->vals[time_id];
645 const cs_real_t *cval_h = NULL, *cval_t = NULL;
646 cs_real_t *_cval_t = NULL;
647
648 cs_real_t tscl_shift = 0;
649
650 /* Initialize pointers (used to simplify future tests) */
651
652 if ( ( cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_HEAT
653 && cs_glob_lagr_specific_physics->itpvar == 1)
654 || cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_COAL) {
655
656 const cs_field_t *f = cs_field_by_name_try("temperature");
657 if (f != NULL)
658 cval_t = f->val;
659 else if (cs_glob_thermal_model->itherm == CS_THERMAL_MODEL_ENTHALPY)
660 cval_h = cs_field_by_name("enthalpy")->val;
661
662 if (cs_glob_thermal_model->itpscl == CS_TEMPERATURE_SCALE_KELVIN)
663 tscl_shift = - cs_physical_constants_celsius_to_kelvin;
664 }
665
666 const cs_real_t pis6 = cs_math_pi / 6.0;
667 const int shape = cs_glob_lagr_model->shape;
668
669 /* Prepare enthalpy to temperature conversion if needed */
670
671 if ( cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_HEAT
672 && cs_glob_lagr_specific_physics->itpvar == 1
673 && cval_t == NULL
674 && cval_h != NULL) {
675
676 BFT_MALLOC(_cval_t, cs_glob_mesh->n_cells_with_ghosts, cs_real_t);
677 cs_ht_convert_h_to_t_cells(cval_h, _cval_t);
678 cval_t = _cval_t;
679
680 }
681
682 /* Loop on zone elements where particles are injected */
683
684 for (cs_lnum_t li = 0; li < n_elts; li++) {
685
686 cs_lnum_t n_e_p = elt_particle_idx[li+1] - elt_particle_idx[li];
687
688 if (n_e_p < 1)
689 continue;
690
691 cs_lnum_t p_s_id = p_set->n_particles + elt_particle_idx[li];
692 cs_lnum_t p_e_id = p_s_id + n_e_p;
693
694 const cs_lnum_t face_id = (face_ids != NULL) ? face_ids[li] : -1;
695
696 /* Loop on particles added for this face */
697
698 for (cs_lnum_t p_id = p_s_id; p_id < p_e_id; p_id++) {
699
700 unsigned char *particle = p_set->p_buffer + p_am->extents * p_id;
701
702 cs_lnum_t cell_id = cs_lagr_particles_get_lnum(p_set, p_id,
703 CS_LAGR_CELL_ID);
704
705 /* Random value associated with each particle */
706
707 cs_real_t part_random = -1;
708 cs_random_uniform(1, &part_random);
709 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_RANDOM_VALUE,
710 part_random);
711
712 /* Particle velocity components */
713
714 cs_real_t *part_vel = cs_lagr_particle_attr(particle, p_am,
715 CS_LAGR_VELOCITY);
716
717 /* prescribed components */
718 if (zis->velocity_profile == CS_LAGR_IN_IMPOSED_COMPONENTS) {
719 for (cs_lnum_t i = 0; i < 3; i++)
720 part_vel[i] = zis->velocity[i];
721 }
722
723 /* prescribed norm */
724 else if (zis->velocity_profile == CS_LAGR_IN_IMPOSED_NORM) {
725 assert(face_id >= 0);
726 for (cs_lnum_t i = 0; i < 3; i++)
727 part_vel[i] = - fvq->b_face_normal[face_id * 3 + i]
728 / fvq->b_face_surf[face_id]
729 * zis->velocity_magnitude;
730 }
731
732 /* velocity as seen from fluid */
733 else if (zis->velocity_profile == CS_LAGR_IN_IMPOSED_FLUID_VALUE) {
734 for (cs_lnum_t i = 0; i < 3; i++)
735 part_vel[i] = vela[cell_id * 3 + i];
736 }
737
738 /* fluid velocity seen */
739 cs_real_t *part_seen_vel = cs_lagr_particle_attr(particle, p_am,
740 CS_LAGR_VELOCITY_SEEN);
741 for (cs_lnum_t i = 0; i < 3; i++)
742 part_seen_vel[i] = vela[cell_id * 3 + i];
743
744 /* Residence time (may be negative to ensure continuous injection) */
745 if (zis->injection_frequency == 1) {
746 cs_real_t res_time = - part_random *cs_glob_lagr_time_step->dtp;
747 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_RESIDENCE_TIME,
748 res_time);
749 }
750 else
751 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_RESIDENCE_TIME,
752 0.0);
753
754 /* Diameter (always set base) */
755
756 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_DIAMETER,
757 zis->diameter);
758
759 /* Shape for spheroids without inertia */
760 if ( shape == CS_LAGR_SHAPE_SPHEROID_STOC_MODEL
761 || shape == CS_LAGR_SHAPE_SPHEROID_JEFFERY_MODEL) {
762
763 /* Spherical radii a b c */
764 cs_real_t *radii = cs_lagr_particle_attr(particle, p_am,
765 CS_LAGR_RADII);
766
767 for (cs_lnum_t i = 0; i < 3; i++) {
768 radii[i] = zis->radii[i];
769 }
770
771 /* Shape parameters */
772 cs_real_t *shape_param = cs_lagr_particle_attr(particle, p_am,
773 CS_LAGR_SHAPE_PARAM);
774
775 /* Compute shape parameters from radii */
776 /* FIXME valid for all spheroids only (a = b, c > a,b ) */
777 cs_real_t lamb = radii[2] / radii[1]; /* FIXME do not divide by 0... */
778 cs_real_t lamb_m1 = (radii[2] - radii[1]) / radii[1];
779 cs_real_t _a2 = radii[0] * radii[0];
780 /* TODO MF shape_param check development in series */
781 cs_real_t aux1 = lamb * lamb;
782 cs_real_t aux2 = aux1 -1;
783 if (lamb_m1 > 1e-10) {
784 cs_real_t aux3 = sqrt(aux2 - 1);
785 cs_real_t kappa = -log(lamb + aux3);
786 shape_param[0] = aux1/aux2 + lamb*kappa/(aux2*aux3);
787 shape_param[1] = shape_param[0];
788 shape_param[2] = -2./aux2 - 2.*lamb*kappa/(aux2*aux3);
789 shape_param[3] = -2. * _a2 *lamb*kappa/aux3;
790 }
791 else if (lamb_m1 < -1e-10) {
792 cs_real_t aux3 = sqrt(1. - aux2);
793 cs_real_t kappa = acos(lamb);
794 shape_param[0] = aux1/aux2+lamb*kappa/(-aux2*aux3);
795 shape_param[1] = shape_param[0];
796 shape_param[2] = -2./aux2 - 2.*lamb*kappa/(-aux2*aux3);
797 shape_param[3] = 2. * _a2 * lamb*kappa/aux3;
798 }
799 else {
800 shape_param[0] = 2.0/3.0;
801 shape_param[1] = 2.0/3.0;
802 shape_param[2] = 2.0/3.0;
803 shape_param[3] = 2. * _a2;
804 }
805
806 if (shape == CS_LAGR_SHAPE_SPHEROID_STOC_MODEL) {
807 /* Orientation */
808 cs_real_t *orientation = cs_lagr_particle_attr(particle, p_am,
809 CS_LAGR_ORIENTATION);
810 cs_real_t *quaternion = cs_lagr_particle_attr(particle, p_am,
811 CS_LAGR_QUATERNION);
812 for (cs_lnum_t i = 0; i < 3; i++) {
813 orientation[i] = zis->orientation[i];
814 }
815
816 /* Compute orientation from uniform orientation on a unit-sphere */
817 cs_real_t theta0;
818 cs_real_t phi0;
819 cs_random_uniform(1, &theta0);
820 cs_random_uniform(1, &phi0);
821 theta0 = acos(2.0*theta0-1.0);
822 phi0 = phi0*2.0*cs_math_pi;
823 orientation[0] = sin(theta0)*cos(phi0);
824 orientation[1] = sin(theta0)*sin(phi0);
825 orientation[2] = cos(theta0);
826 /* Initialize quaternions */
827 quaternion[0] = 1.;
828 quaternion[1] = 0.;
829 quaternion[2] = 0.;
830 quaternion[3] = 0.;
831 /* TODO initialize other things */
832 }
833
834 if (shape == CS_LAGR_SHAPE_SPHEROID_JEFFERY_MODEL) {
835
836 /* Euler parameters */
837 cs_real_t *euler = cs_lagr_particle_attr(particle, p_am,
838 CS_LAGR_EULER);
839
840 for (cs_lnum_t i = 0; i < 4; i++)
841 euler[i] = zis->euler[i];
842
843 /* Compute Euler angles
844 (random orientation with a uniform distribution in [-1;1]) */
845 cs_real_t trans_m[3][3];
846 /* Generate the first two vectors */
847 for (cs_lnum_t id = 0; id < 3; id++) {
848 cs_random_uniform(1, &trans_m[id][0]); /* (?,0) */
849 cs_random_uniform(1, &trans_m[id][1]); /* (?,1) */
850 cs_random_uniform(1, &trans_m[id][2]); /* (?,2) */
851 cs_real_3_t loc_vector = {-1.+2*trans_m[id][0],
852 -1.+2*trans_m[id][1],
853 -1.+2*trans_m[id][2]};
854 cs_real_t norm_trans_m = cs_math_3_norm( loc_vector );
855 while ( norm_trans_m > 1 )
856 {
857 cs_random_uniform(1, &trans_m[id][0]); /* (?,0) */
858 cs_random_uniform(1, &trans_m[id][1]); /* (?,1) */
859 cs_random_uniform(1, &trans_m[id][2]); /* (?,2) */
860 loc_vector[0] = -1.+2*trans_m[id][0];
861 loc_vector[1] = -1.+2*trans_m[id][1];
862 loc_vector[2] = -1.+2*trans_m[id][2];
863 norm_trans_m = cs_math_3_norm( loc_vector );
864 }
865 for (cs_lnum_t id1 = 0; id1 < 3; id1++)
866 trans_m[id][id1] = (-1.+2*trans_m[id][id1]) / norm_trans_m;
867 }
868 /* Correct 2nd vector (for perpendicularity to the 1st) */
869 cs_real_3_t loc_vector0 = {trans_m[0][0],
870 trans_m[0][1],
871 trans_m[0][2]};
872 cs_real_3_t loc_vector1 = {trans_m[1][0],
873 trans_m[1][1],
874 trans_m[1][2]};
875 cs_real_t scal_prod = cs_math_3_dot_product(loc_vector0, loc_vector1);
876 for (cs_lnum_t id = 0; id < 3; id++)
877 trans_m[1][id] -= scal_prod * trans_m[0][id];
878 /* Re-normalize */
879 loc_vector1[0] = trans_m[1][0];
880 loc_vector1[1] = trans_m[1][1];
881 loc_vector1[2] = trans_m[1][2];
882 cs_real_t norm_trans_m = cs_math_3_norm( loc_vector1 );
883 for (cs_lnum_t id = 0; id < 3; id++)
884 trans_m[1][id] /= norm_trans_m;
885
886 /* Compute last vector (cross product of the two others) */
887 loc_vector1[0] = trans_m[1][0];
888 loc_vector1[1] = trans_m[1][1];
889 loc_vector1[2] = trans_m[1][2];
890 cs_real_3_t loc_vector2 = {trans_m[2][0],
891 trans_m[2][1],
892 trans_m[2][2]};
893 cs_math_3_cross_product( loc_vector0, loc_vector1, loc_vector2);
894 for (cs_lnum_t id = 0; id < 3; id++)
895 trans_m[2][id] = loc_vector2[id];
896
897 /* Write Euler angles */
898 cs_real_t random;
899 cs_random_uniform(1, &random);
900 if (random >= 0.5)
901 euler[0] = pow(0.25*(trans_m[0][0]+trans_m[1][1]+trans_m[2][2]+1.),
902 0.5);
903 else
904 euler[0] = -pow(0.25*(trans_m[0][0]+trans_m[1][1]+trans_m[2][2]+1.),
905 0.5);
906 euler[1] = 0.25 * (trans_m[2][1] - trans_m[1][2]) / euler[0];
907 euler[2] = 0.25 * (trans_m[0][2] - trans_m[2][0]) / euler[0];
908 euler[3] = 0.25 * (trans_m[1][0] - trans_m[0][1]) / euler[0];
909
910 /* Compute initial angular velocity */
911 /* Get velocity gradient */
912 cs_lagr_gradients(0, extra->grad_pr, extra->grad_vel);
913
914 /* Local reference frame */
915 cs_real_t grad_vf_r[3][3];
916 cs_math_33_transform_a_to_r(extra->grad_vel[cell_id],
917 trans_m,
918 grad_vf_r);
919
920 cs_real_t *ang_vel = cs_lagr_particle_attr(particle, p_am,
921 CS_LAGR_ANGULAR_VEL);
922
923 ang_vel[0] = 0.5*(grad_vf_r[2][1] - grad_vf_r[1][2]);
924 ang_vel[1] = 0.5*(grad_vf_r[0][2] - grad_vf_r[2][0]);
925 ang_vel[2] = 0.5*(grad_vf_r[0][1] - grad_vf_r[1][0]);
926 }
927
928 }
929
930 if (zis->diameter_variance > 0.0) {
931
932 /* Randomize diameter, ensuring we obtain a
933 positive diameter in the 99,7% range */
934
935 cs_real_t d3 = 3.0 * zis->diameter_variance;
936
937 int i_r = 0; /* avoid infinite loop in case of very improbable
938 random series... */
939
940 for (i_r = 0; i_r < 20; i_r++) {
941 double random;
942 cs_random_normal(1, &random);
943
944 cs_real_t diam = zis->diameter
945 + random * zis->diameter_variance;
946
947 if (diam > 0 && ( diam >= zis->diameter - d3
948 && diam <= zis->diameter + d3)) {
949 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_DIAMETER, diam);
950 break;
951 }
952 }
953
954 }
955
956 /* Other parameters */
957 cs_real_t diam = cs_lagr_particle_get_real(particle, p_am,
958 CS_LAGR_DIAMETER);
959 cs_real_t mporos = cs_glob_lagr_clogging_model->mporos;
960 if (cs_glob_lagr_model->clogging == 1) {
961 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_DIAMETER,
962 diam/(1.-mporos));
963 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_HEIGHT, diam);
964 }
965
966 /* Other variables (mass, ...) depending on physical model */
967 cs_real_t d3 = pow(diam, 3.0);
968
969 if (cs_glob_lagr_model->n_stat_classes > 0)
970 cs_lagr_particle_set_lnum(particle, p_am, CS_LAGR_STAT_CLASS,
971 zis->cluster);
972
973 if ( cs_glob_lagr_model->agglomeration == 1
974 || cs_glob_lagr_model->fragmentation == 1) {
975 cs_lagr_particle_set_lnum(particle, p_am, CS_LAGR_AGGLO_CLASS_ID,
976 zis->aggregat_class_id);
977 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_AGGLO_FRACTAL_DIM,
978 zis->aggregat_fractal_dim);
979 }
980
981 /* used for 2nd order only */
982 if (p_am->displ[0][CS_LAGR_TAUP_AUX] > 0)
983 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_TAUP_AUX, 0.0);
984
985 if ( cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_OFF
986 || cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_HEAT) {
987
988 if (cs_glob_lagr_model->clogging == 0)
989 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_MASS,
990 zis->density * pis6 * d3);
991 else
992 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_MASS,
993 zis->density * pis6 * d3
994 * pow(1.0-mporos, 3));
995
996 if ( cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_HEAT
997 && cs_glob_lagr_specific_physics->itpvar == 1) {
998
999 if (zis->temperature_profile < 1) {
1000 if (cval_t != NULL)
1001 cs_lagr_particle_set_real(particle, p_am,
1002 CS_LAGR_FLUID_TEMPERATURE,
1003 cval_t[cell_id] + tscl_shift);
1004 }
1005
1006 /* constant temperature set, may be modified later by user function */
1007 else if (zis->temperature_profile == 1)
1008 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_TEMPERATURE,
1009 zis->temperature);
1010
1011 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_CP,
1012 zis->cp);
1013 if (extra->radiative_model > 0)
1014 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_EMISSIVITY,
1015 zis->emissivity);
1016
1017 }
1018
1019 }
1020
1021 else if (cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_COAL) {
1022
1023 int coal_id = zis->coal_number - 1;
1024
1025 cs_lagr_particle_set_lnum(particle, p_am, CS_LAGR_COAL_ID, coal_id);
1026 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_FLUID_TEMPERATURE,
1027 cval_t[cell_id] + tscl_shift);
1028
1029 cs_real_t *particle_temp
1030 = cs_lagr_particle_attr(particle, p_am, CS_LAGR_TEMPERATURE);
1031 for (int ilayer = 0;
1032 ilayer < cs_glob_lagr_model->n_temperature_layers;
1033 ilayer++)
1034 particle_temp[ilayer] = zis->temperature;
1035
1036 /* composition from DP_FCP */
1037
1038 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_CP, cp2ch[coal_id]);
1039
1040 cs_real_t mass = rho0ch[coal_id] * pis6 * d3;
1041
1042 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_MASS, mass);
1043 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_WATER_MASS,
1044 xwatch[coal_id] * mass);
1045
1046 cs_real_t *particle_coal_mass
1047 = cs_lagr_particle_attr(particle, p_am, CS_LAGR_COAL_MASS);
1048 cs_real_t *particle_coke_mass
1049 = cs_lagr_particle_attr(particle, p_am, CS_LAGR_COKE_MASS);
1050 for (int ilayer = 0;
1051 ilayer < cs_glob_lagr_model->n_temperature_layers;
1052 ilayer++) {
1053
1054 particle_coal_mass[ilayer]
1055 = (1.0 - xwatch[coal_id]
1056 - xashch[coal_id])
1057 * cs_lagr_particle_get_real(particle, p_am, CS_LAGR_MASS)
1058 / cs_glob_lagr_model->n_temperature_layers;
1059 particle_coke_mass[ilayer] = 0.0;
1060
1061 }
1062
1063 cs_lagr_particle_set_real
1064 (particle, p_am,
1065 CS_LAGR_SHRINKING_DIAMETER,
1066 cs_lagr_particle_get_real(particle, p_am, CS_LAGR_DIAMETER));
1067 cs_lagr_particle_set_real
1068 (particle, p_am,
1069 CS_LAGR_INITIAL_DIAMETER,
1070 cs_lagr_particle_get_real(particle, p_am, CS_LAGR_DIAMETER));
1071
1072 cs_real_t *particle_coal_density
1073 = cs_lagr_particle_attr(particle, p_am, CS_LAGR_COAL_DENSITY);
1074 for (int ilayer = 0;
1075 ilayer < cs_glob_lagr_model->n_temperature_layers;
1076 ilayer++)
1077 particle_coal_density[ilayer] = rho0ch[coal_id];
1078
1079 }
1080
1081 /* statistical weight */
1082 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_STAT_WEIGHT,
1083 zis->stat_weight);
1084
1085 /* Fouling index */
1086 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_FOULING_INDEX,
1087 zis->fouling_index);
1088
1089 /* Initialization of deposition model */
1090
1091 if (cs_glob_lagr_model->deposition == 1) {
1092
1093 cs_real_t random;
1094 cs_random_uniform(1, &random);
1095 cs_lagr_particle_set_real(particle, p_am,
1096 CS_LAGR_INTERF, 5.0 + 15.0 * random);
1097 cs_lagr_particle_set_real(particle, p_am,
1098 CS_LAGR_YPLUS, 1000.0);
1099 cs_lagr_particle_set_lnum(particle, p_am,
1100 CS_LAGR_MARKO_VALUE, -1);
1101 cs_lagr_particle_set_lnum(particle, p_am,
1102 CS_LAGR_NEIGHBOR_FACE_ID, -1);
1103
1104 }
1105
1106 /* Initialization of clogging model */
1107
1108 if (cs_glob_lagr_model->clogging == 1) {
1109
1110 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_DEPO_TIME, 0.0);
1111 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_CONSOL_HEIGHT, 0.0);
1112 cs_lagr_particle_set_real(particle, p_am, CS_LAGR_CLUSTER_NB_PART, 1.0);
1113
1114 }
1115
1116 /* Initialize the additional user variables */
1117
1118 if (cs_glob_lagr_model->n_user_variables > 0) {
1119 cs_real_t *user_attr
1120 = cs_lagr_particle_attr(particle, p_am, CS_LAGR_USER);
1121 for (int i = 0;
1122 i < cs_glob_lagr_model->n_user_variables;
1123 i++)
1124 user_attr[i] = 0.0;
1125 }
1126 }
1127
1128 }
1129
1130 /* Update weights to have the correct flow rate
1131 -------------------------------------------- */
1132
1133 if (zis->flow_rate > 0.0 && zis->n_inject > 0) {
1134
1135 cs_real_t dmass = 0.0;
1136
1137 cs_lnum_t p_s_id = p_set->n_particles;
1138 cs_lnum_t p_e_id = p_s_id + elt_particle_idx[n_elts];
1139
1140 for (cs_lnum_t p_id = p_s_id; p_id < p_e_id; p_id++)
1141 dmass += cs_lagr_particles_get_real(p_set, p_id, CS_LAGR_MASS);
1142
1143 cs_parall_sum(1, CS_REAL_TYPE, &dmass);
1144
1145 /* Compute weights */
1146
1147 if (dmass > 0.0) {
1148 cs_real_t s_weight = zis->flow_rate * cs_glob_lagr_time_step->dtp
1149 / dmass;
1150 for (cs_lnum_t p_id = p_s_id; p_id < p_e_id; p_id++)
1151 cs_lagr_particles_set_real(p_set, p_id, CS_LAGR_STAT_WEIGHT, s_weight);
1152 }
1153
1154 else {
1155
1156 char z_type_name[32] = "unknown";
1157 if (zis->location_id == CS_MESH_LOCATION_BOUNDARY_FACES)
1158 strncpy(z_type_name, _("boundary"), 31);
1159 else if (zis->location_id == CS_MESH_LOCATION_CELLS)
1160 strncpy(z_type_name, _("volume"), 31);
1161 z_type_name[31] = '\0';
1162
1163 bft_error(__FILE__, __LINE__, 0,
1164 _("Lagrangian %s zone %d, set %d:\n"
1165 " imposed flow rate is %g\n"
1166 " while mass of injected particles is 0."),
1167 z_type_name, zis->zone_id, zis->set_id,
1168 (double)zis->flow_rate);
1169
1170 }
1171
1172 }
1173
1174 BFT_FREE(_cval_t);
1175 }
1176
1177 /*----------------------------------------------------------------------------*/
1178 /*!
1179 * \brief Check injected particle values
1180 *
1181 * \param[in,out] p_set particle set
1182 * \param[in] zis injection data for this zone and set
1183 * \param[in] n_elts number of elements in zone
1184 * \param[in] elt_particle_idx starting id of new particles for a given
1185 * element (size: n_elts+1)
1186 */
1187 /*----------------------------------------------------------------------------*/
1188
1189 static void
_check_particles(cs_lagr_particle_set_t * p_set,const cs_lagr_injection_set_t * zis,cs_lnum_t n_elts,const cs_lnum_t elt_particle_idx[])1190 _check_particles(cs_lagr_particle_set_t *p_set,
1191 const cs_lagr_injection_set_t *zis,
1192 cs_lnum_t n_elts,
1193 const cs_lnum_t elt_particle_idx[])
1194 {
1195 const cs_lnum_t s_id = p_set->n_particles;
1196 const cs_lnum_t e_id = s_id + elt_particle_idx[n_elts];
1197
1198 char z_type_name[32] = "unknown";
1199 if (zis->location_id == CS_MESH_LOCATION_BOUNDARY_FACES)
1200 strncpy(z_type_name, _("boundary"), 31);
1201 else if (zis->location_id == CS_MESH_LOCATION_CELLS)
1202 strncpy(z_type_name, _("volume"), 31);
1203 z_type_name[31] = '\0';
1204
1205 int attrs[] = {CS_LAGR_DIAMETER, CS_LAGR_MASS, CS_LAGR_STAT_WEIGHT,
1206 CS_LAGR_CP};
1207
1208 for (cs_lnum_t p_id = s_id; p_id < e_id; p_id++) {
1209
1210 for (int i_attr = 0; i_attr < 4; i_attr++) {
1211
1212 int attr = attrs[i_attr];
1213
1214 if (p_set->p_am->count[1][attr] > 0) {
1215
1216 cs_real_t val = cs_lagr_particles_get_real(p_set, p_id, attr);
1217
1218 if (val <= 0.0)
1219 bft_error(__FILE__, __LINE__, 0,
1220 _("Lagrangian %s zone %d, set %d:\n"
1221 " particle %ld has a negative %s: %g"),
1222 z_type_name, zis->zone_id, zis->set_id,
1223 (long)p_id, cs_lagr_attribute_name[attr], (double)val);
1224
1225 }
1226
1227 }
1228
1229 }
1230
1231 if (cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_COAL) {
1232
1233 int r01_attrs[] = {CS_LAGR_WATER_MASS, CS_LAGR_COAL_MASS, CS_LAGR_COKE_MASS,
1234 CS_LAGR_COAL_DENSITY};
1235 int r00_attrs[] = {CS_LAGR_SHRINKING_DIAMETER, CS_LAGR_INITIAL_DIAMETER};
1236
1237 for (cs_lnum_t p_id = s_id; p_id < e_id; p_id++) {
1238
1239 for (int i_attr = 0; i_attr < 4; i_attr++) {
1240
1241 int attr = r01_attrs[i_attr];
1242 int n_vals = p_set->p_am->count[1][attr];
1243 cs_real_t *vals = cs_lagr_particles_attr(p_set, p_id, attr);
1244
1245 for (int l_id = 0; l_id < n_vals; l_id++) {
1246 if (vals[l_id] < 0.0) {
1247 if (n_vals == 1)
1248 bft_error(__FILE__, __LINE__, 0,
1249 _("Lagrangian %s zone %d, set %d:\n"
1250 " particle %ld has a negative %s: %g"),
1251 z_type_name, zis->zone_id, zis->set_id,
1252 (long)p_id, cs_lagr_attribute_name[attr],
1253 (double)vals[0]);
1254 else
1255 bft_error(__FILE__, __LINE__, 0,
1256 _("Lagrangian %s zone %d, set %d:\n"
1257 " particle %ld has a negative %s\n"
1258 " in layer %d: %g"),
1259 z_type_name, zis->zone_id, zis->set_id,
1260 (long)p_id, cs_lagr_attribute_name[attr], l_id,
1261 (double)vals[l_id]);
1262 }
1263
1264 }
1265
1266 }
1267
1268 for (int i_attr = 0; i_attr < 2; i_attr++) {
1269
1270 int attr = r00_attrs[i_attr];
1271 cs_real_t val = cs_lagr_particles_get_real(p_set, p_id, attr);
1272
1273 if (val < 0) {
1274 bft_error(__FILE__, __LINE__, 0,
1275 _("Lagrangian %s zone %d, set %d:\n"
1276 " particle %ld has a negative %s: %g"),
1277 z_type_name, zis->zone_id, zis->set_id,
1278 (long)p_id, cs_lagr_attribute_name[attr], (double)val);
1279
1280 }
1281
1282 }
1283
1284 }
1285
1286 }
1287 }
1288
1289 /*============================================================================
1290 * Public function definitions
1291 *============================================================================*/
1292
1293 /*----------------------------------------------------------------------------*/
1294 /*!
1295 * \brief Inject particles in the computational domain.
1296 *
1297 * \param[in] time_id time step indicator for fields
1298 * 0: use fields at current time step
1299 * 1: use fields at previous time step
1300 * \param[in] itypfb boundary face types
1301 * \param[in] visc_length viscous layer thickness
1302 * (size: number of mesh boundary faces)
1303 */
1304 /*----------------------------------------------------------------------------*/
1305
1306 void
cs_lagr_injection(int time_id,const int itypfb[],cs_real_t visc_length[])1307 cs_lagr_injection(int time_id,
1308 const int itypfb[],
1309 cs_real_t visc_length[])
1310 {
1311 CS_UNUSED(itypfb);
1312
1313 /* We may be mapped to an auxiliary field with no previous time id */
1314
1315 cs_real_t dnbpnw_preci = 0.;
1316
1317 cs_lagr_extra_module_t *extra = cs_get_lagr_extra_module();
1318 time_id = CS_MIN(time_id, extra->vel->n_time_vals -1);
1319
1320 const cs_mesh_t *mesh = cs_glob_mesh;
1321 const cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
1322
1323 /* Particles management */
1324 cs_lagr_particle_set_t *p_set = cs_glob_lagr_particle_set;
1325
1326 /* Mean fluid velocity field */
1327 cs_real_t *vela = extra->vel->vals[time_id];
1328
1329 cs_lagr_particle_counter_t *pc = cs_lagr_get_particle_counter();
1330 const cs_time_step_t *ts = cs_glob_time_step;
1331
1332 const int n_stats = cs_glob_lagr_model->n_stat_classes + 1;
1333
1334 /* Initialization */
1335
1336 cs_lagr_zone_data_t *zda[2] = {cs_lagr_get_boundary_conditions(),
1337 cs_lagr_get_volume_conditions()};
1338
1339 cs_lagr_get_internal_conditions();
1340
1341 /* Boundary conditions */
1342
1343 {
1344 cs_lagr_zone_data_t *zd = zda[0];
1345
1346 for (int z_id = 0; z_id < zd->n_zones; z_id++) {
1347
1348 if (zd->zone_type[z_id] > CS_LAGR_BC_USER)
1349 bft_error(__FILE__, __LINE__, 0,
1350 _("Lagrangian boundary zone %d nature %d is unknown."),
1351 z_id + 1,
1352 (int)zd->zone_type[z_id]);
1353
1354 if ( zd->zone_type[z_id] == CS_LAGR_FOULING
1355 && cs_glob_lagr_model->physical_model != CS_LAGR_PHYS_COAL)
1356 bft_error
1357 (__FILE__, __LINE__, 0,
1358 _("Lagrangian boundary zone %d nature is of type CS_LAGR_FOULING,\n"
1359 "but cs_glob_lagr_model->physical_model is not equal to CS_LAGR_PHYS_COAL."),
1360 z_id);
1361 if ( zd->zone_type[z_id] == CS_LAGR_FOULING
1362 && cs_glob_lagr_model->fouling != 1)
1363 bft_error
1364 (__FILE__, __LINE__, 0,
1365 _("Lagrangian boundary zone %d nature is of type CS_LAGR_FOULING,\n"
1366 "but fouling is not activated."),
1367 z_id);
1368
1369 }
1370
1371 }
1372
1373 /* Reset some particle counters */
1374
1375 p_set->n_part_new = 0;
1376 p_set->weight_new = 0.0;
1377
1378 for (int i_loc = 0; i_loc < 2; i_loc++) {
1379 cs_lagr_zone_data_t *zd = zda[i_loc];
1380 int fr_size = zd->n_zones * n_stats;
1381 for (int i = 0; i < fr_size; i++)
1382 zd->particle_flow_rate[i] = 0;
1383 }
1384
1385 /* Injection due to precipitation/Dissolution
1386 ------------------------------------------ */
1387
1388 if (cs_glob_lagr_model->precipitation == 1)
1389 cs_lagr_precipitation_injection(vela, &dnbpnw_preci);
1390
1391 /* User-defined injection
1392 ---------------------- */
1393
1394 /* Check various condition types and optional maximum particle limit */
1395
1396 unsigned long long n_g_particles_next = pc->n_g_total;
1397
1398 for (int i_loc = 0; i_loc < 2; i_loc++) {
1399
1400 cs_lagr_zone_data_t *zd = zda[i_loc];
1401
1402 /* compute global number of injected particles */
1403
1404 for (int z_id = 0; z_id < zd->n_zones; z_id++) {
1405 for (int set_id = 0; set_id < zd->n_injection_sets[z_id]; set_id++) {
1406 cs_lagr_injection_set_t *zis
1407 = cs_lagr_get_injection_set(zd, z_id, set_id);
1408 _injection_check(zis);
1409 n_g_particles_next += (unsigned long long) (zis->n_inject);
1410 }
1411 }
1412
1413 }
1414
1415 /* Avoid injection if maximum defined number of particles reached */
1416
1417 if (n_g_particles_next > cs_lagr_get_n_g_particles_max()) {
1418
1419 bft_printf(_("\n Lagrangian module: \n"));
1420 bft_printf
1421 (_(" If particles are injected at time step %d,\n"
1422 " the total number of particles in the domain would increase from\n"
1423 " %llu to %llu, exceeding the maximums set by\n"
1424 " cs_lagr_set_n_g_particles_max. (%llu).\n"
1425 " No particles will be injected for this time step.\n"),
1426 ts->nt_cur,
1427 (unsigned long long)(pc->n_g_total),
1428 (unsigned long long)n_g_particles_next,
1429 (unsigned long long)(cs_lagr_get_n_g_particles_max()));
1430
1431 return;
1432
1433 }
1434
1435 /* Now inject new particles
1436 ------------------------ */
1437
1438 cs_lnum_t n_elts_m = CS_MAX(mesh->n_b_faces, mesh->n_cells);
1439 cs_lnum_t *elt_particle_idx = NULL;
1440 BFT_MALLOC(elt_particle_idx, n_elts_m+1, cs_lnum_t);
1441
1442 /* Loop in injection type (boundary, volume) */
1443
1444 for (int i_loc = 0; i_loc < 2; i_loc++) {
1445
1446 cs_lagr_zone_data_t *zd = zda[i_loc];
1447
1448 int n_zones = 0;
1449
1450 const cs_real_t *elt_weight = NULL;
1451
1452 if (i_loc == 0) { /* boundary */
1453 elt_weight = fvq->b_face_surf;
1454 n_zones = cs_boundary_zone_n_zones();
1455 }
1456 else { /* volume */
1457 elt_weight = fvq->cell_vol;
1458 n_zones = cs_volume_zone_n_zones();
1459 }
1460
1461 /* Loop on injection zones */
1462
1463 for (int z_id = 0; z_id < n_zones; z_id++) {
1464
1465 /* Loop on injected sets */
1466
1467 cs_lnum_t n_z_elts = 0;
1468 const cs_lnum_t *z_elt_ids = NULL;
1469
1470 if (i_loc == 0) {
1471 const cs_zone_t *z = cs_boundary_zone_by_id(z_id);
1472 n_z_elts = z->n_elts;
1473 z_elt_ids = z->elt_ids;
1474 }
1475 else {
1476 const cs_zone_t *z = cs_volume_zone_by_id(z_id);
1477 n_z_elts = z->n_elts;
1478 z_elt_ids = z->elt_ids;
1479 }
1480
1481 for (int set_id = 0;
1482 set_id < zd->n_injection_sets[z_id];
1483 set_id++) {
1484
1485 const cs_lagr_injection_set_t *zis = NULL;
1486
1487 zis = cs_lagr_get_injection_set(zd, z_id, set_id);
1488
1489 int injection_frequency = zis->injection_frequency;
1490
1491 /* Inject only at first time step if injection frequency is zero */
1492
1493 if (injection_frequency <= 0) {
1494 if (ts->nt_cur == ts->nt_prev+1 && pc->n_g_cumulative_total == 0)
1495 injection_frequency = ts->nt_cur;
1496 else
1497 injection_frequency = ts->nt_cur+1;
1498 }
1499
1500 if (ts->nt_cur % injection_frequency != 0)
1501 continue;
1502
1503 cs_real_t *elt_profile = NULL;
1504 if (zis->injection_profile_func != NULL) {
1505 BFT_MALLOC(elt_profile, n_z_elts, cs_real_t);
1506 zis->injection_profile_func(zis->zone_id,
1507 zis->location_id,
1508 zis->injection_profile_input,
1509 n_z_elts,
1510 z_elt_ids,
1511 elt_profile);
1512 }
1513
1514 cs_lnum_t n_inject = _distribute_particles(zis->n_inject,
1515 n_z_elts,
1516 z_elt_ids,
1517 elt_weight,
1518 elt_profile,
1519 elt_particle_idx);
1520
1521 BFT_FREE(elt_profile);
1522
1523 if (cs_lagr_particle_set_resize(p_set->n_particles + n_inject) < 0)
1524 bft_error(__FILE__, __LINE__, 0,
1525 "Lagrangian module internal error: \n"
1526 " resizing of particle set impossible but previous\n"
1527 " size computation did not detect this issue.");
1528
1529 /* Define particle coordinates and place on faces/cells */
1530
1531 if (zis->location_id == CS_MESH_LOCATION_BOUNDARY_FACES)
1532 cs_lagr_new(p_set,
1533 n_z_elts,
1534 z_elt_ids,
1535 elt_particle_idx);
1536 else
1537 cs_lagr_new_v(p_set,
1538 n_z_elts,
1539 z_elt_ids,
1540 elt_particle_idx);
1541
1542 BFT_FREE(elt_profile);
1543
1544 /* Initialize other particle attributes */
1545
1546 _init_particles(p_set,
1547 zis,
1548 time_id,
1549 n_z_elts,
1550 z_elt_ids,
1551 elt_particle_idx);
1552
1553 assert(n_inject == elt_particle_idx[n_z_elts]);
1554
1555 cs_lnum_t particle_range[2] = {p_set->n_particles,
1556 p_set->n_particles + n_inject};
1557
1558 cs_lagr_new_particle_init(particle_range,
1559 time_id,
1560 visc_length);
1561
1562 /* Advanced user modification:
1563
1564 WARNING: the user may change the particle coordinates but is
1565 prevented from changing the previous location (otherwise, if
1566 the particle is not in the same cell anymore, it would be lost).
1567
1568 Moreover, a precaution has to be taken when calling
1569 "current to previous" in the tracking stage.
1570 */
1571
1572 {
1573 cs_lnum_t *particle_face_ids = NULL;
1574
1575 if (zis->location_id == CS_MESH_LOCATION_BOUNDARY_FACES)
1576 particle_face_ids = _get_particle_face_ids(n_z_elts,
1577 z_elt_ids,
1578 elt_particle_idx);
1579
1580 cs_lnum_t *saved_cell_id;
1581 cs_real_3_t *saved_coords;
1582 BFT_MALLOC(saved_cell_id, n_inject, cs_lnum_t);
1583 BFT_MALLOC(saved_coords, n_inject, cs_real_3_t);
1584
1585 for (cs_lnum_t i = 0; i < n_inject; i++) {
1586 cs_lnum_t p_id = particle_range[0] + i;
1587
1588 saved_cell_id[i] = cs_lagr_particles_get_lnum(p_set,
1589 p_id,
1590 CS_LAGR_CELL_ID);
1591 const cs_real_t *p_coords
1592 = cs_lagr_particles_attr_const(p_set,
1593 p_id,
1594 CS_LAGR_COORDS);
1595 for (cs_lnum_t j = 0; j < 3; j++)
1596 saved_coords[i][j] = p_coords[j];
1597 }
1598
1599 cs_user_lagr_in(p_set,
1600 zis,
1601 particle_range,
1602 particle_face_ids,
1603 visc_length);
1604
1605 /* For safety, build values at previous time step, but reset saved values
1606 for previous cell number and particle coordinates */
1607
1608 for (cs_lnum_t i = 0; i < n_inject; i++) {
1609 cs_lnum_t p_id = particle_range[0] + i;
1610
1611 cs_lagr_particles_current_to_previous(p_set, p_id);
1612
1613 cs_lagr_particles_set_lnum_n(p_set,
1614 p_id,
1615 1,
1616 CS_LAGR_CELL_ID,
1617 saved_cell_id[i]);
1618 cs_real_t *p_coords_prev
1619 = cs_lagr_particles_attr_n(p_set,
1620 p_id,
1621 1,
1622 CS_LAGR_COORDS);
1623 for (cs_lnum_t j = 0; j < 3; j++)
1624 p_coords_prev[j] = saved_coords[i][j];
1625
1626 /* Just after injection, compute the next particle position with
1627 * a reduce integration time so as to simulate continuous injection
1628 */
1629 cs_real_t res_time = cs_lagr_particles_get_real(p_set, p_id,
1630 CS_LAGR_RESIDENCE_TIME);
1631
1632 if (res_time < 0) {
1633
1634 cs_real_t *p_coords
1635 = cs_lagr_particles_attr(p_set,
1636 p_id,
1637 CS_LAGR_COORDS);
1638 cs_real_t *p_vel =
1639 cs_lagr_particles_attr(p_set, p_id, CS_LAGR_VELOCITY);
1640 cs_real_t t_fraction = (cs_glob_lagr_time_step->dtp + res_time);
1641
1642 for (cs_lnum_t j = 0; j < 3; j++)
1643 p_coords[j] = p_coords_prev[j] + t_fraction * p_vel[j];
1644
1645 }
1646
1647 }
1648
1649 BFT_FREE(saved_coords);
1650 BFT_FREE(saved_cell_id);
1651
1652 /* Add particle tracking events for boundary injection */
1653
1654 if ( particle_face_ids != NULL
1655 && cs_lagr_stat_is_active(CS_LAGR_STAT_GROUP_TRACKING_EVENT)) {
1656
1657 cs_lagr_event_set_t *events
1658 = cs_lagr_event_set_boundary_interaction();
1659
1660 /* Event set "expected" size: n boundary faces*2 */
1661 cs_lnum_t events_min_size = mesh->n_b_faces * 2;
1662 if (events->n_events_max < events_min_size)
1663 cs_lagr_event_set_resize(events, events_min_size);
1664
1665 for (cs_lnum_t i = 0; i < n_inject; i++) {
1666 cs_lnum_t p_id = particle_range[0] + i;
1667
1668 cs_lnum_t event_id = events->n_events;
1669 events->n_events += 1;
1670
1671 if (event_id >= events->n_events_max) {
1672 /* flush events */
1673 cs_lagr_stat_update_event(events,
1674 CS_LAGR_STAT_GROUP_TRACKING_EVENT);
1675 events->n_events = 0;
1676 event_id = 0;
1677 }
1678
1679 cs_lagr_event_init_from_particle(events, p_set, event_id, p_id);
1680
1681 cs_lnum_t face_id = particle_face_ids[i];
1682 cs_lagr_events_set_lnum(events,
1683 event_id,
1684 CS_LAGR_E_FACE_ID,
1685 face_id);
1686
1687 cs_lnum_t *e_flag = cs_lagr_events_attr(events,
1688 event_id,
1689 CS_LAGR_E_FLAG);
1690
1691 *e_flag = *e_flag | CS_EVENT_INFLOW;
1692
1693 }
1694
1695 }
1696
1697 BFT_FREE(particle_face_ids);
1698
1699 }
1700
1701 /* check some particle attributes consistency */
1702
1703 _check_particles(p_set, zis, n_z_elts, elt_particle_idx);
1704
1705 /* update counters and balances */
1706
1707 cs_real_t z_weight = 0.;
1708
1709 for (cs_lnum_t p_id = particle_range[0];
1710 p_id < particle_range[1];
1711 p_id++) {
1712 cs_real_t s_weight = cs_lagr_particles_get_real(p_set, p_id,
1713 CS_LAGR_STAT_WEIGHT);
1714 cs_real_t flow_rate = ( s_weight
1715 * cs_lagr_particles_get_real(p_set, p_id,
1716 CS_LAGR_MASS));
1717
1718 zd->particle_flow_rate[z_id*n_stats] += flow_rate;
1719
1720 if (n_stats > 1) {
1721 int class_id = cs_lagr_particles_get_lnum(p_set, p_id,
1722 CS_LAGR_STAT_CLASS);
1723 if (class_id > 0 && class_id < n_stats)
1724 zd->particle_flow_rate[z_id*n_stats + class_id] += flow_rate;
1725 }
1726
1727 z_weight += s_weight;
1728 }
1729
1730 p_set->n_particles += n_inject;
1731 p_set->n_part_new += n_inject;
1732 p_set->weight_new += z_weight;
1733
1734 } /* end of loop on sets */
1735
1736 } /* end of loop on zones */
1737
1738 } /* end of loop on zone types (boundary/volume) */
1739
1740 BFT_FREE(elt_particle_idx);
1741
1742 /* Update global particle counters */
1743
1744 pc = cs_lagr_update_particle_counter();
1745 pc->n_g_total += pc->n_g_new;
1746 }
1747
1748 /*----------------------------------------------------------------------------*/
1749
1750 END_C_DECLS
1751