8.3
general documentation
cs_convection_diffusion_priv.h
Go to the documentation of this file.
1#ifndef __CS_CONVECTION_DIFFUSION_PRIV_H__
2#define __CS_CONVECTION_DIFFUSION_PRIV_H__
3
4/*============================================================================
5 * Private functions for convection-diffusion operators.
6 *============================================================================*/
7
8/*
9 This file is part of code_saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-2024 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_math.h"
37#include "cs_parameters.h" // for BC types
38
39#include "cs_convection_diffusion.h"
40
41/*=============================================================================
42 * Macro definitions
43 *============================================================================*/
44
45/*============================================================================
46 * Type definition
47 *============================================================================*/
48
49/*============================================================================
50 * Global variables
51 *============================================================================*/
52
53/*============================================================================
54 * Public inlined function
55 *============================================================================*/
56
57/*----------------------------------------------------------------------------
58 * Synchronize halos for scalar variables.
59 *
60 * parameters:
61 * m <-- pointer to associated mesh structure
62 * halo_type <-> halo type
63 * pvar <-> variable
64 *----------------------------------------------------------------------------*/
65
66inline static void
67cs_sync_scalar_halo(const cs_mesh_t *m,
68 cs_halo_type_t halo_type,
69 cs_real_t pvar[])
70{
71 if (m->halo != NULL)
72 cs_halo_sync_var(m->halo, halo_type, pvar);
73}
74
75/*----------------------------------------------------------------------------*/
86/*----------------------------------------------------------------------------*/
87
88CS_F_HOST_DEVICE inline static cs_real_t
89cs_nvd_scheme_scalar(const cs_nvd_type_t limiter,
90 const cs_real_t nvf_p_c,
91 const cs_real_t nvf_r_f,
92 const cs_real_t nvf_r_c)
93{
94 cs_real_t nvf_p_f;
95
96 cs_real_t beta_m, rfc, r1f, r1, r2, r3, b1, b2;
97
98 switch (limiter) {
99 case CS_NVD_GAMMA: /* Gamma scheme */
100 beta_m = 0.1; /* in [0.1, 0.5] */
101 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
102
103 if (nvf_p_c < beta_m) {
104 nvf_p_f = nvf_p_c*(1.+rfc*(1.-nvf_p_c)/beta_m);
105 } else {
106 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
107
108 nvf_p_f = r1f*nvf_p_c+rfc;
109 }
110
111 break;
112
113 case CS_NVD_SMART: /* SMART scheme */
114 if (nvf_p_c < (nvf_r_c/3.)) {
115 r1 = nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
116 r2 = nvf_r_c*(1.-nvf_r_c);
117
118 nvf_p_f = nvf_p_c*r1/r2;
119 } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
120 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
121 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
122
123 nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c + rfc);
124 } else {
125 nvf_p_f = 1.;
126 }
127
128 break;
129
130 case CS_NVD_CUBISTA: /* CUBISTA scheme */
131 if (nvf_p_c < (3.*nvf_r_c/4.)) {
132 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
133
134 nvf_p_f = nvf_r_f*(1.+rfc/3.)*nvf_p_c/nvf_r_c;
135 } else if (nvf_p_c <= (nvf_r_c*(1.+2.*(nvf_r_f-nvf_r_c))/(2.*nvf_r_f-nvf_r_c))) {
136 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
137 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
138
139 nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c+rfc);
140 } else {
141 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
142
143 nvf_p_f = 1.-.5*r1f*(1.-nvf_p_c);
144 }
145
146 break;
147
148 case CS_NVD_SUPERBEE: /* SuperBee scheme */
149 if (nvf_p_c < (nvf_r_c/(2.-nvf_r_c))) {
150 nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
151 } else if (nvf_p_c < nvf_r_c) {
152 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
153 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
154
155 nvf_p_f = r1f*nvf_p_c+rfc;
156 } else if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
157 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
158 } else {
159 nvf_p_f = 1.;
160 }
161
162 break;
163
164 case CS_NVD_MUSCL: /* MUSCL scheme */
165 if (nvf_p_c < (.5*nvf_r_c)) {
166 nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
167 } else if (nvf_p_c < (1.+nvf_r_c-nvf_r_f)) {
168 nvf_p_f = nvf_p_c+nvf_r_f-nvf_r_c;
169 } else {
170 nvf_p_f = 1.;
171 }
172
173 break;
174
175 case CS_NVD_MINMOD: /* MINMOD scheme */
176 if (nvf_p_c < nvf_r_c) {
177 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
178 } else {
179 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
180 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
181
182 nvf_p_f = r1f*nvf_p_c+rfc;
183 }
184
185 break;
186
187 case CS_NVD_CLAM: /* CLAM scheme */
188 r1 = nvf_r_c*nvf_r_c-nvf_r_f;
189 r2 = nvf_r_c*(nvf_r_c-1.);
190 r3 = nvf_r_f-nvf_r_c;
191
192 nvf_p_f = nvf_p_c*(r1+r3*nvf_p_c)/r2;
193
194 break;
195
196 case CS_NVD_STOIC: /* STOIC scheme */
197 b1 = (nvf_r_c-nvf_r_f)*nvf_r_c;
198 b2 = nvf_r_c+nvf_r_f+2.*nvf_r_f*nvf_r_f-4.*nvf_r_f*nvf_r_c;
199
200 if (nvf_p_c < (b1/b2)) {
201 r1 = -nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
202 r2 = nvf_r_c*(nvf_r_c-1.);
203
204 nvf_p_f = nvf_p_c*r1/r2;
205 } else if (nvf_p_c < nvf_r_c) {
206 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
207 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
208
209 nvf_p_f = r1f*nvf_p_c+rfc;
210 } else if (nvf_p_c < (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
211 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
212 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
213
214 nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
215 } else {
216 nvf_p_f = 1.;
217 }
218
219 break;
220
221 case CS_NVD_OSHER: /* OSHER scheme */
222 if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
223 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
224 } else {
225 nvf_p_f = 1.;
226 }
227
228 break;
229
230 case CS_NVD_WASEB: /* WASEB scheme */
231 r1 = nvf_r_c*nvf_r_f*(nvf_r_f-nvf_r_c);
232 r2 = 2.*nvf_r_c*(1.-nvf_r_c)-nvf_r_f*(1.-nvf_r_f);
233
234 if (nvf_p_c < (r1/r2)) {
235 nvf_p_f = 2.*nvf_p_c;
236 } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
237 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
238 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
239
240 nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
241 } else {
242 nvf_p_f = 1.;
243 }
244
245 break;
246
247 default: /* Upwinding */
248 nvf_p_f = nvf_p_c;
249 break;
250 }
251
252 return nvf_p_f;
253}
254
255/*----------------------------------------------------------------------------*/
271/*----------------------------------------------------------------------------*/
272
273CS_F_HOST_DEVICE inline static cs_real_t
274cs_nvd_vof_scheme_scalar(cs_nvd_type_t limiter,
275 const cs_nreal_t i_face_u_normal[3],
276 cs_real_t nvf_p_c,
277 cs_real_t nvf_r_f,
278 cs_real_t nvf_r_c,
279 const cs_real_t gradv_c[3],
280 const cs_real_t c_courant)
281{
282 assert(limiter >= CS_NVD_VOF_HRIC);
283
284 cs_real_t nvf_p_f;
285 cs_real_t blend, high_order, low_order, ratio;
286
287 /* Compute gradient angle indicator */
288 cs_real_t dotp = cs_math_fabs(cs_math_3_dot_product(gradv_c, i_face_u_normal));
289 cs_real_t sgrad = cs_math_3_norm(gradv_c);
290 cs_real_t denom = sgrad;
291
292 if (limiter == CS_NVD_VOF_HRIC) { /* M-HRIC scheme */
293 /* High order scheme : Bounded Downwind */
294 if (nvf_p_c <= .5) {
295 high_order = 2.*nvf_p_c;
296 } else {
297 high_order = 1.;
298 }
299
300 /* Low order scheme : MUSCL */
301 low_order = cs_nvd_scheme_scalar(CS_NVD_MUSCL,
302 nvf_p_c,
303 nvf_r_f,
304 nvf_r_c);
305
306 /* Compute the blending factor */
307 if (denom <= (cs_math_epzero*dotp)) {
308 blend = 1.;
309 } else {
310 ratio = dotp/denom;
311 blend = cs_math_fmin(1., pow(ratio, .5));
312 }
313
314 /* Blending */
315 nvf_p_f = blend*high_order + (1.-blend)*low_order;
316
317 /* Extra blending due to the cell Courant number */
318 if (c_courant < .7 && c_courant > .3) {
319 nvf_p_f = nvf_p_f + (nvf_p_f - low_order)*(.7 - c_courant )/.4;
320 } else if (c_courant >= .7) {
321 nvf_p_f = low_order;
322 }
323 } else if (limiter == CS_NVD_VOF_CICSAM) { /* M-CICSAM scheme */
324 /* High order scheme : HYPER-C + SUPERBEE */
325 if (c_courant <= .3) {
326 high_order = cs_math_fmin(1., nvf_p_c/(c_courant+cs_math_epzero));
327 } else if (c_courant <= .6) {
328 high_order = cs_math_fmin(1., nvf_p_c/.3);
329 } else if (c_courant <= .7) {
330 cs_real_t superbee = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
331 nvf_p_c,
332 nvf_r_f,
333 nvf_r_c);
334 high_order = 10.*( (.7-c_courant)*cs_math_fmin(1., nvf_p_c/.3)
335 + (c_courant-.6)*superbee);
336 }
337 else {
338 high_order = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
339 nvf_p_c,
340 nvf_r_f,
341 nvf_r_c);
342 }
343
344 /* Low order scheme : MUSCL */
345 low_order = cs_nvd_scheme_scalar(CS_NVD_MUSCL,
346 nvf_p_c,
347 nvf_r_f,
348 nvf_r_c);
349
350 /* Compute the blending factor */
351 if (denom <= (cs_math_epzero*dotp)) {
352 blend = 1.;
353 } else {
354 ratio = dotp/denom;
355 blend = cs_math_fmin(1., pow(ratio, 2.));
356 }
357
358 /* Blending */
359 nvf_p_f = blend*high_order + (1.-blend)*low_order;
360 } else { /* STACS scheme */
361 /* High order scheme : SUPERBEE */
362 high_order = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
363 nvf_p_c,
364 nvf_r_f,
365 nvf_r_c);
366
367 /* Low order scheme : STOIC */
368 low_order = cs_nvd_scheme_scalar(CS_NVD_STOIC,
369 nvf_p_c,
370 nvf_r_f,
371 nvf_r_c);
372
373 /* Compute the blending factor */
374 if (denom <= (cs_math_epzero*dotp)) {
375 blend = 1.;
376 } else {
377 ratio = dotp/denom;
378 blend = cs_math_fmin(1., pow(ratio, 4.));
379 }
380
381 /* Blending */
382 nvf_p_f = blend*high_order + (1.-blend)*low_order;
383 }
384
385 return nvf_p_f;
386}
387
388/*----------------------------------------------------------------------------*/
404/*----------------------------------------------------------------------------*/
405
406CS_F_HOST_DEVICE inline static void
407cs_slope_test(const cs_real_t pi,
408 const cs_real_t pj,
409 const cs_real_t distf,
410 const cs_nreal_t i_face_u_normal[3],
411 const cs_real_t gradi[3],
412 const cs_real_t gradj[3],
413 const cs_real_t grdpai[3],
414 const cs_real_t grdpaj[3],
415 const cs_real_t i_massflux,
416 cs_real_t *testij,
417 cs_real_t *tesqck)
418{
419 cs_real_t dcc, ddi, ddj;
420
421 /* Slope test */
422
423 *testij = grdpai[0]*grdpaj[0]
424 + grdpai[1]*grdpaj[1]
425 + grdpai[2]*grdpaj[2];
426
427 if (i_massflux > 0.) {
428 dcc = gradi[0]*i_face_u_normal[0]
429 + gradi[1]*i_face_u_normal[1]
430 + gradi[2]*i_face_u_normal[2];
431 ddi = grdpai[0]*i_face_u_normal[0]
432 + grdpai[1]*i_face_u_normal[1]
433 + grdpai[2]*i_face_u_normal[2];
434 ddj = (pj-pi)/distf;
435 }
436 else {
437 dcc = gradj[0]*i_face_u_normal[0]
438 + gradj[1]*i_face_u_normal[1]
439 + gradj[2]*i_face_u_normal[2];
440 ddi = (pj-pi)/distf;
441 ddj = grdpaj[0]*i_face_u_normal[0]
442 + grdpaj[1]*i_face_u_normal[1]
443 + grdpaj[2]*i_face_u_normal[2];
444 }
445
446 *tesqck = cs_math_sq(dcc) - cs_math_sq(ddi-ddj);
447}
448
449/*----------------------------------------------------------------------------*/
468/*----------------------------------------------------------------------------*/
469
470template <cs_lnum_t stride>
471CS_F_HOST_DEVICE inline static void
472cs_slope_test_strided(const cs_real_t pi[stride],
473 const cs_real_t pj[stride],
474 const cs_real_t distf,
475 const cs_nreal_t i_face_u_normal[3],
476 const cs_real_t gradi[stride][3],
477 const cs_real_t gradj[stride][3],
478 const cs_real_t gradsti[stride][3],
479 const cs_real_t gradstj[stride][3],
480 const cs_real_t i_massflux,
481 cs_real_t *testij,
482 cs_real_t *tesqck)
483{
484 cs_real_t dcc[stride], ddi[stride], ddj[stride];
485
486 *testij = 0.;
487 *tesqck = 0.;
488
489 /* Slope test */
490
491 for (int i = 0; i < stride; i++) {
492 *testij += gradsti[i][0]*gradstj[i][0]
493 + gradsti[i][1]*gradstj[i][1]
494 + gradsti[i][2]*gradstj[i][2];
495
496 if (i_massflux > 0.) {
497 dcc[i] = gradi[i][0]*i_face_u_normal[0]
498 + gradi[i][1]*i_face_u_normal[1]
499 + gradi[i][2]*i_face_u_normal[2];
500 ddi[i] = gradsti[i][0]*i_face_u_normal[0]
501 + gradsti[i][1]*i_face_u_normal[1]
502 + gradsti[i][2]*i_face_u_normal[2];
503 ddj[i] = (pj[i]-pi[i])/distf;
504 }
505 else {
506 dcc[i] = gradj[i][0]*i_face_u_normal[0]
507 + gradj[i][1]*i_face_u_normal[1]
508 + gradj[i][2]*i_face_u_normal[2];
509 ddi[i] = (pj[i]-pi[i])/distf;
510 ddj[i] = gradstj[i][0]*i_face_u_normal[0]
511 + gradstj[i][1]*i_face_u_normal[1]
512 + gradstj[i][2]*i_face_u_normal[2];
513 }
514
515 *tesqck += cs_math_sq(dcc[i]) - cs_math_sq(ddi[i]-ddj[i]);
516 }
517}
518
519/*----------------------------------------------------------------------------*/
535/*----------------------------------------------------------------------------*/
536
537CS_F_HOST_DEVICE inline static void
538cs_i_compute_quantities(const cs_real_t bldfrp,
539 const cs_rreal_3_t diipf,
540 const cs_rreal_3_t djjpf,
541 const cs_real_3_t gradi,
542 const cs_real_3_t gradj,
543 const cs_real_t pi,
544 const cs_real_t pj,
545 cs_real_t *recoi,
546 cs_real_t *recoj,
547 cs_real_t *pip,
548 cs_real_t *pjp)
549{
550 cs_real_t gradpf[3] = {0.5*(gradi[0] + gradj[0]),
551 0.5*(gradi[1] + gradj[1]),
552 0.5*(gradi[2] + gradj[2])};
553
554 /* reconstruction only if IRCFLP = 1 */
555 *recoi = bldfrp*(cs_math_3_dot_product(gradpf, diipf));
556 *recoj = bldfrp*(cs_math_3_dot_product(gradpf, djjpf));
557 *pip = pi + *recoi;
558 *pjp = pj + *recoj;
559}
560
561/*----------------------------------------------------------------------------*/
580/*----------------------------------------------------------------------------*/
581
582template <cs_lnum_t stride>
583CS_F_HOST_DEVICE inline static void
584cs_i_compute_quantities_strided(const cs_real_t bldfrp,
585 const cs_rreal_t diipf[3],
586 const cs_rreal_t djjpf[3],
587 const cs_real_t gradi[stride][3],
588 const cs_real_t gradj[stride][3],
589 const cs_real_t pi[stride],
590 const cs_real_t pj[stride],
591 cs_real_t recoi[stride],
592 cs_real_t recoj[stride],
593 cs_real_t pip[stride],
594 cs_real_t pjp[stride])
595{
596 cs_real_t dpvf[3];
597
598 /* x-y-z components, p = u, v, w */
599
600 for (int isou = 0; isou < stride; isou++) {
601
602 for (int jsou = 0; jsou < 3; jsou++)
603 dpvf[jsou] = 0.5*( gradi[isou][jsou]
604 + gradj[isou][jsou]);
605
606 /* reconstruction only if IRCFLP = 1 */
607
608 recoi[isou] = bldfrp*(cs_math_3_dot_product(dpvf, diipf));
609 recoj[isou] = bldfrp*(cs_math_3_dot_product(dpvf, djjpf));
610
611 pip[isou] = pi[isou] + recoi[isou];
612 pjp[isou] = pj[isou] + recoj[isou];
613
614 }
615}
616
617/*----------------------------------------------------------------------------*/
633/*----------------------------------------------------------------------------*/
634
635CS_F_HOST_DEVICE inline static void
636cs_i_relax_c_val(const double relaxp,
637 const cs_real_t pia,
638 const cs_real_t pja,
639 const cs_real_t recoi,
640 const cs_real_t recoj,
641 const cs_real_t pi,
642 const cs_real_t pj,
643 cs_real_t *pir,
644 cs_real_t *pjr,
645 cs_real_t *pipr,
646 cs_real_t *pjpr)
647{
648 *pir = pi/relaxp - (1.-relaxp)/relaxp * pia;
649 *pjr = pj/relaxp - (1.-relaxp)/relaxp * pja;
650
651 *pipr = *pir + recoi;
652 *pjpr = *pjr + recoj;
653}
654
655/*----------------------------------------------------------------------------*/
674/*----------------------------------------------------------------------------*/
675
676template <cs_lnum_t stride>
677CS_F_HOST_DEVICE inline static void
678cs_i_relax_c_val_strided(cs_real_t relaxp,
679 const cs_real_t pia[stride],
680 const cs_real_t pja[stride],
681 const cs_real_t recoi[stride],
682 const cs_real_t recoj[stride],
683 const cs_real_t pi[stride],
684 const cs_real_t pj[stride],
685 cs_real_t pir[stride],
686 cs_real_t pjr[stride],
687 cs_real_t pipr[stride],
688 cs_real_t pjpr[stride])
689{
690 for (int isou = 0; isou < stride; isou++) {
691 pir[isou] = pi[isou] /relaxp - (1.-relaxp)/relaxp * pia[isou];
692 pjr[isou] = pj[isou] /relaxp - (1.-relaxp)/relaxp * pja[isou];
693
694 pipr[isou] = pir[isou] + recoi[isou];
695 pjpr[isou] = pjr[isou] + recoj[isou];
696 }
697}
698
699/*----------------------------------------------------------------------------*/
706/*----------------------------------------------------------------------------*/
707
708CS_F_HOST_DEVICE inline static void
709cs_upwind_f_val(const cs_real_t p,
710 cs_real_t *pf)
711{
712 *pf = p;
713}
714
715/*----------------------------------------------------------------------------*/
725/*----------------------------------------------------------------------------*/
726
727template <cs_lnum_t stride>
728CS_F_HOST_DEVICE inline static void
729cs_upwind_f_val_strided(const cs_real_t p[stride],
730 cs_real_t pf[stride])
731{
732 for (int isou = 0; isou < stride; isou++)
733 pf[isou] = p[isou];
734}
735
736/*----------------------------------------------------------------------------*/
745/*----------------------------------------------------------------------------*/
746
747CS_F_HOST_DEVICE inline static void
748cs_centered_f_val(double pnd,
749 cs_real_t pip,
750 cs_real_t pjp,
751 cs_real_t *pf)
752{
753 *pf = pnd*pip + (1.-pnd)*pjp;
754}
755
756/*----------------------------------------------------------------------------*/
768/*----------------------------------------------------------------------------*/
769
770template <cs_lnum_t stride>
771CS_F_HOST_DEVICE inline static void
772cs_centered_f_val_strided(double pnd,
773 const cs_real_t pip[stride],
774 const cs_real_t pjp[stride],
775 cs_real_t pf[stride])
776{
777 for (int isou = 0; isou < stride; isou++)
778 pf[isou] = pnd*pip[isou] + (1.-pnd)*pjp[isou];
779}
780
781/*----------------------------------------------------------------------------*/
791/*----------------------------------------------------------------------------*/
792
793CS_F_HOST_DEVICE inline static void
794cs_solu_f_val(const cs_real_t cell_cen[3],
795 const cs_real_t i_face_cog[3],
796 const cs_real_t grad[3],
797 const cs_real_t p,
798 cs_real_t *pf)
799{
800 cs_real_t df[3];
801
802 df[0] = i_face_cog[0] - cell_cen[0];
803 df[1] = i_face_cog[1] - cell_cen[1];
804 df[2] = i_face_cog[2] - cell_cen[2];
805
806 *pf = p + cs_math_3_dot_product(df, grad);
807}
808
809/*----------------------------------------------------------------------------*/
822/*----------------------------------------------------------------------------*/
823
824template <cs_lnum_t stride>
825CS_F_HOST_DEVICE inline static void
826cs_solu_f_val_strided(const cs_real_t cell_cen[3],
827 const cs_real_t i_face_cog[3],
828 const cs_real_t grad[stride][3],
829 const cs_real_t p[stride],
830 cs_real_t pf[stride])
831{
832 cs_real_t df[3];
833
834 for (cs_lnum_t jsou = 0; jsou < 3; jsou++)
835 df[jsou] = i_face_cog[jsou] - cell_cen[jsou];
836
837 for (cs_lnum_t isou = 0; isou < stride; isou++) {
838 pf[isou] = p[isou] + df[0]*grad[isou][0]
839 + df[1]*grad[isou][1]
840 + df[2]*grad[isou][2];
841 }
842}
843
844/*----------------------------------------------------------------------------*/
854/*----------------------------------------------------------------------------*/
855
856CS_F_HOST_DEVICE inline static void
857cs_blend_f_val(const double blencp,
858 const cs_real_t p,
859 cs_real_t *pf)
860{
861 *pf = blencp * (*pf) + (1. - blencp) * p;
862}
863
864/*----------------------------------------------------------------------------*/
877/*----------------------------------------------------------------------------*/
878
879template <cs_lnum_t stride>
880CS_F_HOST_DEVICE inline static void
881cs_blend_f_val_strided(const double blencp,
882 const cs_real_t p[stride],
883 cs_real_t pf[stride])
884{
885 for (int isou = 0; isou < stride; isou++)
886 pf[isou] = blencp*(pf[isou])+(1.-blencp)*p[isou];
887}
888
889/*----------------------------------------------------------------------------*/
910/*----------------------------------------------------------------------------*/
911
912CS_F_HOST_DEVICE inline static void
913cs_i_conv_flux(const int iconvp,
914 const cs_real_t thetap,
915 const int imasac,
916 const cs_real_t pi,
917 const cs_real_t pj,
918 const cs_real_t pifri,
919 const cs_real_t pifrj,
920 const cs_real_t pjfri,
921 const cs_real_t pjfrj,
922 const cs_real_t i_massflux,
923 const cs_real_t xcppi,
924 const cs_real_t xcppj,
925 cs_real_2_t fluxij)
926{
927 cs_real_t flui, fluj;
928
929 flui = 0.5*(i_massflux + fabs(i_massflux));
930 fluj = 0.5*(i_massflux - fabs(i_massflux));
931
932 fluxij[0] += iconvp*xcppi*(thetap*(flui*pifri + fluj*pjfri) - imasac*i_massflux*pi);
933 fluxij[1] += iconvp*xcppj*(thetap*(flui*pifrj + fluj*pjfrj) - imasac*i_massflux*pj);
934}
935
936/*----------------------------------------------------------------------------*/
957/*----------------------------------------------------------------------------*/
958
959template <cs_lnum_t stride>
960CS_F_HOST_DEVICE inline static void
961cs_i_conv_flux_strided(int iconvp,
962 cs_real_t thetap,
963 int imasac,
964 const cs_real_t pi[stride],
965 const cs_real_t pj[stride],
966 const cs_real_t pifri[stride],
967 const cs_real_t pifrj[stride],
968 const cs_real_t pjfri[stride],
969 const cs_real_t pjfrj[stride],
970 cs_real_t i_massflux,
971 cs_real_t fluxi[stride],
972 cs_real_t fluxj[stride])
973{
974 cs_real_t flui, fluj;
975
976 flui = 0.5*(i_massflux + fabs(i_massflux));
977 fluj = 0.5*(i_massflux - fabs(i_massflux));
978
979 for (int isou = 0; isou < stride; isou++) {
980 fluxi[isou] += iconvp*( thetap*(flui*pifri[isou] + fluj*pjfri[isou])
981 - imasac*i_massflux*pi[isou]);
982 fluxj[isou] += iconvp*( thetap*(flui*pifrj[isou] + fluj*pjfrj[isou])
983 - imasac*i_massflux*pj[isou]);
984 }
985}
986
987/*----------------------------------------------------------------------------*/
1000/*----------------------------------------------------------------------------*/
1001
1002CS_F_HOST_DEVICE inline static void
1003cs_i_diff_flux(const int idiffp,
1004 const cs_real_t thetap,
1005 const cs_real_t pip,
1006 const cs_real_t pjp,
1007 const cs_real_t pipr,
1008 const cs_real_t pjpr,
1009 const cs_real_t i_visc,
1010 cs_real_t fluxij[2])
1011{
1012 fluxij[0] += idiffp*thetap*i_visc*(pipr -pjp);
1013 fluxij[1] += idiffp*thetap*i_visc*(pip -pjpr);
1014}
1015
1016/*----------------------------------------------------------------------------*/
1033/*----------------------------------------------------------------------------*/
1034
1035template <cs_lnum_t stride>
1036CS_F_HOST_DEVICE inline static void
1037cs_i_diff_flux_strided(int idiffp,
1038 cs_real_t thetap,
1039 const cs_real_t pip[stride],
1040 const cs_real_t pjp[stride],
1041 const cs_real_t pipr[stride],
1042 const cs_real_t pjpr[stride],
1043 const cs_real_t i_visc,
1044 cs_real_t fluxi[stride],
1045 cs_real_t fluxj[stride])
1046{
1047 for (int isou = 0; isou < stride; isou++) {
1048 fluxi[isou] += idiffp*thetap*i_visc*(pipr[isou] -pjp[isou]);
1049 fluxj[isou] += idiffp*thetap*i_visc*(pip[isou] -pjpr[isou]);
1050 }
1051}
1052
1053/*----------------------------------------------------------------------------*/
1077/*----------------------------------------------------------------------------*/
1078
1079CS_F_HOST_DEVICE inline static void
1080cs_i_cd_steady_upwind(const cs_real_t bldfrp,
1081 const cs_real_t relaxp,
1082 const cs_rreal_t diipf[3],
1083 const cs_rreal_t djjpf[3],
1084 const cs_real_t gradi[3],
1085 const cs_real_t gradj[3],
1086 const cs_real_t pi,
1087 const cs_real_t pj,
1088 const cs_real_t pia,
1089 const cs_real_t pja,
1090 cs_real_t *pifri,
1091 cs_real_t *pifrj,
1092 cs_real_t *pjfri,
1093 cs_real_t *pjfrj,
1094 cs_real_t *pip,
1095 cs_real_t *pjp,
1096 cs_real_t *pipr,
1097 cs_real_t *pjpr)
1098{
1099 cs_real_t pir, pjr;
1100 cs_real_t recoi, recoj;
1101
1102 cs_i_compute_quantities(bldfrp,
1103 diipf,
1104 djjpf,
1105 gradi,
1106 gradj,
1107 pi,
1108 pj,
1109 &recoi,
1110 &recoj,
1111 pip,
1112 pjp);
1113
1114 cs_i_relax_c_val(relaxp,
1115 pia,
1116 pja,
1117 recoi,
1118 recoj,
1119 pi,
1120 pj,
1121 &pir,
1122 &pjr,
1123 pipr,
1124 pjpr);
1125
1126 cs_upwind_f_val(pi,
1127 pifrj);
1128 cs_upwind_f_val(pir,
1129 pifri);
1130 cs_upwind_f_val(pj,
1131 pjfri);
1132 cs_upwind_f_val(pjr,
1133 pjfrj);
1134}
1135
1136/*----------------------------------------------------------------------------*/
1163/*----------------------------------------------------------------------------*/
1164
1165template <cs_lnum_t stride>
1166CS_F_HOST_DEVICE inline static void
1167cs_i_cd_steady_upwind_strided(const cs_real_t bldfrp,
1168 const cs_real_t relaxp,
1169 const cs_rreal_t diipf[3],
1170 const cs_rreal_t djjpf[3],
1171 const cs_real_t gradi[stride][3],
1172 const cs_real_t gradj[stride][3],
1173 const cs_real_t pi[stride],
1174 const cs_real_t pj[stride],
1175 const cs_real_t pia[stride],
1176 const cs_real_t pja[stride],
1177 cs_real_t pifri[stride],
1178 cs_real_t pifrj[stride],
1179 cs_real_t pjfri[stride],
1180 cs_real_t pjfrj[stride],
1181 cs_real_t pip[stride],
1182 cs_real_t pjp[stride],
1183 cs_real_t pipr[stride],
1184 cs_real_t pjpr[stride])
1185{
1186 cs_real_t pir[stride], pjr[stride];
1187 cs_real_t recoi[stride], recoj[stride];
1188
1189 cs_i_compute_quantities_strided<stride>(bldfrp,
1190 diipf,
1191 djjpf,
1192 gradi,
1193 gradj,
1194 pi,
1195 pj,
1196 recoi,
1197 recoj,
1198 pip,
1199 pjp);
1200
1201 cs_i_relax_c_val_strided<stride>(relaxp,
1202 pia,
1203 pja,
1204 recoi,
1205 recoj,
1206 pi,
1207 pj,
1208 pir,
1209 pjr,
1210 pipr,
1211 pjpr);
1212
1213 cs_upwind_f_val_strided<stride>(pi, pifrj);
1214 cs_upwind_f_val_strided<stride>(pir, pifri);
1215 cs_upwind_f_val_strided<stride>(pj, pjfri);
1216 cs_upwind_f_val_strided<stride>(pjr, pjfrj);
1217}
1218
1219/*----------------------------------------------------------------------------*/
1236/*----------------------------------------------------------------------------*/
1237
1238CS_F_HOST_DEVICE inline static void
1239cs_i_cd_unsteady_upwind(const cs_real_t bldfrp,
1240 const cs_rreal_t diipf[3],
1241 const cs_rreal_t djjpf[3],
1242 const cs_real_t gradi[3],
1243 const cs_real_t gradj[3],
1244 const cs_real_t pi,
1245 const cs_real_t pj,
1246 cs_real_t *pif,
1247 cs_real_t *pjf,
1248 cs_real_t *pip,
1249 cs_real_t *pjp)
1250{
1251 cs_real_t recoi, recoj;
1252
1253 cs_i_compute_quantities(bldfrp,
1254 diipf,
1255 djjpf,
1256 gradi,
1257 gradj,
1258 pi,
1259 pj,
1260 &recoi,
1261 &recoj,
1262 pip,
1263 pjp);
1264
1265 cs_upwind_f_val(pi, pif);
1266 cs_upwind_f_val(pj, pjf);
1267}
1268
1269/*----------------------------------------------------------------------------*/
1289/*----------------------------------------------------------------------------*/
1290
1291template <cs_lnum_t stride>
1292CS_F_HOST_DEVICE inline static void
1293cs_i_cd_unsteady_upwind_strided(const cs_real_t bldfrp,
1294 const cs_rreal_t diipf[3],
1295 const cs_rreal_t djjpf[3],
1296 const cs_real_t gradi[stride][3],
1297 const cs_real_t gradj[stride][3],
1298 const cs_real_t pi[stride],
1299 const cs_real_t pj[stride],
1300 cs_real_t pif[stride],
1301 cs_real_t pjf[stride],
1302 cs_real_t pip[stride],
1303 cs_real_t pjp[stride])
1304{
1305 cs_real_t recoi[stride], recoj[stride];
1306
1307 cs_i_compute_quantities_strided<stride>(bldfrp,
1308 diipf,
1309 djjpf,
1310 gradi,
1311 gradj,
1312 pi,
1313 pj,
1314 recoi,
1315 recoj,
1316 pip,
1317 pjp);
1318
1319 cs_upwind_f_val_strided<stride>(pi, pif);
1320 cs_upwind_f_val_strided<stride>(pj, pjf);
1321}
1322
1323/*----------------------------------------------------------------------------*/
1356/*----------------------------------------------------------------------------*/
1357
1358CS_F_HOST_DEVICE inline static void
1359cs_i_cd_steady(const cs_real_t bldfrp,
1360 const int ischcp,
1361 const double relaxp,
1362 const double blencp,
1363 const cs_real_t weight,
1364 const cs_real_t cell_ceni[3],
1365 const cs_real_t cell_cenj[3],
1366 const cs_real_t i_face_cog[3],
1367 const cs_rreal_t diipf[3],
1368 const cs_rreal_t djjpf[3],
1369 const cs_real_t gradi[3],
1370 const cs_real_t gradj[3],
1371 const cs_real_t gradupi[3],
1372 const cs_real_t gradupj[3],
1373 const cs_real_t pi,
1374 const cs_real_t pj,
1375 const cs_real_t pia,
1376 const cs_real_t pja,
1377 cs_real_t *pifri,
1378 cs_real_t *pifrj,
1379 cs_real_t *pjfri,
1380 cs_real_t *pjfrj,
1381 cs_real_t *pip,
1382 cs_real_t *pjp,
1383 cs_real_t *pipr,
1384 cs_real_t *pjpr)
1385{
1386 cs_real_t pir, pjr;
1387 cs_real_t recoi, recoj;
1388
1389 cs_i_compute_quantities(bldfrp,
1390 diipf,
1391 djjpf,
1392 gradi,
1393 gradj,
1394 pi,
1395 pj,
1396 &recoi,
1397 &recoj,
1398 pip,
1399 pjp);
1400
1401 cs_i_relax_c_val(relaxp,
1402 pia,
1403 pja,
1404 recoi,
1405 recoj,
1406 pi,
1407 pj,
1408 &pir,
1409 &pjr,
1410 pipr,
1411 pjpr);
1412
1413 if (ischcp == 1) {
1414
1415 /* Centered
1416 --------*/
1417
1418 cs_centered_f_val(weight,
1419 *pip,
1420 *pjpr,
1421 pifrj);
1422 cs_centered_f_val(weight,
1423 *pipr,
1424 *pjp,
1425 pifri);
1426 cs_centered_f_val(weight,
1427 *pipr,
1428 *pjp,
1429 pjfri);
1430 cs_centered_f_val(weight,
1431 *pip,
1432 *pjpr,
1433 pjfrj);
1434
1435 } else if (ischcp == 0) {
1436
1437 /* Original SOLU
1438 --------------*/
1439
1440 cs_solu_f_val(cell_ceni,
1441 i_face_cog,
1442 gradi,
1443 pi,
1444 pifrj);
1445 cs_solu_f_val(cell_ceni,
1446 i_face_cog,
1447 gradi,
1448 pir,
1449 pifri);
1450 cs_solu_f_val(cell_cenj,
1451 i_face_cog,
1452 gradj,
1453 pj,
1454 pjfri);
1455 cs_solu_f_val(cell_cenj,
1456 i_face_cog,
1457 gradj,
1458 pjr,
1459 pjfrj);
1460
1461 } else {
1462
1463 /* SOLU
1464 ----*/
1465
1466 cs_solu_f_val(cell_ceni,
1467 i_face_cog,
1468 gradupi,
1469 pi,
1470 pifrj);
1471 cs_solu_f_val(cell_ceni,
1472 i_face_cog,
1473 gradupi,
1474 pir,
1475 pifri);
1476 cs_solu_f_val(cell_cenj,
1477 i_face_cog,
1478 gradupj,
1479 pj,
1480 pjfri);
1481 cs_solu_f_val(cell_cenj,
1482 i_face_cog,
1483 gradupj,
1484 pjr,
1485 pjfrj);
1486
1487 }
1488
1489 /* Blending
1490 --------*/
1491
1492 cs_blend_f_val(blencp,
1493 pi,
1494 pifrj);
1495 cs_blend_f_val(blencp,
1496 pir,
1497 pifri);
1498 cs_blend_f_val(blencp,
1499 pj,
1500 pjfri);
1501 cs_blend_f_val(blencp,
1502 pjr,
1503 pjfrj);
1504}
1505
1506/*----------------------------------------------------------------------------*/
1540/*----------------------------------------------------------------------------*/
1541
1542template <cs_lnum_t stride>
1543CS_F_HOST_DEVICE inline static void
1544cs_i_cd_steady_strided(cs_real_t bldfrp,
1545 int ischcp,
1546 double relaxp,
1547 double blencp,
1548 cs_real_t weight,
1549 const cs_real_t cell_ceni[3],
1550 const cs_real_t cell_cenj[3],
1551 const cs_real_t i_face_cog[3],
1552 const cs_rreal_t diipf[3],
1553 const cs_rreal_t djjpf[3],
1554 const cs_real_t gradi[stride][3],
1555 const cs_real_t gradj[stride][3],
1556 const cs_real_t pi[stride],
1557 const cs_real_t pj[stride],
1558 const cs_real_t pia[stride],
1559 const cs_real_t pja[stride],
1560 cs_real_t pifri[stride],
1561 cs_real_t pifrj[stride],
1562 cs_real_t pjfri[stride],
1563 cs_real_t pjfrj[stride],
1564 cs_real_t pip[stride],
1565 cs_real_t pjp[stride],
1566 cs_real_t pipr[stride],
1567 cs_real_t pjpr[stride])
1568
1569{
1570 cs_real_t pir[stride], pjr[stride];
1571 cs_real_t recoi[stride], recoj[stride];
1572
1573 cs_i_compute_quantities_strided<stride>(bldfrp,
1574 diipf,
1575 djjpf,
1576 gradi,
1577 gradj,
1578 pi,
1579 pj,
1580 recoi,
1581 recoj,
1582 pip,
1583 pjp);
1584
1585 cs_i_relax_c_val_strided<stride>(relaxp,
1586 pia,
1587 pja,
1588 recoi,
1589 recoj,
1590 pi,
1591 pj,
1592 pir,
1593 pjr,
1594 pipr,
1595 pjpr);
1596
1597 if (ischcp == 1) {
1598
1599 /* Centered
1600 --------*/
1601
1602 cs_centered_f_val_strided<stride>(weight, pip, pjpr, pifrj);
1603 cs_centered_f_val_strided<stride>(weight, pipr, pjp, pifri);
1604 cs_centered_f_val_strided<stride>(weight, pipr, pjp, pjfri);
1605 cs_centered_f_val_strided<stride>(weight, pip, pjpr, pjfrj);
1606
1607 }
1608 else {
1609
1610 /* Second order
1611 ------------*/
1612
1613 cs_solu_f_val_strided<stride>(cell_ceni,
1614 i_face_cog,
1615 gradi,
1616 pi,
1617 pifrj);
1618 cs_solu_f_val_strided<stride>(cell_ceni,
1619 i_face_cog,
1620 gradi,
1621 pir,
1622 pifri);
1623 cs_solu_f_val_strided<stride>(cell_cenj,
1624 i_face_cog,
1625 gradj,
1626 pj,
1627 pjfri);
1628 cs_solu_f_val_strided<stride>(cell_cenj,
1629 i_face_cog,
1630 gradj,
1631 pjr,
1632 pjfrj);
1633
1634 }
1635
1636 /* Blending
1637 --------*/
1638
1639 cs_blend_f_val_strided<stride>(blencp, pi, pifrj);
1640 cs_blend_f_val_strided<stride>(blencp, pir, pifri);
1641 cs_blend_f_val_strided<stride>(blencp, pj, pjfri);
1642 cs_blend_f_val_strided<stride>(blencp, pjr, pjfrj);
1643}
1644
1645/*----------------------------------------------------------------------------*/
1675/*----------------------------------------------------------------------------*/
1676
1677CS_F_HOST_DEVICE inline static void
1678cs_i_cd_unsteady(const cs_real_t bldfrp,
1679 const int ischcp,
1680 const double blencp,
1681 const cs_real_t weight,
1682 const cs_real_3_t cell_ceni,
1683 const cs_real_3_t cell_cenj,
1684 const cs_real_3_t i_face_cog,
1685 const cs_real_t hybrid_blend_i,
1686 const cs_real_t hybrid_blend_j,
1687 const cs_rreal_3_t diipf,
1688 const cs_rreal_3_t djjpf,
1689 const cs_real_3_t gradi,
1690 const cs_real_3_t gradj,
1691 const cs_real_3_t gradupi,
1692 const cs_real_3_t gradupj,
1693 const cs_real_t pi,
1694 const cs_real_t pj,
1695 cs_real_t *pif,
1696 cs_real_t *pjf,
1697 cs_real_t *pip,
1698 cs_real_t *pjp)
1699{
1700 cs_real_t recoi, recoj;
1701
1702 cs_i_compute_quantities(bldfrp,
1703 diipf,
1704 djjpf,
1705 gradi,
1706 gradj,
1707 pi,
1708 pj,
1709 &recoi,
1710 &recoj,
1711 pip,
1712 pjp);
1713
1714
1715 if (ischcp == 1) {
1716
1717 /* Centered
1718 --------*/
1719
1720 cs_centered_f_val(weight,
1721 *pip,
1722 *pjp,
1723 pif);
1724 cs_centered_f_val(weight,
1725 *pip,
1726 *pjp,
1727 pjf);
1728
1729 } else if (ischcp == 0) {
1730
1731 /* Legacy SOLU
1732 -----------*/
1733
1734 cs_solu_f_val(cell_ceni,
1735 i_face_cog,
1736 gradi,
1737 pi,
1738 pif);
1739 cs_solu_f_val(cell_cenj,
1740 i_face_cog,
1741 gradj,
1742 pj,
1743 pjf);
1744
1745 } else if (ischcp == 3) {
1746
1747 /* Centered
1748 --------*/
1749
1750 cs_centered_f_val(weight,
1751 *pip,
1752 *pjp,
1753 pif);
1754 cs_centered_f_val(weight,
1755 *pip,
1756 *pjp,
1757 pjf);
1758
1759 /* Legacy SOLU
1760 -----------*/
1761 cs_real_t pif_up, pjf_up;
1762 cs_real_t hybrid_blend_interp;
1763
1764 cs_solu_f_val(cell_ceni,
1765 i_face_cog,
1766 gradi,
1767 pi,
1768 &pif_up);
1769 cs_solu_f_val(cell_cenj,
1770 i_face_cog,
1771 gradj,
1772 pj,
1773 &pjf_up);
1774
1775 hybrid_blend_interp = fmin(hybrid_blend_i,hybrid_blend_j);
1776 *pif = hybrid_blend_interp*(*pif) + (1. - hybrid_blend_interp)*pif_up;
1777 *pjf = hybrid_blend_interp*(*pjf) + (1. - hybrid_blend_interp)*pjf_up;
1778
1779 } else {
1780
1781 /* SOLU
1782 ----*/
1783
1784 cs_solu_f_val(cell_ceni,
1785 i_face_cog,
1786 gradupi,
1787 pi,
1788 pif);
1789 cs_solu_f_val(cell_cenj,
1790 i_face_cog,
1791 gradupj,
1792 pj,
1793 pjf);
1794
1795 }
1796
1797 /* Blending
1798 --------*/
1799
1800 cs_blend_f_val(blencp,
1801 pi,
1802 pif);
1803 cs_blend_f_val(blencp,
1804 pj,
1805 pjf);
1806}
1807
1808/*----------------------------------------------------------------------------*/
1837/*----------------------------------------------------------------------------*/
1838
1839template <cs_lnum_t stride>
1840CS_F_HOST_DEVICE inline static void
1841cs_i_cd_unsteady_strided(cs_real_t bldfrp,
1842 int ischcp,
1843 double blencp,
1844 cs_real_t weight,
1845 const cs_real_t cell_ceni[3],
1846 const cs_real_t cell_cenj[3],
1847 const cs_real_t i_face_cog[3],
1848 const cs_real_t hybrid_blend_i,
1849 const cs_real_t hybrid_blend_j,
1850 const cs_rreal_t diipf[3],
1851 const cs_rreal_t djjpf[3],
1852 const cs_real_t gradi[stride][3],
1853 const cs_real_t gradj[stride][3],
1854 const cs_real_t pi[stride],
1855 const cs_real_t pj[stride],
1856 cs_real_t pif[stride],
1857 cs_real_t pjf[stride],
1858 cs_real_t pip[stride],
1859 cs_real_t pjp[stride])
1860
1861{
1862 cs_real_t recoi[stride], recoj[stride];
1863
1864 cs_i_compute_quantities_strided<stride>(bldfrp,
1865 diipf,
1866 djjpf,
1867 gradi,
1868 gradj,
1869 pi,
1870 pj,
1871 recoi,
1872 recoj,
1873 pip,
1874 pjp);
1875
1876 if (ischcp == 1) {
1877
1878 /* Centered
1879 --------*/
1880
1881 cs_centered_f_val_strided<stride>(weight, pip, pjp, pif);
1882 cs_centered_f_val_strided<stride>(weight, pip, pjp, pjf);
1883
1884 }
1885 else if (ischcp == 3) {
1886
1887 /* Centered
1888 --------*/
1889
1890 cs_centered_f_val_strided<stride>(weight, pip, pjp, pif);
1891 cs_centered_f_val_strided<stride>(weight, pip, pjp, pjf);
1892
1893 /* SOLU
1894 -----*/
1895 cs_real_t pif_up[stride], pjf_up[stride];
1896 cs_real_t hybrid_blend_interp;
1897
1898 cs_solu_f_val_strided<stride>(cell_ceni,
1899 i_face_cog,
1900 gradi,
1901 pi,
1902 pif_up);
1903 cs_solu_f_val_strided<stride>(cell_cenj,
1904 i_face_cog,
1905 gradj,
1906 pj,
1907 pjf_up);
1908
1909 hybrid_blend_interp = fmin(hybrid_blend_i, hybrid_blend_j);
1910 for (int isou = 0; isou < stride; isou++) {
1911 pif[isou] = hybrid_blend_interp *pif[isou]
1912 + (1. - hybrid_blend_interp)*pif_up[isou];
1913 pjf[isou] = hybrid_blend_interp *pjf[isou]
1914 + (1. - hybrid_blend_interp)*pjf_up[isou];
1915 }
1916 }
1917 else {
1918
1919 /* Second order
1920 ------------*/
1921
1922 cs_solu_f_val_strided<stride>(cell_ceni,
1923 i_face_cog,
1924 gradi,
1925 pi,
1926 pif);
1927 cs_solu_f_val_strided<stride>(cell_cenj,
1928 i_face_cog,
1929 gradj,
1930 pj,
1931 pjf);
1932
1933 }
1934
1935 /* Blending
1936 --------*/
1937
1938 cs_blend_f_val_strided<stride>(blencp, pi, pif);
1939 cs_blend_f_val_strided<stride>(blencp, pj, pjf);
1940}
1941
1942/*----------------------------------------------------------------------------*/
1985/*----------------------------------------------------------------------------*/
1986
1987CS_F_HOST_DEVICE inline static void
1988cs_i_cd_steady_slope_test(bool *upwind_switch,
1989 const int iconvp,
1990 const cs_real_t bldfrp,
1991 const int ischcp,
1992 const double relaxp,
1993 const double blencp,
1994 const double blend_st,
1995 const cs_real_t weight,
1996 const cs_real_t i_dist,
1997 const cs_real_t cell_ceni[3],
1998 const cs_real_t cell_cenj[3],
1999 const cs_nreal_t i_face_u_normal[3],
2000 const cs_real_t i_face_cog[3],
2001 const cs_rreal_t diipf[3],
2002 const cs_rreal_t djjpf[3],
2003 const cs_real_t i_massflux,
2004 const cs_real_t gradi[3],
2005 const cs_real_t gradj[3],
2006 const cs_real_t gradupi[3],
2007 const cs_real_t gradupj[3],
2008 const cs_real_t gradsti[3],
2009 const cs_real_t gradstj[3],
2010 const cs_real_t pi,
2011 const cs_real_t pj,
2012 const cs_real_t pia,
2013 const cs_real_t pja,
2014 cs_real_t *pifri,
2015 cs_real_t *pifrj,
2016 cs_real_t *pjfri,
2017 cs_real_t *pjfrj,
2018 cs_real_t *pip,
2019 cs_real_t *pjp,
2020 cs_real_t *pipr,
2021 cs_real_t *pjpr)
2022{
2023 cs_real_t pir, pjr;
2024 cs_real_t recoi, recoj;
2025 cs_real_t testij, tesqck;
2026
2027 *upwind_switch = false;
2028
2029 cs_i_compute_quantities(bldfrp,
2030 diipf,
2031 djjpf,
2032 gradi,
2033 gradj,
2034 pi,
2035 pj,
2036 &recoi,
2037 &recoj,
2038 pip,
2039 pjp);
2040
2041 cs_i_relax_c_val(relaxp,
2042 pia,
2043 pja,
2044 recoi,
2045 recoj,
2046 pi,
2047 pj,
2048 &pir,
2049 &pjr,
2050 pipr,
2051 pjpr);
2052
2053 /* Convection slope test is needed only when iconv >0 */
2054 if (iconvp > 0) {
2055 cs_slope_test(pi,
2056 pj,
2057 i_dist,
2058 i_face_u_normal,
2059 gradi,
2060 gradj,
2061 gradsti,
2062 gradstj,
2063 i_massflux,
2064 &testij,
2065 &tesqck);
2066
2067 if (ischcp==1) {
2068
2069 /* Centered
2070 --------*/
2071
2072 cs_centered_f_val(weight,
2073 *pip,
2074 *pjpr,
2075 pifrj);
2076 cs_centered_f_val(weight,
2077 *pipr,
2078 *pjp,
2079 pifri);
2080 cs_centered_f_val(weight,
2081 *pipr,
2082 *pjp,
2083 pjfri);
2084 cs_centered_f_val(weight,
2085 *pip,
2086 *pjpr,
2087 pjfrj);
2088
2089 } else if (ischcp == 0) {
2090
2091 /* Second order
2092 ------------*/
2093
2094 cs_solu_f_val(cell_ceni,
2095 i_face_cog,
2096 gradi,
2097 pi,
2098 pifrj);
2099 cs_solu_f_val(cell_ceni,
2100 i_face_cog,
2101 gradi,
2102 pir,
2103 pifri);
2104 cs_solu_f_val(cell_cenj,
2105 i_face_cog,
2106 gradj,
2107 pj,
2108 pjfri);
2109 cs_solu_f_val(cell_cenj,
2110 i_face_cog,
2111 gradj,
2112 pjr,
2113 pjfrj);
2114
2115 } else {
2116
2117 /* SOLU
2118 -----*/
2119
2120 cs_solu_f_val(cell_ceni,
2121 i_face_cog,
2122 gradupi,
2123 pi,
2124 pifrj);
2125 cs_solu_f_val(cell_ceni,
2126 i_face_cog,
2127 gradupi,
2128 pir,
2129 pifri);
2130 cs_solu_f_val(cell_cenj,
2131 i_face_cog,
2132 gradupj,
2133 pj,
2134 pjfri);
2135 cs_solu_f_val(cell_cenj,
2136 i_face_cog,
2137 gradupj,
2138 pjr,
2139 pjfrj);
2140 }
2141
2142
2143 /* Slope test: percentage of upwind
2144 ----------------------------------*/
2145
2146 if (tesqck <= 0. || testij <= 0.) {
2147
2148 cs_blend_f_val(blend_st,
2149 pi,
2150 pifrj);
2151 cs_blend_f_val(blend_st,
2152 pir,
2153 pifri);
2154 cs_blend_f_val(blend_st,
2155 pj,
2156 pjfri);
2157 cs_blend_f_val(blend_st,
2158 pjr,
2159 pjfrj);
2160
2161 *upwind_switch = true;
2162
2163 }
2164
2165
2166 /* Blending
2167 --------*/
2168
2169 cs_blend_f_val(blencp,
2170 pi,
2171 pifrj);
2172 cs_blend_f_val(blencp,
2173 pir,
2174 pifri);
2175 cs_blend_f_val(blencp,
2176 pj,
2177 pjfri);
2178 cs_blend_f_val(blencp,
2179 pjr,
2180 pjfrj);
2181
2182 /* If iconv=0 p*fr* are useless */
2183 }
2184 else {
2185 cs_upwind_f_val(pi, pifrj);
2186 cs_upwind_f_val(pir, pifri);
2187 cs_upwind_f_val(pj, pjfri);
2188 cs_upwind_f_val(pjr, pjfrj);
2189 }
2190
2191}
2192
2193/*----------------------------------------------------------------------------*/
2237/*----------------------------------------------------------------------------*/
2238
2239template <cs_lnum_t stride>
2240CS_F_HOST_DEVICE inline static void
2241cs_i_cd_steady_slope_test_strided(bool *upwind_switch,
2242 int iconvp,
2243 cs_real_t bldfrp,
2244 int ischcp,
2245 double relaxp,
2246 double blencp,
2247 double blend_st,
2248 cs_real_t weight,
2249 cs_real_t i_dist,
2250 const cs_real_t cell_ceni[3],
2251 const cs_real_t cell_cenj[3],
2252 const cs_nreal_t i_face_u_normal[3],
2253 const cs_real_t i_face_cog[3],
2254 const cs_rreal_t diipf[3],
2255 const cs_rreal_t djjpf[3],
2256 cs_real_t i_massflux,
2257 const cs_real_t gradi[stride][3],
2258 const cs_real_t gradj[stride][3],
2259 const cs_real_t grdpai[stride][3],
2260 const cs_real_t grdpaj[stride][3],
2261 const cs_real_t pi[stride],
2262 const cs_real_t pj[stride],
2263 const cs_real_t pia[stride],
2264 const cs_real_t pja[stride],
2265 cs_real_t pifri[stride],
2266 cs_real_t pifrj[stride],
2267 cs_real_t pjfri[stride],
2268 cs_real_t pjfrj[stride],
2269 cs_real_t pip[stride],
2270 cs_real_t pjp[stride],
2271 cs_real_t pipr[stride],
2272 cs_real_t pjpr[stride])
2273{
2274 cs_real_t pir[stride], pjr[stride];
2275 cs_real_t recoi[stride], recoj[stride];
2276 cs_real_t testij, tesqck;
2277 int isou;
2278
2279 cs_i_compute_quantities_strided<stride>(bldfrp,
2280 diipf,
2281 djjpf,
2282 gradi,
2283 gradj,
2284 pi,
2285 pj,
2286 recoi,
2287 recoj,
2288 pip,
2289 pjp);
2290
2291 cs_i_relax_c_val_strided<stride>(relaxp,
2292 pia,
2293 pja,
2294 recoi,
2295 recoj,
2296 pi,
2297 pj,
2298 pir,
2299 pjr,
2300 pipr,
2301 pjpr);
2302
2303 /* Convection slope test is needed only when iconv >0 */
2304 if (iconvp > 0) {
2305 cs_slope_test_strided<stride>(pi,
2306 pj,
2307 i_dist,
2308 i_face_u_normal,
2309 gradi,
2310 gradj,
2311 grdpai,
2312 grdpaj,
2313 i_massflux,
2314 &testij,
2315 &tesqck);
2316
2317 for (isou = 0; isou < stride; isou++) {
2318 if (ischcp==1) {
2319
2320 /* Centered
2321 --------*/
2322
2323 cs_centered_f_val(weight, pip[isou], pjpr[isou], &pifrj[isou]);
2324 cs_centered_f_val(weight, pipr[isou], pjp[isou], &pifri[isou]);
2325 cs_centered_f_val(weight, pipr[isou], pjp[isou], &pjfri[isou]);
2326 cs_centered_f_val(weight, pip[isou], pjpr[isou], &pjfrj[isou]);
2327
2328 }
2329 else {
2330
2331 /* Second order
2332 ------------*/
2333
2334 cs_solu_f_val(cell_ceni, i_face_cog, gradi[isou], pi[isou],
2335 &pifrj[isou]);
2336 cs_solu_f_val(cell_ceni, i_face_cog, gradi[isou], pir[isou],
2337 &pifri[isou]);
2338 cs_solu_f_val(cell_cenj, i_face_cog, gradj[isou], pj[isou],
2339 &pjfri[isou]);
2340 cs_solu_f_val(cell_cenj, i_face_cog, gradj[isou], pjr[isou],
2341 &pjfrj[isou]);
2342
2343 }
2344
2345 }
2346
2347 /* Slope test: percentage of upwind
2348 ----------------------------------*/
2349
2350 if (tesqck <= 0. || testij <= 0.) {
2351
2352 cs_blend_f_val_strided<stride>(blend_st, pi, pifrj);
2353 cs_blend_f_val_strided<stride>(blend_st, pir, pifri);
2354 cs_blend_f_val_strided<stride>(blend_st, pj, pjfri);
2355 cs_blend_f_val_strided<stride>(blend_st, pjr, pjfrj);
2356
2357 *upwind_switch = true;
2358
2359 }
2360
2361 /* Blending
2362 --------*/
2363
2364 cs_blend_f_val_strided<stride>(blencp, pi, pifrj);
2365 cs_blend_f_val_strided<stride>(blencp, pir, pifri);
2366 cs_blend_f_val_strided<stride>(blencp, pj, pjfri);
2367 cs_blend_f_val_strided<stride>(blencp, pjr, pjfrj);
2368
2369 /* If iconv=0 p*fr* are useless */
2370 }
2371 else {
2372 for (isou = 0; isou < stride; isou++) {
2373 cs_upwind_f_val(pi[isou], &pifrj[isou]);
2374 cs_upwind_f_val(pir[isou], &pifri[isou]);
2375 cs_upwind_f_val(pj[isou], &pjfri[isou]);
2376 cs_upwind_f_val(pjr[isou], &pjfrj[isou]);
2377 }
2378 }
2379}
2380
2381/*----------------------------------------------------------------------------*/
2417/*----------------------------------------------------------------------------*/
2418
2419CS_F_HOST_DEVICE inline static void
2420cs_i_cd_unsteady_slope_test(bool *upwind_switch,
2421 const int iconvp,
2422 const cs_real_t bldfrp,
2423 const int ischcp,
2424 const double blencp,
2425 const double blend_st,
2426 const cs_real_t weight,
2427 const cs_real_t i_dist,
2428 const cs_real_3_t cell_ceni,
2429 const cs_real_3_t cell_cenj,
2430 const cs_nreal_3_t i_face_u_normal,
2431 const cs_real_3_t i_face_cog,
2432 const cs_rreal_3_t diipf,
2433 const cs_rreal_3_t djjpf,
2434 const cs_real_t i_massflux,
2435 const cs_real_3_t gradi,
2436 const cs_real_3_t gradj,
2437 const cs_real_3_t gradupi,
2438 const cs_real_3_t gradupj,
2439 const cs_real_3_t gradsti,
2440 const cs_real_3_t gradstj,
2441 const cs_real_t pi,
2442 const cs_real_t pj,
2443 cs_real_t *pif,
2444 cs_real_t *pjf,
2445 cs_real_t *pip,
2446 cs_real_t *pjp)
2447{
2448 CS_UNUSED(blend_st);
2449
2450 cs_real_t recoi, recoj;
2451 cs_real_t testij, tesqck;
2452
2453 *upwind_switch = false;
2454
2455 cs_i_compute_quantities(bldfrp,
2456 diipf,
2457 djjpf,
2458 gradi,
2459 gradj,
2460 pi,
2461 pj,
2462 &recoi,
2463 &recoj,
2464 pip,
2465 pjp);
2466
2467 /* Convection slope test is needed only when iconv >0 */
2468 if (iconvp > 0) {
2469 cs_slope_test(pi,
2470 pj,
2471 i_dist,
2472 i_face_u_normal,
2473 gradi,
2474 gradj,
2475 gradsti,
2476 gradstj,
2477 i_massflux,
2478 &testij,
2479 &tesqck);
2480
2481 if (ischcp==1) {
2482
2483 /* Centered
2484 --------*/
2485
2486 cs_centered_f_val(weight,
2487 *pip,
2488 *pjp,
2489 pif);
2490 cs_centered_f_val(weight,
2491 *pip,
2492 *pjp,
2493 pjf);
2494
2495 } else if (ischcp == 0) {
2496
2497 /* Original SOLU
2498 --------------*/
2499
2500 cs_solu_f_val(cell_ceni,
2501 i_face_cog,
2502 gradi,
2503 pi,
2504 pif);
2505 cs_solu_f_val(cell_cenj,
2506 i_face_cog,
2507 gradj,
2508 pj,
2509 pjf);
2510
2511 } else {
2512
2513 /* SOLU
2514 -----*/
2515
2516 cs_solu_f_val(cell_ceni,
2517 i_face_cog,
2518 gradupi,
2519 pi,
2520 pif);
2521 cs_solu_f_val(cell_cenj,
2522 i_face_cog,
2523 gradupj,
2524 pj,
2525 pjf);
2526
2527 }
2528
2529 /* Slope test: percentage of upwind
2530 -------------------------------- */
2531
2532 if (tesqck<=0. || testij<=0.) {
2533
2534 cs_blend_f_val(blend_st,
2535 pi,
2536 pif);
2537 cs_blend_f_val(blend_st,
2538 pj,
2539 pjf);
2540
2541 *upwind_switch = true;
2542
2543 }
2544
2545 /* Blending
2546 --------*/
2547
2548 cs_blend_f_val(blencp,
2549 pi,
2550 pif);
2551 cs_blend_f_val(blencp,
2552 pj,
2553 pjf);
2554
2555 /* If iconv=0 p*f are useless */
2556 } else {
2557 cs_upwind_f_val(pi,
2558 pif);
2559 cs_upwind_f_val(pj,
2560 pjf);
2561 }
2562
2563}
2564
2565/*----------------------------------------------------------------------------*/
2576/*----------------------------------------------------------------------------*/
2577
2578CS_F_HOST_DEVICE inline static void
2579cs_central_downwind_cells(const cs_lnum_t ii,
2580 const cs_lnum_t jj,
2581 const cs_real_t i_massflux,
2582 cs_lnum_t *ic,
2583 cs_lnum_t *id)
2584{
2585 if (i_massflux >= 0.) {
2586 *ic = ii;
2587 *id = jj;
2588 } else {
2589 *ic = jj;
2590 *id = ii;
2591 }
2592}
2593
2594/*----------------------------------------------------------------------------*/
2615/*----------------------------------------------------------------------------*/
2616
2617CS_F_HOST_DEVICE inline static void
2618cs_i_cd_unsteady_nvd(const cs_nvd_type_t limiter,
2619 const double beta,
2620 const cs_real_3_t cell_cen_c,
2621 const cs_real_3_t cell_cen_d,
2622 const cs_nreal_3_t i_face_u_normal,
2623 const cs_real_3_t i_face_cog,
2624 const cs_real_3_t gradv_c,
2625 const cs_real_t p_c,
2626 const cs_real_t p_d,
2627 const cs_real_t local_max_c,
2628 const cs_real_t local_min_c,
2629 const cs_real_t courant_c,
2630 cs_real_t *pif,
2631 cs_real_t *pjf)
2632{
2633 /* distance between face center and central cell center */
2634 /* Distance between face center and central cell center */
2635 cs_real_t dist_fc = cs_math_3_distance(cell_cen_c, i_face_cog);
2636
2637 /* Unit vector and distance between central and downwind cells centers */
2638 cs_real_t diff[3] = {cell_cen_d[0] - cell_cen_c[0],
2639 cell_cen_d[1] - cell_cen_c[1],
2640 cell_cen_d[2] - cell_cen_c[2]};
2641
2642 cs_real_t dist_dc = cs_math_3_norm(diff);
2643 cs_real_t invl = 1./dist_dc;
2644
2645 cs_real_t ndc[3] = {invl*diff[0],
2646 invl*diff[1],
2647 invl*diff[2]};
2648
2649 /* Place the upwind point on the line that joins
2650 the two cells on the upwind side and the same
2651 distance as that between the two cells */
2652 const cs_real_t dist_cu = dist_dc;
2653 const cs_real_t dist_du = dist_dc + dist_cu;
2654
2655 /* Compute the property on the upwind assuming a parabolic
2656 variation of the property between the two cells */
2657 const cs_real_t gradc = cs_math_3_dot_product(gradv_c, ndc);
2658
2659 const cs_real_t grad2c = ((p_d - p_c)/dist_dc - gradc)/dist_dc;
2660
2661 cs_real_t p_u = p_c + (grad2c*dist_cu - gradc)*dist_cu;
2662 p_u = cs_math_fmax(cs_math_fmin(p_u, local_max_c), local_min_c);
2663
2664 /* Compute the normalised distances */
2665 const cs_real_t nvf_r_f = (dist_fc+dist_cu)/dist_du;
2666 const cs_real_t nvf_r_c = dist_cu/dist_du;
2667
2668 /* Check for the bounds of the NVD diagram and compute the face
2669 property according to the selected NVD scheme */
2670 const cs_real_t _small = cs_math_epzero
2671 * (CS_ABS(p_u)+CS_ABS(p_c)+CS_ABS(p_d));
2672
2673 if (CS_ABS(p_d-p_u) <= _small) {
2674 *pif = p_c;
2675 *pjf = p_c;
2676 } else {
2677 const cs_real_t nvf_p_c = (p_c - p_u)/(p_d - p_u);
2678
2679 if (nvf_p_c <= 0. || nvf_p_c >= 1.) {
2680 *pif = p_c;
2681 *pjf = p_c;
2682 } else {
2683 cs_real_t nvf_p_f;
2684
2685 /* Highly compressive NVD scheme for VOF */
2686 if (limiter >= CS_NVD_VOF_HRIC) {
2687 nvf_p_f = cs_nvd_vof_scheme_scalar(limiter,
2688 i_face_u_normal,
2689 nvf_p_c,
2690 nvf_r_f,
2691 nvf_r_c,
2692 gradv_c,
2693 courant_c);
2694 } else { /* Regular NVD scheme */
2695 nvf_p_f = cs_nvd_scheme_scalar(limiter,
2696 nvf_p_c,
2697 nvf_r_f,
2698 nvf_r_c);
2699 }
2700
2701 *pif = p_u + nvf_p_f*(p_d - p_u);
2702 *pif = cs_math_fmax(cs_math_fmin(*pif, local_max_c), local_min_c);
2703
2704 cs_blend_f_val(beta,
2705 p_c,
2706 pif);
2707
2708 *pjf = *pif;
2709 }
2710 }
2711}
2712
2713/*----------------------------------------------------------------------------*/
2750/*----------------------------------------------------------------------------*/
2751
2752template <cs_lnum_t stride>
2753CS_F_HOST_DEVICE inline static void
2754cs_i_cd_unsteady_slope_test_strided(bool *upwind_switch,
2755 int iconvp,
2756 cs_real_t bldfrp,
2757 int ischcp,
2758 double blencp,
2759 double blend_st,
2760 cs_real_t weight,
2761 cs_real_t i_dist,
2762 const cs_real_t cell_ceni[3],
2763 const cs_real_t cell_cenj[3],
2764 const cs_nreal_t i_face_u_normal[3],
2765 const cs_real_t i_face_cog[3],
2766 const cs_rreal_t diipf[3],
2767 const cs_rreal_t djjpf[3],
2768 cs_real_t i_massflux,
2769 const cs_real_t gradi[stride][3],
2770 const cs_real_t gradj[stride][3],
2771 const cs_real_t grdpai[stride][3],
2772 const cs_real_t grdpaj[stride][3],
2773 const cs_real_t pi[stride],
2774 const cs_real_t pj[stride],
2775 cs_real_t pif[stride],
2776 cs_real_t pjf[stride],
2777 cs_real_t pip[stride],
2778 cs_real_t pjp[stride])
2779{
2780 cs_real_t recoi[stride], recoj[stride];
2781 cs_real_t testij, tesqck;
2782 int isou;
2783
2784 cs_i_compute_quantities_strided<stride>(bldfrp,
2785 diipf,
2786 djjpf,
2787 gradi,
2788 gradj,
2789 pi,
2790 pj,
2791 recoi,
2792 recoj,
2793 pip,
2794 pjp);
2795
2796 /* Convection slope test is needed only when iconv >0 */
2797 if (iconvp > 0) {
2798 cs_slope_test_strided<stride>(pi,
2799 pj,
2800 i_dist,
2801 i_face_u_normal,
2802 gradi,
2803 gradj,
2804 grdpai,
2805 grdpaj,
2806 i_massflux,
2807 &testij,
2808 &tesqck);
2809
2810 for (isou = 0; isou < stride; isou++) {
2811
2812 if (ischcp==1) {
2813
2814 /* Centered
2815 --------*/
2816
2817 cs_centered_f_val(weight,
2818 pip[isou],
2819 pjp[isou],
2820 &pif[isou]);
2821 cs_centered_f_val(weight,
2822 pip[isou],
2823 pjp[isou],
2824 &pjf[isou]);
2825
2826 }
2827 else {
2828
2829 /* Second order
2830 ------------*/
2831
2832 cs_solu_f_val(cell_ceni,
2833 i_face_cog,
2834 gradi[isou],
2835 pi[isou],
2836 &pif[isou]);
2837 cs_solu_f_val(cell_cenj,
2838 i_face_cog,
2839 gradj[isou],
2840 pj[isou],
2841 &pjf[isou]);
2842 }
2843
2844 }
2845
2846 /* Slope test activated: percentage of upwind */
2847 if (tesqck <= 0. || testij <= 0.) {
2848
2849 /* Upwind
2850 --------*/
2851
2852 cs_blend_f_val_strided<stride>(blend_st, pi, pif);
2853 cs_blend_f_val_strided<stride>(blend_st, pj, pjf);
2854
2855 *upwind_switch = true;
2856 }
2857
2858 /* Blending
2859 --------*/
2860
2861 cs_blend_f_val_strided<stride>(blencp, pi, pif);
2862 cs_blend_f_val_strided<stride>(blencp, pj, pjf);
2863
2864 /* If iconv=0 p*fr* are useless */
2865 }
2866 else {
2867
2868 for (isou = 0; isou < stride; isou++) {
2869 cs_upwind_f_val(pi[isou], &pif[isou]);
2870 cs_upwind_f_val(pj[isou], &pjf[isou]);
2871 }
2872 }
2873}
2874
2875/*----------------------------------------------------------------------------*/
2884/*----------------------------------------------------------------------------*/
2885
2886CS_F_HOST_DEVICE inline static void
2887cs_b_compute_quantities(const cs_rreal_t diipb[3],
2888 const cs_real_t gradi[3],
2889 const cs_real_t bldfrp,
2890 cs_real_t *recoi)
2891{
2892 *recoi = bldfrp * ( gradi[0]*diipb[0]
2893 + gradi[1]*diipb[1]
2894 + gradi[2]*diipb[2]);
2895}
2896
2897/*----------------------------------------------------------------------------*/
2909/*----------------------------------------------------------------------------*/
2910
2911template <cs_lnum_t stride>
2912CS_F_HOST_DEVICE inline static void
2913cs_b_compute_quantities_strided(const cs_rreal_t diipb[3],
2914 const cs_real_t gradi[stride][3],
2915 const cs_real_t bldfrp,
2916 cs_real_t recoi[stride])
2917{
2918 for (int isou = 0; isou < stride; isou++) {
2919 recoi[isou] = bldfrp * ( gradi[isou][0]*diipb[0]
2920 + gradi[isou][1]*diipb[1]
2921 + gradi[isou][2]*diipb[2]);
2922 }
2923}
2924
2925/*----------------------------------------------------------------------------*/
2936/*----------------------------------------------------------------------------*/
2937
2938CS_F_HOST_DEVICE inline static void
2939cs_b_relax_c_val(const double relaxp,
2940 const cs_real_t pi,
2941 const cs_real_t pia,
2942 const cs_real_t recoi,
2943 cs_real_t *pir,
2944 cs_real_t *pipr)
2945{
2946 *pir = pi/relaxp - (1.-relaxp)/relaxp*pia;
2947 *pipr = *pir + recoi;
2948}
2949
2950/*----------------------------------------------------------------------------*/
2964/*----------------------------------------------------------------------------*/
2965
2966template <cs_lnum_t stride>
2967CS_F_HOST_DEVICE inline static void
2968cs_b_relax_c_val_strided(const double relaxp,
2969 const cs_real_t pi[stride],
2970 const cs_real_t pia[stride],
2971 const cs_real_t recoi[stride],
2972 cs_real_t pir[stride],
2973 cs_real_t pipr[stride])
2974{
2975 for (int isou = 0; isou < stride; isou++) {
2976 pir[isou] = pi[isou]/relaxp - (1.-relaxp)/relaxp*pia[isou];
2977 pipr[isou] = pir[isou] + recoi[isou];
2978 }
2979}
2980
2981/*----------------------------------------------------------------------------*/
3005/*----------------------------------------------------------------------------*/
3006
3007CS_F_HOST_DEVICE inline static void
3008cs_b_imposed_conv_flux(int iconvp,
3009 cs_real_t thetap,
3010 int imasac,
3011 int inc,
3012 int bc_type,
3013 int icvfli,
3014 cs_real_t pi,
3015 cs_real_t pir,
3016 cs_real_t pipr,
3017 cs_real_t coefap,
3018 cs_real_t coefbp,
3019 cs_real_t coface,
3020 cs_real_t cofbce,
3021 cs_real_t b_massflux,
3022 cs_real_t xcpp,
3023 cs_real_t *flux)
3024{
3025 cs_real_t flui, fluj, pfac;
3026
3027 /* Computed convective flux */
3028
3029 if (icvfli == 0) {
3030
3031 /* Remove decentering for coupled faces */
3032 if (bc_type == CS_COUPLED_FD) {
3033 flui = 0.0;
3034 fluj = b_massflux;
3035 } else {
3036 flui = 0.5*(b_massflux +fabs(b_massflux));
3037 fluj = 0.5*(b_massflux -fabs(b_massflux));
3038 }
3039
3040 pfac = inc*coefap + coefbp*pipr;
3041 *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*pfac) -imasac*( b_massflux*pi));
3042
3043 /* Imposed convective flux */
3044
3045 } else {
3046
3047 pfac = inc*coface + cofbce*pipr;
3048 *flux += iconvp*xcpp*(-imasac*(b_massflux*pi) + thetap*(pfac));
3049
3050 }
3051}
3052
3053/*----------------------------------------------------------------------------*/
3078/*----------------------------------------------------------------------------*/
3079
3080template <cs_lnum_t stride>
3081CS_F_HOST_DEVICE inline static void
3082cs_b_imposed_conv_flux_strided(int iconvp,
3083 cs_real_t thetap,
3084 int imasac,
3085 int inc,
3086 int bc_type,
3087 int icvfli,
3088 const cs_real_t pi[stride],
3089 const cs_real_t pir[stride],
3090 const cs_real_t pipr[stride],
3091 const cs_real_t coefap[stride],
3092 const cs_real_t coefbp[stride][stride],
3093 const cs_real_t coface[stride],
3094 const cs_real_t cofbce[stride][stride],
3095 cs_real_t b_massflux,
3096 cs_real_t pfac[stride],
3097 cs_real_t flux[stride])
3098{
3099 cs_real_t flui, fluj;
3100
3101 /* Computed convective flux */
3102
3103 if (icvfli == 0) {
3104
3105 /* Remove decentering for coupled faces */
3106 if (bc_type == CS_COUPLED_FD) {
3107 flui = 0.0;
3108 fluj = b_massflux;
3109 }
3110 else {
3111 flui = 0.5*(b_massflux +fabs(b_massflux));
3112 fluj = 0.5*(b_massflux -fabs(b_massflux));
3113 }
3114 for (int isou = 0; isou < stride; isou++) {
3115 pfac[isou] = inc*coefap[isou];
3116 for (int jsou = 0; jsou < stride; jsou++) {
3117 pfac[isou] += coefbp[isou][jsou]*pipr[jsou];
3118 }
3119 flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac[isou])
3120 - imasac*b_massflux*pi[isou]);
3121 }
3122
3123 /* Imposed convective flux */
3124
3125 }
3126 else {
3127
3128 for (int isou = 0; isou < stride; isou++) {
3129 pfac[isou] = inc*coface[isou];
3130 for (int jsou = 0; jsou < stride; jsou++) {
3131 pfac[isou] += cofbce[isou][jsou]*pipr[jsou];
3132 }
3133 flux[isou] += iconvp*( thetap*pfac[isou]
3134 - imasac*b_massflux*pi[isou]);
3135 }
3136
3137 }
3138}
3139
3140/*----------------------------------------------------------------------------*/
3160/*----------------------------------------------------------------------------*/
3161
3162CS_F_HOST_DEVICE inline static void
3163cs_b_upwind_flux(const int iconvp,
3164 const cs_real_t thetap,
3165 const int imasac,
3166 const int inc,
3167 const int bc_type,
3168 const cs_real_t pi,
3169 const cs_real_t pir,
3170 const cs_real_t pipr,
3171 const cs_real_t coefap,
3172 const cs_real_t coefbp,
3173 const cs_real_t b_massflux,
3174 const cs_real_t xcpp,
3175 cs_real_t *flux)
3176{
3177 cs_real_t flui, fluj, pfac;
3178
3179 /* Remove decentering for coupled faces */
3180 if (bc_type == CS_COUPLED_FD) {
3181 flui = 0.0;
3182 fluj = b_massflux;
3183 } else {
3184 flui = 0.5*(b_massflux +fabs(b_massflux));
3185 fluj = 0.5*(b_massflux -fabs(b_massflux));
3186 }
3187
3188 pfac = inc*coefap + coefbp*pipr;
3189 *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*pfac) -imasac*( b_massflux*pi));
3190}
3191
3192/*----------------------------------------------------------------------------*/
3216/*----------------------------------------------------------------------------*/
3217
3218template <cs_lnum_t stride>
3219CS_F_HOST_DEVICE inline static void
3220cs_b_upwind_flux_strided(int iconvp,
3221 cs_real_t thetap,
3222 int imasac,
3223 int inc,
3224 int bc_type,
3225 const cs_real_t pi[stride],
3226 const cs_real_t pir[stride],
3227 const cs_real_t pipr[stride],
3228 const cs_real_t coefa[stride],
3229 const cs_real_t coefb[stride][stride],
3230 cs_real_t b_massflux,
3231 cs_real_t pfac[stride],
3232 cs_real_t flux[stride])
3233{
3234 cs_real_t flui, fluj;
3235
3236 /* Remove decentering for coupled faces */
3237 if (bc_type == CS_COUPLED_FD) {
3238 flui = 0.0;
3239 fluj = b_massflux;
3240 } else {
3241 flui = 0.5*(b_massflux +fabs(b_massflux));
3242 fluj = 0.5*(b_massflux -fabs(b_massflux));
3243 }
3244 for (int isou = 0; isou < stride; isou++) {
3245 pfac[isou] = inc*coefa[isou];
3246 for (int jsou = 0; jsou < stride; jsou++) {
3247 pfac[isou] += coefb[isou][jsou]*pipr[jsou];
3248 }
3249 flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac[isou])
3250 - imasac*b_massflux*pi[isou]);
3251 }
3252}
3253
3254/*----------------------------------------------------------------------------*/
3267/*----------------------------------------------------------------------------*/
3268
3269CS_F_HOST_DEVICE inline static void
3270cs_b_diff_flux(const int idiffp,
3271 const cs_real_t thetap,
3272 const int inc,
3273 const cs_real_t pipr,
3274 const cs_real_t cofafp,
3275 const cs_real_t cofbfp,
3276 const cs_real_t b_visc,
3277 cs_real_t *flux)
3278{
3279 cs_real_t pfacd = inc*cofafp + cofbfp*pipr;
3280 *flux += idiffp*thetap*b_visc*pfacd;
3281}
3282
3283/*----------------------------------------------------------------------------*/
3299/*----------------------------------------------------------------------------*/
3300
3301template <cs_lnum_t stride>
3302CS_F_HOST_DEVICE inline static void
3303cs_b_diff_flux_strided(int idiffp,
3304 cs_real_t thetap,
3305 int inc,
3306 const cs_real_t pipr[stride],
3307 const cs_real_t cofaf[stride],
3308 const cs_real_t cofbf[stride][stride],
3309 cs_real_t b_visc,
3310 cs_real_t flux[stride])
3311{
3312 cs_real_t pfacd;
3313 for (int isou = 0; isou < stride; isou++) {
3314 pfacd = inc*cofaf[isou];
3315 for (int jsou = 0; jsou < stride; jsou++) {
3316 pfacd += cofbf[isou][jsou]*pipr[jsou];
3317 }
3318 flux[isou] += idiffp*thetap*b_visc*pfacd;
3319 }
3320}
3321
3322/*----------------------------------------------------------------------------*/
3336/*----------------------------------------------------------------------------*/
3337
3338CS_F_HOST_DEVICE inline static void
3339cs_b_cd_steady(const cs_real_t bldfrp,
3340 const double relaxp,
3341 const cs_rreal_3_t diipb,
3342 const cs_real_3_t gradi,
3343 const cs_real_t pi,
3344 const cs_real_t pia,
3345 cs_real_t *pir,
3346 cs_real_t *pipr)
3347{
3348 cs_real_t recoi;
3349
3350 cs_b_compute_quantities(diipb,
3351 gradi,
3352 bldfrp,
3353 &recoi);
3354
3355 cs_b_relax_c_val(relaxp,
3356 pi,
3357 pia,
3358 recoi,
3359 pir,
3360 pipr);
3361}
3362
3363/*----------------------------------------------------------------------------*/
3380/*----------------------------------------------------------------------------*/
3381
3382template <cs_lnum_t stride>
3383CS_F_HOST_DEVICE inline static void
3384cs_b_cd_steady_strided(cs_real_t bldfrp,
3385 double relaxp,
3386 const cs_rreal_t diipb[3],
3387 const cs_real_t gradi[stride][3],
3388 const cs_real_t pi[stride],
3389 const cs_real_t pia[stride],
3390 cs_real_t pir[stride],
3391 cs_real_t pipr[stride])
3392{
3393 cs_real_t recoi[stride];
3394
3395 cs_b_compute_quantities_strided<stride>(diipb,
3396 gradi,
3397 bldfrp,
3398 recoi);
3399
3400 cs_b_relax_c_val_strided<stride>(relaxp,
3401 pi,
3402 pia,
3403 recoi,
3404 pir,
3405 pipr);
3406}
3407
3408/*----------------------------------------------------------------------------*/
3419/*----------------------------------------------------------------------------*/
3420
3421CS_F_HOST_DEVICE inline static void
3422cs_b_cd_unsteady(const cs_real_t bldfrp,
3423 const cs_rreal_t diipb[3],
3424 const cs_real_t gradi[3],
3425 const cs_real_t pi,
3426 cs_real_t *pip)
3427{
3428 cs_real_t recoi;
3429
3430 cs_b_compute_quantities(diipb,
3431 gradi,
3432 bldfrp,
3433 &recoi);
3434
3435 *pip = pi + recoi;
3436}
3437
3438/*----------------------------------------------------------------------------*/
3452/*----------------------------------------------------------------------------*/
3453
3454template <cs_lnum_t stride>
3455CS_F_HOST_DEVICE inline static void
3456cs_b_cd_unsteady_strided(cs_real_t bldfrp,
3457 const cs_rreal_t diipb[3],
3458 const cs_real_t gradi[stride][3],
3459 const cs_real_t pi[stride],
3460 cs_real_t pip[stride])
3461{
3462 cs_real_t recoi[stride];
3463
3464 cs_b_compute_quantities_strided<stride>(diipb,
3465 gradi,
3466 bldfrp,
3467 recoi);
3468
3469 for (int isou = 0; isou< stride; isou++)
3470 pip[isou] = pi[isou] + recoi[isou];
3471}
3472
3473/*----------------------------------------------------------------------------*/
3484/*----------------------------------------------------------------------------*/
3485
3486CS_F_HOST_DEVICE inline static void
3487cs_b_diff_flux_coupling(int idiffp,
3488 cs_real_t pi,
3489 cs_real_t pj,
3490 cs_real_t b_visc,
3491 cs_real_t *fluxi)
3492{
3493 *fluxi += idiffp*b_visc*(pi - pj);
3494}
3495
3496/*----------------------------------------------------------------------------*/
3510/*----------------------------------------------------------------------------*/
3511
3512template <cs_lnum_t stride>
3513CS_F_HOST_DEVICE inline static void
3514cs_b_diff_flux_coupling_strided(int idiffp,
3515 const cs_real_t pi[stride],
3516 const cs_real_t pj[stride],
3517 cs_real_t b_visc,
3518 cs_real_t fluxi[stride])
3519{
3520 for (int k = 0; k < stride; k++)
3521 fluxi[k] += idiffp*b_visc*(pi[k] - pj[k]);
3522}
3523
3524/*----------------------------------------------------------------------------*/
3525
3526#endif /* __CS_CONVECTION_DIFFUSION_PRIV_H__ */
cs_convection_diffusion.h
cs_nvd_type_t
cs_nvd_type_t
Definition: cs_convection_diffusion.h:61
CS_NVD_SUPERBEE
@ CS_NVD_SUPERBEE
Definition: cs_convection_diffusion.h:66
CS_NVD_SMART
@ CS_NVD_SMART
Definition: cs_convection_diffusion.h:64
CS_NVD_CUBISTA
@ CS_NVD_CUBISTA
Definition: cs_convection_diffusion.h:65
CS_NVD_STOIC
@ CS_NVD_STOIC
Definition: cs_convection_diffusion.h:70
CS_NVD_WASEB
@ CS_NVD_WASEB
Definition: cs_convection_diffusion.h:72
CS_NVD_CLAM
@ CS_NVD_CLAM
Definition: cs_convection_diffusion.h:69
CS_NVD_OSHER
@ CS_NVD_OSHER
Definition: cs_convection_diffusion.h:71
CS_NVD_VOF_CICSAM
@ CS_NVD_VOF_CICSAM
Definition: cs_convection_diffusion.h:74
CS_NVD_GAMMA
@ CS_NVD_GAMMA
Definition: cs_convection_diffusion.h:63
CS_NVD_MINMOD
@ CS_NVD_MINMOD
Definition: cs_convection_diffusion.h:68
CS_NVD_VOF_HRIC
@ CS_NVD_VOF_HRIC
Definition: cs_convection_diffusion.h:73
CS_NVD_MUSCL
@ CS_NVD_MUSCL
Definition: cs_convection_diffusion.h:67
cs_b_diff_flux_coupling
static CS_F_HOST_DEVICE 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)
Add diffusive flux to flux at an internal coupling face.
Definition: cs_convection_diffusion_priv.h:3487
cs_i_diff_flux_strided
static CS_F_HOST_DEVICE void cs_i_diff_flux_strided(int idiffp, cs_real_t thetap, const cs_real_t pip[stride], const cs_real_t pjp[stride], const cs_real_t pipr[stride], const cs_real_t pjpr[stride], const cs_real_t i_visc, cs_real_t fluxi[stride], cs_real_t fluxj[stride])
Add diffusive fluxes to fluxes at face ij.
Definition: cs_convection_diffusion_priv.h:1037
cs_nvd_vof_scheme_scalar
static CS_F_HOST_DEVICE cs_real_t cs_nvd_vof_scheme_scalar(cs_nvd_type_t limiter, const cs_nreal_t i_face_u_normal[3], cs_real_t nvf_p_c, cs_real_t nvf_r_f, cs_real_t nvf_r_c, const cs_real_t gradv_c[3], const cs_real_t c_courant)
Compute the normalised face scalar using the specified NVD scheme for the case of a Volume-of-Fluid (...
Definition: cs_convection_diffusion_priv.h:274
cs_centered_f_val_strided
static CS_F_HOST_DEVICE void cs_centered_f_val_strided(double pnd, const cs_real_t pip[stride], const cs_real_t pjp[stride], cs_real_t pf[stride])
Prepare value at face ij by using a centered scheme.
Definition: cs_convection_diffusion_priv.h:772
cs_i_cd_unsteady_slope_test_strided
static CS_F_HOST_DEVICE void cs_i_cd_unsteady_slope_test_strided(bool *upwind_switch, int iconvp, cs_real_t bldfrp, int ischcp, double blencp, double blend_st, cs_real_t weight, cs_real_t i_dist, const cs_real_t cell_ceni[3], const cs_real_t cell_cenj[3], const cs_nreal_t i_face_u_normal[3], const cs_real_t i_face_cog[3], const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], cs_real_t i_massflux, const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t grdpai[stride][3], const cs_real_t grdpaj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], cs_real_t pif[stride], cs_real_t pjf[stride], cs_real_t pip[stride], cs_real_t pjp[stride])
Handle preparation of internal face values for the fluxes computation in case of a unsteady algorithm...
Definition: cs_convection_diffusion_priv.h:2754
cs_i_cd_unsteady
static CS_F_HOST_DEVICE 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_rreal_3_t diipf, const cs_rreal_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)
Handle preparation of internal face values for the fluxes computation in case of a unsteady algorithm...
Definition: cs_convection_diffusion_priv.h:1678
cs_b_diff_flux_strided
static CS_F_HOST_DEVICE void cs_b_diff_flux_strided(int idiffp, cs_real_t thetap, int inc, const cs_real_t pipr[stride], const cs_real_t cofaf[stride], const cs_real_t cofbf[stride][stride], cs_real_t b_visc, cs_real_t flux[stride])
Add diffusive flux to flux at boundary face.
Definition: cs_convection_diffusion_priv.h:3303
cs_i_relax_c_val_strided
static CS_F_HOST_DEVICE void cs_i_relax_c_val_strided(cs_real_t relaxp, const cs_real_t pia[stride], const cs_real_t pja[stride], const cs_real_t recoi[stride], const cs_real_t recoj[stride], const cs_real_t pi[stride], const cs_real_t pj[stride], cs_real_t pir[stride], cs_real_t pjr[stride], cs_real_t pipr[stride], cs_real_t pjpr[stride])
Compute relaxed values at cell i and j.
Definition: cs_convection_diffusion_priv.h:678
cs_b_relax_c_val_strided
static CS_F_HOST_DEVICE void cs_b_relax_c_val_strided(const double relaxp, const cs_real_t pi[stride], const cs_real_t pia[stride], const cs_real_t recoi[stride], cs_real_t pir[stride], cs_real_t pipr[stride])
Compute relaxed values at boundary cell i.
Definition: cs_convection_diffusion_priv.h:2968
cs_i_cd_unsteady_upwind_strided
static CS_F_HOST_DEVICE void cs_i_cd_unsteady_upwind_strided(const cs_real_t bldfrp, const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], cs_real_t pif[stride], cs_real_t pjf[stride], cs_real_t pip[stride], cs_real_t pjp[stride])
Handle preparation of internal face values for the fluxes computation in case of an unsteady algorith...
Definition: cs_convection_diffusion_priv.h:1293
cs_b_upwind_flux
static CS_F_HOST_DEVICE 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)
Add convective flux (substracting the mass accumulation from it) to flux at boundary face....
Definition: cs_convection_diffusion_priv.h:3163
cs_solu_f_val
static CS_F_HOST_DEVICE void cs_solu_f_val(const cs_real_t cell_cen[3], const cs_real_t i_face_cog[3], const cs_real_t grad[3], const cs_real_t p, cs_real_t *pf)
Prepare value at face ij by using a Second Order Linear Upwind scheme.
Definition: cs_convection_diffusion_priv.h:794
cs_sync_scalar_halo
static void cs_sync_scalar_halo(const cs_mesh_t *m, cs_halo_type_t halo_type, cs_real_t pvar[])
Definition: cs_convection_diffusion_priv.h:67
cs_b_cd_steady_strided
static CS_F_HOST_DEVICE void cs_b_cd_steady_strided(cs_real_t bldfrp, double relaxp, const cs_rreal_t diipb[3], const cs_real_t gradi[stride][3], const cs_real_t pi[stride], const cs_real_t pia[stride], cs_real_t pir[stride], cs_real_t pipr[stride])
Handle preparation of boundary face values for the flux computation in case of a steady algorithm.
Definition: cs_convection_diffusion_priv.h:3384
cs_i_relax_c_val
static CS_F_HOST_DEVICE 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)
Compute relaxed values at cell i and j.
Definition: cs_convection_diffusion_priv.h:636
cs_b_upwind_flux_strided
static CS_F_HOST_DEVICE void cs_b_upwind_flux_strided(int iconvp, cs_real_t thetap, int imasac, int inc, int bc_type, const cs_real_t pi[stride], const cs_real_t pir[stride], const cs_real_t pipr[stride], const cs_real_t coefa[stride], const cs_real_t coefb[stride][stride], cs_real_t b_massflux, cs_real_t pfac[stride], cs_real_t flux[stride])
Add convective flux (substracting the mass accumulation from it) to flux at boundary face....
Definition: cs_convection_diffusion_priv.h:3220
cs_i_cd_steady_upwind_strided
static CS_F_HOST_DEVICE void cs_i_cd_steady_upwind_strided(const cs_real_t bldfrp, const cs_real_t relaxp, const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], const cs_real_t pia[stride], const cs_real_t pja[stride], cs_real_t pifri[stride], cs_real_t pifrj[stride], cs_real_t pjfri[stride], cs_real_t pjfrj[stride], cs_real_t pip[stride], cs_real_t pjp[stride], cs_real_t pipr[stride], cs_real_t pjpr[stride])
Handle preparation of internal face values for the fluxes computation in case of a steady algorithm a...
Definition: cs_convection_diffusion_priv.h:1167
cs_i_cd_unsteady_slope_test
static CS_F_HOST_DEVICE 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_3_t cell_ceni, const cs_real_3_t cell_cenj, const cs_nreal_3_t i_face_u_normal, const cs_real_3_t i_face_cog, const cs_rreal_3_t diipf, const cs_rreal_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)
Handle preparation of internal face values for the fluxes computation in case of a steady algorithm a...
Definition: cs_convection_diffusion_priv.h:2420
cs_i_cd_unsteady_nvd
static CS_F_HOST_DEVICE 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_nreal_3_t i_face_u_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)
Handle preparation of internal face values for the convection flux computation in case of an unsteady...
Definition: cs_convection_diffusion_priv.h:2618
cs_b_diff_flux_coupling_strided
static CS_F_HOST_DEVICE void cs_b_diff_flux_coupling_strided(int idiffp, const cs_real_t pi[stride], const cs_real_t pj[stride], cs_real_t b_visc, cs_real_t fluxi[stride])
Add diffusive flux to flux at an internal coupling face for a vector.
Definition: cs_convection_diffusion_priv.h:3514
cs_upwind_f_val_strided
static CS_F_HOST_DEVICE void cs_upwind_f_val_strided(const cs_real_t p[stride], cs_real_t pf[stride])
Prepare value at face ij by using an upwind scheme.
Definition: cs_convection_diffusion_priv.h:729
cs_b_cd_steady
static CS_F_HOST_DEVICE void cs_b_cd_steady(const cs_real_t bldfrp, const double relaxp, const cs_rreal_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)
Handle preparation of boundary face values for the flux computation in case of a steady algorithm.
Definition: cs_convection_diffusion_priv.h:3339
cs_b_cd_unsteady_strided
static CS_F_HOST_DEVICE void cs_b_cd_unsteady_strided(cs_real_t bldfrp, const cs_rreal_t diipb[3], const cs_real_t gradi[stride][3], const cs_real_t pi[stride], cs_real_t pip[stride])
Handle preparation of boundary face values for the flux computation in case of a steady algorithm.
Definition: cs_convection_diffusion_priv.h:3456
cs_b_cd_unsteady
static CS_F_HOST_DEVICE void cs_b_cd_unsteady(const cs_real_t bldfrp, const cs_rreal_t diipb[3], const cs_real_t gradi[3], const cs_real_t pi, cs_real_t *pip)
Handle preparation of boundary face values for the flux computation in case of an unsteady algorithm.
Definition: cs_convection_diffusion_priv.h:3422
cs_b_relax_c_val
static CS_F_HOST_DEVICE 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)
Compute relaxed values at boundary cell i.
Definition: cs_convection_diffusion_priv.h:2939
cs_b_imposed_conv_flux_strided
static CS_F_HOST_DEVICE void cs_b_imposed_conv_flux_strided(int iconvp, cs_real_t thetap, int imasac, int inc, int bc_type, int icvfli, const cs_real_t pi[stride], const cs_real_t pir[stride], const cs_real_t pipr[stride], const cs_real_t coefap[stride], const cs_real_t coefbp[stride][stride], const cs_real_t coface[stride], const cs_real_t cofbce[stride][stride], cs_real_t b_massflux, cs_real_t pfac[stride], cs_real_t flux[stride])
Add convective flux (substracting the mass accumulation from it) to flux at boundary face....
Definition: cs_convection_diffusion_priv.h:3082
cs_i_compute_quantities
static CS_F_HOST_DEVICE void cs_i_compute_quantities(const cs_real_t bldfrp, const cs_rreal_3_t diipf, const cs_rreal_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)
Reconstruct values in I' and J'.
Definition: cs_convection_diffusion_priv.h:538
cs_slope_test
static CS_F_HOST_DEVICE void cs_slope_test(const cs_real_t pi, const cs_real_t pj, const cs_real_t distf, const cs_nreal_t i_face_u_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)
Compute slope test criteria at internal face between cell i and j.
Definition: cs_convection_diffusion_priv.h:407
cs_centered_f_val
static CS_F_HOST_DEVICE void cs_centered_f_val(double pnd, cs_real_t pip, cs_real_t pjp, cs_real_t *pf)
Prepare value at face ij by using a centered scheme.
Definition: cs_convection_diffusion_priv.h:748
cs_i_cd_steady_upwind
static CS_F_HOST_DEVICE void cs_i_cd_steady_upwind(const cs_real_t bldfrp, const cs_real_t relaxp, const cs_rreal_t diipf[3], const cs_rreal_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)
Handle preparation of internal face values for the fluxes computation in case of a steady algorithm a...
Definition: cs_convection_diffusion_priv.h:1080
cs_i_cd_steady_slope_test_strided
static CS_F_HOST_DEVICE void cs_i_cd_steady_slope_test_strided(bool *upwind_switch, int iconvp, cs_real_t bldfrp, int ischcp, double relaxp, double blencp, double blend_st, cs_real_t weight, cs_real_t i_dist, const cs_real_t cell_ceni[3], const cs_real_t cell_cenj[3], const cs_nreal_t i_face_u_normal[3], const cs_real_t i_face_cog[3], const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], cs_real_t i_massflux, const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t grdpai[stride][3], const cs_real_t grdpaj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], const cs_real_t pia[stride], const cs_real_t pja[stride], cs_real_t pifri[stride], cs_real_t pifrj[stride], cs_real_t pjfri[stride], cs_real_t pjfrj[stride], cs_real_t pip[stride], cs_real_t pjp[stride], cs_real_t pipr[stride], cs_real_t pjpr[stride])
Handle preparation of internal face values for the fluxes computation in case of a steady algorithm a...
Definition: cs_convection_diffusion_priv.h:2241
cs_i_cd_unsteady_upwind
static CS_F_HOST_DEVICE void cs_i_cd_unsteady_upwind(const cs_real_t bldfrp, const cs_rreal_t diipf[3], const cs_rreal_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)
Handle preparation of internal face values for the fluxes computation in case of an unsteady algorith...
Definition: cs_convection_diffusion_priv.h:1239
cs_blend_f_val_strided
static CS_F_HOST_DEVICE void cs_blend_f_val_strided(const double blencp, const cs_real_t p[stride], cs_real_t pf[stride])
Blend face values for a centered or SOLU scheme with face values for an upwind scheme.
Definition: cs_convection_diffusion_priv.h:881
cs_i_conv_flux_strided
static CS_F_HOST_DEVICE void cs_i_conv_flux_strided(int iconvp, cs_real_t thetap, int imasac, const cs_real_t pi[stride], const cs_real_t pj[stride], const cs_real_t pifri[stride], const cs_real_t pifrj[stride], const cs_real_t pjfri[stride], const cs_real_t pjfrj[stride], cs_real_t i_massflux, cs_real_t fluxi[stride], cs_real_t fluxj[stride])
Add convective fluxes (substracting the mass accumulation from them) to fluxes at face ij.
Definition: cs_convection_diffusion_priv.h:961
cs_upwind_f_val
static CS_F_HOST_DEVICE void cs_upwind_f_val(const cs_real_t p, cs_real_t *pf)
Prepare value at face ij by using an upwind scheme.
Definition: cs_convection_diffusion_priv.h:709
cs_b_diff_flux
static CS_F_HOST_DEVICE 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)
Add diffusive flux to flux at boundary face.
Definition: cs_convection_diffusion_priv.h:3270
cs_i_diff_flux
static CS_F_HOST_DEVICE 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])
Add diffusive fluxes to fluxes at face ij.
Definition: cs_convection_diffusion_priv.h:1003
cs_i_cd_steady_strided
static CS_F_HOST_DEVICE void cs_i_cd_steady_strided(cs_real_t bldfrp, int ischcp, double relaxp, double blencp, 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_rreal_t diipf[3], const cs_rreal_t djjpf[3], const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], const cs_real_t pia[stride], const cs_real_t pja[stride], cs_real_t pifri[stride], cs_real_t pifrj[stride], cs_real_t pjfri[stride], cs_real_t pjfrj[stride], cs_real_t pip[stride], cs_real_t pjp[stride], cs_real_t pipr[stride], cs_real_t pjpr[stride])
Handle preparation of internal face values for the fluxes computation in case of a steady algorithm a...
Definition: cs_convection_diffusion_priv.h:1544
cs_solu_f_val_strided
static CS_F_HOST_DEVICE void cs_solu_f_val_strided(const cs_real_t cell_cen[3], const cs_real_t i_face_cog[3], const cs_real_t grad[stride][3], const cs_real_t p[stride], cs_real_t pf[stride])
Prepare value at face ij by using a Second Order Linear Upwind scheme.
Definition: cs_convection_diffusion_priv.h:826
cs_i_cd_steady
static CS_F_HOST_DEVICE 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_rreal_t diipf[3], const cs_rreal_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)
Handle preparation of internal face values for the fluxes computation in case of a steady algorithm a...
Definition: cs_convection_diffusion_priv.h:1359
cs_b_compute_quantities_strided
static CS_F_HOST_DEVICE void cs_b_compute_quantities_strided(const cs_rreal_t diipb[3], const cs_real_t gradi[stride][3], const cs_real_t bldfrp, cs_real_t recoi[stride])
Reconstruct values in I' at boundary cell i.
Definition: cs_convection_diffusion_priv.h:2913
cs_central_downwind_cells
static CS_F_HOST_DEVICE 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)
Determine the upwind and downwind sides of an internal face and matching cell indices.
Definition: cs_convection_diffusion_priv.h:2579
cs_i_conv_flux
static CS_F_HOST_DEVICE 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)
Add convective fluxes (substracting the mass accumulation from them) to fluxes at face ij.
Definition: cs_convection_diffusion_priv.h:913
cs_i_compute_quantities_strided
static CS_F_HOST_DEVICE void cs_i_compute_quantities_strided(const cs_real_t bldfrp, const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], cs_real_t recoi[stride], cs_real_t recoj[stride], cs_real_t pip[stride], cs_real_t pjp[stride])
Reconstruct values in I' and J'.
Definition: cs_convection_diffusion_priv.h:584
cs_nvd_scheme_scalar
static CS_F_HOST_DEVICE 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)
Compute the normalized face scalar using the specified NVD scheme.
Definition: cs_convection_diffusion_priv.h:89
cs_b_compute_quantities
static CS_F_HOST_DEVICE void cs_b_compute_quantities(const cs_rreal_t diipb[3], const cs_real_t gradi[3], const cs_real_t bldfrp, cs_real_t *recoi)
Reconstruct values in I' at boundary cell i.
Definition: cs_convection_diffusion_priv.h:2887
cs_i_cd_steady_slope_test
static CS_F_HOST_DEVICE 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 cell_ceni[3], const cs_real_t cell_cenj[3], const cs_nreal_t i_face_u_normal[3], const cs_real_t i_face_cog[3], const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], const cs_real_t i_massflux, 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 gradsti[3], const cs_real_t gradstj[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)
Handle preparation of internal face values for the fluxes computation in case of a steady algorithm a...
Definition: cs_convection_diffusion_priv.h:1988
cs_blend_f_val
static CS_F_HOST_DEVICE void cs_blend_f_val(const double blencp, const cs_real_t p, cs_real_t *pf)
Blend face values for a centered or SOLU scheme with face values for an upwind scheme.
Definition: cs_convection_diffusion_priv.h:857
cs_i_cd_unsteady_strided
static CS_F_HOST_DEVICE void cs_i_cd_unsteady_strided(cs_real_t bldfrp, int ischcp, double blencp, 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 hybrid_blend_i, const cs_real_t hybrid_blend_j, const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], cs_real_t pif[stride], cs_real_t pjf[stride], cs_real_t pip[stride], cs_real_t pjp[stride])
Handle preparation of internal face values for the fluxes computation in case of an unsteady algorith...
Definition: cs_convection_diffusion_priv.h:1841
cs_b_imposed_conv_flux
static CS_F_HOST_DEVICE 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)
Add convective flux (substracting the mass accumulation from it) to flux at boundary face....
Definition: cs_convection_diffusion_priv.h:3008
cs_slope_test_strided
static CS_F_HOST_DEVICE void cs_slope_test_strided(const cs_real_t pi[stride], const cs_real_t pj[stride], const cs_real_t distf, const cs_nreal_t i_face_u_normal[3], const cs_real_t gradi[stride][3], const cs_real_t gradj[stride][3], const cs_real_t gradsti[stride][3], const cs_real_t gradstj[stride][3], const cs_real_t i_massflux, cs_real_t *testij, cs_real_t *tesqck)
Compute slope test criteria at internal face between cell i and j.
Definition: cs_convection_diffusion_priv.h:472
cs_defs.h
CS_F_HOST_DEVICE
#define CS_F_HOST_DEVICE
Definition: cs_defs.h:561
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:342
cs_rreal_3_t
cs_rreal_t cs_rreal_3_t[3]
Definition: cs_defs.h:388
cs_nreal_3_t
cs_nreal_t cs_nreal_3_t[3]
Definition: cs_defs.h:385
CS_ABS
#define CS_ABS(a)
Definition: cs_defs.h:504
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:359
cs_real_2_t
cs_real_t cs_real_2_t[2]
vector of 2 floating-point values
Definition: cs_defs.h:358
CS_UNUSED
#define CS_UNUSED(x)
Definition: cs_defs.h:528
cs_rreal_t
double cs_rreal_t
Definition: cs_defs.h:348
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:335
cs_nreal_t
double cs_nreal_t
Definition: cs_defs.h:346
p
@ p
Definition: cs_field_pointer.h:67
k
@ k
Definition: cs_field_pointer.h:72
cs_halo_sync_var
void cs_halo_sync_var(const cs_halo_t *halo, cs_halo_type_t sync_mode, cs_real_t var[])
Definition: cs_halo.cpp:2084
cs_halo_type_t
cs_halo_type_t
Definition: cs_halo.h:56
cs_math.h
cs_math_3_distance
static CS_F_HOST_DEVICE cs_real_t cs_math_3_distance(const cs_real_t xa[3], const cs_real_t xb[3])
Compute the (euclidean) distance between two points xa and xb in a Cartesian coordinate system of dim...
Definition: cs_math.h:629
cs_math_3_dot_product
static CS_F_HOST_DEVICE cs_real_t cs_math_3_dot_product(const cs_real_t u[3], const cs_real_t v[3])
Compute the dot product of two vectors of 3 real values.
Definition: cs_math.h:696
cs_math_epzero
const cs_real_t cs_math_epzero
cs_math_3_norm
static CS_F_HOST_DEVICE cs_real_t cs_math_3_norm(const cs_real_t v[3])
Compute the euclidean norm of a vector of dimension 3.
Definition: cs_math.h:765
cs_math_fabs
static CS_F_HOST_DEVICE cs_real_t cs_math_fabs(cs_real_t x)
Compute the absolute value of a real value.
Definition: cs_math.h:466
cs_math_fmin
static CS_F_HOST_DEVICE cs_real_t cs_math_fmin(cs_real_t x, cs_real_t y)
Compute the min value of two real values.
Definition: cs_math.h:484
cs_math_sq
static CS_F_HOST_DEVICE cs_real_t cs_math_sq(cs_real_t x)
Compute the square of a real value.
Definition: cs_math.h:545
cs_math_fmax
static CS_F_HOST_DEVICE cs_real_t cs_math_fmax(cs_real_t x, cs_real_t y)
Compute the max value of two real values.
Definition: cs_math.h:503
cs_parameters.h
CS_COUPLED_FD
@ CS_COUPLED_FD
Definition: cs_parameters.h:101
cstnum::pi
double precision pi
value with 16 digits
Definition: cstnum.f90:48
mesh::diipb
double precision, dimension(:,:), pointer diipb
vector II' for interior faces for every boundary face, the three components of the vector ....
Definition: mesh.f90:160
coincl::fmin
real(c_double), pointer, save fmin
Definition: coincl.f90:185
cpincl::b2
real(c_double), dimension(:), pointer, save b2
Definition: cpincl.f90:109
cpincl::b1
real(c_double), dimension(:), pointer, save b1
Definition: cpincl.f90:109
cs_mesh_t
Definition: cs_mesh.h:85
cs_mesh_t::halo
cs_halo_t * halo
Definition: cs_mesh.h:156