1 /*============================================================================
2 * Divergence operators.
3 *============================================================================*/
4
5 /* This file is part of Code_Saturne, a general-purpose CFD tool.
6
7 Copyright (C) 1998-2021 EDF S.A.
8
9 This program is free software; you can redistribute it and/or modify it under
10 the terms of the GNU General Public License as published by the Free Software
11 Foundation; either version 2 of the License, or (at your option) any later
12 version.
13
14 This program is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
17 details.
18
19 You should have received a copy of the GNU General Public License along with
20 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
21 Street, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 /*----------------------------------------------------------------------------*/
24
25 #include "cs_defs.h"
26
27 /*----------------------------------------------------------------------------
28 * Standard C library headers
29 *----------------------------------------------------------------------------*/
30
31 #include <assert.h>
32 #include <errno.h>
33 #include <stdio.h>
34 #include <stdarg.h>
35 #include <string.h>
36 #include <math.h>
37 #include <float.h>
38
39 #if defined(HAVE_MPI)
40 #include <mpi.h>
41 #endif
42
43 /*----------------------------------------------------------------------------
44 * Local headers
45 *----------------------------------------------------------------------------*/
46
47 #include "bft_error.h"
48 #include "bft_mem.h"
49 #include "bft_printf.h"
50
51 #include "cs_blas.h"
52 #include "cs_halo.h"
53 #include "cs_halo_perio.h"
54 #include "cs_log.h"
55 #include "cs_mesh.h"
56 #include "cs_field.h"
57 #include "cs_field_pointer.h"
58 #include "cs_gradient.h"
59 #include "cs_ext_neighborhood.h"
60 #include "cs_mesh_quantities.h"
61 #include "cs_parameters.h"
62 #include "cs_porous_model.h"
63 #include "cs_prototypes.h"
64 #include "cs_timer.h"
65
66 /*----------------------------------------------------------------------------
67 * Header for the current file
68 *----------------------------------------------------------------------------*/
69
70 #include "cs_divergence.h"
71
72 /*----------------------------------------------------------------------------*/
73
74 BEGIN_C_DECLS
75
76 /*=============================================================================
77 * Additional Doxygen documentation
78 *============================================================================*/
79
80 /*! \file cs_divergence.c
81
82 */
83
84 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
85
86 /*=============================================================================
87 * Local Macro Definitions
88 *============================================================================*/
89
90 /*=============================================================================
91 * Local type definitions
92 *============================================================================*/
93
94 /*============================================================================
95 * Private function definitions
96 *============================================================================*/
97
98 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
99
100 /*============================================================================
101 * Public function definitions for Fortran API
102 *============================================================================*/
103
104 /*----------------------------------------------------------------------------
105 * Wrapper to cs_mass_flux
106 *----------------------------------------------------------------------------*/
107
CS_PROCF(inimav,INIMAV)108 void CS_PROCF (inimav, INIMAV)
109 (
110 const int *const f_id,
111 const int *const itypfl,
112 const int *const iflmb0,
113 const int *const init,
114 const int *const inc,
115 const int *const imrgra,
116 const int *const nswrgu,
117 const int *const imligu,
118 const int *const iwarnu,
119 const cs_real_t *const epsrgu,
120 const cs_real_t *const climgu,
121 const cs_real_t rom[],
122 const cs_real_t romb[],
123 const cs_real_3_t vel[],
124 const cs_real_3_t coefav[],
125 const cs_real_33_t coefbv[],
126 cs_real_t i_massflux[],
127 cs_real_t b_massflux[]
128 )
129 {
130 const cs_mesh_t *m = cs_glob_mesh;
131 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
132
133 cs_mass_flux(m,
134 fvq,
135 *f_id,
136 *itypfl,
137 *iflmb0,
138 *init,
139 *inc,
140 *imrgra,
141 *nswrgu,
142 *imligu,
143 *iwarnu,
144 *epsrgu,
145 *climgu,
146 rom,
147 romb,
148 vel,
149 coefav,
150 coefbv,
151 i_massflux,
152 b_massflux);
153 }
154
155 /*----------------------------------------------------------------------------
156 * Wrapper to cs_divergence
157 *----------------------------------------------------------------------------*/
158
CS_PROCF(divmas,DIVMAS)159 void CS_PROCF (divmas, DIVMAS)
160 (
161 const int *const init,
162 const cs_real_t i_massflux[],
163 const cs_real_t b_massflux[],
164 cs_real_t diverg[]
165 )
166 {
167 const cs_mesh_t *m = cs_glob_mesh;
168
169 cs_divergence(m,
170 *init,
171 i_massflux,
172 b_massflux,
173 diverg);
174 }
175
176 /*----------------------------------------------------------------------------
177 * Wrapper to cs_tensor_divergence
178 *----------------------------------------------------------------------------*/
179
CS_PROCF(divmat,DIVMAT)180 void CS_PROCF (divmat, DIVMAT)
181 (
182 const int *const init,
183 const cs_real_3_t i_massflux[],
184 const cs_real_3_t b_massflux[],
185 cs_real_3_t diverg[]
186 )
187 {
188 const cs_mesh_t *m = cs_glob_mesh;
189
190 cs_tensor_divergence(m,
191 *init,
192 i_massflux,
193 b_massflux,
194 diverg);
195 }
196
197 /*----------------------------------------------------------------------------
198 * Wrapper to cs_ext_force_flux
199 *----------------------------------------------------------------------------*/
200
CS_PROCF(projts,PROJTS)201 void CS_PROCF (projts, PROJTS)
202 (
203 const int *const init,
204 const int *const nswrgu,
205 const cs_real_3_t frcxt[],
206 const cs_real_t cofbfp[],
207 cs_real_t i_massflux[],
208 cs_real_t b_massflux[],
209 const cs_real_t i_visc[],
210 const cs_real_t b_visc[],
211 const cs_real_t viselx[],
212 const cs_real_t visely[],
213 const cs_real_t viselz[]
214 )
215 {
216 const cs_mesh_t *m = cs_glob_mesh;
217 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
218
219 cs_ext_force_flux(m,
220 fvq,
221 *init,
222 *nswrgu,
223 frcxt,
224 cofbfp,
225 i_massflux,
226 b_massflux,
227 i_visc,
228 b_visc,
229 viselx,
230 visely,
231 viselz);
232 }
233
234 /*----------------------------------------------------------------------------
235 * Wrapper to cs_ext_force_anisotropic_flux
236 *----------------------------------------------------------------------------*/
237
CS_PROCF(projtv,PROJTV)238 void CS_PROCF (projtv, PROJTV)
239 (
240 const int *const init,
241 const int *const nswrgu,
242 const int *const ircflp,
243 const cs_real_3_t frcxt[],
244 const cs_real_t cofbfp[],
245 const cs_real_t i_visc[],
246 const cs_real_t b_visc[],
247 cs_real_6_t viscel[],
248 const cs_real_2_t weighf[],
249 cs_real_t i_massflux[],
250 cs_real_t b_massflux[])
251 {
252 const cs_mesh_t *m = cs_glob_mesh;
253 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
254
255 cs_ext_force_anisotropic_flux(m,
256 fvq,
257 *init,
258 *nswrgu,
259 *ircflp,
260 frcxt,
261 cofbfp,
262 i_visc,
263 b_visc,
264 viscel,
265 weighf,
266 i_massflux,
267 b_massflux);
268 }
269
270 /*----------------------------------------------------------------------------
271 * Wrapper to cs_tensor_face_flux
272 *----------------------------------------------------------------------------*/
273
CS_PROCF(divrij,DIVRIJ)274 void CS_PROCF (divrij, DIVRIJ)
275 (
276 const int *const f_id,
277 const int *const itypfl,
278 const int *const iflmb0,
279 const int *const init,
280 const int *const inc,
281 const int *const imrgra,
282 const int *const nswrgu,
283 const int *const imligu,
284 const int *const iwarnu,
285 const cs_real_t *const epsrgu,
286 const cs_real_t *const climgu,
287 const cs_real_t rom[],
288 const cs_real_t romb[],
289 const cs_real_6_t tensorvel[],
290 const cs_real_6_t coefat[],
291 const cs_real_66_t coefbt[],
292 cs_real_3_t i_massflux[],
293 cs_real_3_t b_massflux[])
294 {
295 const cs_mesh_t *m = cs_glob_mesh;
296 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
297
298 cs_tensor_face_flux(m,
299 fvq,
300 *f_id,
301 *itypfl,
302 *iflmb0,
303 *init,
304 *inc,
305 *imrgra,
306 *nswrgu,
307 *imligu,
308 *iwarnu,
309 *epsrgu,
310 *climgu,
311 rom,
312 romb,
313 tensorvel,
314 coefat,
315 coefbt,
316 i_massflux,
317 b_massflux);
318 }
319
320 /*============================================================================
321 * Public function definitions
322 *============================================================================*/
323
324 /*----------------------------------------------------------------------------*/
325 /*!
326 * \brief Add \f$ \rho \vect{u} \cdot \vect{s}_\ij\f$ to
327 * the mass flux \f$ \dot{m}_\ij \f$.
328 *
329 * For the reconstruction, \f$ \gradt \left(\rho \vect{u} \right) \f$ is
330 * computed with the following approximated boundary conditions:
331 * - \f$ \vect{a}_{\rho u} = \rho_\fib \vect{a}_u \f$
332 * - \f$ \tens{b}_{\rho u} = \tens{b}_u \f$
333 *
334 * For the mass flux at the boundary we have:
335 * \f[
336 * \dot{m}_\ib = \left[ \rho_\fib \vect{a}_u + \rho_\fib \tens{b}_u \vect{u}
337 * + \tens{b}_u \left(\gradt \vect{u} \cdot \vect{\centi \centip}\right)\right]
338 * \cdot \vect{s}_\ij
339 * \f]
340 * The last equation uses some approximations detailed in the theory guide.
341 *
342 * \param[in] m pointer to mesh
343 * \param[in] fvq pointer to finite volume quantities
344 * \param[in] f_id field id (or -1)
345 * \param[in] itypfl indicator (take rho into account or not)
346 * - 1 compute \f$ \rho\vect{u}\cdot\vect{s} \f$
347 * - 0 compute \f$ \vect{u}\cdot\vect{s} \f$
348 * \param[in] iflmb0 the mass flux is set to 0 on symmetries if = 1
349 * \param[in] init the mass flux is initialized to 0 if > 0
350 * \param[in] inc indicator
351 * - 0 solve an increment
352 * - 1 otherwise
353 * \param[in] imrgra indicator
354 * - 0 iterative gradient
355 * - 1 least square gradient
356 * \param[in] nswrgu number of sweeps for the reconstruction
357 * of the gradients
358 * \param[in] imligu clipping gradient method
359 * - < 0 no clipping
360 * - = 0 thanks to neighbooring gradients
361 * - = 1 thanks to the mean gradient
362 * \param[in] iwarnu verbosity
363 * \param[in] epsrgu relative precision for the gradient
364 * reconstruction
365 * \param[in] climgu clipping coefficient for the computation of
366 * the gradient
367 * \param[in] rom cell density
368 * \param[in] romb density at boundary faces
369 * \param[in] vel vector variable
370 * \param[in] coefav boundary condition array for the variable
371 * (explicit part - vector array )
372 * \param[in] coefbv boundary condition array for the variable
373 * (implicit part - 3x3 tensor array)
374 * \param[in,out] i_massflux mass flux at interior faces \f$ \dot{m}_\fij \f$
375 * \param[in,out] b_massflux mass flux at boundary faces \f$ \dot{m}_\fib \f$
376 */
377 /*----------------------------------------------------------------------------*/
378
379 void
cs_mass_flux(const cs_mesh_t * m,const cs_mesh_quantities_t * fvq,int f_id,int itypfl,int iflmb0,int init,int inc,int imrgra,int nswrgu,int imligu,int iwarnu,double epsrgu,double climgu,const cs_real_t rom[],const cs_real_t romb[],const cs_real_3_t vel[],const cs_real_3_t coefav[],const cs_real_33_t coefbv[],cs_real_t * restrict i_massflux,cs_real_t * restrict b_massflux)380 cs_mass_flux(const cs_mesh_t *m,
381 const cs_mesh_quantities_t *fvq,
382 int f_id,
383 int itypfl,
384 int iflmb0,
385 int init,
386 int inc,
387 int imrgra,
388 int nswrgu,
389 int imligu,
390 int iwarnu,
391 double epsrgu,
392 double climgu,
393 const cs_real_t rom[],
394 const cs_real_t romb[],
395 const cs_real_3_t vel[],
396 const cs_real_3_t coefav[],
397 const cs_real_33_t coefbv[],
398 cs_real_t *restrict i_massflux,
399 cs_real_t *restrict b_massflux)
400 {
401 const cs_halo_t *halo = m->halo;
402
403 const cs_lnum_t n_cells = m->n_cells;
404 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
405 const int n_i_groups = m->i_face_numbering->n_groups;
406 const int n_i_threads = m->i_face_numbering->n_threads;
407 const int n_b_groups = m->b_face_numbering->n_groups;
408 const int n_b_threads = m->b_face_numbering->n_threads;
409 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
410 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
411
412 const cs_lnum_2_t *restrict i_face_cells
413 = (const cs_lnum_2_t *restrict)m->i_face_cells;
414 const cs_lnum_t *restrict b_face_cells
415 = (const cs_lnum_t *restrict)m->b_face_cells;
416 const cs_real_t *restrict weight = fvq->weight;
417 const cs_real_3_t *restrict i_f_face_normal
418 = (const cs_real_3_t *restrict)fvq->i_f_face_normal;
419 const cs_real_3_t *restrict b_f_face_normal
420 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
421 cs_real_2_t *i_f_face_factor;
422 cs_real_t *b_f_face_factor;
423
424 /* Local variables */
425
426 /* Discontinuous porous treatment */
427 cs_real_2_t _i_f_face_factor = {1., 1.};
428 cs_real_t _b_f_face_factor = 1.;
429 int is_p = 0; /* Is porous? */
430
431 if (cs_glob_porous_model == 3) {
432 i_f_face_factor = fvq->i_f_face_factor;
433 b_f_face_factor = fvq->b_f_face_factor;
434 is_p = 1;
435 }
436 else {
437 i_f_face_factor = &_i_f_face_factor;
438 b_f_face_factor = &_b_f_face_factor;
439 }
440
441 const cs_real_3_t *restrict diipb
442 = (const cs_real_3_t *restrict)fvq->diipb;
443 const cs_real_3_t *restrict dofij
444 = (const cs_real_3_t *restrict)fvq->dofij;
445
446 char var_name[64];
447
448 cs_real_3_t *qdm, *f_momentum, *coefaq;
449 cs_real_33_t *grdqdm;
450
451 cs_field_t *f;
452
453 BFT_MALLOC(qdm, n_cells_ext, cs_real_3_t);
454 BFT_MALLOC(f_momentum, m->n_b_faces, cs_real_3_t);
455 BFT_MALLOC(coefaq, m->n_b_faces, cs_real_3_t);
456
457 /*==========================================================================
458 1. Initialization
459 ==========================================================================*/
460
461 /* Choose gradient type */
462
463 cs_halo_type_t halo_type = CS_HALO_STANDARD;
464 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
465
466 cs_gradient_type_by_imrgra(imrgra,
467 &gradient_type,
468 &halo_type);
469
470 if (f_id != -1) {
471 f = cs_field_by_id(f_id);
472 snprintf(var_name, 63, "%s", f->name);
473 }
474 else {
475 strncpy(var_name, "[momentum]", 63);
476 }
477 var_name[63] = '\0';
478
479 /* Momentum computation */
480
481 if (init == 1) {
482 # pragma omp parallel for
483 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
484 i_massflux[face_id] = 0.;
485 }
486 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
487 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
488 b_massflux[face_id] = 0.;
489 }
490
491 } else if (init != 0) {
492 bft_error(__FILE__, __LINE__, 0,
493 _("invalid value of init"));
494 }
495
496 /* Porosity fields */
497 cs_field_t *fporo = cs_field_by_name_try("porosity");
498 cs_field_t *ftporo = cs_field_by_name_try("tensorial_porosity");
499
500 cs_real_t *porosi = NULL;
501 cs_real_6_t *porosf = NULL;
502
503 if (cs_glob_porous_model == 1 || cs_glob_porous_model == 2) {
504 porosi = fporo->val;
505 if (ftporo != NULL) {
506 porosf = (cs_real_6_t *)ftporo->val;
507 }
508 }
509
510 /* Standard mass flux */
511 if (itypfl == 1) {
512
513 /* Without porosity */
514 if (porosi == NULL) {
515 # pragma omp parallel for
516 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
517 for (int isou = 0; isou < 3; isou++) {
518 qdm[cell_id][isou] = rom[cell_id]*vel[cell_id][isou];
519 }
520 }
521 /* With porosity */
522 } else if (porosi != NULL && porosf == NULL) {
523 # pragma omp parallel for
524 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
525 for (int isou = 0; isou < 3; isou++) {
526 qdm[cell_id][isou] = rom[cell_id]*vel[cell_id][isou]*porosi[cell_id];
527 }
528 }
529 /* With anisotropic porosity */
530 } else if (porosi != NULL && porosf != NULL) {
531 # pragma omp parallel for
532 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
533 qdm[cell_id][0] = ( porosf[cell_id][0]*vel[cell_id][0]
534 + porosf[cell_id][3]*vel[cell_id][1]
535 + porosf[cell_id][5]*vel[cell_id][2] )
536 * rom[cell_id];
537 qdm[cell_id][1] = ( porosf[cell_id][3]*vel[cell_id][0]
538 + porosf[cell_id][1]*vel[cell_id][1]
539 + porosf[cell_id][4]*vel[cell_id][2] )
540 * rom[cell_id];
541 qdm[cell_id][2] = ( porosf[cell_id][5]*vel[cell_id][0]
542 + porosf[cell_id][4]*vel[cell_id][1]
543 + porosf[cell_id][2]*vel[cell_id][2] )
544 * rom[cell_id];
545 }
546 }
547
548 /* Velocity flux */
549 } else {
550
551 /* Without porosity */
552 if (porosi == NULL) {
553 # pragma omp parallel for
554 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
555 for (int isou = 0; isou < 3; isou++) {
556 qdm[cell_id][isou] = vel[cell_id][isou];
557 }
558 }
559 /* With porosity */
560 } else if (porosi != NULL && porosf == NULL) {
561 # pragma omp parallel for
562 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
563 for (int isou = 0; isou < 3; isou++) {
564 qdm[cell_id][isou] = vel[cell_id][isou]*porosi[cell_id];
565 }
566 }
567 /* With anisotropic porosity */
568 } else if (porosi != NULL && porosf != NULL) {
569 # pragma omp parallel for
570 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
571 qdm[cell_id][0] = porosf[cell_id][0]*vel[cell_id][0]
572 + porosf[cell_id][3]*vel[cell_id][1]
573 + porosf[cell_id][5]*vel[cell_id][2];
574 qdm[cell_id][1] = porosf[cell_id][3]*vel[cell_id][0]
575 + porosf[cell_id][1]*vel[cell_id][1]
576 + porosf[cell_id][4]*vel[cell_id][2];
577 qdm[cell_id][2] = porosf[cell_id][5]*vel[cell_id][0]
578 + porosf[cell_id][4]*vel[cell_id][1]
579 + porosf[cell_id][2]*vel[cell_id][2];
580 }
581 }
582 }
583
584 /* ---> Periodicity and parallelism treatment */
585
586 if (halo != NULL) {
587 cs_halo_sync_var_strided(halo, halo_type, (cs_real_t *)qdm, 3);
588 if (cs_glob_mesh->n_init_perio > 0)
589 cs_halo_perio_sync_var_vect(halo, halo_type, (cs_real_t *)qdm, 3);
590 }
591
592 /* Standard mass flux */
593 if (itypfl == 1) {
594
595 /* Without porosity */
596 if (porosi == NULL) {
597 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
598 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
599 cs_lnum_t cell_id = b_face_cells[face_id];
600 for (int isou = 0; isou < 3; isou++) {
601 coefaq[face_id][isou] = romb[face_id]*coefav[face_id][isou];
602 f_momentum[face_id][isou] = romb[face_id]*vel[cell_id][isou];
603 }
604 }
605 } /* With porosity */
606 else if (porosi != NULL && porosf == NULL) {
607 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
608 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
609 cs_lnum_t cell_id = b_face_cells[face_id];
610 for (int isou = 0; isou < 3; isou++) {
611 coefaq[face_id][isou] = romb[face_id]
612 *coefav[face_id][isou]*porosi[cell_id];
613 f_momentum[face_id][isou] = romb[face_id]*vel[cell_id][isou]
614 *porosi[cell_id];
615 }
616 }
617
618 } /* With anisotropic porosity */
619 else if (porosi != NULL && porosf != NULL) {
620 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
621 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
622 cs_lnum_t cell_id = b_face_cells[face_id];
623 coefaq[face_id][0] = ( porosf[cell_id][0]*coefav[face_id][0]
624 + porosf[cell_id][3]*coefav[face_id][1]
625 + porosf[cell_id][5]*coefav[face_id][2] )
626 * romb[face_id];
627 coefaq[face_id][1] = ( porosf[cell_id][3]*coefav[face_id][0]
628 + porosf[cell_id][1]*coefav[face_id][1]
629 + porosf[cell_id][4]*coefav[face_id][2] )
630 * romb[face_id];
631 coefaq[face_id][2] = ( porosf[cell_id][5]*coefav[face_id][0]
632 + porosf[cell_id][4]*coefav[face_id][1]
633 + porosf[cell_id][2]*coefav[face_id][2] )
634 * romb[face_id];
635 f_momentum[face_id][0] = ( porosf[cell_id][0]*vel[cell_id][0]
636 + porosf[cell_id][3]*vel[cell_id][1]
637 + porosf[cell_id][5]*vel[cell_id][2] )
638 * romb[face_id];
639 f_momentum[face_id][1] = ( porosf[cell_id][3]*vel[cell_id][0]
640 + porosf[cell_id][1]*vel[cell_id][1]
641 + porosf[cell_id][4]*vel[cell_id][2] )
642 * romb[face_id];
643 f_momentum[face_id][2] = ( porosf[cell_id][5]*vel[cell_id][0]
644 + porosf[cell_id][4]*vel[cell_id][1]
645 + porosf[cell_id][2]*vel[cell_id][2] )
646 * romb[face_id];
647 }
648 }
649
650 /* Velocity flux */
651 } else {
652
653 /* Without porosity */
654 if (porosi == NULL) {
655 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
656 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
657 cs_lnum_t cell_id = b_face_cells[face_id];
658 for (int isou = 0; isou < 3; isou++) {
659 coefaq[face_id][isou] = coefav[face_id][isou];
660 f_momentum[face_id][isou] = vel[cell_id][isou];
661 }
662 }
663 } /* With porosity */
664 else if (porosi != NULL && porosf == NULL) {
665 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
666 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
667 cs_lnum_t cell_id = b_face_cells[face_id];
668 for (int isou = 0; isou < 3; isou++) {
669 coefaq[face_id][isou] = coefav[face_id][isou]*porosi[cell_id];
670 f_momentum[face_id][isou] = vel[cell_id][isou]*porosi[cell_id];
671 }
672 }
673 } /* With anisotropic porosity */
674 else if (porosi != NULL && porosf != NULL) {
675 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
676 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
677 cs_lnum_t cell_id = b_face_cells[face_id];
678 coefaq[face_id][0] = porosf[cell_id][0]*coefav[face_id][0]
679 + porosf[cell_id][3]*coefav[face_id][1]
680 + porosf[cell_id][5]*coefav[face_id][2];
681 coefaq[face_id][1] = porosf[cell_id][3]*coefav[face_id][0]
682 + porosf[cell_id][1]*coefav[face_id][1]
683 + porosf[cell_id][4]*coefav[face_id][2];
684 coefaq[face_id][2] = porosf[cell_id][5]*coefav[face_id][0]
685 + porosf[cell_id][4]*coefav[face_id][1]
686 + porosf[cell_id][2]*coefav[face_id][2];
687 f_momentum[face_id][0] = ( porosf[cell_id][0]*vel[cell_id][0]
688 + porosf[cell_id][3]*vel[cell_id][1]
689 + porosf[cell_id][5]*vel[cell_id][2] );
690 f_momentum[face_id][1] = ( porosf[cell_id][3]*vel[cell_id][0]
691 + porosf[cell_id][1]*vel[cell_id][1]
692 + porosf[cell_id][4]*vel[cell_id][2] );
693 f_momentum[face_id][2] = ( porosf[cell_id][5]*vel[cell_id][0]
694 + porosf[cell_id][4]*vel[cell_id][1]
695 + porosf[cell_id][2]*vel[cell_id][2] );
696 }
697 }
698
699 }
700
701 /*==========================================================================
702 2. Compute mass flux without recontructions
703 ==========================================================================*/
704
705 if (nswrgu <= 1) {
706
707 /* Interior faces */
708
709 for (int g_id = 0; g_id < n_i_groups; g_id++) {
710 # pragma omp parallel for
711 for (int t_id = 0; t_id < n_i_threads; t_id++) {
712 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
713 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
714 face_id++) {
715
716 cs_lnum_t ii = i_face_cells[face_id][0];
717 cs_lnum_t jj = i_face_cells[face_id][1];
718 cs_real_t w_i = weight[face_id] * i_f_face_factor[is_p*face_id][0];
719 cs_real_t w_j = (1. - weight[face_id]) * i_f_face_factor[is_p*face_id][1];
720 /* u, v, w Components */
721 for (int isou = 0; isou < 3; isou++) {
722 i_massflux[face_id] += (w_i * qdm[ii][isou] + w_j * qdm[jj][isou])
723 * i_f_face_normal[face_id][isou];
724 }
725
726 }
727 }
728 }
729
730 /* Boundary faces */
731
732 for (int g_id = 0; g_id < n_b_groups; g_id++) {
733 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
734 for (int t_id = 0; t_id < n_b_threads; t_id++) {
735 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
736 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
737 face_id++) {
738
739 /* u, v, w Components */
740 for (int isou = 0; isou < 3; isou++) {
741 double pfac = inc*coefaq[face_id][isou];
742
743 /* coefbv is a matrix */
744 for (int jsou = 0; jsou < 3; jsou++) {
745 pfac += coefbv[face_id][jsou][isou]*f_momentum[face_id][jsou];
746 }
747 pfac *= b_f_face_factor[is_p*face_id];
748
749 b_massflux[face_id] += pfac*b_f_face_normal[face_id][isou];
750 }
751
752 }
753 }
754 }
755
756 }
757
758 /*==========================================================================
759 4. Compute mass flux with reconstruction method if the mesh is
760 non orthogonal
761 ==========================================================================*/
762
763 if (nswrgu > 1) {
764
765 BFT_MALLOC(grdqdm, n_cells_ext, cs_real_33_t);
766
767 /* Computation of momentum gradient
768 (vectorial gradient, the periodicity has already been treated) */
769
770 cs_gradient_vector(var_name,
771 gradient_type,
772 halo_type,
773 inc,
774 nswrgu,
775 iwarnu,
776 imligu,
777 epsrgu,
778 climgu,
779 (const cs_real_3_t*)coefaq,
780 coefbv,
781 qdm,
782 NULL, /* weighted gradient */
783 NULL, /* cpl */
784 grdqdm);
785
786 /* Mass flow through interior faces */
787
788 for (int g_id = 0; g_id < n_i_groups; g_id++) {
789 # pragma omp parallel for
790 for (int t_id = 0; t_id < n_i_threads; t_id++) {
791 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
792 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
793 face_id++) {
794
795 cs_lnum_t ii = i_face_cells[face_id][0];
796 cs_lnum_t jj = i_face_cells[face_id][1];
797
798 double dofx = dofij[face_id][0];
799 double dofy = dofij[face_id][1];
800 double dofz = dofij[face_id][2];
801
802 cs_real_t w_i = weight[face_id] * i_f_face_factor[is_p*face_id][0];
803 cs_real_t w_j = (1. - weight[face_id]) * i_f_face_factor[is_p*face_id][1];
804
805 /* Terms along U, V, W */
806 for (int isou = 0; isou < 3; isou++) {
807
808 i_massflux[face_id] = i_massflux[face_id]
809 /* Non-reconstructed term */
810 + (w_i * qdm[ii][isou] + w_j * qdm[jj][isou]
811
812 /* ---> -> --> ->
813 (Grad(rho U ) . OFij ) . Sij FIXME for discontinuous porous modelling */
814 + 0.5*(grdqdm[ii][isou][0] +grdqdm[jj][isou][0])*dofx
815 + 0.5*(grdqdm[ii][isou][1] +grdqdm[jj][isou][1])*dofy
816 + 0.5*(grdqdm[ii][isou][2] +grdqdm[jj][isou][2])*dofz
817 )*i_f_face_normal[face_id][isou];
818 }
819
820 }
821 }
822
823 }
824
825 /* Mass flow through boundary faces */
826
827 for (int g_id = 0; g_id < n_b_groups; g_id++) {
828 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
829 for (int t_id = 0; t_id < n_b_threads; t_id++) {
830 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
831 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
832 face_id++) {
833
834 cs_lnum_t ii = b_face_cells[face_id];
835 double diipbx = diipb[face_id][0];
836 double diipby = diipb[face_id][1];
837 double diipbz = diipb[face_id][2];
838
839 /* Terms along U, V, W */
840 for (int isou = 0; isou < 3; isou++) {
841
842 double pfac = inc*coefaq[face_id][isou];
843
844 /* coefu is a matrix */
845 for (int jsou = 0; jsou < 3; jsou++) {
846
847 double pip = f_momentum[face_id][jsou]
848 + grdqdm[ii][jsou][0]*diipbx
849 + grdqdm[ii][jsou][1]*diipby
850 + grdqdm[ii][jsou][2]*diipbz;
851
852 pfac += coefbv[face_id][jsou][isou]*pip;
853
854 }
855
856 pfac *= b_f_face_factor[is_p*face_id];
857
858 b_massflux[face_id] += pfac*b_f_face_normal[face_id][isou];
859
860 }
861
862 }
863 }
864 }
865
866
867 /* Deallocation */
868 BFT_FREE(grdqdm);
869
870 }
871
872 BFT_FREE(qdm);
873 BFT_FREE(coefaq);
874 BFT_FREE(f_momentum);
875
876 /*==========================================================================
877 6. Here, we make sure that the mass flux is null at the boundary faces of
878 type symmetry and coupled walls.
879 ==========================================================================*/
880
881 if (iflmb0 == 1) {
882 /* Force flumab to 0 for velocity */
883 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
884 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
885 if (fvq->b_sym_flag[face_id] == 0) {
886 b_massflux[face_id] = 0.;
887 }
888 }
889 }
890
891 }
892
893 /*----------------------------------------------------------------------------*/
894 /*!
895 * \brief Add \f$ \rho \tens{r} \vect{s}_\ij\f$ to a flux.
896 *
897 * \param[in] m pointer to mesh
898 * \param[in] fvq pointer to finite volume quantities
899 * \param[in] f_id field id (or -1)
900 * \param[in] itypfl indicator (take rho into account or not)
901 * - 1 compute \f$ \rho\vect{u}\cdot\vect{s} \f$
902 * - 0 compute \f$ \vect{u}\cdot\vect{s} \f$
903 * \param[in] iflmb0 the mass flux is set to 0 on symmetries if = 1
904 * \param[in] init the mass flux is initialized to 0 if > 0
905 * \param[in] inc indicator
906 * - 0 solve an increment
907 * - 1 otherwise
908 * \param[in] imrgra indicator
909 * - 0 iterative gradient
910 * - 1 least square gradient
911 * \param[in] nswrgu number of sweeps for the reconstruction
912 * of the gradients
913 * \param[in] imligu clipping gradient method
914 * - < 0 no clipping
915 * - = 0 thanks to neighbooring gradients
916 * - = 1 thanks to the mean gradient
917 * \param[in] iwarnu verbosity
918 * \param[in] epsrgu relative precision for the gradient
919 * reconstruction
920 * \param[in] climgu clipping coefficient for the computation of
921 * the gradient
922 * \param[in] c_rho cell density
923 * \param[in] b_rho density at boundary faces
924 * \param[in] c_var variable
925 * \param[in] coefav boundary condition array for the variable
926 * (explicit part - symmetric tensor array)
927 * \param[in] coefbv boundary condition array for the variable
928 * (implicit part - 6x6 symmetric tensor array)
929 * \param[in,out] i_massflux mass flux at interior faces \f$ \dot{m}_\fij \f$
930 * \param[in,out] b_massflux mass flux at boundary faces \f$ \dot{m}_\fib \f$
931 */
932 /*----------------------------------------------------------------------------*/
933
934 void
cs_tensor_face_flux(const cs_mesh_t * m,cs_mesh_quantities_t * fvq,int f_id,int itypfl,int iflmb0,int init,int inc,int imrgra,int nswrgu,int imligu,int iwarnu,double epsrgu,double climgu,const cs_real_t c_rho[],const cs_real_t b_rho[],const cs_real_6_t c_var[],const cs_real_6_t coefav[],const cs_real_66_t coefbv[],cs_real_3_t * restrict i_massflux,cs_real_3_t * restrict b_massflux)935 cs_tensor_face_flux(const cs_mesh_t *m,
936 cs_mesh_quantities_t *fvq,
937 int f_id,
938 int itypfl,
939 int iflmb0,
940 int init,
941 int inc,
942 int imrgra,
943 int nswrgu,
944 int imligu,
945 int iwarnu,
946 double epsrgu,
947 double climgu,
948 const cs_real_t c_rho[],
949 const cs_real_t b_rho[],
950 const cs_real_6_t c_var[],
951 const cs_real_6_t coefav[],
952 const cs_real_66_t coefbv[],
953 cs_real_3_t *restrict i_massflux,
954 cs_real_3_t *restrict b_massflux)
955 {
956 const cs_halo_t *halo = m->halo;
957
958 const cs_lnum_t n_cells = m->n_cells;
959 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
960 const int n_i_groups = m->i_face_numbering->n_groups;
961 const int n_i_threads = m->i_face_numbering->n_threads;
962 const int n_b_groups = m->b_face_numbering->n_groups;
963 const int n_b_threads = m->b_face_numbering->n_threads;
964 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
965 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
966
967 const cs_lnum_2_t *restrict i_face_cells
968 = (const cs_lnum_2_t *restrict)m->i_face_cells;
969 const cs_lnum_t *restrict b_face_cells
970 = (const cs_lnum_t *restrict)m->b_face_cells;
971 const cs_real_t *restrict weight = fvq->weight;
972 const cs_real_3_t *restrict i_f_face_normal
973 = (const cs_real_3_t *restrict)fvq->i_f_face_normal;
974 const cs_real_3_t *restrict b_f_face_normal
975 = (const cs_real_3_t *restrict)fvq->b_f_face_normal;
976 const cs_real_3_t *restrict diipb
977 = (const cs_real_3_t *restrict)fvq->diipb;
978 const cs_real_3_t *restrict dofij
979 = (const cs_real_3_t *restrict)fvq->dofij;
980
981 /* Local variables */
982
983 char var_name[64];
984
985 cs_real_6_t *c_mass_var, *b_mass_var, *coefaq;
986
987 cs_field_t *f;
988
989 BFT_MALLOC(c_mass_var, n_cells_ext, cs_real_6_t);
990 BFT_MALLOC(b_mass_var, m->n_b_faces, cs_real_6_t);
991 BFT_MALLOC(coefaq, m->n_b_faces, cs_real_6_t);
992
993 /*==========================================================================
994 1. Initialization
995 ==========================================================================*/
996
997 /* Choose gradient type */
998
999 cs_halo_type_t halo_type = CS_HALO_STANDARD;
1000 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
1001
1002 cs_gradient_type_by_imrgra(imrgra,
1003 &gradient_type,
1004 &halo_type);
1005
1006 if (f_id != -1) {
1007 f = cs_field_by_id(f_id);
1008 snprintf(var_name, 63, "%s", f->name);
1009 }
1010 else {
1011 strncpy(var_name, "[tensor face flux]", 63);
1012 }
1013 var_name[63] = '\0';
1014
1015 /* ---> Momentum computation */
1016
1017 if (init == 1) {
1018 # pragma omp parallel for
1019 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
1020 for (int i = 0; i < 3; i++)
1021 i_massflux[face_id][i] = 0.;
1022 }
1023 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1024 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1025 for (int i = 0; i < 3; i++)
1026 b_massflux[face_id][i] = 0.;
1027 }
1028
1029 } else if (init != 0) {
1030 bft_error(__FILE__, __LINE__, 0,
1031 _("invalid value of init"));
1032 }
1033
1034 /* Porosity fields */
1035 cs_field_t *fporo = cs_field_by_name_try("porosity");
1036 cs_field_t *ftporo = cs_field_by_name_try("tensorial_porosity");
1037
1038 cs_real_t *porosi = NULL;
1039 cs_real_6_t *porosf = NULL;
1040
1041 if (cs_glob_porous_model == 1 || cs_glob_porous_model == 2) {
1042 porosi = fporo->val;
1043 if (ftporo != NULL) {
1044 porosf = (cs_real_6_t *)ftporo->val;
1045 }
1046 }
1047
1048 /* Standard mass flux */
1049 if (itypfl == 1) {
1050
1051 /* Without porosity */
1052 if (porosi == NULL) {
1053 # pragma omp parallel for
1054 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1055 for (int isou = 0; isou < 6; isou++) {
1056 c_mass_var[cell_id][isou] = c_rho[cell_id]*c_var[cell_id][isou];
1057 }
1058 }
1059 /* With porosity */
1060 } else if (porosi != NULL && porosf == NULL) {
1061 # pragma omp parallel for
1062 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1063 for (int isou = 0; isou < 6; isou++) {
1064 c_mass_var[cell_id][isou] = c_rho[cell_id]*c_var[cell_id][isou]*porosi[cell_id];
1065 }
1066 }
1067 /* With anisotropic porosity */
1068 } else if (porosi != NULL && porosf != NULL) {
1069 # pragma omp parallel for
1070 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1071 cs_math_sym_33_product(porosf[cell_id],
1072 c_var[cell_id],
1073 c_mass_var[cell_id]);
1074
1075 for (int isou = 0; isou < 6; isou++)
1076 c_mass_var[cell_id][isou] *= c_rho[cell_id];
1077 }
1078 }
1079
1080 /* Velocity flux */
1081 } else {
1082
1083 /* Without porosity */
1084 if (porosi == NULL) {
1085 # pragma omp parallel for
1086 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1087 for (int isou = 0; isou < 6; isou++) {
1088 c_mass_var[cell_id][isou] = c_var[cell_id][isou];
1089 }
1090 }
1091 /* With porosity */
1092 } else if (porosi != NULL && porosf == NULL) {
1093 # pragma omp parallel for
1094 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1095 for (int isou = 0; isou < 6; isou++) {
1096 c_mass_var[cell_id][isou] = c_var[cell_id][isou]*porosi[cell_id];
1097 }
1098 }
1099 /* With anisotropic porosity */
1100 } else if (porosi != NULL && porosf != NULL) {
1101 # pragma omp parallel for
1102 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1103 cs_math_sym_33_product(porosf[cell_id],
1104 c_var[cell_id],
1105 c_mass_var[cell_id]);
1106 }
1107 }
1108 }
1109
1110 /* ---> Periodicity and parallelism treatment */
1111
1112 if (halo != NULL) {
1113 cs_halo_sync_var_strided(halo, halo_type, (cs_real_t *)c_mass_var, 6);
1114 if (cs_glob_mesh->n_init_perio > 0)
1115 cs_halo_perio_sync_var_sym_tens(halo, halo_type, (cs_real_t *)c_mass_var);
1116 }
1117
1118 /* Standard mass flux */
1119 if (itypfl == 1) {
1120
1121 /* Without porosity */
1122 if (porosi == NULL) {
1123 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1124 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1125 cs_lnum_t cell_id = b_face_cells[face_id];
1126 for (int isou = 0; isou < 6; isou++) {
1127 coefaq[face_id][isou] = b_rho[face_id]*coefav[face_id][isou];
1128 b_mass_var[face_id][isou] = b_rho[face_id]*c_var[cell_id][isou];
1129 }
1130 }
1131 /* With porosity */
1132 } else if (porosi != NULL && porosf == NULL) {
1133 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1134 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1135 cs_lnum_t cell_id = b_face_cells[face_id];
1136 for (int isou = 0; isou < 6; isou++) {
1137 coefaq[face_id][isou] = b_rho[face_id]
1138 *coefav[face_id][isou]*porosi[cell_id];
1139 b_mass_var[face_id][isou] = b_rho[face_id]*c_var[cell_id][isou]*porosi[cell_id];
1140 }
1141 }
1142 /* With anisotropic porosity */
1143 } else if (porosi != NULL && porosf != NULL) {
1144 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1145 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1146 cs_lnum_t cell_id = b_face_cells[face_id];
1147
1148 cs_math_sym_33_product(porosf[cell_id],
1149 coefav[face_id],
1150 coefaq[face_id]);
1151
1152 for (int isou = 0; isou < 6; isou++)
1153 coefaq[face_id][isou] *= b_rho[face_id];
1154
1155 cs_math_sym_33_product(porosf[cell_id],
1156 c_var[cell_id],
1157 b_mass_var[face_id]);
1158
1159 for (int isou = 0; isou < 6; isou++)
1160 b_mass_var[face_id][isou] *= b_rho[face_id];
1161
1162 }
1163 }
1164
1165 /* Velocity flux */
1166 } else {
1167
1168 /* Without porosity */
1169 if (porosi == NULL) {
1170 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1171 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1172 cs_lnum_t cell_id = b_face_cells[face_id];
1173 for (int isou = 0; isou < 6; isou++) {
1174 coefaq[face_id][isou] = coefav[face_id][isou];
1175 b_mass_var[face_id][isou] = c_var[cell_id][isou];
1176 }
1177 }
1178 /* With porosity */
1179 } else if (porosi != NULL && porosf == NULL) {
1180 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1181 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1182 cs_lnum_t cell_id = b_face_cells[face_id];
1183 for (int isou = 0; isou < 6; isou++) {
1184 coefaq[face_id][isou] = coefav[face_id][isou]*porosi[cell_id];
1185 b_mass_var[face_id][isou] = c_var[cell_id][isou]*porosi[cell_id];
1186 }
1187 }
1188 /* With anisotropic porosity */
1189 } else if (porosi != NULL && porosf != NULL) {
1190 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1191 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1192 cs_lnum_t cell_id = b_face_cells[face_id];
1193
1194 cs_math_sym_33_product(porosf[cell_id],
1195 coefav[face_id],
1196 coefaq[face_id]);
1197
1198 cs_math_sym_33_product(porosf[cell_id],
1199 c_var[cell_id],
1200 b_mass_var[face_id]);
1201 }
1202 }
1203
1204 }
1205
1206 /*==========================================================================
1207 2. Compute mass flux without recontructions
1208 ==========================================================================*/
1209
1210 if (nswrgu <= 1) {
1211
1212 /* Interior faces */
1213
1214 for (int g_id = 0; g_id < n_i_groups; g_id++) {
1215 # pragma omp parallel for
1216 for (int t_id = 0; t_id < n_i_threads; t_id++) {
1217 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
1218 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
1219 face_id++) {
1220
1221 cs_lnum_t ii = i_face_cells[face_id][0];
1222 cs_lnum_t jj = i_face_cells[face_id][1];
1223
1224 cs_real_t w_i = weight[face_id];
1225 cs_real_t w_j = (1. - weight[face_id]);
1226
1227 cs_real_6_t f_mass_var;
1228
1229 for (int isou = 0; isou < 6; isou++)
1230 f_mass_var[isou] = w_i * c_mass_var[ii][isou] + w_j * c_mass_var[jj][isou];
1231
1232 cs_math_sym_33_3_product_add(f_mass_var,
1233 i_f_face_normal[face_id],
1234 i_massflux[face_id]);
1235
1236 }
1237 }
1238 }
1239
1240 /* Boundary faces */
1241
1242 for (int g_id = 0; g_id < n_b_groups; g_id++) {
1243 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1244 for (int t_id = 0; t_id < n_b_threads; t_id++) {
1245 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
1246 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
1247 face_id++) {
1248
1249 cs_real_6_t f_mass_var;
1250
1251 /* var_f = a + b * var_i */
1252 for (int isou = 0; isou < 6; isou++)
1253 f_mass_var[isou] = inc*coefaq[face_id][isou];
1254
1255 cs_math_66_6_product_add(coefbv[face_id],
1256 b_mass_var[face_id],
1257 f_mass_var);
1258
1259 cs_math_sym_33_3_product_add(f_mass_var,
1260 b_f_face_normal[face_id],
1261 b_massflux[face_id]);
1262
1263 }
1264 }
1265 }
1266
1267 }
1268
1269 /*==========================================================================
1270 4. Compute mass flux with reconstruction technics if the mesh is
1271 non orthogonal
1272 ==========================================================================*/
1273
1274 if (nswrgu > 1) {
1275
1276 cs_real_63_t *c_grad_mvar;
1277 BFT_MALLOC(c_grad_mvar, n_cells_ext, cs_real_63_t);
1278
1279 /* Computation of c_mass_var gradient
1280 (tensor gradient, the periodicity has already been treated) */
1281
1282 cs_gradient_tensor_synced_input(var_name,
1283 gradient_type,
1284 halo_type,
1285 inc,
1286 nswrgu,
1287 iwarnu,
1288 imligu,
1289 epsrgu,
1290 climgu,
1291 (const cs_real_6_t *)coefaq,
1292 (const cs_real_66_t *)coefbv,
1293 (const cs_real_6_t *)c_mass_var,
1294 c_grad_mvar);
1295
1296 /* Mass flow through interior faces */
1297
1298 for (int g_id = 0; g_id < n_i_groups; g_id++) {
1299 # pragma omp parallel for
1300 for (int t_id = 0; t_id < n_i_threads; t_id++) {
1301 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
1302 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
1303 face_id++) {
1304
1305 cs_lnum_t ii = i_face_cells[face_id][0];
1306 cs_lnum_t jj = i_face_cells[face_id][1];
1307
1308 cs_real_t w_i = weight[face_id];
1309 cs_real_t w_j = (1. - weight[face_id]);
1310
1311 cs_real_t f_mass_var[6];
1312
1313 for (cs_lnum_t isou = 0; isou < 6; isou++)
1314 /* Non-reconstructed face value */
1315 f_mass_var[isou] = w_i * c_mass_var[ii][isou] + w_j * c_mass_var[jj][isou]
1316 /* Reconstruction: face gradient times OF */
1317 + 0.5*(c_grad_mvar[ii][isou][0] +c_grad_mvar[jj][isou][0])*dofij[face_id][0]
1318 + 0.5*(c_grad_mvar[ii][isou][1] +c_grad_mvar[jj][isou][1])*dofij[face_id][1]
1319 + 0.5*(c_grad_mvar[ii][isou][2] +c_grad_mvar[jj][isou][2])*dofij[face_id][2];
1320
1321 cs_math_sym_33_3_product_add(f_mass_var,
1322 i_f_face_normal[face_id],
1323 i_massflux[face_id]);
1324
1325 }
1326 }
1327
1328 }
1329
1330 /* Mass flow through boundary faces */
1331
1332 for (int g_id = 0; g_id < n_b_groups; g_id++) {
1333 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1334 for (int t_id = 0; t_id < n_b_threads; t_id++) {
1335 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
1336 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
1337 face_id++) {
1338
1339 cs_lnum_t ii = b_face_cells[face_id];
1340
1341 cs_real_6_t f_mass_var;
1342
1343 /* var_f = a + b * var_I' */
1344 for (int isou = 0; isou < 6; isou++)
1345 f_mass_var[isou] = inc*coefaq[face_id][isou];
1346
1347
1348 /* Add the reconstruction to get value in I' */
1349 for (int jsou = 0; jsou < 6; jsou++)
1350 b_mass_var[face_id][jsou] += c_grad_mvar[ii][jsou][0]*diipb[face_id][0]
1351 + c_grad_mvar[ii][jsou][1]*diipb[face_id][1]
1352 + c_grad_mvar[ii][jsou][2]*diipb[face_id][2];
1353
1354 cs_math_66_6_product_add(coefbv[face_id],
1355 b_mass_var[face_id],
1356 f_mass_var);
1357
1358 cs_math_sym_33_3_product_add(f_mass_var,
1359 b_f_face_normal[face_id],
1360 b_massflux[face_id]);
1361
1362 }
1363 }
1364 }
1365
1366
1367 /* Deallocation */
1368 BFT_FREE(c_grad_mvar);
1369
1370 }
1371
1372 BFT_FREE(c_mass_var);
1373 BFT_FREE(coefaq);
1374 BFT_FREE(b_mass_var);
1375
1376 /*==========================================================================
1377 6. Here, we make sure that the mass flux is null at the boundary faces of
1378 type symmetry and coupled walls.
1379 ==========================================================================*/
1380
1381 if (iflmb0 == 1) {
1382 /* Force flumab to 0 for velocity */
1383 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1384 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1385 if (fvq->b_sym_flag[face_id] == 0) {
1386 for (int isou = 0; isou < 3; isou++)
1387 b_massflux[face_id][isou] = 0.;
1388 }
1389 }
1390 }
1391
1392 }
1393
1394 /*----------------------------------------------------------------------------*/
1395 /*!
1396 * \brief Add the integrated mass flux on the cells.
1397 *
1398 * \f[
1399 * \dot{m}_i = \dot{m}_i + \sum_{\fij \in \Facei{\celli}} \dot{m}_\ij
1400 * \f]
1401 *
1402 * \param[in] m pointer to mesh
1403 * \param[in] init indicator
1404 * - 1 initialize the divergence to 0
1405 * - 0 otherwise
1406 * \param[in] i_massflux mass flux at interior faces
1407 * \param[in] b_massflux mass flux at boundary faces
1408 * \param[in,out] diverg mass flux divergence
1409 */
1410 /*----------------------------------------------------------------------------*/
1411
1412 void
cs_divergence(const cs_mesh_t * m,int init,const cs_real_t i_massflux[],const cs_real_t b_massflux[],cs_real_t * restrict diverg)1413 cs_divergence(const cs_mesh_t *m,
1414 int init,
1415 const cs_real_t i_massflux[],
1416 const cs_real_t b_massflux[],
1417 cs_real_t *restrict diverg)
1418 {
1419 const cs_lnum_t n_cells = m->n_cells;
1420 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1421 const int n_i_groups = m->i_face_numbering->n_groups;
1422 const int n_i_threads = m->i_face_numbering->n_threads;
1423 const int n_b_groups = m->b_face_numbering->n_groups;
1424 const int n_b_threads = m->b_face_numbering->n_threads;
1425 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
1426 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
1427
1428 const cs_lnum_2_t *restrict i_face_cells
1429 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1430 const cs_lnum_t *restrict b_face_cells
1431 = (const cs_lnum_t *restrict)m->b_face_cells;
1432
1433 /*==========================================================================
1434 1. Initialization
1435 ==========================================================================*/
1436
1437 if (init >= 1) {
1438
1439 # pragma omp parallel for
1440 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
1441 diverg[cell_id] = 0.;
1442
1443 }
1444 else if (init == 0 && n_cells_ext > n_cells) {
1445
1446 # pragma omp parallel for if(n_cells_ext - n_cells > CS_THR_MIN)
1447 for (cs_lnum_t cell_id = n_cells+0; cell_id < n_cells_ext; cell_id++)
1448 diverg[cell_id] = 0.;
1449
1450 }
1451 else if (init != 0)
1452 bft_error(__FILE__, __LINE__, 0, _("invalid value of init"));
1453
1454
1455 /*==========================================================================
1456 2. Integration on internal faces
1457 ==========================================================================*/
1458
1459 for (int g_id = 0; g_id < n_i_groups; g_id++) {
1460
1461 # pragma omp parallel for
1462 for (int t_id = 0; t_id < n_i_threads; t_id++) {
1463 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
1464 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
1465 face_id++) {
1466
1467 cs_lnum_t ii = i_face_cells[face_id][0];
1468 cs_lnum_t jj = i_face_cells[face_id][1];
1469
1470 diverg[ii] += i_massflux[face_id];
1471 diverg[jj] -= i_massflux[face_id];
1472
1473 }
1474 }
1475
1476 } /* Loop on openMP groups */
1477
1478
1479 /*==========================================================================
1480 3. Integration on border faces
1481 ==========================================================================*/
1482
1483 for (int g_id = 0; g_id < n_b_groups; g_id++) {
1484
1485 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1486 for (int t_id = 0; t_id < n_b_threads; t_id++) {
1487 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
1488 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
1489 face_id++) {
1490
1491 cs_lnum_t ii = b_face_cells[face_id];
1492 diverg[ii] += b_massflux[face_id];
1493
1494 }
1495 }
1496
1497 } /* Loop on openMP groups */
1498
1499 }
1500
1501 /*----------------------------------------------------------------------------*/
1502 /*!
1503 * \brief Add the integrated mass flux on the cells for a tensor variable.
1504 *
1505 * \f[
1506 * \dot{m}_i = \dot{m}_i + \sum_{\fij \in \Facei{\celli}} \dot{m}_\ij
1507 * \f]
1508 *
1509 * \param[in] m pointer to mesh
1510 * \param[in] init indicator
1511 * - 1 initialize the divergence to 0
1512 * - 0 otherwise
1513 * \param[in] i_massflux mass flux vector at interior faces
1514 * \param[in] b_massflux mass flux vector at boundary faces
1515 * \param[in,out] diverg mass flux divergence vector
1516 */
1517 /*----------------------------------------------------------------------------*/
1518
1519 void
cs_tensor_divergence(const cs_mesh_t * m,int init,const cs_real_3_t i_massflux[],const cs_real_3_t b_massflux[],cs_real_3_t * restrict diverg)1520 cs_tensor_divergence(const cs_mesh_t *m,
1521 int init,
1522 const cs_real_3_t i_massflux[],
1523 const cs_real_3_t b_massflux[],
1524 cs_real_3_t *restrict diverg)
1525 {
1526 const cs_lnum_t n_cells = m->n_cells;
1527 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1528 const int n_i_groups = m->i_face_numbering->n_groups;
1529 const int n_i_threads = m->i_face_numbering->n_threads;
1530 const int n_b_groups = m->b_face_numbering->n_groups;
1531 const int n_b_threads = m->b_face_numbering->n_threads;
1532 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
1533 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
1534
1535 const cs_lnum_2_t *restrict i_face_cells
1536 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1537 const cs_lnum_t *restrict b_face_cells
1538 = (const cs_lnum_t *restrict)m->b_face_cells;
1539
1540 /*==========================================================================
1541 1. Initialization
1542 ==========================================================================*/
1543
1544 if (init >= 1) {
1545 # pragma omp parallel for
1546 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++) {
1547 for (int isou = 0; isou < 3; isou++) {
1548 diverg[cell_id][isou] = 0.;
1549 }
1550 }
1551 } else if (init == 0 && n_cells_ext > n_cells) {
1552 # pragma omp parallel for if(n_cells_ext - n_cells > CS_THR_MIN)
1553 for (cs_lnum_t cell_id = n_cells+0; cell_id < n_cells_ext; cell_id++) {
1554 for (int isou = 0; isou < 3; isou++) {
1555 diverg[cell_id][isou] = 0.;
1556 }
1557 }
1558 } else if (init != 0) {
1559 bft_error(__FILE__, __LINE__, 0,
1560 _("invalid value of init"));
1561 }
1562
1563
1564 /*==========================================================================
1565 2. Integration on internal faces
1566 ==========================================================================*/
1567
1568 for (int g_id = 0; g_id < n_i_groups; g_id++) {
1569 # pragma omp parallel for
1570 for (int t_id = 0; t_id < n_i_threads; t_id++) {
1571 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
1572 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
1573 face_id++) {
1574
1575 cs_lnum_t ii = i_face_cells[face_id][0];
1576 cs_lnum_t jj = i_face_cells[face_id][1];
1577
1578 for (int isou = 0; isou < 3; isou++) {
1579 diverg[ii][isou] += i_massflux[face_id][isou];
1580 diverg[jj][isou] -= i_massflux[face_id][isou];
1581 }
1582
1583 }
1584 }
1585 }
1586
1587
1588 /*==========================================================================
1589 3. Integration on border faces
1590 ==========================================================================*/
1591
1592 for (int g_id = 0; g_id < n_b_groups; g_id++) {
1593 # pragma omp parallel for if(m->n_b_faces > CS_THR_MIN)
1594 for (int t_id = 0; t_id < n_b_threads; t_id++) {
1595 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
1596 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
1597 face_id++) {
1598
1599 cs_lnum_t ii = b_face_cells[face_id];
1600 for (int isou = 0; isou < 3; isou++) {
1601 diverg[ii][isou] += b_massflux[face_id][isou];
1602 }
1603
1604 }
1605 }
1606 }
1607
1608 }
1609
1610 /*----------------------------------------------------------------------------*/
1611 /*!
1612 * \brief Project the external source terms to the faces in coherence with
1613 * cs_face_diffusion_scalar for the improved hydrostatic pressure algorithm
1614 * (iphydr=1).
1615 *
1616 * \param[in] m pointer to mesh
1617 * \param[in] fvq pointer to finite volume quantities
1618 * \param[in] init indicator
1619 * - 1 initialize the mass flux to 0
1620 * - 0 otherwise
1621 * \param[in] nswrgu number of reconstruction sweeps for the
1622 * gradients
1623 * \param[in] frcxt body force creating the hydrostatic pressure
1624 * \param[in] cofbfp boundary condition array for the diffusion
1625 * of the variable (implicit part)
1626 * \param[in,out] i_massflux mass flux at interior faces
1627 * \param[in,out] b_massflux mass flux at boundary faces
1628 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1629 * at interior faces for the r.h.s.
1630 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1631 * at border faces for the r.h.s.
1632 * \param[in] viselx viscosity by cell, dir x
1633 * \param[in] visely viscosity by cell, dir y
1634 * \param[in] viselz viscosity by cell, dir z
1635 */
1636 /*----------------------------------------------------------------------------*/
1637
1638 void
cs_ext_force_flux(const cs_mesh_t * m,cs_mesh_quantities_t * fvq,int init,int nswrgu,const cs_real_3_t frcxt[],const cs_real_t cofbfp[],cs_real_t * restrict i_massflux,cs_real_t * restrict b_massflux,const cs_real_t i_visc[],const cs_real_t b_visc[],const cs_real_t viselx[],const cs_real_t visely[],const cs_real_t viselz[])1639 cs_ext_force_flux(const cs_mesh_t *m,
1640 cs_mesh_quantities_t *fvq,
1641 int init,
1642 int nswrgu,
1643 const cs_real_3_t frcxt[],
1644 const cs_real_t cofbfp[],
1645 cs_real_t *restrict i_massflux,
1646 cs_real_t *restrict b_massflux,
1647 const cs_real_t i_visc[],
1648 const cs_real_t b_visc[],
1649 const cs_real_t viselx[],
1650 const cs_real_t visely[],
1651 const cs_real_t viselz[])
1652 {
1653 const cs_lnum_2_t *restrict i_face_cells
1654 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1655 const cs_lnum_t *restrict b_face_cells
1656 = (const cs_lnum_t *restrict)m->b_face_cells;
1657 const cs_real_t *restrict i_dist = fvq->i_dist;
1658 const cs_real_t *restrict b_dist = fvq->b_dist;
1659 const cs_real_t *restrict i_f_face_surf = fvq->i_f_face_surf;
1660 const cs_real_3_t *restrict cell_cen
1661 = (const cs_real_3_t *restrict)fvq->cell_cen;
1662 const cs_real_3_t *restrict b_face_normal
1663 = (const cs_real_3_t *restrict)fvq->b_face_normal;
1664 const cs_real_3_t *restrict i_face_cog
1665 = (const cs_real_3_t *restrict)fvq->i_face_cog;
1666 const cs_real_3_t *restrict diipf
1667 = (const cs_real_3_t *restrict)fvq->diipf;
1668 const cs_real_3_t *restrict djjpf
1669 = (const cs_real_3_t *restrict)fvq->djjpf;
1670
1671 /*Additional terms due to porosity */
1672
1673 cs_field_t *f_i_poro_duq_0 = cs_field_by_name_try("i_poro_duq_0");
1674
1675 cs_real_t *i_poro_duq_0;
1676 cs_real_t *i_poro_duq_1;
1677 cs_real_t *b_poro_duq;
1678 cs_real_t _f_ext = 0.;
1679
1680 int is_p = 0; /* Is porous ? */
1681 if (f_i_poro_duq_0 != NULL) {
1682
1683 is_p = 1;
1684 i_poro_duq_0 = f_i_poro_duq_0->val;
1685 i_poro_duq_1 = cs_field_by_name("i_poro_duq_1")->val;
1686 b_poro_duq = cs_field_by_name("b_poro_duq")->val;
1687
1688 }
1689 else {
1690
1691 i_poro_duq_0 = &_f_ext;
1692 i_poro_duq_1 = &_f_ext;
1693 b_poro_duq = &_f_ext;
1694
1695 }
1696
1697 /*==========================================================================*/
1698
1699 /*==========================================================================
1700 1. Initialization
1701 ==========================================================================*/
1702
1703 if (init == 1) {
1704
1705 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++)
1706 i_massflux[face_id] = 0.;
1707 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++)
1708 b_massflux[face_id] = 0.;
1709
1710 }
1711 else if (init != 0)
1712 bft_error(__FILE__, __LINE__, 0, _("invalid value of init"));
1713
1714 /*==========================================================================
1715 2. Update mass flux without reconstruction technics
1716 ==========================================================================*/
1717
1718 if (nswrgu <= 1) {
1719
1720 /* Mass flow through interior faces */
1721
1722 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
1723
1724 cs_lnum_t ii = i_face_cells[face_id][0];
1725 cs_lnum_t jj = i_face_cells[face_id][1];
1726
1727 cs_real_2_t poro = { i_poro_duq_0[is_p*face_id],
1728 i_poro_duq_1[is_p*face_id] };
1729
1730 i_massflux[face_id] +=
1731 i_visc[face_id]*( (i_face_cog[face_id][0]-cell_cen[ii][0])*frcxt[ii][0]
1732 + (i_face_cog[face_id][1]-cell_cen[ii][1])*frcxt[ii][1]
1733 + (i_face_cog[face_id][2]-cell_cen[ii][2])*frcxt[ii][2]
1734 + poro[0]
1735 - (i_face_cog[face_id][0]-cell_cen[jj][0])*frcxt[jj][0]
1736 - (i_face_cog[face_id][1]-cell_cen[jj][1])*frcxt[jj][1]
1737 - (i_face_cog[face_id][2]-cell_cen[jj][2])*frcxt[jj][2]
1738 - poro[1]
1739 );
1740
1741 }
1742
1743 /* Mass flux through boundary faces */
1744
1745 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1746
1747 cs_lnum_t ii = b_face_cells[face_id];
1748
1749 /* To avoid division by 0, no division by the fluid surface */
1750
1751 cs_real_3_t normal;
1752 cs_math_3_normalise(b_face_normal[face_id], normal);
1753
1754 cs_real_t poro = b_poro_duq[is_p*face_id];
1755
1756 b_massflux[face_id] += b_visc[face_id] * cofbfp[face_id] *
1757 ( cs_math_3_dot_product(frcxt[ii], normal) * b_dist[face_id] + poro );
1758
1759 }
1760
1761 /*==========================================================================
1762 3. Update mass flux with reconstruction technics
1763 ==========================================================================*/
1764
1765 }
1766 else {
1767
1768 /* Mass flux through interior faces */
1769
1770 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
1771
1772 cs_lnum_t ii = i_face_cells[face_id][0];
1773 cs_lnum_t jj = i_face_cells[face_id][1];
1774
1775 cs_real_2_t poro = { i_poro_duq_0[is_p*face_id],
1776 i_poro_duq_1[is_p*face_id] };
1777
1778 double surfn = i_f_face_surf[face_id];
1779
1780 i_massflux[face_id] += i_visc[face_id]*
1781 ( ( i_face_cog[face_id][0] - cell_cen[ii][0] ) * frcxt[ii][0] +
1782 ( i_face_cog[face_id][1] - cell_cen[ii][1] ) * frcxt[ii][1] +
1783 ( i_face_cog[face_id][2] - cell_cen[ii][2] ) * frcxt[ii][2] +
1784 poro[0] -
1785 ( i_face_cog[face_id][0] - cell_cen[jj][0] ) * frcxt[jj][0] -
1786 ( i_face_cog[face_id][1] - cell_cen[jj][1] ) * frcxt[jj][1] -
1787 ( i_face_cog[face_id][2] - cell_cen[jj][2] ) * frcxt[jj][2] -
1788 poro[1]
1789 )
1790 + surfn/i_dist[face_id]*0.5
1791 *( (djjpf[face_id][0] - diipf[face_id][0])*
1792 (viselx[ii]*frcxt[ii][0] + viselx[jj]*frcxt[jj][0])
1793 + (djjpf[face_id][1] - diipf[face_id][1])*
1794 (visely[ii]*frcxt[ii][1] + visely[jj]*frcxt[jj][1])
1795 + (djjpf[face_id][2] - diipf[face_id][2])*
1796 (viselz[ii]*frcxt[ii][2] + viselz[jj]*frcxt[jj][2])
1797 );
1798
1799 } /* Loop on interior faces */
1800
1801 /* Mass flux through boundary faces */
1802
1803 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1804
1805 cs_lnum_t ii = b_face_cells[face_id];
1806
1807 /* To avoid division by 0, no division by the fluid surface */
1808 cs_real_3_t normal;
1809 cs_math_3_normalise(b_face_normal[face_id], normal);
1810
1811 cs_real_t poro = b_poro_duq[is_p*face_id];
1812
1813 b_massflux[face_id] += b_visc[face_id] * cofbfp[face_id] *
1814 (cs_math_3_dot_product(frcxt[ii], normal) * b_dist[face_id] + poro);
1815
1816 } /* Loop on border faces */
1817
1818 }
1819
1820 }
1821
1822 /*----------------------------------------------------------------------------*/
1823 /*!
1824 * \brief Project the external source terms to the faces in coherence with
1825 * cs_face_anisotropic_diffusion_scalar for the improved hydrostatic pressure
1826 * algorithm (iphydr=1).
1827 *
1828 * \param[in] m pointer to mesh
1829 * \param[in] fvq pointer to finite volume quantities
1830 * \param[in] init indicator
1831 * - 1 initialize the mass flux to 0
1832 * - 0 otherwise
1833 * \param[in] nswrgp number of reconstruction sweeps for the
1834 * gradients
1835 * \param[in] ircflp indicator
1836 * - 1 flux reconstruction,
1837 * - 0 otherwise
1838 * \param[in] frcxt body force creating the hydrostatic pressure
1839 * \param[in] cofbfp boundary condition array for the diffusion
1840 * of the variable (implicit part)
1841 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1842 * at interior faces for the r.h.s.
1843 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1844 * at border faces for the r.h.s.
1845 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
1846 * \param[in] weighf internal face weight between cells i j in case
1847 * of tensor diffusion
1848 * \param[in,out] i_massflux mass flux at interior faces
1849 * \param[in,out] b_massflux mass flux at boundary faces
1850 */
1851 /*----------------------------------------------------------------------------*/
1852
1853 void
cs_ext_force_anisotropic_flux(const cs_mesh_t * m,cs_mesh_quantities_t * fvq,int init,int nswrgp,int ircflp,const cs_real_3_t frcxt[],const cs_real_t cofbfp[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_6_t viscel[],const cs_real_2_t weighf[],cs_real_t * restrict i_massflux,cs_real_t * restrict b_massflux)1854 cs_ext_force_anisotropic_flux(const cs_mesh_t *m,
1855 cs_mesh_quantities_t *fvq,
1856 int init,
1857 int nswrgp,
1858 int ircflp,
1859 const cs_real_3_t frcxt[],
1860 const cs_real_t cofbfp[],
1861 const cs_real_t i_visc[],
1862 const cs_real_t b_visc[],
1863 cs_real_6_t viscel[],
1864 const cs_real_2_t weighf[],
1865 cs_real_t *restrict i_massflux,
1866 cs_real_t *restrict b_massflux)
1867 {
1868 const cs_halo_t *halo = m->halo;
1869
1870 const cs_lnum_t n_cells = m->n_cells;
1871 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1872
1873 const cs_lnum_2_t *restrict i_face_cells
1874 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1875 const cs_lnum_t *restrict b_face_cells
1876 = (const cs_lnum_t *restrict)m->b_face_cells;
1877 const cs_real_t *restrict b_dist = fvq->b_dist;
1878 const cs_real_3_t *restrict cell_cen
1879 = (const cs_real_3_t *restrict)fvq->cell_cen;
1880 const cs_real_3_t *restrict i_f_face_normal
1881 = (const cs_real_3_t *restrict)fvq->i_f_face_normal;
1882 const cs_real_3_t *restrict b_face_normal
1883 = (const cs_real_3_t *restrict)fvq->b_face_normal;
1884 const cs_real_3_t *restrict i_face_cog
1885 = (const cs_real_3_t *restrict)fvq->i_face_cog;
1886
1887 /* Local variables */
1888
1889 double diippf[3], djjppf[3];
1890 double visci[3][3], viscj[3][3];
1891
1892 /* Porosity fields */
1893 cs_field_t *fporo = cs_field_by_name_try("porosity");
1894 cs_field_t *ftporo = cs_field_by_name_try("tensorial_porosity");
1895
1896 cs_real_t *porosi = NULL;
1897 cs_real_6_t *porosf = NULL;
1898
1899 if (cs_glob_porous_model == 1 || cs_glob_porous_model == 2) {
1900 porosi = fporo->val;
1901 if (ftporo != NULL) {
1902 porosf = (cs_real_6_t *)ftporo->val;
1903 }
1904 }
1905
1906 /*==========================================================================*/
1907
1908 /*==========================================================================
1909 1. Initialization
1910 ==========================================================================*/
1911
1912 if (init == 1) {
1913 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
1914 i_massflux[face_id] = 0.;
1915 }
1916 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1917 b_massflux[face_id] = 0.;
1918 }
1919
1920 } else if (init != 0) {
1921 bft_error(__FILE__, __LINE__, 0,
1922 _("invalid value of init"));
1923 }
1924
1925 /*==========================================================================
1926 2. Update mass flux without reconstruction technics
1927 ==========================================================================*/
1928
1929 if (nswrgp <= 1) {
1930
1931 /* ---> Contribution from interior faces */
1932
1933 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
1934
1935 cs_lnum_t ii = i_face_cells[face_id][0];
1936 cs_lnum_t jj = i_face_cells[face_id][1];
1937
1938 i_massflux[face_id] = i_massflux[face_id]+ i_visc[face_id]*(
1939 ( i_face_cog[face_id][0]
1940 -cell_cen[ii][0])*frcxt[ii][0]
1941 +( i_face_cog[face_id][1]
1942 -cell_cen[ii][1])*frcxt[ii][1]
1943 +( i_face_cog[face_id][2]
1944 -cell_cen[ii][2])*frcxt[ii][2]
1945 -( i_face_cog[face_id][0]
1946 -cell_cen[jj][0])*frcxt[jj][0]
1947 -( i_face_cog[face_id][1]
1948 -cell_cen[jj][1])*frcxt[jj][1]
1949 -( i_face_cog[face_id][2]
1950 -cell_cen[jj][2])*frcxt[jj][2]
1951 );
1952
1953 }
1954
1955 /* ---> Contribution from boundary faces */
1956
1957 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
1958
1959 cs_lnum_t ii = b_face_cells[face_id];
1960
1961 /* To avoid division by 0, no division by the fluid surface */
1962 cs_real_3_t normal;
1963 cs_math_3_normalise(b_face_normal[face_id], normal);
1964
1965
1966 b_massflux[face_id] += b_visc[face_id] * b_dist[face_id]
1967 * cofbfp[face_id]
1968 * cs_math_3_dot_product(frcxt[ii], normal);
1969
1970 }
1971
1972 /*========================================================================
1973 3. Update mass flux with reconstruction technics
1974 ========================================================================*/
1975
1976 } else {
1977
1978 cs_real_6_t *viscce = NULL;
1979 cs_real_6_t *w2 = NULL;
1980
1981 /* Without porosity */
1982 if (porosi == NULL) {
1983 viscce = viscel;
1984
1985 /* With porosity */
1986 } else if (porosi != NULL && porosf == NULL) {
1987 BFT_MALLOC(w2, n_cells_ext, cs_real_6_t);
1988 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1989 for (int isou = 0; isou < 6; isou++) {
1990 w2[cell_id][isou] = porosi[cell_id]*viscel[cell_id][isou];
1991 }
1992 }
1993 viscce = w2;
1994
1995 /* With tensorial porosity */
1996 } else if (porosi != NULL && porosf != NULL) {
1997 BFT_MALLOC(w2, n_cells_ext, cs_real_6_t);
1998 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1999 cs_math_sym_33_product(porosf[cell_id],
2000 viscel[cell_id],
2001 w2[cell_id]);
2002 }
2003 viscce = w2;
2004 }
2005
2006 /* ---> Periodicity and parallelism treatment of symmetric tensors */
2007
2008 if (halo != NULL) {
2009 cs_halo_sync_var_strided(halo, CS_HALO_STANDARD, (cs_real_t *)viscce, 6);
2010
2011 if (m->n_init_perio > 0)
2012 cs_halo_perio_sync_var_sym_tens(halo,
2013 CS_HALO_STANDARD,
2014 (cs_real_t *)viscce);
2015 }
2016
2017 /* ---> Contribution from interior faces */
2018
2019 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
2020
2021 cs_lnum_t ii = i_face_cells[face_id][0];
2022 cs_lnum_t jj = i_face_cells[face_id][1];
2023
2024 /* Recompute II' and JJ' at this level */
2025
2026 visci[0][0] = viscce[ii][0];
2027 visci[1][1] = viscce[ii][1];
2028 visci[2][2] = viscce[ii][2];
2029 visci[1][0] = viscce[ii][3];
2030 visci[0][1] = viscce[ii][3];
2031 visci[2][1] = viscce[ii][4];
2032 visci[1][2] = viscce[ii][4];
2033 visci[2][0] = viscce[ii][5];
2034 visci[0][2] = viscce[ii][5];
2035
2036 /* IF.Ki.S / ||Ki.S||^2 */
2037 double fikdvi = weighf[face_id][0];
2038
2039 /* II" = IF + FI" */
2040 for (int i = 0; i < 3; i++) {
2041 diippf[i] = i_face_cog[face_id][i]-cell_cen[ii][i]
2042 - fikdvi*( visci[0][i]*i_f_face_normal[face_id][0]
2043 + visci[1][i]*i_f_face_normal[face_id][1]
2044 + visci[2][i]*i_f_face_normal[face_id][2] );
2045 }
2046
2047 viscj[0][0] = viscce[jj][0];
2048 viscj[1][1] = viscce[jj][1];
2049 viscj[2][2] = viscce[jj][2];
2050 viscj[1][0] = viscce[jj][3];
2051 viscj[0][1] = viscce[jj][3];
2052 viscj[2][1] = viscce[jj][4];
2053 viscj[1][2] = viscce[jj][4];
2054 viscj[2][0] = viscce[jj][5];
2055 viscj[0][2] = viscce[jj][5];
2056
2057 /* FJ.Kj.S / ||Kj.S||^2 */
2058 double fjkdvi = weighf[face_id][1];
2059
2060 /* JJ" = JF + FJ" */
2061 for (int i = 0; i < 3; i++) {
2062 djjppf[i] = i_face_cog[face_id][i]-cell_cen[jj][i]
2063 + fjkdvi*( viscj[0][i]*i_f_face_normal[face_id][0]
2064 + viscj[1][i]*i_f_face_normal[face_id][1]
2065 + viscj[2][i]*i_f_face_normal[face_id][2] );
2066 }
2067
2068 i_massflux[face_id] = i_massflux[face_id]
2069 + i_visc[face_id]
2070 *( frcxt[ii][0]*( i_face_cog[face_id][0]
2071 -cell_cen[ii][0] )
2072 + frcxt[ii][1]*( i_face_cog[face_id][1]
2073 -cell_cen[ii][1] )
2074 + frcxt[ii][2]*( i_face_cog[face_id][2]
2075 -cell_cen[ii][2] )
2076 - frcxt[jj][0]*( i_face_cog[face_id][0]
2077 -cell_cen[jj][0] )
2078 - frcxt[jj][1]*( i_face_cog[face_id][1]
2079 -cell_cen[jj][1] )
2080 - frcxt[jj][2]*( i_face_cog[face_id][2]
2081 -cell_cen[jj][2] )
2082 )
2083 + i_visc[face_id]*ircflp
2084 *(- frcxt[ii][0]*diippf[0]
2085 - frcxt[ii][1]*diippf[1]
2086 - frcxt[ii][2]*diippf[2]
2087 + frcxt[jj][0]*djjppf[0]
2088 + frcxt[jj][1]*djjppf[1]
2089 + frcxt[jj][2]*djjppf[2]
2090 );
2091
2092 }
2093
2094 /* ---> Contribution from boundary faces */
2095
2096 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
2097
2098 cs_lnum_t ii = b_face_cells[face_id];
2099
2100 /* To avoid division by 0, no division by the fluid surface */
2101 cs_real_3_t normal;
2102 cs_math_3_normalise(b_face_normal[face_id], normal);
2103
2104 /* FIXME: wrong if dirichlet and viscce is really a tensor */
2105 b_massflux[face_id] += b_visc[face_id] * b_dist[face_id]
2106 * cofbfp[face_id]
2107 * cs_math_3_dot_product(frcxt[ii], normal);
2108
2109 }
2110
2111 BFT_FREE(w2);
2112 }
2113
2114 }
2115
2116 /*----------------------------------------------------------------------------*/
2117
2118 END_C_DECLS
2119