1 /*============================================================================
2 * Methods for particle localization
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 particle tracking
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_boundary_conditions.h"
55 #include "cs_equation_iterative_solve.h"
56 #include "cs_mesh.h"
57 #include "cs_mesh_quantities.h"
58 #include "cs_parameters.h"
59 #include "cs_gradient.h"
60 #include "cs_face_viscosity.h"
61
62 #include "cs_lagr.h"
63 #include "cs_lagr_tracking.h"
64 #include "cs_lagr_stat.h"
65
66 /*----------------------------------------------------------------------------
67 * Header for the current file
68 *----------------------------------------------------------------------------*/
69
70 #include "cs_lagr_poisson.h"
71
72 /*----------------------------------------------------------------------------*/
73
74 BEGIN_C_DECLS
75
76 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
77
78 /*============================================================================
79 * Private function definitions
80 *============================================================================*/
81
82 /* -------------------------------------------------------
83 * CALCUL DE LA DIVERGENCE D'UN VECTEUR
84 * (On ne s'embete pas, on appelle 3 fois le gradient)
85 * ------------------------------------------------------- */
86
87 static void
diverv(cs_real_t * diverg,cs_real_3_t * u,cs_real_3_t * coefa,cs_real_33_t * coefb)88 diverv (cs_real_t *diverg,
89 cs_real_3_t *u,
90 cs_real_3_t *coefa,
91 cs_real_33_t *coefb)
92 {
93
94 /* ====================================================================
95 * 1. INITIALISATIONS
96 * ====================================================================*/
97 cs_lnum_t ncelet = cs_glob_mesh->n_cells_with_ghosts;
98 cs_lnum_t ncel = cs_glob_mesh->n_cells;
99
100 /* Allocate work arrays */
101 cs_real_33_t *grad;
102 BFT_MALLOC(grad, ncelet, cs_real_33_t);
103
104 /* ====================================================================
105 * Calcul du gradient de U
106 * ====================================================================*/
107
108 cs_halo_type_t halo_type = CS_HALO_STANDARD;
109 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
110
111 cs_gradient_type_by_imrgra(cs_glob_space_disc->imrgra,
112 &gradient_type,
113 &halo_type);
114
115 cs_gradient_vector("Work array",
116 gradient_type,
117 halo_type,
118 1, /* inc */
119 100, /* n_r_sweeps, */
120 2, /* iwarnp */
121 -1, /* imligp */
122 1e-8, /* epsrgp */
123 1.5, /* climgp */
124 (const cs_real_3_t *)coefa,
125 (const cs_real_33_t *)coefb,
126 u,
127 NULL, /* weighted gradient */
128 NULL, /* cpl */
129 grad);
130
131 /* ====================================================================
132 * Calcul de la divergence du vecteur
133 * ====================================================================*/
134
135 for (cs_lnum_t iel = 0; iel < ncel; iel++)
136 diverg[iel] = grad[iel][0][0] + grad[iel][1][1] + grad[iel][2][2];
137
138 /* Free memory */
139 BFT_FREE(grad);
140 }
141
142 /*-------------------------------------------------------------------
143 * RESOLUTION D'UNE EQUATION DE POISSON
144 * div[ALPHA grad(PHI)] = div(ALPHA <Up>)
145 *-------------------------------------------------------------------*/
146
147 static void
_lageqp(cs_real_t * vitessel,cs_real_t * alphal,cs_real_t * phi,const int itypfb[])148 _lageqp(cs_real_t *vitessel,
149 cs_real_t *alphal,
150 cs_real_t *phi,
151 const int itypfb[])
152 {
153 /* ====================================================================
154 * 1. INITIALISATION
155 * ====================================================================*/
156
157 const cs_mesh_t *m = cs_glob_mesh;
158 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
159 cs_lnum_t ncelet = m->n_cells_with_ghosts;
160 cs_lnum_t ncel = m->n_cells;
161 cs_lnum_t nfac = m->n_i_faces;
162 cs_lnum_t nfabor = m->n_b_faces;
163
164 cs_real_t *viscf, *viscb;
165 cs_real_t *smbrs;
166 cs_real_t *rovsdt;
167 cs_real_t *fmala, *fmalb;
168 cs_real_t *phia;
169 cs_real_t *dpvar;
170
171 /* Allocate temporary arrays */
172
173 BFT_MALLOC(viscf , nfac , cs_real_t);
174 BFT_MALLOC(viscb , nfabor, cs_real_t);
175 BFT_MALLOC(smbrs , ncelet, cs_real_t);
176 BFT_MALLOC(rovsdt, ncelet, cs_real_t);
177 BFT_MALLOC(fmala , nfac , cs_real_t);
178 BFT_MALLOC(fmalb , nfabor, cs_real_t);
179 BFT_MALLOC(phia , ncelet, cs_real_t);
180 BFT_MALLOC(dpvar , ncelet, cs_real_t);
181
182 /* Allocate work arrays */
183 cs_real_3_t *w;
184 BFT_MALLOC(w, ncelet, cs_real_3_t);
185
186 bft_printf(_(" ** RESOLUTION POUR LA VARIABLE Pressure correction"));
187
188 /* ====================================================================
189 * 2. TERMES SOURCES
190 * ==================================================================== */
191
192 /* --> Initialisation */
193
194 for (cs_lnum_t iel = 0; iel < ncel; iel++) {
195
196 smbrs[iel] = 0.0;
197 rovsdt[iel] = 0.0;
198 phi[iel] = 0.0;
199 phia[iel] = 0.0;
200
201 }
202
203 /* "VITESSE" DE DIFFUSION FACE */
204 cs_face_viscosity(m,
205 fvq,
206 cs_glob_space_disc->imvisf,
207 alphal,
208 viscf,
209 viscb);
210
211 /* CALCUL de div(Alpha Up) avant correction */
212 for (cs_lnum_t iel = 0; iel < ncel; iel++) {
213
214 for (cs_lnum_t isou = 0; isou < 3; isou++)
215 w[isou][iel] = -vitessel[isou + iel * 3] * alphal[iel];
216
217 }
218
219 /* --> Calcul du gradient de W1 */
220 /* ======================== */
221 /* Allocate temporary arrays */
222 cs_real_3_t *coefaw;
223 cs_real_33_t *coefbw;
224
225 BFT_MALLOC(coefaw, nfabor, cs_real_3_t);
226 BFT_MALLOC(coefbw, nfabor, cs_real_33_t);
227
228 for (cs_lnum_t ifac = 0; ifac < nfabor; ifac++) {
229
230 cs_lnum_t iel = m->b_face_cells[ifac];
231
232 for (cs_lnum_t isou = 0; isou < 3; isou++)
233 coefaw[isou][ifac] = w[isou][iel];
234
235 }
236
237 for (cs_lnum_t ifac = 0; ifac < nfabor; ifac++) {
238
239 for (cs_lnum_t isou = 0; isou < 3; isou++) {
240
241 for (cs_lnum_t jsou = 0; jsou < 3; jsou++)
242 coefbw[jsou][isou][ifac] = 0.0;
243
244 }
245
246 }
247
248 diverv(smbrs, w, coefaw, coefbw);
249
250 /* Free memory */
251 BFT_FREE (coefaw);
252 BFT_FREE (coefbw);
253
254 /* --> Conditions aux limites sur PHI */
255 /* ============================== */
256 /* Allocate temporary arrays */
257 cs_real_t *coefap, *coefbp;
258 cs_real_t *cofafp, *cofbfp;
259 BFT_MALLOC(coefap, nfabor, cs_real_t);
260 BFT_MALLOC(coefbp, nfabor, cs_real_t);
261 BFT_MALLOC(cofafp, nfabor, cs_real_t);
262 BFT_MALLOC(cofbfp, nfabor, cs_real_t);
263
264 for (cs_lnum_t ifac = 0; ifac < nfabor; ifac++) {
265
266 cs_lnum_t iel = m->b_face_cells[ifac];
267 cs_real_t hint = alphal[iel] / fvq->b_dist[ifac];
268
269 if ( itypfb[ifac] == CS_INLET
270 || itypfb[ifac] == CS_SMOOTHWALL
271 || itypfb[ifac] == CS_ROUGHWALL
272 || itypfb[ifac] == CS_SYMMETRY) {
273
274 /* Neumann Boundary Conditions */
275
276 cs_boundary_conditions_set_neumann_scalar(&coefap[ifac],
277 &cofafp[ifac],
278 &coefbp[ifac],
279 &cofbfp[ifac],
280 0.0,
281 hint);
282 coefap[ifac] = 0.0;
283 coefbp[ifac] = 1.0;
284
285 }
286 else if (itypfb[ifac] == CS_OUTLET) {
287
288 /* Dirichlet Boundary Condition */
289
290 cs_boundary_conditions_set_dirichlet_scalar(&coefap[ifac],
291 &cofafp[ifac],
292 &coefbp[ifac],
293 &cofbfp[ifac],
294 phia[iel],
295 hint,
296 -1);
297
298 }
299 else
300 bft_error
301 (__FILE__, __LINE__, 0,
302 _("\n%s (Lagrangian module):\n"
303 " unexpected boundary conditions for Phi."), __func__);
304
305 }
306
307 /* ====================================================================
308 * 3. RESOLUTION
309 * ==================================================================== */
310
311 int idtva0 = 0; /* No steady option here */
312
313 /* Cancel mass fluxes */
314
315 for (cs_lnum_t ifac = 0; ifac < nfac; ifac++) {
316 fmala[ifac] = 0.0;
317 fmalb[ifac] = 0.0;
318 }
319
320 /* In the theta-scheme case, set theta to 1 (order 1) */
321
322 cs_var_cal_opt_t var_cal_opt = cs_parameters_var_cal_opt_default();
323
324 var_cal_opt.verbosity = 2; /* quasi-debug at this stage, TODO clean */
325 var_cal_opt.iconv = 0; /* no convection, pure diffusion here */
326 var_cal_opt.istat = -1;
327 var_cal_opt.ndircl = 1;
328 var_cal_opt.idifft = -1;
329 var_cal_opt.isstpc = 0;
330 var_cal_opt.nswrgr = 10000;
331 var_cal_opt.nswrsm = 2;
332 var_cal_opt.imrgra = cs_glob_space_disc->imrgra;
333 var_cal_opt.imligr = 1;
334
335 cs_equation_iterative_solve_scalar(idtva0,
336 1, /* external sub-iteration? */
337 -1, /* field_id (not a field) */
338 "PoissonL", /* name */
339 0, /* iescap */
340 0, /* imucpp */
341 -1, /* normp */
342 &var_cal_opt,
343 phia, phia,
344 coefap, coefbp,
345 cofafp, cofbfp,
346 fmala, fmalb,
347 viscf, viscb,
348 viscf, viscb,
349 NULL, /* viscel */
350 NULL, /* weighf */
351 NULL, /* weighb */
352 0, /* icvflb (all upwind) */
353 NULL, /* icvfli */
354 rovsdt,
355 smbrs,
356 phi,
357 dpvar,
358 NULL, /* xcpp */
359 NULL); /* eswork */
360
361 /* Free memory */
362
363 BFT_FREE(viscf);
364 BFT_FREE(viscb);
365 BFT_FREE(smbrs);
366 BFT_FREE(rovsdt);
367 BFT_FREE(fmala);
368 BFT_FREE(fmalb);
369 BFT_FREE(coefap);
370 BFT_FREE(coefbp);
371 BFT_FREE(cofafp);
372 BFT_FREE(cofbfp);
373 BFT_FREE(phia);
374 BFT_FREE(w);
375 BFT_FREE(dpvar);
376 }
377
378 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
379
380 /*============================================================================
381 * Public function definitions
382 *============================================================================*/
383
384 /*-----------------------------------------------------------------*/
385 /*! \brief Solve Poisson equation for mean particle velocities
386 * and correct particle instantaneous velocities
387 *
388 * \param[in] itypfb boundary face type
389 */
390 /*-----------------------------------------------------------------*/
391
392 void
cs_lagr_poisson(const int itypfb[])393 cs_lagr_poisson(const int itypfb[])
394 {
395 cs_lnum_t ncel = cs_glob_mesh->n_cells;
396 cs_lnum_t ncelet = cs_glob_mesh->n_cells_with_ghosts;
397 cs_lnum_t nfabor = cs_glob_mesh->n_b_faces;
398
399 /* Allocate a temporary array */
400
401 cs_real_t *phil;
402 BFT_MALLOC(phil, ncelet, cs_real_t);
403
404 /* Initialization */
405
406 cs_lagr_particle_set_t *p_set = cs_lagr_get_particle_set();
407 const cs_lagr_attribute_map_t *p_am = p_set->p_am;
408
409 /* Means of global class */
410
411 int stat_type = cs_lagr_stat_type_from_attr_id(CS_LAGR_VELOCITY);
412
413 cs_field_t *mean_vel
414 = cs_lagr_stat_get_moment(stat_type,
415 CS_LAGR_STAT_GROUP_PARTICLE,
416 CS_LAGR_MOMENT_MEAN,
417 0,
418 -1);
419
420 cs_field_t *mean_fv
421 = cs_lagr_stat_get_moment(CS_LAGR_STAT_VOLUME_FRACTION,
422 CS_LAGR_STAT_GROUP_PARTICLE,
423 CS_LAGR_MOMENT_MEAN,
424 0,
425 -1);
426
427 cs_field_t *stat_weight = cs_lagr_stat_get_stat_weight(0);
428
429 _lageqp(mean_vel->val, mean_fv->val, phil, itypfb);
430
431 /* Compute gradient of phi corrector */
432
433 cs_real_3_t *grad;
434
435 BFT_MALLOC(grad, ncelet, cs_real_3_t);
436
437 cs_real_t *coefap, *coefbp;
438
439 BFT_MALLOC(coefap, nfabor, cs_real_t);
440 BFT_MALLOC(coefbp, nfabor, cs_real_t);
441
442 for (cs_lnum_t ifac = 0; ifac < nfabor; ifac++) {
443 cs_lnum_t iel = cs_glob_mesh->b_face_cells[ifac];
444 coefap[ifac] = phil[iel];
445 coefbp[ifac] = 0.0;
446 }
447
448 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
449 cs_halo_type_t halo_type = CS_HALO_STANDARD;
450
451 cs_gradient_type_by_imrgra(cs_glob_space_disc->imrgra,
452 &gradient_type,
453 &halo_type);
454
455 cs_gradient_scalar("Work array",
456 gradient_type,
457 halo_type,
458 1, /* inc */
459 1, /* recompute_cocg */
460 100, /* n_r_sweeps */
461 0, /* ignored */
462 0, /* hyd_p_flag */
463 1, /* w_stride */
464 2, /* iwarnp */
465 -1, /* imligp */
466 1e-8, /* epsrgp */
467 1.5, /* climgp */
468 NULL, /* f_ext */
469 coefap,
470 coefbp,
471 phil,
472 NULL, /* c_weight */
473 NULL, /* internal coupling */
474 grad);
475
476 BFT_FREE(coefap);
477 BFT_FREE(coefbp);
478 BFT_FREE(phil);
479
480 /* Correct mean velocities */
481
482 for (cs_lnum_t iel = 0; iel < ncel; iel++) {
483
484 if (stat_weight->val[iel] > cs_glob_lagr_stat_options->threshold) {
485 for (cs_lnum_t id = 0; id < 3; id++)
486 mean_vel->val[iel * 3 + id] += - grad[iel][id];
487 }
488
489 }
490
491 /* Correct instant velocities */
492
493 for (cs_lnum_t npt = 0; npt < p_set->n_particles; npt++) {
494
495 unsigned char *part = p_set->p_buffer + p_am->extents * npt;
496 cs_lnum_t iel = cs_lagr_particle_get_lnum(part, p_am, CS_LAGR_CELL_ID);
497
498 if (iel >= 0) {
499
500 cs_real_t *part_vel = cs_lagr_particle_attr(part, p_am, CS_LAGR_VELOCITY);
501
502 for (cs_lnum_t id = 0; id < 3; id++)
503 part_vel[id] += -grad[id][iel];
504
505 }
506
507 }
508
509 /* Free memory */
510
511 BFT_FREE(grad);
512
513 }
514
515 END_C_DECLS
516