7.0
general documentation
cs_convection_diffusion.h
Go to the documentation of this file.
1 #ifndef __CS_CONVECTION_DIFFUSION_H__
2 #define __CS_CONVECTION_DIFFUSION_H__
3 
4 /*============================================================================
5  * Convection-diffusion operators.
6  *============================================================================*/
7 
8 /*
9  This file is part of Code_Saturne, a general-purpose CFD tool.
10 
11  Copyright (C) 1998-2021 EDF S.A.
12 
13  This program is free software; you can redistribute it and/or modify it under
14  the terms of the GNU General Public License as published by the Free Software
15  Foundation; either version 2 of the License, or (at your option) any later
16  version.
17 
18  This program is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21  details.
22 
23  You should have received a copy of the GNU General Public License along with
24  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25  Street, Fifth Floor, Boston, MA 02110-1301, USA.
26 */
27 
28 /*----------------------------------------------------------------------------*/
29 
30 #include "cs_defs.h"
31 
32 /*----------------------------------------------------------------------------
33  * Local headers
34  *----------------------------------------------------------------------------*/
35 
36 #include "cs_base.h"
37 #include "cs_halo.h"
38 #include "cs_math.h"
39 #include "cs_mesh_quantities.h"
40 #include "cs_parameters.h"
41 #include "cs_field_pointer.h"
42 
43 /*----------------------------------------------------------------------------*/
44 
45 BEGIN_C_DECLS
46 
47 /*=============================================================================
48  * Macro definitions
49  *============================================================================*/
50 
51 /*============================================================================
52  * Type definition
53  *============================================================================*/
54 
55 /*----------------------------------------------------------------------------
56  * NVD/TVD Advection Scheme
57  *----------------------------------------------------------------------------*/
58 
59 typedef enum {
60 
61  CS_NVD_GAMMA = 0, /* GAMMA */
62  CS_NVD_SMART = 1, /* SMART */
63  CS_NVD_CUBISTA = 2, /* CUBISTA */
64  CS_NVD_SUPERBEE = 3, /* SUPERBEE */
65  CS_NVD_MUSCL = 4, /* MUSCL */
66  CS_NVD_MINMOD = 5, /* MINMOD */
67  CS_NVD_CLAM = 6, /* CLAM */
68  CS_NVD_STOIC = 7, /* STOIC */
69  CS_NVD_OSHER = 8, /* OSHER */
70  CS_NVD_WASEB = 9, /* WASEB */
71  CS_NVD_VOF_HRIC = 10, /* M-HRIC for VOF */
72  CS_NVD_VOF_CICSAM = 11, /* M-CICSAM for VOF */
73  CS_NVD_VOF_STACS = 12, /* STACS for VOF */
74  CS_NVD_N_TYPES = 13 /* number of NVD schemes */
75 
76 } cs_nvd_type_t;
77 
78 /*============================================================================
79  * Global variables
80  *============================================================================*/
81 
82 /*============================================================================
83  * Public inlined function
84  *============================================================================*/
85 
86 /*----------------------------------------------------------------------------
87  * Synchronize halos for scalar variables.
88  *
89  * parameters:
90  * m <-- pointer to associated mesh structure
91  * tr_dim <-- 0 if pvar does not match a tensor
92  * or there is no periodicity of rotation
93  * 2 for Reynolds stress
94  * pvar <-> variable
95  *----------------------------------------------------------------------------*/
96 
97 inline static void
98 cs_sync_scalar_halo(const cs_mesh_t *m,
99  int tr_dim,
100  cs_real_t pvar[])
101 {
102  if (m->halo != NULL) {
103  if (tr_dim > 0)
104  cs_halo_sync_component(m->halo, CS_HALO_STANDARD, CS_HALO_ROTATION_IGNORE,
105  pvar);
106  else
107  cs_halo_sync_var(m->halo, CS_HALO_STANDARD, pvar);
108  }
109 }
110 
111 /*----------------------------------------------------------------------------
112  * Return pointer to slope test indicator field values if active.
113  *
114  * parameters:
115  * f_id <-- field id (or -1)
116  * var_cal_opt <-- variable calculation options
117  *
118  * return:
119  * pointer to local values array, or NULL;
120  *----------------------------------------------------------------------------*/
121 
122 inline static cs_real_t *
123 cs_get_v_slope_test(int f_id,
124  const cs_var_cal_opt_t var_cal_opt)
125 {
126  const int iconvp = var_cal_opt.iconv;
127  const int isstpp = var_cal_opt.isstpc;
128  const double blencp = var_cal_opt.blencv;
129 
130  cs_real_t *v_slope_test = NULL;
131 
132  if (f_id > -1 && iconvp > 0 && blencp > 0. && isstpp == 0) {
133 
134  static int _k_slope_test_f_id = -1;
135 
136  cs_field_t *f = cs_field_by_id(f_id);
137 
138  int f_track_slope_test_id = -1;
139 
140  if (_k_slope_test_f_id < 0)
141  _k_slope_test_f_id = cs_field_key_id_try("slope_test_upwind_id");
142  if (_k_slope_test_f_id > -1 && isstpp == 0)
143  f_track_slope_test_id = cs_field_get_key_int(f, _k_slope_test_f_id);
144 
145  if (f_track_slope_test_id > -1)
146  v_slope_test = (cs_field_by_id(f_track_slope_test_id))->val;
147 
148  if (v_slope_test != NULL) {
149  const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
150 # pragma omp parallel for
151  for (cs_lnum_t cell_id = 0; cell_id < n_cells_ext; cell_id++)
152  v_slope_test[cell_id] = 0.;
153  }
154 
155  }
156 
157  return v_slope_test;
158 }
159 
160 /*----------------------------------------------------------------------------
161  * Compute the local cell Courant number as the maximum of all cell face based
162  * Courant number at each cell.
163  *
164  * parameters:
165  * f_id <-- field id (or -1)
166  * courant --> cell Courant number
167  */
168 /*----------------------------------------------------------------------------*/
169 
170 inline static void
171 cs_cell_courant_number(const int f_id,
172  cs_real_t *courant)
173 {
174  const cs_mesh_t *m = cs_glob_mesh;
175  const cs_mesh_quantities_t *fvq = cs_glob_mesh_quantities;
176 
177  const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
178  const int n_i_groups = m->i_face_numbering->n_groups;
179  const int n_i_threads = m->i_face_numbering->n_threads;
180  const int n_b_groups = m->b_face_numbering->n_groups;
181  const int n_b_threads = m->b_face_numbering->n_threads;
182  const cs_lnum_t *restrict i_group_index = m->i_face_numbering->group_index;
183  const cs_lnum_t *restrict b_group_index = m->b_face_numbering->group_index;
184 
185  const cs_lnum_2_t *restrict i_face_cells
186  = (const cs_lnum_2_t *restrict)m->i_face_cells;
187  const cs_lnum_t *restrict b_face_cells
188  = (const cs_lnum_t *restrict)m->b_face_cells;
189 
190  const cs_real_t *restrict vol
191  = (cs_real_t *restrict)fvq->cell_vol;
192 
193  cs_field_t *f = cs_field_by_id(f_id);
194  const int kimasf = cs_field_key_id("inner_mass_flux_id");
195  const int kbmasf = cs_field_key_id("boundary_mass_flux_id");
196  const cs_real_t *restrict i_massflux
197  = cs_field_by_id( cs_field_get_key_int(f, kimasf) )->val;
198  const cs_real_t *restrict b_massflux
199  = cs_field_by_id( cs_field_get_key_int(f, kbmasf) )->val;
200 
201  const cs_real_t *restrict dt
202  = (const cs_real_t *restrict)CS_F_(dt)->val;
203 
204  /* Initialisation */
205 
206 # pragma omp parallel for
207  for (cs_lnum_t ii = 0; ii < n_cells_ext; ii++) {
208  courant[ii] = 0.;
209  }
210 
211  /* ---> Contribution from interior faces */
212 
213  cs_real_t cnt;
214 
215  for (int g_id = 0; g_id < n_i_groups; g_id++) {
216 # pragma omp parallel for
217  for (int t_id = 0; t_id < n_i_threads; t_id++) {
218  for (cs_lnum_t face_id = i_group_index[(t_id*n_i_groups + g_id)*2];
219  face_id < i_group_index[(t_id*n_i_groups + g_id)*2 + 1];
220  face_id++) {
221  cs_lnum_t ii = i_face_cells[face_id][0];
222  cs_lnum_t jj = i_face_cells[face_id][1];
223 
224  cnt = CS_ABS(i_massflux[face_id])*dt[ii]/vol[ii];
225  courant[ii] = CS_MAX(courant[ii], cnt);
226 
227  cnt = CS_ABS(i_massflux[face_id])*dt[jj]/vol[jj];
228  courant[jj] = CS_MAX(courant[jj], cnt);
229  }
230  }
231  }
232 
233  /* ---> Contribution from boundary faces */
234 
235  for (int g_id = 0; g_id < n_b_groups; g_id++) {
236 # pragma omp parallel for
237  for (int t_id = 0; t_id < n_b_threads; t_id++) {
238  for (cs_lnum_t face_id = b_group_index[(t_id*n_b_groups + g_id)*2];
239  face_id < b_group_index[(t_id*n_b_groups + g_id)*2 + 1];
240  face_id++) {
241  cs_lnum_t ii = b_face_cells[face_id];
242 
243  cnt = CS_ABS(b_massflux[face_id])*dt[ii]/vol[ii];
244  courant[ii] = CS_MAX(courant[ii], cnt);
245  }
246  }
247  }
248 }
249 
250 /*----------------------------------------------------------------------------*/
261 /*----------------------------------------------------------------------------*/
262 
263 inline static cs_real_t
264 cs_nvd_scheme_scalar(const cs_nvd_type_t limiter,
265  const cs_real_t nvf_p_c,
266  const cs_real_t nvf_r_f,
267  const cs_real_t nvf_r_c)
268 {
269  cs_real_t nvf_p_f;
270 
271  cs_real_t beta_m, rfc, r1f, r1, r2, r3, b1, b2;
272 
273  switch (limiter) {
274  case CS_NVD_GAMMA: /* Gamma scheme */
275  beta_m = 0.1; /* in [0.1, 0.5] */
276  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
277 
278  if (nvf_p_c < beta_m) {
279  nvf_p_f = nvf_p_c*(1.+rfc*(1.-nvf_p_c)/beta_m);
280  } else {
281  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
282 
283  nvf_p_f = r1f*nvf_p_c+rfc;
284  }
285 
286  break;
287 
288  case CS_NVD_SMART: /* SMART scheme */
289  if (nvf_p_c < (nvf_r_c/3.)) {
290  r1 = nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
291  r2 = nvf_r_c*(1.-nvf_r_c);
292 
293  nvf_p_f = nvf_p_c*r1/r2;
294  } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
295  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
296  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
297 
298  nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c + rfc);
299  } else {
300  nvf_p_f = 1.;
301  }
302 
303  break;
304 
305  case CS_NVD_CUBISTA: /* CUBISTA scheme */
306  if (nvf_p_c < (3.*nvf_r_c/4.)) {
307  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
308 
309  nvf_p_f = nvf_r_f*(1.+rfc/3.)*nvf_p_c/nvf_r_c;
310  } else if (nvf_p_c <= (nvf_r_c*(1.+2.*(nvf_r_f-nvf_r_c))/(2.*nvf_r_f-nvf_r_c))) {
311  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
312  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
313 
314  nvf_p_f = nvf_r_f*(r1f*nvf_p_c/nvf_r_c+rfc);
315  } else {
316  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
317 
318  nvf_p_f = 1.-.5*r1f*(1.-nvf_p_c);
319  }
320 
321  break;
322 
323  case CS_NVD_SUPERBEE: /* SuperBee scheme */
324  if (nvf_p_c < (nvf_r_c/(2.-nvf_r_c))) {
325  nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
326  } else if (nvf_p_c < nvf_r_c) {
327  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
328  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
329 
330  nvf_p_f = r1f*nvf_p_c+rfc;
331  } else if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
332  nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
333  } else {
334  nvf_p_f = 1.;
335  }
336 
337  break;
338 
339  case CS_NVD_MUSCL: /* MUSCL scheme */
340  if (nvf_p_c < (.5*nvf_r_c)) {
341  nvf_p_f = (2.*nvf_r_f-nvf_r_c)*nvf_p_c/nvf_r_c;
342  } else if (nvf_p_c < (1.+nvf_r_c-nvf_r_f)) {
343  nvf_p_f = nvf_p_c+nvf_r_f-nvf_r_c;
344  } else {
345  nvf_p_f = 1.;
346  }
347 
348  break;
349 
350  case CS_NVD_MINMOD: /* MINMOD scheme */
351  if (nvf_p_c < nvf_r_c) {
352  nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
353  } else {
354  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
355  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
356 
357  nvf_p_f = r1f*nvf_p_c+rfc;
358  }
359 
360  break;
361 
362  case CS_NVD_CLAM: /* CLAM scheme */
363  r1 = nvf_r_c*nvf_r_c-nvf_r_f;
364  r2 = nvf_r_c*(nvf_r_c-1.);
365  r3 = nvf_r_f-nvf_r_c;
366 
367  nvf_p_f = nvf_p_c*(r1+r3*nvf_p_c)/r2;
368 
369  break;
370 
371  case CS_NVD_STOIC: /* STOIC scheme */
372  b1 = (nvf_r_c-nvf_r_f)*nvf_r_c;
373  b2 = nvf_r_c+nvf_r_f+2.*nvf_r_f*nvf_r_f-4.*nvf_r_f*nvf_r_c;
374 
375  if (nvf_p_c < (b1/b2)) {
376  r1 = -nvf_r_f*(1.-3.*nvf_r_c+2.*nvf_r_f);
377  r2 = nvf_r_c*(nvf_r_c-1.);
378 
379  nvf_p_f = nvf_p_c*r1/r2;
380  } else if (nvf_p_c < nvf_r_c) {
381  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
382  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
383 
384  nvf_p_f = r1f*nvf_p_c+rfc;
385  } else if (nvf_p_c < (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
386  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
387  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
388 
389  nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
390  } else {
391  nvf_p_f = 1.;
392  }
393 
394  break;
395 
396  case CS_NVD_OSHER: /* OSHER scheme */
397  if (nvf_p_c < (nvf_r_c/nvf_r_f)) {
398  nvf_p_f = nvf_r_f*nvf_p_c/nvf_r_c;
399  } else {
400  nvf_p_f = 1.;
401  }
402 
403  break;
404 
405  case CS_NVD_WASEB: /* WASEB scheme */
406  r1 = nvf_r_c*nvf_r_f*(nvf_r_f-nvf_r_c);
407  r2 = 2.*nvf_r_c*(1.-nvf_r_c)-nvf_r_f*(1.-nvf_r_f);
408 
409  if (nvf_p_c < (r1/r2)) {
410  nvf_p_f = 2.*nvf_p_c;
411  } else if (nvf_p_c <= (nvf_r_c*(1.+nvf_r_f-nvf_r_c)/nvf_r_f)) {
412  rfc = (nvf_r_f-nvf_r_c)/(1.-nvf_r_c);
413  r1f = (1.-nvf_r_f)/(1.-nvf_r_c);
414 
415  nvf_p_f = nvf_r_f*(nvf_p_c*r1f/nvf_r_c+rfc);
416  } else {
417  nvf_p_f = 1.;
418  }
419 
420  break;
421 
422  default: /* Upwinding */
423  nvf_p_f = nvf_p_c;
424  break;
425  }
426 
427  return nvf_p_f;
428 }
429 
430 /*----------------------------------------------------------------------------*/
446 /*----------------------------------------------------------------------------*/
447 
448 inline static cs_real_t
449 cs_nvd_vof_scheme_scalar(const cs_nvd_type_t limiter,
450  const cs_real_3_t i_face_normal,
451  const cs_real_t nvf_p_c,
452  const cs_real_t nvf_r_f,
453  const cs_real_t nvf_r_c,
454  const cs_real_3_t gradv_c,
455  const cs_real_t c_courant)
456 {
457  assert(limiter >= CS_NVD_VOF_HRIC);
458 
459  cs_real_t nvf_p_f;
460  cs_real_t blend, high_order, low_order, ratio;
461 
462  /* Compute gradient angle indicator */
463  cs_real_t dotp = CS_ABS(cs_math_3_dot_product(gradv_c, i_face_normal));
464  cs_real_t sgrad = cs_math_3_norm(gradv_c);
465  cs_real_t snorm = cs_math_3_norm(i_face_normal);
466  cs_real_t denom = snorm*sgrad;
467 
468  if (limiter == CS_NVD_VOF_HRIC) { /* M-HRIC scheme */
469  /* High order scheme : Bounded Downwind */
470  if (nvf_p_c <= .5) {
471  high_order = 2.*nvf_p_c;
472  } else {
473  high_order = 1.;
474  }
475 
476  /* Low order scheme : MUSCL */
477  low_order = cs_nvd_scheme_scalar(CS_NVD_MUSCL,
478  nvf_p_c,
479  nvf_r_f,
480  nvf_r_c);
481 
482  /* Compute the blending factor */
483  if (denom <= (cs_math_epzero*dotp)) {
484  blend = 1.;
485  } else {
486  ratio = dotp/denom;
487  blend = CS_MIN(1., pow(ratio, .5));
488  }
489 
490  /* Blending */
491  nvf_p_f = blend*high_order + (1.-blend)*low_order;
492 
493  /* Extra blending due to the cell Courant number */
494  if (c_courant < .7 && c_courant > .3) {
495  nvf_p_f = nvf_p_f + (nvf_p_f - low_order)*(.7 - c_courant )/.4;
496  } else if (c_courant >= .7) {
497  nvf_p_f = low_order;
498  }
499  } else if (limiter == CS_NVD_VOF_CICSAM) { /* M-CICSAM scheme */
500  /* High order scheme : HYPER-C + SUPERBEE */
501  if (c_courant <= .3) {
502  high_order = CS_MIN(1., nvf_p_c/(c_courant+cs_math_epzero));
503  } else if (c_courant <= .6) {
504  high_order = CS_MIN(1., nvf_p_c/.3);
505  } else if (c_courant <= .7) {
506  cs_real_t superbee = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
507  nvf_p_c,
508  nvf_r_f,
509  nvf_r_c);
510  high_order = 10.*( (.7-c_courant)*CS_MIN(1., nvf_p_c/.3)
511  + (c_courant-.6)*superbee);
512  }
513  else {
514  high_order = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
515  nvf_p_c,
516  nvf_r_f,
517  nvf_r_c);
518  }
519 
520  /* Low order scheme : MUSCL */
521  low_order = cs_nvd_scheme_scalar(CS_NVD_MUSCL,
522  nvf_p_c,
523  nvf_r_f,
524  nvf_r_c);
525 
526  /* Compute the blending factor */
527  if (denom <= (cs_math_epzero*dotp)) {
528  blend = 1.;
529  } else {
530  ratio = dotp/denom;
531  blend = CS_MIN(1., pow(ratio, 2.));
532  }
533 
534  /* Blending */
535  nvf_p_f = blend*high_order + (1.-blend)*low_order;
536  } else { /* STACS scheme */
537  /* High order scheme : SUPERBEE */
538  high_order = cs_nvd_scheme_scalar(CS_NVD_SUPERBEE,
539  nvf_p_c,
540  nvf_r_f,
541  nvf_r_c);
542 
543  /* Low order scheme : STOIC */
544  low_order = cs_nvd_scheme_scalar(CS_NVD_STOIC,
545  nvf_p_c,
546  nvf_r_f,
547  nvf_r_c);
548 
549  /* Compute the blending factor */
550  if (denom <= (cs_math_epzero*dotp)) {
551  blend = 1.;
552  } else {
553  ratio = dotp/denom;
554  blend = CS_MIN(1., pow(ratio, 4.));
555  }
556 
557  /* Blending */
558  nvf_p_f = blend*high_order + (1.-blend)*low_order;
559  }
560 
561  return nvf_p_f;
562 }
563 
564 /*----------------------------------------------------------------------------*/
581 /*----------------------------------------------------------------------------*/
582 
583 inline static void
584 cs_slope_test(const cs_real_t pi,
585  const cs_real_t pj,
586  const cs_real_t distf,
587  const cs_real_t srfan,
588  const cs_real_t i_face_normal[3],
589  const cs_real_t gradi[3],
590  const cs_real_t gradj[3],
591  const cs_real_t grdpai[3],
592  const cs_real_t grdpaj[3],
593  const cs_real_t i_massflux,
594  cs_real_t *testij,
595  cs_real_t *tesqck)
596 {
597  cs_real_t testi, testj;
598  cs_real_t dcc, ddi, ddj;
599 
600  /* Slope test */
601 
602  testi = grdpai[0]*i_face_normal[0]
603  + grdpai[1]*i_face_normal[1]
604  + grdpai[2]*i_face_normal[2];
605  testj = grdpaj[0]*i_face_normal[0]
606  + grdpaj[1]*i_face_normal[1]
607  + grdpaj[2]*i_face_normal[2];
608 
609  *testij = grdpai[0]*grdpaj[0]
610  + grdpai[1]*grdpaj[1]
611  + grdpai[2]*grdpaj[2];
612 
613  if (i_massflux>0.) {
614  dcc = gradi[0]*i_face_normal[0]
615  + gradi[1]*i_face_normal[1]
616  + gradi[2]*i_face_normal[2];
617  ddi = testi;
618  ddj = (pj-pi)/distf *srfan;
619  } else {
620  dcc = gradj[0]*i_face_normal[0]
621  + gradj[1]*i_face_normal[1]
622  + gradj[2]*i_face_normal[2];
623  ddi = (pj-pi)/distf *srfan;
624  ddj = testj;
625  }
626 
627  *tesqck = cs_math_sq(dcc) - cs_math_sq(ddi-ddj);
628 }
629 
630 /*----------------------------------------------------------------------------*/
647 /*----------------------------------------------------------------------------*/
648 
649 inline static void
650 cs_slope_test_vector(const cs_real_t pi[3],
651  const cs_real_t pj[3],
652  const cs_real_t distf,
653  const cs_real_t srfan,
654  const cs_real_t i_face_normal[3],
655  const cs_real_t gradi[3][3],
656  const cs_real_t gradj[3][3],
657  const cs_real_t gradsti[3][3],
658  const cs_real_t gradstj[3][3],
659  const cs_real_t i_massflux,
660  cs_real_t *testij,
661  cs_real_t *tesqck)
662 {
663  cs_real_t testi[3], testj[3];
664  cs_real_t dcc[3], ddi[3], ddj[3];
665  *testij = 0.;
666  *tesqck = 0.;
667 
668  /* Slope test */
669 
670  for (int i = 0; i < 3; i++) {
671  *testij += gradsti[i][0]*gradstj[i][0]
672  + gradsti[i][1]*gradstj[i][1]
673  + gradsti[i][2]*gradstj[i][2];
674 
675  testi[i] = gradsti[i][0]*i_face_normal[0]
676  + gradsti[i][1]*i_face_normal[1]
677  + gradsti[i][2]*i_face_normal[2];
678  testj[i] = gradstj[i][0]*i_face_normal[0]
679  + gradstj[i][1]*i_face_normal[1]
680  + gradstj[i][2]*i_face_normal[2];
681 
682  if (i_massflux > 0.) {
683  dcc[i] = gradi[i][0]*i_face_normal[0]
684  + gradi[i][1]*i_face_normal[1]
685  + gradi[i][2]*i_face_normal[2];
686  ddi[i] = testi[i];
687  ddj[i] = (pj[i]-pi[i])/distf *srfan;
688  } else {
689  dcc[i] = gradj[i][0]*i_face_normal[0]
690  + gradj[i][1]*i_face_normal[1]
691  + gradj[i][2]*i_face_normal[2];
692  ddi[i] = (pj[i]-pi[i])/distf *srfan;
693  ddj[i] = testj[i];
694  }
695  }
696 
697  *tesqck = cs_math_3_square_norm(dcc) - cs_math_3_square_distance(ddi, ddj);
698 }
699 
700 /*----------------------------------------------------------------------------*/
717 /*----------------------------------------------------------------------------*/
718 
719 inline static void
720 cs_slope_test_tensor(const cs_real_t pi[6],
721  const cs_real_t pj[6],
722  const cs_real_t distf,
723  const cs_real_t srfan,
724  const cs_real_t i_face_normal[3],
725  const cs_real_t gradi[6][3],
726  const cs_real_t gradj[6][3],
727  const cs_real_t gradsti[6][3],
728  const cs_real_t gradstj[6][3],
729  const cs_real_t i_massflux,
730  cs_real_t *testij,
731  cs_real_t *tesqck)
732 {
733  cs_real_t testi[6], testj[6];
734  cs_real_t dcc[6], ddi[6], ddj[6];
735 
736  *testij = 0.;
737  *tesqck = 0.;
738 
739  /* Slope test */
740 
741  for (int ij = 0; ij < 6; ij++) {
742  *testij += gradsti[ij][0]*gradstj[ij][0]
743  + gradsti[ij][1]*gradstj[ij][1]
744  + gradsti[ij][2]*gradstj[ij][2];
745  testi[ij] = gradsti[ij][0]*i_face_normal[0]
746  + gradsti[ij][1]*i_face_normal[1]
747  + gradsti[ij][2]*i_face_normal[2];
748  testj[ij] = gradstj[ij][0]*i_face_normal[0]
749  + gradstj[ij][1]*i_face_normal[1]
750  + gradstj[ij][2]*i_face_normal[2];
751 
752  if (i_massflux > 0.) {
753  dcc[ij] = gradi[ij][0]*i_face_normal[0]
754  + gradi[ij][1]*i_face_normal[1]
755  + gradi[ij][2]*i_face_normal[2];
756  ddi[ij] = testi[ij];
757  ddj[ij] = (pj[ij]-pi[ij])/distf *srfan;
758  }
759  else {
760  dcc[ij] = gradj[ij][0]*i_face_normal[0]
761  + gradj[ij][1]*i_face_normal[1]
762  + gradj[ij][2]*i_face_normal[2];
763  ddi[ij] = (pj[ij]-pi[ij])/distf *srfan;
764  ddj[ij] = testj[ij];
765  }
766 
767  *tesqck += cs_math_sq(dcc[ij]) - cs_math_sq(ddi[ij]-ddj[ij]);
768  }
769 }
770 
771 /*----------------------------------------------------------------------------*/
787 /*----------------------------------------------------------------------------*/
788 
789 inline static void
790 cs_i_compute_quantities(const cs_real_t bldfrp,
791  const cs_real_3_t diipf,
792  const cs_real_3_t djjpf,
793  const cs_real_3_t gradi,
794  const cs_real_3_t gradj,
795  const cs_real_t pi,
796  const cs_real_t pj,
797  cs_real_t *recoi,
798  cs_real_t *recoj,
799  cs_real_t *pip,
800  cs_real_t *pjp)
801 {
802  cs_real_t gradpf[3] = {0.5*(gradi[0] + gradj[0]),
803  0.5*(gradi[1] + gradj[1]),
804  0.5*(gradi[2] + gradj[2])};
805 
806  /* reconstruction only if IRCFLP = 1 */
807  *recoi = bldfrp*(cs_math_3_dot_product(gradpf, diipf));
808  *recoj = bldfrp*(cs_math_3_dot_product(gradpf, djjpf));
809  *pip = pi + *recoi;
810  *pjp = pj + *recoj;
811 }
812 
813 /*----------------------------------------------------------------------------*/
829 /*----------------------------------------------------------------------------*/
830 
831 inline static void
832 cs_i_compute_quantities_vector(const cs_real_t bldfrp,
833  const cs_real_3_t diipf,
834  const cs_real_3_t djjpf,
835  const cs_real_33_t gradi,
836  const cs_real_33_t gradj,
837  const cs_real_3_t pi,
838  const cs_real_3_t pj,
839  cs_real_t recoi[3],
840  cs_real_t recoj[3],
841  cs_real_t pip[3],
842  cs_real_t pjp[3])
843 {
844  cs_real_3_t dpvf;
845 
846  /* x-y-z components, p = u, v, w */
847 
848  for (int isou = 0; isou < 3; isou++) {
849 
850  for (int jsou = 0; jsou < 3; jsou++)
851  dpvf[jsou] = 0.5*( gradi[isou][jsou]
852  + gradj[isou][jsou]);
853 
854  /* reconstruction only if IRCFLP = 1 */
855 
856  recoi[isou] = bldfrp*(cs_math_3_dot_product(dpvf, diipf));
857  recoj[isou] = bldfrp*(cs_math_3_dot_product(dpvf, djjpf));
858 
859  pip[isou] = pi[isou] + recoi[isou];
860  pjp[isou] = pj[isou] + recoj[isou];
861 
862  }
863 }
864 
865 /*----------------------------------------------------------------------------*/
881 /*----------------------------------------------------------------------------*/
882 
883 inline static void
884 cs_i_compute_quantities_tensor(const cs_real_t bldfrp,
885  const cs_real_3_t diipf,
886  const cs_real_3_t djjpf,
887  const cs_real_63_t gradi,
888  const cs_real_63_t gradj,
889  const cs_real_6_t pi,
890  const cs_real_6_t pj,
891  cs_real_t recoi[6],
892  cs_real_t recoj[6],
893  cs_real_t pip[6],
894  cs_real_t pjp[6])
895 {
896  cs_real_3_t dpvf;
897 
898  /* x-y-z components, p = u, v, w */
899 
900  for (int isou = 0; isou < 6; isou++) {
901 
902  for (int jsou = 0; jsou < 3; jsou++)
903  dpvf[jsou] = 0.5*( gradi[isou][jsou]
904  + gradj[isou][jsou]);
905 
906  /* reconstruction only if IRCFLP = 1 */
907 
908  recoi[isou] = bldfrp*(cs_math_3_dot_product(dpvf, diipf));
909  recoj[isou] = bldfrp*(cs_math_3_dot_product(dpvf, djjpf));
910 
911  pip[isou] = pi[isou] + recoi[isou];
912  pjp[isou] = pj[isou] + recoj[isou];
913 
914  }
915 }
916 
917 /*----------------------------------------------------------------------------*/
933 /*----------------------------------------------------------------------------*/
934 
935 inline static void
936 cs_i_relax_c_val(const double relaxp,
937  const cs_real_t pia,
938  const cs_real_t pja,
939  const cs_real_t recoi,
940  const cs_real_t recoj,
941  const cs_real_t pi,
942  const cs_real_t pj,
943  cs_real_t *pir,
944  cs_real_t *pjr,
945  cs_real_t *pipr,
946  cs_real_t *pjpr)
947 {
948  *pir = pi/relaxp - (1.-relaxp)/relaxp * pia;
949  *pjr = pj/relaxp - (1.-relaxp)/relaxp * pja;
950 
951  *pipr = *pir + recoi;
952  *pjpr = *pjr + recoj;
953 }
954 
955 /*----------------------------------------------------------------------------*/
971 /*----------------------------------------------------------------------------*/
972 
973 inline static void
974 cs_i_relax_c_val_vector(const double relaxp,
975  const cs_real_3_t pia,
976  const cs_real_3_t pja,
977  const cs_real_3_t recoi,
978  const cs_real_3_t recoj,
979  const cs_real_3_t pi,
980  const cs_real_3_t pj,
981  cs_real_t pir[3],
982  cs_real_t pjr[3],
983  cs_real_t pipr[3],
984  cs_real_t pjpr[3])
985 {
986  for (int isou = 0; isou < 3; isou++) {
987  pir[isou] = pi[isou] /relaxp - (1.-relaxp)/relaxp * pia[isou];
988  pjr[isou] = pj[isou] /relaxp - (1.-relaxp)/relaxp * pja[isou];
989 
990  pipr[isou] = pir[isou] + recoi[isou];
991  pjpr[isou] = pjr[isou] + recoj[isou];
992  }
993 }
994 
995 /*----------------------------------------------------------------------------*/
1011 /*----------------------------------------------------------------------------*/
1012 
1013 inline static void
1014 cs_i_relax_c_val_tensor(const cs_real_t relaxp,
1015  const cs_real_t pia[6],
1016  const cs_real_t pja[6],
1017  const cs_real_t recoi[6],
1018  const cs_real_t recoj[6],
1019  const cs_real_t pi[6],
1020  const cs_real_t pj[6],
1021  cs_real_t pir[6],
1022  cs_real_t pjr[6],
1023  cs_real_t pipr[6],
1024  cs_real_t pjpr[6])
1025 {
1026  for (int isou = 0; isou < 6; isou++) {
1027  pir[isou] = pi[isou] /relaxp - (1.-relaxp)/relaxp * pia[isou];
1028  pjr[isou] = pj[isou] /relaxp - (1.-relaxp)/relaxp * pja[isou];
1029 
1030  pipr[isou] = pir[isou] + recoi[isou];
1031  pjpr[isou] = pjr[isou] + recoj[isou];
1032  }
1033 }
1034 
1035 /*----------------------------------------------------------------------------*/
1042 /*----------------------------------------------------------------------------*/
1043 
1044 inline static void
1045 cs_upwind_f_val(const cs_real_t p,
1046  cs_real_t *pf)
1047 {
1048  *pf = p;
1049 }
1050 
1051 /*----------------------------------------------------------------------------*/
1058 /*----------------------------------------------------------------------------*/
1059 
1060 inline static void
1061 cs_upwind_f_val_vector(const cs_real_3_t p,
1062  cs_real_t pf[3])
1063 {
1064  for (int isou = 0; isou < 3; isou++)
1065  pf[isou] = p[isou];
1066 }
1067 
1068 /*----------------------------------------------------------------------------*/
1075 /*----------------------------------------------------------------------------*/
1076 
1077 inline static void
1078 cs_upwind_f_val_tensor(const cs_real_6_t p,
1079  cs_real_t pf[6])
1080 {
1081  for (int isou = 0; isou < 6; isou++)
1082  pf[isou] = p[isou];
1083 }
1084 
1085 /*----------------------------------------------------------------------------*/
1094 /*----------------------------------------------------------------------------*/
1095 
1096 inline static void
1097 cs_centered_f_val(const double pnd,
1098  const cs_real_t pip,
1099  const cs_real_t pjp,
1100  cs_real_t *pf)
1101 {
1102  *pf = pnd*pip + (1.-pnd)*pjp;
1103 }
1104 
1105 /*----------------------------------------------------------------------------*/
1114 /*----------------------------------------------------------------------------*/
1115 
1116 inline static void
1117 cs_centered_f_val_vector(const double pnd,
1118  const cs_real_3_t pip,
1119  const cs_real_3_t pjp,
1120  cs_real_t pf[3])
1121 {
1122  for (int isou = 0; isou < 3; isou++)
1123  pf[isou] = pnd*pip[isou] + (1.-pnd)*pjp[isou];
1124 }
1125 
1126 /*----------------------------------------------------------------------------*/
1135 /*----------------------------------------------------------------------------*/
1136 
1137 inline static void
1138 cs_centered_f_val_tensor(const double pnd,
1139  const cs_real_6_t pip,
1140  const cs_real_6_t pjp,
1141  cs_real_t pf[6])
1142 {
1143  for (int isou = 0; isou < 6; isou++)
1144  pf[isou] = pnd*pip[isou] + (1.-pnd)*pjp[isou];
1145 }
1146 
1147 /*----------------------------------------------------------------------------*/
1157 /*----------------------------------------------------------------------------*/
1158 
1159 inline static void
1160 cs_solu_f_val(const cs_real_3_t cell_cen,
1161  const cs_real_3_t i_face_cog,
1162  const cs_real_3_t grad,
1163  const cs_real_t p,
1164  cs_real_t *pf)
1165 {
1166  cs_real_t df[3];
1167 
1168  df[0] = i_face_cog[0] - cell_cen[0];
1169  df[1] = i_face_cog[1] - cell_cen[1];
1170  df[2] = i_face_cog[2] - cell_cen[2];
1171 
1172  *pf = p + cs_math_3_dot_product(df, grad);
1173 }
1174 
1175 /*----------------------------------------------------------------------------*/
1185 /*----------------------------------------------------------------------------*/
1186 
1187 inline static void
1188 cs_solu_f_val_vector(const cs_real_3_t cell_cen,
1189  const cs_real_3_t i_face_cog,
1190  const cs_real_33_t grad,
1191  const cs_real_3_t p,
1192  cs_real_t pf[3])
1193 {
1194  cs_real_t df[3];
1195 
1196  for (int jsou = 0; jsou < 3; jsou++)
1197  df[jsou] = i_face_cog[jsou] - cell_cen[jsou];
1198 
1199  for (int isou = 0; isou < 3; isou++) {
1200  pf[isou] = p[isou] + df[0]*grad[isou][0]
1201  + df[1]*grad[isou][1]
1202  + df[2]*grad[isou][2];
1203 
1204  }
1205 }
1206 
1207 /*----------------------------------------------------------------------------*/
1217 /*----------------------------------------------------------------------------*/
1218 
1219 inline static void
1220 cs_solu_f_val_tensor(const cs_real_3_t cell_cen,
1221  const cs_real_3_t i_face_cog,
1222  const cs_real_63_t grad,
1223  const cs_real_6_t p,
1224  cs_real_t pf[6])
1225 {
1226  cs_real_t df[3];
1227 
1228  for (int jsou = 0; jsou < 3; jsou++)
1229  df[jsou] = i_face_cog[jsou] - cell_cen[jsou];
1230 
1231  for (int isou = 0; isou < 6; isou++) {
1232  pf[isou] = p[isou] + df[0]*grad[isou][0]
1233  + df[1]*grad[isou][1]
1234  + df[2]*grad[isou][2];
1235  }
1236 }
1237 
1238 /*----------------------------------------------------------------------------*/
1248 /*----------------------------------------------------------------------------*/
1249 
1250 inline static void
1251 cs_blend_f_val(const double blencp,
1252  const cs_real_t p,
1253  cs_real_t *pf)
1254 {
1255  *pf = blencp * (*pf) + (1. - blencp) * p;
1256 }
1257 
1258 /*----------------------------------------------------------------------------*/
1268 /*----------------------------------------------------------------------------*/
1269 
1270 inline static void
1271 cs_blend_f_val_vector(const double blencp,
1272  const cs_real_3_t p,
1273  cs_real_t pf[3])
1274 {
1275  for (int isou = 0; isou < 3; isou++)
1276  pf[isou] = blencp*(pf[isou])+(1.-blencp)*p[isou];
1277 }
1278 
1279 /*----------------------------------------------------------------------------*/
1289 /*----------------------------------------------------------------------------*/
1290 
1291 inline static void
1292 cs_blend_f_val_tensor(const double blencp,
1293  const cs_real_t p[6],
1294  cs_real_t pf[6])
1295 {
1296  for (int isou = 0; isou < 6; isou++)
1297  pf[isou] = blencp*(pf[isou])+(1.-blencp)*p[isou];
1298 }
1299 
1300 /*----------------------------------------------------------------------------*/
1321 /*----------------------------------------------------------------------------*/
1322 
1323 inline static void
1324 cs_i_conv_flux(const int iconvp,
1325  const cs_real_t thetap,
1326  const int imasac,
1327  const cs_real_t pi,
1328  const cs_real_t pj,
1329  const cs_real_t pifri,
1330  const cs_real_t pifrj,
1331  const cs_real_t pjfri,
1332  const cs_real_t pjfrj,
1333  const cs_real_t i_massflux,
1334  const cs_real_t xcppi,
1335  const cs_real_t xcppj,
1336  cs_real_2_t fluxij)
1337 {
1338  cs_real_t flui, fluj;
1339 
1340  flui = 0.5*(i_massflux + fabs(i_massflux));
1341  fluj = 0.5*(i_massflux - fabs(i_massflux));
1342 
1343  fluxij[0] += iconvp*xcppi*(thetap*(flui*pifri + fluj*pjfri) - imasac*i_massflux*pi);
1344  fluxij[1] += iconvp*xcppj*(thetap*(flui*pifrj + fluj*pjfrj) - imasac*i_massflux*pj);
1345 }
1346 
1347 /*----------------------------------------------------------------------------*/
1365 /*----------------------------------------------------------------------------*/
1366 
1367 inline static void
1368 cs_i_conv_flux_vector(const int iconvp,
1369  const cs_real_t thetap,
1370  const int imasac,
1371  const cs_real_t pi[3],
1372  const cs_real_t pj[3],
1373  const cs_real_t pifri[3],
1374  const cs_real_t pifrj[3],
1375  const cs_real_t pjfri[3],
1376  const cs_real_t pjfrj[3],
1377  const cs_real_t i_massflux,
1378  cs_real_t fluxi[3],
1379  cs_real_t fluxj[3])
1380 {
1381  cs_real_t flui, fluj;
1382 
1383  flui = 0.5*(i_massflux + fabs(i_massflux));
1384  fluj = 0.5*(i_massflux - fabs(i_massflux));
1385 
1386  for (int isou = 0; isou < 3; isou++) {
1387 
1388  fluxi[isou] += iconvp*( thetap*(flui*pifri[isou] + fluj*pjfri[isou])
1389  - imasac*i_massflux*pi[isou]);
1390  fluxj[isou] += iconvp*( thetap*(flui*pifrj[isou] + fluj*pjfrj[isou])
1391  - imasac*i_massflux*pj[isou]);
1392  }
1393 }
1394 
1395 /*----------------------------------------------------------------------------*/
1413 /*----------------------------------------------------------------------------*/
1414 
1415 inline static void
1416 cs_i_conv_flux_tensor(const int iconvp,
1417  const cs_real_t thetap,
1418  const int imasac,
1419  const cs_real_t pi[6],
1420  const cs_real_t pj[6],
1421  const cs_real_t pifri[6],
1422  const cs_real_t pifrj[6],
1423  const cs_real_t pjfri[6],
1424  const cs_real_t pjfrj[6],
1425  const cs_real_t i_massflux,
1426  cs_real_t fluxi[6],
1427  cs_real_t fluxj[6])
1428 {
1429  cs_real_t flui, fluj;
1430 
1431  flui = 0.5*(i_massflux + fabs(i_massflux));
1432  fluj = 0.5*(i_massflux - fabs(i_massflux));
1433 
1434  for (int isou = 0; isou < 6; isou++) {
1435  fluxi[isou] += iconvp*( thetap*(flui*pifri[isou] + fluj*pjfri[isou])
1436  - imasac*i_massflux*pi[isou]);
1437  fluxj[isou] += iconvp*( thetap*(flui*pifrj[isou] + fluj*pjfrj[isou])
1438  - imasac*i_massflux*pj[isou]);
1439  }
1440 }
1441 
1442 /*----------------------------------------------------------------------------*/
1455 /*----------------------------------------------------------------------------*/
1456 
1457 inline static void
1458 cs_i_diff_flux(const int idiffp,
1459  const cs_real_t thetap,
1460  const cs_real_t pip,
1461  const cs_real_t pjp,
1462  const cs_real_t pipr,
1463  const cs_real_t pjpr,
1464  const cs_real_t i_visc,
1465  cs_real_t fluxij[2])
1466 {
1467  fluxij[0] += idiffp*thetap*i_visc*(pipr -pjp);
1468  fluxij[1] += idiffp*thetap*i_visc*(pip -pjpr);
1469 }
1470 
1471 /*----------------------------------------------------------------------------*/
1485 /*----------------------------------------------------------------------------*/
1486 
1487 inline static void
1488 cs_i_diff_flux_vector(const int idiffp,
1489  const cs_real_t thetap,
1490  const cs_real_t pip[3],
1491  const cs_real_t pjp[3],
1492  const cs_real_t pipr[3],
1493  const cs_real_t pjpr[3],
1494  const cs_real_t i_visc,
1495  cs_real_t fluxi[3],
1496  cs_real_t fluxj[3])
1497 {
1498  for (int isou = 0; isou < 3; isou++) {
1499  fluxi[isou] += idiffp*thetap*i_visc*(pipr[isou] -pjp[isou]);
1500  fluxj[isou] += idiffp*thetap*i_visc*(pip[isou] -pjpr[isou]);
1501  }
1502 }
1503 
1504 /*----------------------------------------------------------------------------*/
1518 /*----------------------------------------------------------------------------*/
1519 
1520 inline static void
1521 cs_i_diff_flux_tensor(const int idiffp,
1522  const cs_real_t thetap,
1523  const cs_real_t pip[6],
1524  const cs_real_t pjp[6],
1525  const cs_real_t pipr[6],
1526  const cs_real_t pjpr[6],
1527  const cs_real_t i_visc,
1528  cs_real_t fluxi[6],
1529  cs_real_t fluxj[6])
1530 {
1531  for (int isou = 0; isou < 6; isou++) {
1532  fluxi[isou] += idiffp*thetap*i_visc*(pipr[isou] -pjp[isou]);
1533  fluxj[isou] += idiffp*thetap*i_visc*(pip[isou] -pjpr[isou]);
1534  }
1535 }
1536 
1537 /*----------------------------------------------------------------------------*/
1561 /*----------------------------------------------------------------------------*/
1562 
1563 inline static void
1564 cs_i_cd_steady_upwind(const cs_real_t bldfrp,
1565  const cs_real_t relaxp,
1566  const cs_real_t diipf[3],
1567  const cs_real_t djjpf[3],
1568  const cs_real_t gradi[3],
1569  const cs_real_t gradj[3],
1570  const cs_real_t pi,
1571  const cs_real_t pj,
1572  const cs_real_t pia,
1573  const cs_real_t pja,
1574  cs_real_t *pifri,
1575  cs_real_t *pifrj,
1576  cs_real_t *pjfri,
1577  cs_real_t *pjfrj,
1578  cs_real_t *pip,
1579  cs_real_t *pjp,
1580  cs_real_t *pipr,
1581  cs_real_t *pjpr)
1582 {
1583  cs_real_t pir, pjr;
1584  cs_real_t recoi, recoj;
1585 
1586  cs_i_compute_quantities(bldfrp,
1587  diipf,
1588  djjpf,
1589  gradi,
1590  gradj,
1591  pi,
1592  pj,
1593  &recoi,
1594  &recoj,
1595  pip,
1596  pjp);
1597 
1598  cs_i_relax_c_val(relaxp,
1599  pia,
1600  pja,
1601  recoi,
1602  recoj,
1603  pi,
1604  pj,
1605  &pir,
1606  &pjr,
1607  pipr,
1608  pjpr);
1609 
1610  cs_upwind_f_val(pi,
1611  pifrj);
1612  cs_upwind_f_val(pir,
1613  pifri);
1614  cs_upwind_f_val(pj,
1615  pjfri);
1616  cs_upwind_f_val(pjr,
1617  pjfrj);
1618 }
1619 
1620 /*----------------------------------------------------------------------------*/
1644 /*----------------------------------------------------------------------------*/
1645 
1646 inline static void
1647 cs_i_cd_steady_upwind_vector(const cs_real_t bldfrp,
1648  const cs_real_t relaxp,
1649  const cs_real_t diipf[3],
1650  const cs_real_t djjpf[3],
1651  const cs_real_t gradi[3][3],
1652  const cs_real_t gradj[3][3],
1653  const cs_real_t pi[3],
1654  const cs_real_t pj[3],
1655  const cs_real_t pia[3],
1656  const cs_real_t pja[3],
1657  cs_real_t pifri[3],
1658  cs_real_t pifrj[3],
1659  cs_real_t pjfri[3],
1660  cs_real_t pjfrj[3],
1661  cs_real_t pip[3],
1662  cs_real_t pjp[3],
1663  cs_real_t pipr[3],
1664  cs_real_t pjpr[3])
1665 {
1666  cs_real_t pir[3], pjr[3];
1667  cs_real_t recoi[3], recoj[3];
1668 
1669  cs_i_compute_quantities_vector(bldfrp,
1670  diipf,
1671  djjpf,
1672  gradi,
1673  gradj,
1674  pi,
1675  pj,
1676  recoi,
1677  recoj,
1678  pip,
1679  pjp);
1680 
1681  cs_i_relax_c_val_vector(relaxp,
1682  pia,
1683  pja,
1684  recoi,
1685  recoj,
1686  pi,
1687  pj,
1688  pir,
1689  pjr,
1690  pipr,
1691  pjpr);
1692 
1693  cs_upwind_f_val_vector(pi,
1694  pifrj);
1695  cs_upwind_f_val_vector(pir,
1696  pifri);
1697  cs_upwind_f_val_vector(pj,
1698  pjfri);
1699  cs_upwind_f_val_vector(pjr,
1700  pjfrj);
1701 }
1702 
1703 /*----------------------------------------------------------------------------*/
1727 /*----------------------------------------------------------------------------*/
1728 
1729 inline static void
1730 cs_i_cd_steady_upwind_tensor(const cs_real_t bldfrp,
1731  const cs_real_t relaxp,
1732  const cs_real_t diipf[3],
1733  const cs_real_t djjpf[3],
1734  const cs_real_t gradi[6][3],
1735  const cs_real_t gradj[6][3],
1736  const cs_real_t pi[6],
1737  const cs_real_t pj[6],
1738  const cs_real_t pia[6],
1739  const cs_real_t pja[6],
1740  cs_real_t pifri[6],
1741  cs_real_t pifrj[6],
1742  cs_real_t pjfri[6],
1743  cs_real_t pjfrj[6],
1744  cs_real_t pip[6],
1745  cs_real_t pjp[6],
1746  cs_real_t pipr[6],
1747  cs_real_t pjpr[6])
1748 {
1749  cs_real_t pir[6], pjr[6];
1750  cs_real_t recoi[6], recoj[6];
1751 
1752  cs_i_compute_quantities_tensor(bldfrp,
1753  diipf,
1754  djjpf,
1755  gradi,
1756  gradj,
1757  pi,
1758  pj,
1759  recoi,
1760  recoj,
1761  pip,
1762  pjp);
1763 
1764  cs_i_relax_c_val_tensor(relaxp,
1765  pia,
1766  pja,
1767  recoi,
1768  recoj,
1769  pi,
1770  pj,
1771  pir,
1772  pjr,
1773  pipr,
1774  pjpr);
1775 
1776  cs_upwind_f_val_tensor(pi,
1777  pifrj);
1778  cs_upwind_f_val_tensor(pir,
1779  pifri);
1780  cs_upwind_f_val_tensor(pj,
1781  pjfri);
1782  cs_upwind_f_val_tensor(pjr,
1783  pjfrj);
1784 }
1785 
1786 /*----------------------------------------------------------------------------*/
1803 /*----------------------------------------------------------------------------*/
1804 
1805 inline static void
1806 cs_i_cd_unsteady_upwind(const cs_real_t bldfrp,
1807  const cs_real_t diipf[3],
1808  const cs_real_t djjpf[3],
1809  const cs_real_t gradi[3],
1810  const cs_real_t gradj[3],
1811  const cs_real_t pi,
1812  const cs_real_t pj,
1813  cs_real_t *pif,
1814  cs_real_t *pjf,
1815  cs_real_t *pip,
1816  cs_real_t *pjp)
1817 {
1818  cs_real_t recoi, recoj;
1819 
1820  cs_i_compute_quantities(bldfrp,
1821  diipf,
1822  djjpf,
1823  gradi,
1824  gradj,
1825  pi,
1826  pj,
1827  &recoi,
1828  &recoj,
1829  pip,
1830  pjp);
1831 
1832  cs_upwind_f_val(pi, pif);
1833  cs_upwind_f_val(pj, pjf);
1834 }
1835 
1836 /*----------------------------------------------------------------------------*/
1853 /*----------------------------------------------------------------------------*/
1854 
1855 inline static void
1856 cs_i_cd_unsteady_upwind_vector(const cs_real_t bldfrp,
1857  const cs_real_t diipf[3],
1858  const cs_real_t djjpf[3],
1859  const cs_real_t gradi[3][3],
1860  const cs_real_t gradj[3][3],
1861  const cs_real_t pi[3],
1862  const cs_real_t pj[3],
1863  cs_real_t pif[3],
1864  cs_real_t pjf[3],
1865  cs_real_t pip[3],
1866  cs_real_t pjp[3])
1867 {
1868  cs_real_t recoi[3], recoj[3];
1869 
1870  cs_i_compute_quantities_vector(bldfrp,
1871  diipf,
1872  djjpf,
1873  gradi,
1874  gradj,
1875  pi,
1876  pj,
1877  recoi,
1878  recoj,
1879  pip,
1880  pjp);
1881 
1882  cs_upwind_f_val_vector(pi, pif);
1883  cs_upwind_f_val_vector(pj, pjf);
1884 }
1885 
1886 /*----------------------------------------------------------------------------*/
1903 /*----------------------------------------------------------------------------*/
1904 
1905 inline static void
1906 cs_i_cd_unsteady_upwind_tensor(const cs_real_t bldfrp,
1907  const cs_real_t diipf[3],
1908  const cs_real_t djjpf[3],
1909  const cs_real_t gradi[6][3],
1910  const cs_real_t gradj[6][3],
1911  const cs_real_t pi[6],
1912  const cs_real_t pj[6],
1913  cs_real_t pif[6],
1914  cs_real_t pjf[6],
1915  cs_real_t pip[6],
1916  cs_real_t pjp[6])
1917 {
1918  cs_real_t recoi[6], recoj[6];
1919 
1920  cs_i_compute_quantities_tensor(bldfrp,
1921  diipf,
1922  djjpf,
1923  gradi,
1924  gradj,
1925  pi,
1926  pj,
1927  recoi,
1928  recoj,
1929  pip,
1930  pjp);
1931 
1932  cs_upwind_f_val_tensor(pi, pif);
1933  cs_upwind_f_val_tensor(pj, pjf);
1934 
1935 }
1936 
1937 /*----------------------------------------------------------------------------*/
1970 /*----------------------------------------------------------------------------*/
1971 
1972 inline static void
1973 cs_i_cd_steady(const cs_real_t bldfrp,
1974  const int ischcp,
1975  const double relaxp,
1976  const double blencp,
1977  const cs_real_t weight,
1978  const cs_real_t cell_ceni[3],
1979  const cs_real_t cell_cenj[3],
1980  const cs_real_t i_face_cog[3],
1981  const cs_real_t diipf[3],
1982  const cs_real_t djjpf[3],
1983  const cs_real_t gradi[3],
1984  const cs_real_t gradj[3],
1985  const cs_real_t gradupi[3],
1986  const cs_real_t gradupj[3],
1987  const cs_real_t pi,
1988  const cs_real_t pj,
1989  const cs_real_t pia,
1990  const cs_real_t pja,
1991  cs_real_t *pifri,
1992  cs_real_t *pifrj,
1993  cs_real_t *pjfri,
1994  cs_real_t *pjfrj,
1995  cs_real_t *pip,
1996  cs_real_t *pjp,
1997  cs_real_t *pipr,
1998  cs_real_t *pjpr)
1999 {
2000  cs_real_t pir, pjr;
2001  cs_real_t recoi, recoj;
2002 
2003  cs_i_compute_quantities(bldfrp,
2004  diipf,
2005  djjpf,
2006  gradi,
2007  gradj,
2008  pi,
2009  pj,
2010  &recoi,
2011  &recoj,
2012  pip,
2013  pjp);
2014 
2015  cs_i_relax_c_val(relaxp,
2016  pia,
2017  pja,
2018  recoi,
2019  recoj,
2020  pi,
2021  pj,
2022  &pir,
2023  &pjr,
2024  pipr,
2025  pjpr);
2026 
2027  if (ischcp == 1) {
2028 
2029  /* Centered
2030  --------*/
2031 
2032  cs_centered_f_val(weight,
2033  *pip,
2034  *pjpr,
2035  pifrj);
2036  cs_centered_f_val(weight,
2037  *pipr,
2038  *pjp,
2039  pifri);
2040  cs_centered_f_val(weight,
2041  *pipr,
2042  *pjp,
2043  pjfri);
2044  cs_centered_f_val(weight,
2045  *pip,
2046  *pjpr,
2047  pjfrj);
2048 
2049  } else if (ischcp == 0) {
2050 
2051  /* Original SOLU
2052  --------------*/
2053 
2054  cs_solu_f_val(cell_ceni,
2055  i_face_cog,
2056  gradi,
2057  pi,
2058  pifrj);
2059  cs_solu_f_val(cell_ceni,
2060  i_face_cog,
2061  gradi,
2062  pir,
2063  pifri);
2064  cs_solu_f_val(cell_cenj,
2065  i_face_cog,
2066  gradj,
2067  pj,
2068  pjfri);
2069  cs_solu_f_val(cell_cenj,
2070  i_face_cog,
2071  gradj,
2072  pjr,
2073  pjfrj);
2074 
2075  } else {
2076 
2077  /* SOLU
2078  ----*/
2079 
2080  cs_solu_f_val(cell_ceni,
2081  i_face_cog,
2082  gradupi,
2083  pi,
2084  pifrj);
2085  cs_solu_f_val(cell_ceni,
2086  i_face_cog,
2087  gradupi,
2088  pir,
2089  pifri);
2090  cs_solu_f_val(cell_cenj,
2091  i_face_cog,
2092  gradupj,
2093  pj,
2094  pjfri);
2095  cs_solu_f_val(cell_cenj,
2096  i_face_cog,
2097  gradupj,
2098  pjr,
2099  pjfrj);
2100 
2101  }
2102 
2103  /* Blending
2104  --------*/
2105 
2106  cs_blend_f_val(blencp,
2107  pi,
2108  pifrj);
2109  cs_blend_f_val(blencp,
2110  pir,
2111  pifri);
2112  cs_blend_f_val(blencp,
2113  pj,
2114  pjfri);
2115  cs_blend_f_val(blencp,
2116  pjr,
2117  pjfrj);
2118 }
2119 
2120 /*----------------------------------------------------------------------------*/
2151 /*----------------------------------------------------------------------------*/
2152 
2153 inline static void
2154 cs_i_cd_steady_vector(const cs_real_t bldfrp,
2155  const int ischcp,
2156  const double relaxp,
2157  const double blencp,
2158  const cs_real_t weight,
2159  const cs_real_3_t cell_ceni,
2160  const cs_real_3_t cell_cenj,
2161  const cs_real_3_t i_face_cog,
2162  const cs_real_3_t diipf,
2163  const cs_real_3_t djjpf,
2164  const cs_real_33_t gradi,
2165  const cs_real_33_t gradj,
2166  const cs_real_3_t pi,
2167  const cs_real_3_t pj,
2168  const cs_real_3_t pia,
2169  const cs_real_3_t pja,
2170  cs_real_t pifri[3],
2171  cs_real_t pifrj[3],
2172  cs_real_t pjfri[3],
2173  cs_real_t pjfrj[3],
2174  cs_real_t pip[3],
2175  cs_real_t pjp[3],
2176  cs_real_t pipr[3],
2177  cs_real_t pjpr[3])
2178 {
2179  cs_real_t pir[3], pjr[3];
2180  cs_real_t recoi[3], recoj[3];
2181 
2182  cs_i_compute_quantities_vector(bldfrp,
2183  diipf,
2184  djjpf,
2185  gradi,
2186  gradj,
2187  pi,
2188  pj,
2189  recoi,
2190  recoj,
2191  pip,
2192  pjp);
2193 
2194  cs_i_relax_c_val_vector(relaxp,
2195  pia,
2196  pja,
2197  recoi,
2198  recoj,
2199  pi,
2200  pj,
2201  pir,
2202  pjr,
2203  pipr,
2204  pjpr);
2205 
2206  if (ischcp == 1) {
2207 
2208  /* Centered
2209  --------*/
2210 
2211  cs_centered_f_val_vector(weight,
2212  pip,
2213  pjpr,
2214  pifrj);
2215  cs_centered_f_val_vector(weight,
2216  pipr,
2217  pjp,
2218  pifri);
2219  cs_centered_f_val_vector(weight,
2220  pipr,
2221  pjp,
2222  pjfri);
2223  cs_centered_f_val_vector(weight,
2224  pip,
2225  pjpr,
2226  pjfrj);
2227 
2228  } else {
2229 
2230  /* Second order
2231  ------------*/
2232 
2233  cs_solu_f_val_vector(cell_ceni,
2234  i_face_cog,
2235  gradi,
2236  pi,
2237  pifrj);
2238  cs_solu_f_val_vector(cell_ceni,
2239  i_face_cog,
2240  gradi,
2241  pir,
2242  pifri);
2243  cs_solu_f_val_vector(cell_cenj,
2244  i_face_cog,
2245  gradj,
2246  pj,
2247  pjfri);
2248  cs_solu_f_val_vector(cell_cenj,
2249  i_face_cog,
2250  gradj,
2251  pjr,
2252  pjfrj);
2253 
2254  }
2255 
2256  /* Blending
2257  --------*/
2258  cs_blend_f_val_vector(blencp,
2259  pi,
2260  pifrj);
2261  cs_blend_f_val_vector(blencp,
2262  pir,
2263  pifri);
2264  cs_blend_f_val_vector(blencp,
2265  pj,
2266  pjfri);
2267  cs_blend_f_val_vector(blencp,
2268  pjr,
2269  pjfrj);
2270 
2271 }
2272 
2273 /*----------------------------------------------------------------------------*/
2304 /*----------------------------------------------------------------------------*/
2305 
2306 inline static void
2307 cs_i_cd_steady_tensor(const cs_real_t bldfrp,
2308  const int ischcp,
2309  const double relaxp,
2310  const double blencp,
2311  const cs_real_t weight,
2312  const cs_real_3_t cell_ceni,
2313  const cs_real_3_t cell_cenj,
2314  const cs_real_3_t i_face_cog,
2315  const cs_real_3_t diipf,
2316  const cs_real_3_t djjpf,
2317  const cs_real_63_t gradi,
2318  const cs_real_63_t gradj,
2319  const cs_real_6_t pi,
2320  const cs_real_6_t pj,
2321  const cs_real_6_t pia,
2322  const cs_real_6_t pja,
2323  cs_real_t pifri[6],
2324  cs_real_t pifrj[6],
2325  cs_real_t pjfri[6],
2326  cs_real_t pjfrj[6],
2327  cs_real_t pip[6],
2328  cs_real_t pjp[6],
2329  cs_real_t pipr[6],
2330  cs_real_t pjpr[6])
2331 
2332 {
2333  cs_real_t pir[6], pjr[6];
2334  cs_real_t recoi[6], recoj[6];
2335 
2336  cs_i_compute_quantities_tensor(bldfrp,
2337  diipf,
2338  djjpf,
2339  gradi,
2340  gradj,
2341  pi,
2342  pj,
2343  recoi,
2344  recoj,
2345  pip,
2346  pjp);
2347 
2348  cs_i_relax_c_val_tensor(relaxp,
2349  pia,
2350  pja,
2351  recoi,
2352  recoj,
2353  pi,
2354  pj,
2355  pir,
2356  pjr,
2357  pipr,
2358  pjpr);
2359 
2360  if (ischcp == 1) {
2361 
2362  /* Centered
2363  --------*/
2364 
2365  cs_centered_f_val_tensor(weight,
2366  pip,
2367  pjpr,
2368  pifrj);
2369  cs_centered_f_val_tensor(weight,
2370  pipr,
2371  pjp,
2372  pifri);
2373  cs_centered_f_val_tensor(weight,
2374  pipr,
2375  pjp,
2376  pjfri);
2377  cs_centered_f_val_tensor(weight,
2378  pip,
2379  pjpr,
2380  pjfrj);
2381 
2382  } else {
2383 
2384  /* Second order
2385  ------------*/
2386 
2387  cs_solu_f_val_tensor(cell_ceni,
2388  i_face_cog,
2389  gradi,
2390  pi,
2391  pifrj);
2392  cs_solu_f_val_tensor(cell_ceni,
2393  i_face_cog,
2394  gradi,
2395  pir,
2396  pifri);
2397  cs_solu_f_val_tensor(cell_cenj,
2398  i_face_cog,
2399  gradj,
2400  pj,
2401  pjfri);
2402  cs_solu_f_val_tensor(cell_cenj,
2403  i_face_cog,
2404  gradj,
2405  pjr,
2406  pjfrj);
2407 
2408  }
2409 
2410  /* Blending
2411  --------*/
2412 
2413  cs_blend_f_val_tensor(blencp,
2414  pi,
2415  pifrj);
2416  cs_blend_f_val_tensor(blencp,
2417  pir,
2418  pifri);
2419  cs_blend_f_val_tensor(blencp,
2420  pj,
2421  pjfri);
2422  cs_blend_f_val_tensor(blencp,
2423  pjr,
2424  pjfrj);
2425 }
2426 
2427 /*----------------------------------------------------------------------------*/
2457 /*----------------------------------------------------------------------------*/
2458 
2459 inline static void
2460 cs_i_cd_unsteady(const cs_real_t bldfrp,
2461  const int ischcp,
2462  const double blencp,
2463  const cs_real_t weight,
2464  const cs_real_3_t cell_ceni,
2465  const cs_real_3_t cell_cenj,
2466  const cs_real_3_t i_face_cog,
2467  const cs_real_t hybrid_blend_i,
2468  const cs_real_t hybrid_blend_j,
2469  const cs_real_3_t diipf,
2470  const cs_real_3_t djjpf,
2471  const cs_real_3_t gradi,
2472  const cs_real_3_t gradj,
2473  const cs_real_3_t gradupi,
2474  const cs_real_3_t gradupj,
2475  const cs_real_t pi,
2476  const cs_real_t pj,
2477  cs_real_t *pif,
2478  cs_real_t *pjf,
2479  cs_real_t *pip,
2480  cs_real_t *pjp)
2481 {
2482  cs_real_t recoi, recoj;
2483 
2484  cs_i_compute_quantities(bldfrp,
2485  diipf,
2486  djjpf,
2487  gradi,
2488  gradj,
2489  pi,
2490  pj,
2491  &recoi,
2492  &recoj,
2493  pip,
2494  pjp);
2495 
2496 
2497  if (ischcp == 1) {
2498 
2499  /* Centered
2500  --------*/
2501 
2502  cs_centered_f_val(weight,
2503  *pip,
2504  *pjp,
2505  pif);
2506  cs_centered_f_val(weight,
2507  *pip,
2508  *pjp,
2509  pjf);
2510 
2511  } else if (ischcp == 0) {
2512 
2513  /* Legacy SOLU
2514  -----------*/
2515 
2516  cs_solu_f_val(cell_ceni,
2517  i_face_cog,
2518  gradi,
2519  pi,
2520  pif);
2521  cs_solu_f_val(cell_cenj,
2522  i_face_cog,
2523  gradj,
2524  pj,
2525  pjf);
2526 
2527  } else if (ischcp == 3) {
2528 
2529  /* Centered
2530  --------*/
2531 
2532  cs_centered_f_val(weight,
2533  *pip,
2534  *pjp,
2535  pif);
2536  cs_centered_f_val(weight,
2537  *pip,
2538  *pjp,
2539  pjf);
2540 
2541  /* Legacy SOLU
2542  -----------*/
2543  cs_real_t pif_up, pjf_up;
2544  cs_real_t hybrid_blend_interp;
2545 
2546  cs_solu_f_val(cell_ceni,
2547  i_face_cog,
2548  gradi,
2549  pi,
2550  &pif_up);
2551  cs_solu_f_val(cell_cenj,
2552  i_face_cog,
2553  gradj,
2554  pj,
2555  &pjf_up);
2556 
2557  hybrid_blend_interp = fmin(hybrid_blend_i,hybrid_blend_j);
2558  *pif = hybrid_blend_interp*(*pif) + (1. - hybrid_blend_interp)*pif_up;
2559  *pjf = hybrid_blend_interp*(*pjf) + (1. - hybrid_blend_interp)*pjf_up;
2560 
2561  } else {
2562 
2563  /* SOLU
2564  ----*/
2565 
2566  cs_solu_f_val(cell_ceni,
2567  i_face_cog,
2568  gradupi,
2569  pi,
2570  pif);
2571  cs_solu_f_val(cell_cenj,
2572  i_face_cog,
2573  gradupj,
2574  pj,
2575  pjf);
2576 
2577  }
2578 
2579 
2580  /* Blending
2581  --------*/
2582 
2583  cs_blend_f_val(blencp,
2584  pi,
2585  pif);
2586  cs_blend_f_val(blencp,
2587  pj,
2588  pjf);
2589 }
2590 
2591 /*----------------------------------------------------------------------------*/
2617 /*----------------------------------------------------------------------------*/
2618 
2619 inline static void
2620 cs_i_cd_unsteady_vector(const cs_real_t bldfrp,
2621  const int ischcp,
2622  const double blencp,
2623  const cs_real_t weight,
2624  const cs_real_3_t cell_ceni,
2625  const cs_real_3_t cell_cenj,
2626  const cs_real_3_t i_face_cog,
2627  const cs_real_t hybrid_blend_i,
2628  const cs_real_t hybrid_blend_j,
2629  const cs_real_3_t diipf,
2630  const cs_real_3_t djjpf,
2631  const cs_real_33_t gradi,
2632  const cs_real_33_t gradj,
2633  const cs_real_3_t pi,
2634  const cs_real_3_t pj,
2635  cs_real_t pif[3],
2636  cs_real_t pjf[3],
2637  cs_real_t pip[3],
2638  cs_real_t pjp[3])
2639 
2640 {
2641  cs_real_t recoi[3], recoj[3];
2642 
2643  cs_i_compute_quantities_vector(bldfrp,
2644  diipf,
2645  djjpf,
2646  gradi,
2647  gradj,
2648  pi,
2649  pj,
2650  recoi,
2651  recoj,
2652  pip,
2653  pjp);
2654 
2655  if (ischcp == 1) {
2656 
2657  /* Centered
2658  --------*/
2659 
2660  cs_centered_f_val_vector(weight,
2661  pip,
2662  pjp,
2663  pif);
2664  cs_centered_f_val_vector(weight,
2665  pip,
2666  pjp,
2667  pjf);
2668  } else if (ischcp == 3) {
2669 
2670  /* Centered
2671  --------*/
2672 
2673  cs_centered_f_val_vector(weight,
2674  pip,
2675  pjp,
2676  pif);
2677  cs_centered_f_val_vector(weight,
2678  pip,
2679  pjp,
2680  pjf);
2681 
2682  /* SOLU
2683  -----*/
2684  cs_real_t pif_up[3], pjf_up[3];
2685  cs_real_t hybrid_blend_interp;
2686 
2687  cs_solu_f_val_vector(cell_ceni,
2688  i_face_cog,
2689  gradi,
2690  pi,
2691  pif_up);
2692  cs_solu_f_val_vector(cell_cenj,
2693  i_face_cog,
2694  gradj,
2695  pj,
2696  pjf_up);
2697 
2698  hybrid_blend_interp = fmin(hybrid_blend_i,hybrid_blend_j);
2699  for (int isou = 0; isou < 3; isou++) {
2700  pif[isou] = hybrid_blend_interp *pif[isou]
2701  + (1. - hybrid_blend_interp)*pif_up[isou];
2702  pjf[isou] = hybrid_blend_interp *pjf[isou]
2703  + (1. - hybrid_blend_interp)*pjf_up[isou];
2704  }
2705  } else {
2706 
2707  /* Second order
2708  ------------*/
2709 
2710  cs_solu_f_val_vector(cell_ceni,
2711  i_face_cog,
2712  gradi,
2713  pi,
2714  pif);
2715  cs_solu_f_val_vector(cell_cenj,
2716  i_face_cog,
2717  gradj,
2718  pj,
2719  pjf);
2720 
2721  }
2722 
2723  /* Blending
2724  --------*/
2725 
2726  cs_blend_f_val_vector(blencp,
2727  pi,
2728  pif);
2729  cs_blend_f_val_vector(blencp,
2730  pj,
2731  pjf);
2732 
2733 }
2734 
2735 /*----------------------------------------------------------------------------*/
2759 /*----------------------------------------------------------------------------*/
2760 
2761 inline static void
2762 cs_i_cd_unsteady_tensor(const cs_real_t bldfrp,
2763  const int ischcp,
2764  const double blencp,
2765  const cs_real_t weight,
2766  const cs_real_3_t cell_ceni,
2767  const cs_real_3_t cell_cenj,
2768  const cs_real_3_t i_face_cog,
2769  const cs_real_3_t diipf,
2770  const cs_real_3_t djjpf,
2771  const cs_real_63_t gradi,
2772  const cs_real_63_t gradj,
2773  const cs_real_6_t pi,
2774  const cs_real_6_t pj,
2775  cs_real_t pif[6],
2776  cs_real_t pjf[6],
2777  cs_real_t pip[6],
2778  cs_real_t pjp[6])
2779 
2780 {
2781  cs_real_t recoi[6], recoj[6];
2782 
2783  cs_i_compute_quantities_tensor(bldfrp,
2784  diipf,
2785  djjpf,
2786  gradi,
2787  gradj,
2788  pi,
2789  pj,
2790  recoi,
2791  recoj,
2792  pip,
2793  pjp);
2794 
2795  if (ischcp == 1) {
2796 
2797  /* Centered
2798  --------*/
2799 
2800  cs_centered_f_val_tensor(weight,
2801  pip,
2802  pjp,
2803  pif);
2804  cs_centered_f_val_tensor(weight,
2805  pip,
2806  pjp,
2807  pjf);
2808 
2809  } else {
2810 
2811  /* Second order
2812  ------------*/
2813 
2814  cs_solu_f_val_tensor(cell_ceni,
2815  i_face_cog,
2816  gradi,
2817  pi,
2818  pif);
2819  cs_solu_f_val_tensor(cell_cenj,
2820  i_face_cog,
2821  gradj,
2822  pj,
2823  pjf);
2824 
2825  }
2826 
2827  /* Blending
2828  --------*/
2829 
2830  cs_blend_f_val_tensor(blencp,
2831  pi,
2832  pif);
2833  cs_blend_f_val_tensor(blencp,
2834  pj,
2835  pjf);
2836 
2837 }
2838 
2839 /*----------------------------------------------------------------------------*/
2883 /*----------------------------------------------------------------------------*/
2884 
2885 inline static void
2886 cs_i_cd_steady_slope_test(bool *upwind_switch,
2887  const int iconvp,
2888  const cs_real_t bldfrp,
2889  const int ischcp,
2890  const double relaxp,
2891  const double blencp,
2892  const double blend_st,
2893  const cs_real_t weight,
2894  const cs_real_t i_dist,
2895  const cs_real_t i_face_surf,
2896  const cs_real_3_t cell_ceni,
2897  const cs_real_3_t cell_cenj,
2898  const cs_real_3_t i_face_normal,
2899  const cs_real_3_t i_face_cog,
2900  const cs_real_3_t diipf,
2901  const cs_real_3_t djjpf,
2902  const cs_real_t i_massflux,
2903  const cs_real_3_t gradi,
2904  const cs_real_3_t gradj,
2905  const cs_real_3_t gradupi,
2906  const cs_real_3_t gradupj,
2907  const cs_real_3_t gradsti,
2908  const cs_real_3_t gradstj,
2909  const cs_real_t pi,
2910  const cs_real_t pj,
2911  const cs_real_t pia,
2912  const cs_real_t pja,
2913  cs_real_t *pifri,
2914  cs_real_t *pifrj,
2915  cs_real_t *pjfri,
2916  cs_real_t *pjfrj,
2917  cs_real_t *pip,
2918  cs_real_t *pjp,
2919  cs_real_t *pipr,
2920  cs_real_t *pjpr)
2921 {
2922  cs_real_t pir, pjr;
2923  cs_real_t recoi, recoj;
2924  cs_real_t testij, tesqck;
2925 
2926  *upwind_switch = false;
2927 
2928  cs_i_compute_quantities(bldfrp,
2929  diipf,
2930  djjpf,
2931  gradi,
2932  gradj,
2933  pi,
2934  pj,
2935  &recoi,
2936  &recoj,
2937  pip,
2938  pjp);
2939 
2940  cs_i_relax_c_val(relaxp,
2941  pia,
2942  pja,
2943  recoi,
2944  recoj,
2945  pi,
2946  pj,
2947  &pir,
2948  &pjr,
2949  pipr,
2950  pjpr);
2951 
2952  /* Convection slope test is needed only when iconv >0 */
2953  if (iconvp > 0) {
2954  cs_slope_test(pi,
2955  pj,
2956  i_dist,
2957  i_face_surf,
2958  i_face_normal,
2959  gradi,
2960  gradj,
2961  gradsti,
2962  gradstj,
2963  i_massflux,
2964  &testij,
2965  &tesqck);
2966 
2967  if (ischcp==1) {
2968 
2969  /* Centered
2970  --------*/
2971 
2972  cs_centered_f_val(weight,
2973  *pip,
2974  *pjpr,
2975  pifrj);
2976  cs_centered_f_val(weight,
2977  *pipr,
2978  *pjp,
2979  pifri);
2980  cs_centered_f_val(weight,
2981  *pipr,
2982  *pjp,
2983  pjfri);
2984  cs_centered_f_val(weight,
2985  *pip,
2986  *pjpr,
2987  pjfrj);
2988 
2989  } else if (ischcp == 0) {
2990 
2991  /* Second order
2992  ------------*/
2993 
2994  cs_solu_f_val(cell_ceni,
2995  i_face_cog,
2996  gradi,
2997  pi,
2998  pifrj);
2999  cs_solu_f_val(cell_ceni,
3000  i_face_cog,
3001  gradi,
3002  pir,
3003  pifri);
3004  cs_solu_f_val(cell_cenj,
3005  i_face_cog,
3006  gradj,
3007  pj,
3008  pjfri);
3009  cs_solu_f_val(cell_cenj,
3010  i_face_cog,
3011  gradj,
3012  pjr,
3013  pjfrj);
3014 
3015  } else {
3016 
3017  /* SOLU
3018  -----*/
3019 
3020  cs_solu_f_val(cell_ceni,
3021  i_face_cog,
3022  gradupi,
3023  pi,
3024  pifrj);
3025  cs_solu_f_val(cell_ceni,
3026  i_face_cog,
3027  gradupi,
3028  pir,
3029  pifri);
3030  cs_solu_f_val(cell_cenj,
3031  i_face_cog,
3032  gradupj,
3033  pj,
3034  pjfri);
3035  cs_solu_f_val(cell_cenj,
3036  i_face_cog,
3037  gradupj,
3038  pjr,
3039  pjfrj);
3040  }
3041 
3042 
3043  /* Slope test: Pourcentage of upwind
3044  ----------------------------------*/
3045 
3046  if (tesqck <= 0. || testij <= 0.) {
3047 
3048  cs_blend_f_val(blend_st,
3049  pi,
3050  pifrj);
3051  cs_blend_f_val(blend_st,
3052  pir,
3053  pifri);
3054  cs_blend_f_val(blend_st,
3055  pj,
3056  pjfri);
3057  cs_blend_f_val(blend_st,
3058  pjr,
3059  pjfrj);
3060 
3061  *upwind_switch = true;
3062 
3063  }
3064 
3065 
3066  /* Blending
3067  --------*/
3068 
3069  cs_blend_f_val(blencp,
3070  pi,
3071  pifrj);
3072  cs_blend_f_val(blencp,
3073  pir,
3074  pifri);
3075  cs_blend_f_val(blencp,
3076  pj,
3077  pjfri);
3078  cs_blend_f_val(blencp,
3079  pjr,
3080  pjfrj);
3081 
3082  /* If iconv=0 p*fr* are useless */
3083  } else {
3084  cs_upwind_f_val(pi,
3085  pifrj);
3086  cs_upwind_f_val(pir,
3087  pifri);
3088  cs_upwind_f_val(pj,
3089  pjfri);
3090  cs_upwind_f_val(pjr,
3091  pjfrj);
3092  }
3093 
3094 }
3095 
3096 /*----------------------------------------------------------------------------*/
3138 /*----------------------------------------------------------------------------*/
3139 
3140 inline static void
3141 cs_i_cd_steady_slope_test_vector(bool *upwind_switch,
3142  const int iconvp,
3143  const cs_real_t bldfrp,
3144  const int ischcp,
3145  const double relaxp,
3146  const double blencp,
3147  const double blend_st,
3148  const cs_real_t weight,
3149  const cs_real_t i_dist,
3150  const cs_real_t i_face_surf,
3151  const cs_real_3_t cell_ceni,
3152  const cs_real_3_t cell_cenj,
3153  const cs_real_3_t i_face_normal,
3154  const cs_real_3_t i_face_cog,
3155  const cs_real_3_t diipf,
3156  const cs_real_3_t djjpf,
3157  const cs_real_t i_massflux,
3158  const cs_real_33_t gradi,
3159  const cs_real_33_t gradj,
3160  const cs_real_33_t grdpai,
3161  const cs_real_33_t grdpaj,
3162  const cs_real_3_t pi,
3163  const cs_real_3_t pj,
3164  const cs_real_3_t pia,
3165  const cs_real_3_t pja,
3166  cs_real_t pifri[3],
3167  cs_real_t pifrj[3],
3168  cs_real_t pjfri[3],
3169  cs_real_t pjfrj[3],
3170  cs_real_t pip[3],
3171  cs_real_t pjp[3],
3172  cs_real_t pipr[3],
3173  cs_real_t pjpr[3])
3174 {
3175  cs_real_t pir[3], pjr[3];
3176  cs_real_t recoi[3], recoj[3];
3177  cs_real_t testij, tesqck;
3178  int isou;
3179 
3180  cs_i_compute_quantities_vector(bldfrp,
3181  diipf,
3182  djjpf,
3183  gradi,
3184  gradj,
3185  pi,
3186  pj,
3187  recoi,
3188  recoj,
3189  pip,
3190  pjp);
3191 
3192  cs_i_relax_c_val_vector(relaxp,
3193  pia,
3194  pja,
3195  recoi,
3196  recoj,
3197  pi,
3198  pj,
3199  pir,
3200  pjr,
3201  pipr,
3202  pjpr);
3203 
3204  /* Convection slope test is needed only when iconv >0 */
3205  if (iconvp > 0) {
3206  cs_slope_test_vector(pi,
3207  pj,
3208  i_dist,
3209  i_face_surf,
3210  i_face_normal,
3211  gradi,
3212  gradj,
3213  grdpai,
3214  grdpaj,
3215  i_massflux,
3216  &testij,
3217  &tesqck);
3218 
3219  for (isou = 0; isou < 3; isou++) {
3220  if (ischcp==1) {
3221 
3222  /* Centered
3223  --------*/
3224 
3225  cs_centered_f_val(weight,
3226  pip[isou],
3227  pjpr[isou],
3228  &pifrj[isou]);
3229  cs_centered_f_val(weight,
3230  pipr[isou],
3231  pjp[isou],
3232  &pifri[isou]);
3233  cs_centered_f_val(weight,
3234  pipr[isou],
3235  pjp[isou],
3236  &pjfri[isou]);
3237  cs_centered_f_val(weight,
3238  pip[isou],
3239  pjpr[isou],
3240  &pjfrj[isou]);
3241 
3242  } else {
3243 
3244  /* Second order
3245  ------------*/
3246 
3247  cs_solu_f_val(cell_ceni,
3248  i_face_cog,
3249  gradi[isou],
3250  pi[isou],
3251  &pifrj[isou]);
3252  cs_solu_f_val(cell_ceni,
3253  i_face_cog,
3254  gradi[isou],
3255  pir[isou],
3256  &pifri[isou]);
3257  cs_solu_f_val(cell_cenj,
3258  i_face_cog,
3259  gradj[isou],
3260  pj[isou],
3261  &pjfri[isou]);
3262  cs_solu_f_val(cell_cenj,
3263  i_face_cog,
3264  gradj[isou],
3265  pjr[isou],
3266  &pjfrj[isou]);
3267 
3268  }
3269 
3270  }
3271 
3272  /* Slope test: Pourcentage of upwind
3273  ----------------------------------*/
3274 
3275  if (tesqck <= 0. || testij <= 0.) {
3276  cs_blend_f_val_vector(blend_st,
3277  pi,
3278  pifrj);
3279  cs_blend_f_val_vector(blend_st,
3280  pir,
3281  pifri);
3282  cs_blend_f_val_vector(blend_st,
3283  pj,
3284  pjfri);
3285  cs_blend_f_val_vector(blend_st,
3286  pjr,
3287  pjfrj);
3288 
3289  *upwind_switch = true;
3290  }
3291 
3292 
3293  /* Blending
3294  --------*/
3295  cs_blend_f_val_vector(blencp,
3296  pi,
3297  pifrj);
3298  cs_blend_f_val_vector(blencp,
3299  pir,
3300  pifri);
3301  cs_blend_f_val_vector(blencp,
3302  pj,
3303  pjfri);
3304  cs_blend_f_val_vector(blencp,
3305  pjr,
3306  pjfrj);
3307 
3308  /* If iconv=0 p*fr* are useless */
3309  } else {
3310  for (isou = 0; isou < 3; isou++) {
3311  cs_upwind_f_val(pi[isou],
3312  &pifrj[isou]);
3313  cs_upwind_f_val(pir[isou],
3314  &pifri[isou]);
3315  cs_upwind_f_val(pj[isou],
3316  &pjfri[isou]);
3317  cs_upwind_f_val(pjr[isou],
3318  &pjfrj[isou]);
3319  }
3320  }
3321 
3322 }
3323 
3324 /*----------------------------------------------------------------------------*/
3366 /*----------------------------------------------------------------------------*/
3367 
3368 inline static void
3369 cs_i_cd_steady_slope_test_tensor(bool *upwind_switch,
3370  const int iconvp,
3371  const cs_real_t bldfrp,
3372  const int ischcp,
3373  const double relaxp,
3374  const double blencp,
3375  const double blend_st,
3376  const cs_real_t weight,
3377  const cs_real_t i_dist,
3378  const cs_real_t i_face_surf,
3379  const cs_real_3_t cell_ceni,
3380  const cs_real_3_t cell_cenj,
3381  const cs_real_3_t i_face_normal,
3382  const cs_real_3_t i_face_cog,
3383  const cs_real_3_t diipf,
3384  const cs_real_3_t djjpf,
3385  const cs_real_t i_massflux,
3386  const cs_real_63_t gradi,
3387  const cs_real_63_t gradj,
3388  const cs_real_63_t grdpai,
3389  const cs_real_63_t grdpaj,
3390  const cs_real_6_t pi,
3391  const cs_real_6_t pj,
3392  const cs_real_6_t pia,
3393  const cs_real_6_t pja,
3394  cs_real_t pifri[6],
3395  cs_real_t pifrj[6],
3396  cs_real_t pjfri[6],
3397  cs_real_t pjfrj[6],
3398  cs_real_t pip[6],
3399  cs_real_t pjp[6],
3400  cs_real_t pipr[6],
3401  cs_real_t pjpr[6])
3402 {
3403  cs_real_t pir[6], pjr[6];
3404  cs_real_t recoi[6], recoj[6];
3405  cs_real_t testij, tesqck;
3406  int isou;
3407 
3408  cs_i_compute_quantities_tensor(bldfrp,
3409  diipf,
3410  djjpf,
3411  gradi,
3412  gradj,
3413  pi,
3414  pj,
3415  recoi,
3416  recoj,
3417  pip,
3418  pjp);
3419 
3420  cs_i_relax_c_val_tensor(relaxp,
3421  pia,
3422  pja,
3423  recoi,
3424  recoj,
3425  pi,
3426  pj,
3427  pir,
3428  pjr,
3429  pipr,
3430  pjpr);
3431 
3432  /* Convection slope test is needed only when iconv >0 */
3433  if (iconvp > 0) {
3434  cs_slope_test_tensor(pi,
3435  pj,
3436  i_dist,
3437  i_face_surf,
3438  i_face_normal,
3439  gradi,
3440  gradj,
3441  grdpai,
3442  grdpaj,
3443  i_massflux,
3444  &testij,
3445  &tesqck);
3446 
3447  for (isou = 0; isou < 6; isou++) {
3448  if (ischcp==1) {
3449 
3450  /* Centered
3451  --------*/
3452 
3453  cs_centered_f_val(weight,
3454  pip[isou],
3455  pjpr[isou],
3456  &pifrj[isou]);
3457  cs_centered_f_val(weight,
3458  pipr[isou],
3459  pjp[isou],
3460  &pifri[isou]);
3461  cs_centered_f_val(weight,
3462  pipr[isou],
3463  pjp[isou],
3464  &pjfri[isou]);
3465  cs_centered_f_val(weight,
3466  pip[isou],
3467  pjpr[isou],
3468  &pjfrj[isou]);
3469 
3470  } else {
3471 
3472  /* Second order
3473  ------------*/
3474 
3475  cs_solu_f_val(cell_ceni,
3476  i_face_cog,
3477  gradi[isou],
3478  pi[isou],
3479  &pifrj[isou]);
3480  cs_solu_f_val(cell_ceni,
3481  i_face_cog,
3482  gradi[isou],
3483  pir[isou],
3484  &pifri[isou]);
3485  cs_solu_f_val(cell_cenj,
3486  i_face_cog,
3487  gradj[isou],
3488  pj[isou],
3489  &pjfri[isou]);
3490  cs_solu_f_val(cell_cenj,
3491  i_face_cog,
3492  gradj[isou],
3493  pjr[isou],
3494  &pjfrj[isou]);
3495 
3496  }
3497 
3498  }
3499 
3500  /* Slope test: Pourcentage of upwind
3501  ----------------------------------*/
3502 
3503  if (tesqck <= 0. || testij <= 0.) {
3504 
3505  cs_blend_f_val_tensor(blend_st,
3506  pi,
3507  pifrj);
3508  cs_blend_f_val_tensor(blend_st,
3509  pir,
3510  pifri);
3511  cs_blend_f_val_tensor(blend_st,
3512  pj,
3513  pjfri);
3514  cs_blend_f_val_tensor(blend_st,
3515  pjr,
3516  pjfrj);
3517 
3518  *upwind_switch = true;
3519 
3520  }
3521 
3522 
3523  /* Blending
3524  --------*/
3525 
3526  cs_blend_f_val_tensor(blencp,
3527  pi,
3528  pifrj);
3529  cs_blend_f_val_tensor(blencp,
3530  pir,
3531  pifri);
3532  cs_blend_f_val_tensor(blencp,
3533  pj,
3534  pjfri);
3535  cs_blend_f_val_tensor(blencp,
3536  pjr,
3537  pjfrj);
3538 
3539  /* If iconv=0 p*fr* are useless */
3540  } else {
3541  for (isou = 0; isou < 6; isou++) {
3542  cs_upwind_f_val(pi[isou],
3543  &pifrj[isou]);
3544  cs_upwind_f_val(pir[isou],
3545  &pifri[isou]);
3546  cs_upwind_f_val(pj[isou],
3547  &pjfri[isou]);
3548  cs_upwind_f_val(pjr[isou],
3549  &pjfrj[isou]);
3550  }
3551  }
3552 
3553 }
3554 
3555 /*----------------------------------------------------------------------------*/
3592 /*----------------------------------------------------------------------------*/
3593 
3594 inline static void
3595 cs_i_cd_unsteady_slope_test(bool *upwind_switch,
3596  const int iconvp,
3597  const cs_real_t bldfrp,
3598  const int ischcp,
3599  const double blencp,
3600  const double blend_st,
3601  const cs_real_t weight,
3602  const cs_real_t i_dist,
3603  const cs_real_t i_face_surf,
3604  const cs_real_3_t cell_ceni,
3605  const cs_real_3_t cell_cenj,
3606  const cs_real_3_t i_face_normal,
3607  const cs_real_3_t i_face_cog,
3608  const cs_real_3_t diipf,
3609  const cs_real_3_t djjpf,
3610  const cs_real_t i_massflux,
3611  const cs_real_3_t gradi,
3612  const cs_real_3_t gradj,
3613  const cs_real_3_t gradupi,
3614  const cs_real_3_t gradupj,
3615  const cs_real_3_t gradsti,
3616  const cs_real_3_t gradstj,
3617  const cs_real_t pi,
3618  const cs_real_t pj,
3619  cs_real_t *pif,
3620  cs_real_t *pjf,
3621  cs_real_t *pip,
3622  cs_real_t *pjp)
3623 {
3624  CS_UNUSED(blend_st);
3625 
3626  cs_real_t recoi, recoj;
3627  cs_real_t testij, tesqck;
3628 
3629  *upwind_switch = false;
3630 
3631  cs_i_compute_quantities(bldfrp,
3632  diipf,
3633  djjpf,
3634  gradi,
3635  gradj,
3636  pi,
3637  pj,
3638  &recoi,
3639  &recoj,
3640  pip,
3641  pjp);
3642 
3643  /* Convection slope test is needed only when iconv >0 */
3644  if (iconvp > 0) {
3645  cs_slope_test(pi,
3646  pj,
3647  i_dist,
3648  i_face_surf,
3649  i_face_normal,
3650  gradi,
3651  gradj,
3652  gradsti,
3653  gradstj,
3654  i_massflux,
3655  &testij,
3656  &tesqck);
3657 
3658  if (ischcp==1) {
3659 
3660  /* Centered
3661  --------*/
3662 
3663  cs_centered_f_val(weight,
3664  *pip,
3665  *pjp,
3666  pif);
3667  cs_centered_f_val(weight,
3668  *pip,
3669  *pjp,
3670  pjf);
3671 
3672  } else if (ischcp == 0) {
3673 
3674  /* Original SOLU
3675  --------------*/
3676 
3677  cs_solu_f_val(cell_ceni,
3678  i_face_cog,
3679  gradi,
3680  pi,
3681  pif);
3682  cs_solu_f_val(cell_cenj,
3683  i_face_cog,
3684  gradj,
3685  pj,
3686  pjf);
3687 
3688  } else {
3689 
3690  /* SOLU
3691  -----*/
3692 
3693  cs_solu_f_val(cell_ceni,
3694  i_face_cog,
3695  gradupi,
3696  pi,
3697  pif);
3698  cs_solu_f_val(cell_cenj,
3699  i_face_cog,
3700  gradupj,
3701  pj,
3702  pjf);
3703 
3704  }
3705 
3706  /* Slope test: Pourcentage of upwind
3707  ----------------------------------*/
3708 
3709  if (tesqck<=0. || testij<=0.) {
3710 
3711  cs_blend_f_val(blend_st,
3712  pi,
3713  pif);
3714  cs_blend_f_val(blend_st,
3715  pj,
3716  pjf);
3717 
3718  *upwind_switch = true;
3719 
3720  }
3721 
3722  /* Blending
3723  --------*/
3724 
3725  cs_blend_f_val(blencp,
3726  pi,
3727  pif);
3728  cs_blend_f_val(blencp,
3729  pj,
3730  pjf);
3731 
3732  /* If iconv=0 p*f are useless */
3733  } else {
3734  cs_upwind_f_val(pi,
3735  pif);
3736  cs_upwind_f_val(pj,
3737  pjf);
3738  }
3739 
3740 }
3741 
3742 /*----------------------------------------------------------------------------*/
3753 /*----------------------------------------------------------------------------*/
3754 
3755 inline static void
3756 cs_central_downwind_cells(const cs_lnum_t ii,
3757  const cs_lnum_t jj,
3758  const cs_real_t i_massflux,
3759  cs_lnum_t *ic,
3760  cs_lnum_t *id)
3761 {
3762  if (i_massflux >= 0.) {
3763  *ic = ii;
3764  *id = jj;
3765  } else {
3766  *ic = jj;
3767  *id = ii;
3768  }
3769 }
3770 
3771 /*----------------------------------------------------------------------------*/
3793 /*----------------------------------------------------------------------------*/
3794 
3795 inline static void
3796 cs_i_cd_unsteady_nvd(const cs_nvd_type_t limiter,
3797  const double beta,
3798  const cs_real_3_t cell_cen_c,
3799  const cs_real_3_t cell_cen_d,
3800  const cs_real_3_t i_face_normal,
3801  const cs_real_3_t i_face_cog,
3802  const cs_real_3_t gradv_c,
3803  const cs_real_t p_c,
3804  const cs_real_t p_d,
3805  const cs_real_t local_max_c,
3806  const cs_real_t local_min_c,
3807  const cs_real_t courant_c,
3808  cs_real_t *pif,
3809  cs_real_t *pjf)
3810 {
3811  /* distance between face center and central cell center */
3812  cs_real_t dist_fc;
3813  cs_real_3_t nfc;
3814  cs_math_3_length_unitv(cell_cen_c, i_face_cog, &dist_fc, nfc);
3815 
3816  /* unit vector and distance between central and downwind cells centers */
3817  cs_real_t dist_dc;
3818  cs_real_3_t ndc;
3819  cs_math_3_length_unitv(cell_cen_c, cell_cen_d, &dist_dc, ndc);
3820 
3821  /* Place the upwind point on the line that joins
3822  the two cells on the upwind side and the same
3823  distance as that between the two cells */
3824  const cs_real_t dist_cu = dist_dc;
3825  const cs_real_t dist_du = dist_dc + dist_cu;
3826 
3827  /* Compute the property on the upwind assuming a parabolic
3828  variation of the property between the two cells */
3829  const cs_real_t gradc = cs_math_3_dot_product(gradv_c, ndc);
3830 
3831  const cs_real_t grad2c = ((p_d - p_c)/dist_dc - gradc)/dist_dc;
3832 
3833  cs_real_t p_u = p_c + (grad2c*dist_cu - gradc)*dist_cu;
3834  p_u = CS_MAX(CS_MIN(p_u, local_max_c), local_min_c);
3835 
3836  /* Compute the normalised distances */
3837  const cs_real_t nvf_r_f = (dist_fc+dist_cu)/dist_du;
3838  const cs_real_t nvf_r_c = dist_cu/dist_du;
3839 
3840  /* Check for the bounds of the NVD diagram and compute the face
3841  property according to the selected NVD scheme */
3842  const cs_real_t _small = cs_math_epzero
3843  * (CS_ABS(p_u)+CS_ABS(p_c)+CS_ABS(p_d));
3844 
3845  if (CS_ABS(p_d-p_u) <= _small) {
3846  *pif = p_c;
3847  *pjf = p_c;
3848  } else {
3849  const cs_real_t nvf_p_c = (p_c - p_u)/(p_d - p_u);
3850 
3851  if (nvf_p_c <= 0. || nvf_p_c >= 1.) {
3852  *pif = p_c;
3853  *pjf = p_c;
3854  } else {
3855  cs_real_t nvf_p_f;
3856 
3857  /* Highly compressive NVD scheme for VOF */
3858  if (limiter >= CS_NVD_VOF_HRIC) {
3859  nvf_p_f = cs_nvd_vof_scheme_scalar(limiter,
3860  i_face_normal,
3861  nvf_p_c,
3862  nvf_r_f,
3863  nvf_r_c,
3864  gradv_c,
3865  courant_c);
3866  } else { /* Regular NVD scheme */
3867  nvf_p_f = cs_nvd_scheme_scalar(limiter,
3868  nvf_p_c,
3869  nvf_r_f,
3870  nvf_r_c);
3871  }
3872 
3873  *pif = p_u + nvf_p_f*(p_d - p_u);
3874  *pif = CS_MAX(CS_MIN(*pif, local_max_c), local_min_c);
3875 
3876  cs_blend_f_val(beta,
3877  p_c,
3878  pif);
3879 
3880  *pjf = *pif;
3881  }
3882  }
3883 }
3884 
3885 /*----------------------------------------------------------------------------*/
3920 /*----------------------------------------------------------------------------*/
3921 
3922 inline static void
3923 cs_i_cd_unsteady_slope_test_vector(bool *upwind_switch,
3924  const int iconvp,
3925  const cs_real_t bldfrp,
3926  const int ischcp,
3927  const double blencp,
3928  const double blend_st,
3929  const cs_real_t weight,
3930  const cs_real_t i_dist,
3931  const cs_real_t i_face_surf,
3932  const cs_real_3_t cell_ceni,
3933  const cs_real_3_t cell_cenj,
3934  const cs_real_3_t i_face_normal,
3935  const cs_real_3_t i_face_cog,
3936  const cs_real_3_t diipf,
3937  const cs_real_3_t djjpf,
3938  const cs_real_t i_massflux,
3939  const cs_real_33_t gradi,
3940  const cs_real_33_t gradj,
3941  const cs_real_33_t grdpai,
3942  const cs_real_33_t grdpaj,
3943  const cs_real_3_t pi,
3944  const cs_real_3_t pj,
3945  cs_real_t pif[3],
3946  cs_real_t pjf[3],
3947  cs_real_t pip[3],
3948  cs_real_t pjp[3])
3949 {
3950  cs_real_t recoi[3], recoj[3];
3951  cs_real_t testij, tesqck;
3952 
3953  cs_i_compute_quantities_vector(bldfrp,
3954  diipf,
3955  djjpf,
3956  gradi,
3957  gradj,
3958  pi,
3959  pj,
3960  recoi,
3961  recoj,
3962  pip,
3963  pjp);
3964 
3965  /* Convection slope test is needed only when iconv >0 */
3966  if (iconvp > 0) {
3967  cs_slope_test_vector(pi,
3968  pj,
3969  i_dist,
3970  i_face_surf,
3971  i_face_normal,
3972  gradi,
3973  gradj,
3974  grdpai,
3975  grdpaj,
3976  i_massflux,
3977  &testij,
3978  &tesqck);
3979 
3980  for (int isou = 0; isou < 3; isou++) {
3981  if (ischcp == 1) {
3982 
3983  /* Centered
3984  --------*/
3985 
3986  cs_centered_f_val(weight,
3987  pip[isou],
3988  pjp[isou],
3989  &pif[isou]);
3990  cs_centered_f_val(weight,
3991  pip[isou],
3992  pjp[isou],
3993  &pjf[isou]);
3994 
3995  } else {
3996 
3997  /* Second order
3998  ------------*/
3999 
4000  cs_solu_f_val(cell_ceni,
4001  i_face_cog,
4002  gradi[isou],
4003  pi[isou],
4004  &pif[isou]);
4005  cs_solu_f_val(cell_cenj,
4006  i_face_cog,
4007  gradj[isou],
4008  pj[isou],
4009  &pjf[isou]);
4010 
4011  }
4012 
4013  }
4014 
4015  /* Slope test: Pourcentage of upwind
4016  ----------------------------------*/
4017 
4018  if (tesqck <= 0. || testij <= 0.) {
4019 
4020  cs_blend_f_val_vector(blend_st,
4021  pi,
4022  pif);
4023  cs_blend_f_val_vector(blend_st,
4024  pj,
4025  pjf);
4026 
4027  *upwind_switch = true;
4028 
4029  }
4030 
4031 
4032  /* Blending
4033  --------*/
4034  cs_blend_f_val_vector(blencp,
4035  pi,
4036  pif);
4037  cs_blend_f_val_vector(blencp,
4038  pj,
4039  pjf);
4040 
4041  /* If iconv=0 p*f are useless */
4042  } else {
4043 
4044  for (int isou = 0; isou < 3; isou++) {
4045  cs_upwind_f_val(pi[isou],
4046  &pif[isou]);
4047  cs_upwind_f_val(pj[isou],
4048  &pjf[isou]);
4049 
4050  }
4051  }
4052 
4053 }
4054 
4055 /*----------------------------------------------------------------------------*/
4090 /*----------------------------------------------------------------------------*/
4091 
4092 inline static void
4093 cs_i_cd_unsteady_slope_test_tensor(bool *upwind_switch,
4094  const int iconvp,
4095  const cs_real_t bldfrp,
4096  const int ischcp,
4097  const double blencp,
4098  const double blend_st,
4099  const cs_real_t weight,
4100  const cs_real_t i_dist,
4101  const cs_real_t i_face_surf,
4102  const cs_real_3_t cell_ceni,
4103  const cs_real_3_t cell_cenj,
4104  const cs_real_3_t i_face_normal,
4105  const cs_real_3_t i_face_cog,
4106  const cs_real_3_t diipf,
4107  const cs_real_3_t djjpf,
4108  const cs_real_t i_massflux,
4109  const cs_real_63_t gradi,
4110  const cs_real_63_t gradj,
4111  const cs_real_63_t grdpai,
4112  const cs_real_63_t grdpaj,
4113  const cs_real_6_t pi,
4114  const cs_real_6_t pj,
4115  cs_real_t pif[6],
4116  cs_real_t pjf[6],
4117  cs_real_t pip[6],
4118  cs_real_t pjp[6])
4119 {
4120  cs_real_t recoi[6], recoj[6];
4121  cs_real_t testij, tesqck;
4122  int isou;
4123 
4124  cs_i_compute_quantities_tensor(bldfrp,
4125  diipf,
4126  djjpf,
4127  gradi,
4128  gradj,
4129  pi,
4130  pj,
4131  recoi,
4132  recoj,
4133  pip,
4134  pjp);
4135 
4136  /* Convection slope test is needed only when iconv >0 */
4137  if (iconvp > 0) {
4138  cs_slope_test_tensor(pi,
4139  pj,
4140  i_dist,
4141  i_face_surf,
4142  i_face_normal,
4143  gradi,
4144  gradj,
4145  grdpai,
4146  grdpaj,
4147  i_massflux,
4148  &testij,
4149  &tesqck);
4150 
4151  for (isou = 0; isou < 6; isou++) {
4152 
4153  if (ischcp==1) {
4154 
4155  /* Centered
4156  --------*/
4157 
4158  cs_centered_f_val(weight,
4159  pip[isou],
4160  pjp[isou],
4161  &pif[isou]);
4162  cs_centered_f_val(weight,
4163  pip[isou],
4164  pjp[isou],
4165  &pjf[isou]);
4166 
4167  } else {
4168 
4169  /* Second order
4170  ------------*/
4171 
4172  cs_solu_f_val(cell_ceni,
4173  i_face_cog,
4174  gradi[isou],
4175  pi[isou],
4176  &pif[isou]);
4177  cs_solu_f_val(cell_cenj,
4178  i_face_cog,
4179  gradj[isou],
4180  pj[isou],
4181  &pjf[isou]);
4182  }
4183 
4184  }
4185 
4186  /* Slope test activated: poucentage of upwind */
4187  if (tesqck <= 0. || testij <= 0.) {
4188 
4189  /* Upwind
4190  --------*/
4191 
4192  cs_blend_f_val_tensor(blend_st,
4193  pi,
4194  pif);
4195  cs_blend_f_val_tensor(blend_st,
4196  pj,
4197  pjf);
4198 
4199  *upwind_switch = true;
4200  }
4201 
4202 
4203  /* Blending
4204  --------*/
4205 
4206  cs_blend_f_val_tensor(blencp,
4207  pi,
4208  pif);
4209  cs_blend_f_val_tensor(blencp,
4210  pj,
4211  pjf);
4212 
4213  /* If iconv=0 p*fr* are useless */
4214  } else {
4215 
4216  for (isou = 0; isou < 6; isou++) {
4217  cs_upwind_f_val(pi[isou],
4218  &pif[isou]);
4219  cs_upwind_f_val(pj[isou],
4220  &pjf[isou]);
4221  }
4222  }
4223 }
4224 
4225 /*----------------------------------------------------------------------------*/
4234 /*----------------------------------------------------------------------------*/
4235 
4236 inline static void
4237 cs_b_compute_quantities(const cs_real_3_t diipb,
4238  const cs_real_3_t gradi,
4239  const cs_real_t bldfrp,
4240  cs_real_t *recoi)
4241 {
4242  *recoi = bldfrp * ( gradi[0]*diipb[0]
4243  + gradi[1]*diipb[1]
4244  + gradi[2]*diipb[2]);
4245 }
4246 
4247 /*----------------------------------------------------------------------------*/
4256 /*----------------------------------------------------------------------------*/
4257 
4258 inline static void
4259 cs_b_compute_quantities_vector(const cs_real_3_t diipb,
4260  const cs_real_33_t gradi,
4261  const cs_real_t bldfrp,
4262  cs_real_t recoi[3])
4263 {
4264  for (int isou = 0; isou < 3; isou++) {
4265  recoi[isou] = bldfrp * ( gradi[isou][0]*diipb[0]
4266  + gradi[isou][1]*diipb[1]
4267  + gradi[isou][2]*diipb[2]);
4268  }
4269 }
4270 
4271 /*----------------------------------------------------------------------------*/
4280 /*----------------------------------------------------------------------------*/
4281 
4282 inline static void
4283 cs_b_compute_quantities_tensor(const cs_real_3_t diipb,
4284  const cs_real_63_t gradi,
4285  const cs_real_t bldfrp,
4286  cs_real_t recoi[6])
4287 {
4288  for (int isou = 0; isou < 6; isou++) {
4289  recoi[isou] = bldfrp * ( gradi[isou][0]*diipb[0]
4290  + gradi[isou][1]*diipb[1]
4291  + gradi[isou][2]*diipb[2]);
4292  }
4293 }
4294 
4295 /*----------------------------------------------------------------------------*/
4306 /*----------------------------------------------------------------------------*/
4307 
4308 inline static void
4309 cs_b_relax_c_val(const double relaxp,
4310  const cs_real_t pi,
4311  const cs_real_t pia,
4312  const cs_real_t recoi,
4313  cs_real_t *pir,
4314  cs_real_t *pipr)
4315 {
4316  *pir = pi/relaxp - (1.-relaxp)/relaxp*pia;
4317  *pipr = *pir + recoi;
4318 }
4319 
4320 /*----------------------------------------------------------------------------*/
4331 /*----------------------------------------------------------------------------*/
4332 
4333 inline static void
4334 cs_b_relax_c_val_vector(const double relaxp,
4335  const cs_real_3_t pi,
4336  const cs_real_3_t pia,
4337  const cs_real_3_t recoi,
4338  cs_real_t pir[3],
4339  cs_real_t pipr[3])
4340 {
4341  for (int isou = 0; isou < 3; isou++) {
4342  pir[isou] = pi[isou]/relaxp - (1.-relaxp)/relaxp*pia[isou];
4343  pipr[isou] = pir[isou] + recoi[isou];
4344  }
4345 }
4346 
4347 /*----------------------------------------------------------------------------*/
4358 /*----------------------------------------------------------------------------*/
4359 
4360 inline static void
4361 cs_b_relax_c_val_tensor(const double relaxp,
4362  const cs_real_6_t pi,
4363  const cs_real_6_t pia,
4364  const cs_real_6_t recoi,
4365  cs_real_t pir[6],
4366  cs_real_t pipr[6])
4367 {
4368  for (int isou = 0; isou < 6; isou++) {
4369  pir[isou] = pi[isou]/relaxp - (1.-relaxp)/relaxp*pia[isou];
4370  pipr[isou] = pir[isou] + recoi[isou];
4371  }
4372 }
4373 
4374 /*----------------------------------------------------------------------------*/
4398 /*----------------------------------------------------------------------------*/
4399 
4400 inline static void
4401 cs_b_imposed_conv_flux(int iconvp,
4402  cs_real_t thetap,
4403  int imasac,
4404  int inc,
4405  int bc_type,
4406  int icvfli,
4407  cs_real_t pi,
4408  cs_real_t pir,
4409  cs_real_t pipr,
4410  cs_real_t coefap,
4411  cs_real_t coefbp,
4412  cs_real_t coface,
4413  cs_real_t cofbce,
4414  cs_real_t b_massflux,
4415  cs_real_t xcpp,
4416  cs_real_t *flux)
4417 {
4418  cs_real_t flui, fluj, pfac;
4419 
4420  /* Computed convective flux */
4421 
4422  if (icvfli == 0) {
4423 
4424  /* Remove decentering for coupled faces */
4425  if (bc_type == CS_COUPLED_FD) {
4426  flui = 0.0;
4427  fluj = b_massflux;
4428  } else {
4429  flui = 0.5*(b_massflux +fabs(b_massflux));
4430  fluj = 0.5*(b_massflux -fabs(b_massflux));
4431  }
4432 
4433  pfac = inc*coefap + coefbp*pipr;
4434  *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*pfac) -imasac*( b_massflux*pi));
4435 
4436  /* Imposed convective flux */
4437 
4438  } else {
4439 
4440  pfac = inc*coface + cofbce*pipr;
4441  *flux += iconvp*xcpp*(-imasac*(b_massflux*pi) + thetap*(pfac));
4442 
4443  }
4444 }
4445 
4446 /*----------------------------------------------------------------------------*/
4468 /*----------------------------------------------------------------------------*/
4469 
4470 inline static void
4471 cs_b_imposed_conv_flux_vector(int iconvp,
4472  cs_real_t thetap,
4473  int imasac,
4474  int inc,
4475  int bc_type,
4476  int icvfli,
4477  const cs_real_t pi[restrict 3],
4478  const cs_real_t pir[restrict 3],
4479  const cs_real_t pipr[restrict 3],
4480  const cs_real_t coefap[restrict 3],
4481  const cs_real_t coefbp[restrict 3][3],
4482  const cs_real_t coface[restrict 3],
4483  const cs_real_t cofbce[restrict 3][3],
4484  cs_real_t b_massflux,
4485  cs_real_t flux[restrict 3])
4486 {
4487  cs_real_t flui, fluj, pfac;
4488 
4489  /* Computed convective flux */
4490 
4491  if (icvfli == 0) {
4492 
4493  /* Remove decentering for coupled faces */
4494  if (bc_type == CS_COUPLED_FD) {
4495  flui = 0.0;
4496  fluj = b_massflux;
4497  } else {
4498  flui = 0.5*(b_massflux +fabs(b_massflux));
4499  fluj = 0.5*(b_massflux -fabs(b_massflux));
4500  }
4501  for (int isou = 0; isou < 3; isou++) {
4502  pfac = inc*coefap[isou];
4503  for (int jsou = 0; jsou < 3; jsou++) {
4504  pfac += coefbp[isou][jsou]*pipr[jsou];
4505  }
4506  flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac)
4507  - imasac*b_massflux*pi[isou]);
4508  }
4509 
4510  /* Imposed convective flux */
4511 
4512  } else {
4513 
4514  for (int isou = 0; isou < 3; isou++) {
4515  pfac = inc*coface[isou];
4516  for (int jsou = 0; jsou < 3; jsou++) {
4517  pfac += cofbce[isou][jsou]*pipr[jsou];
4518  }
4519  flux[isou] += iconvp*( thetap*pfac
4520  - imasac*b_massflux*pi[isou]);
4521  }
4522 
4523  }
4524 }
4525 
4526 /*----------------------------------------------------------------------------*/
4546 /*----------------------------------------------------------------------------*/
4547 
4548 inline static void
4549 cs_b_upwind_flux(const int iconvp,
4550  const cs_real_t thetap,
4551  const int imasac,
4552  const int inc,
4553  const int bc_type,
4554  const cs_real_t pi,
4555  const cs_real_t pir,
4556  const cs_real_t pipr,
4557  const cs_real_t coefap,
4558  const cs_real_t coefbp,
4559  const cs_real_t b_massflux,
4560  const cs_real_t xcpp,
4561  cs_real_t *flux)
4562 {
4563  cs_real_t flui, fluj, pfac;
4564 
4565  /* Remove decentering for coupled faces */
4566  if (bc_type == CS_COUPLED_FD) {
4567  flui = 0.0;
4568  fluj = b_massflux;
4569  } else {
4570  flui = 0.5*(b_massflux +fabs(b_massflux));
4571  fluj = 0.5*(b_massflux -fabs(b_massflux));
4572  }
4573 
4574  pfac = inc*coefap + coefbp*pipr;
4575  *flux += iconvp*xcpp*(thetap*(flui*pir + fluj*pfac) -imasac*( b_massflux*pi));
4576 }
4577 
4578 /*----------------------------------------------------------------------------*/
4598 /*----------------------------------------------------------------------------*/
4599 
4600 inline static void
4601 cs_b_upwind_flux_vector(const int iconvp,
4602  const cs_real_t thetap,
4603  const int imasac,
4604  const int inc,
4605  const int bc_type,
4606  const cs_real_3_t pi,
4607  const cs_real_3_t pir,
4608  const cs_real_3_t pipr,
4609  const cs_real_3_t coefa,
4610  const cs_real_33_t coefb,
4611  const cs_real_t b_massflux,
4612  cs_real_t flux[3])
4613 {
4614  cs_real_t flui, fluj, pfac;
4615 
4616  /* Remove decentering for coupled faces */
4617  if (bc_type == CS_COUPLED_FD) {
4618  flui = 0.0;
4619  fluj = b_massflux;
4620  } else {
4621  flui = 0.5*(b_massflux +fabs(b_massflux));
4622  fluj = 0.5*(b_massflux -fabs(b_massflux));
4623  }
4624  for (int isou = 0; isou < 3; isou++) {
4625  pfac = inc*coefa[isou];
4626  for (int jsou = 0; jsou < 3; jsou++) {
4627  pfac += coefb[isou][jsou]*pipr[jsou];
4628  }
4629  flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac)
4630  - imasac*b_massflux*pi[isou]);
4631  }
4632 }
4633 
4634 /*----------------------------------------------------------------------------*/
4654 /*----------------------------------------------------------------------------*/
4655 
4656 inline static void
4657 cs_b_upwind_flux_tensor(const int iconvp,
4658  const cs_real_t thetap,
4659  const int imasac,
4660  const int inc,
4661  const int bc_type,
4662  const cs_real_6_t pi,
4663  const cs_real_6_t pir,
4664  const cs_real_6_t pipr,
4665  const cs_real_6_t coefa,
4666  const cs_real_66_t coefb,
4667  const cs_real_t b_massflux,
4668  cs_real_t flux[6])
4669 {
4670  cs_real_t flui, fluj, pfac;
4671 
4672  /* Remove decentering for coupled faces */
4673  if (bc_type == CS_COUPLED_FD) {
4674  flui = 0.0;
4675  fluj = b_massflux;
4676  } else {
4677  flui = 0.5*(b_massflux +fabs(b_massflux));
4678  fluj = 0.5*(b_massflux -fabs(b_massflux));
4679  }
4680  for (int isou = 0; isou < 6; isou++) {
4681  pfac = inc*coefa[isou];
4682  for (int jsou = 0; jsou < 6; jsou++) {
4683  pfac += coefb[isou][jsou]*pipr[jsou];
4684  }
4685  flux[isou] += iconvp*( thetap*(flui*pir[isou] + fluj*pfac)
4686  - imasac*b_massflux*pi[isou]);
4687  }
4688 }
4689 
4690 /*----------------------------------------------------------------------------*/
4703 /*----------------------------------------------------------------------------*/
4704 
4705 inline static void
4706 cs_b_diff_flux(const int idiffp,
4707  const cs_real_t thetap,
4708  const int inc,
4709  const cs_real_t pipr,
4710  const cs_real_t cofafp,
4711  const cs_real_t cofbfp,
4712  const cs_real_t b_visc,
4713  cs_real_t *flux)
4714 {
4715  cs_real_t pfacd = inc*cofafp + cofbfp*pipr;
4716  *flux += idiffp*thetap*b_visc*pfacd;
4717 }
4718 
4719 /*----------------------------------------------------------------------------*/
4732 /*----------------------------------------------------------------------------*/
4733 
4734 inline static void
4735 cs_b_diff_flux_vector(const int idiffp,
4736  const cs_real_t thetap,
4737  const int inc,
4738  const cs_real_3_t pipr,
4739  const cs_real_3_t cofaf,
4740  const cs_real_33_t cofbf,
4741  const cs_real_t b_visc,
4742  cs_real_t flux[3])
4743 {
4744  cs_real_t pfacd ;
4745  for (int isou = 0; isou < 3; isou++) {
4746  pfacd = inc*cofaf[isou];
4747  for (int jsou = 0; jsou < 3; jsou++) {
4748  pfacd += cofbf[isou][jsou]*pipr[jsou];
4749  }
4750  flux[isou] += idiffp*thetap*b_visc*pfacd;
4751  }
4752 }
4753 
4754 /*----------------------------------------------------------------------------*/
4767 /*----------------------------------------------------------------------------*/
4768 
4769 inline static void
4770 cs_b_diff_flux_tensor(const int idiffp,
4771  const cs_real_t thetap,
4772  const int inc,
4773  const cs_real_6_t pipr,
4774  const cs_real_6_t cofaf,
4775  const cs_real_66_t cofbf,
4776  const cs_real_t b_visc,
4777  cs_real_t flux[6])
4778 {
4779  cs_real_t pfacd ;
4780  for (int isou = 0; isou < 6; isou++) {
4781  pfacd = inc*cofaf[isou];
4782  for (int jsou = 0; jsou < 6; jsou++) {
4783  pfacd += cofbf[isou][jsou]*pipr[jsou];
4784  }
4785  flux[isou] += idiffp*thetap*b_visc*pfacd;
4786  }
4787 }
4788 
4789 /*----------------------------------------------------------------------------*/
4803 /*----------------------------------------------------------------------------*/
4804 
4805 inline static void
4806 cs_b_cd_steady(const cs_real_t bldfrp,
4807  const double relaxp,
4808  const cs_real_3_t diipb,
4809  const cs_real_3_t gradi,
4810  const cs_real_t pi,
4811  const cs_real_t pia,
4812  cs_real_t *pir,
4813  cs_real_t *pipr)
4814 {
4815  cs_real_t recoi;
4816 
4817  cs_b_compute_quantities(diipb,
4818  gradi,
4819  bldfrp,
4820  &recoi);
4821 
4822  cs_b_relax_c_val(relaxp,
4823  pi,
4824  pia,
4825  recoi,
4826  pir,
4827  pipr);
4828 }
4829 
4830 /*----------------------------------------------------------------------------*/
4844 /*----------------------------------------------------------------------------*/
4845 
4846 inline static void
4847 cs_b_cd_steady_vector(const cs_real_t bldfrp,
4848  const double relaxp,
4849  const cs_real_3_t diipb,
4850  const cs_real_33_t gradi,
4851  const cs_real_3_t pi,
4852  const cs_real_3_t pia,
4853  cs_real_t pir[3],
4854  cs_real_t pipr[3])
4855 {
4856  cs_real_3_t recoi;
4857 
4858  cs_b_compute_quantities_vector(diipb,
4859  gradi,
4860  bldfrp,
4861  recoi);
4862 
4863  cs_b_relax_c_val_vector(relaxp,
4864  pi,
4865  pia,
4866  recoi,
4867  pir,
4868  pipr);
4869 }
4870 
4871 /*----------------------------------------------------------------------------*/
4885 /*----------------------------------------------------------------------------*/
4886 
4887 inline static void
4888 cs_b_cd_steady_tensor(const cs_real_t bldfrp,
4889  const double relaxp,
4890  const cs_real_3_t diipb,
4891  const cs_real_63_t gradi,
4892  const cs_real_6_t pi,
4893  const cs_real_6_t pia,
4894  cs_real_t pir[6],
4895  cs_real_t pipr[6])
4896 {
4897  cs_real_6_t recoi;
4898 
4899  cs_b_compute_quantities_tensor(diipb,
4900  gradi,
4901  bldfrp,
4902  recoi);
4903 
4904  cs_b_relax_c_val_tensor(relaxp,
4905  pi,
4906  pia,
4907  recoi,
4908  pir,
4909  pipr);
4910 }
4911 
4912 /*----------------------------------------------------------------------------*/
4923 /*----------------------------------------------------------------------------*/
4924 
4925 inline static void
4926 cs_b_cd_unsteady(const cs_real_t bldfrp,
4927  const cs_real_3_t diipb,
4928  const cs_real_3_t gradi,
4929  const cs_real_t pi,
4930  cs_real_t *pip)
4931 {
4932  cs_real_t recoi;
4933 
4934  cs_b_compute_quantities(diipb,
4935  gradi,
4936  bldfrp,
4937  &recoi);
4938 
4939  *pip = pi + recoi;
4940 }
4941 
4942 /*----------------------------------------------------------------------------*/
4953 /*----------------------------------------------------------------------------*/
4954 
4955 inline static void
4956 cs_b_cd_unsteady_vector(const cs_real_t bldfrp,
4957  const cs_real_3_t diipb,
4958  const cs_real_33_t gradi,
4959  const cs_real_3_t pi,
4960  cs_real_t pip[3])
4961 {
4962  cs_real_3_t recoi;
4963 
4964  cs_b_compute_quantities_vector(diipb,
4965  gradi,
4966  bldfrp,
4967  recoi);
4968 
4969  for (int isou = 0; isou < 3; isou++)
4970  pip[isou] = pi[isou] + recoi[isou];
4971 }
4972 
4973 /*----------------------------------------------------------------------------*/
4984 /*----------------------------------------------------------------------------*/
4985 
4986 inline static void
4987 cs_b_cd_unsteady_tensor(const cs_real_t bldfrp,
4988  const cs_real_3_t diipb,
4989  const cs_real_63_t gradi,
4990  const cs_real_6_t pi,
4991  cs_real_t pip[6])
4992 {
4993  cs_real_6_t recoi;
4994 
4995  cs_b_compute_quantities_tensor(diipb,
4996  gradi,
4997  bldfrp,
4998  recoi);
4999 
5000  for(int isou = 0; isou< 6; isou++)
5001  pip[isou] = pi[isou] + recoi[isou];
5002 }
5003 
5004 /*----------------------------------------------------------------------------*/
5015 /*----------------------------------------------------------------------------*/
5016 
5017 inline static void
5018 cs_b_diff_flux_coupling(int idiffp,
5019  cs_real_t pi,
5020  cs_real_t pj,
5021  cs_real_t b_visc,
5022  cs_real_t *fluxi)
5023 {
5024  *fluxi += idiffp*b_visc*(pi - pj);
5025 }
5026 
5027 /*----------------------------------------------------------------------------*/
5038 /*----------------------------------------------------------------------------*/
5039 
5040 inline static void
5041 cs_b_diff_flux_coupling_vector(int idiffp,
5042  const cs_real_t pi[3],
5043  const cs_real_t pj[3],
5044  cs_real_t b_visc,
5045  cs_real_t fluxi[3])
5046 {
5047  for (int k = 0; k < 3; k++)
5048  fluxi[k] += idiffp*b_visc*(pi[k] - pj[k]);
5049 }
5050 
5051 
5052 /*============================================================================
5053  * Public function prototypes for Fortran API
5054  *============================================================================*/
5055 
5056 /*----------------------------------------------------------------------------
5057  * Wrapper to cs_face_diffusion_potential
5058  *----------------------------------------------------------------------------*/
5059 
5060 void CS_PROCF (itrmas, ITRMAS)
5061 (
5062  const int *const f_id,
5063  const int *const init,
5064  const int *const inc,
5065  const int *const imrgra,
5066  const int *const iccocg,
5067  const int *const nswrgp,
5068  const int *const imligp,
5069  const int *const iphydp,
5070  const int *const iwgrp,
5071  const int *const iwarnp,
5072  const cs_real_t *const epsrgp,
5073  const cs_real_t *const climgp,
5074  const cs_real_t *const extrap,
5075  cs_real_3_t frcxt[],
5076  cs_real_t pvar[],
5077  const cs_real_t coefap[],
5078  const cs_real_t coefbp[],
5079  const cs_real_t cofafp[],
5080  const cs_real_t cofbfp[],
5081  const cs_real_t i_visc[],
5082  const cs_real_t b_visc[],
5083  cs_real_t visel[],
5084  cs_real_t i_massflux[],
5085  cs_real_t b_massflux[]
5086 );
5087 
5088 /*----------------------------------------------------------------------------
5089  * Wrapper to cs_face_anisotropic_diffusion_potential
5090  *----------------------------------------------------------------------------*/
5091 
5092 void CS_PROCF (itrmav, ITRMAV)
5093 (
5094  const int *const f_id,
5095  const int *const init,
5096  const int *const inc,
5097  const int *const imrgra,
5098  const int *const iccocg,
5099  const int *const nswrgp,
5100  const int *const imligp,
5101  const int *const ircflp,
5102  const int *const iphydp,
5103  const int *const iwgrp,
5104  const int *const iwarnp,
5105  const cs_real_t *const epsrgp,
5106  const cs_real_t *const climgp,
5107  const cs_real_t *const extrap,
5108  cs_real_3_t frcxt[],
5109  cs_real_t pvar[],
5110  const cs_real_t coefap[],
5111  const cs_real_t coefbp[],
5112  const cs_real_t cofafp[],
5113  const cs_real_t cofbfp[],
5114  const cs_real_t i_visc[],
5115  const cs_real_t b_visc[],
5116  cs_real_6_t viscel[],
5117  const cs_real_2_t weighf[],
5118  const cs_real_t weighb[],
5119  cs_real_t i_massflux[],
5120  cs_real_t b_massflux[]
5121 );
5122 
5123 /*----------------------------------------------------------------------------
5124  * Wrapper to cs_diffusion_potential
5125  *----------------------------------------------------------------------------*/
5126 
5127 void CS_PROCF (itrgrp, ITRGRP)
5128 (
5129  const int *const f_id,
5130  const int *const init,
5131  const int *const inc,
5132  const int *const imrgra,
5133  const int *const iccocg,
5134  const int *const nswrgp,
5135  const int *const imligp,
5136  const int *const iphydp,
5137  const int *const iwgrp,
5138  const int *const iwarnp,
5139  const cs_real_t *const epsrgp,
5140  const cs_real_t *const climgp,
5141  const cs_real_t *const extrap,
5142  cs_real_3_t frcxt[],
5143  cs_real_t pvar[],
5144  const cs_real_t coefap[],
5145  const cs_real_t coefbp[],
5146  const cs_real_t cofafp[],
5147  const cs_real_t cofbfp[],
5148  const cs_real_t i_visc[],
5149  const cs_real_t b_visc[],
5150  cs_real_t visel[],
5151  cs_real_t diverg[]
5152 );
5153 
5154 /*----------------------------------------------------------------------------
5155  * Wrapper to cs_anisotropic_diffusion_potential
5156  *----------------------------------------------------------------------------*/
5157 
5158 void CS_PROCF (itrgrv, ITRGRV)
5159 (
5160  const int *const f_id,
5161  const int *const init,
5162  const int *const inc,
5163  const int *const imrgra,
5164  const int *const iccocg,
5165  const int *const nswrgp,
5166  const int *const imligp,
5167  const int *const ircflp,
5168  const int *const iphydp,
5169  const int *const iwgrp,
5170  const int *const iwarnp,
5171  const cs_real_t *const epsrgp,
5172  const cs_real_t *const climgp,
5173  const cs_real_t *const extrap,
5174  cs_real_3_t frcxt[],
5175  cs_real_t pvar[],
5176  const cs_real_t coefap[],
5177  const cs_real_t coefbp[],
5178  const cs_real_t cofafp[],
5179  const cs_real_t cofbfp[],
5180  const cs_real_t i_visc[],
5181  const cs_real_t b_visc[],
5182  cs_real_6_t viscel[],
5183  const cs_real_2_t weighf[],
5184  const cs_real_t weighb[],
5185  cs_real_t diverg[]
5186 );
5187 
5188 /*=============================================================================
5189  * Public function prototypes
5190  *============================================================================*/
5191 
5192 /*----------------------------------------------------------------------------*/
5211 /*----------------------------------------------------------------------------*/
5212 
5213 void
5214 cs_slope_test_gradient(int f_id,
5215  int inc,
5216  cs_halo_type_t halo_type,
5217  const cs_real_3_t *grad,
5218  cs_real_3_t *grdpa,
5219  const cs_real_t *pvar,
5220  const cs_real_t *coefap,
5221  const cs_real_t *coefbp,
5222  const cs_real_t *i_massflux);
5223 
5224 /*----------------------------------------------------------------------------*/
5240 /*----------------------------------------------------------------------------*/
5241 
5242 void
5243 cs_upwind_gradient(const int f_id,
5244  const int inc,
5245  const cs_halo_type_t halo_type,
5246  const cs_real_t coefap[],
5247  const cs_real_t coefbp[],
5248  const cs_real_t i_massflux[],
5249  const cs_real_t b_massflux[],
5250  const cs_real_t *restrict pvar,
5251  cs_real_3_t *restrict grdpa);
5252 
5253 /*----------------------------------------------------------------------------*/
5271 /*----------------------------------------------------------------------------*/
5272 
5273 void
5274 cs_slope_test_gradient_vector(const int inc,
5275  const cs_halo_type_t halo_type,
5276  const cs_real_33_t *grad,
5277  cs_real_33_t *grdpa,
5278  const cs_real_3_t *pvar,
5279  const cs_real_3_t *coefa,
5280  const cs_real_33_t *coefb,
5281  const cs_real_t *i_massflux);
5282 
5283 /*----------------------------------------------------------------------------*/
5301 /*----------------------------------------------------------------------------*/
5302 
5303 void
5304 cs_slope_test_gradient_tensor(const int inc,
5305  const cs_halo_type_t halo_type,
5306  const cs_real_63_t *grad,
5307  cs_real_63_t *grdpa,
5308  const cs_real_6_t *pvar,
5309  const cs_real_6_t *coefa,
5310  const cs_real_66_t *coefb,
5311  const cs_real_t *i_massflux);
5312 
5313 /*----------------------------------------------------------------------------*/
5322 /*----------------------------------------------------------------------------*/
5323 
5324 void
5325 cs_beta_limiter_building(int f_id,
5326  int inc,
5327  const cs_real_t rovsdt[]);
5328 
5329 /*----------------------------------------------------------------------------*/
5381 /*----------------------------------------------------------------------------*/
5382 
5383 void
5384 cs_convection_diffusion_scalar(int idtvar,
5385  int f_id,
5386  const cs_var_cal_opt_t var_cal_opt,
5387  int icvflb,
5388  int inc,
5389  int iccocg,
5390  int imasac,
5391  cs_real_t *restrict pvar,
5392  const cs_real_t *restrict pvara,
5393  const int icvfli[],
5394  const cs_real_t coefap[],
5395  const cs_real_t coefbp[],
5396  const cs_real_t cofafp[],
5397  const cs_real_t cofbfp[],
5398  const cs_real_t i_massflux[],
5399  const cs_real_t b_massflux[],
5400  const cs_real_t i_visc[],
5401  const cs_real_t b_visc[],
5402  cs_real_t *restrict rhs);
5403 
5404 /*----------------------------------------------------------------------------*/
5443 /*----------------------------------------------------------------------------*/
5444 
5445 void
5446 cs_face_convection_scalar(int idtvar,
5447  int f_id,
5448  const cs_var_cal_opt_t var_cal_opt,
5449  int icvflb,
5450  int inc,
5451  int iccocg,
5452  int imasac,
5453  cs_real_t *restrict pvar,
5454  const cs_real_t *restrict pvara,
5455  const int icvfli[],
5456  const cs_real_t coefap[],
5457  const cs_real_t coefbp[],
5458  const cs_real_t i_massflux[],
5459  const cs_real_t b_massflux[],
5460  cs_real_2_t i_conv_flux[],
5461  cs_real_t b_conv_flux[]);
5462 
5463 /*----------------------------------------------------------------------------*/
5523 /*----------------------------------------------------------------------------*/
5524 
5525 void
5526 cs_convection_diffusion_vector(int idtvar,
5527  int f_id,
5528  const cs_var_cal_opt_t var_cal_opt,
5529  int icvflb,
5530  int inc,
5531  int ivisep,
5532  int imasac,
5533  cs_real_3_t *restrict pvar,
5534  const cs_real_3_t *restrict pvara,
5535  const int icvfli[],
5536  const cs_real_3_t coefav[],
5537  const cs_real_33_t coefbv[],
5538  const cs_real_3_t cofafv[],
5539  const cs_real_33_t cofbfv[],
5540  const cs_real_t i_massflux[],
5541  const cs_real_t b_massflux[],
5542  const cs_real_t i_visc[],
5543  const cs_real_t b_visc[],
5544  const cs_real_t i_secvis[],
5545  const cs_real_t b_secvis[],
5546  cs_real_3_t *restrict rhs);
5547 
5548 /*----------------------------------------------------------------------------*/
5593 /*----------------------------------------------------------------------------*/
5594 
5595 void
5596 cs_convection_diffusion_tensor(int idtvar,
5597  int f_id,
5598  const cs_var_cal_opt_t var_cal_opt,
5599  int icvflb,
5600  int inc,
5601  int imasac,
5602  cs_real_6_t *restrict pvar,
5603  const cs_real_6_t *restrict pvara,
5604  const cs_real_6_t coefa[],
5605  const cs_real_66_t coefb[],
5606  const cs_real_6_t cofaf[],
5607  const cs_real_66_t cofbf[],
5608  const cs_real_t i_massflux[],
5609  const cs_real_t b_massflux[],
5610  const cs_real_t i_visc[],
5611  const cs_real_t b_visc[],
5612  cs_real_6_t *restrict rhs);
5613 
5614 /*----------------------------------------------------------------------------*/
5660 /*----------------------------------------------------------------------------*/
5661 
5662 void
5663 cs_convection_diffusion_thermal(int idtvar,
5664  int f_id,
5665  const cs_var_cal_opt_t var_cal_opt,
5666  int inc,
5667  int iccocg,
5668  int imasac,
5669  cs_real_t *restrict pvar,
5670  const cs_real_t *restrict pvara,
5671  const cs_real_t coefap[],
5672  const cs_real_t coefbp[],
5673  const cs_real_t cofafp[],
5674  const cs_real_t cofbfp[],
5675  const cs_real_t i_massflux[],
5676  const cs_real_t b_massflux[],
5677  const cs_real_t i_visc[],
5678  const cs_real_t b_visc[],
5679  const cs_real_t xcpp[],
5680  cs_real_t *restrict rhs);
5681 
5682 /*----------------------------------------------------------------------------*/
5730 /*----------------------------------------------------------------------------*/
5731 
5732 void
5733 cs_anisotropic_diffusion_scalar(int idtvar,
5734  int f_id,
5735  const cs_var_cal_opt_t var_cal_opt,
5736  int inc,
5737  int iccocg,
5738  cs_real_t *restrict pvar,
5739  const cs_real_t *restrict pvara,
5740  const cs_real_t coefap[],
5741  const cs_real_t coefbp[],
5742  const cs_real_t cofafp[],
5743  const cs_real_t cofbfp[],
5744  const cs_real_t i_visc[],
5745  const cs_real_t b_visc[],
5746  cs_real_6_t *restrict viscel,
5747  const cs_real_2_t weighf[],
5748  const cs_real_t weighb[],
5749  cs_real_t *restrict rhs);
5750 
5751 /*-----------------------------------------------------------------------------*/
5801 /*----------------------------------------------------------------------------*/
5802 
5803 void
5804 cs_anisotropic_left_diffusion_vector(int idtvar,
5805  int f_id,
5806  const cs_var_cal_opt_t var_cal_opt,
5807  int inc,
5808  int ivisep,
5809  cs_real_3_t *restrict pvar,
5810  const cs_real_3_t *restrict pvara,
5811  const cs_real_3_t coefav[],
5812  const cs_real_33_t coefbv[],
5813  const cs_real_3_t cofafv[],
5814  const cs_real_33_t cofbfv[],
5815  const cs_real_33_t i_visc[],
5816  const cs_real_t b_visc[],
5817  const cs_real_t i_secvis[],
5818  cs_real_3_t *restrict rhs);
5819 
5820 /*-----------------------------------------------------------------------------*/
5865 /*----------------------------------------------------------------------------*/
5866 
5867 void
5868 cs_anisotropic_right_diffusion_vector(int idtvar,
5869  int f_id,
5870  const cs_var_cal_opt_t var_cal_opt,
5871  int inc,
5872  cs_real_3_t *restrict pvar,
5873  const cs_real_3_t *restrict pvara,
5874  const cs_real_3_t coefav[],
5875  const cs_real_33_t coefbv[],
5876  const cs_real_3_t cofafv[],
5877  const cs_real_33_t cofbfv[],
5878  const cs_real_t i_visc[],
5879  const cs_real_t b_visc[],
5880  cs_real_6_t *restrict viscel,
5881  const cs_real_2_t weighf[],
5882  const cs_real_t weighb[],
5883  cs_real_3_t *restrict rhs);
5884 
5885 /*----------------------------------------------------------------------------*/
5929 /*----------------------------------------------------------------------------*/
5930 
5931 void
5932 cs_anisotropic_diffusion_tensor(int idtvar,
5933  int f_id,
5934  const cs_var_cal_opt_t var_cal_opt,
5935  int inc,
5936  cs_real_6_t *restrict pvar,
5937  const cs_real_6_t *restrict pvara,
5938  const cs_real_6_t coefa[],
5939  const cs_real_66_t coefb[],
5940  const cs_real_6_t cofaf[],
5941  const cs_real_66_t cofbf[],
5942  const cs_real_t i_visc[],
5943  const cs_real_t b_visc[],
5944  cs_real_6_t *restrict viscel,
5945  const cs_real_2_t weighf[],
5946  const cs_real_t weighb[],
5947  cs_real_6_t *restrict rhs);
5948 
5949 /*----------------------------------------------------------------------------*/
6008 /*----------------------------------------------------------------------------*/
6009 
6010 void
6011 cs_face_diffusion_potential(const int f_id,
6012  const cs_mesh_t *m,
6013  cs_mesh_quantities_t *fvq,
6014  int init,
6015  int inc,
6016  int imrgra,
6017  int iccocg,
6018  int nswrgp,
6019  int imligp,
6020  int iphydp,
6021  int iwgrp,
6022  int iwarnp,
6023  double epsrgp,
6024  double climgp,
6025  cs_real_3_t *restrict frcxt,
6026  cs_real_t *restrict pvar,
6027  const cs_real_t coefap[],
6028  const cs_real_t coefbp[],
6029  const cs_real_t cofafp[],
6030  const cs_real_t cofbfp[],
6031  const cs_real_t i_visc[],
6032  const cs_real_t b_visc[],
6033  cs_real_t *restrict visel,
6034  cs_real_t *restrict i_massflux,
6035  cs_real_t *restrict b_massflux);
6036 
6037 /*----------------------------------------------------------------------------*/
6107 /*----------------------------------------------------------------------------*/
6108 
6109 void
6110 cs_face_anisotropic_diffusion_potential(const int f_id,
6111  const cs_mesh_t *m,
6112  cs_mesh_quantities_t *fvq,
6113  int init,
6114  int inc,
6115  int imrgra,
6116  int iccocg,
6117  int nswrgp,
6118  int imligp,
6119  int ircflp,
6120  int iphydp,
6121  int iwgrp,
6122  int iwarnp,
6123  double epsrgp,
6124  double climgp,
6125  cs_real_3_t *restrict frcxt,
6126  cs_real_t *restrict pvar,
6127  const cs_real_t coefap[],
6128  const cs_real_t coefbp[],
6129  const cs_real_t cofafp[],
6130  const cs_real_t cofbfp[],
6131  const cs_real_t i_visc[],
6132  const cs_real_t b_visc[],
6133  cs_real_6_t *restrict viscel,
6134  const cs_real_2_t weighf[],
6135  const cs_real_t weighb[],
6136  cs_real_t *restrict i_massflux,
6137  cs_real_t *restrict b_massflux);
6138 
6139 /*----------------------------------------------------------------------------*/
6197 /*----------------------------------------------------------------------------*/
6198 
6199 void
6200 cs_diffusion_potential(const int f_id,
6201  const cs_mesh_t *m,
6202  cs_mesh_quantities_t *fvq,
6203  int init,
6204  int inc,
6205  int imrgra,
6206  int iccocg,
6207  int nswrgp,
6208  int imligp,
6209  int iphydp,
6210  int iwgrp,
6211  int iwarnp,
6212  double epsrgp,
6213  double climgp,
6214  cs_real_3_t *restrict frcxt,
6215  cs_real_t *restrict pvar,
6216  const cs_real_t coefap[],
6217  const cs_real_t coefbp[],
6218  const cs_real_t cofafp[],
6219  const cs_real_t cofbfp[],
6220  const cs_real_t i_visc[],
6221  const cs_real_t b_visc[],
6222  cs_real_t visel[],
6223  cs_real_t *restrict diverg);
6224 
6225 /*----------------------------------------------------------------------------*/
6296 /*----------------------------------------------------------------------------*/
6297 
6298 void
6299 cs_anisotropic_diffusion_potential(const int f_id,
6300  const cs_mesh_t *m,
6301  cs_mesh_quantities_t *fvq,
6302  int init,
6303  int inc,
6304  int imrgra,
6305  int iccocg,
6306  int nswrgp,
6307  int imligp,
6308  int ircflp,
6309  int iphydp,
6310  int iwgrp,
6311  int iwarnp,
6312  double epsrgp,
6313  double climgp,
6314  cs_real_3_t *restrict frcxt,
6315  cs_real_t *restrict pvar,
6316  const cs_real_t coefap[],
6317  const cs_real_t coefbp[],
6318  const cs_real_t cofafp[],
6319  const cs_real_t cofbfp[],
6320  const cs_real_t i_visc[],
6321  const cs_real_t b_visc[],
6322  cs_real_6_t *restrict viscel,
6323  const cs_real_2_t weighf[],
6324  const cs_real_t weighb[],
6325  cs_real_t *restrict diverg);
6326 
6327 /*----------------------------------------------------------------------------*/
6328 
6329 END_C_DECLS
6330 
6331 #endif /* __CS_CONVECTION_DIFFUSION_H__ */
k
Definition: cs_field_pointer.h:70
itrgrv
void itrgrv(const int *const f_id, const int *const init, const int *const inc, const int *const imrgra, const int *const iccocg, const int *const nswrgp, const int *const imligp, const int *const ircflp, const int *const iphydp, const int *const iwgrp, const int *const iwarnp, const cs_real_t *const epsrgp, const cs_real_t *const climgp, const cs_real_t *const extrap, cs_real_3_t frcxt[], cs_real_t pvar[], const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t viscel[], const cs_real_2_t weighf[], const cs_real_t weighb[], cs_real_t diverg[])
Definition: cs_convection_diffusion.c:807
cs_convection_diffusion_vector
void cs_convection_diffusion_vector(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int icvflb, int inc, int ivisep, int imasac, cs_real_3_t *restrict pvar, const cs_real_3_t *restrict pvara, const int icvfli[], const cs_real_3_t coefav[], const cs_real_33_t coefbv[], const cs_real_3_t cofafv[], const cs_real_33_t cofbfv[], 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 rhs)
Add the explicit part of the convection/diffusion terms of a transport equation of a vector field ...
Definition: cs_convection_diffusion.c:4174
mesh::diipb
double precision, dimension(:,:), pointer diipb
Definition: mesh.f90:212
restrict
#define restrict
Definition: cs_defs.h:127
cs_real_2_t
cs_real_t cs_real_2_t[2]
vector of 2 floating-point values
Definition: cs_defs.h:319
cs_base.h
cs_real_6_t
cs_real_t cs_real_6_t[6]
vector of 6 floating-point values
Definition: cs_defs.h:322
cs_field_key_id_try
int cs_field_key_id_try(const char *name)
Return an id associated with a given key name if present.
Definition: cs_field.c:2524
cs_anisotropic_diffusion_tensor
void cs_anisotropic_diffusion_tensor(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int inc, cs_real_6_t *restrict pvar, const cs_real_6_t *restrict pvara, const cs_real_6_t coefa[], const cs_real_66_t coefb[], const cs_real_6_t cofaf[], const cs_real_66_t cofbf[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t *restrict viscel, const cs_real_2_t weighf[], const cs_real_t weighb[], cs_real_6_t *restrict rhs)
Add the explicit part of the diffusion terms with a symmetric tensor diffusivity for a transport equa...
Definition: cs_convection_diffusion.c:10184
cs_nvd_type_t
cs_nvd_type_t
Definition: cs_convection_diffusion.h:59
cs_face_diffusion_potential
void cs_face_diffusion_potential(const int f_id, const cs_mesh_t *m, cs_mesh_quantities_t *fvq, int init, int inc, int imrgra, int iccocg, int nswrgp, int imligp, int iphydp, int iwgrp, int iwarnp, double epsrgp, double climgp, cs_real_3_t *restrict frcxt, cs_real_t *restrict pvar, const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_t *restrict visel, cs_real_t *restrict i_massflux, cs_real_t *restrict b_massflux)
Update the face mass flux with the face pressure (or pressure increment, or pressure double increment...
Definition: cs_convection_diffusion.c:10801
cs_anisotropic_diffusion_scalar
void cs_anisotropic_diffusion_scalar(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int inc, int iccocg, cs_real_t *restrict pvar, const cs_real_t *restrict pvara, const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t *restrict viscel, const cs_real_2_t weighf[], const cs_real_t weighb[], cs_real_t *restrict rhs)
Add the explicit part of the diffusion terms with a symmetric tensor diffusivity for a transport equa...
Definition: cs_convection_diffusion.c:8168
cs_field_t
Field descriptor.
Definition: cs_field.h:125
cs_equation_param_t
Set of parameters to handle an unsteady convection-diffusion-reaction equation with term sources...
Definition: cs_equation_param.h:202
CS_ABS
#define CS_ABS(a)
Definition: cs_defs.h:457
cs_field_by_id
cs_field_t * cs_field_by_id(int id)
Return a pointer to a field based on its id.
Definition: cs_field.c:2314
cs_real_66_t
cs_real_t cs_real_66_t[6][6]
6x6 matrix of floating-point values
Definition: cs_defs.h:328
itrgrp
void itrgrp(const int *const f_id, const int *const init, const int *const inc, const int *const imrgra, const int *const iccocg, const int *const nswrgp, const int *const imligp, const int *const iphydp, const int *const iwgrp, const int *const iwarnp, const cs_real_t *const epsrgp, const cs_real_t *const climgp, const cs_real_t *const extrap, cs_real_3_t frcxt[], cs_real_t pvar[], const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_t visel[], cs_real_t diverg[])
Definition: cs_convection_diffusion.c:747
CS_NVD_STOIC
Definition: cs_convection_diffusion.h:68
cs_slope_test_gradient_tensor
void cs_slope_test_gradient_tensor(const int inc, const cs_halo_type_t halo_type, const cs_real_63_t *grad, cs_real_63_t *grdpa, const cs_real_6_t *pvar, const cs_real_6_t *coefa, const cs_real_66_t *coefb, const cs_real_t *i_massflux)
Compute the upwind gradient used in the slope tests.
Definition: cs_convection_diffusion.c:1325
cstnum::pi
double precision pi
value with 16 digits
Definition: cstnum.f90:48
cs_defs.h
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:495
cs_face_convection_scalar
void cs_face_convection_scalar(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int icvflb, int inc, int iccocg, int imasac, cs_real_t *restrict pvar, const cs_real_t *restrict pvara, const int icvfli[], const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t i_massflux[], const cs_real_t b_massflux[], cs_real_2_t i_conv_flux[], cs_real_t b_conv_flux[])
Update face flux with convection contribution of a standard transport equation of a scalar field ...
Definition: cs_convection_diffusion.c:3060
cs_math_epzero
const cs_real_t cs_math_epzero
cs_glob_mesh
cs_mesh_t * cs_glob_mesh
CS_UNUSED
#define CS_UNUSED(x)
Definition: cs_defs.h:481
cs_field_pointer.h
cs_slope_test_gradient_vector
void cs_slope_test_gradient_vector(const int inc, const cs_halo_type_t halo_type, const cs_real_33_t *grad, cs_real_33_t *grdpa, const cs_real_3_t *pvar, const cs_real_3_t *coefa, const cs_real_33_t *coefb, const cs_real_t *i_massflux)
Compute the upwind gradient used in the slope tests.
Definition: cs_convection_diffusion.c:1170
cpincl::beta
double precision, dimension(ncharm), save beta
Definition: cpincl.f90:99
optcal::idtvar
integer(c_int), pointer, save idtvar
option for a variable time step
Definition: optcal.f90:374
cs_anisotropic_left_diffusion_vector
void cs_anisotropic_left_diffusion_vector(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int inc, int ivisep, cs_real_3_t *restrict pvar, const cs_real_3_t *restrict pvara, const cs_real_3_t coefav[], const cs_real_33_t coefbv[], const cs_real_3_t cofafv[], const cs_real_33_t cofbfv[], const cs_real_33_t i_visc[], const cs_real_t b_visc[], const cs_real_t i_secvis[], cs_real_3_t *restrict rhs)
Add explicit part of the terms of diffusion by a left-multiplying symmetric tensorial diffusivity for...
Definition: cs_convection_diffusion.c:8923
p
Definition: cs_field_pointer.h:67
cs_math_3_length_unitv
void cs_math_3_length_unitv(const cs_real_t xa[3], const cs_real_t xb[3], cs_real_t *len, cs_real_3_t unitv)
Compute the length (euclidien norm) between two points xa and xb in a cartesian coordinate system of ...
Definition: cs_math.c:423
CS_NVD_CUBISTA
Definition: cs_convection_diffusion.h:63
cs_field_t::val
cs_real_t * val
Definition: cs_field.h:146
CS_F_
#define CS_F_(e)
Macro used to return a field pointer by its enumerated value.
Definition: cs_field_pointer.h:51
cs_parameters.h
CS_NVD_MINMOD
Definition: cs_convection_diffusion.h:66
cfpoin::icvfli
integer, dimension(:), allocatable icvfli
boundary convection flux indicator of a Rusanov or an analytical flux (some boundary contributions of...
Definition: cfpoin.f90:52
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:307
itrmas
void itrmas(const int *const f_id, const int *const init, const int *const inc, const int *const imrgra, const int *const iccocg, const int *const nswrgp, const int *const imligp, const int *const iphydp, const int *const iwgrp, const int *const iwarnp, const cs_real_t *const epsrgp, const cs_real_t *const climgp, const cs_real_t *const extrap, cs_real_3_t frcxt[], cs_real_t pvar[], const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_t visel[], cs_real_t i_massflux[], cs_real_t b_massflux[])
Definition: cs_convection_diffusion.c:617
cs_mesh_t
Definition: cs_mesh.h:84
cs_field_get_key_int
int cs_field_get_key_int(const cs_field_t *f, int key_id)
Return a integer value for a given key associated with a field.
Definition: cs_field.c:2991
cs_field_key_id
int cs_field_key_id(const char *name)
Return an id associated with a given key name.
Definition: cs_field.c:2497
var_cal_opt
cs_equation_param_t * var_cal_opt
Definition: keywords.h:135
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.c:1434
cs_mesh_quantities.h
cpincl::b2
double precision, dimension(ncharm), save b2
Definition: cpincl.f90:233
cs_equation_param_t::iconv
int iconv
Definition: cs_equation_param.h:453
CS_NVD_MUSCL
Definition: cs_convection_diffusion.h:65
cs_numbering_t::group_index
cs_lnum_t * group_index
Definition: cs_numbering.h:98
cs_anisotropic_right_diffusion_vector
void cs_anisotropic_right_diffusion_vector(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int inc, cs_real_3_t *restrict pvar, const cs_real_3_t *restrict pvara, const cs_real_3_t coefav[], const cs_real_33_t coefbv[], const cs_real_3_t cofafv[], const cs_real_33_t cofbfv[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t *restrict viscel, const cs_real_2_t weighf[], const cs_real_t weighb[], cs_real_3_t *restrict rhs)
Add explicit part of the terms of diffusion by a right-multiplying symmetric tensorial diffusivity fo...
Definition: cs_convection_diffusion.c:9471
dt
Definition: cs_field_pointer.h:65
cs_glob_mesh_quantities
cs_mesh_quantities_t * cs_glob_mesh_quantities
cs_mesh_quantities_t
Definition: cs_mesh_quantities.h:89
cs_equation_param_t::blencv
double blencv
Definition: cs_equation_param.h:472
CS_MIN
#define CS_MIN(a, b)
Definition: cs_defs.h:458
cs_halo_type_t
cs_halo_type_t
Definition: cs_halo.h:50
cs_halo.h
cs_face_anisotropic_diffusion_potential
void cs_face_anisotropic_diffusion_potential(const int f_id, const cs_mesh_t *m, cs_mesh_quantities_t *fvq, int init, int inc, int imrgra, int iccocg, int nswrgp, int imligp, int ircflp, int iphydp, int iwgrp, int iwarnp, double epsrgp, double climgp, cs_real_3_t *restrict frcxt, cs_real_t *restrict pvar, const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t *restrict viscel, const cs_real_2_t weighf[], const cs_real_t weighb[], cs_real_t *restrict i_massflux, cs_real_t *restrict b_massflux)
Add the explicit part of the pressure gradient term to the mass flux in case of anisotropic diffusion...
Definition: cs_convection_diffusion.c:11152
CS_NVD_SMART
Definition: cs_convection_diffusion.h:62
cs_lnum_2_t
cs_lnum_t cs_lnum_2_t[2]
vector of 2 local mesh-entity ids
Definition: cs_defs.h:313
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:320
cs_numbering_t::n_groups
int n_groups
Definition: cs_numbering.h:90
cs_equation_param_t::isstpc
int isstpc
Definition: cs_equation_param.h:463
cs_anisotropic_diffusion_potential
void cs_anisotropic_diffusion_potential(const int f_id, const cs_mesh_t *m, cs_mesh_quantities_t *fvq, int init, int inc, int imrgra, int iccocg, int nswrgp, int imligp, int ircflp, int iphydp, int iwgrp, int iwarnp, double epsrgp, double climgp, cs_real_3_t *restrict frcxt, cs_real_t *restrict pvar, const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t *restrict viscel, const cs_real_2_t weighf[], const cs_real_t weighb[], cs_real_t *restrict diverg)
Add the explicit part of the divergence of the mass flux due to the pressure gradient (routine analog...
Definition: cs_convection_diffusion.c:12022
CS_NVD_N_TYPES
Definition: cs_convection_diffusion.h:74
cs_slope_test_gradient
void cs_slope_test_gradient(int f_id, int inc, cs_halo_type_t halo_type, const cs_real_3_t *grad, cs_real_3_t *grdpa, const cs_real_t *pvar, const cs_real_t *coefap, const cs_real_t *coefbp, const cs_real_t *i_massflux)
Compute the upwind gradient used in the slope tests.
Definition: cs_convection_diffusion.c:894
cs_convection_diffusion_tensor
void cs_convection_diffusion_tensor(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int icvflb, int inc, int imasac, cs_real_6_t *restrict pvar, const cs_real_6_t *restrict pvara, const cs_real_6_t coefa[], const cs_real_66_t coefb[], const cs_real_6_t cofaf[], const cs_real_66_t cofbf[], 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 vector field ...
Definition: cs_convection_diffusion.c:5864
coincl::fmin
double precision, save fmin
Definition: coincl.f90:133
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:301
cs_mesh_t::n_cells_with_ghosts
cs_lnum_t n_cells_with_ghosts
Definition: cs_mesh.h:150
CS_NVD_WASEB
Definition: cs_convection_diffusion.h:70
CS_HALO_STANDARD
Definition: cs_halo.h:52
cs_math.h
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:496
CS_NVD_CLAM
Definition: cs_convection_diffusion.h:67
CS_NVD_VOF_CICSAM
Definition: cs_convection_diffusion.h:72
cs_mesh_t::b_face_numbering
cs_numbering_t * b_face_numbering
Definition: cs_mesh.h:162
cs_mesh_t::i_face_numbering
cs_numbering_t * i_face_numbering
Definition: cs_mesh.h:161
CS_PROCF
#define CS_PROCF(x, y)
Definition: cs_defs.h:509
CS_MAX
#define CS_MAX(a, b)
Definition: cs_defs.h:459
CS_COUPLED_FD
Definition: cs_parameters.h:113
cs_numbering_t::n_threads
int n_threads
Definition: cs_numbering.h:89
cs_diffusion_potential
void cs_diffusion_potential(const int f_id, const cs_mesh_t *m, cs_mesh_quantities_t *fvq, int init, int inc, int imrgra, int iccocg, int nswrgp, int imligp, int iphydp, int iwgrp, int iwarnp, double epsrgp, double climgp, cs_real_3_t *restrict frcxt, cs_real_t *restrict pvar, const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_t visel[], cs_real_t *restrict diverg)
Update the cell mass flux divergence with the face pressure (or pressure increment, or pressure double increment) gradient.
Definition: cs_convection_diffusion.c:11639
numvar::kimasf
integer, save kimasf
interior and boundary convective mass flux key ids of the variables
Definition: numvar.f90:187
cs_real_33_t
cs_real_t cs_real_33_t[3][3]
3x3 matrix of floating-point values
Definition: cs_defs.h:327
CS_NVD_VOF_HRIC
Definition: cs_convection_diffusion.h:71
CS_NVD_GAMMA
Definition: cs_convection_diffusion.h:61
CS_NVD_SUPERBEE
Definition: cs_convection_diffusion.h:64
CS_HALO_ROTATION_IGNORE
Definition: cs_halo.h:64
cs_convection_diffusion_scalar
void cs_convection_diffusion_scalar(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int icvflb, int inc, int iccocg, int imasac, cs_real_t *restrict pvar, const cs_real_t *restrict pvara, const int icvfli[], const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_t *restrict rhs)
Add the explicit part of the convection/diffusion terms of a standard transport equation of a scalar ...
Definition: cs_convection_diffusion.c:1622
cs_real_63_t
cs_real_t cs_real_63_t[6][3]
Definition: cs_defs.h:335
cs_mesh_quantities_t::cell_vol
cs_real_t * cell_vol
Definition: cs_mesh_quantities.h:92
numvar::kbmasf
integer, save kbmasf
Definition: numvar.f90:187
cs_mesh_t::i_face_cells
cs_lnum_2_t * i_face_cells
Definition: cs_mesh.h:110
cs_halo_sync_component
void cs_halo_sync_component(const cs_halo_t *halo, cs_halo_type_t sync_mode, cs_halo_rotation_t rotation_op, cs_real_t var[])
Definition: cs_halo.c:1754
itrmav
void itrmav(const int *const f_id, const int *const init, const int *const inc, const int *const imrgra, const int *const iccocg, const int *const nswrgp, const int *const imligp, const int *const ircflp, const int *const iphydp, const int *const iwgrp, const int *const iwarnp, const cs_real_t *const epsrgp, const cs_real_t *const climgp, const cs_real_t *const extrap, cs_real_3_t frcxt[], cs_real_t pvar[], const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_visc[], const cs_real_t b_visc[], cs_real_6_t viscel[], const cs_real_2_t weighf[], const cs_real_t weighb[], cs_real_t i_massflux[], cs_real_t b_massflux[])
Definition: cs_convection_diffusion.c:679
cs_beta_limiter_building
void cs_beta_limiter_building(int f_id, int inc, const cs_real_t rovsdt[])
Compute the beta blending coefficient of the beta limiter (ensuring preservation of a given min/max p...
Definition: cs_convection_diffusion.c:1471
cs_convection_diffusion_thermal
void cs_convection_diffusion_thermal(int idtvar, int f_id, const cs_var_cal_opt_t var_cal_opt, int inc, int iccocg, int imasac, cs_real_t *restrict pvar, const cs_real_t *restrict pvara, const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t cofafp[], const cs_real_t cofbfp[], const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t i_visc[], const cs_real_t b_visc[], 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.c:6830
cpincl::b1
double precision, dimension(ncharm), save b1
Definition: cpincl.f90:233
cs_upwind_gradient
void cs_upwind_gradient(const int f_id, const int inc, const cs_halo_type_t halo_type, const cs_real_t coefap[], const cs_real_t coefbp[], const cs_real_t i_massflux[], const cs_real_t b_massflux[], const cs_real_t *restrict pvar, cs_real_3_t *restrict grdpa)
Compute the upwind gradient in order to cope with SOLU schemes observed in the litterature.
Definition: cs_convection_diffusion.c:1038
CS_NVD_OSHER
Definition: cs_convection_diffusion.h:69
optcal::imrgra
integer(c_int), pointer, save imrgra
type of gradient reconstruction
Definition: optcal.f90:251
cs_mesh_t::halo
cs_halo_t * halo
Definition: cs_mesh.h:155
cs_mesh_t::b_face_cells
cs_lnum_t * b_face_cells
Definition: cs_mesh.h:111
CS_NVD_VOF_STACS
Definition: cs_convection_diffusion.h:73