1 /*============================================================================
2 * Face viscosity
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_math.h"
52 #include "cs_blas.h"
53 #include "cs_halo.h"
54 #include "cs_halo_perio.h"
55 #include "cs_internal_coupling.h"
56 #include "cs_log.h"
57 #include "cs_mesh.h"
58 #include "cs_field.h"
59 #include "cs_field_pointer.h"
60 #include "cs_gradient.h"
61 #include "cs_ext_neighborhood.h"
62 #include "cs_mesh_quantities.h"
63 #include "cs_parall.h"
64 #include "cs_parameters.h"
65 #include "cs_porous_model.h"
66 #include "cs_prototypes.h"
67 #include "cs_timer.h"
68
69 /*----------------------------------------------------------------------------
70 * Header for the current file
71 *----------------------------------------------------------------------------*/
72
73 #include "cs_face_viscosity.h"
74
75 /*----------------------------------------------------------------------------*/
76
77 BEGIN_C_DECLS
78
79 /*=============================================================================
80 * Additional Doxygen documentation
81 *============================================================================*/
82
83 /*! \file cs_face_viscosity.c
84 *
85 * \brief Face viscosity.
86 *
87 */
88
89 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
90
91 /*=============================================================================
92 * Local type definitions
93 *============================================================================*/
94
95 /*============================================================================
96 * Private function definitions
97 *============================================================================*/
98
99 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
100
101 /*============================================================================
102 * Public function definitions for Fortran API
103 *============================================================================*/
104
105 /*----------------------------------------------------------------------------
106 * Wrapper to cs_face_viscosity
107 *----------------------------------------------------------------------------*/
108
CS_PROCF(viscfa,VISCFA)109 void CS_PROCF (viscfa, VISCFA)
110 (
111 const int *const visc_mean_type,
112 cs_real_t c_visc[],
113 cs_real_t i_visc[],
114 cs_real_t b_visc[]
115 )
116 {
117 const cs_mesh_t *m = cs_glob_mesh;
118 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
119
120 cs_face_viscosity(m,
121 fvq,
122 *visc_mean_type,
123 c_visc,
124 i_visc,
125 b_visc);
126 }
127
128 /*----------------------------------------------------------------------------
129 * Wrapper to cs_face_anisotropic_viscosity_vector
130 *----------------------------------------------------------------------------*/
131
CS_PROCF(vistnv,VISTNV)132 void CS_PROCF (vistnv, VISTNV)
133 (
134 const int *const visc_mean_type,
135 cs_real_6_t c_visc[],
136 cs_real_33_t i_visc[],
137 cs_real_t b_visc[]
138 )
139 {
140 const cs_mesh_t *m = cs_glob_mesh;
141 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
142
143 cs_face_anisotropic_viscosity_vector(m,
144 fvq,
145 *visc_mean_type,
146 c_visc,
147 i_visc,
148 b_visc);
149 }
150
151 /*----------------------------------------------------------------------------
152 * Wrapper to cs_face_anisotropic_viscosity_scalar
153 *----------------------------------------------------------------------------*/
154
CS_PROCF(vitens,VITENS)155 void CS_PROCF (vitens, VITENS)
156 (
157 cs_real_6_t c_visc[],
158 const int *const iwarnp,
159 cs_real_2_t weighf[],
160 cs_real_t weighb[],
161 cs_real_t i_visc[],
162 cs_real_t b_visc[]
163 )
164 {
165 const cs_mesh_t *m = cs_glob_mesh;
166 cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
167
168 cs_face_anisotropic_viscosity_scalar(m,
169 fvq,
170 c_visc,
171 *iwarnp,
172 weighf,
173 weighb,
174 i_visc,
175 b_visc);
176 }
177
178 /*============================================================================
179 * Public function definitions
180 *============================================================================*/
181
182 /*----------------------------------------------------------------------------*/
183 /*!
184 * \brief Compute the diffusion velocity at faces.
185 *
186 * i_visc,b_visc = viscosity*surface/distance, homogeneous to a rate of flow
187 * in kg/s.
188 *
189 * <a name="viscfa"></a>
190 *
191 * Please refer to the
192 * <a href="../../theory.pdf#viscfa"><b>viscfa</b></a> section of the theory
193 * guide for more informations.
194 *
195 * \remark: a priori, no need of reconstruction techniques
196 * (to improve if necessary).
197 *
198 * \param[in] m pointer to mesh
199 * \param[in] fvq pointer to finite volume quantities
200 * \param[in] visc_mean_type method to compute the viscosity at faces:
201 * - 0 arithmetical
202 * - 1 harmonic
203 * \param[in] c_visc cell viscosity (scalar)
204 * \param[out] i_visc inner face viscosity
205 * (times surface divided by distance)
206 * \param[out] b_visc boundary face viscosity
207 * (surface, must be consistent with flux BCs)
208 */
209 /*----------------------------------------------------------------------------*/
210
211 void
cs_face_viscosity(const cs_mesh_t * m,const cs_mesh_quantities_t * fvq,const int visc_mean_type,cs_real_t * restrict c_visc,cs_real_t * restrict i_visc,cs_real_t * restrict b_visc)212 cs_face_viscosity(const cs_mesh_t *m,
213 const cs_mesh_quantities_t *fvq,
214 const int visc_mean_type,
215 cs_real_t *restrict c_visc,
216 cs_real_t *restrict i_visc,
217 cs_real_t *restrict b_visc)
218 {
219 const cs_halo_t *halo = m->halo;
220
221 const cs_lnum_2_t *restrict i_face_cells
222 = (const cs_lnum_2_t *restrict)m->i_face_cells;
223 const cs_lnum_t *restrict b_face_cells
224 = (const cs_lnum_t *restrict)m->b_face_cells;
225 const cs_real_t *restrict weight = fvq->weight;
226 const cs_real_t *restrict i_dist = fvq->i_dist;
227 const cs_real_t *restrict i_f_face_surf = fvq->i_f_face_surf;
228 const cs_real_t *restrict b_f_face_surf = fvq->b_f_face_surf;
229
230 /* Porosity field */
231 cs_field_t *fporo = cs_field_by_name_try("porosity");
232
233 cs_real_t *porosi = NULL;
234
235 if (cs_glob_porous_model == 1 || cs_glob_porous_model == 2) {
236 porosi = fporo->val;
237 }
238
239 /* ---> Periodicity and parallelism treatment */
240 if (halo != NULL) {
241 cs_halo_type_t halo_type = CS_HALO_STANDARD;
242 cs_halo_sync_var(halo, halo_type, c_visc);
243 if (porosi != NULL) cs_halo_sync_var(halo, halo_type, porosi);
244 }
245
246 /* Without porosity */
247 if (porosi == NULL) {
248
249 if (visc_mean_type == 0) {
250
251 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
252
253 cs_lnum_t ii = i_face_cells[face_id][0];
254 cs_lnum_t jj = i_face_cells[face_id][1];
255
256 double visci = c_visc[ii];
257 double viscj = c_visc[jj];
258
259 i_visc[face_id] = 0.5*(visci+viscj)
260 *i_f_face_surf[face_id]/i_dist[face_id];
261
262 }
263
264 } else {
265
266 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
267
268 cs_lnum_t ii = i_face_cells[face_id][0];
269 cs_lnum_t jj = i_face_cells[face_id][1];
270
271 double visci = c_visc[ii];
272 double viscj = c_visc[jj];
273 double pnd = weight[face_id];
274
275 i_visc[face_id] = visci*viscj/CS_MAX(pnd*visci+(1.-pnd)*viscj,
276 DBL_MIN)
277 *i_f_face_surf[face_id]/i_dist[face_id];
278
279 }
280
281 }
282
283 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
284
285 b_visc[face_id] = b_f_face_surf[face_id];
286
287 }
288
289 /* With porosity */
290 } else {
291
292 if (visc_mean_type == 0) {
293
294 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
295
296 cs_lnum_t ii = i_face_cells[face_id][0];
297 cs_lnum_t jj = i_face_cells[face_id][1];
298
299 double visci = c_visc[ii] * porosi[ii];
300 double viscj = c_visc[jj] * porosi[jj];
301
302 i_visc[face_id] = 0.5*(visci+viscj)
303 *i_f_face_surf[face_id]/i_dist[face_id];
304
305 }
306
307 } else {
308
309 for (cs_lnum_t face_id = 0; face_id <m->n_i_faces; face_id++) {
310
311 cs_lnum_t ii = i_face_cells[face_id][0];
312 cs_lnum_t jj = i_face_cells[face_id][1];
313
314 double visci = c_visc[ii] * porosi[ii];
315 double viscj = c_visc[jj] * porosi[jj];
316 double pnd = weight[face_id];
317
318 i_visc[face_id] = visci*viscj/CS_MAX(pnd*visci+(1.-pnd)*viscj,
319 DBL_MIN)
320 *i_f_face_surf[face_id]/i_dist[face_id];
321
322 }
323
324 }
325
326 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
327
328 cs_lnum_t ii = b_face_cells[face_id];
329
330 b_visc[face_id] = b_f_face_surf[face_id]*porosi[ii];
331
332 }
333
334 }
335 }
336
337 /*----------------------------------------------------------------------------*/
338 /*!
339 * \brief Compute the equivalent tensor viscosity at faces for a 3x3 symetric
340 * tensor.
341 *
342 * \param[in] m pointer to mesh
343 * \param[in] fvq pointer to finite volume quantities
344 * \param[in] visc_mean_type method to compute the viscosity at faces:
345 * - 0: arithmetic
346 * - 1: harmonic
347 * \param[in] c_visc cell viscosity symmetric tensor
348 * \param[out] i_visc inner face tensor viscosity
349 * (times surface divided by distance)
350 * \param[out] b_visc boundary face viscosity
351 * (surface, must be consistent with flux BCs)
352 */
353 /*----------------------------------------------------------------------------*/
354
355 void
cs_face_anisotropic_viscosity_vector(const cs_mesh_t * m,const cs_mesh_quantities_t * fvq,const int visc_mean_type,cs_real_6_t * restrict c_visc,cs_real_33_t * restrict i_visc,cs_real_t * restrict b_visc)356 cs_face_anisotropic_viscosity_vector(const cs_mesh_t *m,
357 const cs_mesh_quantities_t *fvq,
358 const int visc_mean_type,
359 cs_real_6_t *restrict c_visc,
360 cs_real_33_t *restrict i_visc,
361 cs_real_t *restrict b_visc)
362 {
363 const cs_halo_t *halo = m->halo;
364
365 const cs_lnum_t n_cells = m->n_cells;
366 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
367
368 const cs_lnum_2_t *restrict i_face_cells
369 = (const cs_lnum_2_t *restrict)m->i_face_cells;
370 const cs_lnum_t *restrict b_face_cells
371 = (const cs_lnum_t *restrict)m->b_face_cells;
372 const cs_real_t *restrict weight = fvq->weight;
373 const cs_real_t *restrict i_dist = fvq->i_dist;
374 const cs_real_t *restrict i_f_face_surf = fvq->i_f_face_surf;
375 const cs_real_t *restrict b_f_face_surf = fvq->b_f_face_surf;
376
377 double visci[3][3], viscj[3][3], s1[6], s2[6];
378
379 cs_real_6_t *c_poro_visc = NULL;
380 cs_real_6_t *w2 = NULL;
381
382 /* Porosity fields */
383 cs_field_t *fporo = cs_field_by_name_try("porosity");
384 cs_field_t *ftporo = cs_field_by_name_try("tensorial_porosity");
385
386 cs_real_t *porosi = NULL;
387 cs_real_6_t *porosf = NULL;
388
389 if (cs_glob_porous_model == 1 || cs_glob_porous_model == 2) {
390 porosi = fporo->val;
391 if (ftporo != NULL) {
392 porosf = (cs_real_6_t *)ftporo->val;
393 }
394 }
395
396 /* Without porosity */
397 if (porosi == NULL) {
398
399 c_poro_visc = c_visc;
400
401 /* With porosity */
402 } else if (porosi != NULL && porosf == NULL) {
403
404 BFT_MALLOC(w2, n_cells_ext, cs_real_6_t);
405 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
406 for (int isou = 0; isou < 6; isou++) {
407 w2[cell_id][isou] = porosi[cell_id]*c_visc[cell_id][isou];
408 }
409 }
410 c_poro_visc = w2;
411
412 /* With tensorial porosity */
413 } else if (porosi != NULL && porosf != NULL) {
414
415 BFT_MALLOC(w2, n_cells_ext, cs_real_6_t);
416 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
417 cs_math_sym_33_product(porosf[cell_id],
418 c_visc[cell_id],
419 w2[cell_id]);
420 }
421 c_poro_visc = w2;
422
423 }
424
425 /* ---> Periodicity and parallelism treatment */
426 if (halo != NULL) {
427 cs_halo_type_t halo_type = CS_HALO_STANDARD;
428 cs_halo_sync_var_strided(halo, halo_type, (cs_real_t *)c_poro_visc, 6);
429 if (m->n_init_perio > 0)
430 cs_halo_perio_sync_var_sym_tens(halo,
431 halo_type,
432 (cs_real_t *)c_poro_visc);
433 }
434
435 /* Arithmetic mean */
436 if (visc_mean_type == 0) {
437
438 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
439
440 cs_lnum_t ii = i_face_cells[face_id][0];
441 cs_lnum_t jj = i_face_cells[face_id][1];
442
443 visci[0][0] = c_poro_visc[ii][0];
444 visci[1][1] = c_poro_visc[ii][1];
445 visci[2][2] = c_poro_visc[ii][2];
446 visci[1][0] = c_poro_visc[ii][3];
447 visci[0][1] = c_poro_visc[ii][3];
448 visci[2][1] = c_poro_visc[ii][4];
449 visci[1][2] = c_poro_visc[ii][4];
450 visci[2][0] = c_poro_visc[ii][5];
451 visci[0][2] = c_poro_visc[ii][5];
452
453 viscj[0][0] = c_poro_visc[jj][0];
454 viscj[1][1] = c_poro_visc[jj][1];
455 viscj[2][2] = c_poro_visc[jj][2];
456 viscj[1][0] = c_poro_visc[jj][3];
457 viscj[0][1] = c_poro_visc[jj][3];
458 viscj[2][1] = c_poro_visc[jj][4];
459 viscj[1][2] = c_poro_visc[jj][4];
460 viscj[2][0] = c_poro_visc[jj][5];
461 viscj[0][2] = c_poro_visc[jj][5];
462
463 for (int isou = 0; isou < 3; isou++) {
464 for (int jsou = 0; jsou < 3; jsou++) {
465 i_visc[face_id][jsou][isou] = 0.5*(visci[jsou][isou]
466 +viscj[jsou][isou])
467 * i_f_face_surf[face_id]/i_dist[face_id];
468 }
469 }
470
471 }
472
473 /* Harmonic mean: Kf = Ki . (pnd Ki +(1-pnd) Kj)^-1 . Kj */
474 } else {
475
476 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
477
478 cs_lnum_t ii = i_face_cells[face_id][0];
479 cs_lnum_t jj = i_face_cells[face_id][1];
480
481 double pnd = weight[face_id];
482
483 for (int isou = 0; isou < 6; isou++) {
484 s1[isou] = pnd*c_poro_visc[ii][isou] + (1.-pnd)*c_poro_visc[jj][isou];
485 }
486
487 cs_math_sym_33_inv_cramer(s1, s2);
488
489 cs_math_sym_33_product(s2, c_poro_visc[jj], s1);
490
491 cs_math_sym_33_product(c_poro_visc[ii], s1, s2);
492
493 double srfddi = i_f_face_surf[face_id]/i_dist[face_id];
494
495 i_visc[face_id][0][0] = s2[0]*srfddi;
496 i_visc[face_id][1][1] = s2[1]*srfddi;
497 i_visc[face_id][2][2] = s2[2]*srfddi;
498 i_visc[face_id][1][0] = s2[3]*srfddi;
499 i_visc[face_id][0][1] = s2[3]*srfddi;
500 i_visc[face_id][2][1] = s2[4]*srfddi;
501 i_visc[face_id][1][2] = s2[4]*srfddi;
502 i_visc[face_id][2][0] = s2[5]*srfddi;
503 i_visc[face_id][0][2] = s2[5]*srfddi;
504
505 }
506
507 }
508
509 /* Without porosity */
510 if (porosi == NULL) {
511
512 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
513
514 b_visc[face_id] = b_f_face_surf[face_id];
515
516 }
517
518 /* With porosity */
519 } else if (porosi != NULL && porosf == NULL) {
520
521 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
522
523 cs_lnum_t ii = b_face_cells[face_id];
524
525 b_visc[face_id] = b_f_face_surf[face_id]*porosi[ii];
526
527 }
528
529 /* With anisotropic porosity */
530 } else if (porosi != NULL && porosf != NULL) {
531
532 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
533
534 cs_lnum_t ii = b_face_cells[face_id];
535
536 b_visc[face_id] = b_f_face_surf[face_id]*porosi[ii];
537
538 }
539
540 }
541
542 BFT_FREE(w2);
543 }
544
545 /*----------------------------------------------------------------------------*/
546 /*!
547 * \brief Compute the equivalent viscosity at faces for a 3x3 symetric tensor,
548 * always using a harmonic mean.
549 *
550 * \param[in] m pointer to mesh
551 * \param[in] fvq pointer to finite volume quantities
552 * \param[in] c_visc cell viscosity symmetric tensor
553 * \param[in] iwarnp verbosity
554 * \param[out] weighf inner face weight between cells i and j
555 * \f$ \frac{\vect{IF} \cdot \tens{K}_\celli}
556 * {\norm{\tens{K}_\celli \cdot \vect{S}}^2} \f$
557 * and
558 * \f$ \frac{\vect{FJ} \cdot \tens{K}_\cellj}
559 * {\norm{\tens{K}_\cellj \cdot \vect{S}}^2} \f$
560 * \param[out] weighb boundary face weight
561 * \f$ \frac{\vect{IF} \cdot \tens{K}_\celli}
562 * {\norm{\tens{K}_\celli \cdot \vect{S}}^2} \f$
563 * \param[out] i_visc inner face viscosity
564 * (times surface divided by distance)
565 * \param[out] b_visc boundary face viscosity
566 * (surface, must be consistent with flux BCs)
567 */
568 /*----------------------------------------------------------------------------*/
569
570 void
cs_face_anisotropic_viscosity_scalar(const cs_mesh_t * m,const cs_mesh_quantities_t * fvq,cs_real_6_t * restrict c_visc,const int iwarnp,cs_real_2_t * restrict weighf,cs_real_t * restrict weighb,cs_real_t * restrict i_visc,cs_real_t * restrict b_visc)571 cs_face_anisotropic_viscosity_scalar(const cs_mesh_t *m,
572 const cs_mesh_quantities_t *fvq,
573 cs_real_6_t *restrict c_visc,
574 const int iwarnp,
575 cs_real_2_t *restrict weighf,
576 cs_real_t *restrict weighb,
577 cs_real_t *restrict i_visc,
578 cs_real_t *restrict b_visc)
579 {
580 const cs_halo_t *halo = m->halo;
581
582 const cs_lnum_t n_cells = m->n_cells;
583 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
584
585 const cs_lnum_2_t *restrict i_face_cells
586 = (const cs_lnum_2_t *restrict)m->i_face_cells;
587 const cs_lnum_t *restrict b_face_cells
588 = (const cs_lnum_t *restrict)m->b_face_cells;
589 const cs_real_t *restrict weight = fvq->weight;
590 const cs_real_t *restrict i_dist = fvq->i_dist;
591 const cs_real_t *restrict b_dist = fvq->b_dist;
592 const cs_real_t *restrict b_f_face_surf = fvq->b_f_face_surf;
593 const cs_real_3_t *restrict cell_cen
594 = (const cs_real_3_t *restrict)fvq->cell_cen;
595 const cs_real_3_t *restrict i_face_normal
596 = (const cs_real_3_t *restrict)fvq->i_face_normal;
597 const cs_real_t *restrict i_face_surf
598 = (const cs_real_t *restrict)fvq->i_face_surf;
599 const cs_real_t *restrict i_f_face_surf
600 = (const cs_real_t *restrict)fvq->i_f_face_surf;
601 const cs_real_3_t *restrict b_face_normal
602 = (const cs_real_3_t *restrict)fvq->b_face_normal;
603 const cs_real_3_t *restrict i_face_cog
604 = (const cs_real_3_t *restrict)fvq->i_face_cog;
605 const cs_real_3_t *restrict b_face_cog
606 = (const cs_real_3_t *restrict)fvq->b_face_cog;
607
608 cs_gnum_t nclipf = 0, nclipb = 0;
609 const double eps = 0.1;
610
611 cs_real_6_t *c_poro_visc = NULL;
612 cs_real_6_t *w2 = NULL;
613
614 /* Porosity fields */
615 cs_field_t *fporo = cs_field_by_name_try("porosity");
616 cs_field_t *ftporo = cs_field_by_name_try("tensorial_porosity");
617
618 cs_real_t *porosi = NULL;
619 cs_real_6_t *porosf = NULL;
620
621 if (cs_glob_porous_model == 1 || cs_glob_porous_model == 2) {
622 porosi = fporo->val;
623 if (ftporo != NULL) {
624 porosf = (cs_real_6_t *)ftporo->val;
625 }
626 }
627
628 /* Without porosity */
629 if (porosi == NULL) {
630
631 c_poro_visc = c_visc;
632
633 /* With porosity */
634 } else if (porosi != NULL && porosf == NULL) {
635
636 BFT_MALLOC(w2, n_cells_ext, cs_real_6_t);
637 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
638 for (int isou = 0; isou < 6; isou++) {
639 w2[cell_id][isou] = porosi[cell_id]*c_visc[cell_id][isou];
640 }
641 }
642 c_poro_visc = w2;
643
644 /* With tensorial porosity */
645 } else if (porosi != NULL && porosf != NULL) {
646
647 BFT_MALLOC(w2, n_cells_ext, cs_real_6_t);
648 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
649 cs_math_sym_33_product(porosf[cell_id],
650 c_visc[cell_id],
651 w2[cell_id]);
652 }
653 c_poro_visc = w2;
654
655 }
656
657 /* ---> Periodicity and parallelism treatment */
658 if (halo != NULL) {
659 cs_halo_type_t halo_type = CS_HALO_STANDARD;
660 cs_halo_sync_var_strided(halo, halo_type, (cs_real_t *)c_poro_visc, 6);
661 if (m->n_init_perio > 0)
662 cs_halo_perio_sync_var_sym_tens(halo,
663 halo_type,
664 (cs_real_t *)c_poro_visc);
665 }
666
667 /* Always Harmonic mean */
668 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++) {
669
670 cs_lnum_t ii = i_face_cells[face_id][0];
671 cs_lnum_t jj = i_face_cells[face_id][1];
672
673 /* ||Ki.S||^2 */
674 cs_real_3_t viscisv;
675 cs_math_sym_33_3_product(c_poro_visc[ii], i_face_normal[face_id], viscisv);
676 cs_real_t viscis = cs_math_3_square_norm(viscisv);
677
678 /* IF */
679 cs_real_3_t fi;
680 for (int kk = 0; kk < 3; kk++)
681 fi[kk] = i_face_cog[face_id][kk]-cell_cen[ii][kk];
682
683 /* IF.Ki.S */
684 cs_real_3_t fiki;
685 cs_math_sym_33_3_product(c_poro_visc[ii], fi, fiki);
686 cs_real_t fikis = cs_math_3_dot_product(fiki, i_face_normal[face_id]);
687
688 double distfi = (1. - weight[face_id])*i_dist[face_id];
689
690 /* Take I" so that I"F= eps*||FI||*Ki.n when I" is in cell rji */
691 double temp = eps*sqrt(viscis)*distfi;
692 if (fikis < temp) {
693 fikis = temp;
694 nclipf++;
695 }
696
697 /* ||Kj.S||^2 */
698 cs_real_3_t viscjsv;
699 cs_math_sym_33_3_product(c_poro_visc[jj], i_face_normal[face_id], viscjsv);
700 cs_real_t viscjs = cs_math_3_square_norm(viscjsv);
701
702 /* FJ */
703 cs_real_3_t fj;
704 for (int kk = 0; kk < 3; kk++)
705 fj[kk] = cell_cen[jj][kk]-i_face_cog[face_id][kk];
706
707 /* FJ.Kj.S */
708 cs_real_3_t fjkj;
709 cs_math_sym_33_3_product(c_poro_visc[jj], fj, fjkj);
710 cs_real_t fjkjs = cs_math_3_dot_product(fjkj, i_face_normal[face_id]);
711
712 double distfj = weight[face_id]*i_dist[face_id];
713
714 /* Take J" so that FJ"= eps*||FJ||*Kj.n when J" is in cell i */
715 temp = eps*sqrt(viscjs)*distfj;
716 if (fjkjs < temp) {
717 fjkjs = temp;
718 nclipf++;
719 }
720
721 weighf[face_id][0] = fikis/viscis;
722 weighf[face_id][1] = fjkjs/viscjs;
723
724 i_visc[face_id] = 1./(weighf[face_id][0] + weighf[face_id][1]);
725
726 }
727
728 /* For the porous modelling based on integral formulation Section and fluid
729 * Section are different */
730 if (cs_glob_porous_model == 3) {
731 for (cs_lnum_t face_id = 0; face_id < m->n_i_faces; face_id++)
732 i_visc[face_id] *= i_f_face_surf[face_id] / i_face_surf[face_id];
733 }
734
735 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
736
737 cs_lnum_t ii = b_face_cells[face_id];
738
739 /* ||Ki.S||^2 */
740 cs_real_3_t viscisv;
741 cs_math_sym_33_3_product(c_poro_visc[ii], b_face_normal[face_id], viscisv);
742 cs_real_t viscis = cs_math_3_square_norm(viscisv);
743
744 /* IF */
745 cs_real_3_t fi;
746 for (int kk = 0; kk < 3; kk++)
747 fi[kk] = b_face_cog[face_id][kk]-cell_cen[ii][kk];
748
749 /* IF.Ki.S */
750 cs_real_3_t fiki;
751 cs_math_sym_33_3_product(c_poro_visc[ii], fi, fiki);
752 cs_real_t fikis = cs_math_3_dot_product(fiki, b_face_normal[face_id]);
753
754 double distfi = b_dist[face_id];
755
756 /* Take I" so that I"F= eps*||FI||*Ki.n when J" is in cell rji */
757 double temp = eps*sqrt(viscis)*distfi;
758 if (fikis < temp) {
759 fikis = temp;
760 nclipb++;
761 }
762
763 weighb[face_id] = fikis/viscis;
764
765 }
766
767 /* Without porosity */
768 if (porosi == NULL) {
769
770 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
771
772 /* Warning: hint must be ||Ki.n||/I"F */
773 b_visc[face_id] = b_f_face_surf[face_id];
774 }
775
776 /* With porosity */
777 } else if (porosi != NULL && porosf == NULL) {
778
779 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
780
781 cs_lnum_t ii = b_face_cells[face_id];
782
783 /* Warning: hint must be ||Ki.n||/I"F */
784 b_visc[face_id] = b_f_face_surf[face_id]*porosi[ii];
785
786 }
787
788 /* With tensorial porosity */
789 } else if (porosi != NULL && porosf != NULL) {
790
791 for (cs_lnum_t face_id = 0; face_id < m->n_b_faces; face_id++) {
792
793 cs_lnum_t ii = b_face_cells[face_id];
794
795 /* Warning: hint must be ||Ki.n||/I"F */
796 b_visc[face_id] = b_f_face_surf[face_id]*porosi[ii];
797
798 }
799
800 }
801
802 if (halo != NULL) {
803 cs_parall_counter(&nclipf, 1);
804 cs_parall_counter(&nclipb, 1);
805 }
806
807 if (iwarnp >= 3) {
808 bft_printf("Computing the face viscosity from the tensorial viscosity:\n"
809 " Number of internal clippings: %lu\n"
810 " Number of boundary clippings: %lu\n",
811 (unsigned long)nclipf, (unsigned long)nclipb);
812 }
813
814 BFT_FREE(w2);
815 }
816
817 /*----------------------------------------------------------------------------*/
818
819 END_C_DECLS
820