1 /*============================================================================
2 * Methods for particle parameters
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 /*============================================================================
28 * Functions dealing with lagrangian initialization
29 *============================================================================*/
30
31 #include "cs_defs.h"
32
33 /*----------------------------------------------------------------------------
34 * Standard C library headers
35 *----------------------------------------------------------------------------*/
36
37 #include <limits.h>
38 #include <stdio.h>
39 #include <stddef.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <math.h>
43 #include <ctype.h>
44 #include <float.h>
45 #include <assert.h>
46
47 /*----------------------------------------------------------------------------
48 * Local headers
49 *----------------------------------------------------------------------------*/
50
51 #include "bft_mem.h"
52 #include "bft_printf.h"
53
54 #include "cs_array.h"
55 #include "cs_base.h"
56
57 #include "cs_field.h"
58 #include "cs_field_pointer.h"
59
60 #include "cs_math.h"
61 #include "cs_mesh_location.h"
62
63 #include "cs_boundary_conditions.h"
64 #include "cs_boundary_zone.h"
65 #include "cs_volume_zone.h"
66
67 #include "cs_parameters.h"
68 #include "cs_prototypes.h"
69 #include "cs_time_step.h"
70 #include "cs_physical_constants.h"
71 #include "cs_thermal_model.h"
72 #include "cs_turbulence_model.h"
73 #include "cs_physical_model.h"
74 #include "cs_parall.h"
75 #include "cs_post.h"
76 #include "cs_post_default.h"
77 #include "cs_prototypes.h"
78
79 #include "cs_gui_particles.h"
80 #include "cs_gui_util.h"
81
82 #include "cs_lagr.h"
83
84 #include "cs_lagr_lec.h"
85 #include "cs_lagr_geom.h"
86 #include "cs_lagr_dlvo.h"
87 #include "cs_lagr_roughness.h"
88 #include "cs_lagr_clogging.h"
89 #include "cs_lagr_injection.h"
90 #include "cs_lagr_gradients.h"
91 #include "cs_lagr_car.h"
92 #include "cs_lagr_coupling.h"
93 #include "cs_lagr_new.h"
94 #include "cs_lagr_particle.h"
95 #include "cs_lagr_resuspension.h"
96 #include "cs_lagr_stat.h"
97 #include "cs_lagr_tracking.h"
98 #include "cs_lagr_print.h"
99 #include "cs_lagr_poisson.h"
100 #include "cs_lagr_post.h"
101 #include "cs_lagr_sde.h"
102 #include "cs_lagr_sde_model.h"
103 #include "cs_lagr_orientation.h"
104 #include "cs_lagr_prototypes.h"
105 #include "cs_lagr_agglo.h"
106 #include "cs_lagr_fragmentation.h"
107
108 #include "cs_random.h"
109
110 /*----------------------------------------------------------------------------
111 * Header for the current file
112 *----------------------------------------------------------------------------*/
113
114 #include "cs_lagr.h"
115
116 /*----------------------------------------------------------------------------*/
117
118 BEGIN_C_DECLS
119
120 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
121
122 /*============================================================================
123 * Type definitions
124 *============================================================================*/
125
126 /*============================================================================
127 * Local structure definitions
128 *============================================================================*/
129
130 /*=============================================================================
131 * Local Enumeration definitions
132 *============================================================================*/
133
134 /*============================================================================
135 * Static global variables
136 *============================================================================*/
137
138 /* Fixed-size parameters */
139
140 cs_lagr_const_dim_t _lagr_const_dim
141 = {.nusbrd = 10,
142 .ndlaim = 10,
143 .ncharm2 = 5,
144 .nlayer = 5};
145
146 const cs_lagr_const_dim_t *cs_glob_lagr_const_dim
147 = (const cs_lagr_const_dim_t *)(&_lagr_const_dim);
148
149 /* General dimensions */
150
151 cs_lagr_dim_t _lagr_dim = {.ntersl = 0,
152 .n_boundary_stats = 0};
153
154 cs_lagr_dim_t *cs_glob_lagr_dim = &_lagr_dim;
155
156 /* Time and Lagrangian-Eulerian coupling scheme */
157
158 static cs_lagr_time_scheme_t _lagr_time_scheme
159 = {.iilagr = CS_LAGR_OFF,
160 .isttio = 0,
161 .isuila = 1,
162 .t_order = 0,
163 .ilapoi = 0,
164 .iadded_mass = 0,
165 .added_mass_const = 0};
166
167 /* Main Lagrangian physical model parameters */
168
169 static cs_lagr_model_t _lagr_model
170 = {.physical_model = CS_LAGR_PHYS_OFF,
171 .n_temperature_layers = 1,
172 .modcpl = 1,
173 .idistu = -1,
174 .idiffl = -1,
175 .deposition = 0,
176 .dlvo = 0,
177 .roughness = 0,
178 .resuspension = 0,
179 .clogging = 0,
180 .shape = 0,
181 .consolidation = 0,
182 .precipitation = 0,
183 .fouling = 0,
184 .n_stat_classes = 0,
185 .agglomeration = 0,
186 .fragmentation = 0,
187 .n_user_variables = 0};
188
189 /* particle counter structure and associated pointer */
190
191 static cs_lagr_particle_counter_t _lagr_particle_counter
192 = {.n_g_cumulative_total = 0,
193 .n_g_cumulative_failed = 0,
194 .n_g_total = 0,
195 .n_g_new = 0,
196 .n_g_exit = 0,
197 .n_g_merged = 0,
198 .n_g_deposited = 0,
199 .n_g_fouling = 0,
200 .n_g_resuspended = 0,
201 .n_g_failed = 0,
202 .w_total = 0.,
203 .w_new = 0.,
204 .w_exit = 0.,
205 .w_merged = 0.,
206 .w_deposited = 0.,
207 .w_fouling = 0.,
208 .w_resuspended = 0.};
209
210 /* lagr specific physics structure and associated pointer */
211 static cs_lagr_specific_physics_t _cs_glob_lagr_specific_physics
212 = {0, 0, 0, 0, 0};
213 cs_lagr_specific_physics_t *cs_glob_lagr_specific_physics
214 = &_cs_glob_lagr_specific_physics;
215
216 /* lagr reentrained model structure and associated pointer */
217 static cs_lagr_reentrained_model_t _lagr_reentrained_model
218 = {.ireent = 0,
219 .iflow = 0,
220 .espasg = 0.,
221 .denasp = 0.,
222 .modyeq = 0.,
223 .rayasp = 0.,
224 .rayasg = 0.};
225 cs_lagr_reentrained_model_t *cs_glob_lagr_reentrained_model
226 = &_lagr_reentrained_model;
227
228 /* lagr precipitation model structure and associated pointer */
229 static cs_lagr_precipitation_model_t _cs_glob_lagr_precipitation_model
230 = {0, 0, 0, NULL, NULL, NULL};
231 cs_lagr_precipitation_model_t *cs_glob_lagr_precipitation_model
232 = &_cs_glob_lagr_precipitation_model;
233
234 /* lagr clogging model structure and associated pointer */
235 static cs_lagr_clogging_model_t _cs_glob_lagr_clogging_model = {0, 0, 0, 0};
236 cs_lagr_clogging_model_t *cs_glob_lagr_clogging_model
237 = &_cs_glob_lagr_clogging_model;
238
239 /* lagr non-spherical model structure and associated pointer */
240 static cs_lagr_shape_model_t _cs_glob_lagr_shape_model = {0};
241 cs_lagr_shape_model_t *cs_glob_lagr_shape_model
242 = &_cs_glob_lagr_shape_model;
243
244 /* lagr agglomeration model structure and associated pointer */
245 static cs_lagr_agglomeration_model_t _cs_glob_lagr_agglomeration_model
246 = {
247 .n_max_classes = 10000,
248 .min_stat_weight = 0.,
249 .max_stat_weight = 0.,
250 .scalar_kernel = 0.,
251 .base_diameter = 0.};
252
253 cs_lagr_agglomeration_model_t *cs_glob_lagr_agglomeration_model
254 = &_cs_glob_lagr_agglomeration_model;
255
256 /* lagr fragmentation model structure and associated pointer */
257 static cs_lagr_fragmentation_model_t _cs_glob_lagr_fragmentation_model
258 = {
259 .scalar_kernel = 0.,
260 .base_diameter = 0.,
261 NULL};
262
263 cs_lagr_fragmentation_model_t *cs_glob_lagr_fragmentation_model
264 = &_cs_glob_lagr_fragmentation_model;
265
266 /* lagr consolidation model structure and associated pointer */
267 static cs_lagr_consolidation_model_t _cs_glob_lagr_consolidation_model
268 = {0, 0, 0, 0};
269 cs_lagr_consolidation_model_t *cs_glob_lagr_consolidation_model
270 = &_cs_glob_lagr_consolidation_model;
271
272 /*! current time step status */
273
274 static cs_lagr_time_step_t _cs_glob_lagr_time_step
275 = {.nor = 0,
276 .dtp = 0.,
277 .ttclag = 0.};
278
279 /* lagr source terms structure and associated pointer */
280 static cs_lagr_source_terms_t _cs_glob_lagr_source_terms
281 = {.ltsdyn = 0,
282 .ltsmas = 0,
283 .ltsthe = 0,
284 .itsli = 0,
285 .itske = 0,
286 .itste = 0,
287 .itsti = 0,
288 .itsmas = 0,
289 .itsmv1 = NULL,
290 .itsmv2 = NULL,
291 .itsco = 0,
292 .itsfp4 = 0,
293 .nstits = 0,
294 .npts = 0,
295 .ntxerr = 0,
296 .vmax = 0,
297 .tmamax = 0,
298 .st_val = NULL};
299
300 cs_lagr_source_terms_t *cs_glob_lagr_source_terms
301 = &_cs_glob_lagr_source_terms;
302
303 /* lagr encrustation structure and associated pointer */
304 static cs_lagr_encrustation_t _cs_glob_lagr_encrustation
305 = {0, 0, NULL, NULL, NULL, NULL, 0} ;
306 cs_lagr_encrustation_t *cs_glob_lagr_encrustation
307 = &_cs_glob_lagr_encrustation;
308
309 /* lagr physico chemical structure and associated pointer */
310 static cs_lagr_physico_chemical_t _cs_glob_lagr_physico_chemical
311 = {0, 0, 0, 0, 0, 0, 0};
312 cs_lagr_physico_chemical_t *cs_glob_lagr_physico_chemical
313 = &_cs_glob_lagr_physico_chemical;
314
315 /* lagr brownian structure and associated pointer */
316 static cs_lagr_brownian_t _cs_glob_lagr_brownian = {0};
317 cs_lagr_brownian_t *cs_glob_lagr_brownian = &_cs_glob_lagr_brownian;
318
319 /* lagr boundary interactions structure and associated pointer */
320 static cs_lagr_boundary_interactions_t _cs_glob_lagr_boundary_interactions
321 = {.npstf = 0,
322 .npstft = 0,
323 .has_part_impact_nbr = 0,
324 .iclgst = 0,
325 .inbr = -1,
326 .inclg = -1,
327 .inclgt = -1,
328 .iclogt = -1,
329 .iscovc = -1,
330 .ihdepm = -1,
331 .ihdiam = -1,
332 .ihsum = -1,
333 .tstatp = 0.,
334 .nombrd = NULL};
335
336 cs_lagr_boundary_interactions_t *cs_glob_lagr_boundary_interactions
337 = &_cs_glob_lagr_boundary_interactions;
338
339 /* lagr extra modules and associated pointer */
340
341 static cs_lagr_extra_module_t _lagr_extra_module
342 = {.iturb = 0,
343 .itytur = 0,
344 .ncharb = 0,
345 .ncharm = 0,
346 .radiative_model = 0,
347 .icp = -1,
348 .diftl0 = 0,
349 .cmu = 0,
350 .visls0 = 0,
351 .ustar = NULL,
352 .cromf = NULL,
353 .pressure = NULL,
354 .scal_t = NULL,
355 .temperature = NULL,
356 .vel = NULL,
357 .viscl = NULL,
358 .cpro_viscls = NULL,
359 .cpro_cp = NULL,
360 .luminance = NULL,
361 .x_oxyd = NULL,
362 .x_eau = NULL,
363 .x_m = NULL,
364 .cvar_k = NULL,
365 .cvar_ep = NULL,
366 .cvar_omg = NULL,
367 .cvar_rij = NULL,
368 .grad_pr = NULL,
369 .grad_vel = NULL};
370
371 cs_lagr_extra_module_t *cs_glob_lagr_extra_module = &_lagr_extra_module;
372
373 /* lagr coal combustion structure and associated pointer */
374
375 static cs_lagr_coal_comb_t _lagr_coal_comb
376 = {0, 0, 0,
377 0, 0.0 , 0.0,
378 0, NULL,
379 0, NULL, NULL,
380 0, NULL, NULL, NULL, NULL,
381 NULL, NULL, NULL, NULL, NULL,
382 NULL, NULL, NULL, NULL, NULL};
383
384 /* boundary and volume condition data */
385
386 static cs_lagr_zone_data_t *_boundary_conditions = NULL;
387 static cs_lagr_zone_data_t *_volume_conditions = NULL;
388
389 /*============================================================================
390 * Global variables
391 *============================================================================*/
392
393 /* Pointers to global structures; some may be const-qualified in the future
394 to set to read-only ouside accessor functions, but this is not done
395 at this stage, as the API is not yet stabilized */
396
397 cs_lagr_time_scheme_t *cs_glob_lagr_time_scheme = &_lagr_time_scheme;
398 cs_lagr_model_t *cs_glob_lagr_model = &_lagr_model;
399
400 const cs_lagr_particle_counter_t *cs_glob_lagr_particle_counter
401 = &_lagr_particle_counter;
402
403 int cs_glob_lagr_log_frequency_n = 1; /* log every frequency_n time steps */
404
405 cs_lagr_time_step_t *cs_glob_lagr_time_step = &_cs_glob_lagr_time_step;
406
407 cs_lagr_coal_comb_t *cs_glob_lagr_coal_comb = &_lagr_coal_comb;
408
409 cs_real_t *bound_stat = NULL;
410
411 const cs_lagr_zone_data_t *cs_glob_lagr_boundary_conditions = NULL;
412 const cs_lagr_zone_data_t *cs_glob_lagr_volume_conditions = NULL;
413
414 cs_lagr_internal_condition_t *cs_glob_lagr_internal_conditions = NULL;
415
416 /* Geometry helper arrays */
417 /*------------------------*/
418
419 /*! Projection matrices for global to local coordinates on boundary faces */
420
421 cs_real_33_t *cs_glob_lagr_b_face_proj = NULL;
422
423 /*============================================================================
424 * Prototypes for functions intended for use only by Fortran wrappers.
425 * (descriptions follow, with function bodies).
426 *============================================================================*/
427
428 void
429 cs_f_lagr_params_pointers(int **p_iilagr,
430 int **p_idepst,
431 int **p_iflow,
432 int **p_ipreci);
433
434 void
435 cs_f_lagr_dim_pointers(int **p_ntersl);
436
437 void
438 cs_f_lagr_clogging_model_pointers(cs_real_t **jamlim,
439 cs_real_t **mporos,
440 cs_real_t **csthpp);
441
442 void
443 cs_f_lagr_shape_model_pointers(cs_real_t **param_chmb);
444
445 void
446 cs_f_lagr_agglomeration_model_pointers( cs_lnum_t **n_max_classes,
447 cs_real_t **min_stat_weight,
448 cs_real_t **max_stat_weight,
449 cs_real_t **scalar_kernel,
450 cs_real_t **base_diameter );
451
452 void
453 cs_f_lagr_consolidation_model_pointers(cs_lnum_t **iconsol,
454 cs_real_t **rate_consol,
455 cs_real_t **slope_consol,
456 cs_real_t **force_consol);
457
458 void
459 cs_f_lagr_source_terms_pointers(int **p_ltsdyn,
460 int **p_ltsmas,
461 int **p_ltsthe,
462 int **p_itsli,
463 int **p_itske,
464 int **p_itste,
465 int **p_itsti,
466 int **p_itsmas,
467 int **p_itsco,
468 int **p_itsmv1,
469 int **p_itsmv2,
470 int *dim_itsmv1,
471 int *dim_itsmv2);
472
473 void
474 cs_f_lagr_specific_physics(int *iirayo,
475 int *ncharb,
476 int *ncharm,
477 cs_real_t *diftl0);
478
479 void
480 cs_f_lagr_coal_comb(int *ih2o,
481 int *io2,
482 int *ico,
483 int *iatc,
484 cs_real_t *prefth,
485 cs_real_t *trefth,
486 int *natom,
487 cs_real_t *wmolat,
488 int *ngazem,
489 cs_real_t *wmole,
490 int *iym1,
491 int *ncharm,
492 cs_real_t *a1ch,
493 cs_real_t *h02ch,
494 cs_real_t *e1ch,
495 cs_real_t *a2ch,
496 cs_real_t *e2ch,
497 cs_real_t *y1ch,
498 cs_real_t *y2ch,
499 cs_real_t *cp2ch,
500 cs_real_t *ahetch,
501 cs_real_t *ehetch,
502 cs_real_t *rho0ch,
503 cs_real_t *xwatch,
504 cs_real_t *xashch,
505 cs_real_t *thcdch);
506
507 /*============================================================================
508 * Fortran wrapper function definitions
509 *============================================================================*/
510
511 /*----------------------------------------------------------------------------
512 * Map some options
513 *
514 * This function is intended for use by Fortran wrappers.
515 *
516 * parameters:
517 * p_iilagr --> lagrangian model type
518 * p_idepo --> deposition option flag
519 * p_iflow -->
520 * p_ipreci --> precipitation option flag
521 *----------------------------------------------------------------------------*/
522
523 void
cs_f_lagr_params_pointers(int ** p_iilagr,int ** p_idepst,int ** p_iflow,int ** p_ipreci)524 cs_f_lagr_params_pointers(int **p_iilagr,
525 int **p_idepst,
526 int **p_iflow,
527 int **p_ipreci)
528 {
529 *p_iilagr = &_lagr_time_scheme.iilagr;
530 *p_idepst = &_lagr_model.deposition;
531 *p_iflow= &_lagr_reentrained_model.iflow;
532 *p_ipreci = &_lagr_model.precipitation;
533 }
534
535 /*----------------------------------------------------------------------------
536 * Get pointers to members of the global lagr dim structure.
537 *
538 * This function is intended for use by Fortran wrappers, and
539 * enables mapping to Fortran global pointers.
540 *
541 * parameters:
542 * ntersl --> pointer to cs_glob_lagr_dim->ntersl
543 *----------------------------------------------------------------------------*/
544
545 void
cs_f_lagr_dim_pointers(int ** p_ntersl)546 cs_f_lagr_dim_pointers(int **p_ntersl)
547 {
548 *p_ntersl = &(_lagr_dim.ntersl);
549 }
550
551 void
cs_f_lagr_clogging_model_pointers(cs_real_t ** jamlim,cs_real_t ** mporos,cs_real_t ** csthpp)552 cs_f_lagr_clogging_model_pointers(cs_real_t **jamlim,
553 cs_real_t **mporos,
554 cs_real_t **csthpp)
555 {
556 *jamlim = &cs_glob_lagr_clogging_model->jamlim;
557 *mporos = &cs_glob_lagr_clogging_model->mporos;
558 *csthpp = &cs_glob_lagr_clogging_model->csthpp;
559 }
560
561 void
cs_f_lagr_shape_model_pointers(cs_real_t ** param_chmb)562 cs_f_lagr_shape_model_pointers(cs_real_t **param_chmb)
563 {
564 *param_chmb = &cs_glob_lagr_shape_model->param_chmb;
565 }
566
567 void
cs_f_lagr_agglomeration_model_pointers(cs_lnum_t ** n_max_classes,cs_real_t ** min_stat_weight,cs_real_t ** max_stat_weight,cs_real_t ** scalar_kernel,cs_real_t ** base_diameter)568 cs_f_lagr_agglomeration_model_pointers(cs_lnum_t **n_max_classes,
569 cs_real_t **min_stat_weight,
570 cs_real_t **max_stat_weight,
571 cs_real_t **scalar_kernel,
572 cs_real_t **base_diameter )
573 {
574 *n_max_classes = &cs_glob_lagr_agglomeration_model->n_max_classes;
575 *min_stat_weight = &cs_glob_lagr_agglomeration_model->min_stat_weight;
576 *max_stat_weight = &cs_glob_lagr_agglomeration_model->max_stat_weight;
577 *scalar_kernel = &cs_glob_lagr_agglomeration_model->scalar_kernel;
578 *base_diameter = &cs_glob_lagr_agglomeration_model->base_diameter;
579 }
580
581
582 void
cs_f_lagr_consolidation_model_pointers(cs_lnum_t ** iconsol,cs_real_t ** rate_consol,cs_real_t ** slope_consol,cs_real_t ** force_consol)583 cs_f_lagr_consolidation_model_pointers(cs_lnum_t **iconsol,
584 cs_real_t **rate_consol,
585 cs_real_t **slope_consol,
586 cs_real_t **force_consol)
587 {
588 *iconsol = &cs_glob_lagr_consolidation_model->iconsol;
589 *rate_consol = &cs_glob_lagr_consolidation_model->rate_consol;
590 *slope_consol = &cs_glob_lagr_consolidation_model->slope_consol;
591 *force_consol = &cs_glob_lagr_consolidation_model->force_consol;
592 }
593
594 void
cs_f_lagr_source_terms_pointers(int ** p_ltsdyn,int ** p_ltsmas,int ** p_ltsthe,int ** p_itsli,int ** p_itske,int ** p_itste,int ** p_itsti,int ** p_itsmas,int ** p_itsco,int ** p_itsmv1,int ** p_itsmv2,int * dim_itsmv1,int * dim_itsmv2)595 cs_f_lagr_source_terms_pointers(int **p_ltsdyn,
596 int **p_ltsmas,
597 int **p_ltsthe,
598 int **p_itsli,
599 int **p_itske,
600 int **p_itste,
601 int **p_itsti,
602 int **p_itsmas,
603 int **p_itsco,
604 int **p_itsmv1,
605 int **p_itsmv2,
606 int *dim_itsmv1,
607 int *dim_itsmv2)
608 {
609 *p_ltsdyn = &cs_glob_lagr_source_terms->ltsdyn;
610 *p_ltsmas = &cs_glob_lagr_source_terms->ltsmas;
611 *p_ltsthe = &cs_glob_lagr_source_terms->ltsthe;
612 *p_itsli = &cs_glob_lagr_source_terms->itsli;
613 *p_itske = &cs_glob_lagr_source_terms->itske;
614 *p_itste = &cs_glob_lagr_source_terms->itste;
615 *p_itsti = &cs_glob_lagr_source_terms->itsti;
616 *p_itsmas = &cs_glob_lagr_source_terms->itsmas;
617 *p_itsco = &cs_glob_lagr_source_terms->itsco;
618
619 if (cs_glob_lagr_source_terms->itsmv1 == NULL)
620 BFT_MALLOC(cs_glob_lagr_source_terms->itsmv1,
621 cs_glob_lagr_const_dim->ncharm2, int);
622 *p_itsmv1 = cs_glob_lagr_source_terms->itsmv1;
623 *dim_itsmv1 = cs_glob_lagr_const_dim->ncharm2;
624
625 if (cs_glob_lagr_source_terms->itsmv2 == NULL)
626 BFT_MALLOC(cs_glob_lagr_source_terms->itsmv2,
627 cs_glob_lagr_const_dim->ncharm2, int);
628 *p_itsmv2 = cs_glob_lagr_source_terms->itsmv2;
629 *dim_itsmv2 = cs_glob_lagr_const_dim->ncharm2;
630 }
631
632 void
cs_f_lagr_specific_physics(int * iirayo,int * ncharb,int * ncharm,cs_real_t * diftl0)633 cs_f_lagr_specific_physics(int *iirayo,
634 int *ncharb,
635 int *ncharm,
636 cs_real_t *diftl0)
637 {
638 cs_turb_model_t *turb_model = cs_get_glob_turb_model();
639
640 if (turb_model == NULL)
641 bft_error(__FILE__, __LINE__, 0,
642 "%s: Turbulence modelling is not set.", __func__);
643
644 _lagr_extra_module.iturb = turb_model->iturb;
645 _lagr_extra_module.itytur = turb_model->itytur;
646 _lagr_extra_module.ncharb = *ncharb;
647 _lagr_extra_module.ncharm = *ncharm;
648 _lagr_extra_module.icp = cs_glob_fluid_properties->icp;
649
650 _lagr_extra_module.radiative_model = *iirayo;
651 _lagr_extra_module.diftl0 = *diftl0;
652 _lagr_extra_module.cmu = cs_turb_cmu;
653 }
654
655 void
cs_f_lagr_coal_comb(int * ih2o,int * io2,int * ico,int * iatc,cs_real_t * prefth,cs_real_t * trefth,int * natom,cs_real_t * wmolat,int * ngazem,cs_real_t * wmole,int * iym1,int * ncharm,cs_real_t * a1ch,cs_real_t * h02ch,cs_real_t * e1ch,cs_real_t * a2ch,cs_real_t * e2ch,cs_real_t * y1ch,cs_real_t * y2ch,cs_real_t * cp2ch,cs_real_t * ahetch,cs_real_t * ehetch,cs_real_t * rho0ch,cs_real_t * xwatch,cs_real_t * xashch,cs_real_t * thcdch)656 cs_f_lagr_coal_comb(int *ih2o,
657 int *io2,
658 int *ico,
659 int *iatc,
660 cs_real_t *prefth,
661 cs_real_t *trefth,
662 int *natom,
663 cs_real_t *wmolat,
664 int *ngazem,
665 cs_real_t *wmole,
666 int *iym1,
667 int *ncharm,
668 cs_real_t *a1ch,
669 cs_real_t *h02ch,
670 cs_real_t *e1ch,
671 cs_real_t *a2ch,
672 cs_real_t *e2ch,
673 cs_real_t *y1ch,
674 cs_real_t *y2ch,
675 cs_real_t *cp2ch,
676 cs_real_t *ahetch,
677 cs_real_t *ehetch,
678 cs_real_t *rho0ch,
679 cs_real_t *xwatch,
680 cs_real_t *xashch,
681 cs_real_t *thcdch)
682 {
683 cs_glob_lagr_coal_comb->ih2o = *ih2o;
684 cs_glob_lagr_coal_comb->io2 = *io2;
685 cs_glob_lagr_coal_comb->ico = *ico;
686
687 cs_glob_lagr_coal_comb->iatc = *iatc;
688 cs_glob_lagr_coal_comb->prefth = *prefth;
689 cs_glob_lagr_coal_comb->trefth = *trefth;
690
691 cs_glob_lagr_coal_comb->natom = *natom;
692 cs_glob_lagr_coal_comb->wmolat = wmolat;
693
694 cs_glob_lagr_coal_comb->ngazem = *ngazem;
695 cs_glob_lagr_coal_comb->wmole = wmole;
696 cs_glob_lagr_coal_comb->iym1 = iym1;
697
698 cs_glob_lagr_coal_comb->ncharm = *ncharm;
699 cs_glob_lagr_coal_comb->a1ch = a1ch;
700 cs_glob_lagr_coal_comb->h02ch = h02ch;
701 cs_glob_lagr_coal_comb->e1ch = e1ch;
702 cs_glob_lagr_coal_comb->a2ch = a2ch;
703 cs_glob_lagr_coal_comb->e2ch = e2ch;
704 cs_glob_lagr_coal_comb->y1ch = y1ch;
705 cs_glob_lagr_coal_comb->y2ch = y2ch;
706 cs_glob_lagr_coal_comb->cp2ch = cp2ch;
707 cs_glob_lagr_coal_comb->ahetch = ahetch;
708 cs_glob_lagr_coal_comb->ehetch = ehetch;
709 cs_glob_lagr_coal_comb->rho0ch = rho0ch;
710 cs_glob_lagr_coal_comb->xwatch = xwatch;
711 cs_glob_lagr_coal_comb->xashch = xashch;
712 cs_glob_lagr_coal_comb->thcdch = thcdch;
713 }
714
715 /*=============================================================================
716 * Private function definitions
717 *============================================================================*/
718
719 static void
_lagr_map_fields_default(void)720 _lagr_map_fields_default(void)
721 {
722 if ( cs_glob_physical_model_flag[CS_COMBUSTION_COAL] >= 0
723 || cs_glob_physical_model_flag[CS_COMBUSTION_FUEL] >= 0) {
724 _lagr_extra_module.cromf = cs_field_by_name_try("rho_gas");
725 }
726 else {
727 _lagr_extra_module.cromf = cs_field_by_name_try("density");
728 }
729
730 _lagr_extra_module.pressure = cs_field_by_name_try("pressure");
731
732 _lagr_extra_module.luminance = cs_field_by_name_try("luminance");
733 if (cs_field_by_name_try("velocity_1") != NULL) {
734 /* we are probably using NEPTUNE_CFD */
735 _lagr_extra_module.vel = cs_field_by_name_try("lagr_velocity");
736
737 _lagr_extra_module.cvar_k = cs_field_by_name_try("lagr_k");
738 _lagr_extra_module.cvar_ep = cs_field_by_name_try("lagr_epsilon");
739 _lagr_extra_module.cvar_omg = NULL;
740 _lagr_extra_module.cvar_rij = cs_field_by_name_try("lagr_rij");
741 _lagr_extra_module.viscl = cs_field_by_name_try
742 ("lagr_molecular_viscosity");
743 _lagr_extra_module.scal_t = cs_field_by_name_try("lagr_enthalpy");
744 _lagr_extra_module.cpro_viscls = cs_field_by_name_try
745 ("lagr_thermal_conductivity");
746 _lagr_extra_module.cpro_cp = cs_field_by_name_try("lagr_specific_heat");
747 _lagr_extra_module.temperature = cs_field_by_name_try("lagr_temperature");
748 _lagr_extra_module.x_oxyd = NULL;
749 _lagr_extra_module.x_eau = NULL;
750 _lagr_extra_module.x_m = NULL;
751 _lagr_extra_module.cromf = cs_field_by_name_try("lagr_density");
752 /* TODO FIXME */
753 _lagr_extra_module.visls0 = 0.;
754
755 _lagr_extra_module.ustar
756 = cs_field_by_name_try("lagr_wall_friction_velocity");
757 }
758 else {
759 /* we use Code_Saturne */
760 _lagr_extra_module.vel = cs_field_by_name_try("velocity");
761 _lagr_extra_module.cvar_k = cs_field_by_name_try("k");
762 _lagr_extra_module.cvar_ep = cs_field_by_name_try("epsilon");
763 _lagr_extra_module.cvar_omg = cs_field_by_name_try("omega");
764 _lagr_extra_module.cvar_rij = cs_field_by_name_try("rij");
765 _lagr_extra_module.viscl = cs_field_by_name_try("molecular_viscosity");
766 _lagr_extra_module.cpro_viscls = NULL;
767
768 _lagr_extra_module.scal_t = cs_thermal_model_field();
769
770 if (_lagr_extra_module.scal_t != NULL) {
771 _lagr_extra_module.visls0
772 = cs_field_get_key_double(_lagr_extra_module.scal_t,
773 cs_field_key_id("diffusivity_ref"));
774
775 int l_id = cs_field_get_key_int(_lagr_extra_module.scal_t,
776 cs_field_key_id("diffusivity_id"));
777 if (l_id >= 0)
778 _lagr_extra_module.cpro_viscls = cs_field_by_id(l_id);
779 }
780
781 _lagr_extra_module.cpro_cp = cs_field_by_name_try("specific_heat");
782 _lagr_extra_module.temperature = cs_field_by_name_try("temperature");
783
784 _lagr_extra_module.x_oxyd = cs_field_by_name_try("ym_o2");
785 _lagr_extra_module.x_eau = cs_field_by_name_try("ym_h2o");
786 _lagr_extra_module.x_m = cs_field_by_name_try("xm");
787
788 _lagr_extra_module.ustar = cs_field_by_name_try("ustar");
789 }
790 }
791
792 /*----------------------------------------------------------------------------*/
793 /*!
794 * \brief Reallocate a zone injection data structure for a given zone and set.
795 *
796 * The data is reallocated (expanded) if needed, and section relative
797 * to the given class id is set to default values.
798 *
799 * Note that the structure is allocated in one block, large enough to contain
800 * the base structure and model-dependent arrays it may point to.
801 *
802 * \param[in] location_id id of associated location
803 * \param[in] zone_id id of associated zone
804 * \param[in] set_id id of requested class
805 * \param[in, out] n_injection_sets number of injection sets for this zone
806 * \param[in, out] injection_sets zone injection data for this zone
807 */
808 /*----------------------------------------------------------------------------*/
809
810 static void
_zone_injection_set_init(int location_id,int zone_id,int set_id,int * n_injection_sets,cs_lagr_injection_set_t ** injection_sets)811 _zone_injection_set_init(int location_id,
812 int zone_id,
813 int set_id,
814 int *n_injection_sets,
815 cs_lagr_injection_set_t **injection_sets)
816 {
817 /* Simply reset if allocation is large enough */
818
819 if (set_id < *n_injection_sets) {
820 cs_lagr_injection_set_t *zis
821 = (cs_lagr_injection_set_t *)((*injection_sets) + set_id);
822 cs_lagr_injection_set_default(zis);
823 }
824
825 /* Reallocate memory so as to maintain
826 sub-arrays at the end of the structure */
827
828 else {
829
830 cs_lagr_injection_set_t *_zis = *injection_sets;
831 BFT_REALLOC(_zis, set_id+1, cs_lagr_injection_set_t);
832
833 for (int i = *n_injection_sets; i <= set_id; i++) {
834
835 cs_lagr_injection_set_t *zis = (cs_lagr_injection_set_t *)(_zis + i);
836
837 memset(zis, 0, sizeof(cs_lagr_injection_set_t));
838
839 zis->location_id = location_id;
840 zis->zone_id = zone_id;
841 zis->set_id = set_id;
842
843 cs_lagr_injection_set_default(zis);
844
845 }
846
847 *n_injection_sets = set_id+1;
848 *injection_sets = _zis;
849
850 }
851 }
852
853 /*----------------------------------------------------------------------------
854 * Initialize or update a zone data structure
855 *
856 * parameters:
857 * zone_data <-> pointer to zone data structure pointer
858 * location_id <-- mesh location id
859 * n_zones <-- number of zones
860 *
861 * returns:
862 * a new defined cs_lagr_injection_sets_t structure
863 *----------------------------------------------------------------------------*/
864
865 static void
_update_zone_data_struct(cs_lagr_zone_data_t ** zone_data,int location_id,int n_zones)866 _update_zone_data_struct(cs_lagr_zone_data_t **zone_data,
867 int location_id,
868 int n_zones)
869 {
870 cs_lagr_zone_data_t *zd = *zone_data;
871
872 if (*zone_data == NULL) {
873 BFT_MALLOC(zd, 1, cs_lagr_zone_data_t);
874 zd->location_id = location_id;
875 zd->n_zones = 0;
876 zd->zone_type = NULL;
877 zd->n_injection_sets = NULL;
878 zd->injection_set = NULL;
879 zd->elt_type = NULL;
880 zd->particle_flow_rate = NULL;
881 *zone_data = zd;
882 }
883
884 if (zd->n_zones < n_zones) {
885 int n_stats = cs_glob_lagr_model->n_stat_classes + 1;
886 BFT_REALLOC(zd->zone_type, n_zones, int);
887 BFT_REALLOC(zd->n_injection_sets, n_zones, int);
888 BFT_REALLOC(zd->injection_set, n_zones, cs_lagr_injection_set_t *);
889 BFT_REALLOC(zd->particle_flow_rate, n_zones*n_stats, cs_real_t);
890 for (int i = zd->n_zones; i < n_zones; i++) {
891 zd->zone_type[i] = -1;
892 zd->n_injection_sets[i] = 0;
893 zd->injection_set[i] = NULL;
894 }
895 for (int i = zd->n_zones*n_stats; i < n_zones*n_stats; i++)
896 zd->particle_flow_rate[i] = 0;
897
898 zd->n_zones = n_zones;
899 }
900 }
901
902 /*----------------------------------------------------------------------------
903 * Initialize a cs_lagr_internal_condition_t structure.
904 *
905 * returns:
906 * a new defined cs_lagr_internal_condition_t structure
907 *----------------------------------------------------------------------------*/
908
909 static cs_lagr_internal_condition_t *
_create_internal_cond_struct(void)910 _create_internal_cond_struct(void)
911 {
912 cs_lagr_internal_condition_t *internal_cond = NULL;
913 cs_mesh_t *mesh = cs_glob_mesh;
914
915 BFT_MALLOC(internal_cond, 1, cs_lagr_internal_condition_t);
916
917 BFT_MALLOC(internal_cond->i_face_zone_id, mesh->n_i_faces, int);
918
919 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++)
920 internal_cond->i_face_zone_id[i] = -1;
921
922 return internal_cond;
923 }
924
925 /*----------------------------------------------------------------------------*/
926 /*!
927 * \brief Update (or build) boundary face types.
928 */
929 /*----------------------------------------------------------------------------*/
930
931 static void
_update_boundary_face_type(void)932 _update_boundary_face_type(void)
933 {
934 cs_lagr_zone_data_t *bcs = cs_lagr_get_boundary_conditions();
935
936 const cs_mesh_t *mesh = cs_glob_mesh;
937
938 BFT_REALLOC(bcs->elt_type, mesh->n_b_faces, char);
939
940 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++)
941 bcs->elt_type[i] = 0;
942
943 for (int z_id = 0; z_id < bcs->n_zones; z_id++) {
944
945 if (bcs->zone_type[z_id] < 0) /* ignore undefined zones */
946 continue;
947
948 char z_type = bcs->zone_type[z_id];
949
950 const cs_zone_t *z = cs_boundary_zone_by_id(z_id);
951 for (cs_lnum_t i = 0; i < z->n_elts; i++)
952 bcs->elt_type[z->elt_ids[i]] = z_type;
953
954 }
955
956 /* Check definitions */
957
958 {
959 int *bc_flag;
960 BFT_MALLOC(bc_flag, mesh->n_b_faces, int);
961
962 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++)
963 bc_flag[i] = bcs->elt_type[i];
964
965 cs_boundary_conditions_error(bc_flag, _("Lagrangian boundary type"));
966
967 BFT_FREE(bc_flag);
968 }
969 }
970
971 /*----------------------------------------------------------------------------*/
972 /*!
973 * \brief Obtain the number of mesh cells occupied by at least one particle.
974 *
975 * \param[in] p_set pointer to particle data structure
976 * \param[in] start start position in p_set
977 * \param[in] end end position in p_set
978 *
979 * \return
980 * integer that gives the number of cells occupied by at least one particle
981 * in a particle sub-set
982 */
983 /*----------------------------------------------------------------------------*/
984
985 static cs_lnum_t
_get_n_occupied_cells(const cs_lagr_particle_set_t * p_set,cs_lnum_t start,cs_lnum_t end)986 _get_n_occupied_cells(const cs_lagr_particle_set_t *p_set,
987 cs_lnum_t start,
988 cs_lnum_t end)
989 {
990 /*Initialization:
991 counter for number of cells that contain particles */
992 cs_lnum_t counter_particle_cells = 0;
993
994 /* Main loop to update the counter */
995 if (end - start > 1) {
996 counter_particle_cells = 1;
997 cs_lnum_t prev_cell_id
998 = cs_lagr_particles_get_lnum(p_set, start, CS_LAGR_CELL_ID);
999 cs_lnum_t curr_cell_id = prev_cell_id;
1000 for (cs_lnum_t p = start + 1; p < end; ++p) {
1001 curr_cell_id = cs_lagr_particles_get_lnum(p_set, p, CS_LAGR_CELL_ID);
1002 if (prev_cell_id != curr_cell_id)
1003 counter_particle_cells++;
1004 prev_cell_id = curr_cell_id;
1005 }
1006 }
1007 else if (p_set->n_particles == 1) {
1008 counter_particle_cells = 1;
1009 }
1010
1011 return counter_particle_cells;
1012 }
1013
1014 /*----------------------------------------------------------------------------*/
1015 /*!
1016 * \brief Obtain the index of cells that are occupied by at least one
1017 * particle and the list of indices of the particle set that contain
1018 * the first element in a cell
1019 *
1020 * The two arrays that contain this information need to be pre-allocated to
1021 * size n_occupied_cells and n_occupied_cells+1 respectively.
1022 *
1023 * The value of particle_gaps at index i contains the first particle in
1024 * in cell[i+1]. The last element of particle_gaps contains the size of p_set.
1025 *
1026 * \param[in] p_set id of associated location
1027 * \param[in] start id of associated zone
1028 * \param[in] end id of requested class
1029 * \param[in] n_occupied_cells number of cells that are occupied by particles
1030 * \param[out] occupied_cell_ids occupied_cell ids
1031 * \param[out] particle_gaps preallocated list of size n_occupied_cells+1
1032 * that will contain the starting indices of
1033 * particles that are in the current cell
1034 */
1035 /*----------------------------------------------------------------------------*/
1036
1037 static void
_occupied_cells(cs_lagr_particle_set_t * p_set,cs_lnum_t start,cs_lnum_t end,cs_lnum_t n_occupied_cells,cs_lnum_t occupied_cell_ids[],cs_lnum_t particle_gaps[])1038 _occupied_cells(cs_lagr_particle_set_t *p_set,
1039 cs_lnum_t start,
1040 cs_lnum_t end,
1041 cs_lnum_t n_occupied_cells,
1042 cs_lnum_t occupied_cell_ids[],
1043 cs_lnum_t particle_gaps[])
1044 {
1045 if (end - start >= 1) {
1046 /* Initialization */
1047 cs_lnum_t prev_cell_id
1048 = cs_lagr_particles_get_lnum(p_set, start, CS_LAGR_CELL_ID);
1049 cs_lnum_t curr_cell_id = prev_cell_id;
1050 cs_lnum_t counter = 0;
1051 occupied_cell_ids[0] = curr_cell_id;
1052 particle_gaps[0] = 0;
1053
1054 counter = 1;
1055
1056 /* Update lists */
1057 for (cs_lnum_t part = start + 1; part < end; ++part) {
1058 curr_cell_id
1059 = cs_lagr_particles_get_lnum(p_set, part, CS_LAGR_CELL_ID);
1060 if (prev_cell_id != curr_cell_id) {
1061 occupied_cell_ids[counter] = curr_cell_id;
1062 particle_gaps[counter] = part;
1063 counter++;
1064 }
1065 prev_cell_id = curr_cell_id;
1066 }
1067 particle_gaps[n_occupied_cells] = p_set->n_particles;
1068 }
1069 }
1070
1071 /*----------------------------------------------------------------------------*/
1072 /*!
1073 * \brief Obtain the number of particles to be deleted
1074 *
1075 * \param[in] p_set pointer to particle data structure
1076 * \param[in] start start position in p_set
1077 * \param[in] end end position in p_se
1078 *
1079 * \return:
1080 * integer that gives the number of particles to be deleted in a particle sub-set
1081 */
1082 /*----------------------------------------------------------------------------*/
1083
1084 static cs_lnum_t
_get_n_deleted(cs_lagr_particle_set_t * p_set,cs_lnum_t start,cs_lnum_t end)1085 _get_n_deleted(cs_lagr_particle_set_t *p_set,
1086 cs_lnum_t start,
1087 cs_lnum_t end)
1088 {
1089 cs_lnum_t res = 0;
1090
1091 for (cs_lnum_t idx = start; idx < end; ++idx) {
1092 if (cs_lagr_particles_get_flag(p_set, idx, CS_LAGR_PART_TO_DELETE))
1093 res++;
1094 }
1095 return res;
1096 }
1097
1098 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1099
1100 /*============================================================================
1101 * Public function definitions
1102 *============================================================================*/
1103
1104 /*----------------------------------------------------------------------------
1105 * Return pointers to lagrangian arrays
1106 *
1107 * This function is intended for use by Fortran wrappers.
1108 *
1109 * parameters:
1110 * dim_bound_stat --> dimensions for bound_stat pointer
1111 * p_bound_stat --> bound_stat pointer
1112 *----------------------------------------------------------------------------*/
1113
1114 void
cs_lagr_init_c_arrays(int dim_cs_glob_lagr_source_terms[2],cs_real_t ** p_cs_glob_lagr_source_terms)1115 cs_lagr_init_c_arrays(int dim_cs_glob_lagr_source_terms[2],
1116 cs_real_t **p_cs_glob_lagr_source_terms)
1117 {
1118 cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
1119 int n_boundary_stats = cs_glob_lagr_dim->n_boundary_stats;
1120
1121 assert(bound_stat == NULL);
1122
1123 if (n_boundary_stats > 0)
1124 BFT_MALLOC(bound_stat, n_b_faces * n_boundary_stats, cs_real_t);
1125
1126 BFT_MALLOC(cs_glob_lagr_source_terms->st_val,
1127 cs_glob_lagr_dim->ntersl * cs_glob_mesh->n_cells_with_ghosts,
1128 cs_real_t);
1129 for (cs_lnum_t i = 0; i < cs_glob_lagr_dim->ntersl; i++) {
1130 cs_real_t *st = cs_glob_lagr_source_terms->st_val
1131 + i*cs_glob_mesh->n_cells_with_ghosts;
1132 cs_array_set_value_real(cs_glob_mesh->n_cells_with_ghosts,
1133 1,
1134 0,
1135 st);
1136 }
1137
1138 *p_cs_glob_lagr_source_terms = cs_glob_lagr_source_terms->st_val;
1139 dim_cs_glob_lagr_source_terms[0] = cs_glob_mesh->n_cells_with_ghosts;
1140 dim_cs_glob_lagr_source_terms[1] = cs_glob_lagr_dim->ntersl;
1141 }
1142
1143 /*----------------------------------------------------------------------------
1144 * Free lagrangian arrays
1145 *
1146 * This function is intended for use by Fortran wrappers.
1147 *----------------------------------------------------------------------------*/
1148
1149 void
cs_lagr_finalize(void)1150 cs_lagr_finalize(void)
1151 {
1152 int n_boundary_stats = cs_glob_lagr_dim->n_boundary_stats;
1153
1154 if (n_boundary_stats > 0) {
1155 assert(bound_stat != NULL);
1156 BFT_FREE(bound_stat);
1157 }
1158
1159 BFT_FREE(cs_glob_lagr_precipitation_model->nbprec);
1160 BFT_FREE(cs_glob_lagr_precipitation_model->solub);
1161
1162 BFT_FREE(cs_glob_lagr_precipitation_model->mp_diss);
1163
1164 BFT_FREE(cs_glob_lagr_source_terms->st_val);
1165
1166 /* geometry */
1167
1168 BFT_FREE(cs_glob_lagr_b_face_proj);
1169
1170 /* encrustation pointers */
1171
1172 BFT_FREE(cs_glob_lagr_encrustation->enc1);
1173 BFT_FREE(cs_glob_lagr_encrustation->enc2);
1174 BFT_FREE(cs_glob_lagr_encrustation->tprenc);
1175 BFT_FREE(cs_glob_lagr_encrustation->visref);
1176
1177 /* boundary interaction pointers */
1178
1179 for (int i = 0; i < cs_glob_lagr_dim->n_boundary_stats; i++)
1180 BFT_FREE(cs_glob_lagr_boundary_interactions->nombrd[i]);
1181 BFT_FREE(cs_glob_lagr_boundary_interactions->nombrd);
1182
1183 /* Source terms */
1184
1185 BFT_FREE(cs_glob_lagr_source_terms->itsmv1);
1186 BFT_FREE(cs_glob_lagr_source_terms->itsmv2);
1187
1188 /* Statistics */
1189
1190 cs_lagr_stat_finalize();
1191
1192 /* Also close log file (TODO move this) */
1193
1194 cs_lagr_print_finalize();
1195
1196 /* Close tracking structures */
1197
1198 cs_lagr_tracking_finalize();
1199
1200 cs_lagr_finalize_zone_conditions();
1201
1202 /* Fluid gradients */
1203 cs_lagr_extra_module_t *extra = cs_glob_lagr_extra_module;
1204 BFT_FREE(extra->grad_pr);
1205 if (extra->grad_vel != NULL)
1206 BFT_FREE(extra->grad_vel);
1207 }
1208
1209 /*----------------------------------------------------------------------------*/
1210 /*!
1211 * \brief Provide access to injection set structure.
1212 *
1213 * This access method ensures the strucure is initialized for the given
1214 * zone and injection set.
1215 *
1216 * \param[in] zone_data pointer to boundary or volume conditions structure
1217 * \param[in] zone_id zone id
1218 * \param[in] set_id injection set id
1219 *
1220 * \return pointer to injection set data structure
1221 */
1222 /*----------------------------------------------------------------------------*/
1223
1224 cs_lagr_injection_set_t *
cs_lagr_get_injection_set(cs_lagr_zone_data_t * zone_data,int zone_id,int set_id)1225 cs_lagr_get_injection_set(cs_lagr_zone_data_t *zone_data,
1226 int zone_id,
1227 int set_id)
1228 {
1229 assert(zone_data != NULL);
1230 assert(zone_id >= 0 && set_id >= 0);
1231 assert(zone_id < zone_data->n_zones);
1232
1233 if (set_id >= zone_data->n_injection_sets[zone_id])
1234 _zone_injection_set_init(zone_data->location_id,
1235 zone_id,
1236 set_id,
1237 &(zone_data->n_injection_sets[zone_id]),
1238 &(zone_data->injection_set[zone_id]));
1239
1240 return &(zone_data->injection_set[zone_id][set_id]);
1241 }
1242
1243 /*----------------------------------------------------------------------------*/
1244 /*!
1245 * \brief Initialize injection set data structure fields to defaults.
1246 *
1247 * \param[in, out] zis pointer to structure to initialize
1248 */
1249 /*----------------------------------------------------------------------------*/
1250
1251 void
cs_lagr_injection_set_default(cs_lagr_injection_set_t * zis)1252 cs_lagr_injection_set_default(cs_lagr_injection_set_t *zis)
1253 {
1254 assert(zis != NULL);
1255
1256 zis->n_inject = 0;
1257 zis->injection_frequency = 0;
1258
1259 zis->injection_profile_func = NULL;
1260 zis->injection_profile_input = NULL;
1261
1262 /* Fluid velocity by default */
1263 zis->velocity_profile = -1;
1264 /* Fluid temperature by default */
1265 zis->temperature_profile = 0;
1266
1267 if (cs_glob_lagr_model->physical_model == CS_LAGR_PHYS_COAL)
1268 zis->coal_number = -2;
1269
1270 zis->cluster = 0;
1271
1272 /* Agglomeration/fragmentation default*/
1273 zis->aggregat_class_id = 1;
1274 zis->aggregat_fractal_dim = 3;
1275
1276 zis->velocity_magnitude = - cs_math_big_r;
1277
1278 for (int i = 0; i < 3; i++)
1279 zis->velocity[i] = - cs_math_big_r;
1280
1281 /* For spheroids without inertia */
1282 /* Default shape: sphere */
1283 zis->shape = CS_LAGR_SHAPE_SPHERE_MODEL;
1284
1285 /* Angular velocity */
1286 for (int i = 0; i < 3; i++)
1287 zis->angular_vel[i] = 0.;
1288
1289 /* Spheroids radii a b c */
1290 for (int i = 0; i < 3; i++)
1291 zis->radii[i] = 0.;
1292
1293 /* Shape parameters */
1294 for (int i = 0; i < 4; i++)
1295 zis->shape_param[i] = 0.;
1296
1297 /* First three Euler parameters */
1298 for (int i = 0; i < 3; i++)
1299 zis->euler[i] = 0.;
1300
1301 zis->euler[3] = 1.;
1302
1303 zis->stat_weight = - cs_math_big_r;
1304 zis->diameter = - cs_math_big_r;
1305 zis->diameter_variance = - cs_math_big_r;
1306 zis->density = - cs_math_big_r;
1307
1308 zis->temperature = - cs_math_big_r;
1309 zis->cp = - cs_math_big_r;
1310 zis->emissivity = - cs_math_big_r;
1311
1312 zis->flow_rate = 0.0;
1313 }
1314
1315 /*----------------------------------------------------------------------------*/
1316 /*!
1317 * \brief Get read/write pointer to global particle counter
1318 *
1319 * \return
1320 * pointer to lagrangian particle counter structure
1321 */
1322 /*----------------------------------------------------------------------------*/
1323
1324 cs_lagr_particle_counter_t *
cs_lagr_get_particle_counter(void)1325 cs_lagr_get_particle_counter(void)
1326 {
1327 return &_lagr_particle_counter;
1328 }
1329
1330 /*----------------------------------------------------------------------------*/
1331 /*!
1332 \brief Update global particle counter
1333 *
1334 * All fields handled in the local particle set are updated relative
1335 * to that data (using global sums).
1336 *
1337 * \return pointer to lagrangian particle counter structure
1338 */
1339 /*----------------------------------------------------------------------------*/
1340
1341 cs_lagr_particle_counter_t *
cs_lagr_update_particle_counter(void)1342 cs_lagr_update_particle_counter(void)
1343 {
1344 cs_lagr_particle_set_t *p_set = cs_glob_lagr_particle_set;
1345 cs_lagr_particle_counter_t *pc = &_lagr_particle_counter;
1346
1347 cs_gnum_t gcount[] = {p_set->n_particles,
1348 p_set->n_part_new,
1349 p_set->n_part_merged,
1350 p_set->n_part_out,
1351 p_set->n_part_dep,
1352 p_set->n_part_fou,
1353 p_set->n_part_resusp,
1354 p_set->n_failed_part};
1355
1356 cs_real_t wsum[] = {p_set->weight,
1357 p_set->weight_new,
1358 p_set->weight_merged,
1359 p_set->weight_out,
1360 p_set->weight_dep,
1361 p_set->weight_fou,
1362 p_set->weight_resusp};
1363
1364 cs_lnum_t size_count = sizeof(gcount) / sizeof(gcount[0]);
1365 cs_lnum_t size_sum = sizeof(wsum) / sizeof(wsum[0]);
1366
1367 cs_parall_counter(gcount, size_count);
1368 cs_parall_sum(size_sum, CS_REAL_TYPE, wsum);
1369
1370 pc->n_g_total = gcount[0];
1371 pc->n_g_new = gcount[1];
1372 pc->n_g_merged = gcount[2];
1373 pc->n_g_exit = gcount[3];
1374 pc->n_g_deposited = gcount[4];
1375 pc->n_g_fouling = gcount[5];
1376 pc->n_g_resuspended = gcount[6];
1377 pc->n_g_failed = gcount[7];
1378
1379 pc->w_total = wsum[0];
1380 pc->w_new = wsum[1];
1381 pc->w_merged = wsum[2];
1382 pc->w_exit = wsum[3];
1383 pc->w_deposited = wsum[4];
1384 pc->w_fouling = wsum[5];
1385 pc->w_resuspended = wsum[6];
1386
1387 return pc;
1388 }
1389
1390 /*----------------------------------------------------------------------------*/
1391 /*! \brief Provide access to cs_lagr_specific_physics_t
1392 *
1393 * \return
1394 * pointer to lagrangian specific physics options
1395 */
1396 /*----------------------------------------------------------------------------*/
1397
1398 cs_lagr_specific_physics_t *
cs_get_lagr_specific_physics(void)1399 cs_get_lagr_specific_physics(void)
1400 {
1401 return &_cs_glob_lagr_specific_physics;
1402 }
1403
1404 /*----------------------------------------------------------------------------*/
1405 /*! \brief Provide access to cs_lagr_reentrained_model_t
1406 *
1407 * \return
1408 * pointer to lagrangian reentrained model options
1409 */
1410 /*----------------------------------------------------------------------------*/
1411
1412 cs_lagr_reentrained_model_t *
cs_get_lagr_reentrained_model(void)1413 cs_get_lagr_reentrained_model(void)
1414 {
1415 return &_lagr_reentrained_model;
1416 }
1417
1418 /*----------------------------------------------------------------------------*/
1419 /*! \brief Provide access to cs_lagr_precipitation_model_t
1420 *
1421 * \return pointer to lagrangian precipitation model options
1422 */
1423 /*----------------------------------------------------------------------------*/
1424
1425 cs_lagr_precipitation_model_t *
cs_get_lagr_precipitation_model(void)1426 cs_get_lagr_precipitation_model(void)
1427 {
1428 return &_cs_glob_lagr_precipitation_model;
1429 }
1430
1431 /*----------------------------------------------------------------------------
1432 * Provide access to cs_lagr_clogging_model_t
1433 *
1434 * needed to initialize structure with GUI
1435 *----------------------------------------------------------------------------*/
1436
1437 cs_lagr_clogging_model_t *
cs_get_lagr_clogging_model(void)1438 cs_get_lagr_clogging_model(void)
1439 {
1440 return &_cs_glob_lagr_clogging_model;
1441 }
1442
1443 /*----------------------------------------------------------------------------
1444 * Provide access to cs_lagr_shape_model_t
1445 *
1446 * needed to initialize structure with GUI
1447 *----------------------------------------------------------------------------*/
1448
1449 cs_lagr_shape_model_t *
cs_get_lagr_shape_model(void)1450 cs_get_lagr_shape_model(void)
1451 {
1452 return &_cs_glob_lagr_shape_model;
1453 }
1454
1455 /*----------------------------------------------------------------------------
1456 * Provide access to cs_lagr_agglomeration_model_t
1457 *
1458 * needed to initialize structure with GUI
1459 *----------------------------------------------------------------------------*/
1460
1461 cs_lagr_agglomeration_model_t *
cs_get_lagr_agglomeration_model(void)1462 cs_get_lagr_agglomeration_model(void)
1463 {
1464 return &_cs_glob_lagr_agglomeration_model;
1465 }
1466
1467 /*----------------------------------------------------------------------------
1468 * Provide access to cs_lagr_consolidation_model_t
1469 *
1470 * needed to initialize structure with GUI
1471 *----------------------------------------------------------------------------*/
1472
1473 cs_lagr_consolidation_model_t *
cs_get_lagr_consolidation_model(void)1474 cs_get_lagr_consolidation_model(void)
1475 {
1476 return &_cs_glob_lagr_consolidation_model;
1477 }
1478
1479 /*----------------------------------------------------------------------------
1480 * Provide access to cs_lagr_time_step_t
1481 *
1482 * needed to initialize structure with GUI
1483 *----------------------------------------------------------------------------*/
1484
1485 cs_lagr_time_step_t *
cs_get_lagr_time_step(void)1486 cs_get_lagr_time_step(void)
1487 {
1488 return &_cs_glob_lagr_time_step;
1489 }
1490
1491 /*----------------------------------------------------------------------------
1492 * Provide access to cs_lagr_source_terms_t
1493 *
1494 * needed to initialize structure with GUI
1495 *----------------------------------------------------------------------------*/
1496
1497 cs_lagr_source_terms_t *
cs_get_lagr_source_terms(void)1498 cs_get_lagr_source_terms(void)
1499 {
1500 return &_cs_glob_lagr_source_terms;
1501 }
1502
1503 /*----------------------------------------------------------------------------
1504 * Provide access to cs_lagr_encrustation_t
1505 *
1506 * needed to initialize structure with GUI
1507 *----------------------------------------------------------------------------*/
1508
1509 cs_lagr_encrustation_t *
cs_get_lagr_encrustation(void)1510 cs_get_lagr_encrustation(void)
1511 {
1512 return &_cs_glob_lagr_encrustation;
1513 }
1514
1515 /*----------------------------------------------------------------------------
1516 * Provide access to cs_lagr_physico_chemical_t
1517 *
1518 * needed to initialize structure with GUI
1519 *----------------------------------------------------------------------------*/
1520
1521 cs_lagr_physico_chemical_t *
cs_get_lagr_physico_chemical(void)1522 cs_get_lagr_physico_chemical(void)
1523 {
1524 return &_cs_glob_lagr_physico_chemical;
1525 }
1526
1527 /*----------------------------------------------------------------------------
1528 * Provide access to cs_lagr_brownian_t
1529 *
1530 * needed to initialize structure with GUI
1531 *----------------------------------------------------------------------------*/
1532
1533 cs_lagr_brownian_t *
cs_get_lagr_brownian(void)1534 cs_get_lagr_brownian(void)
1535 {
1536 return &_cs_glob_lagr_brownian;
1537 }
1538
1539 /*----------------------------------------------------------------------------*/
1540 /*!
1541 * \brief Return pointer to the main internal conditions structure.
1542 *
1543 * \return
1544 * pointer to current internal_contditions or NULL
1545 */
1546 /*----------------------------------------------------------------------------*/
1547
1548 cs_lagr_internal_condition_t *
cs_lagr_get_internal_conditions(void)1549 cs_lagr_get_internal_conditions(void)
1550 {
1551 /* Define a structure with default parameters if not done yet */
1552
1553 if (cs_glob_lagr_internal_conditions == NULL)
1554 cs_glob_lagr_internal_conditions = _create_internal_cond_struct();
1555
1556 if (cs_glob_lagr_internal_conditions->i_face_zone_id == NULL) {
1557 BFT_MALLOC(cs_glob_lagr_internal_conditions->i_face_zone_id,
1558 cs_glob_mesh->n_i_faces,
1559 int);
1560
1561 for (cs_lnum_t i = 0; i < cs_glob_mesh->n_i_faces; i++)
1562 cs_glob_lagr_internal_conditions->i_face_zone_id[i] = -1;
1563 }
1564
1565 return cs_glob_lagr_internal_conditions;
1566 }
1567
1568 /*----------------------------------------------------------------------------*/
1569 /*!
1570 * \brief Return pointer to the main boundary conditions structure.
1571 *
1572 * \return pointer to current boundary zone data structure
1573 */
1574 /*----------------------------------------------------------------------------*/
1575
1576 cs_lagr_zone_data_t *
cs_lagr_get_boundary_conditions(void)1577 cs_lagr_get_boundary_conditions(void)
1578 {
1579 /* Define a structure with default parameters if not done yet */
1580
1581 cs_lnum_t n_zones = cs_boundary_zone_n_zones();
1582
1583 _update_zone_data_struct(&_boundary_conditions,
1584 CS_MESH_LOCATION_BOUNDARY_FACES,
1585 n_zones);
1586
1587 cs_glob_lagr_boundary_conditions = _boundary_conditions;
1588
1589 return _boundary_conditions;
1590 }
1591
1592 /*----------------------------------------------------------------------------*/
1593 /*!
1594 * \brief Return pointer to the main volume conditions structure.
1595 *
1596 * \return pointer to current volume zone data structure
1597 */
1598 /*----------------------------------------------------------------------------*/
1599
1600 cs_lagr_zone_data_t *
cs_lagr_get_volume_conditions(void)1601 cs_lagr_get_volume_conditions(void)
1602 {
1603 /* Define a structure with default parameters if not done yet */
1604
1605 cs_lnum_t n_zones = cs_volume_zone_n_zones();
1606
1607 _update_zone_data_struct(&_volume_conditions,
1608 CS_MESH_LOCATION_CELLS,
1609 n_zones);
1610
1611 cs_glob_lagr_volume_conditions = _volume_conditions;
1612
1613 return _volume_conditions;
1614 }
1615
1616 /*----------------------------------------------------------------------------*/
1617 /*!
1618 * \brief Return pointer to the main internal conditions structure.
1619 *
1620 * \return pointer to current internal conditions structure
1621 */
1622 /*----------------------------------------------------------------------------*/
1623
1624 void
cs_lagr_finalize_zone_conditions(void)1625 cs_lagr_finalize_zone_conditions(void)
1626 {
1627 cs_lagr_zone_data_t *zda[2] = {_boundary_conditions,
1628 _volume_conditions};
1629
1630 for (int i = 0; i < 2; i++) {
1631
1632 cs_lagr_zone_data_t *zd = zda[i];
1633
1634 if (zd != NULL) {
1635
1636 BFT_FREE(zd->zone_type);
1637 for (int j = 0; j < zd->n_zones; j++)
1638 BFT_FREE(zd->injection_set[j]);
1639 BFT_FREE(zd->injection_set);
1640 BFT_FREE(zd->n_injection_sets);
1641
1642 BFT_FREE(zd->elt_type);
1643 BFT_FREE(zd->particle_flow_rate);
1644
1645 BFT_FREE(zda[i]);
1646
1647 }
1648
1649 }
1650 }
1651
1652 /*----------------------------------------------------------------------------
1653 * Destroy finalize the global cs_lagr_internal_condition_t structure.
1654 *----------------------------------------------------------------------------*/
1655
cs_lagr_finalize_internal_cond(void)1656 void cs_lagr_finalize_internal_cond(void)
1657 {
1658 cs_lagr_internal_condition_t *internal_cond
1659 = cs_glob_lagr_internal_conditions;
1660 if (internal_cond != NULL) {
1661 BFT_FREE(internal_cond->i_face_zone_id);
1662 BFT_FREE(internal_cond);
1663 }
1664 }
1665
1666 /*----------------------------------------------------------------------------
1667 * Provide access to cs_lagr_boundary_interactions_t
1668 *
1669 * needed to initialize structure with GUI
1670 *----------------------------------------------------------------------------*/
1671
1672 cs_lagr_boundary_interactions_t *
cs_get_lagr_boundary_interactions(void)1673 cs_get_lagr_boundary_interactions(void)
1674 {
1675 return &_cs_glob_lagr_boundary_interactions;
1676 }
1677
1678 /*----------------------------------------------------------------------------
1679 * Provide access to cs_lagr_extra_module_t
1680 *
1681 *----------------------------------------------------------------------------*/
1682
1683 cs_lagr_extra_module_t *
cs_get_lagr_extra_module(void)1684 cs_get_lagr_extra_module(void)
1685 {
1686 return &_lagr_extra_module;
1687 }
1688
1689 /*----------------------------------------------------------------------------
1690 * Prepare for execution of the Lagrangian model.
1691 *
1692 * This should be called before the fist call to cs_lagr_solve_time_step.
1693 *
1694 * parameters:
1695 * dt <-- time step (per cell)
1696 *----------------------------------------------------------------------------*/
1697
1698 void
cs_lagr_solve_initialize(const cs_real_t * dt)1699 cs_lagr_solve_initialize(const cs_real_t *dt)
1700 {
1701 CS_UNUSED(dt);
1702
1703 /* Allocate pressure and velocity gradients */
1704 cs_lagr_extra_module_t *extra = cs_glob_lagr_extra_module;
1705 cs_lnum_t ncelet = cs_glob_mesh->n_cells_with_ghosts;
1706
1707 BFT_MALLOC(extra->grad_pr, ncelet, cs_real_3_t);
1708 if ( cs_glob_lagr_model->modcpl > 0
1709 || cs_glob_lagr_model->shape > 0)
1710 BFT_MALLOC(extra->grad_vel, ncelet, cs_real_33_t);
1711
1712 /* For frozen field:
1713 values at previous time step = values at current time step */
1714
1715 if (cs_glob_lagr_time_scheme->iilagr == CS_LAGR_FROZEN_CONTINUOUS_PHASE) {
1716
1717 int n_fields = cs_field_n_fields();
1718
1719 for (int f_id = 0; f_id < n_fields; f_id++) {
1720
1721 cs_field_t *f = cs_field_by_id(f_id);
1722 if (f->type & CS_FIELD_VARIABLE)
1723 cs_field_current_to_previous(f);
1724
1725 }
1726
1727 }
1728
1729 /* First call (initializations)
1730 ---------------------------- */
1731
1732 _lagr_map_fields_default();
1733
1734 cs_lagr_tracking_initialize();
1735
1736 /* first initializations */
1737
1738 cs_lagr_post_init();
1739
1740 /* Read particle restart data */
1741
1742 if (cs_glob_lagr_time_scheme->iilagr != CS_LAGR_OFF)
1743 cs_lagr_restart_read_p();
1744
1745 /* Compute gradients of current value fields */
1746 if (cs_glob_lagr_time_scheme->iilagr == CS_LAGR_FROZEN_CONTINUOUS_PHASE)
1747 cs_lagr_gradients(0, extra->grad_pr, extra->grad_vel);
1748
1749 /* Read statistics restart data */
1750
1751 cs_lagr_stat_restart_read();
1752 }
1753
1754 /*----------------------------------------------------------------------------
1755 * Execute one time step of the Lagrangian model.
1756 *
1757 * This is the main function for that model.
1758 *
1759 * parameters:
1760 * itypfb <-- boundary face types
1761 * dt <-- time step (per cell)
1762 *----------------------------------------------------------------------------*/
1763
1764 void
cs_lagr_solve_time_step(const int itypfb[],const cs_real_t * dt)1765 cs_lagr_solve_time_step(const int itypfb[],
1766 const cs_real_t *dt)
1767 {
1768 static int ipass = 0;
1769 cs_time_step_t *ts = cs_get_glob_time_step();
1770
1771 int mode;
1772
1773 cs_lagr_boundary_interactions_t *lag_bdi = cs_glob_lagr_boundary_interactions;
1774 cs_lagr_model_t *lagr_model = cs_glob_lagr_model;
1775 cs_lagr_extra_module_t *extra = cs_glob_lagr_extra_module;
1776 cs_lagr_particle_counter_t *part_c = cs_lagr_get_particle_counter();
1777
1778 cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
1779
1780 cs_lnum_t *ifabor = cs_glob_mesh->b_face_cells;
1781 cs_real_t *surfbo = cs_glob_mesh_quantities->b_face_surf;
1782 cs_real_t *surfbn = cs_glob_mesh_quantities->b_face_normal;
1783
1784 /* Allocate arrays depending on user options */
1785
1786 cs_real_t *tempp = NULL;
1787 if ( lagr_model->clogging == 1
1788 || lagr_model->roughness == 1
1789 || lagr_model->dlvo == 1)
1790 BFT_MALLOC(tempp, cs_glob_mesh->n_cells, cs_real_t);
1791
1792 ipass ++;
1793
1794 static cs_real_t *vislen = NULL;
1795 if ((lagr_model->deposition == 1) && (ipass == 1)) {
1796 BFT_MALLOC(vislen, n_b_faces, cs_real_t);
1797 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++)
1798 vislen[ifac] = -cs_math_big_r;
1799 }
1800
1801 /* mask for for deposition-related states */
1802
1803 int deposition_mask = CS_LAGR_PART_DEPOSITED | CS_LAGR_PART_ROLLING
1804 | CS_LAGR_PART_IMPOSED_MOTION;
1805
1806 /* Initialization
1807 -------------- */
1808
1809 cs_lagr_particle_set_t *p_set = cs_glob_lagr_particle_set;
1810
1811 /* First call (initializations)
1812 ---------------------------- */
1813
1814 if (cs_glob_time_step->nt_cur == cs_glob_time_step->nt_prev + 1) {
1815
1816 /* If the deposition model is activated */
1817
1818 if (lagr_model->deposition > 0) {
1819
1820 cs_real_t ustarmoy = 0.0;
1821 cs_real_t surftot = 0.0;
1822 cs_real_3_t dtmp;
1823 dtmp[0] = 0.0;
1824
1825 /* Boundary faces data */
1826
1827 cs_lagr_geom();
1828
1829 /* Average friction velocity calculation */
1830
1831 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++) {
1832
1833 if (itypfb[ifac] == CS_SMOOTHWALL || itypfb[ifac] == CS_ROUGHWALL) {
1834
1835 cs_lnum_t iel = ifabor[ifac];
1836
1837 /* the density pointer according to the flow location */
1838 cs_real_t romf = extra->cromf->val[iel];
1839 cs_real_t visccf = extra->viscl->val[iel] / romf;
1840
1841 cs_real_t _ustar = CS_MAX(extra->ustar->val[ifac], 1e-15);
1842
1843 cs_real_t surfb = surfbo[ifac];
1844 ustarmoy = ustarmoy + surfb * _ustar;
1845 surftot = surftot + surfb;
1846 vislen[ifac] = visccf / _ustar; // FIXME to be coherent with wall fn: y/y+
1847 dtmp[0] = dtmp[0] + 1.0;
1848
1849 }
1850
1851 }
1852
1853 if (cs_glob_rank_id >= 0) {
1854 dtmp[1] = ustarmoy;
1855 dtmp[2] = surftot;
1856
1857 cs_parall_sum(3, CS_DOUBLE, dtmp);
1858
1859 ustarmoy = dtmp[1];
1860 surftot = dtmp[2];
1861 }
1862
1863 if (dtmp[0] > 0.0)
1864 ustarmoy = ustarmoy / surftot;
1865
1866 /* Average friction velocity display */
1867
1868 if (cs_glob_rank_id <= 0)
1869 bft_printf(_(" ** LAGRANGIAN MODULE:\n"
1870 " ** deposition submodel\n"
1871 "---------------------------------------------\n"
1872 "** Mean friction velocity (ustar) = %7.3f\n"
1873 "---------------------------------------------\n"),
1874 ustarmoy);
1875
1876 }
1877
1878 }
1879
1880 /* Prepare statistics for this time step */
1881
1882 cs_lagr_stat_prepare();
1883
1884 /* Update boundary condition types;
1885 in most cases, this should be useful only at the first iteration,
1886 but we prefer to be safe in case of advanced uses with time-dependent
1887 zones or zone types */
1888
1889 /* User initialization by set and boundary */
1890
1891 if (ts->nt_cur == ts->nt_prev +1)
1892 cs_gui_particles_bcs();
1893
1894 /* User initialization by set and zone */
1895
1896 cs_user_lagr_boundary_conditions(itypfb);
1897 cs_user_lagr_volume_conditions();
1898
1899 _update_boundary_face_type();
1900
1901 /* Initialize counter */
1902
1903 part_c->n_g_total = 0;
1904 part_c->n_g_new = 0;
1905 part_c->n_g_exit = 0;
1906 part_c->n_g_merged = 0;
1907 part_c->n_g_deposited = 0;
1908 part_c->n_g_fouling = 0;
1909 part_c->n_g_resuspended = 0;
1910 part_c->n_g_failed = 0;
1911 part_c->w_total = 0;
1912 part_c->w_new = 0;
1913 part_c->w_exit = 0;
1914 part_c->w_merged = 0;
1915 part_c->w_deposited = 0;
1916 part_c->w_fouling = 0;
1917 part_c->w_resuspended = 0;
1918
1919 /* particles->n_part_new: handled in injection step */
1920 p_set->weight = 0.0;
1921 p_set->n_part_out = 0;
1922 p_set->n_part_dep = 0;
1923 p_set->n_part_fou = 0;
1924 p_set->weight_out = 0.0;
1925 p_set->weight_dep = 0.0;
1926 p_set->weight_fou = 0.0;
1927 p_set->n_failed_part = 0;
1928 p_set->weight_failed = 0.0;
1929
1930 /* Initialization for the dlvo, roughness and clogging model
1931 ---------------------------------------------------------- */
1932
1933 if ( lagr_model->dlvo == 1
1934 || lagr_model->roughness == 1
1935 || lagr_model->clogging == 1) {
1936
1937 if (extra->temperature != NULL) {
1938
1939 if (cs_glob_thermal_model->itpscl == CS_TEMPERATURE_SCALE_CELSIUS) {
1940 for (cs_lnum_t iel = 0; iel < cs_glob_mesh->n_cells; iel++) {
1941 tempp[iel] = extra->temperature->val[iel]
1942 + cs_physical_constants_celsius_to_kelvin;
1943 }
1944 }
1945 else {
1946 for (cs_lnum_t iel = 0; iel < cs_glob_mesh->n_cells; iel++)
1947 tempp[iel] = extra->temperature->val[iel];
1948 }
1949
1950 }
1951
1952 else {
1953 if (cs_glob_thermal_model->itpscl == CS_TEMPERATURE_SCALE_CELSIUS) {
1954 for (cs_lnum_t iel = 0; iel < cs_glob_mesh->n_cells; iel++) {
1955 tempp[iel] = cs_glob_fluid_properties->t0
1956 + cs_physical_constants_celsius_to_kelvin;
1957 }
1958 }
1959 else {
1960 for (cs_lnum_t iel = 0; iel < cs_glob_mesh->n_cells; iel++)
1961 tempp[iel] = cs_glob_fluid_properties->t0;
1962 }
1963 }
1964
1965 }
1966
1967 /* Initialization for the dlvo model
1968 --------------------------------- */
1969
1970 if (lagr_model->dlvo == 1)
1971 cs_lagr_dlvo_init(cs_glob_lagr_physico_chemical->epseau,
1972 cs_glob_lagr_physico_chemical->fion,
1973 tempp,
1974 cs_glob_lagr_physico_chemical->valen,
1975 cs_glob_lagr_physico_chemical->phi_p,
1976 cs_glob_lagr_physico_chemical->phi_s,
1977 cs_glob_lagr_physico_chemical->cstham,
1978 cs_glob_lagr_clogging_model->csthpp,
1979 cs_glob_lagr_physico_chemical->lambda_vdw);
1980
1981 /* Initialization for the roughness surface model
1982 ---------------------------------------------- */
1983
1984 if (lagr_model->roughness == 1)
1985 roughness_init(&cs_glob_lagr_physico_chemical->epseau,
1986 &cs_glob_lagr_physico_chemical->fion,
1987 tempp,
1988 &cs_glob_lagr_physico_chemical->valen,
1989 &cs_glob_lagr_physico_chemical->phi_p,
1990 &cs_glob_lagr_physico_chemical->phi_s,
1991 &cs_glob_lagr_physico_chemical->cstham,
1992 &cs_glob_lagr_physico_chemical->lambda_vdw,
1993 &cs_glob_lagr_reentrained_model->espasg,
1994 &cs_glob_lagr_reentrained_model->denasp,
1995 &cs_glob_lagr_reentrained_model->rayasp,
1996 &cs_glob_lagr_reentrained_model->rayasg);
1997
1998 /* Initialization for the clogging model
1999 ------------------------------------- */
2000
2001 if (lagr_model->clogging == 1)
2002 cloginit(&cs_glob_lagr_physico_chemical->epseau,
2003 &cs_glob_lagr_physico_chemical->fion,
2004 &cs_glob_lagr_clogging_model->jamlim,
2005 &cs_glob_lagr_clogging_model->mporos,
2006 &cs_glob_lagr_clogging_model->diam_mean,
2007 tempp,
2008 &cs_glob_lagr_physico_chemical->valen,
2009 &cs_glob_lagr_physico_chemical->phi_p,
2010 &cs_glob_lagr_physico_chemical->phi_s,
2011 &cs_glob_lagr_physico_chemical->cstham,
2012 &cs_glob_lagr_clogging_model->csthpp,
2013 &cs_glob_lagr_physico_chemical->lambda_vdw);
2014
2015 /* Initialization for the nonsphere model
2016 ------------------------------------- */
2017
2018 if (lagr_model->shape == CS_LAGR_SHAPE_SPHEROID_STOC_MODEL)
2019 cs_glob_lagr_shape_model->param_chmb = 1.0;
2020
2021 /* Update for new particles which entered the domain
2022 ------------------------------------------------- */
2023
2024 /* At the first time step we initialize particles to
2025 values at current time step and not at previous time step, because
2026 values at previous time step = initialization */
2027
2028 int iprev;
2029 if (ts->nt_cur == 1) /* Use fields at current time step */
2030 iprev = 0;
2031
2032 else /* Use fields at previous time step */
2033 iprev = 1;
2034
2035 /* Fields at current time step if using NEPTUNE_CFD */
2036 if (cs_field_by_name_try("velocity_1") != NULL)
2037 iprev = 0;
2038
2039 cs_lagr_injection(iprev, itypfb, vislen);
2040
2041 /* Initialization for the agglomeration/fragmentation models
2042 --------------------------------------------------------- */
2043
2044 /* The algorithms are based on discrete values of particle radii
2045 (CS_LAGR_DIAMETER).
2046 It uses a minimum diameter, corresponding to monomers
2047 (unbreakable particles).
2048 Aggregates are formed of multiples of these monomers. */
2049
2050 /* Evaluation of the minimum diameter */
2051 cs_real_t minimum_particle_diam = 0.;
2052
2053 if (cs_glob_lagr_model->agglomeration == 1)
2054 minimum_particle_diam = cs_glob_lagr_agglomeration_model->base_diameter;
2055
2056 if (cs_glob_lagr_model->fragmentation == 1)
2057 minimum_particle_diam = cs_glob_lagr_fragmentation_model->base_diameter;
2058
2059 /* Management of advancing time
2060 ---------------------------- */
2061
2062 cs_glob_lagr_time_step->dtp = dt[0];
2063 cs_glob_lagr_time_step->ttclag += cs_glob_lagr_time_step->dtp;
2064
2065 part_c = cs_lagr_update_particle_counter();
2066
2067 /* Global number of particles > 0 */
2068 if (part_c->n_g_total > 0) {
2069
2070 /* Record particle's starting cell and rank, and update rebound time id */
2071
2072 for (cs_lnum_t ip = 0; ip < p_set->n_particles; ip++) {
2073
2074 cs_lagr_particles_set_lnum_n
2075 (p_set, ip, 1, CS_LAGR_CELL_ID,
2076 cs_lagr_particles_get_lnum(p_set, ip, CS_LAGR_CELL_ID));
2077
2078 cs_lagr_particles_set_lnum_n(p_set, ip, 1, CS_LAGR_RANK_ID,
2079 cs_glob_rank_id);
2080
2081 cs_lnum_t rebound_id
2082 = cs_lagr_particles_get_lnum(p_set, ip, CS_LAGR_REBOUND_ID);
2083 if (rebound_id >= 0)
2084 cs_lagr_particles_set_lnum(p_set, ip, CS_LAGR_REBOUND_ID,
2085 rebound_id + 1);
2086
2087 }
2088
2089 /* Pressure and fluid velocity gradients
2090 ------------------------------------- */
2091
2092 /* At the first time step we initialize particles to
2093 values at current time step and not at previous time step, because
2094 values at previous time step = initialization (zero gradients) */
2095
2096 if (cs_glob_lagr_time_scheme->iilagr != CS_LAGR_FROZEN_CONTINUOUS_PHASE) {
2097 mode = 1;
2098 if (ts->nt_cur == 1)
2099 mode = 0;
2100
2101 cs_lagr_gradients(mode, extra->grad_pr, extra->grad_vel);
2102 }
2103
2104 /* Particles progression
2105 --------------------- */
2106
2107 bool go_on = true;
2108 cs_lnum_t n_particles_prev = p_set->n_particles - p_set->n_part_new;
2109 while (go_on) {
2110
2111 cs_glob_lagr_time_step->nor
2112 = cs_glob_lagr_time_step->nor % cs_glob_lagr_time_scheme->t_order;
2113 cs_glob_lagr_time_step->nor++;
2114
2115 /* Allocations */
2116
2117 cs_lnum_t nresnew = 0;
2118
2119 cs_real_t *taup;
2120 cs_real_33_t *bx;
2121 cs_real_3_t *tlag, *piil;
2122 BFT_MALLOC(taup, p_set->n_particles, cs_real_t);
2123 BFT_MALLOC(tlag, p_set->n_particles, cs_real_3_t);
2124 BFT_MALLOC(piil, p_set->n_particles, cs_real_3_t);
2125 BFT_MALLOC(bx, p_set->n_particles, cs_real_33_t);
2126
2127 cs_real_t *tsfext = NULL;
2128 if (cs_glob_lagr_time_scheme->iilagr == CS_LAGR_TWOWAY_COUPLING)
2129 BFT_MALLOC(tsfext, p_set->n_particles, cs_real_t);
2130
2131 cs_real_t *cpgd1 = NULL, *cpgd2 = NULL, *cpght = NULL;
2132 if ( cs_glob_lagr_time_scheme->iilagr == CS_LAGR_TWOWAY_COUPLING
2133 && lagr_model->physical_model == CS_LAGR_PHYS_COAL
2134 && cs_glob_lagr_source_terms->ltsthe == 1) {
2135
2136 BFT_MALLOC(cpgd1, p_set->n_particles, cs_real_t);
2137 BFT_MALLOC(cpgd2, p_set->n_particles, cs_real_t);
2138 BFT_MALLOC(cpght, p_set->n_particles, cs_real_t);
2139
2140 }
2141
2142 cs_real_t *tempct = NULL;
2143 if ( ( lagr_model->physical_model == CS_LAGR_PHYS_HEAT
2144 && cs_glob_lagr_specific_physics->itpvar == 1)
2145 || lagr_model->physical_model == CS_LAGR_PHYS_COAL)
2146 BFT_MALLOC(tempct, p_set->n_particles * 2, cs_real_t);
2147
2148 cs_real_t *terbru = NULL;
2149 if (cs_glob_lagr_brownian->lamvbr == 1)
2150 BFT_MALLOC(terbru, p_set->n_particles, cs_real_t);
2151
2152 /* Copy results from previous step */
2153
2154 if (cs_glob_lagr_time_step->nor == 1) {
2155 /* Current to previous but not on new particles at the first time
2156 * because the user may have changed their position */
2157 for (cs_lnum_t ip = 0; ip < n_particles_prev; ip++)
2158 cs_lagr_particles_current_to_previous(p_set, ip);
2159
2160 n_particles_prev = p_set->n_particles;
2161 }
2162
2163 /* Computation of the fluid's pressure and velocity gradient
2164 at n+1 (with values at current time step) */
2165 if ( cs_glob_lagr_time_step->nor == 2
2166 && cs_glob_lagr_time_scheme->iilagr != CS_LAGR_FROZEN_CONTINUOUS_PHASE)
2167 cs_lagr_gradients(0, extra->grad_pr, extra->grad_vel);
2168
2169 /* use fields at previous or current time step */
2170 if (cs_glob_lagr_time_step->nor == 1)
2171 /* Use fields at previous time step */
2172 iprev = 1;
2173
2174 else
2175 iprev = 0;
2176
2177 /* Retrieve bx values associated with particles from previous pass */
2178
2179 if ( cs_glob_lagr_time_scheme->t_order == 2
2180 && cs_glob_lagr_time_step->nor == 2) {
2181
2182 for (cs_lnum_t ip = 0; ip < p_set->n_particles; ip++) {
2183
2184 cs_real_t *jbx1 = cs_lagr_particles_attr(p_set, ip, CS_LAGR_TURB_STATE_1);
2185
2186 for (cs_lnum_t ii = 0; ii < 3; ii++) {
2187
2188 bx[ip][ii][0] = jbx1[ii];
2189
2190 }
2191
2192 }
2193
2194 }
2195
2196 cs_lagr_car(iprev,
2197 dt,
2198 taup,
2199 tlag,
2200 piil,
2201 bx,
2202 tempct,
2203 extra->grad_pr,
2204 extra->grad_vel);
2205
2206 /* Integration of SDEs: position, fluid and particle velocity */
2207
2208 cs_lagr_sde(cs_glob_lagr_time_step->dtp,
2209 (const cs_real_t *)taup,
2210 (const cs_real_3_t *)tlag,
2211 (const cs_real_3_t *)piil,
2212 (const cs_real_33_t *)bx,
2213 tsfext,
2214 (const cs_real_3_t *)extra->grad_pr,
2215 (const cs_real_33_t *)extra->grad_vel,
2216 terbru,
2217 (const cs_real_t *)vislen,
2218 &nresnew);
2219
2220 /* Integration of SDEs for orientation of spheroids without inertia */
2221 if (lagr_model->shape == CS_LAGR_SHAPE_SPHEROID_STOC_MODEL) {
2222 cs_lagr_orientation_dyn_spheroids(iprev,
2223 cs_glob_lagr_time_step->dtp,
2224 (const cs_real_33_t *)extra->grad_vel);
2225 }
2226 /* Integration of Jeffrey equations for ellispoids */
2227 else if (lagr_model->shape == CS_LAGR_SHAPE_SPHEROID_JEFFERY_MODEL) {
2228 cs_lagr_orientation_dyn_jeffery(cs_glob_lagr_time_step->dtp,
2229 (const cs_real_33_t *)extra->grad_vel);
2230 }
2231
2232
2233 /* Save bx values associated with particles for next pass */
2234
2235 if ( cs_glob_lagr_time_scheme->t_order == 2
2236 && cs_glob_lagr_time_step->nor == 1) {
2237
2238 for (cs_lnum_t ip = 0; ip < p_set->n_particles; ip++) {
2239
2240 cs_real_t *jbx1 = cs_lagr_particles_attr(p_set, ip, CS_LAGR_TURB_STATE_1);
2241
2242 for (int ii = 0; ii < 3; ii++)
2243 jbx1[ii] = bx[ip][ii][0];
2244
2245 }
2246
2247 }
2248
2249 /* Integration of SDE related to physical models */
2250
2251 if (lagr_model->physical_model == CS_LAGR_PHYS_HEAT
2252 || lagr_model->physical_model == CS_LAGR_PHYS_COAL) {
2253
2254 if (cs_glob_lagr_time_step->nor == 1)
2255 /* Use fields at previous time step */
2256 iprev = 1;
2257
2258 else
2259 /* Use fields at current time step */
2260 iprev = 0;
2261
2262 cs_lagr_sde_model(tempct, cpgd1, cpgd2, cpght);
2263
2264 }
2265
2266 /* Integration of additional user variables
2267 -----------------------------------------*/
2268
2269 if (cs_glob_lagr_model->n_user_variables > 0)
2270 cs_user_lagr_sde(dt, taup, tlag, tempct);
2271
2272 /* Integration of agglomeration and fragmentation
2273 ----------------------------------------------*/
2274
2275 /* Agglomeration and fragmentation preparation */
2276
2277 /* Preparation: find cells occupied by particles (number)
2278 generate lists of these cells
2279 generate list particles indexes (sublists within a cell) */
2280
2281 cs_lnum_t n_occupied_cells;
2282
2283 cs_lnum_t *occupied_cell_ids = NULL;
2284 cs_lnum_t *particle_list = NULL;
2285
2286 if ( cs_glob_lagr_model->agglomeration == 1
2287 || cs_glob_lagr_model->fragmentation == 1 ) {
2288
2289 n_occupied_cells
2290 = _get_n_occupied_cells(p_set, 0, p_set->n_particles);
2291
2292 BFT_MALLOC(occupied_cell_ids, n_occupied_cells, cs_lnum_t);
2293 BFT_MALLOC(particle_list, n_occupied_cells+1, cs_lnum_t);
2294
2295 _occupied_cells(p_set, 0, p_set->n_particles,
2296 n_occupied_cells,
2297 occupied_cell_ids,
2298 particle_list);
2299
2300 }
2301
2302 /* Compute agglomeration and fragmentation
2303 (avoid second pass if second order scheme is used) */
2304 if ( cs_glob_lagr_time_step->nor == 1
2305 && ((cs_glob_lagr_model->agglomeration == 1) ||
2306 (cs_glob_lagr_model->fragmentation == 1))) {
2307
2308 /* Initialize lists (ids of cells and particles) */
2309 cs_lnum_t *cell_particle_idx;
2310
2311 BFT_MALLOC(cell_particle_idx, n_occupied_cells+1, cs_lnum_t);
2312 cell_particle_idx[0] = 0;
2313
2314 cs_lnum_t enter_parts = p_set->n_particles;
2315
2316 /* Loop on all cells that contain at least one particle */
2317 for (cs_lnum_t icell = 0; icell < n_occupied_cells; ++icell) {
2318
2319 cs_lnum_t cell_id = occupied_cell_ids[icell];
2320
2321 /* Particle indices: between start_part and end_part (list) */
2322 cs_lnum_t start_part = particle_list[icell];
2323 cs_lnum_t end_part = particle_list[icell+1];
2324
2325 cs_lnum_t init_particles = p_set->n_particles;
2326
2327 /* Treat agglomeration */
2328 if (cs_glob_lagr_model->agglomeration == 1) {
2329
2330 cs_lagr_agglomeration(cell_id,
2331 dt[0],
2332 minimum_particle_diam,
2333 start_part,
2334 end_part);
2335 }
2336
2337 /* Save number of created particles */
2338
2339 cs_lnum_t inserted_parts_agglo = p_set->n_particles - init_particles;
2340
2341 /* Create local buffer (deleted particles at the end) */
2342
2343 cs_lnum_t local_size = end_part - start_part;
2344 cs_lnum_t deleted_parts = _get_n_deleted(p_set, start_part, end_part);
2345 size_t swap_buffer_size = p_set->p_am->extents
2346 * (local_size - deleted_parts);
2347 size_t swap_buffer_deleted = p_set->p_am->extents * deleted_parts;
2348
2349 /* Create buffers for deleted particles */
2350 unsigned char * swap_buffer, *deleted_buffer;
2351 BFT_MALLOC(swap_buffer, swap_buffer_size, unsigned char);
2352 BFT_MALLOC(deleted_buffer, swap_buffer_deleted, unsigned char);
2353
2354 /* Update buffer for existing particles */
2355 cs_lnum_t count_del = 0, count_swap = 0;
2356 for (cs_lnum_t i = start_part; i < end_part; ++i) {
2357 if (cs_lagr_particles_get_flag(p_set, i,
2358 CS_LAGR_PART_TO_DELETE)) {
2359 memcpy(deleted_buffer + p_set->p_am->extents * count_del,
2360 p_set->p_buffer + p_set->p_am->extents * i,
2361 p_set->p_am->extents);
2362 count_del++;
2363 }
2364 else {
2365 memcpy(swap_buffer + p_set->p_am->extents * count_swap,
2366 p_set->p_buffer + p_set->p_am->extents * i,
2367 p_set->p_am->extents);
2368 count_swap++;
2369 }
2370 }
2371
2372 memcpy(p_set->p_buffer + p_set->p_am->extents * start_part,
2373 swap_buffer, swap_buffer_size);
2374 memcpy( p_set->p_buffer
2375 + p_set->p_am->extents * (local_size-deleted_parts+start_part),
2376 deleted_buffer, swap_buffer_deleted);
2377
2378 BFT_FREE(deleted_buffer);
2379 BFT_FREE(swap_buffer);
2380
2381 /* Treat fragmentation */
2382 init_particles = p_set->n_particles;
2383
2384 if (cs_glob_lagr_model->fragmentation == 1) {
2385 cs_lagr_fragmentation(dt[0],
2386 minimum_particle_diam,
2387 start_part,
2388 end_part - deleted_parts,
2389 init_particles,
2390 p_set->n_particles);
2391 }
2392 cs_lnum_t inserted_parts_frag = p_set->n_particles - init_particles;
2393
2394 cell_particle_idx[icell+1] = cell_particle_idx[icell]
2395 + inserted_parts_agglo + inserted_parts_frag;
2396 }
2397
2398 p_set->n_particles = enter_parts;
2399
2400 /* Introduce new particles (uniformly in the cell) */
2401 cs_lagr_new_v(p_set,
2402 n_occupied_cells,
2403 occupied_cell_ids,
2404 cell_particle_idx);
2405 p_set->n_particles += cell_particle_idx[n_occupied_cells];
2406
2407 BFT_FREE(cell_particle_idx);
2408 }
2409
2410 BFT_FREE(occupied_cell_ids);
2411 BFT_FREE(particle_list);
2412
2413 /* Reverse coupling: compute source terms
2414 -------------------------------------- */
2415
2416 if ( cs_glob_lagr_time_scheme->iilagr == CS_LAGR_TWOWAY_COUPLING
2417 && cs_glob_lagr_time_step->nor == cs_glob_lagr_time_scheme->t_order)
2418 cs_lagr_coupling(taup, tempct, tsfext, cpgd1, cpgd2, cpght);
2419
2420 /* Deallocate arrays whose size is based on p_set->n_particles
2421 (which may change next) */
2422
2423 BFT_FREE(tlag);
2424 BFT_FREE(taup);
2425 BFT_FREE(piil);
2426 BFT_FREE(bx);
2427
2428 if (cs_glob_lagr_time_scheme->iilagr == CS_LAGR_TWOWAY_COUPLING)
2429 BFT_FREE(tsfext);
2430
2431 if ( cs_glob_lagr_time_scheme->iilagr == CS_LAGR_TWOWAY_COUPLING
2432 && lagr_model->physical_model == CS_LAGR_PHYS_COAL
2433 && cs_glob_lagr_source_terms->ltsthe == 1) {
2434 BFT_FREE(cpgd1);
2435 BFT_FREE(cpgd2);
2436 BFT_FREE(cpght);
2437 }
2438
2439 if ( ( lagr_model->physical_model == CS_LAGR_PHYS_HEAT
2440 && cs_glob_lagr_specific_physics->itpvar == 1)
2441 || lagr_model->physical_model == CS_LAGR_PHYS_COAL)
2442 BFT_FREE(tempct);
2443
2444 if (cs_glob_lagr_brownian->lamvbr == 1)
2445 BFT_FREE(terbru);
2446
2447 p_set->n_particles += nresnew;
2448
2449 /* Location of particles - boundary conditions for particle positions
2450 ------------------------------------------------------------------ */
2451
2452 if (cs_glob_lagr_time_step->nor == 1) {
2453
2454 /* In unsteady case, reset boundary statistics */
2455
2456 if ( cs_glob_lagr_time_scheme->isttio == 0
2457 || ( cs_glob_lagr_time_scheme->isttio == 1
2458 && cs_glob_time_step->nt_cur
2459 <= cs_glob_lagr_stat_options->nstist)) {
2460
2461 lag_bdi->tstatp = 0.0;
2462 lag_bdi->npstf = 0;
2463
2464 for (int ii = 0; ii < cs_glob_lagr_dim->n_boundary_stats; ii++) {
2465
2466 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++)
2467 bound_stat[ii * n_b_faces + ifac] = 0.0;
2468
2469 }
2470
2471 }
2472
2473 lag_bdi->tstatp += cs_glob_lagr_time_step->dtp;
2474 lag_bdi->npstf++;
2475 lag_bdi->npstft++;
2476
2477 cs_lagr_tracking_particle_movement(vislen);
2478
2479 }
2480
2481 /* Update residence time */
2482
2483 if (cs_glob_lagr_time_step->nor == cs_glob_lagr_time_scheme->t_order) {
2484
2485 for (cs_lnum_t npt = 0; npt < p_set->n_particles; npt++) {
2486
2487 if (cs_lagr_particles_get_lnum(p_set, npt, CS_LAGR_CELL_ID) >= 0) {
2488 cs_real_t res_time
2489 = cs_lagr_particles_get_real(p_set, npt,
2490 CS_LAGR_RESIDENCE_TIME)
2491 + cs_glob_lagr_time_step->dtp;
2492 cs_lagr_particles_set_real
2493 (p_set, npt, CS_LAGR_RESIDENCE_TIME, res_time);
2494 }
2495
2496 }
2497
2498 }
2499
2500 /* Compute adhesion for reentrainement model
2501 ----------------------------------------- */
2502
2503 if (lagr_model->resuspension > 0)
2504 cs_lagr_resuspension();
2505
2506 /* Compute statistics
2507 ------------------ */
2508
2509 /* Calculation of consolidation:
2510 * linear increase of the consolidation height with deposit time */
2511
2512 if (cs_glob_lagr_consolidation_model->iconsol > 0) {
2513
2514 for (cs_lnum_t npt = 0; npt < p_set->n_particles; npt++) {
2515
2516 if (cs_lagr_particles_get_flag(p_set, npt, deposition_mask)) {
2517
2518 cs_real_t p_depo_time
2519 = cs_lagr_particles_get_real(p_set, npt, CS_LAGR_DEPO_TIME);
2520 cs_real_t p_consol_height
2521 = CS_MIN(cs_lagr_particles_get_real(p_set, npt, CS_LAGR_HEIGHT),
2522 cs_glob_lagr_consolidation_model->rate_consol * p_depo_time);
2523 cs_lagr_particles_set_real(p_set, npt,
2524 CS_LAGR_CONSOL_HEIGHT, p_consol_height);
2525 cs_lagr_particles_set_real(p_set, npt, CS_LAGR_DEPO_TIME,
2526 p_depo_time + cs_glob_lagr_time_step->dtp);
2527 }
2528 else {
2529 cs_lagr_particles_set_real(p_set, npt, CS_LAGR_CONSOL_HEIGHT, 0.0);
2530 }
2531
2532 }
2533
2534 }
2535
2536 if ( cs_glob_lagr_time_step->nor == cs_glob_lagr_time_scheme->t_order
2537 && cs_glob_time_step->nt_cur >= cs_glob_lagr_stat_options->idstnt)
2538 cs_lagr_stat_update();
2539
2540 /* Statistics for clogging */
2541
2542 if ( cs_glob_lagr_time_step->nor == cs_glob_lagr_time_scheme->t_order
2543 && lagr_model->clogging == 1
2544 && cs_glob_lagr_consolidation_model->iconsol == 1) {
2545
2546 /* Height and time of deposit */
2547
2548 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++) {
2549 bound_stat[lag_bdi->iclogt * n_b_faces + ifac] = 0.0;
2550 bound_stat[lag_bdi->iclogh * n_b_faces + ifac] = 0.0;
2551 bound_stat[lag_bdi->ihdiam * n_b_faces + ifac] = 0.0;
2552 }
2553
2554 for (cs_lnum_t npt = 0; npt < p_set->n_particles; npt++) {
2555
2556 if (cs_lagr_particles_get_flag(p_set, npt, deposition_mask)) {
2557
2558 cs_lnum_t face_id = cs_lagr_particles_get_lnum
2559 (p_set, npt, CS_LAGR_NEIGHBOR_FACE_ID);
2560
2561 cs_real_t p_diam = cs_lagr_particles_get_real
2562 (p_set, npt, CS_LAGR_DIAMETER);
2563
2564 bound_stat[lag_bdi->iclogt * n_b_faces + face_id]
2565 += cs_lagr_particles_get_real(p_set, npt, CS_LAGR_DEPO_TIME);
2566 bound_stat[lag_bdi->iclogh * n_b_faces + face_id]
2567 += cs_lagr_particles_get_real(p_set, npt, CS_LAGR_CONSOL_HEIGHT)
2568 * cs_math_pi * cs_math_sq(p_diam) * 0.25 / surfbn[face_id];
2569
2570 }
2571
2572 }
2573
2574 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++) {
2575
2576 if (bound_stat[lag_bdi->inclg * n_b_faces + ifac] > 0) {
2577
2578 bound_stat[lag_bdi->iclogt * n_b_faces + ifac]
2579 /= bound_stat[lag_bdi->inclg * n_b_faces + ifac];
2580 bound_stat[lag_bdi->ihdiam * n_b_faces + ifac]
2581 = bound_stat[lag_bdi->ihsum * n_b_faces + ifac]
2582 / bound_stat[lag_bdi->inclgt * n_b_faces + ifac];
2583
2584 }
2585 else if (bound_stat[lag_bdi->inclg * n_b_faces + ifac] <= 0) {
2586
2587 bound_stat[lag_bdi->iclogt * n_b_faces + ifac] = 0.0;
2588 bound_stat[lag_bdi->ihdiam * n_b_faces + ifac] = 0.0;
2589
2590 }
2591 else {
2592
2593 /* FIXME */
2594
2595 bft_printf(" ** LAGRANGIAN MODULE:\n"
2596 " ** Error in cs_lagr.c: inclg < 0 ! \n"
2597 "---------------------------------------------\n\n\n"
2598 "** Ifac = %ld and inclg = %g\n"
2599 "-------------------------------------------------\n",
2600 (long)ifac, bound_stat[lag_bdi->inclg * n_b_faces + ifac]);
2601 }
2602
2603 }
2604
2605 }
2606
2607 /* Poisson equation
2608 ---------------- */
2609
2610 if ( cs_glob_lagr_time_step->nor == cs_glob_lagr_time_scheme->t_order
2611 && cs_glob_lagr_time_scheme->ilapoi == 1)
2612 cs_lagr_poisson(itypfb);
2613
2614 /* Loop again ?
2615 ---------- */
2616
2617 if ( cs_glob_lagr_time_scheme->t_order != 2
2618 || cs_glob_lagr_time_step->nor != 1)
2619 go_on = false; //exit the while loop
2620
2621 }
2622
2623 } /* end if number of particles > 0 */
2624
2625 else if (cs_glob_time_step->nt_cur >= cs_glob_lagr_stat_options->idstnt)
2626 cs_lagr_stat_update();
2627
2628 /* Optional user modification of Lagrangian variables at end of iteration */
2629
2630 cs_user_lagr_extra_operations(dt);
2631
2632 /* Update particle counter */
2633 /*-------------------------*/
2634
2635 part_c = cs_lagr_update_particle_counter();
2636 part_c->n_g_cumulative_total += part_c->n_g_new;
2637 part_c->n_g_cumulative_failed += part_c->n_g_failed;
2638
2639 /* Logging
2640 ------- */
2641
2642 int modntl;
2643 if (ipass == 1)
2644 modntl = 0;
2645
2646 else if (cs_glob_lagr_log_frequency_n > 0)
2647 modntl = ts->nt_cur % cs_glob_lagr_log_frequency_n;
2648
2649 else if (cs_glob_lagr_log_frequency_n == -1 && ts->nt_cur == ts->nt_max)
2650 modntl = 0;
2651
2652 else
2653 modntl = 1;
2654
2655 if (modntl == 0)
2656 cs_lagr_print(ts->t_cur);
2657
2658 /* Free memory */
2659
2660 if ( lagr_model->deposition == 1
2661 && ts->nt_cur == ts->nt_max)
2662 BFT_FREE(vislen);
2663
2664 if ( lagr_model->clogging == 1
2665 || lagr_model->roughness == 1
2666 || lagr_model->dlvo == 1)
2667 BFT_FREE(tempp);
2668 }
2669
2670 /*----------------------------------------------------------------------------*/
2671
2672 END_C_DECLS
2673