1 /*============================================================================
2 * Building of the right hand side for a transport of a field.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------*/
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <errno.h>
37 #include <stdio.h>
38 #include <stdarg.h>
39 #include <string.h>
40 #include <math.h>
41 #include <float.h>
42
43 #if defined(HAVE_MPI)
44 #include <mpi.h>
45 #endif
46
47 /*----------------------------------------------------------------------------
48 * Local headers
49 *----------------------------------------------------------------------------*/
50
51 #include "bft_error.h"
52 #include "bft_mem.h"
53 #include "bft_printf.h"
54
55 #include "cs_blas.h"
56 #include "cs_halo.h"
57 #include "cs_halo_perio.h"
58 #include "cs_log.h"
59 #include "cs_math.h"
60 #include "cs_mesh.h"
61 #include "cs_convection_diffusion.h"
62 #include "cs_field.h"
63 #include "cs_field_operator.h"
64 #include "cs_field_pointer.h"
65 #include "cs_gradient.h"
66 #include "cs_ext_neighborhood.h"
67 #include "cs_mesh_quantities.h"
68 #include "cs_parall.h"
69 #include "cs_parameters.h"
70 #include "cs_prototypes.h"
71 #include "cs_timer.h"
72 #include "cs_velocity_pressure.h"
73
74 /*----------------------------------------------------------------------------
75 * Header for the current file
76 *----------------------------------------------------------------------------*/
77
78 #include "cs_balance.h"
79
80 /*----------------------------------------------------------------------------*/
81
82 BEGIN_C_DECLS
83
84 /*=============================================================================
85 * Local macro definitions
86 *============================================================================*/
87
88 /*============================================================================
89 * Public function definitions
90 *============================================================================*/
91
92 /*----------------------------------------------------------------------------*/
93 /*! \file cs_balance.c
94 *
95 * \brief Wrapper to the function which adds the explicit part of the
96 * convection/diffusion terms of a transport equation of a field.
97 */
98 /*----------------------------------------------------------------------------*/
99
100 /*----------------------------------------------------------------------------*/
101 /*!
102 * \brief Wrapper to the function which adds the explicit part of the
103 * convection/diffusion terms of a transport equation of
104 * a scalar field \f$ \varia \f$.
105 *
106 * More precisely, the right hand side \f$ Rhs \f$ is updated as
107 * follows:
108 * \f[
109 * Rhs = Rhs - \sum_{\fij \in \Facei{\celli}} \left(
110 * \dot{m}_\ij \left( \varia_\fij - \varia_\celli \right)
111 * - \mu_\fij \gradv_\fij \varia \cdot \vect{S}_\ij \right)
112 * \f]
113 *
114 * Warning:
115 * - \f$ Rhs \f$ has already been initialized
116 * before calling cs_balance_scalar!
117 * - mind the minus sign
118 *
119 * Options for the convective scheme:
120 * - blencp = 0: upwind scheme for the advection
121 * - blencp = 1: no upwind scheme except in the slope test
122 * - ischcp = 0: SOLU
123 * - ischcp = 1: centered
124 * - ischcp = 2: SOLU with upwind gradient reconstruction
125 * - ischcp = 3: blending SOLU centered
126 * - ischcp = 4: NVD-TVD
127 * - imucpp = 0: do not multiply the convective part by \f$ C_p \f$
128 * - imucpp = 1: multiply the convective part by \f$ C_p \f$
129 *
130 * \param[in] idtvar indicator of the temporal scheme
131 * \param[in] f_id field id (or -1)
132 * \param[in] imucpp indicator
133 * - 0 do not multiply the convectiv term by Cp
134 * - 1 do multiply the convectiv term by Cp
135 * \param[in] imasac take mass accumulation into account?
136 * \param[in] inc indicator
137 * - 0 when solving an increment
138 * - 1 otherwise
139 * \param[in] iccocg indicator
140 * - 1 re-compute cocg matrix
141 * (for iterative gradients)
142 * - 0 otherwise
143 * \param[in] var_cal_opt pointer to a cs_var_cal_opt_t structure which
144 * contains variable calculation options
145 * \param[in] pvar solved variable (current time step)
146 * may be NULL if pvara != NULL
147 * \param[in] pvara solved variable (previous time step)
148 * may be NULL if pvar != NULL
149 * \param[in] coefap boundary condition array for the variable
150 * (explicit part)
151 * \param[in] coefbp boundary condition array for the variable
152 * (implicit part)
153 * \param[in] cofafp boundary condition array for the diffusion
154 * of the variable (explicit part)
155 * \param[in] cofbfp boundary condition array for the diffusion
156 * of the variable (implicit part)
157 * \param[in] i_massflux mass flux at interior faces
158 * \param[in] b_massflux mass flux at boundary faces
159 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
160 * at interior faces for the r.h.s.
161 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
162 * at boundary faces for the r.h.s.
163 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
164 * \param[in] xcpp array of specific heat (Cp)
165 * \param[in] weighf internal face weight between cells i j in case
166 * of tensor diffusion
167 * \param[in] weighb boundary face weight for cells i in case
168 * of tensor diffusion
169 * \param[in] icvflb global indicator of boundary convection flux
170 * - 0 upwind scheme at all boundary faces
171 * - 1 imposed flux at some boundary faces
172 * \param[in] icvfli boundary face indicator array of convection flux
173 * - 0 upwind scheme
174 * - 1 imposed flux
175 * \param[in,out] smbrp right hand side \f$ \vect{Rhs} \f$
176 */
177 /*----------------------------------------------------------------------------*/
178
179 void
cs_balance_scalar(int idtvar,int f_id,int imucpp,int imasac,int inc,int iccocg,cs_var_cal_opt_t * var_cal_opt,cs_real_t pvar[],const cs_real_t 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[],cs_real_6_t viscel[],const cs_real_t xcpp[],const cs_real_2_t weighf[],const cs_real_t weighb[],int icvflb,const int icvfli[],cs_real_t smbrp[])180 cs_balance_scalar(int idtvar,
181 int f_id,
182 int imucpp,
183 int imasac,
184 int inc,
185 int iccocg,
186 cs_var_cal_opt_t *var_cal_opt,
187 cs_real_t pvar[],
188 const cs_real_t pvara[],
189 const cs_real_t coefap[],
190 const cs_real_t coefbp[],
191 const cs_real_t cofafp[],
192 const cs_real_t cofbfp[],
193 const cs_real_t i_massflux[],
194 const cs_real_t b_massflux[],
195 const cs_real_t i_visc[],
196 const cs_real_t b_visc[],
197 cs_real_6_t viscel[],
198 const cs_real_t xcpp[],
199 const cs_real_2_t weighf[],
200 const cs_real_t weighb[],
201 int icvflb,
202 const int icvfli[],
203 cs_real_t smbrp[])
204 {
205 /* Local variables */
206 int iconvp = var_cal_opt->iconv;
207 int idiffp = var_cal_opt->idiff;
208 int idftnp = var_cal_opt->idften;
209 cs_var_cal_opt_t var_cal_opt_loc;
210
211 if (f_id < 0) {
212 var_cal_opt_loc = cs_parameters_var_cal_opt_default();
213 var_cal_opt_loc.verbosity = var_cal_opt->verbosity;
214 var_cal_opt_loc.iconv = var_cal_opt->iconv;
215 var_cal_opt_loc.istat = -1; /* unused in balance */
216 var_cal_opt_loc.idiff = var_cal_opt->idiff;
217 var_cal_opt_loc.idifft = -1; /* unused in balance */
218 var_cal_opt_loc.idften = var_cal_opt->idften;
219 var_cal_opt_loc.iswdyn = -1; /* unused in balance */
220 var_cal_opt_loc.ischcv = var_cal_opt->ischcv;
221 var_cal_opt_loc.isstpc = var_cal_opt->isstpc;
222 var_cal_opt_loc.nswrgr = var_cal_opt->nswrgr;
223 var_cal_opt_loc.nswrsm = -1; /* unused in balance */
224 var_cal_opt_loc.imrgra = var_cal_opt->imrgra;
225 var_cal_opt_loc.imligr = var_cal_opt->imligr;
226 var_cal_opt_loc.ircflu = var_cal_opt->ircflu;
227 var_cal_opt_loc.iwgrec = 0; /* require field id */
228 var_cal_opt_loc.icoupl = -1; /* require field id */
229 var_cal_opt_loc.thetav = var_cal_opt->thetav;
230 var_cal_opt_loc.blencv = var_cal_opt->blencv;
231 var_cal_opt_loc.blend_st = var_cal_opt->blend_st;
232 var_cal_opt_loc.epsilo = -1.; /* unused in balance */
233 var_cal_opt_loc.epsrsm = -1.; /* unused in balance */
234 var_cal_opt_loc.epsrgr = var_cal_opt->epsrgr;
235 var_cal_opt_loc.climgr = var_cal_opt->climgr;
236 var_cal_opt_loc.relaxv = var_cal_opt->relaxv;
237 } else {
238 cs_field_t *f = cs_field_by_id(f_id);
239 int k_id = cs_field_key_id("var_cal_opt");
240 cs_field_get_key_struct(f, k_id, &var_cal_opt_loc);
241 var_cal_opt_loc.thetav = var_cal_opt->thetav;
242 }
243
244 /* Scalar diffusivity */
245 if (idftnp & CS_ISOTROPIC_DIFFUSION) {
246 if (imucpp == 0) {
247 cs_convection_diffusion_scalar(idtvar,
248 f_id,
249 var_cal_opt_loc,
250 icvflb,
251 inc,
252 iccocg,
253 imasac,
254 pvar,
255 pvara,
256 icvfli,
257 coefap,
258 coefbp,
259 cofafp,
260 cofbfp,
261 i_massflux,
262 b_massflux,
263 i_visc,
264 b_visc,
265 smbrp);
266 } else {
267 /* The convective part is multiplied by Cp for the temperature */
268 cs_convection_diffusion_thermal(idtvar,
269 f_id,
270 var_cal_opt_loc,
271 inc,
272 iccocg,
273 imasac,
274 pvar,
275 pvara,
276 coefap,
277 coefbp,
278 cofafp,
279 cofbfp,
280 i_massflux,
281 b_massflux,
282 i_visc,
283 b_visc,
284 xcpp,
285 smbrp);
286 }
287 }
288 /* Symmetric tensor diffusivity */
289 else if (idftnp & CS_ANISOTROPIC_DIFFUSION) {
290 var_cal_opt_loc.idiff = 0;
291 /* Convective part */
292 if (imucpp == 0 && iconvp == 1) {
293 cs_convection_diffusion_scalar(idtvar,
294 f_id,
295 var_cal_opt_loc,
296 icvflb,
297 inc,
298 iccocg,
299 imasac,
300 pvar,
301 pvara,
302 icvfli,
303 coefap,
304 coefbp,
305 cofafp,
306 cofbfp,
307 i_massflux,
308 b_massflux,
309 i_visc,
310 b_visc,
311 smbrp);
312 }
313 /* The convective part is multiplied by Cp for the temperature */
314 else if (imucpp == 1 && iconvp == 1) {
315 cs_convection_diffusion_thermal(idtvar,
316 f_id,
317 var_cal_opt_loc,
318 inc,
319 iccocg,
320 imasac,
321 pvar,
322 pvara,
323 coefap,
324 coefbp,
325 cofafp,
326 cofbfp,
327 i_massflux,
328 b_massflux,
329 i_visc,
330 b_visc,
331 xcpp,
332 smbrp);
333 }
334
335 /* Diffusive part */
336 if (idiffp == 1) {
337 cs_anisotropic_diffusion_scalar(idtvar,
338 f_id,
339 var_cal_opt_loc,
340 inc,
341 iccocg,
342 pvar,
343 pvara,
344 coefap,
345 coefbp,
346 cofafp,
347 cofbfp,
348 i_visc,
349 b_visc,
350 viscel,
351 weighf,
352 weighb,
353 smbrp);
354 }
355 }
356 }
357
358 /*----------------------------------------------------------------------------*/
359 /*!
360 * \brief Wrapper to the function which adds the explicit part of the
361 * convection/diffusion
362 * terms of a transport equation of a vector field \f$ \vect{\varia} \f$.
363 *
364 * More precisely, the right hand side \f$ \vect{Rhs} \f$ is updated as
365 * follows:
366 * \f[
367 * \vect{Rhs} = \vect{Rhs} - \sum_{\fij \in \Facei{\celli}} \left(
368 * \dot{m}_\ij \left( \vect{\varia}_\fij - \vect{\varia}_\celli \right)
369 * - \mu_\fij \gradt_\fij \vect{\varia} \cdot \vect{S}_\ij \right)
370 * \f]
371 *
372 * Remark:
373 * if ivisep = 1, then we also take \f$ \mu \transpose{\gradt\vect{\varia}}
374 * + \lambda \trace{\gradt\vect{\varia}} \f$, where \f$ \lambda \f$ is
375 * the secondary viscosity, i.e. usually \f$ -\frac{2}{3} \mu \f$.
376 *
377 * Warning:
378 * - \f$ \vect{Rhs} \f$ has already been initialized
379 * before calling cs_balance_vector!
380 * - mind the sign minus
381 *
382 * Options for the convective scheme:
383 * - blencp = 0: upwind scheme for the advection
384 * - blencp = 1: no upwind scheme except in the slope test
385 * - ischcp = 0: SOLU
386 * - ischcp = 1: centered
387 * - ischcp = 2: SOLU with upwind gradient reconstruction
388 * - ischcp = 3: blending SOLU centered
389 * - ischcp = 4: NVD-TVD
390 *
391 * \param[in] idtvar indicator of the temporal scheme
392 * \param[in] f_id field id (or -1)
393 * \param[in] imasac take mass accumulation into account?
394 * \param[in] inc indicator
395 * - 0 when solving an increment
396 * - 1 otherwise
397 * \param[in] ivisep indicator to take \f$ \divv
398 * \left(\mu \gradt \transpose{\vect{a}} \right)
399 * -2/3 \grad\left( \mu \dive \vect{a} \right)\f$
400 * - 1 take into account,
401 * - 0 otherwise
402 * \param[in] var_cal_opt pointer to a cs_var_cal_opt_t structure which
403 * contains variable calculation options
404 * \param[in] pvar solved velocity (current time step)
405 * \param[in] pvara solved velocity (previous time step)
406 * \param[in] coefav boundary condition array for the variable
407 * (explicit part)
408 * \param[in] coefbv boundary condition array for the variable
409 * (implicit part)
410 * \param[in] cofafv boundary condition array for the diffusion
411 * of the variable (explicit part)
412 * \param[in] cofbfv boundary condition array for the diffusion
413 * of the variable (implicit part)
414 * \param[in] i_massflux mass flux at interior faces
415 * \param[in] b_massflux mass flux at boundary faces
416 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
417 * at interior faces for the r.h.s.
418 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
419 * at boundary faces for the r.h.s.
420 * \param[in] secvif secondary viscosity at interior faces
421 * \param[in] secvib secondary viscosity at boundary faces
422 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
423 * \param[in] weighf internal face weight between cells i j in case
424 * of tensor diffusion
425 * \param[in] weighb boundary face weight for cells i in case
426 * of tensor diffusion
427 * \param[in] icvflb global indicator of boundary convection flux
428 * - 0 upwind scheme at all boundary faces
429 * - 1 imposed flux at some boundary faces
430 * \param[in] icvfli boundary face indicator array of convection flux
431 * - 0 upwind scheme
432 * - 1 imposed flux
433 * \param[in,out] smbr right hand side \f$ \vect{Rhs} \f$
434 */
435 /*----------------------------------------------------------------------------*/
436
437 void
cs_balance_vector(int idtvar,int f_id,int imasac,int inc,int ivisep,cs_var_cal_opt_t * var_cal_opt,cs_real_3_t pvar[],const cs_real_3_t 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_massflux[],const cs_real_t b_massflux[],const cs_real_t i_visc[],const cs_real_t b_visc[],const cs_real_t secvif[],const cs_real_t secvib[],cs_real_6_t viscel[],const cs_real_2_t weighf[],const cs_real_t weighb[],int icvflb,const int icvfli[],cs_real_3_t smbr[])438 cs_balance_vector(int idtvar,
439 int f_id,
440 int imasac,
441 int inc,
442 int ivisep,
443 cs_var_cal_opt_t *var_cal_opt,
444 cs_real_3_t pvar[],
445 const cs_real_3_t pvara[],
446 const cs_real_3_t coefav[],
447 const cs_real_33_t coefbv[],
448 const cs_real_3_t cofafv[],
449 const cs_real_33_t cofbfv[],
450 const cs_real_t i_massflux[],
451 const cs_real_t b_massflux[],
452 const cs_real_t i_visc[],
453 const cs_real_t b_visc[],
454 const cs_real_t secvif[],
455 const cs_real_t secvib[],
456 cs_real_6_t viscel[],
457 const cs_real_2_t weighf[],
458 const cs_real_t weighb[],
459 int icvflb,
460 const int icvfli[],
461 cs_real_3_t smbr[])
462 {
463 /* Local variables */
464 int iconvp = var_cal_opt->iconv;
465 int idiffp = var_cal_opt->idiff;
466 int idftnp = var_cal_opt->idften;
467 cs_var_cal_opt_t var_cal_opt_loc;
468
469 if (f_id < 0) {
470 var_cal_opt_loc = cs_parameters_var_cal_opt_default();
471 var_cal_opt_loc.verbosity = var_cal_opt->verbosity;
472 var_cal_opt_loc.iconv = var_cal_opt->iconv;
473 var_cal_opt_loc.istat = -1; /* unused in balance */
474 var_cal_opt_loc.idiff = var_cal_opt->idiff;
475 var_cal_opt_loc.idifft = -1; /* unused in balance */
476 var_cal_opt_loc.idften = var_cal_opt->idften;
477 var_cal_opt_loc.iswdyn = -1; /* unused in balance */
478 var_cal_opt_loc.ischcv = var_cal_opt->ischcv;
479 var_cal_opt_loc.isstpc = var_cal_opt->isstpc;
480 var_cal_opt_loc.nswrgr = var_cal_opt->nswrgr;
481 var_cal_opt_loc.nswrsm = -1; /* unused in balance */
482 var_cal_opt_loc.imrgra = var_cal_opt->imrgra;
483 var_cal_opt_loc.imligr = var_cal_opt->imligr;
484 var_cal_opt_loc.ircflu = var_cal_opt->ircflu;
485 var_cal_opt_loc.iwgrec = 0; /* require field id */
486 var_cal_opt_loc.icoupl = -1; /* require field id */
487 var_cal_opt_loc.thetav = var_cal_opt->thetav;
488 var_cal_opt_loc.blencv = var_cal_opt->blencv;
489 var_cal_opt_loc.blend_st = var_cal_opt->blend_st;
490 var_cal_opt_loc.epsilo = -1.; /* unused in balance */
491 var_cal_opt_loc.epsrsm = -1.; /* unused in balance */
492 var_cal_opt_loc.epsrgr = var_cal_opt->epsrgr;
493 var_cal_opt_loc.climgr = var_cal_opt->climgr;
494 var_cal_opt_loc.relaxv = var_cal_opt->relaxv;
495 } else {
496 cs_field_t *f = cs_field_by_id(f_id);
497 int k_id = cs_field_key_id("var_cal_opt");
498 cs_field_get_key_struct(f, k_id, &var_cal_opt_loc);
499 var_cal_opt_loc.thetav = var_cal_opt->thetav;
500 }
501
502 /* Scalar diffusivity */
503 if (idftnp & CS_ISOTROPIC_DIFFUSION) {
504 cs_convection_diffusion_vector(idtvar,
505 f_id,
506 var_cal_opt_loc,
507 icvflb,
508 inc,
509 ivisep,
510 imasac,
511 pvar,
512 pvara,
513 icvfli,
514 coefav,
515 coefbv,
516 cofafv,
517 cofbfv,
518 i_massflux,
519 b_massflux,
520 i_visc,
521 b_visc,
522 secvif,
523 secvib,
524 smbr);
525 }
526 /* Symmetric tensor diffusivity */
527 else if (idftnp & CS_ANISOTROPIC_DIFFUSION) {
528 /* ! Nor diffusive part neither secondary viscosity or transpose of gradient */
529 var_cal_opt_loc.idiff = 0;
530 /* Convective part */
531 if (iconvp == 1) {
532 cs_convection_diffusion_vector(idtvar,
533 f_id,
534 var_cal_opt_loc,
535 icvflb,
536 inc,
537 ivisep,
538 imasac,
539 pvar,
540 pvara,
541 icvfli,
542 coefav,
543 coefbv,
544 cofafv,
545 cofbfv,
546 i_massflux,
547 b_massflux,
548 i_visc,
549 b_visc,
550 secvif,
551 secvib,
552 smbr);
553
554 }
555
556 /* Diffusive part (with a 3x3 symmetric diffusivity) */
557
558 /* either Daly-Harlow type i.e. Nabla(v).K */
559 if (idiffp == 1 && idftnp & CS_ANISOTROPIC_RIGHT_DIFFUSION) {
560 /* ! Neither diffusive part neither secondary viscosity
561 nor transpose of gradient */
562 cs_anisotropic_right_diffusion_vector(idtvar,
563 f_id,
564 var_cal_opt_loc,
565 inc,
566 pvar,
567 pvara,
568 coefav,
569 coefbv,
570 cofafv,
571 cofbfv,
572 i_visc,
573 b_visc,
574 viscel,
575 weighf,
576 weighb,
577 smbr);
578 }
579
580 /* or K.Nabla(v) ) */
581 else if (idiffp == 1 && idftnp & CS_ANISOTROPIC_LEFT_DIFFUSION) {
582 cs_anisotropic_left_diffusion_vector(idtvar,
583 f_id,
584 var_cal_opt_loc,
585 inc,
586 ivisep,
587 pvar,
588 pvara,
589 coefav,
590 coefbv,
591 cofafv,
592 cofbfv,
593 (const cs_real_33_t *)i_visc,
594 b_visc,
595 secvif,
596 smbr);
597 }
598 }
599 }
600
601 /*----------------------------------------------------------------------------*/
602 /*!
603 * \brief Wrapper to the function which adds the explicit part of the
604 * convection/diffusion
605 * terms of a transport equation of a tensor field \f$ \tens{\varia} \f$.
606 *
607 * More precisely, the right hand side \f$ \vect{Rhs} \f$ is updated as
608 * follows:
609 * \f[
610 * \tens{Rhs} = \tens{Rhs} - \sum_{\fij \in \Facei{\celli}} \left(
611 * \dot{m}_\ij \left( \tens{\varia}_\fij - \tens{\varia}_\celli \right)
612 * - \mu_\fij \gradt_\fij \tens{\varia} \cdot \tens{S}_\ij \right)
613 * \f]
614 *
615 * Warning:
616 * - \f$ \tens{Rhs} \f$ has already been initialized before calling bilscts!
617 * - mind the sign minus
618 *
619 * Options for the convective scheme:
620 * - blencp = 0: upwind scheme for the advection
621 * - blencp = 1: no upwind scheme except in the slope test
622 * - ischcp = 0: SOLU
623 * - ischcp = 1: centered
624 * - ischcp = 2: SOLU with upwind gradient reconstruction
625 * - ischcp = 3: blending SOLU centered
626 * - ischcp = 4: NVD-TVD
627 *
628 * \param[in] idtvar indicator of the temporal scheme
629 * \param[in] f_id field id (or -1)
630 * \param[in] imasac take mass accumulation into account?
631 * \param[in] inc indicator
632 * \param[in] var_cal_opt pointer to a cs_var_cal_opt_t structure which
633 * contains variable calculation options
634 * \param[in] pvar solved velocity (current time step)
635 * \param[in] pvara solved velocity (previous time step)
636 * \param[in] coefa boundary condition array for the variable
637 * (Explicit part)
638 * \param[in] coefb boundary condition array for the variable
639 * (Impplicit part)
640 * \param[in] cofaf boundary condition array for the diffusion
641 * of the variable (Explicit part)
642 * \param[in] cofbf boundary condition array for the diffusion
643 * of the variable (Implicit part)
644 * \param[in] i_massflux mass flux at interior faces
645 * \param[in] b_massflux mass flux at boundary faces
646 * \param[in] i_visc \f$ \mu_\fij \dfrac{S_\fij}{\ipf \jpf} \f$
647 * at interior faces for the r.h.s.
648 * \param[in] b_visc \f$ \mu_\fib \dfrac{S_\fib}{\ipf \centf} \f$
649 * at boundary faces for the r.h.s.
650 * \param[in] viscel symmetric cell tensor \f$ \tens{\mu}_\celli \f$
651 * \param[in] weighf internal face weight between cells i j in case
652 * of tensor diffusion
653 * \param[in] weighb boundary face weight for cells i in case
654 * of tensor diffusion
655 * \param[in] icvflb global indicator of boundary convection flux
656 * - 0 upwind scheme at all boundary faces
657 * - 1 imposed flux at some boundary faces
658 * \param[in] icvfli boundary face indicator array of convection flux
659 * - 0 upwind scheme
660 * - 1 imposed flux
661 * \param[in,out] smbrp right hand side \f$ \vect{Rhs} \f$
662 */
663 /*----------------------------------------------------------------------------*/
664
665 void
cs_balance_tensor(int idtvar,int f_id,int imasac,int inc,cs_var_cal_opt_t * var_cal_opt,cs_real_6_t pvar[],const cs_real_6_t 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 viscel[],const cs_real_2_t weighf[],const cs_real_t weighb[],int icvflb,const int icvfli[],cs_real_6_t smbrp[])666 cs_balance_tensor(int idtvar,
667 int f_id,
668 int imasac,
669 int inc,
670 cs_var_cal_opt_t *var_cal_opt,
671 cs_real_6_t pvar[],
672 const cs_real_6_t pvara[],
673 const cs_real_6_t coefa[],
674 const cs_real_66_t coefb[],
675 const cs_real_6_t cofaf[],
676 const cs_real_66_t cofbf[],
677 const cs_real_t i_massflux[],
678 const cs_real_t b_massflux[],
679 const cs_real_t i_visc[],
680 const cs_real_t b_visc[],
681 cs_real_6_t viscel[],
682 const cs_real_2_t weighf[],
683 const cs_real_t weighb[],
684 int icvflb,
685 const int icvfli[],
686 cs_real_6_t smbrp[])
687 {
688 CS_UNUSED(icvfli);
689
690 /* Local variables */
691 int iconvp = var_cal_opt->iconv;
692 int idiffp = var_cal_opt->idiff;
693 int idftnp = var_cal_opt->idften;
694 cs_var_cal_opt_t var_cal_opt_loc;
695
696 if (f_id < 0) {
697 var_cal_opt_loc = cs_parameters_var_cal_opt_default();
698 var_cal_opt_loc.verbosity = var_cal_opt->verbosity;
699 var_cal_opt_loc.iconv = var_cal_opt->iconv;
700 var_cal_opt_loc.istat = -1; /* unused in balance */
701 var_cal_opt_loc.idiff = var_cal_opt->idiff;
702 var_cal_opt_loc.idifft = -1; /* unused in balance */
703 var_cal_opt_loc.idften = var_cal_opt->idften;
704 var_cal_opt_loc.iswdyn = -1; /* unused in balance */
705 var_cal_opt_loc.ischcv = var_cal_opt->ischcv;
706 var_cal_opt_loc.isstpc = var_cal_opt->isstpc;
707 var_cal_opt_loc.nswrgr = var_cal_opt->nswrgr;
708 var_cal_opt_loc.nswrsm = -1; /* unused in balance */
709 var_cal_opt_loc.imrgra = var_cal_opt->imrgra;
710 var_cal_opt_loc.imligr = var_cal_opt->imligr;
711 var_cal_opt_loc.ircflu = var_cal_opt->ircflu;
712 var_cal_opt_loc.iwgrec = 0; /* require field id */
713 var_cal_opt_loc.icoupl = -1; /* require field id */
714 var_cal_opt_loc.thetav = var_cal_opt->thetav;
715 var_cal_opt_loc.blencv = var_cal_opt->blencv;
716 var_cal_opt_loc.blend_st = var_cal_opt->blend_st;
717 var_cal_opt_loc.epsilo = -1.; /* unused in balance */
718 var_cal_opt_loc.epsrsm = -1.; /* unused in balance */
719 var_cal_opt_loc.epsrgr = var_cal_opt->epsrgr;
720 var_cal_opt_loc.climgr = var_cal_opt->climgr;
721 var_cal_opt_loc.relaxv = var_cal_opt->relaxv;
722 } else {
723 cs_field_t *f = cs_field_by_id(f_id);
724 int k_id = cs_field_key_id("var_cal_opt");
725 cs_field_get_key_struct(f, k_id, &var_cal_opt_loc);
726 var_cal_opt_loc.thetav = var_cal_opt->thetav;
727 }
728
729 /* Scalar diffusivity */
730 if (idftnp & CS_ISOTROPIC_DIFFUSION) {
731 cs_convection_diffusion_tensor(idtvar,
732 f_id,
733 var_cal_opt_loc,
734 icvflb,
735 inc,
736 imasac,
737 pvar,
738 pvara,
739 coefa,
740 coefb,
741 cofaf,
742 cofbf,
743 i_massflux,
744 b_massflux,
745 i_visc,
746 b_visc,
747 smbrp);
748 }
749 /* Symmetric tensor diffusivity */
750 else if (idftnp & CS_ANISOTROPIC_RIGHT_DIFFUSION) {
751 /* No diffusive part */
752 var_cal_opt_loc.idiff = 0;
753 /* Convective part */
754 if (iconvp == 1) {
755 cs_convection_diffusion_tensor(idtvar,
756 f_id,
757 var_cal_opt_loc,
758 icvflb,
759 inc,
760 imasac,
761 pvar,
762 pvara,
763 coefa,
764 coefb,
765 cofaf,
766 cofbf,
767 i_massflux,
768 b_massflux,
769 i_visc,
770 b_visc,
771 smbrp);
772 }
773 /* Diffusive part (with a 6x6 symmetric diffusivity) */
774 if (idiffp == 1) {
775 cs_anisotropic_diffusion_tensor(idtvar,
776 f_id,
777 var_cal_opt_loc,
778 inc,
779 pvar,
780 pvara,
781 coefa,
782 coefb,
783 cofaf,
784 cofbf,
785 i_visc,
786 b_visc,
787 viscel,
788 weighf,
789 weighb,
790 smbrp);
791 }
792 }
793 }
794
795 /*----------------------------------------------------------------------------*/
796
797 END_C_DECLS
798