9.2
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-2026 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 "base/cs_defs.h"
31
32/*----------------------------------------------------------------------------
33 * Local headers
34 *----------------------------------------------------------------------------*/
35
36#include "base/cs_math.h"
37#include "base/cs_field.h"
38#include "base/cs_parameters.h" // for BC types
39#include "base/cs_reducers.h"
40
42
43/*=============================================================================
44 * Macro definitions
45 *============================================================================*/
46
47/*============================================================================
48 * Type definition
49 *============================================================================*/
50
51/*============================================================================
52 * Global variables
53 *============================================================================*/
54
55/*============================================================================
56 * Public inlined function
57 *============================================================================*/
58
59/*----------------------------------------------------------------------------
60 * Synchronize halos for scalar variables.
61 *
62 * parameters:
63 * m <-- pointer to associated mesh structure
64 * halo_type <-> halo type
65 * pvar <-> variable
66 *----------------------------------------------------------------------------*/
67
68inline static void
70 cs_halo_type_t halo_type,
71 cs_real_t pvar[])
72{
73 if (m->halo != NULL)
74 cs_halo_sync_var(m->halo, halo_type, pvar);
75}
76
77/*----------------------------------------------------------------------------*/
88/*----------------------------------------------------------------------------*/
89
90CS_F_HOST_DEVICE inline static cs_real_t
92 const cs_real_t nvf_p_c,
93 const cs_real_t nvf_r_f,
94 const cs_real_t nvf_r_c)
95{
96 cs_real_t nvf_p_f;
97
98 cs_real_t beta_m, rfc, r1f, r1, r2, r3, b1, b2;
99
100 switch (limiter) {
101 case CS_NVD_GAMMA: /* Gamma scheme */
102 beta_m = 0.1; /* in [0.1, 0.5] */
103 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
104
105 if (nvf_p_c < beta_m) {
106 nvf_p_f = nvf_p_c*(1.+rfc*(1.-nvf_p_c)/beta_m);
107 } else {
108 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
109
110 nvf_p_f = r1f*nvf_p_c+rfc;
111 }
112
113 break;
114
115 case CS_NVD_SMART: /* SMART scheme */
116 if (nvf_p_c < (nvf_r_c/3.)) {
117 r1 = nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
118 r2 = nvf_r_c*(1.-nvf_r_c);
119
120 nvf_p_f = nvf_p_c*r1/r2;
121 } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
122 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
123 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
124
125 nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c + rfc);
126 } else {
127 nvf_p_f = 1.;
128 }
129
130 break;
131
132 case CS_NVD_CUBISTA: /* CUBISTA scheme */
133 if (nvf_p_c < (3.*nvf_r_c/4.)) {
134 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
135
136 nvf_p_f = nvf_r_f*(1.+rfc/3.)*nvf_p_c/nvf_r_c;
137 } else if (nvf_p_c <= (nvf_r_c*(1.+2.*(nvf_r_f-nvf_r_c))/(2.*nvf_r_f-nvf_r_c))) {
138 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
139 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
140
141 nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c+rfc);
142 } else {
143 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
144
145 nvf_p_f = 1.-.5*r1f*(1.-nvf_p_c);
146 }
147
148 break;
149
150 case CS_NVD_SUPERBEE: /* SuperBee scheme */
151 if (nvf_p_c < (nvf_r_c/(2.-nvf_r_c))) {
152 nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
153 } else if (nvf_p_c < nvf_r_c) {
154 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
155 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
156
157 nvf_p_f = r1f*nvf_p_c+rfc;
158 } else if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
159 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
160 } else {
161 nvf_p_f = 1.;
162 }
163
164 break;
165
166 case CS_NVD_MUSCL: /* MUSCL scheme */
167 if (nvf_p_c < (.5*nvf_r_c)) {
168 nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
169 } else if (nvf_p_c < (1.+nvf_r_c-nvf_r_f)) {
170 nvf_p_f = nvf_p_c+nvf_r_f-nvf_r_c;
171 } else {
172 nvf_p_f = 1.;
173 }
174
175 break;
176
177 case CS_NVD_MINMOD: /* MINMOD scheme */
178 if (nvf_p_c < nvf_r_c) {
179 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
180 } else {
181 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
182 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
183
184 nvf_p_f = r1f*nvf_p_c+rfc;
185 }
186
187 break;
188
189 case CS_NVD_CLAM: /* CLAM scheme */
190 r1 = nvf_r_c*nvf_r_c-nvf_r_f;
191 r2 = nvf_r_c*(nvf_r_c-1.);
192 r3 = nvf_r_f-nvf_r_c;
193
194 nvf_p_f = nvf_p_c*(r1+r3*nvf_p_c)/r2;
195
196 break;
197
198 case CS_NVD_STOIC: /* STOIC scheme */
199 b1 = (nvf_r_c-nvf_r_f)*nvf_r_c;
200 b2 = nvf_r_c+nvf_r_f+2.*nvf_r_f*nvf_r_f-4.*nvf_r_f*nvf_r_c;
201
202 if (nvf_p_c < (b1/b2)) {
203 r1 = -nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
204 r2 = nvf_r_c*(nvf_r_c-1.);
205
206 nvf_p_f = nvf_p_c*r1/r2;
207 } else if (nvf_p_c < nvf_r_c) {
208 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
209 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
210
211 nvf_p_f = r1f*nvf_p_c+rfc;
212 } else if (nvf_p_c < (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
213 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
214 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
215
216 nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
217 } else {
218 nvf_p_f = 1.;
219 }
220
221 break;
222
223 case CS_NVD_OSHER: /* OSHER scheme */
224 if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
225 nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
226 } else {
227 nvf_p_f = 1.;
228 }
229
230 break;
231
232 case CS_NVD_WASEB: /* WASEB scheme */
233 r1 = nvf_r_c*nvf_r_f*(nvf_r_f-nvf_r_c);
234 r2 = 2.*nvf_r_c*(1.-nvf_r_c)-nvf_r_f*(1.-nvf_r_f);
235
236 if (nvf_p_c < (r1/r2)) {
237 nvf_p_f = 2.*nvf_p_c;
238 } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
239 rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
240 r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
241
242 nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
243 } else {
244 nvf_p_f = 1.;
245 }
246
247 break;
248
249 default: /* Upwinding */
250 nvf_p_f = nvf_p_c;
251 break;
252 }
253
254 return nvf_p_f;
255}
256
257/*----------------------------------------------------------------------------*/
272/*----------------------------------------------------------------------------*/
273
274template <typename T>
275CS_F_HOST_DEVICE inline static cs_real_t
277 const cs_nreal_t i_face_u_normal[3],
278 cs_real_t nvf_p_c,
279 cs_real_t nvf_r_f,
280 cs_real_t nvf_r_c,
281 const T gradv_c[3],
282 const cs_real_t c_courant)
283{
284 assert(limiter >= CS_NVD_VOF_HRIC);
285
286 cs_real_t nvf_p_f;
287 cs_real_t blend, high_order, low_order, ratio;
288
289 /* Compute gradient angle indicator */
290 cs_real_t dotp = cs::abs(cs_math_3_dot_product(gradv_c, i_face_u_normal));
291 cs_real_t sgrad = cs_math_3_norm(gradv_c);
292 cs_real_t denom = sgrad;
293
294 if (limiter == CS_NVD_VOF_HRIC) { /* M-HRIC scheme */
295 /* High order scheme : Bounded Downwind */
296 if (nvf_p_c <= .5) {
297 high_order = 2.*nvf_p_c;
298 } else {
299 high_order = 1.;
300 }
301
302 /* Low order scheme : MUSCL */
304 nvf_p_c,
305 nvf_r_f,
306 nvf_r_c);
307
308 /* Compute the blending factor */
309 if (denom <= (cs_math_epzero*dotp)) {
310 blend = 1.;
311 } else {
312 ratio = dotp/denom;
313 blend = cs_math_fmin(1., pow(ratio, .5));
314 }
315
316 /* Blending */
317 nvf_p_f = blend*high_order + (1.-blend)*low_order;
318
319 /* Extra blending due to the cell Courant number */
320 if (c_courant < .7 && c_courant > .3) {
321 nvf_p_f = nvf_p_f + (nvf_p_f - low_order)*(.7 - c_courant )/.4;
322 } else if (c_courant >= .7) {
323 nvf_p_f = low_order;
324 }
325 } else if (limiter == CS_NVD_VOF_CICSAM) { /* M-CICSAM scheme */
326 /* High order scheme : HYPER-C + SUPERBEE */
327 if (c_courant <= .3) {
328 high_order = cs_math_fmin(1., nvf_p_c/(c_courant+cs_math_epzero));
329 } else if (c_courant <= .6) {
330 high_order = cs_math_fmin(1., nvf_p_c/.3);
331 } else if (c_courant <= .7) {
333 nvf_p_c,
334 nvf_r_f,
335 nvf_r_c);
336 high_order = 10.*( (.7-c_courant)*cs_math_fmin(1., nvf_p_c/.3)
337 + (c_courant-.6)*superbee);
338 }
339 else {
341 nvf_p_c,
342 nvf_r_f,
343 nvf_r_c);
344 }
345
346 /* Low order scheme : MUSCL */
348 nvf_p_c,
349 nvf_r_f,
350 nvf_r_c);
351
352 /* Compute the blending factor */
353 if (denom <= (cs_math_epzero*dotp)) {
354 blend = 1.;
355 } else {
356 ratio = dotp/denom;
357 blend = cs_math_fmin(1., pow(ratio, 2.));
358 }
359
360 /* Blending */
361 nvf_p_f = blend*high_order + (1.-blend)*low_order;
362 } else { /* STACS scheme */
363 /* High order scheme : SUPERBEE */
365 nvf_p_c,
366 nvf_r_f,
367 nvf_r_c);
368
369 /* Low order scheme : STOIC */
371 nvf_p_c,
372 nvf_r_f,
373 nvf_r_c);
374
375 /* Compute the blending factor */
376 if (denom <= (cs_math_epzero*dotp)) {
377 blend = 1.;
378 } else {
379 ratio = dotp/denom;
380 blend = cs_math_fmin(1., pow(ratio, 4.));
381 }
382
383 /* Blending */
384 nvf_p_f = blend*high_order + (1.-blend)*low_order;
385 }
386
387 return nvf_p_f;
388}
389
390/*----------------------------------------------------------------------------*/
406/*----------------------------------------------------------------------------*/
407
408template <typename T>
409CS_F_HOST_DEVICE inline static void
411 const cs_real_t pj,
412 const cs_real_t distf,
413 const cs_nreal_t i_face_u_normal[3],
414 const T gradi[3],
415 const T gradj[3],
416 const T grdpai[3],
417 const T grdpaj[3],
418 const cs_real_t i_massflux,
419 T *testij,
420 T *tesqck)
421{
422 T dcc, ddi, ddj;
423
424 /* Slope test */
425
426 *testij = grdpai[0]*grdpaj[0]
427 + grdpai[1]*grdpaj[1]
428 + grdpai[2]*grdpaj[2];
429
430 if (i_massflux > 0.) {
431 dcc = gradi[0]*i_face_u_normal[0]
432 + gradi[1]*i_face_u_normal[1]
433 + gradi[2]*i_face_u_normal[2];
434 ddi = grdpai[0]*i_face_u_normal[0]
435 + grdpai[1]*i_face_u_normal[1]
436 + grdpai[2]*i_face_u_normal[2];
437 ddj = (pj-pi)/distf;
438 }
439 else {
440 dcc = gradj[0]*i_face_u_normal[0]
441 + gradj[1]*i_face_u_normal[1]
442 + gradj[2]*i_face_u_normal[2];
443 ddi = (pj-pi)/distf;
444 ddj = grdpaj[0]*i_face_u_normal[0]
445 + grdpaj[1]*i_face_u_normal[1]
446 + grdpaj[2]*i_face_u_normal[2];
447 }
448
449 T diff = ddi-ddj;
450
451 *tesqck = cs_math_sq(dcc) - cs_math_sq(diff);
452}
453
454/*----------------------------------------------------------------------------*/
473/*----------------------------------------------------------------------------*/
474
475template <cs_lnum_t stride, typename T>
476CS_F_HOST_DEVICE inline static void
478 const cs_real_t pj[stride],
479 const cs_real_t distf,
480 const cs_nreal_t i_face_u_normal[3],
481 const T gradi[stride][3],
482 const T gradj[stride][3],
483 const T gradsti[stride][3],
484 const T gradstj[stride][3],
485 const cs_real_t i_massflux,
486 T *testij,
487 T *tesqck)
488{
489 *testij = 0.;
490 *tesqck = 0.;
491
492 /* Slope test */
493
494 for (int i = 0; i < stride; i++) {
495
496 T dcc, ddi, ddj;
497
498 *testij += gradsti[i][0]*gradstj[i][0]
499 + gradsti[i][1]*gradstj[i][1]
500 + gradsti[i][2]*gradstj[i][2];
501
502 if (i_massflux > 0.) {
503 dcc = gradi[i][0]*i_face_u_normal[0]
504 + gradi[i][1]*i_face_u_normal[1]
505 + gradi[i][2]*i_face_u_normal[2];
506 ddi = gradsti[i][0]*i_face_u_normal[0]
507 + gradsti[i][1]*i_face_u_normal[1]
508 + gradsti[i][2]*i_face_u_normal[2];
509 ddj = (pj[i]-pi[i])/distf;
510 }
511 else {
512 dcc = gradj[i][0]*i_face_u_normal[0]
513 + gradj[i][1]*i_face_u_normal[1]
514 + gradj[i][2]*i_face_u_normal[2];
515 ddi = (pj[i]-pi[i])/distf;
516 ddj = gradstj[i][0]*i_face_u_normal[0]
517 + gradstj[i][1]*i_face_u_normal[1]
518 + gradstj[i][2]*i_face_u_normal[2];
519 }
520
521 *tesqck += cs_math_sq(dcc) - cs_math_sq(ddi-ddj);
522 }
523}
524
525/*----------------------------------------------------------------------------*/
541/*----------------------------------------------------------------------------*/
542
543template <typename T>
544CS_F_HOST_DEVICE inline static void
546 const cs_rreal_3_t diipf,
547 const cs_rreal_3_t djjpf,
548 const T gradi[3],
549 const T gradj[3],
550 const cs_real_t pi,
551 const cs_real_t pj,
552 T *recoi,
553 T *recoj,
554 cs_real_t *pip,
555 cs_real_t *pjp)
556{
557 T gradpf[3] = {(T)0.5*(gradi[0] + gradj[0]),
558 (T)0.5*(gradi[1] + gradj[1]),
559 (T)0.5*(gradi[2] + gradj[2])};
560
561 /* reconstruction only if IRCFLP = 1 */
562 *recoi = bldfrp*(cs_math_3_dot_product(gradpf, diipf));
563 *recoj = bldfrp*(cs_math_3_dot_product(gradpf, djjpf));
564 *pip = pi + *recoi;
565 *pjp = pj + *recoj;
566}
567
568/*----------------------------------------------------------------------------*/
587/*----------------------------------------------------------------------------*/
588
589template <cs_lnum_t stride, typename T>
590CS_F_HOST_DEVICE inline static void
592 const cs_rreal_t diipf[3],
593 const cs_rreal_t djjpf[3],
594 const T gradi[stride][3],
595 const T gradj[stride][3],
596 const cs_real_t pi[stride],
597 const cs_real_t pj[stride],
598 T recoi[stride],
599 T recoj[stride],
600 cs_real_t pip[stride],
601 cs_real_t pjp[stride])
602{
603 T dpvf[3];
604
605 /* x-y-z components, p = u, v, w */
606
607 for (int isou = 0; isou < stride; isou++) {
608
609 for (int jsou = 0; jsou < 3; jsou++)
610 dpvf[jsou] = 0.5*( gradi[isou][jsou]
611 + gradj[isou][jsou]);
612
613 /* reconstruction only if IRCFLP = 1 */
614
615 recoi[isou] = bldfrp*(cs_math_3_dot_product(dpvf, diipf));
616 recoj[isou] = bldfrp*(cs_math_3_dot_product(dpvf, djjpf));
617
618 pip[isou] = pi[isou] + recoi[isou];
619 pjp[isou] = pj[isou] + recoj[isou];
620
621 }
622}
623
624/*----------------------------------------------------------------------------*/
631/*----------------------------------------------------------------------------*/
632
633CS_F_HOST_DEVICE inline static void
635 cs_real_t &pip)
636{
637 if (pip < bounds[0]) {
638 pip = bounds[0];
639 }
640 else if (pip > bounds[1]) {
641 pip = bounds[1];
642 }
643}
644
645/*----------------------------------------------------------------------------*/
656/*----------------------------------------------------------------------------*/
657
658template <cs_lnum_t stride>
659CS_F_HOST_DEVICE inline static void
661 const cs_real_t pi[stride],
662 cs_real_t pip[stride])
663{
664 cs_real_t v_r[stride];
665 for (cs_lnum_t isou = 0; isou < stride; isou++)
666 v_r[isou] = pip[isou] - pi[isou];
667
668 cs_real_t d2 = cs_math_square_norm<stride>(v_r);
669 if (d2 > bounds) {
670 cs_real_t s = sqrt(bounds / d2); // scaling factor
671
672 for (cs_lnum_t isou = 0; isou < stride; isou++)
673 pip[isou] = pip[isou] + v_r[isou]*s;
674 }
675}
676
677/*----------------------------------------------------------------------------*/
693/*----------------------------------------------------------------------------*/
694
695CS_F_HOST_DEVICE inline static void
696cs_i_relax_c_val(const double relaxp,
697 const cs_real_t pia,
698 const cs_real_t pja,
699 const cs_real_t recoi,
700 const cs_real_t recoj,
701 const cs_real_t pi,
702 const cs_real_t pj,
703 cs_real_t *pir,
704 cs_real_t *pjr,
705 cs_real_t *pipr,
706 cs_real_t *pjpr)
707{
708 *pir = pi/relaxp - (1.-relaxp)/relaxp * pia;
709 *pjr = pj/relaxp - (1.-relaxp)/relaxp * pja;
710
711 *pipr = *pir + recoi;
712 *pjpr = *pjr + recoj;
713}
714
715/*----------------------------------------------------------------------------*/
734/*----------------------------------------------------------------------------*/
735
736template <cs_lnum_t stride, typename T>
737CS_F_HOST_DEVICE inline static void
739 const cs_real_t pia[stride],
740 const cs_real_t pja[stride],
741 const T recoi[stride],
742 const T recoj[stride],
743 const cs_real_t pi[stride],
744 const cs_real_t pj[stride],
745 cs_real_t pir[stride],
746 cs_real_t pjr[stride],
747 cs_real_t pipr[stride],
748 cs_real_t pjpr[stride])
749{
750 for (int isou = 0; isou < stride; isou++) {
751 pir[isou] = pi[isou] /relaxp - (1.-relaxp)/relaxp * pia[isou];
752 pjr[isou] = pj[isou] /relaxp - (1.-relaxp)/relaxp * pja[isou];
753
754 pipr[isou] = pir[isou] + recoi[isou];
755 pjpr[isou] = pjr[isou] + recoj[isou];
756 }
757}
758
759/*----------------------------------------------------------------------------*/
768/*----------------------------------------------------------------------------*/
769
770CS_F_HOST_DEVICE inline static void
772 cs_real_t pip,
773 cs_real_t pjp,
774 cs_real_t *pf)
775{
776 *pf = pnd*pip + (1.-pnd)*pjp;
777}
778
779/*----------------------------------------------------------------------------*/
791/*----------------------------------------------------------------------------*/
792
793template <cs_lnum_t stride>
794CS_F_HOST_DEVICE inline static void
796 const cs_real_t pip[stride],
797 const cs_real_t pjp[stride],
798 cs_real_t pf[stride])
799{
800 for (int isou = 0; isou < stride; isou++)
801 pf[isou] = pnd*pip[isou] + (1.-pnd)*pjp[isou];
802}
803
804/*----------------------------------------------------------------------------*/
814/*----------------------------------------------------------------------------*/
815
816template <typename T>
817CS_F_HOST_DEVICE inline static void
818cs_solu_f_val(const cs_real_t cell_cen[3],
819 const cs_real_t i_face_cog[3],
820 const T grad[3],
821 const cs_real_t p,
822 cs_real_t *pf)
823{
824 cs_real_t df[3];
825
826 df[0] = i_face_cog[0] - cell_cen[0];
827 df[1] = i_face_cog[1] - cell_cen[1];
828 df[2] = i_face_cog[2] - cell_cen[2];
829
830 *pf = p + cs_math_3_dot_product(df, grad);
831}
832
833/*----------------------------------------------------------------------------*/
846/*----------------------------------------------------------------------------*/
847
848template <cs_lnum_t stride, typename T>
849CS_F_HOST_DEVICE inline static void
851 const cs_real_t i_face_cog[3],
852 const T grad[stride][3],
853 const cs_real_t p[stride],
854 cs_real_t pf[stride])
855{
856 cs_real_t df[3];
857
858 for (cs_lnum_t jsou = 0; jsou < 3; jsou++)
859 df[jsou] = i_face_cog[jsou] - cell_cen[jsou];
860
861 for (cs_lnum_t isou = 0; isou < stride; isou++) {
862 pf[isou] = p[isou] + df[0]*grad[isou][0]
863 + df[1]*grad[isou][1]
864 + df[2]*grad[isou][2];
865 }
866}
867
868/*----------------------------------------------------------------------------*/
878/*----------------------------------------------------------------------------*/
879
880CS_F_HOST_DEVICE inline static void
881cs_blend_f_val(const double blencp,
882 const cs_real_t p,
883 cs_real_t *pf)
884{
885 *pf = blencp * (*pf) + (1. - blencp) * p;
886}
887
888/*----------------------------------------------------------------------------*/
901/*----------------------------------------------------------------------------*/
902
903template <cs_lnum_t stride>
904CS_F_HOST_DEVICE inline static void
905cs_blend_f_val_strided(const double blencp,
906 const cs_real_t p[stride],
907 cs_real_t pf[stride])
908{
909 for (int isou = 0; isou < stride; isou++)
910 pf[isou] = blencp*(pf[isou])+(1.-blencp)*p[isou];
911}
912
913/*----------------------------------------------------------------------------*/
934/*----------------------------------------------------------------------------*/
935
936CS_F_HOST_DEVICE inline static void
937cs_i_conv_flux(const int iconvp,
938 const cs_real_t thetap,
939 const int imasac,
940 const cs_real_t pi,
941 const cs_real_t pj,
942 const cs_real_t pifri,
943 const cs_real_t pifrj,
944 const cs_real_t pjfri,
945 const cs_real_t pjfrj,
946 const cs_real_t i_massflux,
947 const cs_real_t xcppi,
948 const cs_real_t xcppj,
949 cs_real_2_t fluxij)
950{
951 cs_real_t flui, fluj;
952
953 flui = 0.5*(i_massflux + fabs(i_massflux));
954 fluj = 0.5*(i_massflux - fabs(i_massflux));
955
956 fluxij[0] += iconvp*xcppi*(thetap*(flui*pifri + fluj*pjfri) - imasac*i_massflux*pi);
957 fluxij[1] += iconvp*xcppj*(thetap*(flui*pifrj + fluj*pjfrj) - imasac*i_massflux*pj);
958}
959
960/*----------------------------------------------------------------------------*/
981/*----------------------------------------------------------------------------*/
982
983template <cs_lnum_t stride>
984CS_F_HOST_DEVICE inline static void
986 cs_real_t thetap,
987 int imasac,
988 const cs_real_t pi[stride],
989 const cs_real_t pj[stride],
990 const cs_real_t pifri[stride],
991 const cs_real_t pifrj[stride],
992 const cs_real_t pjfri[stride],
993 const cs_real_t pjfrj[stride],
994 cs_real_t i_massflux,
995 cs_real_t fluxi[stride],
996 cs_real_t fluxj[stride])
997{
998 cs_real_t flui, fluj;
999
1000 flui = 0.5*(i_massflux + fabs(i_massflux));
1001 fluj = 0.5*(i_massflux - fabs(i_massflux));
1002
1003 for (int isou = 0; isou < stride; isou++) {
1004 fluxi[isou] += iconvp*( thetap*(flui*pifri[isou] + fluj*pjfri[isou])
1005 - imasac*i_massflux*pi[isou]);
1006 fluxj[isou] += iconvp*( thetap*(flui*pifrj[isou] + fluj*pjfrj[isou])
1007 - imasac*i_massflux*pj[isou]);
1008 }
1009}
1010
1011/*----------------------------------------------------------------------------*/
1024/*----------------------------------------------------------------------------*/
1025
1026CS_F_HOST_DEVICE inline static void
1027cs_i_diff_flux(const int idiffp,
1028 const cs_real_t thetap,
1029 const cs_real_t pip,
1030 const cs_real_t pjp,
1031 const cs_real_t pipr,
1032 const cs_real_t pjpr,
1033 const cs_real_t i_visc,
1034 cs_real_t fluxij[2])
1035{
1036 fluxij[0] += idiffp*thetap*i_visc*(pipr -pjp);
1037 fluxij[1] += idiffp*thetap*i_visc*(pip -pjpr);
1038}
1039
1040/*----------------------------------------------------------------------------*/
1057/*----------------------------------------------------------------------------*/
1058
1059template <cs_lnum_t stride>
1060CS_F_HOST_DEVICE inline static void
1062 cs_real_t thetap,
1063 const cs_real_t pip[stride],
1064 const cs_real_t pjp[stride],
1065 const cs_real_t pipr[stride],
1066 const cs_real_t pjpr[stride],
1067 const cs_real_t i_visc,
1068 cs_real_t fluxi[stride],
1069 cs_real_t fluxj[stride])
1070{
1071 for (int isou = 0; isou < stride; isou++) {
1072 fluxi[isou] += idiffp*thetap*i_visc*(pipr[isou] -pjp[isou]);
1073 fluxj[isou] += idiffp*thetap*i_visc*(pip[isou] -pjpr[isou]);
1074 }
1075}
1076
1077/*----------------------------------------------------------------------------*/
1094/*----------------------------------------------------------------------------*/
1095
1096CS_F_HOST_DEVICE inline static void
1098 const cs_rreal_t diipf[3],
1099 const cs_rreal_t djjpf[3],
1100 const cs_real_t gradi[3],
1101 const cs_real_t gradj[3],
1102 const cs_real_t pi,
1103 const cs_real_t pj,
1104 cs_real_t *pif,
1105 cs_real_t *pjf,
1106 cs_real_t *pip,
1107 cs_real_t *pjp)
1108{
1109 cs_real_t recoi, recoj;
1110
1112 diipf,
1113 djjpf,
1114 gradi,
1115 gradj,
1116 pi,
1117 pj,
1118 &recoi,
1119 &recoj,
1120 pip,
1121 pjp);
1122
1123 *pif = pi;
1124 *pjf = pj;
1125}
1126
1127/*----------------------------------------------------------------------------*/
1147/*----------------------------------------------------------------------------*/
1148
1149template <cs_lnum_t stride>
1150CS_F_HOST_DEVICE inline static void
1152 const cs_rreal_t diipf[3],
1153 const cs_rreal_t djjpf[3],
1154 const cs_real_t gradi[stride][3],
1155 const cs_real_t gradj[stride][3],
1156 const cs_real_t pi[stride],
1157 const cs_real_t pj[stride],
1158 cs_real_t pif[stride],
1159 cs_real_t pjf[stride],
1160 cs_real_t pip[stride],
1161 cs_real_t pjp[stride])
1162{
1163 cs_real_t recoi[stride], recoj[stride];
1164
1165 cs_i_compute_quantities_strided<stride>(bldfrp,
1166 diipf,
1167 djjpf,
1168 gradi,
1169 gradj,
1170 pi,
1171 pj,
1172 recoi,
1173 recoj,
1174 pip,
1175 pjp);
1176
1177 for (int l = 0; l < stride; l++) {
1178 pif[l] = pi[l];
1179 pjf[l] = pj[l];
1180 }
1181}
1182
1183/*----------------------------------------------------------------------------*/
1211/*----------------------------------------------------------------------------*/
1212
1213CS_F_HOST_DEVICE inline static void
1215 const int ischcp,
1216 const double blencp,
1217 const cs_real_t weight,
1218 const cs_real_3_t cell_ceni,
1219 const cs_real_3_t cell_cenj,
1220 const cs_real_3_t i_face_cog,
1221 const cs_real_t hybrid_blend_i,
1222 const cs_real_t hybrid_blend_j,
1223 const cs_rreal_3_t diipf,
1224 const cs_rreal_3_t djjpf,
1225 const cs_real_3_t gradi,
1226 const cs_real_3_t gradj,
1227 const cs_real_3_t gradupi,
1228 const cs_real_3_t gradupj,
1229 const cs_real_t pi,
1230 const cs_real_t pj,
1231 cs_real_t *pif,
1232 cs_real_t *pjf,
1233 cs_real_t *pip,
1234 cs_real_t *pjp)
1235{
1236 cs_real_t recoi, recoj;
1237
1239 diipf,
1240 djjpf,
1241 gradi,
1242 gradj,
1243 pi,
1244 pj,
1245 &recoi,
1246 &recoj,
1247 pip,
1248 pjp);
1249
1250
1251 if (ischcp == 1) {
1252
1253 /* Centered
1254 --------*/
1255
1256 cs_centered_f_val(weight,
1257 *pip,
1258 *pjp,
1259 pif);
1260 cs_centered_f_val(weight,
1261 *pip,
1262 *pjp,
1263 pjf);
1264
1265 } else if (ischcp == 0) {
1266
1267 /* Legacy SOLU
1268 -----------*/
1269
1270 cs_solu_f_val(cell_ceni,
1271 i_face_cog,
1272 gradi,
1273 pi,
1274 pif);
1275 cs_solu_f_val(cell_cenj,
1276 i_face_cog,
1277 gradj,
1278 pj,
1279 pjf);
1280
1281 } else if (ischcp == 3) {
1282
1283 /* Centered
1284 --------*/
1285
1286 cs_centered_f_val(weight,
1287 *pip,
1288 *pjp,
1289 pif);
1290 cs_centered_f_val(weight,
1291 *pip,
1292 *pjp,
1293 pjf);
1294
1295 /* Legacy SOLU
1296 -----------*/
1297 cs_real_t pif_up, pjf_up;
1298 cs_real_t hybrid_blend_interp;
1299
1300 cs_solu_f_val(cell_ceni,
1301 i_face_cog,
1302 gradi,
1303 pi,
1304 &pif_up);
1305 cs_solu_f_val(cell_cenj,
1306 i_face_cog,
1307 gradj,
1308 pj,
1309 &pjf_up);
1310
1311 hybrid_blend_interp = fmin(hybrid_blend_i,hybrid_blend_j);
1312 *pif = hybrid_blend_interp*(*pif) + (1. - hybrid_blend_interp)*pif_up;
1313 *pjf = hybrid_blend_interp*(*pjf) + (1. - hybrid_blend_interp)*pjf_up;
1314
1315 } else {
1316
1317 /* SOLU
1318 ----*/
1319
1320 cs_solu_f_val(cell_ceni,
1321 i_face_cog,
1322 gradupi,
1323 pi,
1324 pif);
1325 cs_solu_f_val(cell_cenj,
1326 i_face_cog,
1327 gradupj,
1328 pj,
1329 pjf);
1330
1331 }
1332
1333 /* Blending
1334 --------*/
1335
1336 cs_blend_f_val(blencp,
1337 pi,
1338 pif);
1339 cs_blend_f_val(blencp,
1340 pj,
1341 pjf);
1342}
1343
1344/*----------------------------------------------------------------------------*/
1373/*----------------------------------------------------------------------------*/
1374
1375template <cs_lnum_t stride>
1376CS_F_HOST_DEVICE inline static void
1378 int ischcp,
1379 double blencp,
1380 cs_real_t weight,
1381 const cs_real_t cell_ceni[3],
1382 const cs_real_t cell_cenj[3],
1383 const cs_real_t i_face_cog[3],
1384 const cs_real_t hybrid_blend_i,
1385 const cs_real_t hybrid_blend_j,
1386 const cs_rreal_t diipf[3],
1387 const cs_rreal_t djjpf[3],
1388 const cs_real_t gradi[stride][3],
1389 const cs_real_t gradj[stride][3],
1390 const cs_real_t pi[stride],
1391 const cs_real_t pj[stride],
1392 cs_real_t pif[stride],
1393 cs_real_t pjf[stride],
1394 cs_real_t pip[stride],
1395 cs_real_t pjp[stride])
1396
1397{
1398 cs_real_t recoi[stride], recoj[stride];
1399
1400 cs_i_compute_quantities_strided<stride>(bldfrp,
1401 diipf,
1402 djjpf,
1403 gradi,
1404 gradj,
1405 pi,
1406 pj,
1407 recoi,
1408 recoj,
1409 pip,
1410 pjp);
1411
1412 if (ischcp == 1) {
1413
1414 /* Centered
1415 --------*/
1416
1417 cs_centered_f_val_strided<stride>(weight, pip, pjp, pif);
1418 cs_centered_f_val_strided<stride>(weight, pip, pjp, pjf);
1419
1420 }
1421 else if (ischcp == 3) {
1422
1423 /* Centered
1424 --------*/
1425
1426 cs_centered_f_val_strided<stride>(weight, pip, pjp, pif);
1427 cs_centered_f_val_strided<stride>(weight, pip, pjp, pjf);
1428
1429 /* SOLU
1430 -----*/
1431 cs_real_t pif_up[stride], pjf_up[stride];
1432 cs_real_t hybrid_blend_interp;
1433
1434 cs_solu_f_val_strided<stride>(cell_ceni,
1435 i_face_cog,
1436 gradi,
1437 pi,
1438 pif_up);
1439 cs_solu_f_val_strided<stride>(cell_cenj,
1440 i_face_cog,
1441 gradj,
1442 pj,
1443 pjf_up);
1444
1445 hybrid_blend_interp = fmin(hybrid_blend_i, hybrid_blend_j);
1446 for (int isou = 0; isou < stride; isou++) {
1447 pif[isou] = hybrid_blend_interp *pif[isou]
1448 + (1. - hybrid_blend_interp)*pif_up[isou];
1449 pjf[isou] = hybrid_blend_interp *pjf[isou]
1450 + (1. - hybrid_blend_interp)*pjf_up[isou];
1451 }
1452 }
1453 else {
1454
1455 /* Second order
1456 ------------*/
1457
1458 cs_solu_f_val_strided<stride>(cell_ceni,
1459 i_face_cog,
1460 gradi,
1461 pi,
1462 pif);
1463 cs_solu_f_val_strided<stride>(cell_cenj,
1464 i_face_cog,
1465 gradj,
1466 pj,
1467 pjf);
1468
1469 }
1470
1471 /* Blending
1472 --------*/
1473
1474 cs_blend_f_val_strided<stride>(blencp, pi, pif);
1475 cs_blend_f_val_strided<stride>(blencp, pj, pjf);
1476}
1477
1478/*----------------------------------------------------------------------------*/
1514/*----------------------------------------------------------------------------*/
1515
1516CS_F_HOST_DEVICE inline static void
1518 const int iconvp,
1519 const cs_real_t bldfrp,
1520 const int ischcp,
1521 const double blencp,
1522 const double blend_st,
1523 const cs_real_t weight,
1524 const cs_real_t i_dist,
1525 const cs_real_3_t cell_ceni,
1526 const cs_real_3_t cell_cenj,
1527 const cs_nreal_3_t i_face_u_normal,
1528 const cs_real_3_t i_face_cog,
1529 const cs_rreal_3_t diipf,
1530 const cs_rreal_3_t djjpf,
1531 const cs_real_t i_massflux,
1532 const cs_real_3_t gradi,
1533 const cs_real_3_t gradj,
1534 const cs_real_3_t gradupi,
1535 const cs_real_3_t gradupj,
1536 const cs_real_3_t gradsti,
1537 const cs_real_3_t gradstj,
1538 const cs_real_t pi,
1539 const cs_real_t pj,
1540 cs_real_t *pif,
1541 cs_real_t *pjf,
1542 cs_real_t *pip,
1543 cs_real_t *pjp)
1544{
1545 CS_UNUSED(blend_st);
1546
1547 cs_real_t recoi, recoj;
1548 cs_real_t testij, tesqck;
1549
1550 *upwind_switch = false;
1551
1553 diipf,
1554 djjpf,
1555 gradi,
1556 gradj,
1557 pi,
1558 pj,
1559 &recoi,
1560 &recoj,
1561 pip,
1562 pjp);
1563
1564 /* Convection slope test is needed only when iconv >0 */
1565 if (iconvp > 0) {
1567 pj,
1568 i_dist,
1569 i_face_u_normal,
1570 gradi,
1571 gradj,
1572 gradsti,
1573 gradstj,
1574 i_massflux,
1575 &testij,
1576 &tesqck);
1577
1578 if (ischcp==1) {
1579
1580 /* Centered
1581 --------*/
1582
1583 cs_centered_f_val(weight,
1584 *pip,
1585 *pjp,
1586 pif);
1587 cs_centered_f_val(weight,
1588 *pip,
1589 *pjp,
1590 pjf);
1591
1592 } else if (ischcp == 0) {
1593
1594 /* Original SOLU
1595 --------------*/
1596
1597 cs_solu_f_val(cell_ceni,
1598 i_face_cog,
1599 gradi,
1600 pi,
1601 pif);
1602 cs_solu_f_val(cell_cenj,
1603 i_face_cog,
1604 gradj,
1605 pj,
1606 pjf);
1607
1608 } else {
1609
1610 /* SOLU
1611 -----*/
1612
1613 cs_solu_f_val(cell_ceni,
1614 i_face_cog,
1615 gradupi,
1616 pi,
1617 pif);
1618 cs_solu_f_val(cell_cenj,
1619 i_face_cog,
1620 gradupj,
1621 pj,
1622 pjf);
1623
1624 }
1625
1626 /* Slope test: percentage of upwind
1627 -------------------------------- */
1628
1629 if (tesqck<=0. || testij<=0.) {
1630
1631 cs_blend_f_val(blend_st,
1632 pi,
1633 pif);
1634 cs_blend_f_val(blend_st,
1635 pj,
1636 pjf);
1637
1638 *upwind_switch = true;
1639
1640 }
1641
1642 /* Blending
1643 --------*/
1644
1645 cs_blend_f_val(blencp,
1646 pi,
1647 pif);
1648 cs_blend_f_val(blencp,
1649 pj,
1650 pjf);
1651
1652 /* If iconv=0 p*f are useless */
1653 }
1654 else {
1655 *pif = pi;
1656 *pjf = pj;
1657 }
1658
1659}
1660
1661/*----------------------------------------------------------------------------*/
1672/*----------------------------------------------------------------------------*/
1673
1674CS_F_HOST_DEVICE inline static void
1676 const cs_lnum_t jj,
1677 const cs_real_t i_massflux,
1678 cs_lnum_t *ic,
1679 cs_lnum_t *id)
1680{
1681 if (i_massflux >= 0.) {
1682 *ic = ii;
1683 *id = jj;
1684 } else {
1685 *ic = jj;
1686 *id = ii;
1687 }
1688}
1689
1690/*----------------------------------------------------------------------------*/
1711/*----------------------------------------------------------------------------*/
1712
1713template <typename T>
1714CS_F_HOST_DEVICE inline static void
1716 const double beta,
1717 const cs_real_3_t cell_cen_c,
1718 const cs_real_3_t cell_cen_d,
1719 const cs_nreal_3_t i_face_u_normal,
1720 const cs_real_3_t i_face_cog,
1721 const T gradv_c,
1722 const cs_real_t p_c,
1723 const cs_real_t p_d,
1724 const cs_real_t local_max_c,
1725 const cs_real_t local_min_c,
1726 const cs_real_t courant_c,
1727 cs_real_t *pif,
1728 cs_real_t *pjf)
1729{
1730 /* distance between face center and central cell center */
1731 /* Distance between face center and central cell center */
1732 cs_real_t dist_fc = cs_math_3_distance(cell_cen_c, i_face_cog);
1733
1734 /* Unit vector and distance between central and downwind cells centers */
1735 cs_real_t diff[3] = {cell_cen_d[0] - cell_cen_c[0],
1736 cell_cen_d[1] - cell_cen_c[1],
1737 cell_cen_d[2] - cell_cen_c[2]};
1738
1739 cs_real_t dist_dc = cs_math_3_norm(diff);
1740 cs_real_t invl = 1./dist_dc;
1741
1742 cs_real_t ndc[3] = {invl*diff[0],
1743 invl*diff[1],
1744 invl*diff[2]};
1745
1746 /* Place the upwind point on the line that joins
1747 the two cells on the upwind side and the same
1748 distance as that between the two cells */
1749 const cs_real_t dist_cu = dist_dc;
1750 const cs_real_t dist_du = dist_dc + dist_cu;
1751
1752 /* Compute the property on the upwind assuming a parabolic
1753 variation of the property between the two cells */
1754 const cs_real_t gradc = cs_math_3_dot_product(gradv_c, ndc);
1755
1756 const cs_real_t grad2c = ((p_d - p_c)/dist_dc - gradc)/dist_dc;
1757
1758 cs_real_t p_u = p_c + (grad2c*dist_cu - gradc)*dist_cu;
1759 p_u = cs::max(cs::min(p_u, local_max_c), local_min_c);
1760
1761 /* Compute the normalised distances */
1762 const cs_real_t nvf_r_f = (dist_fc+dist_cu)/dist_du;
1763 const cs_real_t nvf_r_c = dist_cu/dist_du;
1764
1765 /* Check for the bounds of the NVD diagram and compute the face
1766 property according to the selected NVD scheme */
1767 const cs_real_t _small
1768 = cs_math_epzero * (cs::abs(p_u) + cs::abs(p_c) + cs::abs(p_d));
1769
1770 if (cs::abs(p_d-p_u) <= _small) {
1771 *pif = p_c;
1772 *pjf = p_c;
1773 }
1774 else {
1775 const cs_real_t nvf_p_c = (p_c - p_u)/(p_d - p_u);
1776
1777 if (nvf_p_c <= 0. || nvf_p_c >= 1.) {
1778 *pif = p_c;
1779 *pjf = p_c;
1780 }
1781 else {
1782 cs_real_t nvf_p_f;
1783
1784 /* Highly compressive NVD scheme for VOF */
1785 if (limiter >= CS_NVD_VOF_HRIC) {
1786 nvf_p_f = cs_nvd_vof_scheme_scalar(limiter,
1787 i_face_u_normal,
1788 nvf_p_c,
1789 nvf_r_f,
1790 nvf_r_c,
1791 gradv_c,
1792 courant_c);
1793 } else { /* Regular NVD scheme */
1794 nvf_p_f = cs_nvd_scheme_scalar(limiter,
1795 nvf_p_c,
1796 nvf_r_f,
1797 nvf_r_c);
1798 }
1799
1800 *pif = p_u + nvf_p_f*(p_d - p_u);
1801 *pif = cs::max(cs::min(*pif, local_max_c), local_min_c);
1802
1803 cs_blend_f_val(beta,
1804 p_c,
1805 pif);
1806
1807 *pjf = *pif;
1808 }
1809 }
1810}
1811
1812/*----------------------------------------------------------------------------*/
1849/*----------------------------------------------------------------------------*/
1850
1851template <cs_lnum_t stride>
1852CS_F_HOST_DEVICE inline static void
1854 int iconvp,
1855 cs_real_t bldfrp,
1856 int ischcp,
1857 double blencp,
1858 double blend_st,
1859 cs_real_t weight,
1860 cs_real_t i_dist,
1861 const cs_real_t cell_ceni[3],
1862 const cs_real_t cell_cenj[3],
1863 const cs_nreal_t i_face_u_normal[3],
1864 const cs_real_t i_face_cog[3],
1865 const cs_rreal_t diipf[3],
1866 const cs_rreal_t djjpf[3],
1867 cs_real_t i_massflux,
1868 const cs_real_t gradi[stride][3],
1869 const cs_real_t gradj[stride][3],
1870 const cs_real_t grdpai[stride][3],
1871 const cs_real_t grdpaj[stride][3],
1872 const cs_real_t pi[stride],
1873 const cs_real_t pj[stride],
1874 cs_real_t pif[stride],
1875 cs_real_t pjf[stride],
1876 cs_real_t pip[stride],
1877 cs_real_t pjp[stride])
1878{
1879 cs_real_t recoi[stride], recoj[stride];
1880 cs_real_t testij, tesqck;
1881 int isou;
1882
1883 cs_i_compute_quantities_strided<stride>(bldfrp,
1884 diipf,
1885 djjpf,
1886 gradi,
1887 gradj,
1888 pi,
1889 pj,
1890 recoi,
1891 recoj,
1892 pip,
1893 pjp);
1894
1895 /* Convection slope test is needed only when iconv >0 */
1896 if (iconvp > 0) {
1897 cs_slope_test_strided<stride>(pi,
1898 pj,
1899 i_dist,
1900 i_face_u_normal,
1901 gradi,
1902 gradj,
1903 grdpai,
1904 grdpaj,
1905 i_massflux,
1906 &testij,
1907 &tesqck);
1908
1909 for (isou = 0; isou < stride; isou++) {
1910
1911 if (ischcp==1) {
1912
1913 /* Centered
1914 --------*/
1915
1916 cs_centered_f_val(weight,
1917 pip[isou],
1918 pjp[isou],
1919 &pif[isou]);
1920 cs_centered_f_val(weight,
1921 pip[isou],
1922 pjp[isou],
1923 &pjf[isou]);
1924
1925 }
1926 else {
1927
1928 /* Second order
1929 ------------*/
1930
1931 cs_solu_f_val(cell_ceni,
1932 i_face_cog,
1933 gradi[isou],
1934 pi[isou],
1935 &pif[isou]);
1936 cs_solu_f_val(cell_cenj,
1937 i_face_cog,
1938 gradj[isou],
1939 pj[isou],
1940 &pjf[isou]);
1941 }
1942
1943 }
1944
1945 /* Slope test activated: percentage of upwind */
1946 if (tesqck <= 0. || testij <= 0.) {
1947
1948 /* Upwind
1949 --------*/
1950
1951 cs_blend_f_val_strided<stride>(blend_st, pi, pif);
1952 cs_blend_f_val_strided<stride>(blend_st, pj, pjf);
1953
1954 *upwind_switch = true;
1955 }
1956
1957 /* Blending
1958 --------*/
1959
1960 cs_blend_f_val_strided<stride>(blencp, pi, pif);
1961 cs_blend_f_val_strided<stride>(blencp, pj, pjf);
1962
1963 /* If iconv=0 p*fr* are useless */
1964 }
1965 else {
1966
1967 for (int l = 0; l < stride; l++) {
1968 pif[l] = pi[l];
1969 pjf[l] = pj[l];
1970 }
1971
1972 }
1973}
1974
1975/*----------------------------------------------------------------------------*/
1984/*----------------------------------------------------------------------------*/
1985
1986template <typename T>
1987CS_F_HOST_DEVICE inline static void
1989 const T gradi[3],
1990 const T bldfrp,
1991 T *recoi)
1992{
1993 *recoi = bldfrp * ( gradi[0]*diipb[0]
1994 + gradi[1]*diipb[1]
1995 + gradi[2]*diipb[2]);
1996}
1997
1998/*----------------------------------------------------------------------------*/
2010/*----------------------------------------------------------------------------*/
2011
2012template <cs_lnum_t stride, typename T>
2013CS_F_HOST_DEVICE inline static void
2015 const T gradi[stride][3],
2016 const T bldfrp,
2017 T recoi[stride])
2018{
2019 for (int isou = 0; isou < stride; isou++) {
2020 recoi[isou] = bldfrp * ( gradi[isou][0]*diipb[0]
2021 + gradi[isou][1]*diipb[1]
2022 + gradi[isou][2]*diipb[2]);
2023 }
2024}
2025
2026/*----------------------------------------------------------------------------*/
2037/*----------------------------------------------------------------------------*/
2038
2039CS_F_HOST_DEVICE inline static void
2040cs_b_relax_c_val(const double relaxp,
2041 const cs_real_t pi,
2042 const cs_real_t pia,
2043 const cs_real_t recoi,
2044 cs_real_t *pir,
2045 cs_real_t *pipr)
2046{
2047 *pir = pi/relaxp - (1.-relaxp)/relaxp*pia;
2048 *pipr = *pir + recoi;
2049}
2050
2051/*----------------------------------------------------------------------------*/
2065/*----------------------------------------------------------------------------*/
2066
2067template <cs_lnum_t stride>
2068CS_F_HOST_DEVICE inline static void
2069cs_b_relax_c_val_strided(const double relaxp,
2070 const cs_real_t pi[stride],
2071 const cs_real_t pia[stride],
2072 const cs_real_t recoi[stride],
2073 cs_real_t pir[stride],
2074 cs_real_t pipr[stride])
2075{
2076 for (int isou = 0; isou < stride; isou++) {
2077 pir[isou] = pi[isou]/relaxp - (1.-relaxp)/relaxp*pia[isou];
2078 pipr[isou] = pir[isou] + recoi[isou];
2079 }
2080}
2081
2082/*----------------------------------------------------------------------------*/
2105/*----------------------------------------------------------------------------*/
2106
2107CS_F_HOST_DEVICE inline static void
2109 cs_real_t thetap,
2110 int imasac,
2111 int inc,
2112 int bc_type,
2113 int icvfli,
2114 cs_real_t pi,
2115 cs_real_t pir,
2116 cs_real_t pipr,
2117 cs_real_t coface,
2118 cs_real_t cofbce,
2119 cs_real_t b_massflux,
2120 cs_real_t xcpp,
2121 cs_real_t val_f,
2122 cs_real_t *flux)
2123{
2124 cs_real_t flui, fluj, pfac;
2125
2126 /* Computed convective flux */
2127
2128 if (icvfli == 0) {
2129
2130 /* Remove decentering for coupled faces */
2131 if (bc_type == CS_COUPLED_FD) {
2132 flui = 0.0;
2133 fluj = b_massflux;
2134 } else {
2135 flui = 0.5*(b_massflux +fabs(b_massflux));
2136 fluj = 0.5*(b_massflux -fabs(b_massflux));
2137 }
2138
2139 *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*val_f) -imasac*( b_massflux*pi));
2140
2141 /* Imposed convective flux */
2142
2143 } else { // TODO: store val_c = inc*coface + cofbce*pipr
2144
2145 pfac = inc*coface + cofbce*pipr;
2146 *flux += iconvp*xcpp*(-imasac*(b_massflux*pi) + thetap*(pfac));
2147
2148 }
2149}
2150
2151/*----------------------------------------------------------------------------*/
2175/*----------------------------------------------------------------------------*/
2176
2177template <cs_lnum_t stride>
2178CS_F_HOST_DEVICE inline static void
2180 cs_real_t thetap,
2181 int imasac,
2182 int inc,
2183 int bc_type,
2184 int icvfli,
2185 const cs_real_t pi[stride],
2186 const cs_real_t pir[stride],
2187 const cs_real_t pipr[stride],
2188 const cs_real_t coface[stride],
2189 const cs_real_t cofbce[stride][stride],
2190 cs_real_t b_massflux,
2191 cs_real_t pfac[stride],
2192 cs_real_t flux[stride])
2193{
2194 cs_real_t flui, fluj;
2195
2196 /* Computed convective flux */
2197
2198 if (icvfli == 0) {
2199
2200 /* Remove decentering for coupled faces */
2201 if (bc_type == CS_COUPLED_FD) {
2202 flui = 0.0;
2203 fluj = b_massflux;
2204 }
2205 else {
2206 flui = 0.5*(b_massflux +fabs(b_massflux));
2207 fluj = 0.5*(b_massflux -fabs(b_massflux));
2208 }
2209
2210 for (int isou = 0; isou < stride; isou++)
2211 flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac[isou])
2212 - imasac*b_massflux*pi[isou]);
2213
2214 /* Imposed convective flux */
2215
2216 }
2217 else {
2218
2219 for (int isou = 0; isou < stride; isou++) {
2220 pfac[isou] = inc*coface[isou];
2221 for (int jsou = 0; jsou < stride; jsou++) {
2222 pfac[isou] += cofbce[isou][jsou]*pipr[jsou];
2223 }
2224 flux[isou] += iconvp*( thetap*pfac[isou]
2225 - imasac*b_massflux*pi[isou]);
2226 }
2227
2228 }
2229}
2230
2231/*----------------------------------------------------------------------------*/
2248/*----------------------------------------------------------------------------*/
2249
2250CS_F_HOST_DEVICE inline static void
2251cs_b_upwind_flux(const int iconvp,
2252 const cs_real_t thetap,
2253 const int imasac,
2254 const int bc_type,
2255 const cs_real_t pi,
2256 const cs_real_t pir,
2257 const cs_real_t val_f,
2258 const cs_real_t b_massflux,
2259 const cs_real_t xcpp,
2260 cs_real_t *flux)
2261{
2262 cs_real_t flui, fluj;
2263
2264 /* Remove decentering for coupled faces */
2265 if (bc_type == CS_COUPLED_FD) {
2266 flui = 0.0;
2267 fluj = b_massflux;
2268 } else {
2269 flui = 0.5*(b_massflux +fabs(b_massflux));
2270 fluj = 0.5*(b_massflux -fabs(b_massflux));
2271 }
2272
2273 *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*val_f) -imasac*( b_massflux*pi));
2274}
2275
2276/*----------------------------------------------------------------------------*/
2296/*----------------------------------------------------------------------------*/
2297
2298CS_F_HOST_DEVICE inline static void
2299cs_b_upwind_flux(const int iconvp,
2300 const cs_real_t thetap,
2301 const int imasac,
2302 const int inc,
2303 const int bc_type,
2304 const cs_real_t pi,
2305 const cs_real_t pir,
2306 const cs_real_t pipr,
2307 const cs_real_t coefap,
2308 const cs_real_t coefbp,
2309 const cs_real_t b_massflux,
2310 const cs_real_t xcpp,
2311 cs_real_t *flux)
2312{
2313 cs_real_t flui, fluj, pfac;
2314
2315 /* Remove decentering for coupled faces */
2316 if (bc_type == CS_COUPLED_FD) {
2317 flui = 0.0;
2318 fluj = b_massflux;
2319 } else {
2320 flui = 0.5*(b_massflux +fabs(b_massflux));
2321 fluj = 0.5*(b_massflux -fabs(b_massflux));
2322 }
2323
2324 pfac = inc*coefap + coefbp*pipr;
2325 *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*pfac) -imasac*( b_massflux*pi));
2326}
2327
2328/*----------------------------------------------------------------------------*/
2346/*----------------------------------------------------------------------------*/
2347
2348template <cs_lnum_t stride>
2349CS_F_HOST_DEVICE inline static void
2351 cs_real_t thetap,
2352 int imasac,
2353 int bc_type,
2354 const cs_real_t pi[stride],
2355 const cs_real_t pir[stride],
2356 const cs_real_t b_massflux,
2357 const cs_real_t pfac[stride],
2358 cs_real_t flux[stride])
2359{
2360 cs_real_t flui, fluj;
2361
2362 /* Remove decentering for coupled faces */
2363 if (bc_type == CS_COUPLED_FD) {
2364 flui = 0.0;
2365 fluj = b_massflux;
2366 } else {
2367 flui = 0.5*(b_massflux +fabs(b_massflux));
2368 fluj = 0.5*(b_massflux -fabs(b_massflux));
2369 }
2370
2371 for (int isou = 0; isou < stride; isou++)
2372 flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac[isou])
2373 - imasac*b_massflux*pi[isou]);
2374}
2375
2376/*----------------------------------------------------------------------------*/
2386/*----------------------------------------------------------------------------*/
2387
2388CS_F_HOST_DEVICE inline static void
2389cs_b_diff_flux(const int idiffp,
2390 const cs_real_t thetap,
2391 const cs_real_t pfacd,
2392 const cs_real_t b_visc,
2393 cs_real_t *flux)
2394{
2395 *flux += idiffp*thetap*b_visc*pfacd;
2396}
2397
2398/*----------------------------------------------------------------------------*/
2411/*----------------------------------------------------------------------------*/
2412
2413CS_F_HOST_DEVICE inline static void
2414cs_b_diff_flux(const int idiffp,
2415 const cs_real_t thetap,
2416 const int inc,
2417 const cs_real_t pipr,
2418 const cs_real_t cofafp,
2419 const cs_real_t cofbfp,
2420 const cs_real_t b_visc,
2421 cs_real_t *flux)
2422{
2423 cs_real_t pfacd = inc*cofafp + cofbfp*pipr;
2424 *flux += idiffp*thetap*b_visc*pfacd;
2425}
2426
2427/*----------------------------------------------------------------------------*/
2440/*----------------------------------------------------------------------------*/
2441
2442template <cs_lnum_t stride>
2443CS_F_HOST_DEVICE inline static void
2445 cs_real_t thetap,
2446 cs_real_t b_visc,
2447 const cs_real_t pfacd[stride],
2448 cs_real_t flux[stride])
2449{
2450 for (int isou = 0; isou < stride; isou++)
2451 flux[isou] += idiffp*thetap*b_visc*pfacd[isou];
2452}
2453
2454/*----------------------------------------------------------------------------*/
2465/*----------------------------------------------------------------------------*/
2466
2467template <typename T>
2468CS_F_HOST_DEVICE inline static void
2469cs_b_cd_unsteady(const T bldfrp,
2470 const cs_rreal_t diipb[3],
2471 const T gradi[3],
2472 const cs_real_t pi,
2473 cs_real_t *pip)
2474{
2475 T recoi;
2476
2478 gradi,
2479 bldfrp,
2480 &recoi);
2481
2482 *pip = pi + recoi;
2483}
2484
2485/*----------------------------------------------------------------------------*/
2499/*----------------------------------------------------------------------------*/
2500
2501template <cs_lnum_t stride, typename T>
2502CS_F_HOST_DEVICE inline static void
2504 const cs_rreal_t diipb[3],
2505 const T gradi[stride][3],
2506 const cs_real_t pi[stride],
2507 cs_real_t pip[stride])
2508{
2509 T recoi[stride];
2510
2511 cs_b_compute_quantities_strided<stride>(diipb,
2512 gradi,
2513 bldfrp,
2514 recoi);
2515
2516 for (int isou = 0; isou < stride; isou++)
2517 pip[isou] = pi[isou] + recoi[isou];
2518}
2519
2520/*----------------------------------------------------------------------------*/
2531/*----------------------------------------------------------------------------*/
2532
2533CS_F_HOST_DEVICE inline static void
2535 cs_real_t pi,
2536 cs_real_t pj,
2537 cs_real_t b_visc,
2538 cs_real_t *fluxi)
2539{
2540 *fluxi += idiffp*b_visc*(pi - pj);
2541}
2542
2543/*----------------------------------------------------------------------------*/
2557/*----------------------------------------------------------------------------*/
2558
2559template <cs_lnum_t stride>
2560CS_F_HOST_DEVICE inline static void
2562 const cs_real_t pi[stride],
2563 const cs_real_t pj[stride],
2564 cs_real_t b_visc,
2565 cs_real_t fluxi[stride])
2566{
2567 for (int k = 0; k < stride; k++)
2568 fluxi[k] += idiffp*b_visc*(pi[k] - pj[k]);
2569}
2570
2571/*============================================================================
2572 * Semi private functions
2573 *============================================================================*/
2574
2575/*----------------------------------------------------------------------------
2576 * Adust bounds if rc_clip_factor < or > 1, and check bounds
2577 * if required by verbosity level or postprocessing.
2578 *----------------------------------------------------------------------------*/
2579
2580template <typename T>
2581void
2583(
2585 const char *var_name,
2586 const cs_equation_param_t &eqp,
2587 bool face_gradient,
2588 int ircflb,
2589 const cs_mesh_t *m,
2590 const cs_mesh_quantities_t *fvq,
2591 const cs_real_t *restrict pvar,
2592 const T (*restrict grad)[3],
2593 const cs_real_t *df_limiter,
2594 cs_real_t (*restrict bounds)[2]
2595);
2596
2597/*----------------------------------------------------------------------------
2598 * Adust bounds if rc_clip_factor < or > 1, and check bounds
2599 *----------------------------------------------------------------------------*/
2600
2601template <cs_lnum_t stride, typename T>
2602void
2604(
2606 const char *var_name,
2607 const cs_equation_param_t &eqp,
2608 bool face_gradient,
2609 int ircflb,
2610 const cs_mesh_t *m,
2611 const cs_mesh_quantities_t *fvq,
2612 const T (*restrict grad)[stride][3],
2613 const cs_real_t *df_limiter,
2614 cs_real_t *restrict bounds
2615);
2616
2617/*----------------------------------------------------------------------------*
2618 * Add the explicit part of the convection/diffusion terms of a
2619 * standard transport equation of a scalar field \f$ \varia \f$.
2620 *
2621 * Implemented and documented in cs_convection_diffusion_steady.cpp.
2622 *----------------------------------------------------------------------------*/
2623
2624void
2626(
2627 const cs_field_t *f,
2628 const cs_equation_param_t &eqp,
2629 bool icvflb,
2630 int inc,
2631 const cs_real_t *restrict pvar,
2632 const cs_real_t *restrict pvara,
2633 const int icvfli[],
2634 const cs_field_bc_coeffs_t *bc_coeffs,
2635 const cs_real_t i_massflux[],
2636 const cs_real_t b_massflux[],
2637 const cs_real_t i_visc[],
2638 const cs_real_t b_visc[],
2639 const cs_real_t xcpp[],
2640 cs_real_t *restrict rhs
2641);
2642
2643/*----------------------------------------------------------------------------
2644 * \brief Add the explicit part of the convection/diffusion terms of a transport
2645 * equation of a vector or tensor field.
2646 *
2647 * Implemented and documented in cs_convection_diffusion_steady.cpp.
2648 *----------------------------------------------------------------------------*/
2649
2650template <cs_lnum_t stride>
2651void
2653(
2654 cs_field_t *f,
2655 const char *var_name,
2656 const cs_equation_param_t &eqp,
2657 int icvflb,
2658 int inc,
2659 cs_real_t (*pvar)[stride],
2660 const cs_real_t (*pvara)[stride],
2661 const int icvfli[],
2662 const cs_field_bc_coeffs_t *bc_coeffs,
2663 const cs_real_t i_massflux[],
2664 const cs_real_t b_massflux[],
2665 const cs_real_t i_visc[],
2666 const cs_real_t b_visc[],
2667 cs_real_t (*restrict grad)[stride][3],
2668 cs_real_t (*restrict rhs)[stride]
2669);
2670
2671/*----------------------------------------------------------------------------
2672 * Update face flux with convection contribution of a standard transport
2673 * equation of a scalar field \f$ \varia \f$.
2674 *
2675 * Implemented and documented in cs_convection_diffusion_steady.cpp.
2676 *----------------------------------------------------------------------------*/
2677
2678void
2680(
2681 const cs_field_t *f,
2682 const cs_equation_param_t eqp,
2683 int icvflb,
2684 int inc,
2685 cs_real_t *restrict pvar,
2686 const cs_real_t *restrict pvara,
2687 const int icvfli[],
2688 const cs_field_bc_coeffs_t *bc_coeffs,
2689 const cs_real_t i_massflux[],
2690 const cs_real_t b_massflux[],
2691 cs_real_t i_conv_flux[][2],
2692 cs_real_t b_conv_flux[]
2693);
2694
2695/*----------------------------------------------------------------------------
2696 * Add the explicit part of the convection/diffusion terms of a
2697 * standard transport equation of a scalar field \f$ \varia \f$.
2698 *----------------------------------------------------------------------------*/
2699
2700bool
2702 const cs_equation_param_t eqp,
2703 int icvflb,
2704 int inc,
2705 int imasac,
2706 const cs_real_t *restrict pvar,
2707 const int icvfli[],
2708 cs_field_bc_coeffs_t *bc_coeffs,
2709 const cs_real_t i_massflux[],
2710 const cs_real_t b_massflux[],
2711 const cs_real_t i_visc[],
2712 const cs_real_t b_visc[],
2713 const cs_real_t *c_weight,
2714 cs_real_t *rhs,
2715 cs_real_2_t i_flux[],
2716 cs_real_t b_flux[]);
2717
2718/*----------------------------------------------------------------------------
2719 * Add the explicit part of the convection/diffusion terms of a transport
2720 * equation of a scalar field \f$ \varia \f$ such as the temperature.
2721 *----------------------------------------------------------------------------*/
2722
2723bool
2725 const cs_equation_param_t eqp,
2726 int inc,
2727 int imasac,
2728 const cs_real_t * pvar,
2729 cs_field_bc_coeffs_t *bc_coeffs,
2730 const cs_real_t i_massflux[],
2731 const cs_real_t b_massflux[],
2732 const cs_real_t i_visc[],
2733 const cs_real_t b_visc[],
2734 const cs_real_t *c_weight,
2735 const cs_real_t xcpp[],
2736 cs_real_t *restrict rhs);
2737
2738/*----------------------------------------------------------------------------
2739 * Add the explicit part of the convection/diffusion terms of a transport
2740 * equation of a vector field \f$ \vect{\varia} \f$.
2741 *----------------------------------------------------------------------------*/
2742
2743bool
2745 int f_id,
2746 const cs_equation_param_t eqp,
2747 int icvflb,
2748 int inc,
2749 int ivisep,
2750 int imasac,
2751 cs_real_3_t *restrict pvar,
2752 const cs_real_3_t *restrict pvara,
2753 const int icvfli[],
2754 cs_field_bc_coeffs_t *bc_coeffs,
2755 const cs_real_t i_massflux[],
2756 const cs_real_t b_massflux[],
2757 const cs_real_t i_visc[],
2758 const cs_real_t b_visc[],
2759 const cs_real_t i_secvis[],
2760 const cs_real_t b_secvis[],
2761 cs_real_3_t *restrict i_pvar,
2762 cs_real_3_t *restrict b_pvar,
2763 cs_real_3_t *restrict rhs);
2764
2765/*----------------------------------------------------------------------------
2766 * Add the explicit part of the convection/diffusion terms of a transport
2767 * equation of a tensor field \f$ \tens{\varia} \f$.
2768 *----------------------------------------------------------------------------*/
2769
2770bool
2772 int f_id,
2773 const cs_equation_param_t eqp,
2774 int icvflb,
2775 int inc,
2776 int imasac,
2777 cs_real_6_t *restrict pvar,
2778 const cs_real_6_t *restrict pvara,
2779 cs_field_bc_coeffs_t *bc_coeffs,
2780 const cs_real_t i_massflux[],
2781 const cs_real_t b_massflux[],
2782 const cs_real_t i_visc[],
2783 const cs_real_t b_visc[],
2784 cs_real_6_t *restrict rhs);
2785
2786/*----------------------------------------------------------------------------*/
2787
2788#endif /* CS_CONVECTION_DIFFUSION_PRIV_H */
Definition: cs_dispatch.h:2288
Field boundary condition descriptor (for variables)
Definition: cs_field.h:107
Field descriptor.
Definition: cs_field.h:275
cs_nvd_type_t
Definition: cs_convection_diffusion.h:55
@ CS_NVD_SUPERBEE
Definition: cs_convection_diffusion.h:60
@ CS_NVD_SMART
Definition: cs_convection_diffusion.h:58
@ CS_NVD_CUBISTA
Definition: cs_convection_diffusion.h:59
@ CS_NVD_STOIC
Definition: cs_convection_diffusion.h:64
@ CS_NVD_WASEB
Definition: cs_convection_diffusion.h:66
@ CS_NVD_CLAM
Definition: cs_convection_diffusion.h:63
@ CS_NVD_OSHER
Definition: cs_convection_diffusion.h:65
@ CS_NVD_VOF_CICSAM
Definition: cs_convection_diffusion.h:68
@ CS_NVD_GAMMA
Definition: cs_convection_diffusion.h:57
@ CS_NVD_MINMOD
Definition: cs_convection_diffusion.h:62
@ CS_NVD_VOF_HRIC
Definition: cs_convection_diffusion.h:67
@ CS_NVD_MUSCL
Definition: cs_convection_diffusion.h:61
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:2534
static CS_F_HOST_DEVICE void cs_i_compute_quantities(const T bldfrp, const cs_rreal_3_t diipf, const cs_rreal_3_t djjpf, const T gradi[3], const T gradj[3], const cs_real_t pi, const cs_real_t pj, T *recoi, T *recoj, cs_real_t *pip, cs_real_t *pjp)
Reconstruct values in I' and J'.
Definition: cs_convection_diffusion_priv.h:545
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:1061
void cs_convection_diffusion_steady_scalar(const cs_field_t *f, const cs_equation_param_t &eqp, bool icvflb, int inc, const cs_real_t *restrict pvar, const cs_real_t *restrict pvara, const int icvfli[], const cs_field_bc_coeffs_t *bc_coeffs, const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], const cs_real_t xcpp[], cs_real_t *restrict rhs)
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:795
bool cs_convection_diffusion_vector_v9(int idtvar, int f_id, const cs_equation_param_t eqp, int icvflb, int inc, int ivisep, int imasac, cs_real_3_t *restrict pvar, const cs_real_3_t *restrict pvara, const int icvfli[], cs_field_bc_coeffs_t *bc_coeffs, const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], const cs_real_t i_secvis[], const cs_real_t b_secvis[], cs_real_3_t *restrict i_pvar, cs_real_3_t *restrict b_pvar, cs_real_3_t *restrict rhs)
Add the explicit part of the convection/diffusion terms of a transport equation of a vector field .
Definition: cs_convection_diffusion_v9.cpp:2284
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:1853
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:1214
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 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:1715
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:2069
void cs_convection_diffusion_adjust_and_check_bounds_scalar(cs_dispatch_context &ctx, const char *var_name, const cs_equation_param_t &eqp, bool face_gradient, int ircflb, const cs_mesh_t *m, const cs_mesh_quantities_t *fvq, const cs_real_t *restrict pvar, const T(*restrict grad)[3], const cs_real_t *df_limiter, cs_real_t(*restrict bounds)[2])
Adust bounds if rc_clip_factor < or > 1, and check bounds if required by verbosity level or postproce...
Definition: cs_convection_diffusion.cpp:4646
bool cs_convection_diffusion_thermal_v9(const cs_field_t *f, const cs_equation_param_t eqp, int inc, int imasac, const cs_real_t *pvar, cs_field_bc_coeffs_t *bc_coeffs, const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], const cs_real_t *c_weight, const cs_real_t xcpp[], cs_real_t *restrict rhs)
Add the explicit part of the convection/diffusion terms of a transport equation of a scalar field su...
Definition: cs_convection_diffusion_v9.cpp:2187
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:1151
void cs_convection_diffusion_adjust_and_check_bounds_strided(cs_dispatch_context &ctx, const char *var_name, const cs_equation_param_t &eqp, bool face_gradient, int ircflb, const cs_mesh_t *m, const cs_mesh_quantities_t *fvq, const T(*restrict grad)[stride][3], const cs_real_t *df_limiter, cs_real_t *restrict bounds)
Adust bounds if rc_clip_factor < or > 1, and check bounds if required by verbosity level or postproce...
Definition: cs_convection_diffusion.cpp:4734
void cs_face_convection_steady_scalar(const cs_field_t *f, const cs_equation_param_t eqp, int icvflb, int inc, cs_real_t *restrict pvar, const cs_real_t *restrict pvara, const int icvfli[], const cs_field_bc_coeffs_t *bc_coeffs, const cs_real_t i_massflux[], const cs_real_t b_massflux[], cs_real_t i_conv_flux[][2], cs_real_t b_conv_flux[])
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:69
bool cs_convection_diffusion_tensor_v9(int idtvar, int f_id, const cs_equation_param_t eqp, int icvflb, int inc, int imasac, cs_real_6_t *restrict pvar, const cs_real_6_t *restrict pvara, cs_field_bc_coeffs_t *bc_coeffs, const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t *restrict rhs)
Add the explicit part of the convection/diffusion terms of a transport equation of a tensor field .
Definition: cs_convection_diffusion_v9.cpp:2543
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:696
static CS_F_HOST_DEVICE void cs_b_upwind_flux(const int iconvp, const cs_real_t thetap, const int imasac, const int bc_type, const cs_real_t pi, const cs_real_t pir, const cs_real_t val_f, 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:2251
static CS_F_HOST_DEVICE void cs_i_compute_quantities_strided(const T bldfrp, const cs_rreal_t diipf[3], const cs_rreal_t djjpf[3], const T gradi[stride][3], const T gradj[stride][3], const cs_real_t pi[stride], const cs_real_t pj[stride], T recoi[stride], 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:591
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:1517
static CS_F_HOST_DEVICE void cs_b_upwind_flux_strided(int iconvp, cs_real_t thetap, int imasac, int bc_type, const cs_real_t pi[stride], const cs_real_t pir[stride], const cs_real_t b_massflux, const 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:2350
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:2561
static CS_F_HOST_DEVICE void cs_b_cd_unsteady_strided(T bldfrp, const cs_rreal_t diipb[3], const 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:2503
static CS_F_HOST_DEVICE void cs_clip_quantity_strided(const cs_real_t &bounds, const cs_real_t pi[stride], cs_real_t pip[stride])
Clip reconstructed value in strided case.
Definition: cs_convection_diffusion_priv.h:660
static CS_F_HOST_DEVICE void cs_b_diff_flux(const int idiffp, const cs_real_t thetap, const cs_real_t pfacd, const cs_real_t b_visc, cs_real_t *flux)
Add diffusive flux to flux at boundary face.
Definition: cs_convection_diffusion_priv.h:2389
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:2040
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:771
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 T gradi[3], const T gradj[3], const T grdpai[3], const T grdpaj[3], const cs_real_t i_massflux, T *testij, T *tesqck)
Compute slope test criteria at internal face between cell i and j.
Definition: cs_convection_diffusion_priv.h:410
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 coface, cs_real_t cofbce, cs_real_t b_massflux, cs_real_t xcpp, cs_real_t val_f, cs_real_t *flux)
Add convective flux (substracting the mass accumulation from it) to flux at boundary face....
Definition: cs_convection_diffusion_priv.h:2108
static CS_F_HOST_DEVICE void cs_clip_quantity(const cs_real_t bounds[2], cs_real_t &pip)
Clip reconstructed value.
Definition: cs_convection_diffusion_priv.h:634
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 T recoi[stride], const 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:738
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 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:2179
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:1097
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:905
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:985
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 T gradi[stride][3], const T gradj[stride][3], const T gradsti[stride][3], const T gradstj[stride][3], const cs_real_t i_massflux, T *testij, T *tesqck)
Compute slope test criteria at internal face between cell i and j.
Definition: cs_convection_diffusion_priv.h:477
static CS_F_HOST_DEVICE void cs_b_compute_quantities_strided(const cs_rreal_t diipb[3], const T gradi[stride][3], const T bldfrp, T recoi[stride])
Reconstruct values in I' at boundary cell i.
Definition: cs_convection_diffusion_priv.h:2014
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:1027
bool cs_convection_diffusion_scalar_v9(const cs_field_t *f, const cs_equation_param_t eqp, int icvflb, int inc, int imasac, const cs_real_t *restrict pvar, const int icvfli[], cs_field_bc_coeffs_t *bc_coeffs, const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], const cs_real_t *c_weight, cs_real_t *rhs, cs_real_2_t i_flux[], cs_real_t b_flux[])
Add the explicit part of the convection/diffusion terms of a standard transport equation of a scalar ...
Definition: cs_convection_diffusion_v9.cpp:2081
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 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:276
static CS_F_HOST_DEVICE void cs_b_compute_quantities(const cs_rreal_t diipb[3], const T gradi[3], const T bldfrp, T *recoi)
Reconstruct values in I' at boundary cell i.
Definition: cs_convection_diffusion_priv.h:1988
static CS_F_HOST_DEVICE void cs_b_diff_flux_strided(int idiffp, cs_real_t thetap, cs_real_t b_visc, const cs_real_t pfacd[stride], cs_real_t flux[stride])
Add diffusive flux to flux at boundary face.
Definition: cs_convection_diffusion_priv.h:2444
static CS_F_HOST_DEVICE void cs_b_cd_unsteady(const T bldfrp, const cs_rreal_t diipb[3], const 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:2469
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:1675
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:937
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 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:818
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:91
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:881
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:1377
void cs_convection_diffusion_steady_strided(cs_field_t *f, const char *var_name, const cs_equation_param_t &eqp, int icvflb, int inc, cs_real_t(*pvar)[stride], const cs_real_t(*pvara)[stride], const int icvfli[], const cs_field_bc_coeffs_t *bc_coeffs, const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_t(*restrict grad)[stride][3], cs_real_t(*restrict rhs)[stride])
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 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:850
#define restrict
Definition: cs_defs.h:148
#define CS_F_HOST_DEVICE
Definition: cs_defs.h:555
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
cs_rreal_t cs_rreal_3_t[3]
Definition: cs_defs.h:378
cs_nreal_t cs_nreal_3_t[3]
Definition: cs_defs.h:375
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:349
cs_real_t cs_real_2_t[2]
vector of 2 floating-point values
Definition: cs_defs.h:348
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:351
#define CS_UNUSED(x)
Definition: cs_defs.h:518
double cs_rreal_t
Definition: cs_defs.h:338
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
double cs_nreal_t
Definition: cs_defs.h:336
@ p
Definition: cs_field_pointer.h:63
@ k
Definition: cs_field_pointer.h:68
void cs_halo_sync_var(const cs_halo_t *halo, cs_halo_type_t sync_mode, cs_real_t var[])
Definition: cs_halo.cpp:2372
cs_halo_type_t
Definition: cs_halo.h:53
CS_F_HOST_DEVICE cs_real_t cs_math_3_dot_product(const T u[3], const U v[3])
Compute the dot product of two vectors of 3 real values.
Definition: cs_math.h:194
static constexpr cs_real_t cs_math_epzero
Definition: cs_math.h:102
CS_F_HOST_DEVICE T cs_math_sq(T x)
Compute the square of a real value.
Definition: cs_math.h:149
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:970
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:825
CS_F_HOST_DEVICE T cs_math_3_norm(const T v[3])
Compute the euclidean norm of a vector of dimension 3.
Definition: cs_math.h:240
@ CS_COUPLED_FD
Definition: cs_parameters.h:99
double precision pi
value with 16 digits
Definition: cstnum.f90:42
integer(c_int), pointer, save idtvar
option for a variable time step
Definition: optcal.f90:70
CS_F_HOST_DEVICE T max(const T a, const T b)
Definition: cs_defs.h:735
CS_F_HOST_DEVICE T min(const T a, const T b)
Definition: cs_defs.h:712
CS_F_HOST_DEVICE T abs(const T a)
Definition: cs_defs.h:690
Set of parameters to handle an unsteady convection-diffusion-reaction equation with term sources.
Definition: cs_equation_param.h:190
Definition: cs_mesh_quantities.h:88
Definition: cs_mesh.h:85
cs_halo_t * halo
Definition: cs_mesh.h:156