1 /*============================================================================
2 * Matrix building
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_gradient.h"
58 #include "cs_ext_neighborhood.h"
59 #include "cs_mesh_quantities.h"
60 #include "cs_parameters.h"
61 #include "cs_porous_model.h"
62 #include "cs_prototypes.h"
63 #include "cs_timer.h"
64
65 #include "cs_parall.h"
66
67 /*----------------------------------------------------------------------------
68 * Header for the current file
69 *----------------------------------------------------------------------------*/
70
71 #include "cs_matrix_building.h"
72
73 /*----------------------------------------------------------------------------*/
74
75 BEGIN_C_DECLS
76
77 /*=============================================================================
78 * Additional Doxygen documentation
79 *============================================================================*/
80
81 /*! \file cs_matrix_building.c
82
83 */
84
85 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
86
87 /*=============================================================================
88 * Local Macro Definitions
89 *============================================================================*/
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 * Fortran wrapper to cs_matrix_scalar (or its counterpart for
107 * symmetric matrices)
108 *----------------------------------------------------------------------------*/
109
CS_PROCF(matrix,MATRIX)110 void CS_PROCF (matrix, MATRIX)
111 (
112 const int *iconvp,
113 const int *idiffp,
114 const int *ndircp,
115 const int *isym,
116 const cs_real_t *thetap,
117 const int *imucpp,
118 const cs_real_t coefbp[],
119 const cs_real_t cofbfp[],
120 const cs_real_t rovsdt[],
121 const cs_real_t i_massflux[],
122 const cs_real_t b_massflux[],
123 const cs_real_t i_visc[],
124 const cs_real_t b_visc[],
125 const cs_real_t xcpp[],
126 cs_real_t da[],
127 cs_real_t xa[]
128 )
129 {
130 cs_matrix_wrapper_scalar(*iconvp,
131 *idiffp,
132 *ndircp,
133 *isym,
134 *thetap,
135 *imucpp,
136 coefbp,
137 cofbfp,
138 rovsdt,
139 i_massflux,
140 b_massflux,
141 i_visc,
142 b_visc,
143 xcpp,
144 da,
145 xa);
146 }
147
148 /*----------------------------------------------------------------------------
149 * Fortran wrapper to cs_matrix_time_step
150 *----------------------------------------------------------------------------*/
151
CS_PROCF(matrdt,MATRDT)152 void CS_PROCF (matrdt, MATRDT)
153 (
154 const int *const iconvp,
155 const int *const idiffp,
156 const int *const isym,
157 const cs_real_t coefbp[],
158 const cs_real_t cofbfp[],
159 const cs_real_t i_massflux[],
160 const cs_real_t b_massflux[],
161 const cs_real_t i_visc[],
162 const cs_real_t b_visc[],
163 cs_real_t da[]
164 )
165 {
166 const cs_mesh_t *m = cs_glob_mesh;
167
168 cs_matrix_time_step(m,
169 *iconvp,
170 *idiffp,
171 *isym,
172 coefbp,
173 cofbfp,
174 i_massflux,
175 b_massflux,
176 i_visc,
177 b_visc,
178 da);
179 }
180
181 /*============================================================================
182 * Public function definitions
183 *============================================================================*/
184
185 /*----------------------------------------------------------------------------
186 * Wrapper to cs_matrix_scalar (or its counterpart for
187 * symmetric matrices)
188 *----------------------------------------------------------------------------*/
189
190 void
cs_matrix_wrapper_scalar(int iconvp,int idiffp,int ndircp,int isym,double thetap,int imucpp,const cs_real_t coefbp[],const cs_real_t cofbfp[],const cs_real_t rovsdt[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],const cs_real_t xcpp[],cs_real_t da[],cs_real_t xa[])191 cs_matrix_wrapper_scalar(int iconvp,
192 int idiffp,
193 int ndircp,
194 int isym,
195 double thetap,
196 int imucpp,
197 const cs_real_t coefbp[],
198 const cs_real_t cofbfp[],
199 const cs_real_t rovsdt[],
200 const cs_real_t i_massflux[],
201 const cs_real_t b_massflux[],
202 const cs_real_t i_visc[],
203 const cs_real_t b_visc[],
204 const cs_real_t xcpp[],
205 cs_real_t da[],
206 cs_real_t xa[])
207 {
208 const cs_mesh_t *m = cs_glob_mesh;
209 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
210 const cs_lnum_t n_cells = m->n_cells;
211
212 if (isym != 1 && isym != 2) {
213 bft_error(__FILE__, __LINE__, 0,
214 _("invalid value of isym"));
215 }
216
217 /* Symmetric matrix */
218 if (isym == 1) {
219 cs_sym_matrix_scalar(m,
220 idiffp,
221 thetap,
222 cofbfp,
223 rovsdt,
224 i_visc,
225 b_visc,
226 da,
227 xa);
228 }
229
230 /* Non-symmetric matrix */
231 else {
232 cs_matrix_scalar(m,
233 iconvp,
234 idiffp,
235 thetap,
236 imucpp,
237 coefbp,
238 cofbfp,
239 rovsdt,
240 i_massflux,
241 b_massflux,
242 i_visc,
243 b_visc,
244 xcpp,
245 da,
246 (cs_real_2_t*) xa);
247 }
248
249 /* Penalization if non invertible matrix */
250
251 /* If no Dirichlet condition, the diagonal is slightly increased in order
252 to shift the eigenvalues spectrum (if IDIRCL=0, we force NDIRCP to be at
253 least 1 in order not to shift the diagonal). */
254
255 if (ndircp <= 0) {
256 const double epsi = 1.e-7;
257
258 # pragma omp parallel for
259 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
260 da[cell_id] = (1.+epsi)*da[cell_id];
261 }
262 }
263
264 /* If a whole line of the matrix is 0, the diagonal is set to 1 */
265 if (mq->has_disable_flag == 1) {
266 # pragma omp parallel for
267 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
268 da[cell_id] += mq->c_disable_flag[cell_id];
269 }
270 }
271
272 }
273
274 /*----------------------------------------------------------------------------
275 * Wrapper to cs_matrix_scalar for convection/diffusion multigrid
276 *----------------------------------------------------------------------------*/
277
278 void
cs_matrix_wrapper_scalar_conv_diff(int iconvp,int idiffp,int ndircp,double thetap,int imucpp,const cs_real_t coefbp[],const cs_real_t cofbfp[],const cs_real_t rovsdt[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],const cs_real_t xcpp[],cs_real_t da[],cs_real_t xa[],cs_real_t da_conv[],cs_real_t xa_conv[],cs_real_t da_diff[],cs_real_t xa_diff[])279 cs_matrix_wrapper_scalar_conv_diff(int iconvp,
280 int idiffp,
281 int ndircp,
282 double thetap,
283 int imucpp,
284 const cs_real_t coefbp[],
285 const cs_real_t cofbfp[],
286 const cs_real_t rovsdt[],
287 const cs_real_t i_massflux[],
288 const cs_real_t b_massflux[],
289 const cs_real_t i_visc[],
290 const cs_real_t b_visc[],
291 const cs_real_t xcpp[],
292 cs_real_t da[],
293 cs_real_t xa[],
294 cs_real_t da_conv[],
295 cs_real_t xa_conv[],
296 cs_real_t da_diff[],
297 cs_real_t xa_diff[])
298 {
299 const cs_mesh_t *m = cs_glob_mesh;
300 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
301 const cs_lnum_t n_cells = m->n_cells;
302 const cs_lnum_t n_i_faces = m->n_i_faces;
303 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
304
305 /* Diffusion matrix */
306 cs_sym_matrix_scalar(m,
307 idiffp,
308 thetap,
309 cofbfp,
310 rovsdt,
311 i_visc,
312 b_visc,
313 da_diff,
314 xa_diff);
315
316 /* Convection matrix */
317 cs_matrix_scalar(m,
318 iconvp,
319 0.,
320 thetap,
321 imucpp,
322 coefbp,
323 cofbfp,
324 rovsdt,
325 i_massflux,
326 b_massflux,
327 i_visc,
328 b_visc,
329 xcpp,
330 da_conv,
331 (cs_real_2_t*)xa_conv);
332
333 /* Global matrix */
334
335 # pragma omp parallel for
336 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
337 da_conv[cell_id] -= rovsdt[cell_id];
338 da_diff[cell_id] -= rovsdt[cell_id];
339 da[cell_id] = rovsdt[cell_id] + da_diff[cell_id] + da_conv[cell_id];
340 }
341 if (n_cells_ext > n_cells) {
342 # pragma omp parallel for if (n_cells_ext - n_cells > CS_THR_MIN)
343 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
344 da[cell_id] = 0.;
345 }
346 }
347
348 # pragma omp parallel for
349 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
350 xa[2*face_id] = xa_diff[face_id] + xa_conv[2*face_id];
351 xa[2*face_id +1] = xa_diff[face_id] + xa_conv[2*face_id +1];
352 }
353
354 /* Penalization if non invertible matrix */
355
356 /* If no Dirichlet condition, the diagonal is slightly increased in order
357 to shift the eigenvalues spectrum (if IDIRCL=0, we force NDIRCP to be at
358 least 1 in order not to shift the diagonal). */
359
360 if (ndircp <= 0) {
361 const double epsi = 1.e-7;
362
363 # pragma omp parallel for
364 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
365 da[cell_id] = (1.+epsi)*da[cell_id];
366 }
367 }
368
369 /* If a whole line of the matrix is 0, the diagonal is set to 1 */
370 if (mq->has_disable_flag == 1) {
371 # pragma omp parallel for
372 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
373 da[cell_id] += mq->c_disable_flag[cell_id];
374 }
375 }
376
377 }
378
379 /*----------------------------------------------------------------------------
380 * Wrapper to cs_matrix_vector (or its counterpart for
381 * symmetric matrices)
382 *----------------------------------------------------------------------------*/
383
384 void
cs_matrix_wrapper_vector(int iconvp,int idiffp,int tensorial_diffusion,int ndircp,int isym,cs_lnum_t eb_size[4],double thetap,const cs_real_33_t coefbp[],const cs_real_33_t cofbfp[],const cs_real_33_t fimp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_33_t da[],cs_real_t xa[])385 cs_matrix_wrapper_vector(int iconvp,
386 int idiffp,
387 int tensorial_diffusion,
388 int ndircp,
389 int isym,
390 cs_lnum_t eb_size[4],
391 double thetap,
392 const cs_real_33_t coefbp[],
393 const cs_real_33_t cofbfp[],
394 const cs_real_33_t fimp[],
395 const cs_real_t i_massflux[],
396 const cs_real_t b_massflux[],
397 const cs_real_t i_visc[],
398 const cs_real_t b_visc[],
399 cs_real_33_t da[],
400 cs_real_t xa[])
401 {
402 const cs_mesh_t *m = cs_glob_mesh;
403 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
404 const cs_lnum_t n_cells = m->n_cells;
405
406 if (isym != 1 && isym != 2) {
407 bft_error(__FILE__, __LINE__, 0,
408 _("invalid value of isym"));
409 }
410
411 /* scalar diffusion or right anisotropic diffusion */
412 if (tensorial_diffusion == 1) {
413 /* Symmetric matrix */
414 if (isym == 1) {
415 assert(eb_size[0] == 1);
416 cs_sym_matrix_vector(m,
417 idiffp,
418 thetap,
419 cofbfp,
420 fimp,
421 i_visc,
422 b_visc,
423 da,
424 xa);
425
426 /* Non-symmetric matrix */
427 }
428 else {
429 cs_matrix_vector(m,
430 mq,
431 iconvp,
432 idiffp,
433 eb_size,
434 thetap,
435 coefbp,
436 cofbfp,
437 fimp,
438 i_massflux,
439 b_massflux,
440 i_visc,
441 b_visc,
442 da,
443 (cs_real_2_t*) xa);
444 }
445 }
446 /* left tensor diffusion */
447 else {
448
449 /* Symmetric matrix */
450 if (isym == 1) {
451 cs_sym_matrix_anisotropic_diffusion(m,
452 idiffp,
453 thetap,
454 cofbfp,
455 fimp,
456 (const cs_real_33_t *)i_visc,
457 b_visc,
458 da,
459 (cs_real_33_t *) xa);
460
461 /* Non-symmetric matrix */
462 } else {
463 cs_matrix_anisotropic_diffusion(m,
464 mq,
465 iconvp,
466 idiffp,
467 thetap,
468 coefbp,
469 cofbfp,
470 fimp,
471 i_massflux,
472 b_massflux,
473 (const cs_real_33_t *)i_visc,
474 b_visc,
475 da,
476 (cs_real_332_t *) xa);
477 }
478
479 }
480
481 /* Penalization if non invertible matrix */
482
483 /* If no Dirichlet condition, the diagonal is slightly increased in order
484 to shift the eigenvalues spectrum. */
485
486 if (ndircp <= 0) {
487 const double epsi = 1.e-7;
488 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
489 for (int isou = 0; isou < 3; isou++) {
490 da[cell_id][isou][isou] = (1.+epsi)*da[cell_id][isou][isou];
491 }
492 }
493 }
494
495 /* If a whole line of the matrix is 0, the diagonal is set to 1 */
496 if (mq->has_disable_flag == 1) {
497 # pragma omp parallel for
498 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
499 for (int isou = 0; isou < 3; isou++)
500 da[cell_id][isou][isou]
501 += (cs_real_t)(mq->c_disable_flag[cell_id]);
502 }
503 }
504
505 }
506
507 /*----------------------------------------------------------------------------
508 * Wrapper to cs_matrix_tensor (or its counterpart for
509 * symmetric matrices)
510 *----------------------------------------------------------------------------*/
511
512 void
cs_matrix_wrapper_tensor(int iconvp,int idiffp,int tensorial_diffusion,int ndircp,int isym,double thetap,const cs_real_66_t coefbts[],const cs_real_66_t cofbfts[],const cs_real_66_t fimp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_66_t da[],cs_real_t xa[])513 cs_matrix_wrapper_tensor(int iconvp,
514 int idiffp,
515 int tensorial_diffusion,
516 int ndircp,
517 int isym,
518 double thetap,
519 const cs_real_66_t coefbts[],
520 const cs_real_66_t cofbfts[],
521 const cs_real_66_t fimp[],
522 const cs_real_t i_massflux[],
523 const cs_real_t b_massflux[],
524 const cs_real_t i_visc[],
525 const cs_real_t b_visc[],
526 cs_real_66_t da[],
527 cs_real_t xa[])
528 {
529 const cs_mesh_t *m = cs_glob_mesh;
530 const cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
531 const cs_lnum_t n_cells = m->n_cells;
532
533 if (isym != 1 && isym != 2) {
534 bft_error(__FILE__, __LINE__, 0,
535 _("invalid value of isym"));
536 }
537
538 /* scalar diffusion or right anisotropic diffusion */
539 if (tensorial_diffusion == 1) {
540 /* Symmetric matrix */
541 if (isym == 1) {
542 cs_sym_matrix_tensor(m,
543 idiffp,
544 thetap,
545 cofbfts,
546 fimp,
547 i_visc,
548 b_visc,
549 da,
550 xa);
551
552 /* Non-symmetric matrix */
553 } else {
554 cs_matrix_tensor(m,
555 iconvp,
556 idiffp,
557 thetap,
558 coefbts,
559 cofbfts,
560 fimp,
561 i_massflux,
562 b_massflux,
563 i_visc,
564 b_visc,
565 da,
566 (cs_real_2_t*) xa);
567 }
568 }
569 /* left tensor diffusion */
570 else {
571 /* Symmetric matrix */
572 if (isym == 1) {
573 cs_sym_matrix_anisotropic_diffusion_tensor(m,
574 idiffp,
575 thetap,
576 cofbfts,
577 fimp,
578 (const cs_real_66_t *)i_visc,
579 b_visc,
580 da,
581 (cs_real_66_t *)xa);
582
583 /* Non-symmetric matrix */
584 }
585 else {
586 cs_matrix_anisotropic_diffusion_tensor(m,
587 iconvp,
588 idiffp,
589 thetap,
590 coefbts,
591 cofbfts,
592 fimp,
593 i_massflux,
594 b_massflux,
595 (cs_real_66_t *)i_visc,
596 b_visc,
597 da,
598 (cs_real_662_t *)xa);
599 }
600 }
601
602 /* Penalization if non invertible matrix */
603
604 /* If no Dirichlet condition, the diagonal is slightly increased in order
605 to shift the eigenvalues spectrum (if IDIRCL=0, we force NDIRCP to be at
606 least 1 in order not to shift the diagonal). */
607
608 if (ndircp <= 0) {
609 const double epsi = 1.e-7;
610
611 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
612 for (int isou = 0; isou < 6; isou++) {
613 da[cell_id][isou][isou] = (1.+epsi)*da[cell_id][isou][isou];
614 }
615 }
616 }
617
618 /* If a whole line of the matrix is 0, the diagonal is set to 1 */
619 if (mq->has_disable_flag == 1) {
620 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
621 for (int isou = 0; isou < 6; isou++) {
622 da[cell_id][isou][isou] += mq->c_disable_flag[cell_id];
623 }
624 }
625 }
626
627 }
628
629 /*----------------------------------------------------------------------------*/
630 /*!
631 * \brief Build the diffusion matrix for a scalar field.
632 * (symmetric matrix).
633 *
634 * The diffusion is not reconstructed.
635 * The matrix is split into a diagonal block (number of cells)
636 * and an extra diagonal part (of dimension the number of internal
637 * faces).
638 *
639 * \param[in] m pointer to mesh structure
640 * \param[in] idiffp indicator
641 * - 1 diffusion
642 * - 0 otherwise
643 * \param[in] thetap weighting coefficient for the theta-scheme,
644 * - thetap = 0: explicit scheme
645 * - thetap = 0.5: time-centered
646 * scheme (mix between Crank-Nicolson and
647 * Adams-Bashforth)
648 * - thetap = 1: implicit scheme
649 * \param[in] cofbfp boundary condition array for the variable flux
650 * (implicit part)
651 * \param[in] rovsdt working array
652 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
653 * at interior faces for the matrix
654 * \param[in] b_visc \f$ S_\fib \f$
655 * at border faces for the matrix
656 * \param[out] da diagonal part of the matrix
657 * \param[out] xa extra diagonal part of the matrix
658 */
659 /*----------------------------------------------------------------------------*/
660
661 void
cs_sym_matrix_scalar(const cs_mesh_t * m,int idiffp,double thetap,const cs_real_t cofbfp[],const cs_real_t rovsdt[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_t * restrict da,cs_real_t * restrict xa)662 cs_sym_matrix_scalar(const cs_mesh_t *m,
663 int idiffp,
664 double thetap,
665 const cs_real_t cofbfp[],
666 const cs_real_t rovsdt[],
667 const cs_real_t i_visc[],
668 const cs_real_t b_visc[],
669 cs_real_t *restrict da,
670 cs_real_t *restrict xa)
671 {
672 const cs_lnum_t n_cells = m->n_cells;
673 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
674 const cs_lnum_t n_i_faces = m->n_i_faces;
675 const int n_i_groups = m->i_face_numbering->n_groups;
676 const int n_i_threads = m->i_face_numbering->n_threads;
677 const int n_b_groups = m->b_face_numbering->n_groups;
678 const int n_b_threads = m->b_face_numbering->n_threads;
679 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
680 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
681
682 const cs_lnum_2_t *restrict i_face_cells
683 = (const cs_lnum_2_t *restrict)m->i_face_cells;
684 const cs_lnum_t *restrict b_face_cells
685 = (const cs_lnum_t *restrict)m->b_face_cells;
686
687 /* 1. Initialization */
688
689 # pragma omp parallel for
690 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
691 da[cell_id] = rovsdt[cell_id];
692 }
693 if (n_cells_ext > n_cells) {
694 # pragma omp parallel for if (n_cells_ext - n_cells > CS_THR_MIN)
695 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
696 da[cell_id] = 0.;
697 }
698 }
699
700 /* 2. Computation of extradiagonal terms and contribution to the diagonal */
701
702 if (idiffp) {
703
704 for (int g_id = 0; g_id < n_i_groups; g_id++) {
705 # pragma omp parallel for firstprivate(thetap)
706 for (int t_id = 0; t_id < n_i_threads; t_id++) {
707 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
708 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
709 face_id++) {
710
711 cs_lnum_t ii = i_face_cells[face_id][0];
712 cs_lnum_t jj = i_face_cells[face_id][1];
713
714 cs_real_t aij = -thetap*i_visc[face_id];
715
716 xa[face_id] = aij;
717 da[ii] -= aij;
718 da[jj] -= aij;
719
720 }
721 }
722 }
723
724 }
725
726 else {
727
728 # pragma omp parallel for
729 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
730 xa[face_id] = 0.;
731 }
732
733 }
734
735 /* 3. Contribution of border faces to the diagonal */
736
737 if (idiffp) {
738
739 for (int g_id = 0; g_id < n_b_groups; g_id++) {
740 # pragma omp parallel for firstprivate(thetap, idiffp) \
741 if (m->n_b_faces > CS_THR_MIN)
742 for (int t_id = 0; t_id < n_b_threads; t_id++) {
743 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
744 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
745 face_id++) {
746
747 cs_lnum_t ii = b_face_cells[face_id];
748
749 da[ii] += thetap*b_visc[face_id]*cofbfp[face_id];
750
751 }
752 }
753 }
754
755 }
756
757 }
758
759 /*----------------------------------------------------------------------------*/
760 /*!
761 * \brief Build the advection/diffusion matrix for a scalar field
762 * (non-symmetric matrix).
763 *
764 * The advection is upwind, the diffusion is not reconstructed.
765 * The matrix is split into a diagonal block (number of cells)
766 * and an extra diagonal part (of dimension 2 time the number of internal
767 * faces).
768 *
769 * \param[in] m pointer to mesh structure
770 * \param[in] iconvp indicator
771 * - 1 advection
772 * - 0 otherwise
773 * \param[in] idiffp indicator
774 * - 1 diffusion
775 * - 0 otherwise
776 * \param[in] thetap weighting coefficient for the theta-scheme,
777 * - thetap = 0: explicit scheme
778 * - thetap = 0.5: time-centered
779 * scheme (mix between Crank-Nicolson and
780 * Adams-Bashforth)
781 * - thetap = 1: implicit scheme
782 * \param[in] imucpp indicator
783 * - 0 do not multiply the convective term by Cp
784 * - 1 do multiply the convective term by Cp
785 * \param[in] coefbp boundary condition array for the variable
786 * (implicit part)
787 * \param[in] cofbfp boundary condition array for the variable flux
788 * (implicit part)
789 * \param[in] rovsdt working array
790 * \param[in] i_massflux mass flux at interior faces
791 * \param[in] b_massflux mass flux at border faces
792 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
793 * at interior faces for the matrix
794 * \param[in] b_visc \f$ S_\fib \f$
795 * at border faces for the matrix
796 * \param[in] xcpp array of specific heat (Cp)
797 * \param[out] da diagonal part of the matrix
798 * \param[out] xa extra interleaved diagonal part of the matrix
799 */
800 /*----------------------------------------------------------------------------*/
801
802 void
cs_matrix_scalar(const cs_mesh_t * m,int iconvp,int idiffp,double thetap,int imucpp,const cs_real_t coefbp[],const cs_real_t cofbfp[],const cs_real_t rovsdt[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],const cs_real_t xcpp[],cs_real_t * restrict da,cs_real_2_t * restrict xa)803 cs_matrix_scalar(const cs_mesh_t *m,
804 int iconvp,
805 int idiffp,
806 double thetap,
807 int imucpp,
808 const cs_real_t coefbp[],
809 const cs_real_t cofbfp[],
810 const cs_real_t rovsdt[],
811 const cs_real_t i_massflux[],
812 const cs_real_t b_massflux[],
813 const cs_real_t i_visc[],
814 const cs_real_t b_visc[],
815 const cs_real_t xcpp[],
816 cs_real_t *restrict da,
817 cs_real_2_t *restrict xa)
818 {
819 const cs_lnum_t n_cells = m->n_cells;
820 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
821 const cs_lnum_t n_i_faces = m->n_i_faces;
822 const int n_i_groups = m->i_face_numbering->n_groups;
823 const int n_i_threads = m->i_face_numbering->n_threads;
824 const int n_b_groups = m->b_face_numbering->n_groups;
825 const int n_b_threads = m->b_face_numbering->n_threads;
826 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
827 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
828
829 const cs_lnum_2_t *restrict i_face_cells
830 = (const cs_lnum_2_t *restrict)m->i_face_cells;
831 const cs_lnum_t *restrict b_face_cells
832 = (const cs_lnum_t *restrict)m->b_face_cells;
833
834 /* 1. Initialization */
835
836 # pragma omp parallel for
837 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
838 da[cell_id] = rovsdt[cell_id];
839 }
840 if (n_cells_ext > n_cells) {
841 # pragma omp parallel for if (n_cells_ext - n_cells > CS_THR_MIN)
842 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
843 da[cell_id] = 0.;
844 }
845 }
846
847 # pragma omp parallel for
848 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
849 xa[face_id][0] = 0.;
850 xa[face_id][1] = 0.;
851 }
852
853 /* When solving the temperature, the convective part is multiplied by Cp */
854 if (imucpp == 0) {
855
856 /* 2. Computation of extradiagonal terms */
857
858 # pragma omp parallel for firstprivate(thetap, iconvp, idiffp)
859 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
860
861 double flui = 0.5*(i_massflux[face_id] -fabs(i_massflux[face_id]));
862 double fluj =-0.5*(i_massflux[face_id] +fabs(i_massflux[face_id]));
863
864 xa[face_id][0] = thetap*(iconvp*flui -idiffp*i_visc[face_id]);
865 xa[face_id][1] = thetap*(iconvp*fluj -idiffp*i_visc[face_id]);
866
867 }
868
869 /* 3. Contribution of the extra-diagonal terms to the diagonal */
870
871 for (int g_id = 0; g_id < n_i_groups; g_id++) {
872 # pragma omp parallel for
873 for (int t_id = 0; t_id < n_i_threads; t_id++) {
874 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
875 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
876 face_id++) {
877
878 cs_lnum_t ii = i_face_cells[face_id][0];
879 cs_lnum_t jj = i_face_cells[face_id][1];
880
881 /* D_ii = theta (m_ij)^+ - m_ij
882 * = -X_ij - (1-theta)*m_ij
883 * D_jj = -theta (m_ij)^- + m_ij
884 * = -X_ji + (1-theta)*m_ij
885 */
886 da[ii] -= xa[face_id][0] + iconvp*(1. - thetap)*i_massflux[face_id];
887 da[jj] -= xa[face_id][1] - iconvp*(1. - thetap)*i_massflux[face_id];
888
889 }
890 }
891 }
892
893 /* 4. Contribution of border faces to the diagonal */
894
895 for (int g_id = 0; g_id < n_b_groups; g_id++) {
896 # pragma omp parallel for firstprivate(thetap, iconvp, idiffp) \
897 if (m->n_b_faces > CS_THR_MIN)
898 for (int t_id = 0; t_id < n_b_threads; t_id++) {
899 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
900 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
901 face_id++) {
902
903 cs_lnum_t ii = b_face_cells[face_id];
904
905 double flui = 0.5*(b_massflux[face_id] - fabs(b_massflux[face_id]));
906
907 /* D_ii = theta (m_f)^+ + theta B (m_f)^- - m_f
908 * = (theta B -1)*(m_f)^- - (1-theta)*(m_f)^+
909 * = theta*(B -1)*(m_f)^- - (1-theta)*m_f
910 */
911 da[ii] += iconvp*(flui*thetap*(coefbp[face_id]-1.)
912 -(1.-thetap)*b_massflux[face_id])
913 + idiffp*thetap*b_visc[face_id]*cofbfp[face_id];
914 }
915 }
916 }
917
918 /* When solving the temperature, the convective part is multiplied by Cp */
919 } else {
920
921 /* 2. Computation of extradiagonal terms */
922
923 # pragma omp parallel for firstprivate(thetap, iconvp, idiffp)
924 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
925
926 double flui = 0.5*(i_massflux[face_id] -fabs(i_massflux[face_id]));
927 double fluj =-0.5*(i_massflux[face_id] +fabs(i_massflux[face_id]));
928
929 cs_lnum_t ii = i_face_cells[face_id][0];
930 cs_lnum_t jj = i_face_cells[face_id][1];
931
932 xa[face_id][0] = thetap*( iconvp*xcpp[ii]*flui
933 -idiffp*i_visc[face_id]);
934 xa[face_id][1] = thetap*( iconvp*xcpp[jj]*fluj
935 -idiffp*i_visc[face_id]);
936
937 }
938
939 /* 3. Contribution of the extra-diagonal terms to the diagonal */
940
941 for (int g_id = 0; g_id < n_i_groups; g_id++) {
942 # pragma omp parallel for
943 for (int t_id = 0; t_id < n_i_threads; t_id++) {
944 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
945 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
946 face_id++) {
947
948 cs_lnum_t ii = i_face_cells[face_id][0];
949 cs_lnum_t jj = i_face_cells[face_id][1];
950
951 /* D_ii = theta (m_ij)^+ - m_ij
952 * = -X_ij - (1-theta)*m_ij
953 * D_jj = -theta (m_ij)^- + m_ij
954 * = -X_ji + (1-theta)*m_ij
955 */
956 da[ii] -= xa[face_id][0]
957 + iconvp*(1. - thetap)*xcpp[ii]*i_massflux[face_id];
958 da[jj] -= xa[face_id][1]
959 - iconvp*(1. - thetap)*xcpp[jj]*i_massflux[face_id];
960
961 }
962 }
963 }
964
965 /* 4. Contribution of boundary faces to the diagonal */
966
967 for (int g_id = 0; g_id < n_b_groups; g_id++) {
968 # pragma omp parallel for firstprivate(thetap, iconvp, idiffp) \
969 if (m->n_b_faces > CS_THR_MIN)
970 for (int t_id = 0; t_id < n_b_threads; t_id++) {
971 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
972 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
973 face_id++) {
974
975 cs_lnum_t ii = b_face_cells[face_id];
976
977 double flui = 0.5*(b_massflux[face_id] - fabs(b_massflux[face_id]));
978
979 /* D_ii = theta (m_f)^+ + theta B (m_f)^- - m_f
980 * = (theta B -1)*(m_f)^- - (1-theta)*(m_f)^+
981 * = theta*(B -1)*(m_f)^- - (1-theta)*m_f
982 */
983 da[ii] += iconvp*xcpp[ii]*(flui*thetap*(coefbp[face_id]-1.)
984 -(1.-thetap)*b_massflux[face_id])
985 + idiffp*thetap*b_visc[face_id]*cofbfp[face_id];
986 }
987 }
988 }
989
990 }
991
992 }
993
994 /*----------------------------------------------------------------------------*/
995 /*!
996 * \brief Build the diffusion matrix for a vector field
997 * (symmetric matrix).
998 *
999 * The diffusion is not reconstructed.
1000 * The matrix is split into a diagonal block (3x3 times number of cells)
1001 * and an extra diagonal part (of dimension the number of internal
1002 * faces).
1003 *
1004 * \param[in] m pointer to mesh structure
1005 * \param[in] idiffp indicator
1006 * - 1 diffusion
1007 * - 0 otherwise
1008 * \param[in] thetap weighting coefficient for the theta-scheme,
1009 * - thetap = 0: explicit scheme
1010 * - thetap = 0.5: time-centered
1011 * scheme (mix between Crank-Nicolson and
1012 * Adams-Bashforth)
1013 * - thetap = 1: implicit scheme
1014 * \param[in] cofbfp boundary condition array for the variable flux
1015 * (implicit part - 3x3 tensor array)
1016 * \param[in] fimp \f$ \tens{f_s}^{imp} \f$
1017 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1018 * at interior faces for the matrix
1019 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1020 * at border faces for the matrix
1021 * \param[out] da diagonal part of the matrix
1022 * \param[out] xa extra interleaved diagonal part of the matrix
1023 */
1024 /*----------------------------------------------------------------------------*/
1025
1026 void
cs_sym_matrix_vector(const cs_mesh_t * m,int idiffp,double thetap,const cs_real_33_t cofbfp[],const cs_real_33_t fimp[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_33_t * restrict da,cs_real_t * restrict xa)1027 cs_sym_matrix_vector(const cs_mesh_t *m,
1028 int idiffp,
1029 double thetap,
1030 const cs_real_33_t cofbfp[],
1031 const cs_real_33_t fimp[],
1032 const cs_real_t i_visc[],
1033 const cs_real_t b_visc[],
1034 cs_real_33_t *restrict da,
1035 cs_real_t *restrict xa)
1036 {
1037 const cs_lnum_t n_cells = m->n_cells;
1038 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1039 const cs_lnum_t n_i_faces = m->n_i_faces;
1040 const cs_lnum_t n_b_faces = m->n_b_faces;
1041
1042 const cs_lnum_2_t *restrict i_face_cells
1043 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1044 const cs_lnum_t *restrict b_face_cells
1045 = (const cs_lnum_t *restrict)m->b_face_cells;
1046
1047 /* 1. Initialization */
1048
1049 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1050
1051 for (int isou = 0; isou < 3; isou++) {
1052 for (int jsou = 0; jsou < 3; jsou++) {
1053 da[cell_id][jsou][isou] = fimp[cell_id][jsou][isou];
1054 }
1055 }
1056
1057 }
1058
1059 if (n_cells_ext > n_cells) {
1060 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
1061
1062 for (int isou = 0; isou < 3; isou++) {
1063 for (int jsou = 0; jsou < 3; jsou++) {
1064 da[cell_id][jsou][isou] = 0.;
1065 }
1066 }
1067
1068 }
1069 }
1070
1071 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1072 xa[face_id] = 0.;
1073 }
1074
1075 /* 2. Computation of extradiagonal terms */
1076
1077 for (cs_lnum_t face_id = 0; face_id <n_i_faces; face_id++) {
1078
1079 xa[face_id] = -thetap*idiffp*i_visc[face_id];
1080
1081 }
1082
1083 /* 3. Contribution of the extra-diagonal terms to the diagonal */
1084
1085 for (cs_lnum_t face_id = 0; face_id <n_i_faces; face_id++) {
1086
1087 cs_lnum_t ii = i_face_cells[face_id][0];
1088 cs_lnum_t jj = i_face_cells[face_id][1];
1089
1090 for (int isou = 0; isou < 3; isou++) {
1091 da[ii][isou][isou] -= xa[face_id];
1092 da[jj][isou][isou] -= xa[face_id];
1093 }
1094
1095 }
1096
1097 /* 4. Contribution of border faces to the diagonal */
1098
1099 for (cs_lnum_t face_id = 0; face_id <n_b_faces; face_id++) {
1100
1101 cs_lnum_t ii = b_face_cells[face_id];
1102
1103 for (int isou = 0; isou < 3; isou++) {
1104 for (int jsou = 0; jsou < 3; jsou++) {
1105 da[ii][jsou][isou] += thetap*idiffp*b_visc[face_id]
1106 *cofbfp[face_id][jsou][isou];
1107 }
1108 }
1109
1110 }
1111 }
1112
1113 /*----------------------------------------------------------------------------*/
1114 /*!
1115 * \brief Build the diffusion matrix for a tensor field
1116 * (symmetric matrix).
1117 *
1118 * The diffusion is not reconstructed.
1119 * The matrix is split into a diagonal block (6x6 times number of cells)
1120 * and an extra diagonal part (of dimension the number of internal
1121 * faces).
1122 *
1123 * \param[in] m pointer to mesh structure
1124 * \param[in] idiffp indicator
1125 * - 1 diffusion
1126 * - 0 otherwise
1127 * \param[in] thetap weighting coefficient for the theta-scheme,
1128 * - thetap = 0: explicit scheme
1129 * - thetap = 0.5: time-centered
1130 * scheme (mix between Crank-Nicolson and
1131 * Adams-Bashforth)
1132 * - thetap = 1: implicit scheme
1133 * \param[in] cofbfts boundary condition array for the variable flux
1134 * (Implicit part - 6x6 tensor array)
1135 * \param[in] fimp part of the diagonal
1136 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1137 * at interior faces for the matrix
1138 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1139 * at border faces for the matrix
1140 * \param[out] da diagonal part of the matrix
1141 * \param[out] xa extra interleaved diagonal part of the matrix
1142 */
1143 /*----------------------------------------------------------------------------*/
1144
1145 void
cs_sym_matrix_tensor(const cs_mesh_t * m,int idiffp,double thetap,const cs_real_66_t cofbfts[],const cs_real_66_t fimp[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_66_t * restrict da,cs_real_t * restrict xa)1146 cs_sym_matrix_tensor(const cs_mesh_t *m,
1147 int idiffp,
1148 double thetap,
1149 const cs_real_66_t cofbfts[],
1150 const cs_real_66_t fimp[],
1151 const cs_real_t i_visc[],
1152 const cs_real_t b_visc[],
1153 cs_real_66_t *restrict da,
1154 cs_real_t *restrict xa)
1155 {
1156 const cs_lnum_t n_cells = m->n_cells;
1157 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1158 const cs_lnum_t n_i_faces = m->n_i_faces;
1159 const cs_lnum_t n_b_faces = m->n_b_faces;
1160
1161 const cs_lnum_2_t *restrict i_face_cells
1162 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1163 const cs_lnum_t *restrict b_face_cells
1164 = (const cs_lnum_t *restrict)m->b_face_cells;
1165
1166 /* 1. Initialization */
1167
1168 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1169
1170 for (int isou = 0; isou < 6; isou++) {
1171 for (int jsou = 0; jsou < 6; jsou++) {
1172 da[cell_id][jsou][isou] = fimp[cell_id][jsou][isou];
1173 }
1174 }
1175
1176 }
1177
1178 if (n_cells_ext > n_cells) {
1179 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
1180
1181 for (int isou = 0; isou < 6; isou++) {
1182 for (int jsou = 0; jsou < 6; jsou++) {
1183 da[cell_id][jsou][isou] = 0.;
1184 }
1185 }
1186
1187 }
1188 }
1189
1190 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1191 xa[face_id] = 0.;
1192 }
1193
1194 /* 2. Computation of extradiagonal terms */
1195
1196 for (cs_lnum_t face_id = 0; face_id <n_i_faces; face_id++) {
1197
1198 xa[face_id] = -thetap*idiffp*i_visc[face_id];
1199
1200 }
1201
1202 /* 3. Contribution of the extra-diagonal terms to the diagonal */
1203
1204 for (cs_lnum_t face_id = 0; face_id <n_i_faces; face_id++) {
1205
1206 cs_lnum_t ii = i_face_cells[face_id][0];
1207 cs_lnum_t jj = i_face_cells[face_id][1];
1208
1209 for (int isou = 0; isou < 6; isou++) {
1210 da[ii][isou][isou] -= xa[face_id];
1211 da[jj][isou][isou] -= xa[face_id];
1212 }
1213
1214 }
1215
1216 /* 4. Contribution of border faces to the diagonal */
1217
1218 for (cs_lnum_t face_id = 0; face_id <n_b_faces; face_id++) {
1219
1220 cs_lnum_t ii = b_face_cells[face_id];
1221
1222 for (int isou = 0; isou < 6; isou++) {
1223 for (int jsou = 0; jsou < 6; jsou++) {
1224 da[ii][jsou][isou] += thetap*idiffp*b_visc[face_id]
1225 *cofbfts[face_id][jsou][isou];
1226 }
1227 }
1228
1229 }
1230 }
1231
1232 /*----------------------------------------------------------------------------*/
1233 /*!
1234 * \brief Build the advection/diffusion matrix for a vector field
1235 * (non-symmetric matrix).
1236 *
1237 * The advection is upwind, the diffusion is not reconstructed.
1238 * The matrix is split into a diagonal block (3x3 times number of cells)
1239 * and an extra diagonal part (of dimension 2 time the number of internal
1240 * faces).
1241 *
1242 * \param[in] m pointer to mesh structure
1243 * \param[in] mq pointer to mesh quantities structure
1244 * \param[in] iconvp indicator
1245 * - 1 advection
1246 * - 0 otherwise
1247 * \param[in] idiffp indicator
1248 * - 1 diffusion
1249 * - 0 otherwise
1250 * \param[in] thetap weighting coefficient for the theta-scheme,
1251 * - thetap = 0: explicit scheme
1252 * - thetap = 0.5: time-centered
1253 * scheme (mix between Crank-Nicolson and
1254 * Adams-Bashforth)
1255 * - thetap = 1: implicit scheme
1256 * \param[in] coefbp boundary condition array for the variable
1257 * (implicit part - 3x3 tensor array)
1258 * \param[in] cofbfp boundary condition array for the variable flux
1259 * (implicit part - 3x3 tensor array)
1260 * \param[in] fimp part of the diagonal
1261 * \param[in] i_massflux mass flux at interior faces
1262 * \param[in] b_massflux mass flux at border faces
1263 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1264 * at interior faces for the matrix
1265 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1266 * at border faces for the matrix
1267 * \param[out] da diagonal part of the matrix
1268 * \param[out] xa extra interleaved diagonal part of the matrix
1269 */
1270 /*----------------------------------------------------------------------------*/
1271
1272 void
cs_matrix_vector(const cs_mesh_t * m,const cs_mesh_quantities_t * mq,int iconvp,int idiffp,cs_lnum_t eb_size[4],double thetap,const cs_real_33_t coefbp[],const cs_real_33_t cofbfp[],const cs_real_33_t fimp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_33_t * restrict da,cs_real_2_t * restrict xa)1273 cs_matrix_vector(const cs_mesh_t *m,
1274 const cs_mesh_quantities_t *mq,
1275 int iconvp,
1276 int idiffp,
1277 cs_lnum_t eb_size[4],
1278 double thetap,
1279 const cs_real_33_t coefbp[],
1280 const cs_real_33_t cofbfp[],
1281 const cs_real_33_t fimp[],
1282 const cs_real_t i_massflux[],
1283 const cs_real_t b_massflux[],
1284 const cs_real_t i_visc[],
1285 const cs_real_t b_visc[],
1286 cs_real_33_t *restrict da,
1287 cs_real_2_t *restrict xa)
1288 {
1289 const cs_lnum_t n_cells = m->n_cells;
1290 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1291 const cs_lnum_t n_i_faces = m->n_i_faces;
1292 const cs_lnum_t n_b_faces = m->n_b_faces;
1293
1294 const cs_lnum_2_t *restrict i_face_cells
1295 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1296 const cs_lnum_t *restrict b_face_cells
1297 = (const cs_lnum_t *restrict)m->b_face_cells;
1298
1299 const cs_real_2_t *i_f_face_factor;
1300 const cs_real_t *b_f_face_factor;
1301
1302 /* Discontinuous porous treatment */
1303 cs_real_2_t _i_f_face_factor = {1., 1.};
1304 cs_real_t _b_f_face_factor = 1.;
1305 int is_p = 0;
1306
1307 if (cs_glob_porous_model == 3) {
1308 i_f_face_factor = (const cs_real_2_t *)(mq->i_f_face_factor);
1309 b_f_face_factor = mq->b_f_face_factor;
1310 is_p = 1;
1311 }
1312 else {
1313 i_f_face_factor = &_i_f_face_factor;
1314 b_f_face_factor = &_b_f_face_factor;
1315 }
1316
1317 const cs_real_3_t *restrict i_face_normal
1318 = (const cs_real_3_t *restrict)mq->i_face_normal;
1319 const cs_real_3_t *restrict b_face_normal
1320 = (const cs_real_3_t *restrict)mq->b_face_normal;
1321
1322
1323 cs_real_332_t *_xa = (cs_real_332_t *) xa; //FIXME why 332 and use 233...
1324
1325 /* 1. Initialization */
1326
1327 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1328
1329 for (int isou = 0; isou < 3; isou++) {
1330 for (int jsou = 0; jsou < 3; jsou++) {
1331 da[cell_id][jsou][isou] = fimp[cell_id][jsou][isou];
1332 }
1333 }
1334
1335 }
1336
1337 if (n_cells_ext > n_cells) {
1338 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
1339
1340 for (int isou = 0; isou < 3; isou++) {
1341 for (int jsou = 0; jsou < 3; jsou++) {
1342 da[cell_id][jsou][isou] = 0.;
1343 }
1344 }
1345
1346 }
1347 }
1348
1349 if (eb_size[0] == 1) {
1350 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1351 xa[face_id][0] = 0.;
1352 xa[face_id][1] = 0.;
1353 }
1354 } else {
1355 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1356 for (cs_lnum_t i = 0; i < eb_size[0]; i++) {
1357 for (cs_lnum_t j = 0; j < eb_size[1]; j++) {
1358 _xa[face_id][0][i][j] = 0.;
1359 _xa[face_id][1][i][j] = 0.;
1360 }
1361 }
1362 }
1363 }
1364
1365 /* 2. Computation of extradiagonal terms */
1366
1367 if (eb_size[0] == 1) {
1368 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1369
1370 /*
1371 * X_ij = - theta f_j (m_ij)^-
1372 * X_ji = - theta f_i (m_ij)^+
1373 */
1374 cs_real_2_t flu = {
1375 0.5 * iconvp * (i_massflux[face_id] - fabs(i_massflux[face_id])),
1376 -0.5 * iconvp * (i_massflux[face_id] + fabs(i_massflux[face_id]))
1377 };
1378
1379 xa[face_id][0] = thetap*(flu[0] -idiffp*i_visc[face_id])
1380 * i_f_face_factor[is_p*face_id][1];//FIXME also diffusion? MF thinks so
1381 xa[face_id][1] = thetap*(flu[1] -idiffp*i_visc[face_id])
1382 * i_f_face_factor[is_p*face_id][0];
1383
1384 }
1385 }
1386 else {
1387 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1388
1389 /*
1390 * X_ij = - theta f_j (m_ij)^-
1391 * X_ji = - theta f_i (m_ij)^+
1392 */
1393 cs_real_2_t flu = {
1394 0.5 * iconvp * (i_massflux[face_id] - fabs(i_massflux[face_id]))
1395 - idiffp*i_visc[face_id],
1396 -0.5 * iconvp * (i_massflux[face_id] + fabs(i_massflux[face_id]))
1397 - idiffp*i_visc[face_id]
1398 };
1399
1400 cs_real_3_t normal;
1401 cs_math_3_normalise(i_face_normal[face_id], normal);
1402 /* Diagonal part:
1403 * the n(x)n term is multiplied by i_f_face_factor and (1 - n(x)n) by 1
1404 * XA_ij <= XA_ik n_k n_j (factor - 1) + XA_ij
1405 * XA_ij used to be diagonal: XA_ik n_k n_j = XA_ii n_i n_j*/
1406 for (cs_lnum_t i = 0; i < eb_size[0]; i++) {
1407 _xa[face_id][0][i][i] = flu[0];
1408 _xa[face_id][1][i][i] = flu[1];
1409 for (cs_lnum_t j = 0; j < eb_size[1]; j++) {
1410 _xa[face_id][0][i][j] = thetap*(_xa[face_id][0][i][j]
1411 + flu[0]
1412 * (i_f_face_factor[is_p*face_id][1] - 1.) * normal[i] * normal[j]);
1413 //FIXME also diffusion? MF thinks so
1414 _xa[face_id][1][i][j] = thetap*(_xa[face_id][1][i][j]
1415 + flu[1]
1416 * (i_f_face_factor[is_p*face_id][0] - 1.) * normal[i] * normal[j]);
1417 }
1418 }
1419
1420 }
1421 }
1422
1423 /* 3. Contribution of the extra-diagonal terms to the diagonal */
1424
1425 if (eb_size[0] == 1) {
1426 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1427
1428 cs_lnum_t ii = i_face_cells[face_id][0];
1429 cs_lnum_t jj = i_face_cells[face_id][1];
1430
1431 /* D_ii = theta f_i (m_ij)^+ - m_ij
1432 * = -X_ij - (1-theta)*m_ij
1433 * = -X_ji - m_ij
1434 * D_jj = -theta f_j (m_ij)^- + m_ij
1435 * = -X_ji + (1-theta)*m_ij
1436 * = -X_ij + m_ij
1437 */
1438 for (int i = 0; i < 3; i++) {
1439 da[ii][i][i] -= xa[face_id][1]
1440 + iconvp*i_massflux[face_id];
1441 da[jj][i][i] -= xa[face_id][0]
1442 - iconvp*i_massflux[face_id];
1443 }
1444 }
1445 }
1446 else {
1447 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1448
1449 cs_lnum_t ii = i_face_cells[face_id][0];
1450 cs_lnum_t jj = i_face_cells[face_id][1];
1451
1452 /* D_ii = theta (m_ij)^+ - m_ij
1453 * = -X_ij - (1-theta)*m_ij
1454 * = -X_ji - m_ij
1455 * D_jj = -theta (m_ij)^- + m_ij
1456 * = -X_ji + (1-theta)*m_ij
1457 * = -X_ij + m_ij
1458 */
1459 for (cs_lnum_t i = 0; i < eb_size[0]; i++) {
1460 da[ii][i][i] -= iconvp * i_massflux[face_id];
1461 da[jj][i][i] += iconvp * i_massflux[face_id];
1462
1463 for (cs_lnum_t j = 0; j < eb_size[1]; j++) {
1464 da[ii][i][j] -= _xa[face_id][1][i][j];
1465 da[jj][i][j] -= _xa[face_id][0][i][j];
1466 }
1467 }
1468
1469 }
1470 }
1471
1472 /* 4. Contribution of border faces to the diagonal */
1473
1474 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1475
1476 cs_lnum_t ii = b_face_cells[face_id];
1477
1478 cs_real_2_t flu = {
1479 /* (m_ij)^+ */
1480 iconvp * 0.5 * (b_massflux[face_id] + fabs(b_massflux[face_id])),
1481 /* (m_ij)^- */
1482 iconvp * 0.5 * (b_massflux[face_id] - fabs(b_massflux[face_id]))
1483 };
1484
1485 cs_real_3_t normal;
1486 cs_math_3_normalise(b_face_normal[face_id], normal);
1487 cs_real_t n_b_n
1488 = cs_math_3_33_3_dot_product(normal, coefbp[face_id], normal);
1489 cs_real_t n_bf_n
1490 = cs_math_3_33_3_dot_product(normal, cofbfp[face_id], normal);
1491
1492 for (int i = 0; i < 3; i++) {
1493 for (int j = 0; j < 3; j++) {
1494 cs_real_t d_ij = ((i == j) ? 1. : 0.);
1495 /* D = theta (m_f)^+.1 + theta B (m_f)^- - m_f.1
1496 * NB: stop here because the first two terms maybe scaled
1497 * = (theta B -1)*(m_f)^- - (1-theta)*(m_f)^+
1498 * = theta*(B -1)*(m_f)^- - (1-theta)*m_f
1499 */
1500 da[ii][i][j] +=
1501 thetap * (
1502 d_ij * flu[0]
1503 + flu[1] * coefbp[face_id][i][j]
1504 + idiffp * b_visc[face_id] * cofbfp[face_id][i][j]
1505 + (flu[0] + flu[1] * n_b_n + idiffp * b_visc[face_id] * n_bf_n)
1506 * (b_f_face_factor[is_p*face_id] - 1.) * normal[i] * normal[j]
1507 )
1508 - iconvp * d_ij * b_massflux[face_id];
1509 }
1510 }
1511
1512 }
1513 }
1514
1515 /*----------------------------------------------------------------------------*/
1516 /*!
1517 * \brief Build the advection/diffusion matrix for a tensor field
1518 * (non-symmetric matrix).
1519 *
1520 * The advection is upwind, the diffusion is not reconstructed.
1521 * The matrix is split into a diagonal block (6x6 times number of cells)
1522 * and an extra diagonal part (of dimension 2 time the number of internal
1523 * faces).
1524 *
1525 * \param[in] m pointer to mesh structure
1526 * \param[in] iconvp indicator
1527 * - 1 advection
1528 * - 0 otherwise
1529 * \param[in] idiffp indicator
1530 * - 1 diffusion
1531 * - 0 otherwise
1532 * \param[in] thetap weighting coefficient for the theta-scheme,
1533 * - thetap = 0: explicit scheme
1534 * - thetap = 0.5: time-centered
1535 * scheme (mix between Crank-Nicolson and
1536 * Adams-Bashforth)
1537 * - thetap = 1: implicit scheme
1538 * \param[in] coefbts boundary condition array for the variable
1539 * (Implicit part - 6x6 tensor array)
1540 * \param[in] cofbfts boundary condition array for the variable flux
1541 * (Implicit part - 6x6 tensor array)
1542 * \param[in] fimp part of the diagonal
1543 * \param[in] i_massflux mass flux at interior faces
1544 * \param[in] b_massflux mass flux at border faces
1545 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1546 * at interior faces for the matrix
1547 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1548 * at border faces for the matrix
1549 * \param[out] da diagonal part of the matrix
1550 * \param[out] xa extra interleaved diagonal part of the matrix
1551 */
1552 /*----------------------------------------------------------------------------*/
1553
1554 void
cs_matrix_tensor(const cs_mesh_t * m,int iconvp,int idiffp,double thetap,const cs_real_66_t coefbts[],const cs_real_66_t cofbfts[],const cs_real_66_t fimp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_66_t * restrict da,cs_real_2_t * restrict xa)1555 cs_matrix_tensor(const cs_mesh_t *m,
1556 int iconvp,
1557 int idiffp,
1558 double thetap,
1559 const cs_real_66_t coefbts[],
1560 const cs_real_66_t cofbfts[],
1561 const cs_real_66_t fimp[],
1562 const cs_real_t i_massflux[],
1563 const cs_real_t b_massflux[],
1564 const cs_real_t i_visc[],
1565 const cs_real_t b_visc[],
1566 cs_real_66_t *restrict da,
1567 cs_real_2_t *restrict xa)
1568 {
1569 const cs_lnum_t n_cells = m->n_cells;
1570 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1571 const cs_lnum_t n_i_faces = m->n_i_faces;
1572 const cs_lnum_t n_b_faces = m->n_b_faces;
1573
1574 const cs_lnum_2_t *restrict i_face_cells
1575 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1576 const cs_lnum_t *restrict b_face_cells
1577 = (const cs_lnum_t *restrict)m->b_face_cells;
1578
1579 /* 1. Initialization */
1580
1581 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1582
1583 for (int i = 0; i < 6; i++) {
1584 for (int j = 0; j < 6; j++) {
1585 da[cell_id][i][j] = fimp[cell_id][i][j];
1586 }
1587 }
1588
1589 }
1590
1591 if (n_cells_ext > n_cells) {
1592 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
1593
1594 for (int i = 0; i < 6; i++) {
1595 for (int j = 0; j < 6; j++) {
1596 da[cell_id][i][j] = 0.;
1597 }
1598 }
1599
1600 }
1601 }
1602
1603 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1604 xa[face_id][0] = 0.;
1605 xa[face_id][1] = 0.;
1606 }
1607
1608 /* 2. Computation of extradiagonal terms */
1609
1610 for (cs_lnum_t face_id = 0; face_id <n_i_faces; face_id++) {
1611
1612 double flui = 0.5*( i_massflux[face_id] -fabs(i_massflux[face_id]) );
1613 double fluj =-0.5*( i_massflux[face_id] +fabs(i_massflux[face_id]) );
1614
1615 xa[face_id][0] = thetap*(iconvp*flui -idiffp*i_visc[face_id]);
1616 xa[face_id][1] = thetap*(iconvp*fluj -idiffp*i_visc[face_id]);
1617
1618 }
1619
1620 /* 3. Contribution of the extra-diagonal terms to the diagonal */
1621
1622 for (cs_lnum_t face_id = 0; face_id <n_i_faces; face_id++) {
1623
1624 cs_lnum_t ii = i_face_cells[face_id][0];
1625 cs_lnum_t jj = i_face_cells[face_id][1];
1626
1627 /* D_ii = theta (m_ij)^+ - m_ij
1628 * = -X_ij - (1-theta)*m_ij
1629 * D_jj = -theta (m_ij)^- + m_ij
1630 * = -X_ji + (1-theta)*m_ij
1631 */
1632 for (int isou = 0; isou < 6; isou++) {
1633 da[ii][isou][isou] -= xa[face_id][0]
1634 + iconvp*(1. - thetap)*i_massflux[face_id];
1635 da[jj][isou][isou] -= xa[face_id][1]
1636 - iconvp*(1. - thetap)*i_massflux[face_id];
1637 }
1638
1639 }
1640
1641 /* 4. Contribution of border faces to the diagonal */
1642
1643 for (cs_lnum_t face_id = 0; face_id <n_b_faces; face_id++) {
1644
1645 cs_lnum_t ii = b_face_cells[face_id];
1646 double flui = 0.5*(b_massflux[face_id] -fabs(b_massflux[face_id]));
1647
1648 for (int isou = 0; isou < 6; isou++) {
1649 for (int jsou = 0; jsou < 6; jsou++) {
1650 /* D_ii = theta (m_f)^+ + theta B (m_f)^- - m_f
1651 * = (theta B -1)*(m_f)^- - (1-theta)*(m_f)^+
1652 * = theta*(B -1)*(m_f)^- - (1-theta)*m_f
1653 */
1654 if (isou == jsou) {
1655 da[ii][jsou][isou] += iconvp*( thetap*flui
1656 *(coefbts[face_id][jsou][isou]-1.)
1657 - (1. - thetap)*b_massflux[face_id])
1658 + idiffp*thetap*b_visc[face_id]
1659 *cofbfts[face_id][jsou][isou];
1660 } else {
1661 da[ii][jsou][isou] += thetap*( iconvp*flui*coefbts[face_id][jsou][isou]
1662 +idiffp*b_visc[face_id]
1663 *cofbfts[face_id][jsou][isou] );
1664 }
1665 }
1666 }
1667
1668 }
1669 }
1670
1671 /*----------------------------------------------------------------------------*/
1672 /*!
1673 * \brief Build the diagonal of the advection/diffusion matrix
1674 * for determining the variable time step, flow, Fourier.
1675 *
1676 * \param[in] m pointer to mesh structure
1677 * \param[in] iconvp indicator
1678 * - 1 advection
1679 * - 0 otherwise
1680 * \param[in] idiffp indicator
1681 * - 1 diffusion
1682 * - 0 otherwise
1683 * \param[in] isym indicator
1684 * - 1 symmetric matrix
1685 * - 2 non symmmetric matrix
1686 * \param[in] coefbp boundary condition array for the variable
1687 * (implicit part)
1688 * \param[in] cofbfp boundary condition array for the variable flux
1689 * (implicit part)
1690 * \param[in] i_massflux mass flux at interior faces
1691 * \param[in] b_massflux mass flux at border faces
1692 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1693 * at interior faces for the matrix
1694 * \param[in] b_visc \f$ S_\fib \f$
1695 * at border faces for the matrix
1696 * \param[out] da diagonal part of the matrix
1697 */
1698 /*----------------------------------------------------------------------------*/
1699
1700 void
cs_matrix_time_step(const cs_mesh_t * m,int iconvp,int idiffp,int isym,const cs_real_t coefbp[],const cs_real_t cofbfp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],cs_real_t * restrict da)1701 cs_matrix_time_step(const cs_mesh_t *m,
1702 int iconvp,
1703 int idiffp,
1704 int isym,
1705 const cs_real_t coefbp[],
1706 const cs_real_t cofbfp[],
1707 const cs_real_t i_massflux[],
1708 const cs_real_t b_massflux[],
1709 const cs_real_t i_visc[],
1710 const cs_real_t b_visc[],
1711 cs_real_t *restrict da)
1712 {
1713 const int n_cells = m->n_cells;
1714 const int n_cells_ext = m->n_cells_with_ghosts;
1715 const int n_i_groups = m->i_face_numbering->n_groups;
1716 const int n_i_threads = m->i_face_numbering->n_threads;
1717 const int n_b_groups = m->b_face_numbering->n_groups;
1718 const int n_b_threads = m->b_face_numbering->n_threads;
1719 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
1720 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
1721
1722 const cs_lnum_2_t *restrict i_face_cells
1723 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1724 const cs_lnum_t *restrict b_face_cells
1725 = (const cs_lnum_t *restrict)m->b_face_cells;
1726
1727 /* 1. Initialization */
1728
1729 if (isym != 1 && isym != 2) {
1730 bft_error(__FILE__, __LINE__, 0,
1731 _("invalid value of isym"));
1732 }
1733
1734 # pragma omp parallel for
1735 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1736 da[cell_id] = 0.;
1737 }
1738 if (n_cells_ext > n_cells) {
1739 # pragma omp parallel for if (n_cells_ext - n_cells > CS_THR_MIN)
1740 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
1741 da[cell_id] = 0.;
1742 }
1743 }
1744
1745 /* 2. Computation of extradiagonal terms unnecessary */
1746
1747 /* 3. Contribution of the extra-diagonal terms to the diagonal */
1748
1749 if (isym == 2) {
1750
1751 for (int g_id = 0; g_id < n_i_groups; g_id++) {
1752 # pragma omp parallel for
1753 for (int t_id = 0; t_id < n_i_threads; t_id++) {
1754 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
1755 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
1756 face_id++) {
1757
1758 cs_lnum_t ii = i_face_cells[face_id][0];
1759 cs_lnum_t jj = i_face_cells[face_id][1];
1760
1761 double fluj =-0.5*(i_massflux[face_id] + fabs(i_massflux[face_id]));
1762 double flui = 0.5*(i_massflux[face_id] - fabs(i_massflux[face_id]));
1763
1764 double xaifa2 = iconvp*fluj -idiffp*i_visc[face_id];
1765 double xaifa1 = iconvp*flui -idiffp*i_visc[face_id];
1766 da[ii] -= xaifa2;
1767 da[jj] -= xaifa1;
1768
1769 }
1770 }
1771 }
1772
1773 } else {
1774
1775 for (int g_id = 0; g_id < n_i_groups; g_id++) {
1776 # pragma omp parallel for
1777 for (int t_id = 0; t_id < n_i_threads; t_id++) {
1778 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
1779 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
1780 face_id++) {
1781
1782 cs_lnum_t ii = i_face_cells[face_id][0];
1783 cs_lnum_t jj = i_face_cells[face_id][1];
1784
1785 double flui = 0.5*(i_massflux[face_id] - fabs(i_massflux[face_id]));
1786
1787 double xaifa1 = iconvp*flui -idiffp*i_visc[face_id];
1788 da[ii] -= xaifa1;
1789 da[jj] -= xaifa1;
1790
1791 }
1792 }
1793 }
1794
1795 }
1796
1797 /* 4. Contribution of border faces to the diagonal */
1798
1799 for (int g_id = 0; g_id < n_b_groups; g_id++) {
1800 # pragma omp parallel for if (m->n_b_faces > CS_THR_MIN)
1801 for (int t_id = 0; t_id < n_b_threads; t_id++) {
1802 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
1803 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
1804 face_id++) {
1805
1806 cs_lnum_t ii = b_face_cells[face_id];
1807
1808 double flui = 0.5*(b_massflux[face_id] - fabs(b_massflux[face_id]));
1809 double fluj = -0.5*(b_massflux[face_id] + fabs(b_massflux[face_id]));
1810
1811 da[ii] += iconvp*(-fluj + flui*coefbp[face_id])
1812 + idiffp*b_visc[face_id]*cofbfp[face_id];
1813
1814 }
1815 }
1816 }
1817 }
1818
1819 /*----------------------------------------------------------------------------*/
1820 /*!
1821 * \brief Build the advection/diffusion matrix for a vector field with a
1822 * tensorial diffusivity.
1823 *
1824 * The advection is upwind, the diffusion is not reconstructed.
1825 * The matrix is split into a diagonal block (3x3 times number of cells)
1826 * and an extra diagonal part (of dimension 2 times 3x3 the number of internal
1827 * faces).
1828 *
1829 * \param[in] m pointer to mesh structure
1830 * \param[in] mq pointer to mesh quantities structure
1831 * \param[in] iconvp indicator
1832 * - 1 advection
1833 * - 0 otherwise
1834 * \param[in] idiffp indicator
1835 * - 1 diffusion
1836 * - 0 otherwise
1837 * \param[in] thetap weighting coefficient for the theta-scheme,
1838 * - thetap = 0: explicit scheme
1839 * - thetap = 0.5: time-centered
1840 * scheme (mix between Crank-Nicolson and
1841 * Adams-Bashforth)
1842 * - thetap = 1: implicit scheme
1843 * \param[in] coefbp boundary condition array for the variable
1844 * (implicit part - 3x3 tensor array)
1845 * \param[in] cofbfp boundary condition array for the variable flux
1846 * (implicit part - 3x3 tensor array)
1847 * \param[in] fimp part of the diagonal
1848 * \param[in] i_massflux mass flux at interior faces
1849 * \param[in] b_massflux mass flux at border faces
1850 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
1851 * at interior faces for the matrix
1852 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
1853 * at border faces for the matrix
1854 * \param[out] da diagonal part of the matrix
1855 * \param[out] xa extra interleaved diagonal part of the matrix
1856 */
1857 /*----------------------------------------------------------------------------*/
1858
1859 void
cs_matrix_anisotropic_diffusion(const cs_mesh_t * m,const cs_mesh_quantities_t * mq,int iconvp,int idiffp,double thetap,const cs_real_33_t coefbp[],const cs_real_33_t cofbfp[],const cs_real_33_t fimp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_33_t i_visc[],const cs_real_t b_visc[],cs_real_33_t * restrict da,cs_real_332_t * restrict xa)1860 cs_matrix_anisotropic_diffusion(const cs_mesh_t *m,
1861 const cs_mesh_quantities_t *mq,
1862 int iconvp,
1863 int idiffp,
1864 double thetap,
1865 const cs_real_33_t coefbp[],
1866 const cs_real_33_t cofbfp[],
1867 const cs_real_33_t fimp[],
1868 const cs_real_t i_massflux[],
1869 const cs_real_t b_massflux[],
1870 const cs_real_33_t i_visc[],
1871 const cs_real_t b_visc[],
1872 cs_real_33_t *restrict da,
1873 cs_real_332_t *restrict xa)
1874 {
1875 const cs_lnum_t n_cells = m->n_cells;
1876 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
1877 const cs_lnum_t n_i_faces = m->n_i_faces;
1878 const cs_lnum_t n_b_faces = m->n_b_faces;
1879
1880 const cs_lnum_2_t *restrict i_face_cells
1881 = (const cs_lnum_2_t *restrict)m->i_face_cells;
1882 const cs_lnum_t *restrict b_face_cells
1883 = (const cs_lnum_t *restrict)m->b_face_cells;
1884
1885 const cs_real_2_t *i_f_face_factor;
1886 const cs_real_t *b_f_face_factor;
1887
1888 const cs_real_3_t *restrict i_face_normal
1889 = (const cs_real_3_t *restrict)mq->i_face_normal;
1890
1891 const cs_real_3_t *restrict b_face_normal
1892 = (const cs_real_3_t *restrict)mq->b_face_normal;
1893
1894 /* Discontinuous porous treatment */
1895 cs_real_2_t _i_f_face_factor = {1., 1.};
1896 cs_real_t _b_f_face_factor = 1.;
1897 int is_p = 0;
1898
1899 /* Is it porous? */
1900 if (cs_glob_porous_model == 3) {
1901 i_f_face_factor = (const cs_real_2_t *)(mq->i_f_face_factor);
1902 b_f_face_factor = mq->b_f_face_factor;
1903 is_p = 1;
1904 }
1905 else {
1906 i_f_face_factor = &_i_f_face_factor;
1907 b_f_face_factor = &_b_f_face_factor;
1908 }
1909
1910 /* 1. Initialization */
1911
1912 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
1913 for (int i = 0; i < 3; i++) {
1914 for (int j = 0; j < 3; j++) {
1915 da[cell_id][j][i] = fimp[cell_id][j][i];
1916 }
1917 }
1918 }
1919 if (n_cells_ext > n_cells) {
1920 for (cs_lnum_t cell_id = n_cells; cell_id < n_cells_ext; cell_id++) {
1921 for (int i = 0; i < 3; i++) {
1922 for (int j = 0; j < 3; j++) {
1923 da[cell_id][j][i] = 0.;
1924 }
1925 }
1926 }
1927 }
1928
1929 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1930 for (int i = 0; i < 3; i++) {
1931 for (int j = 0; j < 3; j++) {
1932 xa[face_id][0][j][i] = 0.;
1933 xa[face_id][1][j][i] = 0.;
1934 }
1935 }
1936 }
1937
1938 /* 2. Computation of extradiagonal terms */
1939
1940 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1941
1942 cs_real_2_t flu = {
1943 iconvp * 0.5*(i_massflux[face_id] -fabs(i_massflux[face_id])),
1944 -iconvp * 0.5*(i_massflux[face_id] +fabs(i_massflux[face_id]))
1945 };
1946
1947 cs_real_3_t normal;
1948 cs_math_3_normalise(i_face_normal[face_id], normal);
1949
1950 for (int i = 0; i < 3; i++) {
1951 xa[face_id][0][i][i] = flu[0];
1952 xa[face_id][1][i][i] = flu[1];
1953 for (int j = 0; j < 3; j++) {
1954 xa[face_id][0][i][j] = thetap*( xa[face_id][0][i][j]
1955 + (i_f_face_factor[is_p*face_id][0] - 1.)
1956 * normal[i] * normal[j] * flu[0]
1957 - idiffp*i_visc[face_id][i][j]);//FIXME also diffusion? MF thinks so
1958 xa[face_id][1][i][j] = thetap*( xa[face_id][1][i][j]
1959 + (i_f_face_factor[is_p*face_id][1] - 1.)
1960 * normal[i] * normal[j] * flu[1]
1961 - idiffp*i_visc[face_id][i][j]);
1962 }
1963 }
1964
1965 }
1966
1967 /* 3. Contribution of the extra-diagonal terms to the diagonal */
1968
1969 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
1970
1971 cs_lnum_t ii = i_face_cells[face_id][0];
1972 cs_lnum_t jj = i_face_cells[face_id][1];
1973
1974 /* D_ii = theta (m_ij)^+ - m_ij
1975 * = -X_ij - (1-theta)*m_ij
1976 * D_jj = -theta (m_ij)^- + m_ij
1977 * = -X_ji + (1-theta)*m_ij
1978 */
1979 for (int i = 0; i < 3; i++) {
1980 da[ii][i][i] -= iconvp*(1. - thetap)*i_massflux[face_id];
1981 da[jj][i][i] += iconvp*(1. - thetap)*i_massflux[face_id];
1982 for (int j = 0; j < 3; j++) {
1983 da[ii][i][j] -= xa[face_id][0][i][j];
1984 da[jj][i][j] -= xa[face_id][1][i][j];
1985 }
1986 }
1987 }
1988
1989 /* 4. Contribution of border faces to the diagonal */
1990
1991 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
1992
1993 cs_lnum_t ii = b_face_cells[face_id];
1994
1995 cs_real_2_t flu = {
1996 /* (m_ij)^+ */
1997 iconvp * 0.5 * (b_massflux[face_id] + fabs(b_massflux[face_id])),
1998 /* (m_ij)^- */
1999 iconvp * 0.5 * (b_massflux[face_id] - fabs(b_massflux[face_id]))
2000 };
2001
2002 cs_real_3_t normal;
2003 cs_math_3_normalise(b_face_normal[face_id], normal);
2004 cs_real_t n_b_n
2005 = cs_math_3_33_3_dot_product(normal, coefbp[face_id], normal);
2006 cs_real_t n_bf_n
2007 = cs_math_3_33_3_dot_product(normal, cofbfp[face_id], normal);
2008
2009 for (int i = 0; i < 3; i++) {
2010 for (int j = 0; j < 3; j++) {
2011 cs_real_t d_ij = ((i == j) ? 1. : 0.);
2012 /* D = theta (m_f)^+.1 + theta B (m_f)^- - m_f.1
2013 * NB: stop here because the first two terms maybe scaled
2014 * = (theta B -1)*(m_f)^- - (1-theta)*(m_f)^+
2015 * = theta*(B -1)*(m_f)^- - (1-theta)*m_f
2016 */
2017 da[ii][i][j] +=
2018 thetap * (
2019 d_ij * flu[0]
2020 + flu[1] * coefbp[face_id][i][j]
2021 + idiffp * b_visc[face_id] * cofbfp[face_id][i][j]
2022 + (flu[0] + flu[1] * n_b_n + idiffp * b_visc[face_id] * n_bf_n)
2023 * (b_f_face_factor[is_p*face_id] - 1.) * normal[i] * normal[j]
2024 )
2025 - iconvp * d_ij * b_massflux[face_id];
2026 }
2027 }
2028
2029 }
2030 }
2031
2032 /*----------------------------------------------------------------------------*/
2033 /*!
2034 * \brief Build the advection/diffusion matrix for a tensor field with a
2035 * tensorial diffusivity.
2036 *
2037 * The advection is upwind, the diffusion is not reconstructed.
2038 * The matrix is split into a diagonal block (3x3 times number of cells)
2039 * and an extra diagonal part (of dimension 2 times 3x3 the number of internal
2040 * faces).
2041 *
2042 * \param[in] m pointer to mesh structure
2043 * \param[in] iconvp indicator
2044 * - 1 advection
2045 * - 0 otherwise
2046 * \param[in] idiffp indicator
2047 * - 1 diffusion
2048 * - 0 otherwise
2049 * \param[in] thetap weighting coefficient for the theta-scheme,
2050 * - thetap = 0: explicit scheme
2051 * - thetap = 0.5: time-centered
2052 * scheme (mix between Crank-Nicolson and
2053 * Adams-Bashforth)
2054 * - thetap = 1: implicit scheme
2055 * \param[in] coefbp boundary condition array for the variable
2056 * (implicit part - 3x3 tensor array)
2057 * \param[in] cofbfp boundary condition array for the variable flux
2058 * (implicit part - 3x3 tensor array)
2059 * \param[in] fimp part of the diagonal
2060 * \param[in] i_massflux mass flux at interior faces
2061 * \param[in] b_massflux mass flux at border faces
2062 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
2063 * at interior faces for the matrix
2064 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
2065 * at border faces for the matrix
2066 * \param[out] da diagonal part of the matrix
2067 * \param[out] xa extra interleaved diagonal part of the matrix
2068 */
2069 /*----------------------------------------------------------------------------*/
2070
2071 void
cs_matrix_anisotropic_diffusion_tensor(const cs_mesh_t * m,int iconvp,int idiffp,double thetap,const cs_real_66_t coefbp[],const cs_real_66_t cofbfp[],const cs_real_66_t fimp[],const cs_real_t i_massflux[],const cs_real_t b_massflux[],const cs_real_66_t i_visc[],const cs_real_t b_visc[],cs_real_66_t * restrict da,cs_real_662_t * restrict xa)2072 cs_matrix_anisotropic_diffusion_tensor(const cs_mesh_t *m,
2073 int iconvp,
2074 int idiffp,
2075 double thetap,
2076 const cs_real_66_t coefbp[],
2077 const cs_real_66_t cofbfp[],
2078 const cs_real_66_t fimp[],
2079 const cs_real_t i_massflux[],
2080 const cs_real_t b_massflux[],
2081 const cs_real_66_t i_visc[],
2082 const cs_real_t b_visc[],
2083 cs_real_66_t *restrict da,
2084 cs_real_662_t *restrict xa)
2085 {
2086 const cs_lnum_t n_cells = m->n_cells;
2087 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
2088 const cs_lnum_t n_i_faces = m->n_i_faces;
2089 const cs_lnum_t n_b_faces = m->n_b_faces;
2090
2091 const cs_lnum_2_t *restrict i_face_cells
2092 = (const cs_lnum_2_t *restrict)m->i_face_cells;
2093 const cs_lnum_t *restrict b_face_cells
2094 = (const cs_lnum_t *restrict)m->b_face_cells;
2095
2096 /* 1. Initialization */
2097
2098 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
2099 for (int j = 0; j < 6; j++) {
2100 for (int i = 0; i < 6; i++) {
2101 da[cell_id][i][j] = fimp[cell_id][i][j];
2102 }
2103 }
2104 }
2105 if (n_cells_ext > n_cells) {
2106 for (cs_lnum_t cell_id = n_cells+0; cell_id < n_cells_ext; cell_id++) {
2107 for (int j = 0; j < 6; j++) {
2108 for (int i = 0; i < 6; i++) {
2109 da[cell_id][i][j] = 0.;
2110 }
2111 }
2112 }
2113 }
2114
2115 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
2116 for (int j = 0; j < 6; j++) {
2117 for (int i = 0; i < 6; i++) {
2118 xa[face_id][0][j][i] = 0.;
2119 xa[face_id][1][j][i] = 0.;
2120 }
2121 }
2122 }
2123
2124 /* 2. Computation of extradiagonal terms */
2125
2126 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
2127
2128 double flui = 0.5*( i_massflux[face_id] -fabs(i_massflux[face_id]) );
2129 double fluj =-0.5*( i_massflux[face_id] +fabs(i_massflux[face_id]) );
2130
2131 for (int i = 0; i < 6; i++) {
2132 xa[face_id][0][i][i] = iconvp*flui;
2133 xa[face_id][1][i][i] = iconvp*fluj;
2134 for (int j = 0; j < 6; j++) {
2135 xa[face_id][0][i][j] = thetap*( xa[face_id][0][i][j]
2136 - idiffp*i_visc[face_id][i][j]);
2137 xa[face_id][1][i][j] = thetap*( xa[face_id][1][i][j]
2138 - idiffp*i_visc[face_id][i][j]);
2139 }
2140 }
2141
2142 }
2143
2144 /* 3. Contribution of the extra-diagonal terms to the diagonal */
2145
2146 for (cs_lnum_t face_id = 0; face_id < n_i_faces; face_id++) {
2147
2148 cs_lnum_t ii = i_face_cells[face_id][0];
2149 cs_lnum_t jj = i_face_cells[face_id][1];
2150
2151 /* D_ii = theta (m_ij)^+ - m_ij
2152 * = -X_ij - (1-theta)*m_ij
2153 * D_jj = -theta (m_ij)^- + m_ij
2154 * = -X_ji + (1-theta)*m_ij
2155 */
2156 for (int i = 0; i < 6; i++) {
2157 da[ii][i][i] -= iconvp*(1. - thetap)*i_massflux[face_id];
2158 da[jj][i][i] += iconvp*(1. - thetap)*i_massflux[face_id];
2159
2160 for (int j = 0; j < 6; j++) {
2161 da[ii][i][j] -= xa[face_id][0][i][j];
2162 da[jj][i][j] -= xa[face_id][1][i][j];
2163 }
2164 }
2165 }
2166
2167 /* 4. Contribution of border faces to the diagonal */
2168
2169 for (cs_lnum_t face_id = 0; face_id < n_b_faces; face_id++) {
2170
2171 cs_lnum_t ii = b_face_cells[face_id];
2172 double flui = 0.5*( b_massflux[face_id] -fabs(b_massflux[face_id]) );
2173
2174 for (int i = 0; i < 6; i++) {
2175 for (int j = 0; j < 6; j++) {
2176 /* D_ii = theta (m_f)^+ + theta B (m_f)^- - m_f
2177 * = (theta B -1)*(m_f)^- - (1-theta)*(m_f)^+
2178 * = theta*(B -1)*(m_f)^- - (1-theta)*m_f
2179 */
2180 if (j == i) {
2181 da[ii][i][j] += iconvp*( thetap*flui
2182 *(coefbp[face_id][i][j]-1.)
2183 - (1. - thetap)*b_massflux[face_id])
2184 + idiffp*thetap*b_visc[face_id]
2185 *cofbfp[face_id][i][j];
2186 } else {
2187 da[ii][i][j] += thetap*( iconvp*flui*coefbp[face_id][i][j]
2188 +idiffp*b_visc[face_id]
2189 *cofbfp[face_id][i][j] );
2190 }
2191 }
2192 }
2193
2194 }
2195 }
2196
2197 /*----------------------------------------------------------------------------*/
2198 /*!
2199 * \brief Build the diffusion matrix for a vector field with a
2200 * tensorial diffusivity (symmetric matrix).
2201 *
2202 * The diffusion is not reconstructed.
2203 * The matrix is split into a diagonal block (3x3 times number of cells)
2204 * and an extra diagonal part (of dimension 3x3 the number of internal
2205 * faces).
2206 *
2207 * \param[in] m pointer to mesh structure
2208 * \param[in] idiffp indicator
2209 * - 1 diffusion
2210 * - 0 otherwise
2211 * \param[in] thetap weighting coefficient for the theta-scheme,
2212 * - thetap = 0: explicit scheme
2213 * - thetap = 0.5: time-centered
2214 * scheme (mix between Crank-Nicolson and
2215 * Adams-Bashforth)
2216 * - thetap = 1: implicit scheme
2217 * \param[in] cofbfp boundary condition array for the variable flux
2218 * (implicit part - 3x3 tensor array)
2219 * \param[in] fimp \f$ \tens{f_s}^{imp} \f$
2220 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
2221 * at interior faces for the matrix
2222 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
2223 * at border faces for the matrix
2224 * \param[out] da diagonal part of the matrix
2225 * \param[out] xa extra interleaved diagonal part of the matrix
2226 */
2227 /*----------------------------------------------------------------------------*/
2228
2229 void
cs_sym_matrix_anisotropic_diffusion(const cs_mesh_t * m,int idiffp,double thetap,const cs_real_33_t cofbfp[],const cs_real_33_t fimp[],const cs_real_33_t i_visc[],const cs_real_t b_visc[],cs_real_33_t * restrict da,cs_real_33_t * restrict xa)2230 cs_sym_matrix_anisotropic_diffusion(const cs_mesh_t *m,
2231 int idiffp,
2232 double thetap,
2233 const cs_real_33_t cofbfp[],
2234 const cs_real_33_t fimp[],
2235 const cs_real_33_t i_visc[],
2236 const cs_real_t b_visc[],
2237 cs_real_33_t *restrict da,
2238 cs_real_33_t *restrict xa)
2239 {
2240 const cs_lnum_t n_cells = m->n_cells;
2241 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
2242 const cs_lnum_t n_i_faces = m->n_i_faces;
2243 const cs_lnum_t n_b_faces = m->n_b_faces;
2244
2245 const cs_lnum_2_t *restrict i_face_cells
2246 = (const cs_lnum_2_t *restrict)m->i_face_cells;
2247 const cs_lnum_t *restrict b_face_cells
2248 = (const cs_lnum_t *restrict)m->b_face_cells;
2249
2250 /* 1. Initialization */
2251
2252 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
2253 for (int isou = 0; isou < 3; isou++) {
2254 for (int jsou = 0; jsou < 3; jsou++) {
2255 da[cell_id][jsou][isou] = fimp[cell_id][jsou][isou];
2256 }
2257 }
2258 }
2259 if (n_cells_ext > n_cells) {
2260 for (cs_lnum_t cell_id = n_cells+0; cell_id < n_cells_ext; cell_id++) {
2261 for (int isou = 0; isou < 3; isou++) {
2262 for (int jsou = 0; jsou < 3; jsou++) {
2263 da[cell_id][jsou][isou] = 0.;
2264 }
2265 }
2266 }
2267 }
2268
2269 for (cs_lnum_t ifac = 0; ifac < n_i_faces; ifac++) {
2270 for (int isou = 0; isou < 3; isou++) {
2271 for (int jsou = 0; jsou < 3; jsou++) {
2272 xa[ifac][jsou][isou] = 0.;
2273 }
2274 }
2275 }
2276
2277 /* 2. Computation of extradiagonal terms */
2278
2279 for (cs_lnum_t ifac = 0; ifac < n_i_faces; ifac++) {
2280
2281 for (int isou = 0; isou < 3; isou++) {
2282 for (int jsou = 0; jsou < 3; jsou++) {
2283 xa[ifac][jsou][isou] = -thetap*idiffp*i_visc[ifac][jsou][isou];
2284 }
2285 }
2286
2287 }
2288
2289 /* 3. Contribution of the extra-diagonal terms to the diagonal */
2290
2291 for (cs_lnum_t ifac = 0; ifac < n_i_faces; ifac++) {
2292
2293 cs_lnum_t ii = i_face_cells[ifac][0];
2294 cs_lnum_t jj = i_face_cells[ifac][1];
2295
2296 for (int isou = 0; isou < 3; isou++) {
2297 for (int jsou = 0; jsou < 3; jsou++) {
2298 da[ii][jsou][isou] -= xa[ifac][jsou][isou];
2299 da[jj][jsou][isou] -= xa[ifac][jsou][isou];
2300 }
2301 }
2302 }
2303
2304 /* 4. Contribution of border faces to the diagonal */
2305
2306 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++) {
2307
2308 cs_lnum_t ii = b_face_cells[ifac];
2309
2310 for (int isou = 0; isou < 3; isou++) {
2311 for (int jsou = 0; jsou < 3; jsou++) {
2312 if (isou == jsou) {
2313 da[ii][jsou][isou] += thetap*idiffp*b_visc[ifac]
2314 *cofbfp[ifac][jsou][isou];
2315 }
2316 else {
2317 da[ii][jsou][isou] += thetap*idiffp*b_visc[ifac]
2318 *cofbfp[ifac][jsou][isou];
2319 }
2320 }
2321 }
2322
2323 }
2324 }
2325
2326 /*----------------------------------------------------------------------------*/
2327 /*!
2328 * \brief Build the diffusion matrix for a tensor field with a
2329 * tensorial diffusivity (symmetric matrix).
2330 *
2331 * The diffusion is not reconstructed.
2332 * The matrix is split into a diagonal block (3x3 times number of cells)
2333 * and an extra diagonal part (of dimension 3x3 the number of internal
2334 * faces).
2335 *
2336 * \param[in] m pointer to mesh structure
2337 * \param[in] idiffp indicator
2338 * - 1 diffusion
2339 * - 0 otherwise
2340 * \param[in] thetap weighting coefficient for the theta-scheme,
2341 * - thetap = 0: explicit scheme
2342 * - thetap = 0.5: time-centered
2343 * scheme (mix between Crank-Nicolson and
2344 * Adams-Bashforth)
2345 * - thetap = 1: implicit scheme
2346 * \param[in] cofbfp boundary condition array for the variable flux
2347 * (implicit part - 3x3 tensor array)
2348 * \param[in] fimp \f$ \tens{f_s}^{imp} \f$
2349 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
2350 * at interior faces for the matrix
2351 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
2352 * at border faces for the matrix
2353 * \param[out] da diagonal part of the matrix
2354 * \param[out] xa extra interleaved diagonal part of the matrix
2355 */
2356 /*----------------------------------------------------------------------------*/
2357
2358 void
cs_sym_matrix_anisotropic_diffusion_tensor(const cs_mesh_t * m,int idiffp,double thetap,const cs_real_66_t cofbfp[],const cs_real_66_t fimp[],const cs_real_66_t i_visc[],const cs_real_t b_visc[],cs_real_66_t * restrict da,cs_real_66_t * restrict xa)2359 cs_sym_matrix_anisotropic_diffusion_tensor(const cs_mesh_t *m,
2360 int idiffp,
2361 double thetap,
2362 const cs_real_66_t cofbfp[],
2363 const cs_real_66_t fimp[],
2364 const cs_real_66_t i_visc[],
2365 const cs_real_t b_visc[],
2366 cs_real_66_t *restrict da,
2367 cs_real_66_t *restrict xa)
2368 {
2369 const cs_lnum_t n_cells = m->n_cells;
2370 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
2371 const cs_lnum_t n_i_faces = m->n_i_faces;
2372 const cs_lnum_t n_b_faces = m->n_b_faces;
2373
2374 const cs_lnum_2_t *restrict i_face_cells
2375 = (const cs_lnum_2_t *restrict)m->i_face_cells;
2376 const cs_lnum_t *restrict b_face_cells
2377 = (const cs_lnum_t *restrict)m->b_face_cells;
2378
2379 /* 1. Initialization */
2380
2381 for (cs_lnum_t cell_id = 0; cell_id < n_cells; cell_id++) {
2382 for (int isou = 0; isou < 6; isou++) {
2383 for (int jsou = 0; jsou < 6; jsou++) {
2384 da[cell_id][jsou][isou] = fimp[cell_id][jsou][isou];
2385 }
2386 }
2387 }
2388 if (n_cells_ext > n_cells) {
2389 for (cs_lnum_t cell_id = n_cells+0; cell_id < n_cells_ext; cell_id++) {
2390 for (int isou = 0; isou < 6; isou++) {
2391 for (int jsou = 0; jsou < 6; jsou++) {
2392 da[cell_id][jsou][isou] = 0.;
2393 }
2394 }
2395 }
2396 }
2397
2398 for (cs_lnum_t ifac = 0; ifac < n_i_faces; ifac++) {
2399 for (int isou = 0; isou < 6; isou++) {
2400 for (int jsou = 0; jsou < 6; jsou++) {
2401 xa[ifac][jsou][isou] = 0.;
2402 }
2403 }
2404 }
2405
2406 /* 2. Computation of extradiagonal terms */
2407
2408 for (cs_lnum_t ifac = 0; ifac < n_i_faces; ifac++) {
2409
2410 for (int isou = 0; isou < 6; isou++) {
2411 for (int jsou = 0; jsou < 6; jsou++) {
2412 xa[ifac][jsou][isou] = -thetap*idiffp*i_visc[ifac][jsou][isou];
2413 }
2414 }
2415
2416 }
2417
2418 /* 3. Contribution of the extra-diagonal terms to the diagonal */
2419
2420 for (cs_lnum_t ifac = 0; ifac < n_i_faces; ifac++) {
2421
2422 cs_lnum_t ii = i_face_cells[ifac][0];
2423 cs_lnum_t jj = i_face_cells[ifac][1];
2424
2425 for (int isou = 0; isou < 6; isou++) {
2426 for (int jsou = 0; jsou < 6; jsou++) {
2427 da[ii][jsou][isou] -= xa[ifac][jsou][isou];
2428 da[jj][jsou][isou] -= xa[ifac][jsou][isou];
2429 }
2430 }
2431 }
2432
2433 /* 4. Contribution of border faces to the diagonal */
2434
2435 for (cs_lnum_t ifac = 0; ifac < n_b_faces; ifac++) {
2436
2437 cs_lnum_t ii = b_face_cells[ifac];
2438
2439 for (int isou = 0; isou < 6; isou++) {
2440 for (int jsou = 0; jsou < 6; jsou++) {
2441 da[ii][jsou][isou] += thetap*idiffp*b_visc[ifac]
2442 *cofbfp[ifac][jsou][isou];
2443 }
2444 }
2445
2446 }
2447 }
2448
2449 /*----------------------------------------------------------------------------*/
2450
2451 END_C_DECLS
2452