1 #ifndef __CS_CONVECTION_DIFFUSION_H__
2 #define __CS_CONVECTION_DIFFUSION_H__
3
4 /*============================================================================
5 * Convection-diffusion operators.
6 *============================================================================*/
7
8 /*
9 This file is part of Code_Saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-2021 EDF S.A.
12
13 This program is free software; you can redistribute it and/or modify it under
14 the terms of the GNU General Public License as published by the Free Software
15 Foundation; either version 2 of the License, or (at your option) any later
16 version.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21 details.
22
23 You should have received a copy of the GNU General Public License along with
24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25 Street, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27
28 /*----------------------------------------------------------------------------*/
29
30 #include "cs_defs.h"
31
32 /*----------------------------------------------------------------------------
33 * Local headers
34 *----------------------------------------------------------------------------*/
35
36 #include "cs_base.h"
37 #include "cs_halo.h"
38 #include "cs_math.h"
39 #include "cs_mesh_quantities.h"
40 #include "cs_parameters.h"
41 #include "cs_field_pointer.h"
42
43 /*----------------------------------------------------------------------------*/
44
45 BEGIN_C_DECLS
46
47 /*=============================================================================
48 * Macro definitions
49 *============================================================================*/
50
51 /*============================================================================
52 * Type definition
53 *============================================================================*/
54
55 /*----------------------------------------------------------------------------
56 * NVD/TVD Advection Scheme
57 *----------------------------------------------------------------------------*/
58
59 typedef enum {
60
61 CS_NVD_GAMMA = 0, /* GAMMA */
62 CS_NVD_SMART = 1, /* SMART */
63 CS_NVD_CUBISTA = 2, /* CUBISTA */
64 CS_NVD_SUPERBEE = 3, /* SUPERBEE */
65 CS_NVD_MUSCL = 4, /* MUSCL */
66 CS_NVD_MINMOD = 5, /* MINMOD */
67 CS_NVD_CLAM = 6, /* CLAM */
68 CS_NVD_STOIC = 7, /* STOIC */
69 CS_NVD_OSHER = 8, /* OSHER */
70 CS_NVD_WASEB = 9, /* WASEB */
71 CS_NVD_VOF_HRIC = 10, /* M-HRIC for VOF */
72 CS_NVD_VOF_CICSAM = 11, /* M-CICSAM for VOF */
73 CS_NVD_VOF_STACS = 12, /* STACS for VOF */
74 CS_NVD_N_TYPES = 13 /* number of NVD schemes */
75
76 } cs_nvd_type_t;
77
78 /*============================================================================
79 * Global variables
80 *============================================================================*/
81
82 /*============================================================================
83 * Public inlined function
84 *============================================================================*/
85
86 /*----------------------------------------------------------------------------
87 * Synchronize halos for scalar variables.
88 *
89 * parameters:
90 * m <-- pointer to associated mesh structure
91 * pvar <-> variable
92 *----------------------------------------------------------------------------*/
93
94 inline static void
cs_sync_scalar_halo(const cs_mesh_t * m,cs_real_t pvar[])95 cs_sync_scalar_halo(const cs_mesh_t *m,
96 cs_real_t pvar[])
97 {
98 if (m->halo != NULL)
99 cs_halo_sync_var(m->halo, CS_HALO_STANDARD, pvar);
100 }
101
102 /*----------------------------------------------------------------------------
103 * Return pointer to slope test indicator field values if active.
104 *
105 * parameters:
106 * f_id <-- field id (or -1)
107 * var_cal_opt <-- variable calculation options
108 *
109 * return:
110 * pointer to local values array, or NULL;
111 *----------------------------------------------------------------------------*/
112
113 inline static cs_real_t *
cs_get_v_slope_test(int f_id,const cs_var_cal_opt_t var_cal_opt)114 cs_get_v_slope_test(int f_id,
115 const cs_var_cal_opt_t var_cal_opt)
116 {
117 const int iconvp = var_cal_opt.iconv;
118 const int isstpp = var_cal_opt.isstpc;
119 const double blencp = var_cal_opt.blencv;
120
121 cs_real_t *v_slope_test = NULL;
122
123 if (f_id > -1 && iconvp > 0 && blencp > 0. && isstpp == 0) {
124
125 static int _k_slope_test_f_id = -1;
126
127 cs_field_t *f = cs_field_by_id(f_id);
128
129 int f_track_slope_test_id = -1;
130
131 if (_k_slope_test_f_id < 0)
132 _k_slope_test_f_id = cs_field_key_id_try("slope_test_upwind_id");
133 if (_k_slope_test_f_id > -1 && isstpp == 0)
134 f_track_slope_test_id = cs_field_get_key_int(f, _k_slope_test_f_id);
135
136 if (f_track_slope_test_id > -1)
137 v_slope_test = (cs_field_by_id(f_track_slope_test_id))->val;
138
139 if (v_slope_test != NULL) {
140 const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
141 # pragma omp parallel for
142 for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
143 v_slope_test[cell_id] = 0.;
144 }
145
146 }
147
148 return v_slope_test;
149 }
150
151 /*----------------------------------------------------------------------------
152 * Compute the local cell Courant number as the maximum of all cell face based
153 * Courant number at each cell.
154 *
155 * parameters:
156 * f_id <-- field id (or -1)
157 * courant --> cell Courant number
158 */
159 /*----------------------------------------------------------------------------*/
160
161 inline static void
cs_cell_courant_number(const int f_id,cs_real_t * courant)162 cs_cell_courant_number(const int f_id,
163 cs_real_t *courant)
164 {
165 const cs_mesh_t *m = cs_glob_mesh;
166 const cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
167
168 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
169 const int n_i_groups = m->i_face_numbering->n_groups;
170 const int n_i_threads = m->i_face_numbering->n_threads;
171 const int n_b_groups = m->b_face_numbering->n_groups;
172 const int n_b_threads = m->b_face_numbering->n_threads;
173 const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
174 const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
175
176 const cs_lnum_2_t *restrict i_face_cells
177 = (const cs_lnum_2_t *restrict)m->i_face_cells;
178 const cs_lnum_t *restrict b_face_cells
179 = (const cs_lnum_t *restrict)m->b_face_cells;
180
181 const cs_real_t *restrict vol
182 = (cs_real_t *restrict)fvq->cell_vol;
183
184 cs_field_t *f = cs_field_by_id(f_id);
185 const int kimasf = cs_field_key_id("inner_mass_flux_id");
186 const int kbmasf = cs_field_key_id("boundary_mass_flux_id");
187 const cs_real_t *restrict i_massflux
188 = cs_field_by_id( cs_field_get_key_int(f, kimasf) )->val;
189 const cs_real_t *restrict b_massflux
190 = cs_field_by_id( cs_field_get_key_int(f, kbmasf) )->val;
191
192 const cs_real_t *restrict dt
193 = (const cs_real_t *restrict)CS_F_(dt)->val;
194
195 /* Initialisation */
196
197 # pragma omp parallel for
198 for (cs_lnum_t ii = 0; ii < n_cells_ext; ii++) {
199 courant[ii] = 0.;
200 }
201
202 /* ---> Contribution from interior faces */
203
204 cs_real_t cnt;
205
206 for (int g_id = 0; g_id < n_i_groups; g_id++) {
207 # pragma omp parallel for
208 for (int t_id = 0; t_id < n_i_threads; t_id++) {
209 for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
210 face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
211 face_id++) {
212 cs_lnum_t ii = i_face_cells[face_id][0];
213 cs_lnum_t jj = i_face_cells[face_id][1];
214
215 cnt = CS_ABS(i_massflux[face_id])*dt[ii]/vol[ii];
216 courant[ii] = CS_MAX(courant[ii], cnt);
217
218 cnt = CS_ABS(i_massflux[face_id])*dt[jj]/vol[jj];
219 courant[jj] = CS_MAX(courant[jj], cnt);
220 }
221 }
222 }
223
224 /* ---> Contribution from boundary faces */
225
226 for (int g_id = 0; g_id < n_b_groups; g_id++) {
227 # pragma omp parallel for
228 for (int t_id = 0; t_id < n_b_threads; t_id++) {
229 for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
230 face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
231 face_id++) {
232 cs_lnum_t ii = b_face_cells[face_id];
233
234 cnt = CS_ABS(b_massflux[face_id])*dt[ii]/vol[ii];
235 courant[ii] = CS_MAX(courant[ii], cnt);
236 }
237 }
238 }
239 }
240
241 /*----------------------------------------------------------------------------*/
242 /*!
243 * \brief Compute the normalised face scalar using the specified NVD scheme.
244 *
245 * \param[in] scheme choice of the NVD scheme
246 * \param[in] nvf_p_c normalised property of the current cell
247 * \param[in] nvf_r_f normalised distance from the face
248 * \param[in] nvf_r_c normalised distance from the current cell
249 *
250 * \return normalised face scalar value
251 */
252 /*----------------------------------------------------------------------------*/
253
254 inline static cs_real_t
cs_nvd_scheme_scalar(const cs_nvd_type_t limiter,const cs_real_t nvf_p_c,const cs_real_t nvf_r_f,const cs_real_t nvf_r_c)255 cs_nvd_scheme_scalar(const cs_nvd_type_t limiter,
256 const cs_real_t nvf_p_c,
257 const cs_real_t nvf_r_f,
258 const cs_real_t nvf_r_c)
259 {
260 cs_real_t nvf_p_f;
261
262 cs_real_t beta_m, rfc, r1f, r1, r2, r3, b1, b2;
263
264 switch (limiter) {
265 case CS_NVD_GAMMA: /* Gamma scheme */
266 beta_m = 0.1; /* in [0.1, 0.5] */
267 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
268
269 if (nvf_p_c < beta_m) {
270 nvf_p_f = nvf_p_c*(1.+rfc*(1.-nvf_p_c)/beta_m);
271 } else {
272 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
273
274 nvf_p_f = r1f*nvf_p_c+rfc;
275 }
276
277 break;
278
279 case CS_NVD_SMART: /* SMART scheme */
280 if (nvf_p_c < (nvf_r_c/3.)) {
281 r1 = nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
282 r2 = nvf_r_c*(1.-nvf_r_c);
283
284 nvf_p_f = nvf_p_c*r1/r2;
285 } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
286 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
287 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
288
289 nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c + rfc);
290 } else {
291 nvf_p_f = 1.;
292 }
293
294 break;
295
296 case CS_NVD_CUBISTA: /* CUBISTA scheme */
297 if (nvf_p_c < (3.*nvf_r_c/4.)) {
298 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
299
300 nvf_p_f = nvf_r_f*(1.+rfc/3.)*nvf_p_c/nvf_r_c;
301 } else if (nvf_p_c <= (nvf_r_c*(1.+2.*(nvf_r_f-nvf_r_c))/(2.*nvf_r_f-nvf_r_c))) {
302 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
303 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
304
305 nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c+rfc);
306 } else {
307 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
308
309 nvf_p_f = 1.-.5*r1f*(1.-nvf_p_c);
310 }
311
312 break;
313
314 case CS_NVD_SUPERBEE: /* SuperBee scheme */
315 if (nvf_p_c < (nvf_r_c/(2.-nvf_r_c))) {
316 nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
317 } else if (nvf_p_c < nvf_r_c) {
318 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
319 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
320
321 nvf_p_f = r1f*nvf_p_c+rfc;
322 } else if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
323 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
324 } else {
325 nvf_p_f = 1.;
326 }
327
328 break;
329
330 case CS_NVD_MUSCL: /* MUSCL scheme */
331 if (nvf_p_c < (.5*nvf_r_c)) {
332 nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
333 } else if (nvf_p_c < (1.+nvf_r_c-nvf_r_f)) {
334 nvf_p_f = nvf_p_c+nvf_r_f-nvf_r_c;
335 } else {
336 nvf_p_f = 1.;
337 }
338
339 break;
340
341 case CS_NVD_MINMOD: /* MINMOD scheme */
342 if (nvf_p_c < nvf_r_c) {
343 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
344 } else {
345 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
346 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
347
348 nvf_p_f = r1f*nvf_p_c+rfc;
349 }
350
351 break;
352
353 case CS_NVD_CLAM: /* CLAM scheme */
354 r1 = nvf_r_c*nvf_r_c-nvf_r_f;
355 r2 = nvf_r_c*(nvf_r_c-1.);
356 r3 = nvf_r_f-nvf_r_c;
357
358 nvf_p_f = nvf_p_c*(r1+r3*nvf_p_c)/r2;
359
360 break;
361
362 case CS_NVD_STOIC: /* STOIC scheme */
363 b1 = (nvf_r_c-nvf_r_f)*nvf_r_c;
364 b2 = nvf_r_c+nvf_r_f+2.*nvf_r_f*nvf_r_f-4.*nvf_r_f*nvf_r_c;
365
366 if (nvf_p_c < (b1/b2)) {
367 r1 = -nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
368 r2 = nvf_r_c*(nvf_r_c-1.);
369
370 nvf_p_f = nvf_p_c*r1/r2;
371 } else if (nvf_p_c < nvf_r_c) {
372 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
373 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
374
375 nvf_p_f = r1f*nvf_p_c+rfc;
376 } else if (nvf_p_c < (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
377 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
378 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
379
380 nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
381 } else {
382 nvf_p_f = 1.;
383 }
384
385 break;
386
387 case CS_NVD_OSHER: /* OSHER scheme */
388 if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
389 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
390 } else {
391 nvf_p_f = 1.;
392 }
393
394 break;
395
396 case CS_NVD_WASEB: /* WASEB scheme */
397 r1 = nvf_r_c*nvf_r_f*(nvf_r_f-nvf_r_c);
398 r2 = 2.*nvf_r_c*(1.-nvf_r_c)-nvf_r_f*(1.-nvf_r_f);
399
400 if (nvf_p_c < (r1/r2)) {
401 nvf_p_f = 2.*nvf_p_c;
402 } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
403 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
404 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
405
406 nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
407 } else {
408 nvf_p_f = 1.;
409 }
410
411 break;
412
413 default: /* Upwinding */
414 nvf_p_f = nvf_p_c;
415 break;
416 }
417
418 return nvf_p_f;
419 }
420
421 /*----------------------------------------------------------------------------*/
422 /*!
423 * \brief Compute the normalised face scalar using the specified NVD scheme
424 * for the case of a Volume-of-Fluid (VOF) transport equation.
425 *
426 * \param[in] limiter choice of the NVD scheme
427 * \param[in] i_face_normal normal of face ij
428 * \param[in] face_id the current cell face
429 * \param[in] nvf_p_c normalised property of the current cell
430 * \param[in] nvf_r_f normalised distance from the face
431 * \param[in] nvf_r_c normalised distance from the current cell
432 * \param[in] gradv_c gradient at central cell
433 * \param[in] courant_c courant at central cell
434 *
435 * \return normalised face scalar
436 */
437 /*----------------------------------------------------------------------------*/
438
439 inline static cs_real_t
cs_nvd_vof_scheme_scalar(const cs_nvd_type_t limiter,const cs_real_3_t i_face_normal,const cs_real_t nvf_p_c,const cs_real_t nvf_r_f,const cs_real_t nvf_r_c,const cs_real_3_t gradv_c,const cs_real_t c_courant)440 cs_nvd_vof_scheme_scalar(const cs_nvd_type_t limiter,
441 const cs_real_3_t i_face_normal,
442 const cs_real_t nvf_p_c,
443 const cs_real_t nvf_r_f,
444 const cs_real_t nvf_r_c,
445 const cs_real_3_t gradv_c,
446 const cs_real_t c_courant)
447 {
448 assert(limiter >= CS_NVD_VOF_HRIC);
449
450 cs_real_t nvf_p_f;
451 cs_real_t blend, high_order, low_order, ratio;
452
453 /* Compute gradient angle indicator */
454 cs_real_t dotp = CS_ABS(cs_math_3_dot_product(gradv_c, i_face_normal));
455 cs_real_t sgrad = cs_math_3_norm(gradv_c);
456 cs_real_t snorm = cs_math_3_norm(i_face_normal);
457 cs_real_t denom = snorm*sgrad;
458
459 if (limiter == CS_NVD_VOF_HRIC) { /* M-HRIC scheme */
460 /* High order scheme : Bounded Downwind */
461 if (nvf_p_c <= .5) {
462 high_order = 2.*nvf_p_c;
463 } else {
464 high_order = 1.;
465 }
466
467 /* Low order scheme : MUSCL */
468 low_order = cs_nvd_scheme_scalar(CS_NVD_MUSCL,
469 nvf_p_c,
470 nvf_r_f,
471 nvf_r_c);
472
473 /* Compute the blending factor */
474 if (denom <= (cs_math_epzero*dotp)) {
475 blend = 1.;
476 } else {
477 ratio = dotp/denom;
478 blend = CS_MIN(1., pow(ratio, .5));
479 }
480
481 /* Blending */
482 nvf_p_f = blend*high_order + (1.-blend)*low_order;
483
484 /* Extra blending due to the cell Courant number */
485 if (c_courant < .7 && c_courant > .3) {
486 nvf_p_f = nvf_p_f + (nvf_p_f - low_order)*(.7 - c_courant )/.4;
487 } else if (c_courant >= .7) {
488 nvf_p_f = low_order;
489 }
490 } else if (limiter == CS_NVD_VOF_CICSAM) { /* M-CICSAM scheme */
491 /* High order scheme : HYPER-C + SUPERBEE */
492 if (c_courant <= .3) {
493 high_order = CS_MIN(1., nvf_p_c/(c_courant+cs_math_epzero));
494 } else if (c_courant <= .6) {
495 high_order = CS_MIN(1., nvf_p_c/.3);
496 } else if (c_courant <= .7) {
497 cs_real_t superbee = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
498 nvf_p_c,
499 nvf_r_f,
500 nvf_r_c);
501 high_order = 10.*( (.7-c_courant)*CS_MIN(1., nvf_p_c/.3)
502 + (c_courant-.6)*superbee);
503 }
504 else {
505 high_order = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
506 nvf_p_c,
507 nvf_r_f,
508 nvf_r_c);
509 }
510
511 /* Low order scheme : MUSCL */
512 low_order = cs_nvd_scheme_scalar(CS_NVD_MUSCL,
513 nvf_p_c,
514 nvf_r_f,
515 nvf_r_c);
516
517 /* Compute the blending factor */
518 if (denom <= (cs_math_epzero*dotp)) {
519 blend = 1.;
520 } else {
521 ratio = dotp/denom;
522 blend = CS_MIN(1., pow(ratio, 2.));
523 }
524
525 /* Blending */
526 nvf_p_f = blend*high_order + (1.-blend)*low_order;
527 } else { /* STACS scheme */
528 /* High order scheme : SUPERBEE */
529 high_order = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
530 nvf_p_c,
531 nvf_r_f,
532 nvf_r_c);
533
534 /* Low order scheme : STOIC */
535 low_order = cs_nvd_scheme_scalar(CS_NVD_STOIC,
536 nvf_p_c,
537 nvf_r_f,
538 nvf_r_c);
539
540 /* Compute the blending factor */
541 if (denom <= (cs_math_epzero*dotp)) {
542 blend = 1.;
543 } else {
544 ratio = dotp/denom;
545 blend = CS_MIN(1., pow(ratio, 4.));
546 }
547
548 /* Blending */
549 nvf_p_f = blend*high_order + (1.-blend)*low_order;
550 }
551
552 return nvf_p_f;
553 }
554
555 /*----------------------------------------------------------------------------*/
556 /*!
557 * \brief Compute slope test criteria at internal face between cell i and j.
558 *
559 * \param[in] pi value at cell i
560 * \param[in] pj value at cell j
561 * \param[in] distf distance IJ.Nij
562 * \param[in] srfan face surface
563 * \param[in] i_face_normal face normal
564 * \param[in] gradi gradient at cell i
565 * \param[in] gradj gradient at cell j
566 * \param[in] grdpai upwind gradient at cell i
567 * \param[in] grdpaj upwind gradient at cell j
568 * \param[in] i_massflux mass flux at face (from i to j)
569 * \param[out] testij value of slope test first criterion
570 * \param[out] tesqck value of slope test second criterion
571 */
572 /*----------------------------------------------------------------------------*/
573
574 inline static void
cs_slope_test(const cs_real_t pi,const cs_real_t pj,const cs_real_t distf,const cs_real_t srfan,const cs_real_t i_face_normal[3],const cs_real_t gradi[3],const cs_real_t gradj[3],const cs_real_t grdpai[3],const cs_real_t grdpaj[3],const cs_real_t i_massflux,cs_real_t * testij,cs_real_t * tesqck)575 cs_slope_test(const cs_real_t pi,
576 const cs_real_t pj,
577 const cs_real_t distf,
578 const cs_real_t srfan,
579 const cs_real_t i_face_normal[3],
580 const cs_real_t gradi[3],
581 const cs_real_t gradj[3],
582 const cs_real_t grdpai[3],
583 const cs_real_t grdpaj[3],
584 const cs_real_t i_massflux,
585 cs_real_t *testij,
586 cs_real_t *tesqck)
587 {
588 cs_real_t testi, testj;
589 cs_real_t dcc, ddi, ddj;
590
591 /* Slope test */
592
593 testi = grdpai[0]*i_face_normal[0]
594 + grdpai[1]*i_face_normal[1]
595 + grdpai[2]*i_face_normal[2];
596 testj = grdpaj[0]*i_face_normal[0]
597 + grdpaj[1]*i_face_normal[1]
598 + grdpaj[2]*i_face_normal[2];
599
600 *testij = grdpai[0]*grdpaj[0]
601 + grdpai[1]*grdpaj[1]
602 + grdpai[2]*grdpaj[2];
603
604 if (i_massflux>0.) {
605 dcc = gradi[0]*i_face_normal[0]
606 + gradi[1]*i_face_normal[1]
607 + gradi[2]*i_face_normal[2];
608 ddi = testi;
609 ddj = (pj-pi)/distf *srfan;
610 } else {
611 dcc = gradj[0]*i_face_normal[0]
612 + gradj[1]*i_face_normal[1]
613 + gradj[2]*i_face_normal[2];
614 ddi = (pj-pi)/distf *srfan;
615 ddj = testj;
616 }
617
618 *tesqck = cs_math_sq(dcc) - cs_math_sq(ddi-ddj);
619 }
620
621 /*----------------------------------------------------------------------------*/
622 /*!
623 * \brief Compute slope test criteria at internal face between cell i and j.
624 *
625 * \param[in] pi value at cell i
626 * \param[in] pj value at cell j
627 * \param[in] distf distance IJ.Nij
628 * \param[in] srfan face surface
629 * \param[in] i_face_normal face normal
630 * \param[in] gradi gradient at cell i
631 * \param[in] gradj gradient at cell j
632 * \param[in] grdpai upwind gradient at cell i
633 * \param[in] grdpaj upwind gradient at cell j
634 * \param[in] i_massflux mass flux at face (from i to j)
635 * \param[out] testij value of slope test first criterion
636 * \param[out] tesqck value of slope test second criterion
637 */
638 /*----------------------------------------------------------------------------*/
639
640 inline static void
cs_slope_test_vector(const cs_real_t pi[3],const cs_real_t pj[3],const cs_real_t distf,const cs_real_t srfan,const cs_real_t i_face_normal[3],const cs_real_t gradi[3][3],const cs_real_t gradj[3][3],const cs_real_t gradsti[3][3],const cs_real_t gradstj[3][3],const cs_real_t i_massflux,cs_real_t * testij,cs_real_t * tesqck)641 cs_slope_test_vector(const cs_real_t pi[3],
642 const cs_real_t pj[3],
643 const cs_real_t distf,
644 const cs_real_t srfan,
645 const cs_real_t i_face_normal[3],
646 const cs_real_t gradi[3][3],
647 const cs_real_t gradj[3][3],
648 const cs_real_t gradsti[3][3],
649 const cs_real_t gradstj[3][3],
650 const cs_real_t i_massflux,
651 cs_real_t *testij,
652 cs_real_t *tesqck)
653 {
654 cs_real_t testi[3], testj[3];
655 cs_real_t dcc[3], ddi[3], ddj[3];
656 *testij = 0.;
657 *tesqck = 0.;
658
659 /* Slope test */
660
661 for (int i = 0; i < 3; i++) {
662 *testij += gradsti[i][0]*gradstj[i][0]
663 + gradsti[i][1]*gradstj[i][1]
664 + gradsti[i][2]*gradstj[i][2];
665
666 testi[i] = gradsti[i][0]*i_face_normal[0]
667 + gradsti[i][1]*i_face_normal[1]
668 + gradsti[i][2]*i_face_normal[2];
669 testj[i] = gradstj[i][0]*i_face_normal[0]
670 + gradstj[i][1]*i_face_normal[1]
671 + gradstj[i][2]*i_face_normal[2];
672
673 if (i_massflux > 0.) {
674 dcc[i] = gradi[i][0]*i_face_normal[0]
675 + gradi[i][1]*i_face_normal[1]
676 + gradi[i][2]*i_face_normal[2];
677 ddi[i] = testi[i];
678 ddj[i] = (pj[i]-pi[i])/distf *srfan;
679 } else {
680 dcc[i] = gradj[i][0]*i_face_normal[0]
681 + gradj[i][1]*i_face_normal[1]
682 + gradj[i][2]*i_face_normal[2];
683 ddi[i] = (pj[i]-pi[i])/distf *srfan;
684 ddj[i] = testj[i];
685 }
686 }
687
688 *tesqck = cs_math_3_square_norm(dcc) - cs_math_3_square_distance(ddi, ddj);
689 }
690
691 /*----------------------------------------------------------------------------*/
692 /*!
693 * \brief Compute slope test criteria at internal face between cell i and j.
694 *
695 * \param[in] pi value at cell i
696 * \param[in] pj value at cell j
697 * \param[in] distf distance IJ.Nij
698 * \param[in] srfan face surface
699 * \param[in] i_face_normal face normal
700 * \param[in] gradi gradient at cell i
701 * \param[in] gradj gradient at cell j
702 * \param[in] grdpai upwind gradient at cell i
703 * \param[in] grdpaj upwind gradient at cell j
704 * \param[in] i_massflux mass flux at face (from i to j)
705 * \param[out] testij value of slope test first criterion
706 * \param[out] tesqck value of slope test second criterion
707 */
708 /*----------------------------------------------------------------------------*/
709
710 inline static void
cs_slope_test_tensor(const cs_real_t pi[6],const cs_real_t pj[6],const cs_real_t distf,const cs_real_t srfan,const cs_real_t i_face_normal[3],const cs_real_t gradi[6][3],const cs_real_t gradj[6][3],const cs_real_t gradsti[6][3],const cs_real_t gradstj[6][3],const cs_real_t i_massflux,cs_real_t * testij,cs_real_t * tesqck)711 cs_slope_test_tensor(const cs_real_t pi[6],
712 const cs_real_t pj[6],
713 const cs_real_t distf,
714 const cs_real_t srfan,
715 const cs_real_t i_face_normal[3],
716 const cs_real_t gradi[6][3],
717 const cs_real_t gradj[6][3],
718 const cs_real_t gradsti[6][3],
719 const cs_real_t gradstj[6][3],
720 const cs_real_t i_massflux,
721 cs_real_t *testij,
722 cs_real_t *tesqck)
723 {
724 cs_real_t testi[6], testj[6];
725 cs_real_t dcc[6], ddi[6], ddj[6];
726
727 *testij = 0.;
728 *tesqck = 0.;
729
730 /* Slope test */
731
732 for (int ij = 0; ij < 6; ij++) {
733 *testij += gradsti[ij][0]*gradstj[ij][0]
734 + gradsti[ij][1]*gradstj[ij][1]
735 + gradsti[ij][2]*gradstj[ij][2];
736 testi[ij] = gradsti[ij][0]*i_face_normal[0]
737 + gradsti[ij][1]*i_face_normal[1]
738 + gradsti[ij][2]*i_face_normal[2];
739 testj[ij] = gradstj[ij][0]*i_face_normal[0]
740 + gradstj[ij][1]*i_face_normal[1]
741 + gradstj[ij][2]*i_face_normal[2];
742
743 if (i_massflux > 0.) {
744 dcc[ij] = gradi[ij][0]*i_face_normal[0]
745 + gradi[ij][1]*i_face_normal[1]
746 + gradi[ij][2]*i_face_normal[2];
747 ddi[ij] = testi[ij];
748 ddj[ij] = (pj[ij]-pi[ij])/distf *srfan;
749 }
750 else {
751 dcc[ij] = gradj[ij][0]*i_face_normal[0]
752 + gradj[ij][1]*i_face_normal[1]
753 + gradj[ij][2]*i_face_normal[2];
754 ddi[ij] = (pj[ij]-pi[ij])/distf *srfan;
755 ddj[ij] = testj[ij];
756 }
757
758 *tesqck += cs_math_sq(dcc[ij]) - cs_math_sq(ddi[ij]-ddj[ij]);
759 }
760 }
761
762 /*----------------------------------------------------------------------------*/
763 /*!
764 * \brief Reconstruct values in I' and J'.
765 *
766 * \param[in] bldfrp reconstruction blending factor
767 * \param[in] diipf distance II'
768 * \param[in] djjpf distance JJ'
769 * \param[in] gradi gradient at cell i
770 * \param[in] gradj gradient at cell j
771 * \param[in] pi value at cell i
772 * \param[in] pj value at cell j
773 * \param[out] recoi reconstruction at cell i
774 * \param[out] recoj reconstruction at cell j
775 * \param[out] pip reconstructed value at cell i
776 * \param[out] pjp reconstructed value at cell j
777 */
778 /*----------------------------------------------------------------------------*/
779
780 inline static void
cs_i_compute_quantities(const cs_real_t bldfrp,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_3_t gradi,const cs_real_3_t gradj,const cs_real_t pi,const cs_real_t pj,cs_real_t * recoi,cs_real_t * recoj,cs_real_t * pip,cs_real_t * pjp)781 cs_i_compute_quantities(const cs_real_t bldfrp,
782 const cs_real_3_t diipf,
783 const cs_real_3_t djjpf,
784 const cs_real_3_t gradi,
785 const cs_real_3_t gradj,
786 const cs_real_t pi,
787 const cs_real_t pj,
788 cs_real_t *recoi,
789 cs_real_t *recoj,
790 cs_real_t *pip,
791 cs_real_t *pjp)
792 {
793 cs_real_t gradpf[3] = {0.5*(gradi[0] + gradj[0]),
794 0.5*(gradi[1] + gradj[1]),
795 0.5*(gradi[2] + gradj[2])};
796
797 /* reconstruction only if IRCFLP = 1 */
798 *recoi = bldfrp*(cs_math_3_dot_product(gradpf, diipf));
799 *recoj = bldfrp*(cs_math_3_dot_product(gradpf, djjpf));
800 *pip = pi + *recoi;
801 *pjp = pj + *recoj;
802 }
803
804 /*----------------------------------------------------------------------------*/
805 /*!
806 * \brief Reconstruct values in I' and J'.
807 *
808 * \param[in] bldfrp reconstruction blending factor
809 * \param[in] diipf distance II'
810 * \param[in] djjpf distance JJ'
811 * \param[in] gradi gradient at cell i
812 * \param[in] gradj gradient at cell j
813 * \param[in] pi value at cell i
814 * \param[in] pj value at cell j
815 * \param[out] recoi reconstruction at cell i
816 * \param[out] recoj reconstruction at cell j
817 * \param[out] pip reconstructed value at cell i
818 * \param[out] pjp reconstructed value at cell j
819 */
820 /*----------------------------------------------------------------------------*/
821
822 inline static void
cs_i_compute_quantities_vector(const cs_real_t bldfrp,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_33_t gradi,const cs_real_33_t gradj,const cs_real_3_t pi,const cs_real_3_t pj,cs_real_t recoi[3],cs_real_t recoj[3],cs_real_t pip[3],cs_real_t pjp[3])823 cs_i_compute_quantities_vector(const cs_real_t bldfrp,
824 const cs_real_3_t diipf,
825 const cs_real_3_t djjpf,
826 const cs_real_33_t gradi,
827 const cs_real_33_t gradj,
828 const cs_real_3_t pi,
829 const cs_real_3_t pj,
830 cs_real_t recoi[3],
831 cs_real_t recoj[3],
832 cs_real_t pip[3],
833 cs_real_t pjp[3])
834 {
835 cs_real_3_t dpvf;
836
837 /* x-y-z components, p = u, v, w */
838
839 for (int isou = 0; isou < 3; isou++) {
840
841 for (int jsou = 0; jsou < 3; jsou++)
842 dpvf[jsou] = 0.5*( gradi[isou][jsou]
843 + gradj[isou][jsou]);
844
845 /* reconstruction only if IRCFLP = 1 */
846
847 recoi[isou] = bldfrp*(cs_math_3_dot_product(dpvf, diipf));
848 recoj[isou] = bldfrp*(cs_math_3_dot_product(dpvf, djjpf));
849
850 pip[isou] = pi[isou] + recoi[isou];
851 pjp[isou] = pj[isou] + recoj[isou];
852
853 }
854 }
855
856 /*----------------------------------------------------------------------------*/
857 /*!
858 * \brief Reconstruct values in I' and J'.
859 *
860 * \param[in] bldfrp reconstruction blending factor
861 * \param[in] diipf distance II'
862 * \param[in] djjpf distance JJ'
863 * \param[in] gradi gradient at cell i
864 * \param[in] gradj gradient at cell j
865 * \param[in] pi value at cell i
866 * \param[in] pj value at cell j
867 * \param[out] recoi reconstruction at cell i
868 * \param[out] recoj reconstruction at cell j
869 * \param[out] pip reconstructed value at cell i
870 * \param[out] pjp reconstructed value at cell j
871 */
872 /*----------------------------------------------------------------------------*/
873
874 inline static void
cs_i_compute_quantities_tensor(const cs_real_t bldfrp,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_63_t gradi,const cs_real_63_t gradj,const cs_real_6_t pi,const cs_real_6_t pj,cs_real_t recoi[6],cs_real_t recoj[6],cs_real_t pip[6],cs_real_t pjp[6])875 cs_i_compute_quantities_tensor(const cs_real_t bldfrp,
876 const cs_real_3_t diipf,
877 const cs_real_3_t djjpf,
878 const cs_real_63_t gradi,
879 const cs_real_63_t gradj,
880 const cs_real_6_t pi,
881 const cs_real_6_t pj,
882 cs_real_t recoi[6],
883 cs_real_t recoj[6],
884 cs_real_t pip[6],
885 cs_real_t pjp[6])
886 {
887 cs_real_3_t dpvf;
888
889 /* x-y-z components, p = u, v, w */
890
891 for (int isou = 0; isou < 6; isou++) {
892
893 for (int jsou = 0; jsou < 3; jsou++)
894 dpvf[jsou] = 0.5*( gradi[isou][jsou]
895 + gradj[isou][jsou]);
896
897 /* reconstruction only if IRCFLP = 1 */
898
899 recoi[isou] = bldfrp*(cs_math_3_dot_product(dpvf, diipf));
900 recoj[isou] = bldfrp*(cs_math_3_dot_product(dpvf, djjpf));
901
902 pip[isou] = pi[isou] + recoi[isou];
903 pjp[isou] = pj[isou] + recoj[isou];
904
905 }
906 }
907
908 /*----------------------------------------------------------------------------*/
909 /*!
910 * \brief Compute relaxed values at cell i and j.
911 *
912 * \param[in] relaxp relaxation coefficient
913 * \param[in] pia old value at cell i
914 * \param[in] pja old value at cell j
915 * \param[in] recoi reconstruction at cell i
916 * \param[in] recoj reconstruction at cell j
917 * \param[in] pi value at cell i
918 * \param[in] pj value at cell j
919 * \param[out] pir relaxed value at cell i
920 * \param[out] pjr relaxed value at cell j
921 * \param[out] pipr relaxed reconstructed value at cell i
922 * \param[out] pjpr relaxed reconstructed value at cell j
923 */
924 /*----------------------------------------------------------------------------*/
925
926 inline static void
cs_i_relax_c_val(const double relaxp,const cs_real_t pia,const cs_real_t pja,const cs_real_t recoi,const cs_real_t recoj,const cs_real_t pi,const cs_real_t pj,cs_real_t * pir,cs_real_t * pjr,cs_real_t * pipr,cs_real_t * pjpr)927 cs_i_relax_c_val(const double relaxp,
928 const cs_real_t pia,
929 const cs_real_t pja,
930 const cs_real_t recoi,
931 const cs_real_t recoj,
932 const cs_real_t pi,
933 const cs_real_t pj,
934 cs_real_t *pir,
935 cs_real_t *pjr,
936 cs_real_t *pipr,
937 cs_real_t *pjpr)
938 {
939 *pir = pi/relaxp - (1.-relaxp)/relaxp * pia;
940 *pjr = pj/relaxp - (1.-relaxp)/relaxp * pja;
941
942 *pipr = *pir + recoi;
943 *pjpr = *pjr + recoj;
944 }
945
946 /*----------------------------------------------------------------------------*/
947 /*!
948 * \brief Compute relaxed values at cell i and j.
949 *
950 * \param[in] relaxp relaxation coefficient
951 * \param[in] pia old value at cell i
952 * \param[in] pja old value at cell j
953 * \param[in] recoi reconstruction at cell i
954 * \param[in] recoj reconstruction at cell j
955 * \param[in] pi value at cell i
956 * \param[in] pj value at cell j
957 * \param[out] pir relaxed value at cell i
958 * \param[out] pjr relaxed value at cell j
959 * \param[out] pipr relaxed reconstructed value at cell i
960 * \param[out] pjpr relaxed reconstructed value at cell j
961 */
962 /*----------------------------------------------------------------------------*/
963
964 inline static void
cs_i_relax_c_val_vector(const double relaxp,const cs_real_3_t pia,const cs_real_3_t pja,const cs_real_3_t recoi,const cs_real_3_t recoj,const cs_real_3_t pi,const cs_real_3_t pj,cs_real_t pir[3],cs_real_t pjr[3],cs_real_t pipr[3],cs_real_t pjpr[3])965 cs_i_relax_c_val_vector(const double relaxp,
966 const cs_real_3_t pia,
967 const cs_real_3_t pja,
968 const cs_real_3_t recoi,
969 const cs_real_3_t recoj,
970 const cs_real_3_t pi,
971 const cs_real_3_t pj,
972 cs_real_t pir[3],
973 cs_real_t pjr[3],
974 cs_real_t pipr[3],
975 cs_real_t pjpr[3])
976 {
977 for (int isou = 0; isou < 3; isou++) {
978 pir[isou] = pi[isou] /relaxp - (1.-relaxp)/relaxp * pia[isou];
979 pjr[isou] = pj[isou] /relaxp - (1.-relaxp)/relaxp * pja[isou];
980
981 pipr[isou] = pir[isou] + recoi[isou];
982 pjpr[isou] = pjr[isou] + recoj[isou];
983 }
984 }
985
986 /*----------------------------------------------------------------------------*/
987 /*!
988 * \brief Compute relaxed values at cell i and j.
989 *
990 * \param[in] relaxp relaxation coefficient
991 * \param[in] pia old value at cell i
992 * \param[in] pja old value at cell j
993 * \param[in] recoi reconstruction at cell i
994 * \param[in] recoj reconstruction at cell j
995 * \param[in] pi value at cell i
996 * \param[in] pj value at cell j
997 * \param[out] pir relaxed value at cell i
998 * \param[out] pjr relaxed value at cell j
999 * \param[out] pipr relaxed reconstructed value at cell i
1000 * \param[out] pjpr relaxed reconstructed value at cell j
1001 */
1002 /*----------------------------------------------------------------------------*/
1003
1004 inline static void
cs_i_relax_c_val_tensor(const cs_real_t relaxp,const cs_real_t pia[6],const cs_real_t pja[6],const cs_real_t recoi[6],const cs_real_t recoj[6],const cs_real_t pi[6],const cs_real_t pj[6],cs_real_t pir[6],cs_real_t pjr[6],cs_real_t pipr[6],cs_real_t pjpr[6])1005 cs_i_relax_c_val_tensor(const cs_real_t relaxp,
1006 const cs_real_t pia[6],
1007 const cs_real_t pja[6],
1008 const cs_real_t recoi[6],
1009 const cs_real_t recoj[6],
1010 const cs_real_t pi[6],
1011 const cs_real_t pj[6],
1012 cs_real_t pir[6],
1013 cs_real_t pjr[6],
1014 cs_real_t pipr[6],
1015 cs_real_t pjpr[6])
1016 {
1017 for (int isou = 0; isou < 6; isou++) {
1018 pir[isou] = pi[isou] /relaxp - (1.-relaxp)/relaxp * pia[isou];
1019 pjr[isou] = pj[isou] /relaxp - (1.-relaxp)/relaxp * pja[isou];
1020
1021 pipr[isou] = pir[isou] + recoi[isou];
1022 pjpr[isou] = pjr[isou] + recoj[isou];
1023 }
1024 }
1025
1026 /*----------------------------------------------------------------------------*/
1027 /*!
1028 * \brief Prepare value at face ij by using an upwind scheme.
1029 *
1030 * \param[in] p value at cell
1031 * \param[out] pf value at face
1032 */
1033 /*----------------------------------------------------------------------------*/
1034
1035 inline static void
cs_upwind_f_val(const cs_real_t p,cs_real_t * pf)1036 cs_upwind_f_val(const cs_real_t p,
1037 cs_real_t *pf)
1038 {
1039 *pf = p;
1040 }
1041
1042 /*----------------------------------------------------------------------------*/
1043 /*!
1044 * \brief Prepare value at face ij by using an upwind scheme.
1045 *
1046 * \param[in] p value at cell
1047 * \param[out] pf value at face
1048 */
1049 /*----------------------------------------------------------------------------*/
1050
1051 inline static void
cs_upwind_f_val_vector(const cs_real_3_t p,cs_real_t pf[3])1052 cs_upwind_f_val_vector(const cs_real_3_t p,
1053 cs_real_t pf[3])
1054 {
1055 for (int isou = 0; isou < 3; isou++)
1056 pf[isou] = p[isou];
1057 }
1058
1059 /*----------------------------------------------------------------------------*/
1060 /*!
1061 * \brief Prepare value at face ij by using an upwind scheme.
1062 *
1063 * \param[in] p value at cell
1064 * \param[out] pf value at face
1065 */
1066 /*----------------------------------------------------------------------------*/
1067
1068 inline static void
cs_upwind_f_val_tensor(const cs_real_6_t p,cs_real_t pf[6])1069 cs_upwind_f_val_tensor(const cs_real_6_t p,
1070 cs_real_t pf[6])
1071 {
1072 for (int isou = 0; isou < 6; isou++)
1073 pf[isou] = p[isou];
1074 }
1075
1076 /*----------------------------------------------------------------------------*/
1077 /*!
1078 * \brief Prepare value at face ij by using a centered scheme.
1079 *
1080 * \param[in] pnd weight
1081 * \param[in] pip (relaxed) reconstructed value at cell i
1082 * \param[in] pjp (relaxed) reconstructed value at cell j
1083 * \param[out] pf face value
1084 */
1085 /*----------------------------------------------------------------------------*/
1086
1087 inline static void
cs_centered_f_val(const double pnd,const cs_real_t pip,const cs_real_t pjp,cs_real_t * pf)1088 cs_centered_f_val(const double pnd,
1089 const cs_real_t pip,
1090 const cs_real_t pjp,
1091 cs_real_t *pf)
1092 {
1093 *pf = pnd*pip + (1.-pnd)*pjp;
1094 }
1095
1096 /*----------------------------------------------------------------------------*/
1097 /*!
1098 * \brief Prepare value at face ij by using a centered scheme.
1099 *
1100 * \param[in] pnd weight
1101 * \param[in] pip (relaxed) reconstructed value at cell i
1102 * \param[in] pjp (relaxed) reconstructed value at cell j
1103 * \param[out] pf face value
1104 */
1105 /*----------------------------------------------------------------------------*/
1106
1107 inline static void
cs_centered_f_val_vector(const double pnd,const cs_real_3_t pip,const cs_real_3_t pjp,cs_real_t pf[3])1108 cs_centered_f_val_vector(const double pnd,
1109 const cs_real_3_t pip,
1110 const cs_real_3_t pjp,
1111 cs_real_t pf[3])
1112 {
1113 for (int isou = 0; isou < 3; isou++)
1114 pf[isou] = pnd*pip[isou] + (1.-pnd)*pjp[isou];
1115 }
1116
1117 /*----------------------------------------------------------------------------*/
1118 /*!
1119 * \brief Prepare value at face ij by using a centered scheme.
1120 *
1121 * \param[in] pnd weight
1122 * \param[in] pip (relaxed) reconstructed value at cell i
1123 * \param[in] pjp (relaxed) reconstructed value at cell j
1124 * \param[out] pf face value
1125 */
1126 /*----------------------------------------------------------------------------*/
1127
1128 inline static void
cs_centered_f_val_tensor(const double pnd,const cs_real_6_t pip,const cs_real_6_t pjp,cs_real_t pf[6])1129 cs_centered_f_val_tensor(const double pnd,
1130 const cs_real_6_t pip,
1131 const cs_real_6_t pjp,
1132 cs_real_t pf[6])
1133 {
1134 for (int isou = 0; isou < 6; isou++)
1135 pf[isou] = pnd*pip[isou] + (1.-pnd)*pjp[isou];
1136 }
1137
1138 /*----------------------------------------------------------------------------*/
1139 /*!
1140 * \brief Prepare value at face ij by using a Second Order Linear Upwind scheme.
1141 *
1142 * \param[in] cell_cen center of gravity coordinates of cell
1143 * \param[in] i_face_cog center of gravity coordinates of face
1144 * \param[in] grad gradient at cell
1145 * \param[in] p (relaced) value at cell
1146 * \param[out] pf face value
1147 */
1148 /*----------------------------------------------------------------------------*/
1149
1150 inline static void
cs_solu_f_val(const cs_real_3_t cell_cen,const cs_real_3_t i_face_cog,const cs_real_3_t grad,const cs_real_t p,cs_real_t * pf)1151 cs_solu_f_val(const cs_real_3_t cell_cen,
1152 const cs_real_3_t i_face_cog,
1153 const cs_real_3_t grad,
1154 const cs_real_t p,
1155 cs_real_t *pf)
1156 {
1157 cs_real_t df[3];
1158
1159 df[0] = i_face_cog[0] - cell_cen[0];
1160 df[1] = i_face_cog[1] - cell_cen[1];
1161 df[2] = i_face_cog[2] - cell_cen[2];
1162
1163 *pf = p + cs_math_3_dot_product(df, grad);
1164 }
1165
1166 /*----------------------------------------------------------------------------*/
1167 /*!
1168 * \brief Prepare value at face ij by using a Second Order Linear Upwind scheme.
1169 *
1170 * \param[in] cell_cen center of gravity coordinates of cell
1171 * \param[in] i_face_cog center of gravity coordinates of face
1172 * \param[in] grad gradient at cell
1173 * \param[in] p (relaced) value at cell
1174 * \param[out] pf face value
1175 */
1176 /*----------------------------------------------------------------------------*/
1177
1178 inline static void
cs_solu_f_val_vector(const cs_real_3_t cell_cen,const cs_real_3_t i_face_cog,const cs_real_33_t grad,const cs_real_3_t p,cs_real_t pf[3])1179 cs_solu_f_val_vector(const cs_real_3_t cell_cen,
1180 const cs_real_3_t i_face_cog,
1181 const cs_real_33_t grad,
1182 const cs_real_3_t p,
1183 cs_real_t pf[3])
1184 {
1185 cs_real_t df[3];
1186
1187 for (int jsou = 0; jsou < 3; jsou++)
1188 df[jsou] = i_face_cog[jsou] - cell_cen[jsou];
1189
1190 for (int isou = 0; isou < 3; isou++) {
1191 pf[isou] = p[isou] + df[0]*grad[isou][0]
1192 + df[1]*grad[isou][1]
1193 + df[2]*grad[isou][2];
1194
1195 }
1196 }
1197
1198 /*----------------------------------------------------------------------------*/
1199 /*!
1200 * \brief Prepare value at face ij by using a Second Order Linear Upwind scheme.
1201 *
1202 * \param[in] cell_cen center of gravity coordinates of cell
1203 * \param[in] i_face_cog center of gravity coordinates of face
1204 * \param[in] grad gradient at cell
1205 * \param[in] p (relaced) value at cell
1206 * \param[out] pf face value
1207 */
1208 /*----------------------------------------------------------------------------*/
1209
1210 inline static void
cs_solu_f_val_tensor(const cs_real_3_t cell_cen,const cs_real_3_t i_face_cog,const cs_real_63_t grad,const cs_real_6_t p,cs_real_t pf[6])1211 cs_solu_f_val_tensor(const cs_real_3_t cell_cen,
1212 const cs_real_3_t i_face_cog,
1213 const cs_real_63_t grad,
1214 const cs_real_6_t p,
1215 cs_real_t pf[6])
1216 {
1217 cs_real_t df[3];
1218
1219 for (int jsou = 0; jsou < 3; jsou++)
1220 df[jsou] = i_face_cog[jsou] - cell_cen[jsou];
1221
1222 for (int isou = 0; isou < 6; isou++) {
1223 pf[isou] = p[isou] + df[0]*grad[isou][0]
1224 + df[1]*grad[isou][1]
1225 + df[2]*grad[isou][2];
1226 }
1227 }
1228
1229 /*----------------------------------------------------------------------------*/
1230 /*!
1231 * \brief Blend face values for a centered or SOLU scheme with face values for
1232 * an upwind scheme.
1233 *
1234 * \param[in] blencp proportion of second order scheme,
1235 * (1-blencp) is the proportion of upwind.
1236 * \param[in] p (relaxed) value at cell
1237 * \param[out] pf face value
1238 */
1239 /*----------------------------------------------------------------------------*/
1240
1241 inline static void
cs_blend_f_val(const double blencp,const cs_real_t p,cs_real_t * pf)1242 cs_blend_f_val(const double blencp,
1243 const cs_real_t p,
1244 cs_real_t *pf)
1245 {
1246 *pf = blencp * (*pf) + (1. - blencp) * p;
1247 }
1248
1249 /*----------------------------------------------------------------------------*/
1250 /*!
1251 * \brief Blend face values for a centered or SOLU scheme with face values for
1252 * an upwind scheme.
1253 *
1254 * \param[in] blencp proportion of second order scheme,
1255 * (1-blencp) is the proportion of upwind.
1256 * \param[in] p (relaxed) value at cell
1257 * \param[in,out] pf face value
1258 */
1259 /*----------------------------------------------------------------------------*/
1260
1261 inline static void
cs_blend_f_val_vector(const double blencp,const cs_real_3_t p,cs_real_t pf[3])1262 cs_blend_f_val_vector(const double blencp,
1263 const cs_real_3_t p,
1264 cs_real_t pf[3])
1265 {
1266 for (int isou = 0; isou < 3; isou++)
1267 pf[isou] = blencp*(pf[isou])+(1.-blencp)*p[isou];
1268 }
1269
1270 /*----------------------------------------------------------------------------*/
1271 /*!
1272 * \brief Blend face values for a centered or SOLU scheme with face values for
1273 * an upwind scheme.
1274 *
1275 * \param[in] blencp proportion of second order scheme,
1276 * (1-blencp) is the proportion of upwind.
1277 * \param[in] p (relaxed) value at cell
1278 * \param[out] pf face value
1279 */
1280 /*----------------------------------------------------------------------------*/
1281
1282 inline static void
cs_blend_f_val_tensor(const double blencp,const cs_real_t p[6],cs_real_t pf[6])1283 cs_blend_f_val_tensor(const double blencp,
1284 const cs_real_t p[6],
1285 cs_real_t pf[6])
1286 {
1287 for (int isou = 0; isou < 6; isou++)
1288 pf[isou] = blencp*(pf[isou])+(1.-blencp)*p[isou];
1289 }
1290
1291 /*----------------------------------------------------------------------------*/
1292 /*!
1293 * \brief Add convective fluxes (substracting the mass accumulation from them)
1294 * to fluxes at face ij.
1295 *
1296 * \param[in] iconvp convection flag
1297 * \param[in] thetap weighting coefficient for the theta-schema,
1298 * \param[in] imasac take mass accumulation into account?
1299 * \param[in] pi value at cell i
1300 * \param[in] pj value at cell j
1301 * \param[in] pifri contribution of i to flux from i to j
1302 * \param[in] pifrj contribution of i to flux from j to i
1303 * \param[in] pjfri contribution of j to flux from i to j
1304 * \param[in] pjfrj contribution of j to flux from j to i
1305 * \param[in] i_massflux mass flux at face ij
1306 * \param[in] xcppi specific heat value if the scalar is the temperature,
1307 * 1 otherwise at cell i
1308 * \param[in] xcppj specific heat value if the scalar is the temperature,
1309 * 1 otherwise at cell j
1310 * \param[in,out] fluxij fluxes at face ij
1311 */
1312 /*----------------------------------------------------------------------------*/
1313
1314 inline static void
cs_i_conv_flux(const int iconvp,const cs_real_t thetap,const int imasac,const cs_real_t pi,const cs_real_t pj,const cs_real_t pifri,const cs_real_t pifrj,const cs_real_t pjfri,const cs_real_t pjfrj,const cs_real_t i_massflux,const cs_real_t xcppi,const cs_real_t xcppj,cs_real_2_t fluxij)1315 cs_i_conv_flux(const int iconvp,
1316 const cs_real_t thetap,
1317 const int imasac,
1318 const cs_real_t pi,
1319 const cs_real_t pj,
1320 const cs_real_t pifri,
1321 const cs_real_t pifrj,
1322 const cs_real_t pjfri,
1323 const cs_real_t pjfrj,
1324 const cs_real_t i_massflux,
1325 const cs_real_t xcppi,
1326 const cs_real_t xcppj,
1327 cs_real_2_t fluxij)
1328 {
1329 cs_real_t flui, fluj;
1330
1331 flui = 0.5*(i_massflux + fabs(i_massflux));
1332 fluj = 0.5*(i_massflux - fabs(i_massflux));
1333
1334 fluxij[0] += iconvp*xcppi*(thetap*(flui*pifri + fluj*pjfri) - imasac*i_massflux*pi);
1335 fluxij[1] += iconvp*xcppj*(thetap*(flui*pifrj + fluj*pjfrj) - imasac*i_massflux*pj);
1336 }
1337
1338 /*----------------------------------------------------------------------------*/
1339 /*!
1340 * \brief Add convective fluxes (substracting the mass accumulation from them)
1341 * to fluxes at face ij.
1342 *
1343 * \param[in] iconvp convection flag
1344 * \param[in] thetap weighting coefficient for the theta-schema,
1345 * \param[in] imasac take mass accumulation into account?
1346 * \param[in] pi value at cell i
1347 * \param[in] pj value at cell j
1348 * \param[out] pifri contribution of i to flux from i to j
1349 * \param[out] pifrj contribution of i to flux from j to i
1350 * \param[out] pjfri contribution of j to flux from i to j
1351 * \param[out] pjfrj contribution of j to flux from j to i
1352 * \param[in] i_massflux mass flux at face ij
1353 * \param[in,out] fluxi fluxes at face i
1354 * \param[in,out] fluxj fluxes at face j
1355 */
1356 /*----------------------------------------------------------------------------*/
1357
1358 inline static void
cs_i_conv_flux_vector(const int iconvp,const cs_real_t thetap,const int imasac,const cs_real_t pi[3],const cs_real_t pj[3],const cs_real_t pifri[3],const cs_real_t pifrj[3],const cs_real_t pjfri[3],const cs_real_t pjfrj[3],const cs_real_t i_massflux,cs_real_t fluxi[3],cs_real_t fluxj[3])1359 cs_i_conv_flux_vector(const int iconvp,
1360 const cs_real_t thetap,
1361 const int imasac,
1362 const cs_real_t pi[3],
1363 const cs_real_t pj[3],
1364 const cs_real_t pifri[3],
1365 const cs_real_t pifrj[3],
1366 const cs_real_t pjfri[3],
1367 const cs_real_t pjfrj[3],
1368 const cs_real_t i_massflux,
1369 cs_real_t fluxi[3],
1370 cs_real_t fluxj[3])
1371 {
1372 cs_real_t flui, fluj;
1373
1374 flui = 0.5*(i_massflux + fabs(i_massflux));
1375 fluj = 0.5*(i_massflux - fabs(i_massflux));
1376
1377 for (int isou = 0; isou < 3; isou++) {
1378
1379 fluxi[isou] += iconvp*( thetap*(flui*pifri[isou] + fluj*pjfri[isou])
1380 - imasac*i_massflux*pi[isou]);
1381 fluxj[isou] += iconvp*( thetap*(flui*pifrj[isou] + fluj*pjfrj[isou])
1382 - imasac*i_massflux*pj[isou]);
1383 }
1384 }
1385
1386 /*----------------------------------------------------------------------------*/
1387 /*!
1388 * \brief Add convective fluxes (substracting the mass accumulation from them)
1389 * to fluxes at face ij.
1390 *
1391 * \param[in] iconvp convection flag
1392 * \param[in] thetap weighting coefficient for the theta-schema,
1393 * \param[in] imasac take mass accumulation into account?
1394 * \param[in] pi value at cell i
1395 * \param[in] pj value at cell j
1396 * \param[out] pifri contribution of i to flux from i to j
1397 * \param[out] pifrj contribution of i to flux from j to i
1398 * \param[out] pjfri contribution of j to flux from i to j
1399 * \param[out] pjfrj contribution of j to flux from j to i
1400 * \param[in] i_massflux mass flux at face ij
1401 * \param[in,out] fluxi fluxes at face i
1402 * \param[in,out] fluxj fluxes at face j
1403 */
1404 /*----------------------------------------------------------------------------*/
1405
1406 inline static void
cs_i_conv_flux_tensor(const int iconvp,const cs_real_t thetap,const int imasac,const cs_real_t pi[6],const cs_real_t pj[6],const cs_real_t pifri[6],const cs_real_t pifrj[6],const cs_real_t pjfri[6],const cs_real_t pjfrj[6],const cs_real_t i_massflux,cs_real_t fluxi[6],cs_real_t fluxj[6])1407 cs_i_conv_flux_tensor(const int iconvp,
1408 const cs_real_t thetap,
1409 const int imasac,
1410 const cs_real_t pi[6],
1411 const cs_real_t pj[6],
1412 const cs_real_t pifri[6],
1413 const cs_real_t pifrj[6],
1414 const cs_real_t pjfri[6],
1415 const cs_real_t pjfrj[6],
1416 const cs_real_t i_massflux,
1417 cs_real_t fluxi[6],
1418 cs_real_t fluxj[6])
1419 {
1420 cs_real_t flui, fluj;
1421
1422 flui = 0.5*(i_massflux + fabs(i_massflux));
1423 fluj = 0.5*(i_massflux - fabs(i_massflux));
1424
1425 for (int isou = 0; isou < 6; isou++) {
1426 fluxi[isou] += iconvp*( thetap*(flui*pifri[isou] + fluj*pjfri[isou])
1427 - imasac*i_massflux*pi[isou]);
1428 fluxj[isou] += iconvp*( thetap*(flui*pifrj[isou] + fluj*pjfrj[isou])
1429 - imasac*i_massflux*pj[isou]);
1430 }
1431 }
1432
1433 /*----------------------------------------------------------------------------*/
1434 /*!
1435 * \brief Add diffusive fluxes to fluxes at face ij.
1436 *
1437 * \param[in] idiffp diffusion flag
1438 * \param[in] thetap weighting coefficient for the theta-schema,
1439 * \param[in] pip reconstructed value at cell i
1440 * \param[in] pjp reconstructed value at cell j
1441 * \param[in] pipr relaxed reconstructed value at cell i
1442 * \param[in] pjpr relaxed reconstructed value at cell j
1443 * \param[in] i_visc diffusion coefficient (divided by IJ) at face ij
1444 * \param[in,out] fluxij fluxes at face ij
1445 */
1446 /*----------------------------------------------------------------------------*/
1447
1448 inline static void
cs_i_diff_flux(const int idiffp,const cs_real_t thetap,const cs_real_t pip,const cs_real_t pjp,const cs_real_t pipr,const cs_real_t pjpr,const cs_real_t i_visc,cs_real_t fluxij[2])1449 cs_i_diff_flux(const int idiffp,
1450 const cs_real_t thetap,
1451 const cs_real_t pip,
1452 const cs_real_t pjp,
1453 const cs_real_t pipr,
1454 const cs_real_t pjpr,
1455 const cs_real_t i_visc,
1456 cs_real_t fluxij[2])
1457 {
1458 fluxij[0] += idiffp*thetap*i_visc*(pipr -pjp);
1459 fluxij[1] += idiffp*thetap*i_visc*(pip -pjpr);
1460 }
1461
1462 /*----------------------------------------------------------------------------*/
1463 /*!
1464 * \brief Add diffusive fluxes to fluxes at face ij.
1465 *
1466 * \param[in] idiffp diffusion flag
1467 * \param[in] thetap weighting coefficient for the theta-schema,
1468 * \param[in] pip reconstructed value at cell i
1469 * \param[in] pjp reconstructed value at cell j
1470 * \param[in] pipr relaxed reconstructed value at cell i
1471 * \param[in] pjpr relaxed reconstructed value at cell j
1472 * \param[in] i_visc diffusion coefficient (divided by IJ) at face ij
1473 * \param[in,out] fluxi fluxes at face i
1474 * \param[in,out] fluxj fluxes at face j
1475 */
1476 /*----------------------------------------------------------------------------*/
1477
1478 inline static void
cs_i_diff_flux_vector(const int idiffp,const cs_real_t thetap,const cs_real_t pip[3],const cs_real_t pjp[3],const cs_real_t pipr[3],const cs_real_t pjpr[3],const cs_real_t i_visc,cs_real_t fluxi[3],cs_real_t fluxj[3])1479 cs_i_diff_flux_vector(const int idiffp,
1480 const cs_real_t thetap,
1481 const cs_real_t pip[3],
1482 const cs_real_t pjp[3],
1483 const cs_real_t pipr[3],
1484 const cs_real_t pjpr[3],
1485 const cs_real_t i_visc,
1486 cs_real_t fluxi[3],
1487 cs_real_t fluxj[3])
1488 {
1489 for (int isou = 0; isou < 3; isou++) {
1490 fluxi[isou] += idiffp*thetap*i_visc*(pipr[isou] -pjp[isou]);
1491 fluxj[isou] += idiffp*thetap*i_visc*(pip[isou] -pjpr[isou]);
1492 }
1493 }
1494
1495 /*----------------------------------------------------------------------------*/
1496 /*!
1497 * \brief Add diffusive fluxes to fluxes at face ij.
1498 *
1499 * \param[in] idiffp diffusion flag
1500 * \param[in] thetap weighting coefficient for the theta-schema,
1501 * \param[in] pip reconstructed value at cell i
1502 * \param[in] pjp reconstructed value at cell j
1503 * \param[in] pipr relaxed reconstructed value at cell i
1504 * \param[in] pjpr relaxed reconstructed value at cell j
1505 * \param[in] i_visc diffusion coefficient (divided by IJ) at face ij
1506 * \param[in,out] fluxi fluxes at face i
1507 * \param[in,out] fluxj fluxes at face j
1508 */
1509 /*----------------------------------------------------------------------------*/
1510
1511 inline static void
cs_i_diff_flux_tensor(const int idiffp,const cs_real_t thetap,const cs_real_t pip[6],const cs_real_t pjp[6],const cs_real_t pipr[6],const cs_real_t pjpr[6],const cs_real_t i_visc,cs_real_t fluxi[6],cs_real_t fluxj[6])1512 cs_i_diff_flux_tensor(const int idiffp,
1513 const cs_real_t thetap,
1514 const cs_real_t pip[6],
1515 const cs_real_t pjp[6],
1516 const cs_real_t pipr[6],
1517 const cs_real_t pjpr[6],
1518 const cs_real_t i_visc,
1519 cs_real_t fluxi[6],
1520 cs_real_t fluxj[6])
1521 {
1522 for (int isou = 0; isou < 6; isou++) {
1523 fluxi[isou] += idiffp*thetap*i_visc*(pipr[isou] -pjp[isou]);
1524 fluxj[isou] += idiffp*thetap*i_visc*(pip[isou] -pjpr[isou]);
1525 }
1526 }
1527
1528 /*----------------------------------------------------------------------------*/
1529 /*!
1530 * \brief Handle preparation of internal face values for the fluxes computation
1531 * in case of a steady algorithm and a pure upwind flux.
1532 *
1533 * \param[in] bldfrp reconstruction blending factor
1534 * \param[in] relaxp relaxation coefficient
1535 * \param[in] diipf distance II'
1536 * \param[in] djjpf distance JJ'
1537 * \param[in] gradi gradient at cell i
1538 * \param[in] gradj gradient at cell j
1539 * \param[in] pi value at cell i
1540 * \param[in] pj value at cell j
1541 * \param[in] pia old value at cell i
1542 * \param[in] pja old value at cell j
1543 * \param[out] pifri contribution of i to flux from i to j
1544 * \param[out] pifrj contribution of i to flux from j to i
1545 * \param[out] pjfri contribution of j to flux from i to j
1546 * \param[out] pjfrj contribution of j to flux from j to i
1547 * \param[out] pip reconstructed value at cell i
1548 * \param[out] pjp reconstructed value at cell j
1549 * \param[out] pipr relaxed reconstructed value at cell i
1550 * \param[out] pjpr relaxed reconstructed value at cell j
1551 */
1552 /*----------------------------------------------------------------------------*/
1553
1554 inline static void
cs_i_cd_steady_upwind(const cs_real_t bldfrp,const cs_real_t relaxp,const cs_real_t diipf[3],const cs_real_t djjpf[3],const cs_real_t gradi[3],const cs_real_t gradj[3],const cs_real_t pi,const cs_real_t pj,const cs_real_t pia,const cs_real_t pja,cs_real_t * pifri,cs_real_t * pifrj,cs_real_t * pjfri,cs_real_t * pjfrj,cs_real_t * pip,cs_real_t * pjp,cs_real_t * pipr,cs_real_t * pjpr)1555 cs_i_cd_steady_upwind(const cs_real_t bldfrp,
1556 const cs_real_t relaxp,
1557 const cs_real_t diipf[3],
1558 const cs_real_t djjpf[3],
1559 const cs_real_t gradi[3],
1560 const cs_real_t gradj[3],
1561 const cs_real_t pi,
1562 const cs_real_t pj,
1563 const cs_real_t pia,
1564 const cs_real_t pja,
1565 cs_real_t *pifri,
1566 cs_real_t *pifrj,
1567 cs_real_t *pjfri,
1568 cs_real_t *pjfrj,
1569 cs_real_t *pip,
1570 cs_real_t *pjp,
1571 cs_real_t *pipr,
1572 cs_real_t *pjpr)
1573 {
1574 cs_real_t pir, pjr;
1575 cs_real_t recoi, recoj;
1576
1577 cs_i_compute_quantities(bldfrp,
1578 diipf,
1579 djjpf,
1580 gradi,
1581 gradj,
1582 pi,
1583 pj,
1584 &recoi,
1585 &recoj,
1586 pip,
1587 pjp);
1588
1589 cs_i_relax_c_val(relaxp,
1590 pia,
1591 pja,
1592 recoi,
1593 recoj,
1594 pi,
1595 pj,
1596 &pir,
1597 &pjr,
1598 pipr,
1599 pjpr);
1600
1601 cs_upwind_f_val(pi,
1602 pifrj);
1603 cs_upwind_f_val(pir,
1604 pifri);
1605 cs_upwind_f_val(pj,
1606 pjfri);
1607 cs_upwind_f_val(pjr,
1608 pjfrj);
1609 }
1610
1611 /*----------------------------------------------------------------------------*/
1612 /*!
1613 * \brief Handle preparation of internal face values for the fluxes computation
1614 * in case of a steady algorithm and a pure upwind flux.
1615 *
1616 * \param[in] bldfrp reconstruction blending factor
1617 * \param[in] relaxp relaxation coefficient
1618 * \param[in] diipf distance II'
1619 * \param[in] djjpf distance JJ'
1620 * \param[in] gradi gradient at cell i
1621 * \param[in] gradj gradient at cell j
1622 * \param[in] pi value at cell i
1623 * \param[in] pj value at cell j
1624 * \param[in] pia old value at cell i
1625 * \param[in] pja old value at cell j
1626 * \param[out] pifri contribution of i to flux from i to j
1627 * \param[out] pifrj contribution of i to flux from j to i
1628 * \param[out] pjfri contribution of j to flux from i to j
1629 * \param[out] pjfrj contribution of j to flux from j to i
1630 * \param[out] pip reconstructed value at cell i
1631 * \param[out] pjp reconstructed value at cell j
1632 * \param[out] pipr relaxed reconstructed value at cell i
1633 * \param[out] pjpr relaxed reconstructed value at cell j
1634 */
1635 /*----------------------------------------------------------------------------*/
1636
1637 inline static void
cs_i_cd_steady_upwind_vector(const cs_real_t bldfrp,const cs_real_t relaxp,const cs_real_t diipf[3],const cs_real_t djjpf[3],const cs_real_t gradi[3][3],const cs_real_t gradj[3][3],const cs_real_t pi[3],const cs_real_t pj[3],const cs_real_t pia[3],const cs_real_t pja[3],cs_real_t pifri[3],cs_real_t pifrj[3],cs_real_t pjfri[3],cs_real_t pjfrj[3],cs_real_t pip[3],cs_real_t pjp[3],cs_real_t pipr[3],cs_real_t pjpr[3])1638 cs_i_cd_steady_upwind_vector(const cs_real_t bldfrp,
1639 const cs_real_t relaxp,
1640 const cs_real_t diipf[3],
1641 const cs_real_t djjpf[3],
1642 const cs_real_t gradi[3][3],
1643 const cs_real_t gradj[3][3],
1644 const cs_real_t pi[3],
1645 const cs_real_t pj[3],
1646 const cs_real_t pia[3],
1647 const cs_real_t pja[3],
1648 cs_real_t pifri[3],
1649 cs_real_t pifrj[3],
1650 cs_real_t pjfri[3],
1651 cs_real_t pjfrj[3],
1652 cs_real_t pip[3],
1653 cs_real_t pjp[3],
1654 cs_real_t pipr[3],
1655 cs_real_t pjpr[3])
1656 {
1657 cs_real_t pir[3], pjr[3];
1658 cs_real_t recoi[3], recoj[3];
1659
1660 cs_i_compute_quantities_vector(bldfrp,
1661 diipf,
1662 djjpf,
1663 gradi,
1664 gradj,
1665 pi,
1666 pj,
1667 recoi,
1668 recoj,
1669 pip,
1670 pjp);
1671
1672 cs_i_relax_c_val_vector(relaxp,
1673 pia,
1674 pja,
1675 recoi,
1676 recoj,
1677 pi,
1678 pj,
1679 pir,
1680 pjr,
1681 pipr,
1682 pjpr);
1683
1684 cs_upwind_f_val_vector(pi,
1685 pifrj);
1686 cs_upwind_f_val_vector(pir,
1687 pifri);
1688 cs_upwind_f_val_vector(pj,
1689 pjfri);
1690 cs_upwind_f_val_vector(pjr,
1691 pjfrj);
1692 }
1693
1694 /*----------------------------------------------------------------------------*/
1695 /*!
1696 * \brief Handle preparation of internal face values for the fluxes computation
1697 * in case of a steady algorithm and a pure upwind flux.
1698 *
1699 * \param[in] bldfrp reconstruction blending factor
1700 * \param[in] relaxp relaxation coefficient
1701 * \param[in] diipf distance II'
1702 * \param[in] djjpf distance JJ'
1703 * \param[in] gradi gradient at cell i
1704 * \param[in] gradj gradient at cell j
1705 * \param[in] pi value at cell i
1706 * \param[in] pj value at cell j
1707 * \param[in] pia old value at cell i
1708 * \param[in] pja old value at cell j
1709 * \param[out] pifri contribution of i to flux from i to j
1710 * \param[out] pifrj contribution of i to flux from j to i
1711 * \param[out] pjfri contribution of j to flux from i to j
1712 * \param[out] pjfrj contribution of j to flux from j to i
1713 * \param[out] pip reconstructed value at cell i
1714 * \param[out] pjp reconstructed value at cell j
1715 * \param[out] pipr relaxed reconstructed value at cell i
1716 * \param[out] pjpr relaxed reconstructed value at cell j
1717 */
1718 /*----------------------------------------------------------------------------*/
1719
1720 inline static void
cs_i_cd_steady_upwind_tensor(const cs_real_t bldfrp,const cs_real_t relaxp,const cs_real_t diipf[3],const cs_real_t djjpf[3],const cs_real_t gradi[6][3],const cs_real_t gradj[6][3],const cs_real_t pi[6],const cs_real_t pj[6],const cs_real_t pia[6],const cs_real_t pja[6],cs_real_t pifri[6],cs_real_t pifrj[6],cs_real_t pjfri[6],cs_real_t pjfrj[6],cs_real_t pip[6],cs_real_t pjp[6],cs_real_t pipr[6],cs_real_t pjpr[6])1721 cs_i_cd_steady_upwind_tensor(const cs_real_t bldfrp,
1722 const cs_real_t relaxp,
1723 const cs_real_t diipf[3],
1724 const cs_real_t djjpf[3],
1725 const cs_real_t gradi[6][3],
1726 const cs_real_t gradj[6][3],
1727 const cs_real_t pi[6],
1728 const cs_real_t pj[6],
1729 const cs_real_t pia[6],
1730 const cs_real_t pja[6],
1731 cs_real_t pifri[6],
1732 cs_real_t pifrj[6],
1733 cs_real_t pjfri[6],
1734 cs_real_t pjfrj[6],
1735 cs_real_t pip[6],
1736 cs_real_t pjp[6],
1737 cs_real_t pipr[6],
1738 cs_real_t pjpr[6])
1739 {
1740 cs_real_t pir[6], pjr[6];
1741 cs_real_t recoi[6], recoj[6];
1742
1743 cs_i_compute_quantities_tensor(bldfrp,
1744 diipf,
1745 djjpf,
1746 gradi,
1747 gradj,
1748 pi,
1749 pj,
1750 recoi,
1751 recoj,
1752 pip,
1753 pjp);
1754
1755 cs_i_relax_c_val_tensor(relaxp,
1756 pia,
1757 pja,
1758 recoi,
1759 recoj,
1760 pi,
1761 pj,
1762 pir,
1763 pjr,
1764 pipr,
1765 pjpr);
1766
1767 cs_upwind_f_val_tensor(pi,
1768 pifrj);
1769 cs_upwind_f_val_tensor(pir,
1770 pifri);
1771 cs_upwind_f_val_tensor(pj,
1772 pjfri);
1773 cs_upwind_f_val_tensor(pjr,
1774 pjfrj);
1775 }
1776
1777 /*----------------------------------------------------------------------------*/
1778 /*!
1779 * \brief Handle preparation of internal face values for the fluxes computation
1780 * in case of an unsteady algorithm and a pure upwind flux.
1781 *
1782 * \param[in] bldfrp reconstruction blending factor
1783 * \param[in] diipf distance II'
1784 * \param[in] djjpf distance JJ'
1785 * \param[in] gradi gradient at cell i
1786 * \param[in] gradj gradient at cell j
1787 * \param[in] pi value at cell i
1788 * \param[in] pj value at cell j
1789 * \param[out] pif contribution of i to flux from i to j
1790 * \param[out] pjf contribution of j to flux from i to j
1791 * \param[out] pip reconstructed value at cell i
1792 * \param[out] pjp reconstructed value at cell j
1793 */
1794 /*----------------------------------------------------------------------------*/
1795
1796 inline static void
cs_i_cd_unsteady_upwind(const cs_real_t bldfrp,const cs_real_t diipf[3],const cs_real_t djjpf[3],const cs_real_t gradi[3],const cs_real_t gradj[3],const cs_real_t pi,const cs_real_t pj,cs_real_t * pif,cs_real_t * pjf,cs_real_t * pip,cs_real_t * pjp)1797 cs_i_cd_unsteady_upwind(const cs_real_t bldfrp,
1798 const cs_real_t diipf[3],
1799 const cs_real_t djjpf[3],
1800 const cs_real_t gradi[3],
1801 const cs_real_t gradj[3],
1802 const cs_real_t pi,
1803 const cs_real_t pj,
1804 cs_real_t *pif,
1805 cs_real_t *pjf,
1806 cs_real_t *pip,
1807 cs_real_t *pjp)
1808 {
1809 cs_real_t recoi, recoj;
1810
1811 cs_i_compute_quantities(bldfrp,
1812 diipf,
1813 djjpf,
1814 gradi,
1815 gradj,
1816 pi,
1817 pj,
1818 &recoi,
1819 &recoj,
1820 pip,
1821 pjp);
1822
1823 cs_upwind_f_val(pi, pif);
1824 cs_upwind_f_val(pj, pjf);
1825 }
1826
1827 /*----------------------------------------------------------------------------*/
1828 /*!
1829 * \brief Handle preparation of internal face values for the fluxes computation
1830 * in case of an unsteady algorithm and a pure upwind flux.
1831 *
1832 * \param[in] bldfrp reconstruction blending factor
1833 * \param[in] diipf distance II'
1834 * \param[in] djjpf distance JJ'
1835 * \param[in] gradi gradient at cell i
1836 * \param[in] gradj gradient at cell j
1837 * \param[in] pi value at cell i
1838 * \param[in] pj value at cell j
1839 * \param[out] pif contribution of i to flux from i to j
1840 * \param[out] pjf contribution of j to flux from i to j
1841 * \param[out] pip reconstructed value at cell i
1842 * \param[out] pjp reconstructed value at cell j
1843 */
1844 /*----------------------------------------------------------------------------*/
1845
1846 inline static void
cs_i_cd_unsteady_upwind_vector(const cs_real_t bldfrp,const cs_real_t diipf[3],const cs_real_t djjpf[3],const cs_real_t gradi[3][3],const cs_real_t gradj[3][3],const cs_real_t pi[3],const cs_real_t pj[3],cs_real_t pif[3],cs_real_t pjf[3],cs_real_t pip[3],cs_real_t pjp[3])1847 cs_i_cd_unsteady_upwind_vector(const cs_real_t bldfrp,
1848 const cs_real_t diipf[3],
1849 const cs_real_t djjpf[3],
1850 const cs_real_t gradi[3][3],
1851 const cs_real_t gradj[3][3],
1852 const cs_real_t pi[3],
1853 const cs_real_t pj[3],
1854 cs_real_t pif[3],
1855 cs_real_t pjf[3],
1856 cs_real_t pip[3],
1857 cs_real_t pjp[3])
1858 {
1859 cs_real_t recoi[3], recoj[3];
1860
1861 cs_i_compute_quantities_vector(bldfrp,
1862 diipf,
1863 djjpf,
1864 gradi,
1865 gradj,
1866 pi,
1867 pj,
1868 recoi,
1869 recoj,
1870 pip,
1871 pjp);
1872
1873 cs_upwind_f_val_vector(pi, pif);
1874 cs_upwind_f_val_vector(pj, pjf);
1875 }
1876
1877 /*----------------------------------------------------------------------------*/
1878 /*!
1879 * \brief Handle preparation of internal face values for the fluxes computation
1880 * in case of an unsteady algorithm and a pure upwind flux.
1881 *
1882 * \param[in] bldfrp reconstruction blending factor
1883 * \param[in] diipf distance II'
1884 * \param[in] djjpf distance JJ'
1885 * \param[in] gradi gradient at cell i
1886 * \param[in] gradj gradient at cell j
1887 * \param[in] pi value at cell i
1888 * \param[in] pj value at cell j
1889 * \param[out] pif contribution of i to flux from i to j
1890 * \param[out] pjf contribution of j to flux from i to j
1891 * \param[out] pip reconstructed value at cell i
1892 * \param[out] pjp reconstructed value at cell j
1893 */
1894 /*----------------------------------------------------------------------------*/
1895
1896 inline static void
cs_i_cd_unsteady_upwind_tensor(const cs_real_t bldfrp,const cs_real_t diipf[3],const cs_real_t djjpf[3],const cs_real_t gradi[6][3],const cs_real_t gradj[6][3],const cs_real_t pi[6],const cs_real_t pj[6],cs_real_t pif[6],cs_real_t pjf[6],cs_real_t pip[6],cs_real_t pjp[6])1897 cs_i_cd_unsteady_upwind_tensor(const cs_real_t bldfrp,
1898 const cs_real_t diipf[3],
1899 const cs_real_t djjpf[3],
1900 const cs_real_t gradi[6][3],
1901 const cs_real_t gradj[6][3],
1902 const cs_real_t pi[6],
1903 const cs_real_t pj[6],
1904 cs_real_t pif[6],
1905 cs_real_t pjf[6],
1906 cs_real_t pip[6],
1907 cs_real_t pjp[6])
1908 {
1909 cs_real_t recoi[6], recoj[6];
1910
1911 cs_i_compute_quantities_tensor(bldfrp,
1912 diipf,
1913 djjpf,
1914 gradi,
1915 gradj,
1916 pi,
1917 pj,
1918 recoi,
1919 recoj,
1920 pip,
1921 pjp);
1922
1923 cs_upwind_f_val_tensor(pi, pif);
1924 cs_upwind_f_val_tensor(pj, pjf);
1925
1926 }
1927
1928 /*----------------------------------------------------------------------------*/
1929 /*!
1930 * \brief Handle preparation of internal face values for the fluxes computation
1931 * in case of a steady algorithm and without enabling slope tests.
1932 *
1933 * \param[in] bldfrp reconstruction blending factor
1934 * \param[in] ischcp second order convection scheme flag
1935 * \param[in] relaxp relaxation coefficient
1936 * \param[in] blencp proportion of second order scheme,
1937 * (1-blencp) is the proportion of upwind.
1938 * \param[in] weight geometrical weight
1939 * \param[in] cell_ceni center of gravity coordinates of cell i
1940 * \param[in] cell_cenj center of gravity coordinates of cell i
1941 * \param[in] i_face_cog center of gravity coordinates of face ij
1942 * \param[in] diipf distance II'
1943 * \param[in] djjpf distance JJ'
1944 * \param[in] gradi gradient at cell i
1945 * \param[in] gradj gradient at cell j
1946 * \param[in] gradupi gradient upwind at cell i
1947 * \param[in] gradupj gradient upwind at cell j
1948 * \param[in] pi value at cell i
1949 * \param[in] pj value at cell j
1950 * \param[in] pia old value at cell i
1951 * \param[in] pja old value at cell j
1952 * \param[out] pifri contribution of i to flux from i to j
1953 * \param[out] pifrj contribution of i to flux from j to i
1954 * \param[out] pjfri contribution of j to flux from i to j
1955 * \param[out] pjfrj contribution of j to flux from j to i
1956 * \param[out] pip reconstructed value at cell i
1957 * \param[out] pjp reconstructed value at cell j
1958 * \param[out] pipr relaxed reconstructed value at cell i
1959 * \param[out] pjpr relaxed reconstructed value at cell j
1960 */
1961 /*----------------------------------------------------------------------------*/
1962
1963 inline static void
cs_i_cd_steady(const cs_real_t bldfrp,const int ischcp,const double relaxp,const double blencp,const cs_real_t weight,const cs_real_t cell_ceni[3],const cs_real_t cell_cenj[3],const cs_real_t i_face_cog[3],const cs_real_t diipf[3],const cs_real_t djjpf[3],const cs_real_t gradi[3],const cs_real_t gradj[3],const cs_real_t gradupi[3],const cs_real_t gradupj[3],const cs_real_t pi,const cs_real_t pj,const cs_real_t pia,const cs_real_t pja,cs_real_t * pifri,cs_real_t * pifrj,cs_real_t * pjfri,cs_real_t * pjfrj,cs_real_t * pip,cs_real_t * pjp,cs_real_t * pipr,cs_real_t * pjpr)1964 cs_i_cd_steady(const cs_real_t bldfrp,
1965 const int ischcp,
1966 const double relaxp,
1967 const double blencp,
1968 const cs_real_t weight,
1969 const cs_real_t cell_ceni[3],
1970 const cs_real_t cell_cenj[3],
1971 const cs_real_t i_face_cog[3],
1972 const cs_real_t diipf[3],
1973 const cs_real_t djjpf[3],
1974 const cs_real_t gradi[3],
1975 const cs_real_t gradj[3],
1976 const cs_real_t gradupi[3],
1977 const cs_real_t gradupj[3],
1978 const cs_real_t pi,
1979 const cs_real_t pj,
1980 const cs_real_t pia,
1981 const cs_real_t pja,
1982 cs_real_t *pifri,
1983 cs_real_t *pifrj,
1984 cs_real_t *pjfri,
1985 cs_real_t *pjfrj,
1986 cs_real_t *pip,
1987 cs_real_t *pjp,
1988 cs_real_t *pipr,
1989 cs_real_t *pjpr)
1990 {
1991 cs_real_t pir, pjr;
1992 cs_real_t recoi, recoj;
1993
1994 cs_i_compute_quantities(bldfrp,
1995 diipf,
1996 djjpf,
1997 gradi,
1998 gradj,
1999 pi,
2000 pj,
2001 &recoi,
2002 &recoj,
2003 pip,
2004 pjp);
2005
2006 cs_i_relax_c_val(relaxp,
2007 pia,
2008 pja,
2009 recoi,
2010 recoj,
2011 pi,
2012 pj,
2013 &pir,
2014 &pjr,
2015 pipr,
2016 pjpr);
2017
2018 if (ischcp == 1) {
2019
2020 /* Centered
2021 --------*/
2022
2023 cs_centered_f_val(weight,
2024 *pip,
2025 *pjpr,
2026 pifrj);
2027 cs_centered_f_val(weight,
2028 *pipr,
2029 *pjp,
2030 pifri);
2031 cs_centered_f_val(weight,
2032 *pipr,
2033 *pjp,
2034 pjfri);
2035 cs_centered_f_val(weight,
2036 *pip,
2037 *pjpr,
2038 pjfrj);
2039
2040 } else if (ischcp == 0) {
2041
2042 /* Original SOLU
2043 --------------*/
2044
2045 cs_solu_f_val(cell_ceni,
2046 i_face_cog,
2047 gradi,
2048 pi,
2049 pifrj);
2050 cs_solu_f_val(cell_ceni,
2051 i_face_cog,
2052 gradi,
2053 pir,
2054 pifri);
2055 cs_solu_f_val(cell_cenj,
2056 i_face_cog,
2057 gradj,
2058 pj,
2059 pjfri);
2060 cs_solu_f_val(cell_cenj,
2061 i_face_cog,
2062 gradj,
2063 pjr,
2064 pjfrj);
2065
2066 } else {
2067
2068 /* SOLU
2069 ----*/
2070
2071 cs_solu_f_val(cell_ceni,
2072 i_face_cog,
2073 gradupi,
2074 pi,
2075 pifrj);
2076 cs_solu_f_val(cell_ceni,
2077 i_face_cog,
2078 gradupi,
2079 pir,
2080 pifri);
2081 cs_solu_f_val(cell_cenj,
2082 i_face_cog,
2083 gradupj,
2084 pj,
2085 pjfri);
2086 cs_solu_f_val(cell_cenj,
2087 i_face_cog,
2088 gradupj,
2089 pjr,
2090 pjfrj);
2091
2092 }
2093
2094 /* Blending
2095 --------*/
2096
2097 cs_blend_f_val(blencp,
2098 pi,
2099 pifrj);
2100 cs_blend_f_val(blencp,
2101 pir,
2102 pifri);
2103 cs_blend_f_val(blencp,
2104 pj,
2105 pjfri);
2106 cs_blend_f_val(blencp,
2107 pjr,
2108 pjfrj);
2109 }
2110
2111 /*----------------------------------------------------------------------------*/
2112 /*!
2113 * \brief Handle preparation of internal face values for the fluxes computation
2114 * in case of a steady algorithm and without enabling slope tests.
2115 *
2116 * \param[in] bldfrp reconstruction blending factor
2117 * \param[in] ischcp second order convection scheme flag
2118 * \param[in] relaxp relaxation coefficient
2119 * \param[in] blencp proportion of second order scheme,
2120 * (1-blencp) is the proportion of upwind.
2121 * \param[in] weight geometrical weight
2122 * \param[in] cell_ceni center of gravity coordinates of cell i
2123 * \param[in] cell_cenj center of gravity coordinates of cell i
2124 * \param[in] i_face_cog center of gravity coordinates of face ij
2125 * \param[in] diipf distance II'
2126 * \param[in] djjpf distance JJ'
2127 * \param[in] gradi gradient at cell i
2128 * \param[in] gradj gradient at cell j
2129 * \param[in] pi value at cell i
2130 * \param[in] pj value at cell j
2131 * \param[in] pia old value at cell i
2132 * \param[in] pja old value at cell j
2133 * \param[out] pifri contribution of i to flux from i to j
2134 * \param[out] pifrj contribution of i to flux from j to i
2135 * \param[out] pjfri contribution of j to flux from i to j
2136 * \param[out] pjfrj contribution of j to flux from j to i
2137 * \param[out] pip reconstructed value at cell i
2138 * \param[out] pjp reconstructed value at cell j
2139 * \param[out] pipr relaxed reconstructed value at cell i
2140 * \param[out] pjpr relaxed reconstructed value at cell j
2141 */
2142 /*----------------------------------------------------------------------------*/
2143
2144 inline static void
cs_i_cd_steady_vector(const cs_real_t bldfrp,const int ischcp,const double relaxp,const double blencp,const cs_real_t weight,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_33_t gradi,const cs_real_33_t gradj,const cs_real_3_t pi,const cs_real_3_t pj,const cs_real_3_t pia,const cs_real_3_t pja,cs_real_t pifri[3],cs_real_t pifrj[3],cs_real_t pjfri[3],cs_real_t pjfrj[3],cs_real_t pip[3],cs_real_t pjp[3],cs_real_t pipr[3],cs_real_t pjpr[3])2145 cs_i_cd_steady_vector(const cs_real_t bldfrp,
2146 const int ischcp,
2147 const double relaxp,
2148 const double blencp,
2149 const cs_real_t weight,
2150 const cs_real_3_t cell_ceni,
2151 const cs_real_3_t cell_cenj,
2152 const cs_real_3_t i_face_cog,
2153 const cs_real_3_t diipf,
2154 const cs_real_3_t djjpf,
2155 const cs_real_33_t gradi,
2156 const cs_real_33_t gradj,
2157 const cs_real_3_t pi,
2158 const cs_real_3_t pj,
2159 const cs_real_3_t pia,
2160 const cs_real_3_t pja,
2161 cs_real_t pifri[3],
2162 cs_real_t pifrj[3],
2163 cs_real_t pjfri[3],
2164 cs_real_t pjfrj[3],
2165 cs_real_t pip[3],
2166 cs_real_t pjp[3],
2167 cs_real_t pipr[3],
2168 cs_real_t pjpr[3])
2169 {
2170 cs_real_t pir[3], pjr[3];
2171 cs_real_t recoi[3], recoj[3];
2172
2173 cs_i_compute_quantities_vector(bldfrp,
2174 diipf,
2175 djjpf,
2176 gradi,
2177 gradj,
2178 pi,
2179 pj,
2180 recoi,
2181 recoj,
2182 pip,
2183 pjp);
2184
2185 cs_i_relax_c_val_vector(relaxp,
2186 pia,
2187 pja,
2188 recoi,
2189 recoj,
2190 pi,
2191 pj,
2192 pir,
2193 pjr,
2194 pipr,
2195 pjpr);
2196
2197 if (ischcp == 1) {
2198
2199 /* Centered
2200 --------*/
2201
2202 cs_centered_f_val_vector(weight,
2203 pip,
2204 pjpr,
2205 pifrj);
2206 cs_centered_f_val_vector(weight,
2207 pipr,
2208 pjp,
2209 pifri);
2210 cs_centered_f_val_vector(weight,
2211 pipr,
2212 pjp,
2213 pjfri);
2214 cs_centered_f_val_vector(weight,
2215 pip,
2216 pjpr,
2217 pjfrj);
2218
2219 } else {
2220
2221 /* Second order
2222 ------------*/
2223
2224 cs_solu_f_val_vector(cell_ceni,
2225 i_face_cog,
2226 gradi,
2227 pi,
2228 pifrj);
2229 cs_solu_f_val_vector(cell_ceni,
2230 i_face_cog,
2231 gradi,
2232 pir,
2233 pifri);
2234 cs_solu_f_val_vector(cell_cenj,
2235 i_face_cog,
2236 gradj,
2237 pj,
2238 pjfri);
2239 cs_solu_f_val_vector(cell_cenj,
2240 i_face_cog,
2241 gradj,
2242 pjr,
2243 pjfrj);
2244
2245 }
2246
2247 /* Blending
2248 --------*/
2249 cs_blend_f_val_vector(blencp,
2250 pi,
2251 pifrj);
2252 cs_blend_f_val_vector(blencp,
2253 pir,
2254 pifri);
2255 cs_blend_f_val_vector(blencp,
2256 pj,
2257 pjfri);
2258 cs_blend_f_val_vector(blencp,
2259 pjr,
2260 pjfrj);
2261
2262 }
2263
2264 /*----------------------------------------------------------------------------*/
2265 /*!
2266 * \brief Handle preparation of internal face values for the fluxes computation
2267 * in case of a steady algorithm and without enabling slope tests.
2268 *
2269 * \param[in] bldfrp reconstruction blending factor
2270 * \param[in] ischcp second order convection scheme flag
2271 * \param[in] relaxp relaxation coefficient
2272 * \param[in] blencp proportion of second order scheme,
2273 * (1-blencp) is the proportion of upwind.
2274 * \param[in] weight geometrical weight
2275 * \param[in] cell_ceni center of gravity coordinates of cell i
2276 * \param[in] cell_cenj center of gravity coordinates of cell i
2277 * \param[in] i_face_cog center of gravity coordinates of face ij
2278 * \param[in] diipf distance II'
2279 * \param[in] djjpf distance JJ'
2280 * \param[in] gradi gradient at cell i
2281 * \param[in] gradj gradient at cell j
2282 * \param[in] pi value at cell i
2283 * \param[in] pj value at cell j
2284 * \param[in] pia old value at cell i
2285 * \param[in] pja old value at cell j
2286 * \param[out] pifri contribution of i to flux from i to j
2287 * \param[out] pifrj contribution of i to flux from j to i
2288 * \param[out] pjfri contribution of j to flux from i to j
2289 * \param[out] pjfrj contribution of j to flux from j to i
2290 * \param[out] pip reconstructed value at cell i
2291 * \param[out] pjp reconstructed value at cell j
2292 * \param[out] pipr relaxed reconstructed value at cell i
2293 * \param[out] pjpr relaxed reconstructed value at cell j
2294 */
2295 /*----------------------------------------------------------------------------*/
2296
2297 inline static void
cs_i_cd_steady_tensor(const cs_real_t bldfrp,const int ischcp,const double relaxp,const double blencp,const cs_real_t weight,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_63_t gradi,const cs_real_63_t gradj,const cs_real_6_t pi,const cs_real_6_t pj,const cs_real_6_t pia,const cs_real_6_t pja,cs_real_t pifri[6],cs_real_t pifrj[6],cs_real_t pjfri[6],cs_real_t pjfrj[6],cs_real_t pip[6],cs_real_t pjp[6],cs_real_t pipr[6],cs_real_t pjpr[6])2298 cs_i_cd_steady_tensor(const cs_real_t bldfrp,
2299 const int ischcp,
2300 const double relaxp,
2301 const double blencp,
2302 const cs_real_t weight,
2303 const cs_real_3_t cell_ceni,
2304 const cs_real_3_t cell_cenj,
2305 const cs_real_3_t i_face_cog,
2306 const cs_real_3_t diipf,
2307 const cs_real_3_t djjpf,
2308 const cs_real_63_t gradi,
2309 const cs_real_63_t gradj,
2310 const cs_real_6_t pi,
2311 const cs_real_6_t pj,
2312 const cs_real_6_t pia,
2313 const cs_real_6_t pja,
2314 cs_real_t pifri[6],
2315 cs_real_t pifrj[6],
2316 cs_real_t pjfri[6],
2317 cs_real_t pjfrj[6],
2318 cs_real_t pip[6],
2319 cs_real_t pjp[6],
2320 cs_real_t pipr[6],
2321 cs_real_t pjpr[6])
2322
2323 {
2324 cs_real_t pir[6], pjr[6];
2325 cs_real_t recoi[6], recoj[6];
2326
2327 cs_i_compute_quantities_tensor(bldfrp,
2328 diipf,
2329 djjpf,
2330 gradi,
2331 gradj,
2332 pi,
2333 pj,
2334 recoi,
2335 recoj,
2336 pip,
2337 pjp);
2338
2339 cs_i_relax_c_val_tensor(relaxp,
2340 pia,
2341 pja,
2342 recoi,
2343 recoj,
2344 pi,
2345 pj,
2346 pir,
2347 pjr,
2348 pipr,
2349 pjpr);
2350
2351 if (ischcp == 1) {
2352
2353 /* Centered
2354 --------*/
2355
2356 cs_centered_f_val_tensor(weight,
2357 pip,
2358 pjpr,
2359 pifrj);
2360 cs_centered_f_val_tensor(weight,
2361 pipr,
2362 pjp,
2363 pifri);
2364 cs_centered_f_val_tensor(weight,
2365 pipr,
2366 pjp,
2367 pjfri);
2368 cs_centered_f_val_tensor(weight,
2369 pip,
2370 pjpr,
2371 pjfrj);
2372
2373 } else {
2374
2375 /* Second order
2376 ------------*/
2377
2378 cs_solu_f_val_tensor(cell_ceni,
2379 i_face_cog,
2380 gradi,
2381 pi,
2382 pifrj);
2383 cs_solu_f_val_tensor(cell_ceni,
2384 i_face_cog,
2385 gradi,
2386 pir,
2387 pifri);
2388 cs_solu_f_val_tensor(cell_cenj,
2389 i_face_cog,
2390 gradj,
2391 pj,
2392 pjfri);
2393 cs_solu_f_val_tensor(cell_cenj,
2394 i_face_cog,
2395 gradj,
2396 pjr,
2397 pjfrj);
2398
2399 }
2400
2401 /* Blending
2402 --------*/
2403
2404 cs_blend_f_val_tensor(blencp,
2405 pi,
2406 pifrj);
2407 cs_blend_f_val_tensor(blencp,
2408 pir,
2409 pifri);
2410 cs_blend_f_val_tensor(blencp,
2411 pj,
2412 pjfri);
2413 cs_blend_f_val_tensor(blencp,
2414 pjr,
2415 pjfrj);
2416 }
2417
2418 /*----------------------------------------------------------------------------*/
2419 /*!
2420 * \brief Handle preparation of internal face values for the fluxes computation
2421 * in case of a unsteady algorithm and without enabling slope tests.
2422 *
2423 * \param[in] bldfrp reconstruction blending factor
2424 * \param[in] ischcp second order convection scheme flag
2425 * \param[in] blencp proportion of second order scheme,
2426 * (1-blencp) is the proportion of upwind.
2427 * \param[in] weight geometrical weight
2428 * \param[in] cell_ceni center of gravity coordinates of cell i
2429 * \param[in] cell_cenj center of gravity coordinates of cell i
2430 * \param[in] i_face_cog center of gravity coordinates of face ij
2431 * \param[in] hybrid_blend_i blending factor between SOLU and centered
2432 * \param[in] hybrid_blend_j blending factor between SOLU and centered
2433 * \param[in] diipf distance II'
2434 * \param[in] djjpf distance JJ'
2435 * \param[in] gradi gradient at cell i
2436 * \param[in] gradj gradient at cell j
2437 * \param[in] gradupi upwind gradient at cell i
2438 * \param[in] gradupj upwind gradient at cell j
2439 * \param[in] gradi gradient at cell i
2440 * \param[in] gradj gradient at cell j
2441 * \param[in] pi value at cell i
2442 * \param[in] pj value at cell j
2443 * \param[out] pif contribution of i to flux from i to j
2444 * \param[out] pjf contribution of j to flux from i to j
2445 * \param[out] pip reconstructed value at cell i
2446 * \param[out] pjp reconstructed value at cell j
2447 */
2448 /*----------------------------------------------------------------------------*/
2449
2450 inline static void
cs_i_cd_unsteady(const cs_real_t bldfrp,const int ischcp,const double blencp,const cs_real_t weight,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_cog,const cs_real_t hybrid_blend_i,const cs_real_t hybrid_blend_j,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_3_t gradi,const cs_real_3_t gradj,const cs_real_3_t gradupi,const cs_real_3_t gradupj,const cs_real_t pi,const cs_real_t pj,cs_real_t * pif,cs_real_t * pjf,cs_real_t * pip,cs_real_t * pjp)2451 cs_i_cd_unsteady(const cs_real_t bldfrp,
2452 const int ischcp,
2453 const double blencp,
2454 const cs_real_t weight,
2455 const cs_real_3_t cell_ceni,
2456 const cs_real_3_t cell_cenj,
2457 const cs_real_3_t i_face_cog,
2458 const cs_real_t hybrid_blend_i,
2459 const cs_real_t hybrid_blend_j,
2460 const cs_real_3_t diipf,
2461 const cs_real_3_t djjpf,
2462 const cs_real_3_t gradi,
2463 const cs_real_3_t gradj,
2464 const cs_real_3_t gradupi,
2465 const cs_real_3_t gradupj,
2466 const cs_real_t pi,
2467 const cs_real_t pj,
2468 cs_real_t *pif,
2469 cs_real_t *pjf,
2470 cs_real_t *pip,
2471 cs_real_t *pjp)
2472 {
2473 cs_real_t recoi, recoj;
2474
2475 cs_i_compute_quantities(bldfrp,
2476 diipf,
2477 djjpf,
2478 gradi,
2479 gradj,
2480 pi,
2481 pj,
2482 &recoi,
2483 &recoj,
2484 pip,
2485 pjp);
2486
2487
2488 if (ischcp == 1) {
2489
2490 /* Centered
2491 --------*/
2492
2493 cs_centered_f_val(weight,
2494 *pip,
2495 *pjp,
2496 pif);
2497 cs_centered_f_val(weight,
2498 *pip,
2499 *pjp,
2500 pjf);
2501
2502 } else if (ischcp == 0) {
2503
2504 /* Legacy SOLU
2505 -----------*/
2506
2507 cs_solu_f_val(cell_ceni,
2508 i_face_cog,
2509 gradi,
2510 pi,
2511 pif);
2512 cs_solu_f_val(cell_cenj,
2513 i_face_cog,
2514 gradj,
2515 pj,
2516 pjf);
2517
2518 } else if (ischcp == 3) {
2519
2520 /* Centered
2521 --------*/
2522
2523 cs_centered_f_val(weight,
2524 *pip,
2525 *pjp,
2526 pif);
2527 cs_centered_f_val(weight,
2528 *pip,
2529 *pjp,
2530 pjf);
2531
2532 /* Legacy SOLU
2533 -----------*/
2534 cs_real_t pif_up, pjf_up;
2535 cs_real_t hybrid_blend_interp;
2536
2537 cs_solu_f_val(cell_ceni,
2538 i_face_cog,
2539 gradi,
2540 pi,
2541 &pif_up);
2542 cs_solu_f_val(cell_cenj,
2543 i_face_cog,
2544 gradj,
2545 pj,
2546 &pjf_up);
2547
2548 hybrid_blend_interp = fmin(hybrid_blend_i,hybrid_blend_j);
2549 *pif = hybrid_blend_interp*(*pif) + (1. - hybrid_blend_interp)*pif_up;
2550 *pjf = hybrid_blend_interp*(*pjf) + (1. - hybrid_blend_interp)*pjf_up;
2551
2552 } else {
2553
2554 /* SOLU
2555 ----*/
2556
2557 cs_solu_f_val(cell_ceni,
2558 i_face_cog,
2559 gradupi,
2560 pi,
2561 pif);
2562 cs_solu_f_val(cell_cenj,
2563 i_face_cog,
2564 gradupj,
2565 pj,
2566 pjf);
2567
2568 }
2569
2570
2571 /* Blending
2572 --------*/
2573
2574 cs_blend_f_val(blencp,
2575 pi,
2576 pif);
2577 cs_blend_f_val(blencp,
2578 pj,
2579 pjf);
2580 }
2581
2582 /*----------------------------------------------------------------------------*/
2583 /*!
2584 * \brief Handle preparation of internal face values for the fluxes computation
2585 * in case of an unsteady algorithm and without enabling slope tests.
2586 *
2587 * \param[in] bldfrp reconstruction blending factor
2588 * \param[in] ischcp second order convection scheme flag
2589 * \param[in] blencp proportion of second order scheme,
2590 * (1-blencp) is the proportion of upwind.
2591 * \param[in] weight geometrical weight
2592 * \param[in] cell_ceni center of gravity coordinates of cell i
2593 * \param[in] cell_cenj center of gravity coordinates of cell i
2594 * \param[in] i_face_cog center of gravity coordinates of face ij
2595 * \param[in] hybrid_blend_i blending factor between SOLU and centered
2596 * \param[in] hybrid_blend_j blending factor between SOLU and centered
2597 * \param[in] diipf distance II'
2598 * \param[in] djjpf distance JJ'
2599 * \param[in] gradi gradient at cell i
2600 * \param[in] gradj gradient at cell j
2601 * \param[in] pi value at cell i
2602 * \param[in] pj value at cell j
2603 * \param[out] pif contribution of i to flux from i to j
2604 * \param[out] pjf contribution of j to flux from j to i
2605 * \param[out] pip reconstructed value at cell i
2606 * \param[out] pjp reconstructed value at cell j
2607 */
2608 /*----------------------------------------------------------------------------*/
2609
2610 inline static void
cs_i_cd_unsteady_vector(const cs_real_t bldfrp,const int ischcp,const double blencp,const cs_real_t weight,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_cog,const cs_real_t hybrid_blend_i,const cs_real_t hybrid_blend_j,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_33_t gradi,const cs_real_33_t gradj,const cs_real_3_t pi,const cs_real_3_t pj,cs_real_t pif[3],cs_real_t pjf[3],cs_real_t pip[3],cs_real_t pjp[3])2611 cs_i_cd_unsteady_vector(const cs_real_t bldfrp,
2612 const int ischcp,
2613 const double blencp,
2614 const cs_real_t weight,
2615 const cs_real_3_t cell_ceni,
2616 const cs_real_3_t cell_cenj,
2617 const cs_real_3_t i_face_cog,
2618 const cs_real_t hybrid_blend_i,
2619 const cs_real_t hybrid_blend_j,
2620 const cs_real_3_t diipf,
2621 const cs_real_3_t djjpf,
2622 const cs_real_33_t gradi,
2623 const cs_real_33_t gradj,
2624 const cs_real_3_t pi,
2625 const cs_real_3_t pj,
2626 cs_real_t pif[3],
2627 cs_real_t pjf[3],
2628 cs_real_t pip[3],
2629 cs_real_t pjp[3])
2630
2631 {
2632 cs_real_t recoi[3], recoj[3];
2633
2634 cs_i_compute_quantities_vector(bldfrp,
2635 diipf,
2636 djjpf,
2637 gradi,
2638 gradj,
2639 pi,
2640 pj,
2641 recoi,
2642 recoj,
2643 pip,
2644 pjp);
2645
2646 if (ischcp == 1) {
2647
2648 /* Centered
2649 --------*/
2650
2651 cs_centered_f_val_vector(weight,
2652 pip,
2653 pjp,
2654 pif);
2655 cs_centered_f_val_vector(weight,
2656 pip,
2657 pjp,
2658 pjf);
2659 } else if (ischcp == 3) {
2660
2661 /* Centered
2662 --------*/
2663
2664 cs_centered_f_val_vector(weight,
2665 pip,
2666 pjp,
2667 pif);
2668 cs_centered_f_val_vector(weight,
2669 pip,
2670 pjp,
2671 pjf);
2672
2673 /* SOLU
2674 -----*/
2675 cs_real_t pif_up[3], pjf_up[3];
2676 cs_real_t hybrid_blend_interp;
2677
2678 cs_solu_f_val_vector(cell_ceni,
2679 i_face_cog,
2680 gradi,
2681 pi,
2682 pif_up);
2683 cs_solu_f_val_vector(cell_cenj,
2684 i_face_cog,
2685 gradj,
2686 pj,
2687 pjf_up);
2688
2689 hybrid_blend_interp = fmin(hybrid_blend_i,hybrid_blend_j);
2690 for (int isou = 0; isou < 3; isou++) {
2691 pif[isou] = hybrid_blend_interp *pif[isou]
2692 + (1. - hybrid_blend_interp)*pif_up[isou];
2693 pjf[isou] = hybrid_blend_interp *pjf[isou]
2694 + (1. - hybrid_blend_interp)*pjf_up[isou];
2695 }
2696 } else {
2697
2698 /* Second order
2699 ------------*/
2700
2701 cs_solu_f_val_vector(cell_ceni,
2702 i_face_cog,
2703 gradi,
2704 pi,
2705 pif);
2706 cs_solu_f_val_vector(cell_cenj,
2707 i_face_cog,
2708 gradj,
2709 pj,
2710 pjf);
2711
2712 }
2713
2714 /* Blending
2715 --------*/
2716
2717 cs_blend_f_val_vector(blencp,
2718 pi,
2719 pif);
2720 cs_blend_f_val_vector(blencp,
2721 pj,
2722 pjf);
2723
2724 }
2725
2726 /*----------------------------------------------------------------------------*/
2727 /*!
2728 * \brief Handle preparation of internal face values for the fluxes computation
2729 * in case of an unsteady algorithm and without enabling slope tests.
2730 *
2731 * \param[in] bldfrp reconstruction blending factor
2732 * \param[in] ischcp second order convection scheme flag
2733 * \param[in] blencp proportion of second order scheme,
2734 * (1-blencp) is the proportion of upwind.
2735 * \param[in] weight geometrical weight
2736 * \param[in] cell_ceni center of gravity coordinates of cell i
2737 * \param[in] cell_cenj center of gravity coordinates of cell i
2738 * \param[in] i_face_cog center of gravity coordinates of face ij
2739 * \param[in] diipf distance II'
2740 * \param[in] djjpf distance JJ'
2741 * \param[in] gradi gradient at cell i
2742 * \param[in] gradj gradient at cell j
2743 * \param[in] pi value at cell i
2744 * \param[in] pj value at cell j
2745 * \param[out] pif contribution of i to flux from i to j
2746 * \param[out] pjf contribution of j to flux from i to j
2747 * \param[out] pip reconstructed value at cell i
2748 * \param[out] pjp reconstructed value at cell j
2749 */
2750 /*----------------------------------------------------------------------------*/
2751
2752 inline static void
cs_i_cd_unsteady_tensor(const cs_real_t bldfrp,const int ischcp,const double blencp,const cs_real_t weight,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_63_t gradi,const cs_real_63_t gradj,const cs_real_6_t pi,const cs_real_6_t pj,cs_real_t pif[6],cs_real_t pjf[6],cs_real_t pip[6],cs_real_t pjp[6])2753 cs_i_cd_unsteady_tensor(const cs_real_t bldfrp,
2754 const int ischcp,
2755 const double blencp,
2756 const cs_real_t weight,
2757 const cs_real_3_t cell_ceni,
2758 const cs_real_3_t cell_cenj,
2759 const cs_real_3_t i_face_cog,
2760 const cs_real_3_t diipf,
2761 const cs_real_3_t djjpf,
2762 const cs_real_63_t gradi,
2763 const cs_real_63_t gradj,
2764 const cs_real_6_t pi,
2765 const cs_real_6_t pj,
2766 cs_real_t pif[6],
2767 cs_real_t pjf[6],
2768 cs_real_t pip[6],
2769 cs_real_t pjp[6])
2770
2771 {
2772 cs_real_t recoi[6], recoj[6];
2773
2774 cs_i_compute_quantities_tensor(bldfrp,
2775 diipf,
2776 djjpf,
2777 gradi,
2778 gradj,
2779 pi,
2780 pj,
2781 recoi,
2782 recoj,
2783 pip,
2784 pjp);
2785
2786 if (ischcp == 1) {
2787
2788 /* Centered
2789 --------*/
2790
2791 cs_centered_f_val_tensor(weight,
2792 pip,
2793 pjp,
2794 pif);
2795 cs_centered_f_val_tensor(weight,
2796 pip,
2797 pjp,
2798 pjf);
2799
2800 } else {
2801
2802 /* Second order
2803 ------------*/
2804
2805 cs_solu_f_val_tensor(cell_ceni,
2806 i_face_cog,
2807 gradi,
2808 pi,
2809 pif);
2810 cs_solu_f_val_tensor(cell_cenj,
2811 i_face_cog,
2812 gradj,
2813 pj,
2814 pjf);
2815
2816 }
2817
2818 /* Blending
2819 --------*/
2820
2821 cs_blend_f_val_tensor(blencp,
2822 pi,
2823 pif);
2824 cs_blend_f_val_tensor(blencp,
2825 pj,
2826 pjf);
2827
2828 }
2829
2830 /*----------------------------------------------------------------------------*/
2831 /*!
2832 * \brief Handle preparation of internal face values for the fluxes computation
2833 * in case of a steady algorithm and using slope tests.
2834 *
2835 * \param[out] upwind_switch slope test result
2836 * \param[in] iconvp convection flag
2837 * \param[in] bldfrp reconstruction blending factor
2838 * \param[in] ischcp second order convection scheme flag
2839 * \param[in] relaxp relaxation coefficient
2840 * \param[in] blencp proportion of second order scheme,
2841 * (1-blencp) is the proportion of upwind.
2842 * \param[in] blend_st proportion of centered or SOLU scheme,
2843 * when the slope test is activated
2844 * (1-blend_st) is the proportion of upwind.
2845 * \param[in] weight geometrical weight
2846 * \param[in] i_dist distance IJ.Nij
2847 * \param[in] i_face_surf face surface
2848 * \param[in] cell_ceni center of gravity coordinates of cell i
2849 * \param[in] cell_cenj center of gravity coordinates of cell j
2850 * \param[in] i_face_normal face normal
2851 * \param[in] i_face_cog center of gravity coordinates of face ij
2852 * \param[in] diipf distance II'
2853 * \param[in] djjpf distance JJ'
2854 * \param[in] i_massflux mass flux at face ij
2855 * \param[in] gradi gradient at cell i
2856 * \param[in] gradj gradient at cell j
2857 * \param[in] gradupi upwind gradient at cell i
2858 * \param[in] gradupj upwind gradient at cell j
2859 * \param[in] gradsti slope test gradient at cell i
2860 * \param[in] gradstj slope test gradient at cell j
2861 * \param[in] pi value at cell i
2862 * \param[in] pj value at cell j
2863 * \param[in] pia old value at cell i
2864 * \param[in] pja old value at cell j
2865 * \param[out] pifri contribution of i to flux from i to j
2866 * \param[out] pifrj contribution of i to flux from j to i
2867 * \param[out] pjfri contribution of j to flux from i to j
2868 * \param[out] pjfrj contribution of j to flux from j to i
2869 * \param[out] pip reconstructed value at cell i
2870 * \param[out] pjp reconstructed value at cell j
2871 * \param[out] pipr relaxed reconstructed value at cell i
2872 * \param[out] pjpr relaxed reconstructed value at cell j
2873 */
2874 /*----------------------------------------------------------------------------*/
2875
2876 inline static void
cs_i_cd_steady_slope_test(bool * upwind_switch,const int iconvp,const cs_real_t bldfrp,const int ischcp,const double relaxp,const double blencp,const double blend_st,const cs_real_t weight,const cs_real_t i_dist,const cs_real_t i_face_surf,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_normal,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_t i_massflux,const cs_real_3_t gradi,const cs_real_3_t gradj,const cs_real_3_t gradupi,const cs_real_3_t gradupj,const cs_real_3_t gradsti,const cs_real_3_t gradstj,const cs_real_t pi,const cs_real_t pj,const cs_real_t pia,const cs_real_t pja,cs_real_t * pifri,cs_real_t * pifrj,cs_real_t * pjfri,cs_real_t * pjfrj,cs_real_t * pip,cs_real_t * pjp,cs_real_t * pipr,cs_real_t * pjpr)2877 cs_i_cd_steady_slope_test(bool *upwind_switch,
2878 const int iconvp,
2879 const cs_real_t bldfrp,
2880 const int ischcp,
2881 const double relaxp,
2882 const double blencp,
2883 const double blend_st,
2884 const cs_real_t weight,
2885 const cs_real_t i_dist,
2886 const cs_real_t i_face_surf,
2887 const cs_real_3_t cell_ceni,
2888 const cs_real_3_t cell_cenj,
2889 const cs_real_3_t i_face_normal,
2890 const cs_real_3_t i_face_cog,
2891 const cs_real_3_t diipf,
2892 const cs_real_3_t djjpf,
2893 const cs_real_t i_massflux,
2894 const cs_real_3_t gradi,
2895 const cs_real_3_t gradj,
2896 const cs_real_3_t gradupi,
2897 const cs_real_3_t gradupj,
2898 const cs_real_3_t gradsti,
2899 const cs_real_3_t gradstj,
2900 const cs_real_t pi,
2901 const cs_real_t pj,
2902 const cs_real_t pia,
2903 const cs_real_t pja,
2904 cs_real_t *pifri,
2905 cs_real_t *pifrj,
2906 cs_real_t *pjfri,
2907 cs_real_t *pjfrj,
2908 cs_real_t *pip,
2909 cs_real_t *pjp,
2910 cs_real_t *pipr,
2911 cs_real_t *pjpr)
2912 {
2913 cs_real_t pir, pjr;
2914 cs_real_t recoi, recoj;
2915 cs_real_t testij, tesqck;
2916
2917 *upwind_switch = false;
2918
2919 cs_i_compute_quantities(bldfrp,
2920 diipf,
2921 djjpf,
2922 gradi,
2923 gradj,
2924 pi,
2925 pj,
2926 &recoi,
2927 &recoj,
2928 pip,
2929 pjp);
2930
2931 cs_i_relax_c_val(relaxp,
2932 pia,
2933 pja,
2934 recoi,
2935 recoj,
2936 pi,
2937 pj,
2938 &pir,
2939 &pjr,
2940 pipr,
2941 pjpr);
2942
2943 /* Convection slope test is needed only when iconv >0 */
2944 if (iconvp > 0) {
2945 cs_slope_test(pi,
2946 pj,
2947 i_dist,
2948 i_face_surf,
2949 i_face_normal,
2950 gradi,
2951 gradj,
2952 gradsti,
2953 gradstj,
2954 i_massflux,
2955 &testij,
2956 &tesqck);
2957
2958 if (ischcp==1) {
2959
2960 /* Centered
2961 --------*/
2962
2963 cs_centered_f_val(weight,
2964 *pip,
2965 *pjpr,
2966 pifrj);
2967 cs_centered_f_val(weight,
2968 *pipr,
2969 *pjp,
2970 pifri);
2971 cs_centered_f_val(weight,
2972 *pipr,
2973 *pjp,
2974 pjfri);
2975 cs_centered_f_val(weight,
2976 *pip,
2977 *pjpr,
2978 pjfrj);
2979
2980 } else if (ischcp == 0) {
2981
2982 /* Second order
2983 ------------*/
2984
2985 cs_solu_f_val(cell_ceni,
2986 i_face_cog,
2987 gradi,
2988 pi,
2989 pifrj);
2990 cs_solu_f_val(cell_ceni,
2991 i_face_cog,
2992 gradi,
2993 pir,
2994 pifri);
2995 cs_solu_f_val(cell_cenj,
2996 i_face_cog,
2997 gradj,
2998 pj,
2999 pjfri);
3000 cs_solu_f_val(cell_cenj,
3001 i_face_cog,
3002 gradj,
3003 pjr,
3004 pjfrj);
3005
3006 } else {
3007
3008 /* SOLU
3009 -----*/
3010
3011 cs_solu_f_val(cell_ceni,
3012 i_face_cog,
3013 gradupi,
3014 pi,
3015 pifrj);
3016 cs_solu_f_val(cell_ceni,
3017 i_face_cog,
3018 gradupi,
3019 pir,
3020 pifri);
3021 cs_solu_f_val(cell_cenj,
3022 i_face_cog,
3023 gradupj,
3024 pj,
3025 pjfri);
3026 cs_solu_f_val(cell_cenj,
3027 i_face_cog,
3028 gradupj,
3029 pjr,
3030 pjfrj);
3031 }
3032
3033
3034 /* Slope test: Pourcentage of upwind
3035 ----------------------------------*/
3036
3037 if (tesqck <= 0. || testij <= 0.) {
3038
3039 cs_blend_f_val(blend_st,
3040 pi,
3041 pifrj);
3042 cs_blend_f_val(blend_st,
3043 pir,
3044 pifri);
3045 cs_blend_f_val(blend_st,
3046 pj,
3047 pjfri);
3048 cs_blend_f_val(blend_st,
3049 pjr,
3050 pjfrj);
3051
3052 *upwind_switch = true;
3053
3054 }
3055
3056
3057 /* Blending
3058 --------*/
3059
3060 cs_blend_f_val(blencp,
3061 pi,
3062 pifrj);
3063 cs_blend_f_val(blencp,
3064 pir,
3065 pifri);
3066 cs_blend_f_val(blencp,
3067 pj,
3068 pjfri);
3069 cs_blend_f_val(blencp,
3070 pjr,
3071 pjfrj);
3072
3073 /* If iconv=0 p*fr* are useless */
3074 } else {
3075 cs_upwind_f_val(pi,
3076 pifrj);
3077 cs_upwind_f_val(pir,
3078 pifri);
3079 cs_upwind_f_val(pj,
3080 pjfri);
3081 cs_upwind_f_val(pjr,
3082 pjfrj);
3083 }
3084
3085 }
3086
3087 /*----------------------------------------------------------------------------*/
3088 /*!
3089 * \brief Handle preparation of internal face values for the fluxes computation
3090 * in case of a steady algorithm and using slope tests.
3091 *
3092 * \param[out] upwind_switch slope test result
3093 * \param[in] iconvp convection flag
3094 * \param[in] bldfrp reconstruction blending factor
3095 * \param[in] ischcp second order convection scheme flag
3096 * \param[in] relaxp relaxation coefficient
3097 * \param[in] blencp proportion of second order scheme,
3098 * (1-blencp) is the proportion of upwind.
3099 * \param[in] blend_st proportion of centered or SOLU scheme,
3100 * when the slope test is activated
3101 * (1-blend_st) is the proportion of upwind.
3102 * \param[in] weight geometrical weight
3103 * \param[in] i_dist distance IJ.Nij
3104 * \param[in] i_face_surf face surface
3105 * \param[in] cell_ceni center of gravity coordinates of cell i
3106 * \param[in] cell_cenj center of gravity coordinates of cell i
3107 * \param[in] i_face_normal face normal
3108 * \param[in] i_face_cog center of gravity coordinates of face ij
3109 * \param[in] diipf distance II'
3110 * \param[in] djjpf distance JJ'
3111 * \param[in] i_massflux mass flux at face ij
3112 * \param[in] gradi gradient at cell i
3113 * \param[in] gradj gradient at cell j
3114 * \param[in] grdpai upwind gradient at cell i
3115 * \param[in] grdpaj upwind gradient at cell j
3116 * \param[in] pi value at cell i
3117 * \param[in] pj value at cell j
3118 * \param[in] pia old value at cell i
3119 * \param[in] pja old value at cell j
3120 * \param[out] pifri contribution of i to flux from i to j
3121 * \param[out] pifrj contribution of i to flux from j to i
3122 * \param[out] pjfri contribution of j to flux from i to j
3123 * \param[out] pjfrj contribution of j to flux from j to i
3124 * \param[out] pip reconstructed value at cell i
3125 * \param[out] pjp reconstructed value at cell j
3126 * \param[out] pipr relaxed reconstructed value at cell i
3127 * \param[out] pjpr relaxed reconstructed value at cell j
3128 */
3129 /*----------------------------------------------------------------------------*/
3130
3131 inline static void
cs_i_cd_steady_slope_test_vector(bool * upwind_switch,const int iconvp,const cs_real_t bldfrp,const int ischcp,const double relaxp,const double blencp,const double blend_st,const cs_real_t weight,const cs_real_t i_dist,const cs_real_t i_face_surf,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_normal,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_t i_massflux,const cs_real_33_t gradi,const cs_real_33_t gradj,const cs_real_33_t grdpai,const cs_real_33_t grdpaj,const cs_real_3_t pi,const cs_real_3_t pj,const cs_real_3_t pia,const cs_real_3_t pja,cs_real_t pifri[3],cs_real_t pifrj[3],cs_real_t pjfri[3],cs_real_t pjfrj[3],cs_real_t pip[3],cs_real_t pjp[3],cs_real_t pipr[3],cs_real_t pjpr[3])3132 cs_i_cd_steady_slope_test_vector(bool *upwind_switch,
3133 const int iconvp,
3134 const cs_real_t bldfrp,
3135 const int ischcp,
3136 const double relaxp,
3137 const double blencp,
3138 const double blend_st,
3139 const cs_real_t weight,
3140 const cs_real_t i_dist,
3141 const cs_real_t i_face_surf,
3142 const cs_real_3_t cell_ceni,
3143 const cs_real_3_t cell_cenj,
3144 const cs_real_3_t i_face_normal,
3145 const cs_real_3_t i_face_cog,
3146 const cs_real_3_t diipf,
3147 const cs_real_3_t djjpf,
3148 const cs_real_t i_massflux,
3149 const cs_real_33_t gradi,
3150 const cs_real_33_t gradj,
3151 const cs_real_33_t grdpai,
3152 const cs_real_33_t grdpaj,
3153 const cs_real_3_t pi,
3154 const cs_real_3_t pj,
3155 const cs_real_3_t pia,
3156 const cs_real_3_t pja,
3157 cs_real_t pifri[3],
3158 cs_real_t pifrj[3],
3159 cs_real_t pjfri[3],
3160 cs_real_t pjfrj[3],
3161 cs_real_t pip[3],
3162 cs_real_t pjp[3],
3163 cs_real_t pipr[3],
3164 cs_real_t pjpr[3])
3165 {
3166 cs_real_t pir[3], pjr[3];
3167 cs_real_t recoi[3], recoj[3];
3168 cs_real_t testij, tesqck;
3169 int isou;
3170
3171 cs_i_compute_quantities_vector(bldfrp,
3172 diipf,
3173 djjpf,
3174 gradi,
3175 gradj,
3176 pi,
3177 pj,
3178 recoi,
3179 recoj,
3180 pip,
3181 pjp);
3182
3183 cs_i_relax_c_val_vector(relaxp,
3184 pia,
3185 pja,
3186 recoi,
3187 recoj,
3188 pi,
3189 pj,
3190 pir,
3191 pjr,
3192 pipr,
3193 pjpr);
3194
3195 /* Convection slope test is needed only when iconv >0 */
3196 if (iconvp > 0) {
3197 cs_slope_test_vector(pi,
3198 pj,
3199 i_dist,
3200 i_face_surf,
3201 i_face_normal,
3202 gradi,
3203 gradj,
3204 grdpai,
3205 grdpaj,
3206 i_massflux,
3207 &testij,
3208 &tesqck);
3209
3210 for (isou = 0; isou < 3; isou++) {
3211 if (ischcp==1) {
3212
3213 /* Centered
3214 --------*/
3215
3216 cs_centered_f_val(weight,
3217 pip[isou],
3218 pjpr[isou],
3219 &pifrj[isou]);
3220 cs_centered_f_val(weight,
3221 pipr[isou],
3222 pjp[isou],
3223 &pifri[isou]);
3224 cs_centered_f_val(weight,
3225 pipr[isou],
3226 pjp[isou],
3227 &pjfri[isou]);
3228 cs_centered_f_val(weight,
3229 pip[isou],
3230 pjpr[isou],
3231 &pjfrj[isou]);
3232
3233 } else {
3234
3235 /* Second order
3236 ------------*/
3237
3238 cs_solu_f_val(cell_ceni,
3239 i_face_cog,
3240 gradi[isou],
3241 pi[isou],
3242 &pifrj[isou]);
3243 cs_solu_f_val(cell_ceni,
3244 i_face_cog,
3245 gradi[isou],
3246 pir[isou],
3247 &pifri[isou]);
3248 cs_solu_f_val(cell_cenj,
3249 i_face_cog,
3250 gradj[isou],
3251 pj[isou],
3252 &pjfri[isou]);
3253 cs_solu_f_val(cell_cenj,
3254 i_face_cog,
3255 gradj[isou],
3256 pjr[isou],
3257 &pjfrj[isou]);
3258
3259 }
3260
3261 }
3262
3263 /* Slope test: Pourcentage of upwind
3264 ----------------------------------*/
3265
3266 if (tesqck <= 0. || testij <= 0.) {
3267 cs_blend_f_val_vector(blend_st,
3268 pi,
3269 pifrj);
3270 cs_blend_f_val_vector(blend_st,
3271 pir,
3272 pifri);
3273 cs_blend_f_val_vector(blend_st,
3274 pj,
3275 pjfri);
3276 cs_blend_f_val_vector(blend_st,
3277 pjr,
3278 pjfrj);
3279
3280 *upwind_switch = true;
3281 }
3282
3283
3284 /* Blending
3285 --------*/
3286 cs_blend_f_val_vector(blencp,
3287 pi,
3288 pifrj);
3289 cs_blend_f_val_vector(blencp,
3290 pir,
3291 pifri);
3292 cs_blend_f_val_vector(blencp,
3293 pj,
3294 pjfri);
3295 cs_blend_f_val_vector(blencp,
3296 pjr,
3297 pjfrj);
3298
3299 /* If iconv=0 p*fr* are useless */
3300 } else {
3301 for (isou = 0; isou < 3; isou++) {
3302 cs_upwind_f_val(pi[isou],
3303 &pifrj[isou]);
3304 cs_upwind_f_val(pir[isou],
3305 &pifri[isou]);
3306 cs_upwind_f_val(pj[isou],
3307 &pjfri[isou]);
3308 cs_upwind_f_val(pjr[isou],
3309 &pjfrj[isou]);
3310 }
3311 }
3312
3313 }
3314
3315 /*----------------------------------------------------------------------------*/
3316 /*!
3317 * \brief Handle preparation of internal face values for the fluxes computation
3318 * in case of a steady algorithm and using slope tests.
3319 *
3320 * \param[out] upwind_switch slope test result
3321 * \param[in] iconvp convection flag
3322 * \param[in] bldfrp reconstruction blending factor
3323 * \param[in] ischcp second order convection scheme flag
3324 * \param[in] relaxp relaxation coefficient
3325 * \param[in] blencp proportion of second order scheme,
3326 * (1-blencp) is the proportion of upwind.
3327 * \param[in] blend_st proportion of centered or SOLU scheme,
3328 * when the slope test is activated
3329 * (1-blend_st) is the proportion of upwind.
3330 * \param[in] weight geometrical weight
3331 * \param[in] i_dist distance IJ.Nij
3332 * \param[in] i_face_surf face surface
3333 * \param[in] cell_ceni center of gravity coordinates of cell i
3334 * \param[in] cell_cenj center of gravity coordinates of cell i
3335 * \param[in] i_face_normal face normal
3336 * \param[in] i_face_cog center of gravity coordinates of face ij
3337 * \param[in] diipf distance II'
3338 * \param[in] djjpf distance JJ'
3339 * \param[in] i_massflux mass flux at face ij
3340 * \param[in] gradi gradient at cell i
3341 * \param[in] gradj gradient at cell j
3342 * \param[in] grdpai upwind gradient at cell i
3343 * \param[in] grdpaj upwind gradient at cell j
3344 * \param[in] pi value at cell i
3345 * \param[in] pj value at cell j
3346 * \param[in] pia old value at cell i
3347 * \param[in] pja old value at cell j
3348 * \param[out] pifri contribution of i to flux from i to j
3349 * \param[out] pifrj contribution of i to flux from j to i
3350 * \param[out] pjfri contribution of j to flux from i to j
3351 * \param[out] pjfrj contribution of j to flux from j to i
3352 * \param[out] pip reconstructed value at cell i
3353 * \param[out] pjp reconstructed value at cell j
3354 * \param[out] pipr relaxed reconstructed value at cell i
3355 * \param[out] pjpr relaxed reconstructed value at cell j
3356 */
3357 /*----------------------------------------------------------------------------*/
3358
3359 inline static void
cs_i_cd_steady_slope_test_tensor(bool * upwind_switch,const int iconvp,const cs_real_t bldfrp,const int ischcp,const double relaxp,const double blencp,const double blend_st,const cs_real_t weight,const cs_real_t i_dist,const cs_real_t i_face_surf,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_normal,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_t i_massflux,const cs_real_63_t gradi,const cs_real_63_t gradj,const cs_real_63_t grdpai,const cs_real_63_t grdpaj,const cs_real_6_t pi,const cs_real_6_t pj,const cs_real_6_t pia,const cs_real_6_t pja,cs_real_t pifri[6],cs_real_t pifrj[6],cs_real_t pjfri[6],cs_real_t pjfrj[6],cs_real_t pip[6],cs_real_t pjp[6],cs_real_t pipr[6],cs_real_t pjpr[6])3360 cs_i_cd_steady_slope_test_tensor(bool *upwind_switch,
3361 const int iconvp,
3362 const cs_real_t bldfrp,
3363 const int ischcp,
3364 const double relaxp,
3365 const double blencp,
3366 const double blend_st,
3367 const cs_real_t weight,
3368 const cs_real_t i_dist,
3369 const cs_real_t i_face_surf,
3370 const cs_real_3_t cell_ceni,
3371 const cs_real_3_t cell_cenj,
3372 const cs_real_3_t i_face_normal,
3373 const cs_real_3_t i_face_cog,
3374 const cs_real_3_t diipf,
3375 const cs_real_3_t djjpf,
3376 const cs_real_t i_massflux,
3377 const cs_real_63_t gradi,
3378 const cs_real_63_t gradj,
3379 const cs_real_63_t grdpai,
3380 const cs_real_63_t grdpaj,
3381 const cs_real_6_t pi,
3382 const cs_real_6_t pj,
3383 const cs_real_6_t pia,
3384 const cs_real_6_t pja,
3385 cs_real_t pifri[6],
3386 cs_real_t pifrj[6],
3387 cs_real_t pjfri[6],
3388 cs_real_t pjfrj[6],
3389 cs_real_t pip[6],
3390 cs_real_t pjp[6],
3391 cs_real_t pipr[6],
3392 cs_real_t pjpr[6])
3393 {
3394 cs_real_t pir[6], pjr[6];
3395 cs_real_t recoi[6], recoj[6];
3396 cs_real_t testij, tesqck;
3397 int isou;
3398
3399 cs_i_compute_quantities_tensor(bldfrp,
3400 diipf,
3401 djjpf,
3402 gradi,
3403 gradj,
3404 pi,
3405 pj,
3406 recoi,
3407 recoj,
3408 pip,
3409 pjp);
3410
3411 cs_i_relax_c_val_tensor(relaxp,
3412 pia,
3413 pja,
3414 recoi,
3415 recoj,
3416 pi,
3417 pj,
3418 pir,
3419 pjr,
3420 pipr,
3421 pjpr);
3422
3423 /* Convection slope test is needed only when iconv >0 */
3424 if (iconvp > 0) {
3425 cs_slope_test_tensor(pi,
3426 pj,
3427 i_dist,
3428 i_face_surf,
3429 i_face_normal,
3430 gradi,
3431 gradj,
3432 grdpai,
3433 grdpaj,
3434 i_massflux,
3435 &testij,
3436 &tesqck);
3437
3438 for (isou = 0; isou < 6; isou++) {
3439 if (ischcp==1) {
3440
3441 /* Centered
3442 --------*/
3443
3444 cs_centered_f_val(weight,
3445 pip[isou],
3446 pjpr[isou],
3447 &pifrj[isou]);
3448 cs_centered_f_val(weight,
3449 pipr[isou],
3450 pjp[isou],
3451 &pifri[isou]);
3452 cs_centered_f_val(weight,
3453 pipr[isou],
3454 pjp[isou],
3455 &pjfri[isou]);
3456 cs_centered_f_val(weight,
3457 pip[isou],
3458 pjpr[isou],
3459 &pjfrj[isou]);
3460
3461 } else {
3462
3463 /* Second order
3464 ------------*/
3465
3466 cs_solu_f_val(cell_ceni,
3467 i_face_cog,
3468 gradi[isou],
3469 pi[isou],
3470 &pifrj[isou]);
3471 cs_solu_f_val(cell_ceni,
3472 i_face_cog,
3473 gradi[isou],
3474 pir[isou],
3475 &pifri[isou]);
3476 cs_solu_f_val(cell_cenj,
3477 i_face_cog,
3478 gradj[isou],
3479 pj[isou],
3480 &pjfri[isou]);
3481 cs_solu_f_val(cell_cenj,
3482 i_face_cog,
3483 gradj[isou],
3484 pjr[isou],
3485 &pjfrj[isou]);
3486
3487 }
3488
3489 }
3490
3491 /* Slope test: Pourcentage of upwind
3492 ----------------------------------*/
3493
3494 if (tesqck <= 0. || testij <= 0.) {
3495
3496 cs_blend_f_val_tensor(blend_st,
3497 pi,
3498 pifrj);
3499 cs_blend_f_val_tensor(blend_st,
3500 pir,
3501 pifri);
3502 cs_blend_f_val_tensor(blend_st,
3503 pj,
3504 pjfri);
3505 cs_blend_f_val_tensor(blend_st,
3506 pjr,
3507 pjfrj);
3508
3509 *upwind_switch = true;
3510
3511 }
3512
3513
3514 /* Blending
3515 --------*/
3516
3517 cs_blend_f_val_tensor(blencp,
3518 pi,
3519 pifrj);
3520 cs_blend_f_val_tensor(blencp,
3521 pir,
3522 pifri);
3523 cs_blend_f_val_tensor(blencp,
3524 pj,
3525 pjfri);
3526 cs_blend_f_val_tensor(blencp,
3527 pjr,
3528 pjfrj);
3529
3530 /* If iconv=0 p*fr* are useless */
3531 } else {
3532 for (isou = 0; isou < 6; isou++) {
3533 cs_upwind_f_val(pi[isou],
3534 &pifrj[isou]);
3535 cs_upwind_f_val(pir[isou],
3536 &pifri[isou]);
3537 cs_upwind_f_val(pj[isou],
3538 &pjfri[isou]);
3539 cs_upwind_f_val(pjr[isou],
3540 &pjfrj[isou]);
3541 }
3542 }
3543
3544 }
3545
3546 /*----------------------------------------------------------------------------*/
3547 /*!
3548 * \brief Handle preparation of internal face values for the fluxes computation
3549 * in case of a steady algorithm and using slope tests.
3550 *
3551 * \param[out] upwind_switch slope test result
3552 * \param[in] iconvp convection flag
3553 * \param[in] bldfrp reconstruction blending factor
3554 * \param[in] ischcp second order convection scheme flag
3555 * \param[in] blencp proportion of second order scheme,
3556 * (1-blencp) is the proportion of upwind.
3557 * \param[in] blend_st proportion of centered or SOLU scheme,
3558 * when the slope test is activated
3559 * (1-blend_st) is the proportion of upwind.
3560 * \param[in] weight geometrical weight
3561 * \param[in] i_dist distance IJ.Nij
3562 * \param[in] i_face_surf face surface
3563 * \param[in] cell_ceni center of gravity coordinates of cell i
3564 * \param[in] cell_cenj center of gravity coordinates of cell i
3565 * \param[in] i_face_normal face normal
3566 * \param[in] i_face_cog center of gravity coordinates of face ij
3567 * \param[in] diipf distance II'
3568 * \param[in] djjpf distance JJ'
3569 * \param[in] i_massflux mass flux at face ij
3570 * \param[in] gradi gradient at cell i
3571 * \param[in] gradj gradient at cell j
3572 * \param[in] gradupi upwind gradient at cell i
3573 * \param[in] gradupj upwind gradient at cell j
3574 * \param[in] gradsti slope test gradient at cell i
3575 * \param[in] gradstj slope test gradient at cell j
3576 * \param[in] pi value at cell i
3577 * \param[in] pj value at cell j
3578 * \param[out] pif contribution of i to flux from i to j
3579 * \param[out] pjf contribution of j to flux from i to j
3580 * \param[out] pip reconstructed value at cell i
3581 * \param[out] pjp reconstructed value at cell j
3582 */
3583 /*----------------------------------------------------------------------------*/
3584
3585 inline static void
cs_i_cd_unsteady_slope_test(bool * upwind_switch,const int iconvp,const cs_real_t bldfrp,const int ischcp,const double blencp,const double blend_st,const cs_real_t weight,const cs_real_t i_dist,const cs_real_t i_face_surf,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_normal,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_t i_massflux,const cs_real_3_t gradi,const cs_real_3_t gradj,const cs_real_3_t gradupi,const cs_real_3_t gradupj,const cs_real_3_t gradsti,const cs_real_3_t gradstj,const cs_real_t pi,const cs_real_t pj,cs_real_t * pif,cs_real_t * pjf,cs_real_t * pip,cs_real_t * pjp)3586 cs_i_cd_unsteady_slope_test(bool *upwind_switch,
3587 const int iconvp,
3588 const cs_real_t bldfrp,
3589 const int ischcp,
3590 const double blencp,
3591 const double blend_st,
3592 const cs_real_t weight,
3593 const cs_real_t i_dist,
3594 const cs_real_t i_face_surf,
3595 const cs_real_3_t cell_ceni,
3596 const cs_real_3_t cell_cenj,
3597 const cs_real_3_t i_face_normal,
3598 const cs_real_3_t i_face_cog,
3599 const cs_real_3_t diipf,
3600 const cs_real_3_t djjpf,
3601 const cs_real_t i_massflux,
3602 const cs_real_3_t gradi,
3603 const cs_real_3_t gradj,
3604 const cs_real_3_t gradupi,
3605 const cs_real_3_t gradupj,
3606 const cs_real_3_t gradsti,
3607 const cs_real_3_t gradstj,
3608 const cs_real_t pi,
3609 const cs_real_t pj,
3610 cs_real_t *pif,
3611 cs_real_t *pjf,
3612 cs_real_t *pip,
3613 cs_real_t *pjp)
3614 {
3615 CS_UNUSED(blend_st);
3616
3617 cs_real_t recoi, recoj;
3618 cs_real_t testij, tesqck;
3619
3620 *upwind_switch = false;
3621
3622 cs_i_compute_quantities(bldfrp,
3623 diipf,
3624 djjpf,
3625 gradi,
3626 gradj,
3627 pi,
3628 pj,
3629 &recoi,
3630 &recoj,
3631 pip,
3632 pjp);
3633
3634 /* Convection slope test is needed only when iconv >0 */
3635 if (iconvp > 0) {
3636 cs_slope_test(pi,
3637 pj,
3638 i_dist,
3639 i_face_surf,
3640 i_face_normal,
3641 gradi,
3642 gradj,
3643 gradsti,
3644 gradstj,
3645 i_massflux,
3646 &testij,
3647 &tesqck);
3648
3649 if (ischcp==1) {
3650
3651 /* Centered
3652 --------*/
3653
3654 cs_centered_f_val(weight,
3655 *pip,
3656 *pjp,
3657 pif);
3658 cs_centered_f_val(weight,
3659 *pip,
3660 *pjp,
3661 pjf);
3662
3663 } else if (ischcp == 0) {
3664
3665 /* Original SOLU
3666 --------------*/
3667
3668 cs_solu_f_val(cell_ceni,
3669 i_face_cog,
3670 gradi,
3671 pi,
3672 pif);
3673 cs_solu_f_val(cell_cenj,
3674 i_face_cog,
3675 gradj,
3676 pj,
3677 pjf);
3678
3679 } else {
3680
3681 /* SOLU
3682 -----*/
3683
3684 cs_solu_f_val(cell_ceni,
3685 i_face_cog,
3686 gradupi,
3687 pi,
3688 pif);
3689 cs_solu_f_val(cell_cenj,
3690 i_face_cog,
3691 gradupj,
3692 pj,
3693 pjf);
3694
3695 }
3696
3697 /* Slope test: Pourcentage of upwind
3698 ----------------------------------*/
3699
3700 if (tesqck<=0. || testij<=0.) {
3701
3702 cs_blend_f_val(blend_st,
3703 pi,
3704 pif);
3705 cs_blend_f_val(blend_st,
3706 pj,
3707 pjf);
3708
3709 *upwind_switch = true;
3710
3711 }
3712
3713 /* Blending
3714 --------*/
3715
3716 cs_blend_f_val(blencp,
3717 pi,
3718 pif);
3719 cs_blend_f_val(blencp,
3720 pj,
3721 pjf);
3722
3723 /* If iconv=0 p*f are useless */
3724 } else {
3725 cs_upwind_f_val(pi,
3726 pif);
3727 cs_upwind_f_val(pj,
3728 pjf);
3729 }
3730
3731 }
3732
3733 /*----------------------------------------------------------------------------*/
3734 /*!
3735 * \brief Determine the upwind and downwind sides of an internal face and
3736 * matching cell indices.
3737 *
3738 * \param[in] ii index of cell (0)
3739 * \param[in] jj index of cell (1)
3740 * \param[in] i_massflux mass flux at face ij
3741 * \param[out] ic index of central cell (upwind w.r.t. the face)
3742 * \param[out] id index of downwind cell
3743 */
3744 /*----------------------------------------------------------------------------*/
3745
3746 inline static void
cs_central_downwind_cells(const cs_lnum_t ii,const cs_lnum_t jj,const cs_real_t i_massflux,cs_lnum_t * ic,cs_lnum_t * id)3747 cs_central_downwind_cells(const cs_lnum_t ii,
3748 const cs_lnum_t jj,
3749 const cs_real_t i_massflux,
3750 cs_lnum_t *ic,
3751 cs_lnum_t *id)
3752 {
3753 if (i_massflux >= 0.) {
3754 *ic = ii;
3755 *id = jj;
3756 } else {
3757 *ic = jj;
3758 *id = ii;
3759 }
3760 }
3761
3762 /*----------------------------------------------------------------------------*/
3763 /*!
3764 * \brief Handle preparation of internal face values for the convection flux
3765 * computation in case of an unsteady algorithm and using NVD schemes.
3766 *
3767 * \param[in] limiter choice of the NVD scheme
3768 * \param[in] beta proportion of second order scheme,
3769 * (1-blencp) is the proportion of upwind.
3770 * \param[in] cell_cen_c center of gravity coordinates of central cell
3771 * \param[in] cell_cen_d center of gravity coordinates of downwind cell
3772 * \param[in] i_face_normal normal of face ij
3773 * \param[in] i_face_cog center of gravity coordinates of face ij
3774 * \param[in] i_massflux mass flux at face ij
3775 * \param[in] gradv_c gradient at central cell
3776 * \param[in] p_c value at central cell
3777 * \param[in] p_d value at downwind cell
3778 * \param[in] local_max_c local maximum of variable
3779 * \param[in] local_min_c local minimum of variable
3780 * \param[in] courant_c central cell courant number
3781 * \param[out] pif contribution of i to flux from i to j
3782 * \param[out] pjf contribution of j to flux from i to j
3783 */
3784 /*----------------------------------------------------------------------------*/
3785
3786 inline static void
cs_i_cd_unsteady_nvd(const cs_nvd_type_t limiter,const double beta,const cs_real_3_t cell_cen_c,const cs_real_3_t cell_cen_d,const cs_real_3_t i_face_normal,const cs_real_3_t i_face_cog,const cs_real_3_t gradv_c,const cs_real_t p_c,const cs_real_t p_d,const cs_real_t local_max_c,const cs_real_t local_min_c,const cs_real_t courant_c,cs_real_t * pif,cs_real_t * pjf)3787 cs_i_cd_unsteady_nvd(const cs_nvd_type_t limiter,
3788 const double beta,
3789 const cs_real_3_t cell_cen_c,
3790 const cs_real_3_t cell_cen_d,
3791 const cs_real_3_t i_face_normal,
3792 const cs_real_3_t i_face_cog,
3793 const cs_real_3_t gradv_c,
3794 const cs_real_t p_c,
3795 const cs_real_t p_d,
3796 const cs_real_t local_max_c,
3797 const cs_real_t local_min_c,
3798 const cs_real_t courant_c,
3799 cs_real_t *pif,
3800 cs_real_t *pjf)
3801 {
3802 /* distance between face center and central cell center */
3803 cs_real_t dist_fc;
3804 cs_real_3_t nfc;
3805 cs_math_3_length_unitv(cell_cen_c, i_face_cog, &dist_fc, nfc);
3806
3807 /* unit vector and distance between central and downwind cells centers */
3808 cs_real_t dist_dc;
3809 cs_real_3_t ndc;
3810 cs_math_3_length_unitv(cell_cen_c, cell_cen_d, &dist_dc, ndc);
3811
3812 /* Place the upwind point on the line that joins
3813 the two cells on the upwind side and the same
3814 distance as that between the two cells */
3815 const cs_real_t dist_cu = dist_dc;
3816 const cs_real_t dist_du = dist_dc + dist_cu;
3817
3818 /* Compute the property on the upwind assuming a parabolic
3819 variation of the property between the two cells */
3820 const cs_real_t gradc = cs_math_3_dot_product(gradv_c, ndc);
3821
3822 const cs_real_t grad2c = ((p_d - p_c)/dist_dc - gradc)/dist_dc;
3823
3824 cs_real_t p_u = p_c + (grad2c*dist_cu - gradc)*dist_cu;
3825 p_u = CS_MAX(CS_MIN(p_u, local_max_c), local_min_c);
3826
3827 /* Compute the normalised distances */
3828 const cs_real_t nvf_r_f = (dist_fc+dist_cu)/dist_du;
3829 const cs_real_t nvf_r_c = dist_cu/dist_du;
3830
3831 /* Check for the bounds of the NVD diagram and compute the face
3832 property according to the selected NVD scheme */
3833 const cs_real_t _small = cs_math_epzero
3834 * (CS_ABS(p_u)+CS_ABS(p_c)+CS_ABS(p_d));
3835
3836 if (CS_ABS(p_d-p_u) <= _small) {
3837 *pif = p_c;
3838 *pjf = p_c;
3839 } else {
3840 const cs_real_t nvf_p_c = (p_c - p_u)/(p_d - p_u);
3841
3842 if (nvf_p_c <= 0. || nvf_p_c >= 1.) {
3843 *pif = p_c;
3844 *pjf = p_c;
3845 } else {
3846 cs_real_t nvf_p_f;
3847
3848 /* Highly compressive NVD scheme for VOF */
3849 if (limiter >= CS_NVD_VOF_HRIC) {
3850 nvf_p_f = cs_nvd_vof_scheme_scalar(limiter,
3851 i_face_normal,
3852 nvf_p_c,
3853 nvf_r_f,
3854 nvf_r_c,
3855 gradv_c,
3856 courant_c);
3857 } else { /* Regular NVD scheme */
3858 nvf_p_f = cs_nvd_scheme_scalar(limiter,
3859 nvf_p_c,
3860 nvf_r_f,
3861 nvf_r_c);
3862 }
3863
3864 *pif = p_u + nvf_p_f*(p_d - p_u);
3865 *pif = CS_MAX(CS_MIN(*pif, local_max_c), local_min_c);
3866
3867 cs_blend_f_val(beta,
3868 p_c,
3869 pif);
3870
3871 *pjf = *pif;
3872 }
3873 }
3874 }
3875
3876 /*----------------------------------------------------------------------------*/
3877 /*!
3878 * \brief Handle preparation of internal face values for the fluxes computation
3879 * in case of a unsteady algorithm and using slope tests.
3880 *
3881 * \param[out] upwind_switch slope test result
3882 * \param[in] iconvp convection flag
3883 * \param[in] bldfrp reconstruction blending factor
3884 * \param[in] ischcp second order convection scheme flag
3885 * \param[in] blencp proportion of second order scheme,
3886 * (1-blencp) is the proportion of upwind.
3887 * \param[in] blend_st proportion of centered or SOLU scheme,
3888 * when the slope test is activated
3889 * (1-blend_st) is the proportion of upwind.
3890 * \param[in] weight geometrical weight
3891 * \param[in] i_dist distance IJ.Nij
3892 * \param[in] i_face_surf face surface
3893 * \param[in] cell_ceni center of gravity coordinates of cell i
3894 * \param[in] cell_cenj center of gravity coordinates of cell i
3895 * \param[in] i_face_normal face normal
3896 * \param[in] i_face_cog center of gravity coordinates of face ij
3897 * \param[in] diipf distance II'
3898 * \param[in] djjpf distance JJ'
3899 * \param[in] i_massflux mass flux at face ij
3900 * \param[in] gradi gradient at cell i
3901 * \param[in] gradj gradient at cell j
3902 * \param[in] grdpai upwind gradient at cell i
3903 * \param[in] grdpaj upwind gradient at cell j
3904 * \param[in] pi value at cell i
3905 * \param[in] pj value at cell j
3906 * \param[out] pif contribution of i to flux from i to j
3907 * \param[out] pjf contribution of j to flux from i to j
3908 * \param[out] pip reconstructed value at cell i
3909 * \param[out] pjp reconstructed value at cell j
3910 */
3911 /*----------------------------------------------------------------------------*/
3912
3913 inline static void
cs_i_cd_unsteady_slope_test_vector(bool * upwind_switch,const int iconvp,const cs_real_t bldfrp,const int ischcp,const double blencp,const double blend_st,const cs_real_t weight,const cs_real_t i_dist,const cs_real_t i_face_surf,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_normal,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_t i_massflux,const cs_real_33_t gradi,const cs_real_33_t gradj,const cs_real_33_t grdpai,const cs_real_33_t grdpaj,const cs_real_3_t pi,const cs_real_3_t pj,cs_real_t pif[3],cs_real_t pjf[3],cs_real_t pip[3],cs_real_t pjp[3])3914 cs_i_cd_unsteady_slope_test_vector(bool *upwind_switch,
3915 const int iconvp,
3916 const cs_real_t bldfrp,
3917 const int ischcp,
3918 const double blencp,
3919 const double blend_st,
3920 const cs_real_t weight,
3921 const cs_real_t i_dist,
3922 const cs_real_t i_face_surf,
3923 const cs_real_3_t cell_ceni,
3924 const cs_real_3_t cell_cenj,
3925 const cs_real_3_t i_face_normal,
3926 const cs_real_3_t i_face_cog,
3927 const cs_real_3_t diipf,
3928 const cs_real_3_t djjpf,
3929 const cs_real_t i_massflux,
3930 const cs_real_33_t gradi,
3931 const cs_real_33_t gradj,
3932 const cs_real_33_t grdpai,
3933 const cs_real_33_t grdpaj,
3934 const cs_real_3_t pi,
3935 const cs_real_3_t pj,
3936 cs_real_t pif[3],
3937 cs_real_t pjf[3],
3938 cs_real_t pip[3],
3939 cs_real_t pjp[3])
3940 {
3941 cs_real_t recoi[3], recoj[3];
3942 cs_real_t testij, tesqck;
3943
3944 cs_i_compute_quantities_vector(bldfrp,
3945 diipf,
3946 djjpf,
3947 gradi,
3948 gradj,
3949 pi,
3950 pj,
3951 recoi,
3952 recoj,
3953 pip,
3954 pjp);
3955
3956 /* Convection slope test is needed only when iconv >0 */
3957 if (iconvp > 0) {
3958 cs_slope_test_vector(pi,
3959 pj,
3960 i_dist,
3961 i_face_surf,
3962 i_face_normal,
3963 gradi,
3964 gradj,
3965 grdpai,
3966 grdpaj,
3967 i_massflux,
3968 &testij,
3969 &tesqck);
3970
3971 for (int isou = 0; isou < 3; isou++) {
3972 if (ischcp == 1) {
3973
3974 /* Centered
3975 --------*/
3976
3977 cs_centered_f_val(weight,
3978 pip[isou],
3979 pjp[isou],
3980 &pif[isou]);
3981 cs_centered_f_val(weight,
3982 pip[isou],
3983 pjp[isou],
3984 &pjf[isou]);
3985
3986 } else {
3987
3988 /* Second order
3989 ------------*/
3990
3991 cs_solu_f_val(cell_ceni,
3992 i_face_cog,
3993 gradi[isou],
3994 pi[isou],
3995 &pif[isou]);
3996 cs_solu_f_val(cell_cenj,
3997 i_face_cog,
3998 gradj[isou],
3999 pj[isou],
4000 &pjf[isou]);
4001
4002 }
4003
4004 }
4005
4006 /* Slope test: Pourcentage of upwind
4007 ----------------------------------*/
4008
4009 if (tesqck <= 0. || testij <= 0.) {
4010
4011 cs_blend_f_val_vector(blend_st,
4012 pi,
4013 pif);
4014 cs_blend_f_val_vector(blend_st,
4015 pj,
4016 pjf);
4017
4018 *upwind_switch = true;
4019
4020 }
4021
4022
4023 /* Blending
4024 --------*/
4025 cs_blend_f_val_vector(blencp,
4026 pi,
4027 pif);
4028 cs_blend_f_val_vector(blencp,
4029 pj,
4030 pjf);
4031
4032 /* If iconv=0 p*f are useless */
4033 } else {
4034
4035 for (int isou = 0; isou < 3; isou++) {
4036 cs_upwind_f_val(pi[isou],
4037 &pif[isou]);
4038 cs_upwind_f_val(pj[isou],
4039 &pjf[isou]);
4040
4041 }
4042 }
4043
4044 }
4045
4046 /*----------------------------------------------------------------------------*/
4047 /*!
4048 * \brief Handle preparation of internal face values for the fluxes computation
4049 * in case of a unsteady algorithm and using slope tests.
4050 *
4051 * \param[out] upwind_switch slope test result
4052 * \param[in] iconvp convection flag
4053 * \param[in] bldfrp reconstruction blending factor
4054 * \param[in] ischcp second order convection scheme flag
4055 * \param[in] blencp proportion of second order scheme,
4056 * (1-blencp) is the proportion of upwind.
4057 * \param[in] blend_st proportion of centered or SOLU scheme,
4058 * when the slope test is activated
4059 * (1-blend_st) is the proportion of upwind.
4060 * \param[in] weight geometrical weight
4061 * \param[in] i_dist distance IJ.Nij
4062 * \param[in] i_face_surf face surface
4063 * \param[in] cell_ceni center of gravity coordinates of cell i
4064 * \param[in] cell_cenj center of gravity coordinates of cell i
4065 * \param[in] i_face_normal face normal
4066 * \param[in] i_face_cog center of gravity coordinates of face ij
4067 * \param[in] diipf distance II'
4068 * \param[in] djjpf distance JJ'
4069 * \param[in] i_massflux mass flux at face ij
4070 * \param[in] gradi gradient at cell i
4071 * \param[in] gradj gradient at cell j
4072 * \param[in] grdpai upwind gradient at cell i
4073 * \param[in] grdpaj upwind gradient at cell j
4074 * \param[in] pi value at cell i
4075 * \param[in] pj value at cell j
4076 * \param[out] pif contribution of i to flux from i to j
4077 * \param[out] pjf contribution of j to flux from i to j
4078 * \param[out] pip reconstructed value at cell i
4079 * \param[out] pjp reconstructed value at cell j
4080 */
4081 /*----------------------------------------------------------------------------*/
4082
4083 inline static void
cs_i_cd_unsteady_slope_test_tensor(bool * upwind_switch,const int iconvp,const cs_real_t bldfrp,const int ischcp,const double blencp,const double blend_st,const cs_real_t weight,const cs_real_t i_dist,const cs_real_t i_face_surf,const cs_real_3_t cell_ceni,const cs_real_3_t cell_cenj,const cs_real_3_t i_face_normal,const cs_real_3_t i_face_cog,const cs_real_3_t diipf,const cs_real_3_t djjpf,const cs_real_t i_massflux,const cs_real_63_t gradi,const cs_real_63_t gradj,const cs_real_63_t grdpai,const cs_real_63_t grdpaj,const cs_real_6_t pi,const cs_real_6_t pj,cs_real_t pif[6],cs_real_t pjf[6],cs_real_t pip[6],cs_real_t pjp[6])4084 cs_i_cd_unsteady_slope_test_tensor(bool *upwind_switch,
4085 const int iconvp,
4086 const cs_real_t bldfrp,
4087 const int ischcp,
4088 const double blencp,
4089 const double blend_st,
4090 const cs_real_t weight,
4091 const cs_real_t i_dist,
4092 const cs_real_t i_face_surf,
4093 const cs_real_3_t cell_ceni,
4094 const cs_real_3_t cell_cenj,
4095 const cs_real_3_t i_face_normal,
4096 const cs_real_3_t i_face_cog,
4097 const cs_real_3_t diipf,
4098 const cs_real_3_t djjpf,
4099 const cs_real_t i_massflux,
4100 const cs_real_63_t gradi,
4101 const cs_real_63_t gradj,
4102 const cs_real_63_t grdpai,
4103 const cs_real_63_t grdpaj,
4104 const cs_real_6_t pi,
4105 const cs_real_6_t pj,
4106 cs_real_t pif[6],
4107 cs_real_t pjf[6],
4108 cs_real_t pip[6],
4109 cs_real_t pjp[6])
4110 {
4111 cs_real_t recoi[6], recoj[6];
4112 cs_real_t testij, tesqck;
4113 int isou;
4114
4115 cs_i_compute_quantities_tensor(bldfrp,
4116 diipf,
4117 djjpf,
4118 gradi,
4119 gradj,
4120 pi,
4121 pj,
4122 recoi,
4123 recoj,
4124 pip,
4125 pjp);
4126
4127 /* Convection slope test is needed only when iconv >0 */
4128 if (iconvp > 0) {
4129 cs_slope_test_tensor(pi,
4130 pj,
4131 i_dist,
4132 i_face_surf,
4133 i_face_normal,
4134 gradi,
4135 gradj,
4136 grdpai,
4137 grdpaj,
4138 i_massflux,
4139 &testij,
4140 &tesqck);
4141
4142 for (isou = 0; isou < 6; isou++) {
4143
4144 if (ischcp==1) {
4145
4146 /* Centered
4147 --------*/
4148
4149 cs_centered_f_val(weight,
4150 pip[isou],
4151 pjp[isou],
4152 &pif[isou]);
4153 cs_centered_f_val(weight,
4154 pip[isou],
4155 pjp[isou],
4156 &pjf[isou]);
4157
4158 } else {
4159
4160 /* Second order
4161 ------------*/
4162
4163 cs_solu_f_val(cell_ceni,
4164 i_face_cog,
4165 gradi[isou],
4166 pi[isou],
4167 &pif[isou]);
4168 cs_solu_f_val(cell_cenj,
4169 i_face_cog,
4170 gradj[isou],
4171 pj[isou],
4172 &pjf[isou]);
4173 }
4174
4175 }
4176
4177 /* Slope test activated: poucentage of upwind */
4178 if (tesqck <= 0. || testij <= 0.) {
4179
4180 /* Upwind
4181 --------*/
4182
4183 cs_blend_f_val_tensor(blend_st,
4184 pi,
4185 pif);
4186 cs_blend_f_val_tensor(blend_st,
4187 pj,
4188 pjf);
4189
4190 *upwind_switch = true;
4191 }
4192
4193
4194 /* Blending
4195 --------*/
4196
4197 cs_blend_f_val_tensor(blencp,
4198 pi,
4199 pif);
4200 cs_blend_f_val_tensor(blencp,
4201 pj,
4202 pjf);
4203
4204 /* If iconv=0 p*fr* are useless */
4205 } else {
4206
4207 for (isou = 0; isou < 6; isou++) {
4208 cs_upwind_f_val(pi[isou],
4209 &pif[isou]);
4210 cs_upwind_f_val(pj[isou],
4211 &pjf[isou]);
4212 }
4213 }
4214 }
4215
4216 /*----------------------------------------------------------------------------*/
4217 /*!
4218 * \brief Reconstruct values in I' at boundary cell i.
4219 *
4220 * \param[in] diipb distance I'I'
4221 * \param[in] gradi gradient at cell i
4222 * \param[in] bldfrp reconstruction blending factor
4223 * \param[out] recoi reconstruction at cell i
4224 */
4225 /*----------------------------------------------------------------------------*/
4226
4227 inline static void
cs_b_compute_quantities(const cs_real_3_t diipb,const cs_real_3_t gradi,const cs_real_t bldfrp,cs_real_t * recoi)4228 cs_b_compute_quantities(const cs_real_3_t diipb,
4229 const cs_real_3_t gradi,
4230 const cs_real_t bldfrp,
4231 cs_real_t *recoi)
4232 {
4233 *recoi = bldfrp * ( gradi[0]*diipb[0]
4234 + gradi[1]*diipb[1]
4235 + gradi[2]*diipb[2]);
4236 }
4237
4238 /*----------------------------------------------------------------------------*/
4239 /*!
4240 * \brief Reconstruct values in I' at boundary cell i.
4241 *
4242 * \param[in] diipb distance I'I'
4243 * \param[in] gradi gradient at cell i
4244 * \param[in] bldfrp reconstruction blending factor
4245 * \param[out] recoi reconstruction at cell i
4246 */
4247 /*----------------------------------------------------------------------------*/
4248
4249 inline static void
cs_b_compute_quantities_vector(const cs_real_3_t diipb,const cs_real_33_t gradi,const cs_real_t bldfrp,cs_real_t recoi[3])4250 cs_b_compute_quantities_vector(const cs_real_3_t diipb,
4251 const cs_real_33_t gradi,
4252 const cs_real_t bldfrp,
4253 cs_real_t recoi[3])
4254 {
4255 for (int isou = 0; isou < 3; isou++) {
4256 recoi[isou] = bldfrp * ( gradi[isou][0]*diipb[0]
4257 + gradi[isou][1]*diipb[1]
4258 + gradi[isou][2]*diipb[2]);
4259 }
4260 }
4261
4262 /*----------------------------------------------------------------------------*/
4263 /*!
4264 * \brief Reconstruct values in I' at boundary cell i.
4265 *
4266 * \param[in] diipb distance I'I'
4267 * \param[in] gradi gradient at cell i
4268 * \param[in] bldfrp reconstruction blending factor
4269 * \param[out] recoi reconstruction at cell i
4270 */
4271 /*----------------------------------------------------------------------------*/
4272
4273 inline static void
cs_b_compute_quantities_tensor(const cs_real_3_t diipb,const cs_real_63_t gradi,const cs_real_t bldfrp,cs_real_t recoi[6])4274 cs_b_compute_quantities_tensor(const cs_real_3_t diipb,
4275 const cs_real_63_t gradi,
4276 const cs_real_t bldfrp,
4277 cs_real_t recoi[6])
4278 {
4279 for (int isou = 0; isou < 6; isou++) {
4280 recoi[isou] = bldfrp * ( gradi[isou][0]*diipb[0]
4281 + gradi[isou][1]*diipb[1]
4282 + gradi[isou][2]*diipb[2]);
4283 }
4284 }
4285
4286 /*----------------------------------------------------------------------------*/
4287 /*!
4288 * \brief Compute relaxed values at boundary cell i.
4289 *
4290 * \param[in] relaxp relaxation coefficient
4291 * \param[in] pi value at cell i
4292 * \param[in] pia old value at cell i
4293 * \param[in] recoi reconstruction at cell i
4294 * \param[out] pir relaxed value at cell i
4295 * \param[out] pipr relaxed reconstructed value at cell i
4296 */
4297 /*----------------------------------------------------------------------------*/
4298
4299 inline static void
cs_b_relax_c_val(const double relaxp,const cs_real_t pi,const cs_real_t pia,const cs_real_t recoi,cs_real_t * pir,cs_real_t * pipr)4300 cs_b_relax_c_val(const double relaxp,
4301 const cs_real_t pi,
4302 const cs_real_t pia,
4303 const cs_real_t recoi,
4304 cs_real_t *pir,
4305 cs_real_t *pipr)
4306 {
4307 *pir = pi/relaxp - (1.-relaxp)/relaxp*pia;
4308 *pipr = *pir + recoi;
4309 }
4310
4311 /*----------------------------------------------------------------------------*/
4312 /*!
4313 * \brief Compute relaxed values at boundary cell i.
4314 *
4315 * \param[in] relaxp relaxation coefficient
4316 * \param[in] pi value at cell i
4317 * \param[in] pia old value at cell i
4318 * \param[in] recoi reconstruction at cell i
4319 * \param[out] pir relaxed value at cell i
4320 * \param[out] pipr relaxed reconstructed value at cell i
4321 */
4322 /*----------------------------------------------------------------------------*/
4323
4324 inline static void
cs_b_relax_c_val_vector(const double relaxp,const cs_real_3_t pi,const cs_real_3_t pia,const cs_real_3_t recoi,cs_real_t pir[3],cs_real_t pipr[3])4325 cs_b_relax_c_val_vector(const double relaxp,
4326 const cs_real_3_t pi,
4327 const cs_real_3_t pia,
4328 const cs_real_3_t recoi,
4329 cs_real_t pir[3],
4330 cs_real_t pipr[3])
4331 {
4332 for (int isou = 0; isou < 3; isou++) {
4333 pir[isou] = pi[isou]/relaxp - (1.-relaxp)/relaxp*pia[isou];
4334 pipr[isou] = pir[isou] + recoi[isou];
4335 }
4336 }
4337
4338 /*----------------------------------------------------------------------------*/
4339 /*!
4340 * \brief Compute relaxed values at boundary cell i.
4341 *
4342 * \param[in] relaxp relaxation coefficient
4343 * \param[in] pi value at cell i
4344 * \param[in] pia old value at cell i
4345 * \param[in] recoi reconstruction at cell i
4346 * \param[out] pir relaxed value at cell i
4347 * \param[out] pipr relaxed reconstructed value at cell i
4348 */
4349 /*----------------------------------------------------------------------------*/
4350
4351 inline static void
cs_b_relax_c_val_tensor(const double relaxp,const cs_real_6_t pi,const cs_real_6_t pia,const cs_real_6_t recoi,cs_real_t pir[6],cs_real_t pipr[6])4352 cs_b_relax_c_val_tensor(const double relaxp,
4353 const cs_real_6_t pi,
4354 const cs_real_6_t pia,
4355 const cs_real_6_t recoi,
4356 cs_real_t pir[6],
4357 cs_real_t pipr[6])
4358 {
4359 for (int isou = 0; isou < 6; isou++) {
4360 pir[isou] = pi[isou]/relaxp - (1.-relaxp)/relaxp*pia[isou];
4361 pipr[isou] = pir[isou] + recoi[isou];
4362 }
4363 }
4364
4365 /*----------------------------------------------------------------------------*/
4366 /*!
4367 * \brief Add convective flux (substracting the mass accumulation from it) to
4368 * flux at boundary face. The convective flux can be either an upwind flux or an
4369 * imposed value.
4370 *
4371 * \param[in] iconvp convection flag
4372 * \param[in] thetap weighting coefficient for the theta-schema,
4373 * \param[in] imasac take mass accumulation into account?
4374 * \param[in] inc Not an increment flag
4375 * \param[in] bc_type type of boundary face
4376 * \param[in] icvfli imposed convective flux flag
4377 * \param[in] pi value at cell i
4378 * \param[in] pir relaxed value at cell i
4379 * \param[in] pipr relaxed reconstructed value at cell i
4380 * \param[in] coefap explicit boundary coefficient for convection operator
4381 * \param[in] coefbp implicit boundary coefficient for convection operator
4382 * \param[in] coface explicit imposed convective flux value (0 otherwise).
4383 * \param[in] cofbce implicit part of imp. conv. flux value
4384 * \param[in] b_massflux mass flux at boundary face
4385 * \param[in] xcpp specific heat value if the scalar is the temperature,
4386 * 1 otherwise
4387 * \param[in,out] flux flux at boundary face
4388 */
4389 /*----------------------------------------------------------------------------*/
4390
4391 inline static void
cs_b_imposed_conv_flux(int iconvp,cs_real_t thetap,int imasac,int inc,int bc_type,int icvfli,cs_real_t pi,cs_real_t pir,cs_real_t pipr,cs_real_t coefap,cs_real_t coefbp,cs_real_t coface,cs_real_t cofbce,cs_real_t b_massflux,cs_real_t xcpp,cs_real_t * flux)4392 cs_b_imposed_conv_flux(int iconvp,
4393 cs_real_t thetap,
4394 int imasac,
4395 int inc,
4396 int bc_type,
4397 int icvfli,
4398 cs_real_t pi,
4399 cs_real_t pir,
4400 cs_real_t pipr,
4401 cs_real_t coefap,
4402 cs_real_t coefbp,
4403 cs_real_t coface,
4404 cs_real_t cofbce,
4405 cs_real_t b_massflux,
4406 cs_real_t xcpp,
4407 cs_real_t *flux)
4408 {
4409 cs_real_t flui, fluj, pfac;
4410
4411 /* Computed convective flux */
4412
4413 if (icvfli == 0) {
4414
4415 /* Remove decentering for coupled faces */
4416 if (bc_type == CS_COUPLED_FD) {
4417 flui = 0.0;
4418 fluj = b_massflux;
4419 } else {
4420 flui = 0.5*(b_massflux +fabs(b_massflux));
4421 fluj = 0.5*(b_massflux -fabs(b_massflux));
4422 }
4423
4424 pfac = inc*coefap + coefbp*pipr;
4425 *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*pfac) -imasac*( b_massflux*pi));
4426
4427 /* Imposed convective flux */
4428
4429 } else {
4430
4431 pfac = inc*coface + cofbce*pipr;
4432 *flux += iconvp*xcpp*(-imasac*(b_massflux*pi) + thetap*(pfac));
4433
4434 }
4435 }
4436
4437 /*----------------------------------------------------------------------------*/
4438 /*!
4439 * \brief Add convective flux (substracting the mass accumulation from it) to
4440 * flux at boundary face. The convective flux can be either an upwind flux or an
4441 * imposed value.
4442 *
4443 * \param[in] iconvp convection flag
4444 * \param[in] thetap weighting coefficient for the theta-schema,
4445 * \param[in] imasac take mass accumulation into account?
4446 * \param[in] inc Not an increment flag
4447 * \param[in] bc_type type of boundary face
4448 * \param[in] icvfli imposed convective flux flag
4449 * \param[in] pi value at cell i
4450 * \param[in] pir relaxed value at cell i
4451 * \param[in] pipr relaxed reconstructed value at cell i
4452 * \param[in] coefap explicit boundary coefficient for convection operator
4453 * \param[in] coefbp implicit boundary coefficient for convection operator
4454 * \param[in] coface explicit imposed convective flux value (0 otherwise).
4455 * \param[in] cofbce implicit part of imp. conv. flux value
4456 * \param[in] b_massflux mass flux at boundary face
4457 * \param[in,out] flux flux at boundary face
4458 */
4459 /*----------------------------------------------------------------------------*/
4460
4461 inline static void
cs_b_imposed_conv_flux_vector(int iconvp,cs_real_t thetap,int imasac,int inc,int bc_type,int icvfli,const cs_real_t pi[restrict3],const cs_real_t pir[restrict3],const cs_real_t pipr[restrict3],const cs_real_t coefap[restrict3],const cs_real_t coefbp[restrict3][3],const cs_real_t coface[restrict3],const cs_real_t cofbce[restrict3][3],cs_real_t b_massflux,cs_real_t flux[restrict3])4462 cs_b_imposed_conv_flux_vector(int iconvp,
4463 cs_real_t thetap,
4464 int imasac,
4465 int inc,
4466 int bc_type,
4467 int icvfli,
4468 const cs_real_t pi[restrict 3],
4469 const cs_real_t pir[restrict 3],
4470 const cs_real_t pipr[restrict 3],
4471 const cs_real_t coefap[restrict 3],
4472 const cs_real_t coefbp[restrict 3][3],
4473 const cs_real_t coface[restrict 3],
4474 const cs_real_t cofbce[restrict 3][3],
4475 cs_real_t b_massflux,
4476 cs_real_t flux[restrict 3])
4477 {
4478 cs_real_t flui, fluj, pfac;
4479
4480 /* Computed convective flux */
4481
4482 if (icvfli == 0) {
4483
4484 /* Remove decentering for coupled faces */
4485 if (bc_type == CS_COUPLED_FD) {
4486 flui = 0.0;
4487 fluj = b_massflux;
4488 } else {
4489 flui = 0.5*(b_massflux +fabs(b_massflux));
4490 fluj = 0.5*(b_massflux -fabs(b_massflux));
4491 }
4492 for (int isou = 0; isou < 3; isou++) {
4493 pfac = inc*coefap[isou];
4494 for (int jsou = 0; jsou < 3; jsou++) {
4495 pfac += coefbp[isou][jsou]*pipr[jsou];
4496 }
4497 flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac)
4498 - imasac*b_massflux*pi[isou]);
4499 }
4500
4501 /* Imposed convective flux */
4502
4503 } else {
4504
4505 for (int isou = 0; isou < 3; isou++) {
4506 pfac = inc*coface[isou];
4507 for (int jsou = 0; jsou < 3; jsou++) {
4508 pfac += cofbce[isou][jsou]*pipr[jsou];
4509 }
4510 flux[isou] += iconvp*( thetap*pfac
4511 - imasac*b_massflux*pi[isou]);
4512 }
4513
4514 }
4515 }
4516
4517 /*----------------------------------------------------------------------------*/
4518 /*!
4519 * \brief Add convective flux (substracting the mass accumulation from it) to
4520 * flux at boundary face. The convective flux is a pure upwind flux.
4521 *
4522 * \param[in] iconvp convection flag
4523 * \param[in] thetap weighting coefficient for the theta-schema,
4524 * \param[in] imasac take mass accumulation into account?
4525 * \param[in] inc Not an increment flag
4526 * \param[in] bc_type type of boundary face
4527 * \param[in] pi value at cell i
4528 * \param[in] pir relaxed value at cell i
4529 * \param[in] pipr relaxed reconstructed value at cell i
4530 * \param[in] coefap explicit boundary coefficient for convection operator
4531 * \param[in] coefbp implicit boundary coefficient for convection operator
4532 * \param[in] b_massflux mass flux at boundary face
4533 * \param[in] xcpp specific heat value if the scalar is the
4534 * temperature, 1 otherwise
4535 * \param[in,out] flux flux at boundary face
4536 */
4537 /*----------------------------------------------------------------------------*/
4538
4539 inline static void
cs_b_upwind_flux(const int iconvp,const cs_real_t thetap,const int imasac,const int inc,const int bc_type,const cs_real_t pi,const cs_real_t pir,const cs_real_t pipr,const cs_real_t coefap,const cs_real_t coefbp,const cs_real_t b_massflux,const cs_real_t xcpp,cs_real_t * flux)4540 cs_b_upwind_flux(const int iconvp,
4541 const cs_real_t thetap,
4542 const int imasac,
4543 const int inc,
4544 const int bc_type,
4545 const cs_real_t pi,
4546 const cs_real_t pir,
4547 const cs_real_t pipr,
4548 const cs_real_t coefap,
4549 const cs_real_t coefbp,
4550 const cs_real_t b_massflux,
4551 const cs_real_t xcpp,
4552 cs_real_t *flux)
4553 {
4554 cs_real_t flui, fluj, pfac;
4555
4556 /* Remove decentering for coupled faces */
4557 if (bc_type == CS_COUPLED_FD) {
4558 flui = 0.0;
4559 fluj = b_massflux;
4560 } else {
4561 flui = 0.5*(b_massflux +fabs(b_massflux));
4562 fluj = 0.5*(b_massflux -fabs(b_massflux));
4563 }
4564
4565 pfac = inc*coefap + coefbp*pipr;
4566 *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*pfac) -imasac*( b_massflux*pi));
4567 }
4568
4569 /*----------------------------------------------------------------------------*/
4570 /*!
4571 * \brief Add convective flux (substracting the mass accumulation from it) to
4572 * flux at boundary face. The convective flux is a pure upwind flux.
4573 *
4574 * \param[in] iconvp convection flag
4575 * \param[in] thetap weighting coefficient for the theta-schema,
4576 * \param[in] imasac take mass accumulation into account?
4577 * \param[in] inc Not an increment flag
4578 * \param[in] bc_type type of boundary face
4579 * \param[in] pi value at cell i
4580 * \param[in] pir relaxed value at cell i
4581 * \param[in] pipr relaxed reconstructed value at cell i
4582 * \param[in] coefa explicit boundary coefficient for convection
4583 * operator
4584 * \param[in] coefb implicit boundary coefficient for convection
4585 * operator
4586 * \param[in] b_massflux mass flux at boundary face
4587 * \param[in,out] flux flux at boundary face
4588 */
4589 /*----------------------------------------------------------------------------*/
4590
4591 inline static void
cs_b_upwind_flux_vector(const int iconvp,const cs_real_t thetap,const int imasac,const int inc,const int bc_type,const cs_real_3_t pi,const cs_real_3_t pir,const cs_real_3_t pipr,const cs_real_3_t coefa,const cs_real_33_t coefb,const cs_real_t b_massflux,cs_real_t flux[3])4592 cs_b_upwind_flux_vector(const int iconvp,
4593 const cs_real_t thetap,
4594 const int imasac,
4595 const int inc,
4596 const int bc_type,
4597 const cs_real_3_t pi,
4598 const cs_real_3_t pir,
4599 const cs_real_3_t pipr,
4600 const cs_real_3_t coefa,
4601 const cs_real_33_t coefb,
4602 const cs_real_t b_massflux,
4603 cs_real_t flux[3])
4604 {
4605 cs_real_t flui, fluj, pfac;
4606
4607 /* Remove decentering for coupled faces */
4608 if (bc_type == CS_COUPLED_FD) {
4609 flui = 0.0;
4610 fluj = b_massflux;
4611 } else {
4612 flui = 0.5*(b_massflux +fabs(b_massflux));
4613 fluj = 0.5*(b_massflux -fabs(b_massflux));
4614 }
4615 for (int isou = 0; isou < 3; isou++) {
4616 pfac = inc*coefa[isou];
4617 for (int jsou = 0; jsou < 3; jsou++) {
4618 pfac += coefb[isou][jsou]*pipr[jsou];
4619 }
4620 flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac)
4621 - imasac*b_massflux*pi[isou]);
4622 }
4623 }
4624
4625 /*----------------------------------------------------------------------------*/
4626 /*!
4627 * \brief Add convective flux (substracting the mass accumulation from it) to
4628 * flux at boundary face. The convective flux is a pure upwind flux.
4629 *
4630 * \param[in] iconvp convection flag
4631 * \param[in] thetap weighting coefficient for the theta-schema,
4632 * \param[in] imasac take mass accumulation into account?
4633 * \param[in] inc Not an increment flag
4634 * \param[in] bc_type type of boundary face
4635 * \param[in] pi value at cell i
4636 * \param[in] pir relaxed value at cell i
4637 * \param[in] pipr relaxed reconstructed value at cell i
4638 * \param[in] coefa explicit boundary coefficient for convection
4639 * operator
4640 * \param[in] coefb implicit boundary coefficient for convection
4641 * operator
4642 * \param[in] b_massflux mass flux at boundary face
4643 * \param[in,out] flux flux at boundary face
4644 */
4645 /*----------------------------------------------------------------------------*/
4646
4647 inline static void
cs_b_upwind_flux_tensor(const int iconvp,const cs_real_t thetap,const int imasac,const int inc,const int bc_type,const cs_real_6_t pi,const cs_real_6_t pir,const cs_real_6_t pipr,const cs_real_6_t coefa,const cs_real_66_t coefb,const cs_real_t b_massflux,cs_real_t flux[6])4648 cs_b_upwind_flux_tensor(const int iconvp,
4649 const cs_real_t thetap,
4650 const int imasac,
4651 const int inc,
4652 const int bc_type,
4653 const cs_real_6_t pi,
4654 const cs_real_6_t pir,
4655 const cs_real_6_t pipr,
4656 const cs_real_6_t coefa,
4657 const cs_real_66_t coefb,
4658 const cs_real_t b_massflux,
4659 cs_real_t flux[6])
4660 {
4661 cs_real_t flui, fluj, pfac;
4662
4663 /* Remove decentering for coupled faces */
4664 if (bc_type == CS_COUPLED_FD) {
4665 flui = 0.0;
4666 fluj = b_massflux;
4667 } else {
4668 flui = 0.5*(b_massflux +fabs(b_massflux));
4669 fluj = 0.5*(b_massflux -fabs(b_massflux));
4670 }
4671 for (int isou = 0; isou < 6; isou++) {
4672 pfac = inc*coefa[isou];
4673 for (int jsou = 0; jsou < 6; jsou++) {
4674 pfac += coefb[isou][jsou]*pipr[jsou];
4675 }
4676 flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac)
4677 - imasac*b_massflux*pi[isou]);
4678 }
4679 }
4680
4681 /*----------------------------------------------------------------------------*/
4682 /*!
4683 * \brief Add diffusive flux to flux at boundary face.
4684 *
4685 * \param[in] idiffp diffusion flag
4686 * \param[in] thetap weighting coefficient for the theta-schema,
4687 * \param[in] inc Not an increment flag
4688 * \param[in] pipr relaxed reconstructed value at cell i
4689 * \param[in] cofafp explicit boundary coefficient for diffusion operator
4690 * \param[in] cofbfp implicit boundary coefficient for diffusion operator
4691 * \param[in] b_visc boundary face surface
4692 * \param[in,out] flux flux at boundary face
4693 */
4694 /*----------------------------------------------------------------------------*/
4695
4696 inline static void
cs_b_diff_flux(const int idiffp,const cs_real_t thetap,const int inc,const cs_real_t pipr,const cs_real_t cofafp,const cs_real_t cofbfp,const cs_real_t b_visc,cs_real_t * flux)4697 cs_b_diff_flux(const int idiffp,
4698 const cs_real_t thetap,
4699 const int inc,
4700 const cs_real_t pipr,
4701 const cs_real_t cofafp,
4702 const cs_real_t cofbfp,
4703 const cs_real_t b_visc,
4704 cs_real_t *flux)
4705 {
4706 cs_real_t pfacd = inc*cofafp + cofbfp*pipr;
4707 *flux += idiffp*thetap*b_visc*pfacd;
4708 }
4709
4710 /*----------------------------------------------------------------------------*/
4711 /*!
4712 * \brief Add diffusive flux to flux at boundary face.
4713 *
4714 * \param[in] idiffp diffusion flag
4715 * \param[in] thetap weighting coefficient for the theta-schema,
4716 * \param[in] inc Not an increment flag
4717 * \param[in] pipr relaxed reconstructed value at cell i
4718 * \param[in] cofaf explicit boundary coefficient for diffusion operator
4719 * \param[in] cofbf implicit boundary coefficient for diffusion operator
4720 * \param[in] b_visc boundary face surface
4721 * \param[in,out] flux flux at boundary face
4722 */
4723 /*----------------------------------------------------------------------------*/
4724
4725 inline static void
cs_b_diff_flux_vector(const int idiffp,const cs_real_t thetap,const int inc,const cs_real_3_t pipr,const cs_real_3_t cofaf,const cs_real_33_t cofbf,const cs_real_t b_visc,cs_real_t flux[3])4726 cs_b_diff_flux_vector(const int idiffp,
4727 const cs_real_t thetap,
4728 const int inc,
4729 const cs_real_3_t pipr,
4730 const cs_real_3_t cofaf,
4731 const cs_real_33_t cofbf,
4732 const cs_real_t b_visc,
4733 cs_real_t flux[3])
4734 {
4735 cs_real_t pfacd ;
4736 for (int isou = 0; isou < 3; isou++) {
4737 pfacd = inc*cofaf[isou];
4738 for (int jsou = 0; jsou < 3; jsou++) {
4739 pfacd += cofbf[isou][jsou]*pipr[jsou];
4740 }
4741 flux[isou] += idiffp*thetap*b_visc*pfacd;
4742 }
4743 }
4744
4745 /*----------------------------------------------------------------------------*/
4746 /*!
4747 * \brief Add diffusive flux to flux at boundary face.
4748 *
4749 * \param[in] idiffp diffusion flag
4750 * \param[in] thetap weighting coefficient for the theta-schema,
4751 * \param[in] inc Not an increment flag
4752 * \param[in] pipr relaxed reconstructed value at cell i
4753 * \param[in] cofaf explicit boundary coefficient for diffusion operator
4754 * \param[in] cofbf implicit boundary coefficient for diffusion operator
4755 * \param[in] b_visc boundary face surface
4756 * \param[in,out] flux flux at boundary face
4757 */
4758 /*----------------------------------------------------------------------------*/
4759
4760 inline static void
cs_b_diff_flux_tensor(const int idiffp,const cs_real_t thetap,const int inc,const cs_real_6_t pipr,const cs_real_6_t cofaf,const cs_real_66_t cofbf,const cs_real_t b_visc,cs_real_t flux[6])4761 cs_b_diff_flux_tensor(const int idiffp,
4762 const cs_real_t thetap,
4763 const int inc,
4764 const cs_real_6_t pipr,
4765 const cs_real_6_t cofaf,
4766 const cs_real_66_t cofbf,
4767 const cs_real_t b_visc,
4768 cs_real_t flux[6])
4769 {
4770 cs_real_t pfacd ;
4771 for (int isou = 0; isou < 6; isou++) {
4772 pfacd = inc*cofaf[isou];
4773 for (int jsou = 0; jsou < 6; jsou++) {
4774 pfacd += cofbf[isou][jsou]*pipr[jsou];
4775 }
4776 flux[isou] += idiffp*thetap*b_visc*pfacd;
4777 }
4778 }
4779
4780 /*----------------------------------------------------------------------------*/
4781 /*!
4782 * \brief Handle preparation of boundary face values for the flux computation in
4783 * case of a steady algorithm.
4784 *
4785 * \param[in] bldfrp reconstruction blending factor
4786 * \param[in] relaxp relaxation coefficient
4787 * \param[in] diipb distance I'I'
4788 * \param[in] gradi gradient at cell i
4789 * \param[in] pi value at cell i
4790 * \param[in] pia old value at cell i
4791 * \param[out] pir relaxed value at cell i
4792 * \param[out] pipr relaxed reconstructed value at cell i
4793 */
4794 /*----------------------------------------------------------------------------*/
4795
4796 inline static void
cs_b_cd_steady(const cs_real_t bldfrp,const double relaxp,const cs_real_3_t diipb,const cs_real_3_t gradi,const cs_real_t pi,const cs_real_t pia,cs_real_t * pir,cs_real_t * pipr)4797 cs_b_cd_steady(const cs_real_t bldfrp,
4798 const double relaxp,
4799 const cs_real_3_t diipb,
4800 const cs_real_3_t gradi,
4801 const cs_real_t pi,
4802 const cs_real_t pia,
4803 cs_real_t *pir,
4804 cs_real_t *pipr)
4805 {
4806 cs_real_t recoi;
4807
4808 cs_b_compute_quantities(diipb,
4809 gradi,
4810 bldfrp,
4811 &recoi);
4812
4813 cs_b_relax_c_val(relaxp,
4814 pi,
4815 pia,
4816 recoi,
4817 pir,
4818 pipr);
4819 }
4820
4821 /*----------------------------------------------------------------------------*/
4822 /*!
4823 * \brief Handle preparation of boundary face values for the flux computation in
4824 * case of a steady algorithm.
4825 *
4826 * \param[in] bldfrp reconstruction blending factor
4827 * \param[in] relaxp relaxation coefficient
4828 * \param[in] diipb distance I'I'
4829 * \param[in] gradi gradient at cell i
4830 * \param[in] pi value at cell i
4831 * \param[in] pia old value at cell i
4832 * \param[out] pir relaxed value at cell i
4833 * \param[out] pipr relaxed reconstructed value at cell i
4834 */
4835 /*----------------------------------------------------------------------------*/
4836
4837 inline static void
cs_b_cd_steady_vector(const cs_real_t bldfrp,const double relaxp,const cs_real_3_t diipb,const cs_real_33_t gradi,const cs_real_3_t pi,const cs_real_3_t pia,cs_real_t pir[3],cs_real_t pipr[3])4838 cs_b_cd_steady_vector(const cs_real_t bldfrp,
4839 const double relaxp,
4840 const cs_real_3_t diipb,
4841 const cs_real_33_t gradi,
4842 const cs_real_3_t pi,
4843 const cs_real_3_t pia,
4844 cs_real_t pir[3],
4845 cs_real_t pipr[3])
4846 {
4847 cs_real_3_t recoi;
4848
4849 cs_b_compute_quantities_vector(diipb,
4850 gradi,
4851 bldfrp,
4852 recoi);
4853
4854 cs_b_relax_c_val_vector(relaxp,
4855 pi,
4856 pia,
4857 recoi,
4858 pir,
4859 pipr);
4860 }
4861
4862 /*----------------------------------------------------------------------------*/
4863 /*!
4864 * \brief Handle preparation of boundary face values for the flux computation in
4865 * case of a steady algorithm.
4866 *
4867 * \param[in] bldfrp reconstruction blending factor
4868 * \param[in] relaxp relaxation coefficient
4869 * \param[in] diipb distance I'I'
4870 * \param[in] gradi gradient at cell i
4871 * \param[in] pi value at cell i
4872 * \param[in] pia old value at cell i
4873 * \param[out] pir relaxed value at cell i
4874 * \param[out] pipr relaxed reconstructed value at cell i
4875 */
4876 /*----------------------------------------------------------------------------*/
4877
4878 inline static void
cs_b_cd_steady_tensor(const cs_real_t bldfrp,const double relaxp,const cs_real_3_t diipb,const cs_real_63_t gradi,const cs_real_6_t pi,const cs_real_6_t pia,cs_real_t pir[6],cs_real_t pipr[6])4879 cs_b_cd_steady_tensor(const cs_real_t bldfrp,
4880 const double relaxp,
4881 const cs_real_3_t diipb,
4882 const cs_real_63_t gradi,
4883 const cs_real_6_t pi,
4884 const cs_real_6_t pia,
4885 cs_real_t pir[6],
4886 cs_real_t pipr[6])
4887 {
4888 cs_real_6_t recoi;
4889
4890 cs_b_compute_quantities_tensor(diipb,
4891 gradi,
4892 bldfrp,
4893 recoi);
4894
4895 cs_b_relax_c_val_tensor(relaxp,
4896 pi,
4897 pia,
4898 recoi,
4899 pir,
4900 pipr);
4901 }
4902
4903 /*----------------------------------------------------------------------------*/
4904 /*!
4905 * \brief Handle preparation of boundary face values for the flux computation in
4906 * case of an unsteady algorithm.
4907 *
4908 * \param[in] bldfrp reconstruction blending factor
4909 * \param[in] diipb distance I'I'
4910 * \param[in] gradi gradient at cell i
4911 * \param[in] pi value at cell i
4912 * \param[out] pip reconstructed value at cell i
4913 */
4914 /*----------------------------------------------------------------------------*/
4915
4916 inline static void
cs_b_cd_unsteady(const cs_real_t bldfrp,const cs_real_3_t diipb,const cs_real_3_t gradi,const cs_real_t pi,cs_real_t * pip)4917 cs_b_cd_unsteady(const cs_real_t bldfrp,
4918 const cs_real_3_t diipb,
4919 const cs_real_3_t gradi,
4920 const cs_real_t pi,
4921 cs_real_t *pip)
4922 {
4923 cs_real_t recoi;
4924
4925 cs_b_compute_quantities(diipb,
4926 gradi,
4927 bldfrp,
4928 &recoi);
4929
4930 *pip = pi + recoi;
4931 }
4932
4933 /*----------------------------------------------------------------------------*/
4934 /*!
4935 * \brief Handle preparation of boundary face values for the flux computation in
4936 * case of a steady algorithm.
4937 *
4938 * \param[in] bldfrp reconstruction blending factor
4939 * \param[in] diipb distance I'I'
4940 * \param[in] gradi gradient at cell i
4941 * \param[in] pi value at cell i
4942 * \param[out] pip reconstructed value at cell i
4943 */
4944 /*----------------------------------------------------------------------------*/
4945
4946 inline static void
cs_b_cd_unsteady_vector(const cs_real_t bldfrp,const cs_real_3_t diipb,const cs_real_33_t gradi,const cs_real_3_t pi,cs_real_t pip[3])4947 cs_b_cd_unsteady_vector(const cs_real_t bldfrp,
4948 const cs_real_3_t diipb,
4949 const cs_real_33_t gradi,
4950 const cs_real_3_t pi,
4951 cs_real_t pip[3])
4952 {
4953 cs_real_3_t recoi;
4954
4955 cs_b_compute_quantities_vector(diipb,
4956 gradi,
4957 bldfrp,
4958 recoi);
4959
4960 for (int isou = 0; isou < 3; isou++)
4961 pip[isou] = pi[isou] + recoi[isou];
4962 }
4963
4964 /*----------------------------------------------------------------------------*/
4965 /*!
4966 * \brief Handle preparation of boundary face values for the flux computation in
4967 * case of a steady algorithm.
4968 *
4969 * \param[in] bldfrp reconstruction blending factor
4970 * \param[in] diipb distance I'I'
4971 * \param[in] gradi gradient at cell i
4972 * \param[in] pi value at cell i
4973 * \param[out] pip reconstructed value at cell i
4974 */
4975 /*----------------------------------------------------------------------------*/
4976
4977 inline static void
cs_b_cd_unsteady_tensor(const cs_real_t bldfrp,const cs_real_3_t diipb,const cs_real_63_t gradi,const cs_real_6_t pi,cs_real_t pip[6])4978 cs_b_cd_unsteady_tensor(const cs_real_t bldfrp,
4979 const cs_real_3_t diipb,
4980 const cs_real_63_t gradi,
4981 const cs_real_6_t pi,
4982 cs_real_t pip[6])
4983 {
4984 cs_real_6_t recoi;
4985
4986 cs_b_compute_quantities_tensor(diipb,
4987 gradi,
4988 bldfrp,
4989 recoi);
4990
4991 for(int isou = 0; isou< 6; isou++)
4992 pip[isou] = pi[isou] + recoi[isou];
4993 }
4994
4995 /*----------------------------------------------------------------------------*/
4996 /*!
4997 * \brief Add diffusive flux to flux at an internal coupling face.
4998 *
4999 * \param[in] idiffp diffusion flag
5000 * \param[in] pi value at cell i
5001 * \param[in] pj value at cell j
5002 * \param[in] b_visc equivalent exchange coefficient at an internal
5003 * coupling face
5004 * \param[in,out] fluxi flux at internal coupling face
5005 */
5006 /*----------------------------------------------------------------------------*/
5007
5008 inline static void
cs_b_diff_flux_coupling(int idiffp,cs_real_t pi,cs_real_t pj,cs_real_t b_visc,cs_real_t * fluxi)5009 cs_b_diff_flux_coupling(int idiffp,
5010 cs_real_t pi,
5011 cs_real_t pj,
5012 cs_real_t b_visc,
5013 cs_real_t *fluxi)
5014 {
5015 *fluxi += idiffp*b_visc*(pi - pj);
5016 }
5017
5018 /*----------------------------------------------------------------------------*/
5019 /*!
5020 * \brief Add diffusive flux to flux at an internal coupling face for a vector.
5021 *
5022 * \param[in] idiffp diffusion flag
5023 * \param[in] pi value at cell i
5024 * \param[in] pj value at cell j
5025 * \param[in] b_visc equivalent exchange coefficient at an internal
5026 * coupling face
5027 * \param[in,out] fluxi flux at internal coupling face
5028 */
5029 /*----------------------------------------------------------------------------*/
5030
5031 inline static void
cs_b_diff_flux_coupling_vector(int idiffp,const cs_real_t pi[3],const cs_real_t pj[3],cs_real_t b_visc,cs_real_t fluxi[3])5032 cs_b_diff_flux_coupling_vector(int idiffp,
5033 const cs_real_t pi[3],
5034 const cs_real_t pj[3],
5035 cs_real_t b_visc,
5036 cs_real_t fluxi[3])
5037 {
5038 for (int k = 0; k < 3; k++)
5039 fluxi[k] += idiffp*b_visc*(pi[k] - pj[k]);
5040 }
5041
5042
5043 /*============================================================================
5044 * Public function prototypes for Fortran API
5045 *============================================================================*/
5046
5047 /*----------------------------------------------------------------------------
5048 * Wrapper to cs_face_diffusion_potential
5049 *----------------------------------------------------------------------------*/
5050
5051 void CS_PROCF (itrmas, ITRMAS)
5052 (
5053 const int *const f_id,
5054 const int *const init,
5055 const int *const inc,
5056 const int *const imrgra,
5057 const int *const iccocg,
5058 const int *const nswrgp,
5059 const int *const imligp,
5060 const int *const iphydp,
5061 const int *const iwgrp,
5062 const int *const iwarnp,
5063 const cs_real_t *const epsrgp,
5064 const cs_real_t *const climgp,
5065 const cs_real_t *const extrap,
5066 cs_real_3_t frcxt[],
5067 cs_real_t pvar[],
5068 const cs_real_t coefap[],
5069 const cs_real_t coefbp[],
5070 const cs_real_t cofafp[],
5071 const cs_real_t cofbfp[],
5072 const cs_real_t i_visc[],
5073 const cs_real_t b_visc[],
5074 cs_real_t visel[],
5075 cs_real_t i_massflux[],
5076 cs_real_t b_massflux[]
5077 );
5078
5079 /*----------------------------------------------------------------------------
5080 * Wrapper to cs_face_anisotropic_diffusion_potential
5081 *----------------------------------------------------------------------------*/
5082
5083 void CS_PROCF (itrmav, ITRMAV)
5084 (
5085 const int *const f_id,
5086 const int *const init,
5087 const int *const inc,
5088 const int *const imrgra,
5089 const int *const iccocg,
5090 const int *const nswrgp,
5091 const int *const imligp,
5092 const int *const ircflp,
5093 const int *const iphydp,
5094 const int *const iwgrp,
5095 const int *const iwarnp,
5096 const cs_real_t *const epsrgp,
5097 const cs_real_t *const climgp,
5098 const cs_real_t *const extrap,
5099 cs_real_3_t frcxt[],
5100 cs_real_t pvar[],
5101 const cs_real_t coefap[],
5102 const cs_real_t coefbp[],
5103 const cs_real_t cofafp[],
5104 const cs_real_t cofbfp[],
5105 const cs_real_t i_visc[],
5106 const cs_real_t b_visc[],
5107 cs_real_6_t viscel[],
5108 const cs_real_2_t weighf[],
5109 const cs_real_t weighb[],
5110 cs_real_t i_massflux[],
5111 cs_real_t b_massflux[]
5112 );
5113
5114 /*----------------------------------------------------------------------------
5115 * Wrapper to cs_diffusion_potential
5116 *----------------------------------------------------------------------------*/
5117
5118 void CS_PROCF (itrgrp, ITRGRP)
5119 (
5120 const int *const f_id,
5121 const int *const init,
5122 const int *const inc,
5123 const int *const imrgra,
5124 const int *const iccocg,
5125 const int *const nswrgp,
5126 const int *const imligp,
5127 const int *const iphydp,
5128 const int *const iwgrp,
5129 const int *const iwarnp,
5130 const cs_real_t *const epsrgp,
5131 const cs_real_t *const climgp,
5132 const cs_real_t *const extrap,
5133 cs_real_3_t frcxt[],
5134 cs_real_t pvar[],
5135 const cs_real_t coefap[],
5136 const cs_real_t coefbp[],
5137 const cs_real_t cofafp[],
5138 const cs_real_t cofbfp[],
5139 const cs_real_t i_visc[],
5140 const cs_real_t b_visc[],
5141 cs_real_t visel[],
5142 cs_real_t diverg[]
5143 );
5144
5145 /*----------------------------------------------------------------------------
5146 * Wrapper to cs_anisotropic_diffusion_potential
5147 *----------------------------------------------------------------------------*/
5148
5149 void CS_PROCF (itrgrv, ITRGRV)
5150 (
5151 const int *const f_id,
5152 const int *const init,
5153 const int *const inc,
5154 const int *const imrgra,
5155 const int *const iccocg,
5156 const int *const nswrgp,
5157 const int *const imligp,
5158 const int *const ircflp,
5159 const int *const iphydp,
5160 const int *const iwgrp,
5161 const int *const iwarnp,
5162 const cs_real_t *const epsrgp,
5163 const cs_real_t *const climgp,
5164 const cs_real_t *const extrap,
5165 cs_real_3_t frcxt[],
5166 cs_real_t pvar[],
5167 const cs_real_t coefap[],
5168 const cs_real_t coefbp[],
5169 const cs_real_t cofafp[],
5170 const cs_real_t cofbfp[],
5171 const cs_real_t i_visc[],
5172 const cs_real_t b_visc[],
5173 cs_real_6_t viscel[],
5174 const cs_real_2_t weighf[],
5175 const cs_real_t weighb[],
5176 cs_real_t diverg[]
5177 );
5178
5179 /*=============================================================================
5180 * Public function prototypes
5181 *============================================================================*/
5182
5183 /*----------------------------------------------------------------------------*/
5184 /*!
5185 * \brief Compute the upwind gradient used in the slope tests.
5186 *
5187 * This function assumes the input gradient and pvar values have already
5188 * been synchronized.
5189 *
5190 * \param[in] f_id field id
5191 * \param[in] inc Not an increment flag
5192 * \param[in] halo_type halo type
5193 * \param[in] grad standard gradient
5194 * \param[out] grdpa upwind gradient
5195 * \param[in] pvar values
5196 * \param[in] coefap boundary condition array for the variable
5197 * (explicit part)
5198 * \param[in] coefbp boundary condition array for the variable
5199 * (implicit part)
5200 * \param[in] i_massflux mass flux at interior faces
5201 */
5202 /*----------------------------------------------------------------------------*/
5203
5204 void
5205 cs_slope_test_gradient(int f_id,
5206 int inc,
5207 cs_halo_type_t halo_type,
5208 const cs_real_3_t *grad,
5209 cs_real_3_t *grdpa,
5210 const cs_real_t *pvar,
5211 const cs_real_t *coefap,
5212 const cs_real_t *coefbp,
5213 const cs_real_t *i_massflux);
5214
5215 /*----------------------------------------------------------------------------*/
5216 /*!
5217 * \brief Compute the upwind gradient in order to cope with SOLU schemes
5218 * observed in the litterature.
5219 *
5220 * \param[in] f_id field index
5221 * \param[in] inc Not an increment flag
5222 * \param[in] halo_type halo type
5223 * \param[out] grdpa upwind gradient
5224 * \param[in] pvar values
5225 * \param[in] coefap boundary condition array for the variable
5226 * (explicit part)
5227 * \param[in] coefbp boundary condition array for the variable
5228 * (implicit part)
5229 * \param[in] i_massflux mass flux at interior faces
5230 */
5231 /*----------------------------------------------------------------------------*/
5232
5233 void
5234 cs_upwind_gradient(const int f_id,
5235 const int inc,
5236 const cs_halo_type_t halo_type,
5237 const cs_real_t coefap[],
5238 const cs_real_t coefbp[],
5239 const cs_real_t i_massflux[],
5240 const cs_real_t b_massflux[],
5241 const cs_real_t *restrict pvar,
5242 cs_real_3_t *restrict grdpa);
5243
5244 /*----------------------------------------------------------------------------*/
5245 /*!
5246 * \brief Compute the upwind gradient used in the slope tests.
5247 *
5248 * This function assumes the input gradient and pvar values have already
5249 * been synchronized.
5250 *
5251 * \param[in] inc Not an increment flag
5252 * \param[in] halo_type halo type
5253 * \param[in] grad standard gradient
5254 * \param[out] grdpa upwind gradient
5255 * \param[in] pvar values
5256 * \param[in] coefa boundary condition array for the variable
5257 * (explicit part)
5258 * \param[in] coefb boundary condition array for the variable
5259 * (implicit part)
5260 * \param[in] i_massflux mass flux at interior faces
5261 */
5262 /*----------------------------------------------------------------------------*/
5263
5264 void
5265 cs_slope_test_gradient_vector(const int inc,
5266 const cs_halo_type_t halo_type,
5267 const cs_real_33_t *grad,
5268 cs_real_33_t *grdpa,
5269 const cs_real_3_t *pvar,
5270 const cs_real_3_t *coefa,
5271 const cs_real_33_t *coefb,
5272 const cs_real_t *i_massflux);
5273
5274 /*----------------------------------------------------------------------------*/
5275 /*!
5276 * \brief Compute the upwind gradient used in the slope tests.
5277 *
5278 * This function assumes the input gradient and pvar values have already
5279 * been synchronized.
5280 *
5281 * \param[in] inc Not an increment flag
5282 * \param[in] halo_type halo type
5283 * \param[in] grad standard gradient
5284 * \param[out] grdpa upwind gradient
5285 * \param[in] pvar values
5286 * \param[in] coefa boundary condition array for the variable
5287 * (explicit part)
5288 * \param[in] coefb boundary condition array for the variable
5289 * (implicit part)
5290 * \param[in] i_massflux mass flux at interior faces
5291 */
5292 /*----------------------------------------------------------------------------*/
5293
5294 void
5295 cs_slope_test_gradient_tensor(const int inc,
5296 const cs_halo_type_t halo_type,
5297 const cs_real_63_t *grad,
5298 cs_real_63_t *grdpa,
5299 const cs_real_6_t *pvar,
5300 const cs_real_6_t *coefa,
5301 const cs_real_66_t *coefb,
5302 const cs_real_t *i_massflux);
5303
5304 /*----------------------------------------------------------------------------*/
5305 /*!
5306 * \brief Compute the beta blending coefficient of the
5307 * beta limiter (ensuring preservation of a given min/max pair of values).
5308 *
5309 * \param[in] f_id field id
5310 * \param[in] inc "not an increment" flag
5311 * \param[in] rovsdt rho * volume / dt
5312 */
5313 /*----------------------------------------------------------------------------*/
5314
5315 void
5316 cs_beta_limiter_building(int f_id,
5317 int inc,
5318 const cs_real_t rovsdt[]);
5319
5320 /*----------------------------------------------------------------------------*/
5321 /*!
5322 * \brief Add the explicit part of the convection/diffusion terms of a
5323 * standard transport equation of a scalar field \f$ \varia \f$.
5324 *
5325 * More precisely, the right hand side \f$ Rhs \f$ is updated as
5326 * follows:
5327 * \f[
5328 * Rhs = Rhs - \sum_{\fij \in \Facei{\celli}} \left(
5329 * \dot{m}_\ij \left( \varia_\fij - \varia_\celli \right)
5330 * - \mu_\fij \gradv_\fij \varia \cdot \vect{S}_\ij \right)
5331 * \f]
5332 *
5333 * Warning:
5334 * - \f$ Rhs \f$ has already been initialized before calling bilsc2!
5335 * - mind the sign minus
5336 *
5337 * \param[in] idtvar indicator of the temporal scheme
5338 * \param[in] f_id field id (or -1)
5339 * \param[in] var_cal_opt variable calculation options
5340 * \param[in] icvflb global indicator of boundary convection flux
5341 * - 0 upwind scheme at all boundary faces
5342 * - 1 imposed flux at some boundary faces
5343 * \param[in] inc indicator
5344 * - 0 when solving an increment
5345 * - 1 otherwise
5346 * \param[in] iccocg indicator
5347 * - 1 re-compute cocg matrix
5348 * (for iterative gradients)
5349 * - 0 otherwise
5350 * \param[in] imasac take mass accumulation into account?
5351 * \param[in] pvar solved variable (current time step)
5352 * \param[in] pvara solved variable (previous time step)
5353 * \param[in] icvfli boundary face indicator array of convection flux
5354 * - 0 upwind scheme
5355 * - 1 imposed flux
5356 * \param[in] coefap boundary condition array for the variable
5357 * (explicit part)
5358 * \param[in] coefbp boundary condition array for the variable
5359 * (implicit part)
5360 * \param[in] cofafp boundary condition array for the diffusion
5361 * of the variable (explicit part)
5362 * \param[in] cofbfp boundary condition array for the diffusion
5363 * of the variable (implicit part)
5364 * \param[in] i_massflux mass flux at interior faces
5365 * \param[in] b_massflux mass flux at boundary faces
5366 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
5367 * at interior faces for the r.h.s.
5368 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
5369 * at border faces for the r.h.s.
5370 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5371 */
5372 /*----------------------------------------------------------------------------*/
5373
5374 void
5375 cs_convection_diffusion_scalar(int idtvar,
5376 int f_id,
5377 const cs_var_cal_opt_t var_cal_opt,
5378 int icvflb,
5379 int inc,
5380 int iccocg,
5381 int imasac,
5382 cs_real_t *restrict pvar,
5383 const cs_real_t *restrict pvara,
5384 const int icvfli[],
5385 const cs_real_t coefap[],
5386 const cs_real_t coefbp[],
5387 const cs_real_t cofafp[],
5388 const cs_real_t cofbfp[],
5389 const cs_real_t i_massflux[],
5390 const cs_real_t b_massflux[],
5391 const cs_real_t i_visc[],
5392 const cs_real_t b_visc[],
5393 cs_real_t *restrict rhs);
5394
5395 /*----------------------------------------------------------------------------*/
5396 /*!
5397 * <a name="cs_face_convection_scalar"></a>
5398 *
5399 * \brief Update face flux with convection contribution of a standard transport
5400 * equation of a scalar field \f$ \varia \f$.
5401 *
5402 * \f[
5403 * C_\ij = \dot{m}_\ij \left( \varia_\fij - \varia_\celli \right)
5404 * \f]
5405 *
5406 * \param[in] idtvar indicator of the temporal scheme
5407 * \param[in] f_id field id (or -1)
5408 * \param[in] var_cal_opt variable calculation options
5409 * \param[in] icvflb global indicator of boundary convection flux
5410 * - 0 upwind scheme at all boundary faces
5411 * - 1 imposed flux at some boundary faces
5412 * \param[in] inc indicator
5413 * - 0 when solving an increment
5414 * - 1 otherwise
5415 * \param[in] iccocg indicator
5416 * - 1 re-compute cocg matrix
5417 * (for iterative gradients)
5418 * - 0 otherwise
5419 * \param[in] imasac take mass accumulation into account?
5420 * \param[in] pvar solved variable (current time step)
5421 * \param[in] pvara solved variable (previous time step)
5422 * \param[in] icvfli boundary face indicator array of convection flux
5423 * - 0 upwind scheme
5424 * - 1 imposed flux
5425 * \param[in] coefap boundary condition array for the variable
5426 * (explicit part)
5427 * \param[in] coefbp boundary condition array for the variable
5428 * (implicit part)
5429 * \param[in] i_massflux mass flux at interior faces
5430 * \param[in] b_massflux mass flux at boundary faces
5431 * \param[in,out] i_conv_flux scalar convection flux at interior faces
5432 * \param[in,out] b_conv_flux scalar convection flux at boundary faces
5433 */
5434 /*----------------------------------------------------------------------------*/
5435
5436 void
5437 cs_face_convection_scalar(int idtvar,
5438 int f_id,
5439 const cs_var_cal_opt_t var_cal_opt,
5440 int icvflb,
5441 int inc,
5442 int iccocg,
5443 int imasac,
5444 cs_real_t *restrict pvar,
5445 const cs_real_t *restrict pvara,
5446 const int icvfli[],
5447 const cs_real_t coefap[],
5448 const cs_real_t coefbp[],
5449 const cs_real_t i_massflux[],
5450 const cs_real_t b_massflux[],
5451 cs_real_2_t i_conv_flux[],
5452 cs_real_t b_conv_flux[]);
5453
5454 /*----------------------------------------------------------------------------*/
5455 /*!
5456 * \brief Add the explicit part of the convection/diffusion terms of a transport
5457 * equation of a vector field \f$ \vect{\varia} \f$.
5458 *
5459 * More precisely, the right hand side \f$ \vect{Rhs} \f$ is updated as
5460 * follows:
5461 * \f[
5462 * \vect{Rhs} = \vect{Rhs} - \sum_{\fij \in \Facei{\celli}} \left(
5463 * \dot{m}_\ij \left( \vect{\varia}_\fij - \vect{\varia}_\celli \right)
5464 * - \mu_\fij \gradt_\fij \vect{\varia} \cdot \vect{S}_\ij \right)
5465 * \f]
5466 *
5467 * Remark:
5468 * if ivisep = 1, then we also take \f$ \mu \transpose{\gradt\vect{\varia}}
5469 * + \lambda \trace{\gradt\vect{\varia}} \f$, where \f$ \lambda \f$ is
5470 * the secondary viscosity, i.e. usually \f$ -\frac{2}{3} \mu \f$.
5471 *
5472 * Warning:
5473 * - \f$ \vect{Rhs} \f$ has already been initialized before calling bilsc!
5474 * - mind the sign minus
5475 *
5476 * \param[in] idtvar indicator of the temporal scheme
5477 * \param[in] f_id index of the current variable
5478 * \param[in] var_cal_opt variable calculation options
5479 * \param[in] icvflb global indicator of boundary convection flux
5480 * - 0 upwind scheme at all boundary faces
5481 * - 1 imposed flux at some boundary faces
5482 * \param[in] inc indicator
5483 * - 0 when solving an increment
5484 * - 1 otherwise
5485 * \param[in] ivisep indicator to take \f$ \divv
5486 * \left(\mu \gradt \transpose{\vect{a}} \right)
5487 * -2/3 \grad\left( \mu \dive \vect{a} \right)\f$
5488 * - 1 take into account,
5489 * - 0 otherwise
5490 * \param[in] imasac take mass accumulation into account?
5491 * \param[in] pvar solved velocity (current time step)
5492 * \param[in] pvara solved velocity (previous time step)
5493 * \param[in] icvfli boundary face indicator array of convection flux
5494 * - 0 upwind scheme
5495 * - 1 imposed flux
5496 * \param[in] coefav boundary condition array for the variable
5497 * (explicit part)
5498 * \param[in] coefbv boundary condition array for the variable
5499 * (implicit part)
5500 * \param[in] cofafv boundary condition array for the diffusion
5501 * of the variable (explicit part)
5502 * \param[in] cofbfv boundary condition array for the diffusion
5503 * of the variable (implicit part)
5504 * \param[in] i_massflux mass flux at interior faces
5505 * \param[in] b_massflux mass flux at boundary faces
5506 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
5507 * at interior faces for the r.h.s.
5508 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
5509 * at border faces for the r.h.s.
5510 * \param[in] i_secvis secondary viscosity at interior faces
5511 * \param[in] b_secvis secondary viscosity at boundary faces
5512 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5513 */
5514 /*----------------------------------------------------------------------------*/
5515
5516 void
5517 cs_convection_diffusion_vector(int idtvar,
5518 int f_id,
5519 const cs_var_cal_opt_t var_cal_opt,
5520 int icvflb,
5521 int inc,
5522 int ivisep,
5523 int imasac,
5524 cs_real_3_t *restrict pvar,
5525 const cs_real_3_t *restrict pvara,
5526 const int icvfli[],
5527 const cs_real_3_t coefav[],
5528 const cs_real_33_t coefbv[],
5529 const cs_real_3_t cofafv[],
5530 const cs_real_33_t cofbfv[],
5531 const cs_real_t i_massflux[],
5532 const cs_real_t b_massflux[],
5533 const cs_real_t i_visc[],
5534 const cs_real_t b_visc[],
5535 const cs_real_t i_secvis[],
5536 const cs_real_t b_secvis[],
5537 cs_real_3_t *restrict rhs);
5538
5539 /*----------------------------------------------------------------------------*/
5540 /*!
5541 * \brief Add the explicit part of the convection/diffusion terms of a transport
5542 * equation of a vector field \f$ \vect{\varia} \f$.
5543 *
5544 * More precisely, the right hand side \f$ \vect{Rhs} \f$ is updated as
5545 * follows:
5546 * \f[
5547 * \vect{Rhs} = \vect{Rhs} - \sum_{\fij \in \Facei{\celli}} \left(
5548 * \dot{m}_\ij \left( \vect{\varia}_\fij - \vect{\varia}_\celli \right)
5549 * - \mu_\fij \gradt_\fij \vect{\varia} \cdot \vect{S}_\ij \right)
5550 * \f]
5551 *
5552 * Warning:
5553 * - \f$ \vect{Rhs} \f$ has already been initialized before calling bilsc!
5554 * - mind the sign minus
5555 *
5556 * \param[in] idtvar indicator of the temporal scheme
5557 * \param[in] f_id index of the current variable
5558 * \param[in] var_cal_opt variable calculation options
5559 * \param[in] icvflb global indicator of boundary convection flux
5560 * - 0 upwind scheme at all boundary faces
5561 * - 1 imposed flux at some boundary faces
5562 * \param[in] inc indicator
5563 * - 0 when solving an increment
5564 * - 1 otherwise
5565 * \param[in] imasac take mass accumulation into account?
5566 * \param[in] pvar solved velocity (current time step)
5567 * \param[in] pvara solved velocity (previous time step)
5568 * \param[in] coefa boundary condition array for the variable
5569 * (Explicit part)
5570 * \param[in] coefb boundary condition array for the variable
5571 * (Implicit part)
5572 * \param[in] cofaf boundary condition array for the diffusion
5573 * of the variable (Explicit part)
5574 * \param[in] cofbf boundary condition array for the diffusion
5575 * of the variable (Implicit part)
5576 * \param[in] i_massflux mass flux at interior faces
5577 * \param[in] b_massflux mass flux at boundary faces
5578 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
5579 * at interior faces for the r.h.s.
5580 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
5581 * at border faces for the r.h.s.
5582 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5583 */
5584 /*----------------------------------------------------------------------------*/
5585
5586 void
5587 cs_convection_diffusion_tensor(int idtvar,
5588 int f_id,
5589 const cs_var_cal_opt_t var_cal_opt,
5590 int icvflb,
5591 int inc,
5592 int imasac,
5593 cs_real_6_t *restrict pvar,
5594 const cs_real_6_t *restrict pvara,
5595 const cs_real_6_t coefa[],
5596 const cs_real_66_t coefb[],
5597 const cs_real_6_t cofaf[],
5598 const cs_real_66_t cofbf[],
5599 const cs_real_t i_massflux[],
5600 const cs_real_t b_massflux[],
5601 const cs_real_t i_visc[],
5602 const cs_real_t b_visc[],
5603 cs_real_6_t *restrict rhs);
5604
5605 /*----------------------------------------------------------------------------*/
5606 /*!
5607 * \brief Add the explicit part of the convection/diffusion terms of a transport
5608 * equation of a scalar field \f$ \varia \f$ such as the temperature.
5609 *
5610 * More precisely, the right hand side \f$ Rhs \f$ is updated as
5611 * follows:
5612 * \f[
5613 * Rhs = Rhs + \sum_{\fij \in \Facei{\celli}} \left(
5614 * C_p\dot{m}_\ij \varia_\fij
5615 * - \lambda_\fij \gradv_\fij \varia \cdot \vect{S}_\ij \right)
5616 * \f]
5617 *
5618 * Warning:
5619 * \f$ Rhs \f$ has already been initialized before calling bilsct!
5620 *
5621 * \param[in] idtvar indicator of the temporal scheme
5622 * \param[in] f_id index of the current variable
5623 * \param[in] var_cal_opt variable calculation options
5624 * \param[in] inc indicator
5625 * - 0 when solving an increment
5626 * - 1 otherwise
5627 * \param[in] iccocg indicator
5628 * - 1 re-compute cocg matrix
5629 * (for iterative gradients)
5630 * - 0 otherwise
5631 * \param[in] imasac take mass accumulation into account?
5632 * \param[in] pvar solved variable (current time step)
5633 * \param[in] pvara solved variable (previous time step)
5634 * \param[in] coefap boundary condition array for the variable
5635 * (explicit part)
5636 * \param[in] coefbp boundary condition array for the variable
5637 * (implicit part)
5638 * \param[in] cofafp boundary condition array for the diffusion
5639 * of the variable (explicit part)
5640 * \param[in] cofbfp boundary condition array for the diffusion
5641 * of the variable (implicit part)
5642 * \param[in] i_massflux mass flux at interior faces
5643 * \param[in] b_massflux mass flux at boundary faces
5644 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
5645 * at interior faces for the r.h.s.
5646 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
5647 * at border faces for the r.h.s.
5648 * \param[in] xcpp array of specific heat (\f$ C_p \f$)
5649 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5650 */
5651 /*----------------------------------------------------------------------------*/
5652
5653 void
5654 cs_convection_diffusion_thermal(int idtvar,
5655 int f_id,
5656 const cs_var_cal_opt_t var_cal_opt,
5657 int inc,
5658 int iccocg,
5659 int imasac,
5660 cs_real_t *restrict pvar,
5661 const cs_real_t *restrict pvara,
5662 const cs_real_t coefap[],
5663 const cs_real_t coefbp[],
5664 const cs_real_t cofafp[],
5665 const cs_real_t cofbfp[],
5666 const cs_real_t i_massflux[],
5667 const cs_real_t b_massflux[],
5668 const cs_real_t i_visc[],
5669 const cs_real_t b_visc[],
5670 const cs_real_t xcpp[],
5671 cs_real_t *restrict rhs);
5672
5673 /*----------------------------------------------------------------------------*/
5674 /*!
5675 * \brief Add the explicit part of the diffusion terms with a symmetric tensor
5676 * diffusivity for a transport equation of a scalar field \f$ \varia \f$.
5677 *
5678 * More precisely, the right hand side \f$ Rhs \f$ is updated as
5679 * follows:
5680 * \f[
5681 * Rhs = Rhs - \sum_{\fij \in \Facei{\celli}} \left(
5682 * - \tens{\mu}_\fij \gradv_\fij \varia \cdot \vect{S}_\ij \right)
5683 * \f]
5684 *
5685 * Warning:
5686 * - \f$ Rhs \f$ has already been initialized before
5687 * calling cs_anisotropic_diffusion_scalar!
5688 * - mind the sign minus
5689 *
5690 * \param[in] idtvar indicator of the temporal scheme
5691 * \param[in] f_id index of the current variable
5692 * \param[in] var_cal_opt variable calculation options
5693 * \param[in] inc indicator
5694 * - 0 when solving an increment
5695 * - 1 otherwise
5696 * \param[in] iccocg indicator
5697 * - 1 re-compute cocg matrix
5698 (for iterativ gradients)
5699 * - 0 otherwise
5700 * \param[in] pvar solved variable (current time step)
5701 * \param[in] pvara solved variable (previous time step)
5702 * \param[in] coefap boundary condition array for the variable
5703 * (explicit part)
5704 * \param[in] coefbp boundary condition array for the variable
5705 * (implicit part)
5706 * \param[in] cofafp boundary condition array for the diffusion
5707 * of the variable (explicit part)
5708 * \param[in] cofbfp boundary condition array for the diffusion
5709 * of the variable (implicit part)
5710 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
5711 * at interior faces for the r.h.s.
5712 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
5713 * at border faces for the r.h.s.
5714 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
5715 * \param[in] weighf internal face weight between cells i j in case
5716 * of tensor diffusion
5717 * \param[in] weighb boundary face weight for cells i in case
5718 * of tensor diffusion
5719 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5720 */
5721 /*----------------------------------------------------------------------------*/
5722
5723 void
5724 cs_anisotropic_diffusion_scalar(int idtvar,
5725 int f_id,
5726 const cs_var_cal_opt_t var_cal_opt,
5727 int inc,
5728 int iccocg,
5729 cs_real_t *restrict pvar,
5730 const cs_real_t *restrict pvara,
5731 const cs_real_t coefap[],
5732 const cs_real_t coefbp[],
5733 const cs_real_t cofafp[],
5734 const cs_real_t cofbfp[],
5735 const cs_real_t i_visc[],
5736 const cs_real_t b_visc[],
5737 cs_real_6_t *restrict viscel,
5738 const cs_real_2_t weighf[],
5739 const cs_real_t weighb[],
5740 cs_real_t *restrict rhs);
5741
5742 /*-----------------------------------------------------------------------------*/
5743 /*!
5744 * \brief Add explicit part of the terms of diffusion by a left-multiplying
5745 * symmetric tensorial diffusivity for a transport equation of a vector field
5746 * \f$ \vect{\varia} \f$.
5747 *
5748 * More precisely, the right hand side \f$ \vect{Rhs} \f$ is updated as
5749 * follows:
5750 * \f[
5751 * \vect{Rhs} = \vect{Rhs} - \sum_{\fij \in \Facei{\celli}} \left(
5752 * - \gradt_\fij \vect{\varia} \tens{\mu}_\fij \cdot \vect{S}_\ij \right)
5753 * \f]
5754 *
5755 * Remark:
5756 * if ivisep = 1, then we also take \f$ \mu \transpose{\gradt\vect{\varia}}
5757 * + \lambda \trace{\gradt\vect{\varia}} \f$, where \f$ \lambda \f$ is
5758 * the secondary viscosity, i.e. usually \f$ -\frac{2}{3} \mu \f$.
5759 *
5760 * Warning:
5761 * - \f$ \vect{Rhs} \f$ has already been initialized before calling the present
5762 * function
5763 * - mind the sign minus
5764 *
5765 * \param[in] idtvar indicator of the temporal scheme
5766 * \param[in] f_id index of the current variable
5767 * \param[in] var_cal_opt variable calculation options
5768 * \param[in] inc indicator
5769 * - 0 when solving an increment
5770 * - 1 otherwise
5771 * \param[in] ivisep indicator to take \f$ \divv
5772 * \left(\mu \gradt \transpose{\vect{a}} \right)
5773 * -2/3 \grad\left( \mu \dive \vect{a} \right)\f$
5774 * - 1 take into account,
5775 * \param[in] pvar solved variable (current time step)
5776 * \param[in] pvara solved variable (previous time step)
5777 * \param[in] coefav boundary condition array for the variable
5778 * (explicit part)
5779 * \param[in] coefbv boundary condition array for the variable
5780 * (implicit part)
5781 * \param[in] cofafv boundary condition array for the diffusion
5782 * of the variable (explicit part)
5783 * \param[in] cofbfv boundary condition array for the diffusion
5784 * of the variable (implicit part)
5785 * \param[in] i_visc \f$ \tens{\mu}_\fij \dfrac{S_\fij}{\ipf\jpf} \f$
5786 * at interior faces for the r.h.s.
5787 * \param[in] b_visc \f$ \dfrac{S_\fib}{\ipf \centf} \f$
5788 * at border faces for the r.h.s.
5789 * \param[in] i_secvis secondary viscosity at interior faces
5790 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5791 */
5792 /*----------------------------------------------------------------------------*/
5793
5794 void
5795 cs_anisotropic_left_diffusion_vector(int idtvar,
5796 int f_id,
5797 const cs_var_cal_opt_t var_cal_opt,
5798 int inc,
5799 int ivisep,
5800 cs_real_3_t *restrict pvar,
5801 const cs_real_3_t *restrict pvara,
5802 const cs_real_3_t coefav[],
5803 const cs_real_33_t coefbv[],
5804 const cs_real_3_t cofafv[],
5805 const cs_real_33_t cofbfv[],
5806 const cs_real_33_t i_visc[],
5807 const cs_real_t b_visc[],
5808 const cs_real_t i_secvis[],
5809 cs_real_3_t *restrict rhs);
5810
5811 /*-----------------------------------------------------------------------------*/
5812 /*!
5813 * \brief Add explicit part of the terms of diffusion by a right-multiplying
5814 * symmetric tensorial diffusivity for a transport equation of a vector field
5815 * \f$ \vect{\varia} \f$.
5816 *
5817 * More precisely, the right hand side \f$ \vect{Rhs} \f$ is updated as
5818 * follows:
5819 * \f[
5820 * \vect{Rhs} = \vect{Rhs} - \sum_{\fij \in \Facei{\celli}} \left(
5821 * - \gradt_\fij \vect{\varia} \tens{\mu}_\fij \cdot \vect{S}_\ij \right)
5822 * \f]
5823 *
5824 * Warning:
5825 * - \f$ \vect{Rhs} \f$ has already been initialized before calling the present
5826 * function
5827 * - mind the sign minus
5828 *
5829 * \param[in] idtvar indicator of the temporal scheme
5830 * \param[in] f_id index of the current variable
5831 * \param[in] var_cal_opt variable calculation options
5832 * \param[in] inc indicator
5833 * - 0 when solving an increment
5834 * - 1 otherwise
5835 * \param[in] pvar solved variable (current time step)
5836 * \param[in] pvara solved variable (previous time step)
5837 * \param[in] coefav boundary condition array for the variable
5838 * (explicit part)
5839 * \param[in] coefbv boundary condition array for the variable
5840 * (implicit part)
5841 * \param[in] cofafv boundary condition array for the diffusion
5842 * of the variable (explicit part)
5843 * \param[in] cofbfv boundary condition array for the diffusion
5844 * of the variable (implicit part)
5845 * \param[in] i_visc \f$ \tens{\mu}_\fij \dfrac{S_\fij}{\ipf\jpf} \f$
5846 * at interior faces for the r.h.s.
5847 * \param[in] b_visc \f$ \dfrac{S_\fib}{\ipf \centf} \f$
5848 * at border faces for the r.h.s.
5849 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
5850 * \param[in] weighf internal face weight between cells i j in case
5851 * of tensor diffusion
5852 * \param[in] weighb boundary face weight for cells i in case
5853 * of tensor diffusion
5854 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5855 */
5856 /*----------------------------------------------------------------------------*/
5857
5858 void
5859 cs_anisotropic_right_diffusion_vector(int idtvar,
5860 int f_id,
5861 const cs_var_cal_opt_t var_cal_opt,
5862 int inc,
5863 cs_real_3_t *restrict pvar,
5864 const cs_real_3_t *restrict pvara,
5865 const cs_real_3_t coefav[],
5866 const cs_real_33_t coefbv[],
5867 const cs_real_3_t cofafv[],
5868 const cs_real_33_t cofbfv[],
5869 const cs_real_t i_visc[],
5870 const cs_real_t b_visc[],
5871 cs_real_6_t *restrict viscel,
5872 const cs_real_2_t weighf[],
5873 const cs_real_t weighb[],
5874 cs_real_3_t *restrict rhs);
5875
5876 /*----------------------------------------------------------------------------*/
5877 /*!
5878 * \brief Add the explicit part of the diffusion terms with a symmetric tensor
5879 * diffusivity for a transport equation of a scalar field \f$ \varia \f$.
5880 *
5881 * More precisely, the right hand side \f$ Rhs \f$ is updated as
5882 * follows:
5883 * \f[
5884 * Rhs = Rhs - \sum_{\fij \in \Facei{\celli}} \left(
5885 * - \tens{\mu}_\fij \gradv_\fij \varia \cdot \vect{S}_\ij \right)
5886 * \f]
5887 *
5888 * Warning:
5889 * - \f$ Rhs \f$ has already been initialized before
5890 * calling cs_anisotropic_diffusion_scalar!
5891 * - mind the sign minus
5892 *
5893 * \param[in] idtvar indicator of the temporal scheme
5894 * \param[in] f_id index of the current variable
5895 * \param[in] var_cal_opt variable calculation options
5896 * \param[in] inc indicator
5897 * - 0 when solving an increment
5898 * - 1 otherwise
5899 * \param[in] pvar solved variable (current time step)
5900 * \param[in] pvara solved variable (previous time step)
5901 * \param[in] coefa boundary condition array for the variable
5902 * (explicit part)
5903 * \param[in] coefb boundary condition array for the variable
5904 * (implicit part)
5905 * \param[in] cofaf boundary condition array for the diffusion
5906 * of the variable (explicit part)
5907 * \param[in] cofbf boundary condition array for the diffusion
5908 * of the variable (implicit part)
5909 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
5910 * at interior faces for the r.h.s.
5911 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
5912 * at border faces for the r.h.s.
5913 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
5914 * \param[in] weighf internal face weight between cells i j in case
5915 * of tensor diffusion
5916 * \param[in] weighb boundary face weight for cells i in case
5917 * of tensor diffusion
5918 * \param[in,out] rhs right hand side \f$ \vect{Rhs} \f$
5919 */
5920 /*----------------------------------------------------------------------------*/
5921
5922 void
5923 cs_anisotropic_diffusion_tensor(int idtvar,
5924 int f_id,
5925 const cs_var_cal_opt_t var_cal_opt,
5926 int inc,
5927 cs_real_6_t *restrict pvar,
5928 const cs_real_6_t *restrict pvara,
5929 const cs_real_6_t coefa[],
5930 const cs_real_66_t coefb[],
5931 const cs_real_6_t cofaf[],
5932 const cs_real_66_t cofbf[],
5933 const cs_real_t i_visc[],
5934 const cs_real_t b_visc[],
5935 cs_real_6_t *restrict viscel,
5936 const cs_real_2_t weighf[],
5937 const cs_real_t weighb[],
5938 cs_real_6_t *restrict rhs);
5939
5940 /*----------------------------------------------------------------------------*/
5941 /*!
5942 * \brief Update the face mass flux with the face pressure (or pressure
5943 * increment, or pressure double increment) gradient.
5944 *
5945 * \f[
5946 * \dot{m}_\ij = \dot{m}_\ij
5947 * - \Delta t \grad_\fij \delta p \cdot \vect{S}_\ij
5948 * \f]
5949 *
5950 * \param[in] f_id field id (or -1)
5951 * \param[in] m pointer to mesh
5952 * \param[in] fvq pointer to finite volume quantities
5953 * \param[in] init indicator
5954 * - 1 initialize the mass flux to 0
5955 * - 0 otherwise
5956 * \param[in] inc indicator
5957 * - 0 when solving an increment
5958 * - 1 otherwise
5959 * \param[in] imrgra indicator
5960 * - 0 iterative gradient
5961 * - 1 least square gradient
5962 * \param[in] iccocg indicator
5963 * - 1 re-compute cocg matrix
5964 * (for iterative gradients)
5965 * - 0 otherwise
5966 * \param[in] nswrgp number of reconstruction sweeps for the
5967 * gradients
5968 * \param[in] imligp clipping gradient method
5969 * - < 0 no clipping
5970 * - = 0 thank to neighbooring gradients
5971 * - = 1 thank to the mean gradient
5972 * \param[in] iphydp hydrostatic pressure indicator
5973 * \param[in] iwgrp indicator
5974 * - 1 weight gradient by vicosity*porosity
5975 * - weighting determined by field options
5976 * \param[in] iwarnp verbosity
5977 * \param[in] epsrgp relative precision for the gradient
5978 * reconstruction
5979 * \param[in] climgp clipping coeffecient for the computation of
5980 * the gradient
5981 * \param[in] frcxt body force creating the hydrostatic pressure
5982 * \param[in] pvar solved variable (current time step)
5983 * \param[in] coefap boundary condition array for the variable
5984 * (explicit part)
5985 * \param[in] coefbp boundary condition array for the variable
5986 * (implicit part)
5987 * \param[in] cofafp boundary condition array for the diffusion
5988 * of the variable (explicit part)
5989 * \param[in] cofbfp boundary condition array for the diffusion
5990 * of the variable (implicit part)
5991 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
5992 * at interior faces for the r.h.s.
5993 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
5994 * at border faces for the r.h.s.
5995 * \param[in] visel viscosity by cell
5996 * \param[in,out] i_massflux mass flux at interior faces
5997 * \param[in,out] b_massflux mass flux at boundary faces
5998 */
5999 /*----------------------------------------------------------------------------*/
6000
6001 void
6002 cs_face_diffusion_potential(const int f_id,
6003 const cs_mesh_t *m,
6004 cs_mesh_quantities_t *fvq,
6005 int init,
6006 int inc,
6007 int imrgra,
6008 int iccocg,
6009 int nswrgp,
6010 int imligp,
6011 int iphydp,
6012 int iwgrp,
6013 int iwarnp,
6014 double epsrgp,
6015 double climgp,
6016 cs_real_3_t *restrict frcxt,
6017 cs_real_t *restrict pvar,
6018 const cs_real_t coefap[],
6019 const cs_real_t coefbp[],
6020 const cs_real_t cofafp[],
6021 const cs_real_t cofbfp[],
6022 const cs_real_t i_visc[],
6023 const cs_real_t b_visc[],
6024 cs_real_t *restrict visel,
6025 cs_real_t *restrict i_massflux,
6026 cs_real_t *restrict b_massflux);
6027
6028 /*----------------------------------------------------------------------------*/
6029 /*!
6030 * \brief Add the explicit part of the pressure gradient term to the mass flux
6031 * in case of anisotropic diffusion of the pressure field \f$ P \f$.
6032 *
6033 * More precisely, the mass flux side \f$ \dot{m}_\fij \f$ is updated as
6034 * follows:
6035 * \f[
6036 * \dot{m}_\fij = \dot{m}_\fij -
6037 * \left( \tens{\mu}_\fij \gradv_\fij P \cdot \vect{S}_\ij \right)
6038 * \f]
6039 *
6040 * \param[in] f_id field id (or -1)
6041 * \param[in] m pointer to mesh
6042 * \param[in] fvq pointer to finite volume quantities
6043 * \param[in] init indicator
6044 * - 1 initialize the mass flux to 0
6045 * - 0 otherwise
6046 * \param[in] inc indicator
6047 * - 0 when solving an increment
6048 * - 1 otherwise
6049 * \param[in] imrgra indicator
6050 * - 0 iterative gradient
6051 * - 1 least square gradient
6052 * \param[in] iccocg indicator
6053 * - 1 re-compute cocg matrix
6054 (for iterativ gradients)
6055 * - 0 otherwise
6056 * \param[in] nswrgp number of reconstruction sweeps for the
6057 * gradients
6058 * \param[in] imligp clipping gradient method
6059 * - < 0 no clipping
6060 * - = 0 thank to neighbooring gradients
6061 * - = 1 thank to the mean gradient
6062 * \param[in] ircflp indicator
6063 * - 1 flux reconstruction,
6064 * - 0 otherwise
6065 * \param[in] iphydp indicator
6066 * - 1 hydrostatic pressure taken into account
6067 * - 0 otherwise
6068 * \param[in] iwgrp indicator
6069 * - 1 weight gradient by vicosity*porosity
6070 * - weighting determined by field options
6071 * \param[in] iwarnp verbosity
6072 * \param[in] epsrgp relative precision for the gradient
6073 * reconstruction
6074 * \param[in] climgp clipping coeffecient for the computation of
6075 * the gradient
6076 * \param[in] frcxt body force creating the hydrostatic pressure
6077 * \param[in] pvar solved variable (pressure)
6078 * \param[in] coefap boundary condition array for the variable
6079 * (explicit part)
6080 * \param[in] coefbp boundary condition array for the variable
6081 * (implicit part)
6082 * \param[in] cofafp boundary condition array for the diffusion
6083 * of the variable (explicit part)
6084 * \param[in] cofbfp boundary condition array for the diffusion
6085 * of the variable (implicit part)
6086 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
6087 * at interior faces for the r.h.s.
6088 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
6089 * at border faces for the r.h.s.
6090 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
6091 * \param[in] weighf internal face weight between cells i j in case
6092 * of tensor diffusion
6093 * \param[in] weighb boundary face weight for cells i in case
6094 * of tensor diffusion
6095 * \param[in,out] i_massflux mass flux at interior faces
6096 * \param[in,out] b_massflux mass flux at boundary faces
6097 */
6098 /*----------------------------------------------------------------------------*/
6099
6100 void
6101 cs_face_anisotropic_diffusion_potential(const int f_id,
6102 const cs_mesh_t *m,
6103 cs_mesh_quantities_t *fvq,
6104 int init,
6105 int inc,
6106 int imrgra,
6107 int iccocg,
6108 int nswrgp,
6109 int imligp,
6110 int ircflp,
6111 int iphydp,
6112 int iwgrp,
6113 int iwarnp,
6114 double epsrgp,
6115 double climgp,
6116 cs_real_3_t *restrict frcxt,
6117 cs_real_t *restrict pvar,
6118 const cs_real_t coefap[],
6119 const cs_real_t coefbp[],
6120 const cs_real_t cofafp[],
6121 const cs_real_t cofbfp[],
6122 const cs_real_t i_visc[],
6123 const cs_real_t b_visc[],
6124 cs_real_6_t *restrict viscel,
6125 const cs_real_2_t weighf[],
6126 const cs_real_t weighb[],
6127 cs_real_t *restrict i_massflux,
6128 cs_real_t *restrict b_massflux);
6129
6130 /*----------------------------------------------------------------------------*/
6131 /*!
6132 * \brief Update the cell mass flux divergence with the face pressure (or
6133 * pressure increment, or pressure double increment) gradient.
6134 *
6135 * \f[
6136 * \dot{m}_\ij = \dot{m}_\ij
6137 * - \sum_j \Delta t \grad_\fij p \cdot \vect{S}_\ij
6138 * \f]
6139 *
6140 * \param[in] f_id field id (or -1)
6141 * \param[in] m pointer to mesh
6142 * \param[in] fvq pointer to finite volume quantities
6143 * \param[in] init indicator
6144 * - 1 initialize the mass flux to 0
6145 * - 0 otherwise
6146 * \param[in] inc indicator
6147 * - 0 when solving an increment
6148 * - 1 otherwise
6149 * \param[in] imrgra indicator
6150 * - 0 iterative gradient
6151 * - 1 least square gradient
6152 * \param[in] iccocg indicator
6153 * - 1 re-compute cocg matrix
6154 * (for iterative gradients)
6155 * - 0 otherwise
6156 * \param[in] nswrgp number of reconstruction sweeps for the
6157 * gradients
6158 * \param[in] imligp clipping gradient method
6159 * - < 0 no clipping
6160 * - = 0 thank to neighbooring gradients
6161 * - = 1 thank to the mean gradient
6162 * \param[in] iphydp hydrostatic pressure indicator
6163 * \param[in] iwgrp indicator
6164 * - 1 weight gradient by vicosity*porosity
6165 * - weighting determined by field options
6166 * \param[in] iwarnp verbosity
6167 * \param[in] epsrgp relative precision for the gradient
6168 * reconstruction
6169 * \param[in] climgp clipping coeffecient for the computation of
6170 * the gradient
6171 * \param[in] frcxt body force creating the hydrostatic pressure
6172 * \param[in] pvar solved variable (current time step)
6173 * \param[in] coefap boundary condition array for the variable
6174 * (explicit part)
6175 * \param[in] coefbp boundary condition array for the variable
6176 * (implicit part)
6177 * \param[in] cofafp boundary condition array for the diffusion
6178 * of the variable (explicit part)
6179 * \param[in] cofbfp boundary condition array for the diffusion
6180 * of the variable (implicit part)
6181 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
6182 * at interior faces for the r.h.s.
6183 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
6184 * at border faces for the r.h.s.
6185 * \param[in] visel viscosity by cell
6186 * \param[in,out] diverg mass flux divergence
6187 */
6188 /*----------------------------------------------------------------------------*/
6189
6190 void
6191 cs_diffusion_potential(const int f_id,
6192 const cs_mesh_t *m,
6193 cs_mesh_quantities_t *fvq,
6194 int init,
6195 int inc,
6196 int imrgra,
6197 int iccocg,
6198 int nswrgp,
6199 int imligp,
6200 int iphydp,
6201 int iwgrp,
6202 int iwarnp,
6203 double epsrgp,
6204 double climgp,
6205 cs_real_3_t *restrict frcxt,
6206 cs_real_t *restrict pvar,
6207 const cs_real_t coefap[],
6208 const cs_real_t coefbp[],
6209 const cs_real_t cofafp[],
6210 const cs_real_t cofbfp[],
6211 const cs_real_t i_visc[],
6212 const cs_real_t b_visc[],
6213 cs_real_t visel[],
6214 cs_real_t *restrict diverg);
6215
6216 /*----------------------------------------------------------------------------*/
6217 /*!
6218 * \brief Add the explicit part of the divergence of the mass flux due to the
6219 * pressure gradient (routine analog to cs_anisotropic_diffusion_scalar).
6220 *
6221 * More precisely, the divergence of the mass flux side
6222 * \f$ \sum_{\fij \in \Facei{\celli}} \dot{m}_\fij \f$ is updated as follows:
6223 * \f[
6224 * \sum_{\fij \in \Facei{\celli}} \dot{m}_\fij
6225 * = \sum_{\fij \in \Facei{\celli}} \dot{m}_\fij
6226 * - \sum_{\fij \in \Facei{\celli}}
6227 * \left( \tens{\mu}_\fij \gradv_\fij P \cdot \vect{S}_\ij \right)
6228 * \f]
6229 *
6230 * \param[in] f_id field id (or -1)
6231 * \param[in] m pointer to mesh
6232 * \param[in] fvq pointer to finite volume quantities
6233 * \param[in] init indicator
6234 * - 1 initialize the mass flux to 0
6235 * - 0 otherwise
6236 * \param[in] inc indicator
6237 * - 0 when solving an increment
6238 * - 1 otherwise
6239 * \param[in] imrgra indicator
6240 * - 0 iterative gradient
6241 * - 1 least square gradient
6242 * \param[in] iccocg indicator
6243 * - 1 re-compute cocg matrix
6244 (for iterativ gradients)
6245 * - 0 otherwise
6246 * \param[in] nswrgp number of reconstruction sweeps for the
6247 * gradients
6248 * \param[in] imligp clipping gradient method
6249 * - < 0 no clipping
6250 * - = 0 thank to neighbooring gradients
6251 * - = 1 thank to the mean gradient
6252 * \param[in] ircflp indicator
6253 * - 1 flux reconstruction,
6254 * - 0 otherwise
6255 * \param[in] iphydp indicator
6256 * - 1 hydrostatic pressure taken into account
6257 * - 0 otherwise
6258 * \param[in] iwgrp indicator
6259 * - 1 weight gradient by vicosity*porosity
6260 * - weighting determined by field options
6261 * \param[in] iwarnp verbosity
6262 * \param[in] epsrgp relative precision for the gradient
6263 * reconstruction
6264 * \param[in] climgp clipping coeffecient for the computation of
6265 * the gradient
6266 * \param[in] frcxt body force creating the hydrostatic pressure
6267 * \param[in] pvar solved variable (pressure)
6268 * \param[in] coefap boundary condition array for the variable
6269 * (explicit part)
6270 * \param[in] coefbp boundary condition array for the variable
6271 * (implicit part)
6272 * \param[in] cofafp boundary condition array for the diffusion
6273 * of the variable (explicit part)
6274 * \param[in] cofbfp boundary condition array for the diffusion
6275 * of the variable (implicit part)
6276 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
6277 * at interior faces for the r.h.s.
6278 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
6279 * at border faces for the r.h.s.
6280 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
6281 * \param[in] weighf internal face weight between cells i j in case
6282 * of tensor diffusion
6283 * \param[in] weighb boundary face weight for cells i in case
6284 * of tensor diffusion
6285 * \param[in,out] diverg divergence of the mass flux
6286 */
6287 /*----------------------------------------------------------------------------*/
6288
6289 void
6290 cs_anisotropic_diffusion_potential(const int f_id,
6291 const cs_mesh_t *m,
6292 cs_mesh_quantities_t *fvq,
6293 int init,
6294 int inc,
6295 int imrgra,
6296 int iccocg,
6297 int nswrgp,
6298 int imligp,
6299 int ircflp,
6300 int iphydp,
6301 int iwgrp,
6302 int iwarnp,
6303 double epsrgp,
6304 double climgp,
6305 cs_real_3_t *restrict frcxt,
6306 cs_real_t *restrict pvar,
6307 const cs_real_t coefap[],
6308 const cs_real_t coefbp[],
6309 const cs_real_t cofafp[],
6310 const cs_real_t cofbfp[],
6311 const cs_real_t i_visc[],
6312 const cs_real_t b_visc[],
6313 cs_real_6_t *restrict viscel,
6314 const cs_real_2_t weighf[],
6315 const cs_real_t weighb[],
6316 cs_real_t *restrict diverg);
6317
6318 /*----------------------------------------------------------------------------*/
6319
6320 END_C_DECLS
6321
6322 #endif /* __CS_CONVECTION_DIFFUSION_H__ */
6323