1 /*============================================================================
2 * Build discrete convection operators for CDO schemes
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 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <stdio.h>
34 #include <stdlib.h>
35 #include <math.h>
36 #include <float.h>
37 #include <limits.h>
38 #include <assert.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include <bft_mem.h>
45
46 #include "cs_cdo_bc.h"
47 #include "cs_flag.h"
48 #include "cs_property.h"
49 #include "cs_math.h"
50 #include "cs_scheme_geometry.h"
51
52 #if defined(DEBUG) && !defined(NDEBUG) /* For debugging purpose */
53 #include "cs_dbg.h"
54 #include "cs_log.h"
55 #endif
56
57 /*----------------------------------------------------------------------------
58 * Header for the current file
59 *----------------------------------------------------------------------------*/
60
61 #include "cs_cdo_advection.h"
62
63 /*----------------------------------------------------------------------------*/
64
65 BEGIN_C_DECLS
66
67 /*=============================================================================
68 * Additional doxygen documentation
69 *============================================================================*/
70
71 /*!
72 \file cs_cdo_advection.c
73
74 \brief Build discrete advection operators for CDO vertex-based schemes
75
76 */
77
78 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
79
80 /*=============================================================================
81 * Local Macro definitions
82 *============================================================================*/
83
84 #define CS_CDO_ADVECTION_DBG 0
85
86 /* Redefined the name of functions from cs_math to get shorter names */
87 #define _dp3 cs_math_3_dot_product
88
89 /*=============================================================================
90 * Local structure definitions
91 *============================================================================*/
92
93 /*----------------------------------------------------------------------------*/
94 /*!
95 * \brief Compute the value of the weighting function related to upwinding
96 *
97 * \param[in] criterion dot product between advection and normal vectors or
98 * estimation of a local Peclet number
99 *
100 * \return the weight value
101 */
102 /*----------------------------------------------------------------------------*/
103
104 typedef double
105 (_upwind_weight_t)(double criterion);
106
107 /*----------------------------------------------------------------------------*/
108 /*!
109 * \brief Update the system by taking into account the boundary conditions
110 *
111 * \param[in] csys pointer to a cellwise view of the system
112 * \param[in] bf local face number id (0..n_bc_faces)
113 * \param[in] flx advective flux across the triangle "tef"
114 * \param[in] v1 first vertex to consider in the cell numbering
115 * \param[in] v2 second vertex to consider in the cell numbering
116 * \param[in, out] rhs rhs of the local system
117 * \param[in, out] diag diagonal of the local system
118 */
119 /*----------------------------------------------------------------------------*/
120
121 typedef void
122 (_update_vb_system_with_bc_t)(const cs_cell_sys_t *csys,
123 const short int bf,
124 const double flx,
125 const short int v1,
126 const short int v2,
127 double *rhs,
128 double *diag);
129
130 /*============================================================================
131 * Local variables
132 *============================================================================*/
133
134 /*============================================================================
135 * Private constant variables
136 *============================================================================*/
137
138 /* Advanced developer parameters (stabilization coefficient) */
139 static double cs_cip_stab_coef = 1e-2;
140
141 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
142
143 /*============================================================================
144 * Private function prototypes
145 *============================================================================*/
146
147 /*----------------------------------------------------------------------------*/
148 /*!
149 * \brief Compute the weight of upwinding for an upwind scheme
150 *
151 * \param[in] criterion dot product between advection and normal vectors or
152 * estimation of a local Peclet number
153 *
154 * \return the weight value
155 */
156 /*----------------------------------------------------------------------------*/
157
158 inline static double
_get_upwind_weight(double criterion)159 _get_upwind_weight(double criterion)
160 {
161 if (criterion > cs_math_zero_threshold)
162 return 1;
163 else if (criterion < -cs_math_zero_threshold)
164 return 0;
165 else
166 return 0.5;
167 }
168
169 /*----------------------------------------------------------------------------*/
170 /*!
171 * \brief Compute the weight of upwinding for an Samarskii scheme
172 *
173 * \param[in] criterion dot product between advection and normal vectors or
174 * estimation of a local Peclet number
175 *
176 * \return the weight value
177 */
178 /*----------------------------------------------------------------------------*/
179
180 inline static double
_get_samarskii_weight(double criterion)181 _get_samarskii_weight(double criterion)
182 {
183 double weight = 1.0;
184 if (criterion < 0)
185 weight = 1./(2 - criterion);
186 else
187 weight = (1 + criterion)/(2 + criterion);
188
189 return weight;
190 }
191
192 /*----------------------------------------------------------------------------*/
193 /*!
194 * \brief Compute the weight of upwinding for an Sharfetter-Gummel scheme
195 *
196 * \param[in] criterion dot product between advection and normal vectors or
197 * estimation of a local Peclet number
198 *
199 * \return the weight value
200 */
201 /*----------------------------------------------------------------------------*/
202
203 inline static double
_get_sg_weight(double criterion)204 _get_sg_weight(double criterion)
205 {
206 double weight = 1.0;
207
208 if (criterion < 0)
209 weight = 0.5*exp(criterion);
210 else
211 weight = 1 - 0.5*exp(-criterion);
212
213 return weight;
214 }
215
216 /*----------------------------------------------------------------------------*/
217 /*!
218 * \brief Retrieve the function related to the kind of upwind weighting
219 *
220 * \param[in] scheme scheme used for discretizing the advection term
221 *
222 * \return a function pointer
223 */
224 /*----------------------------------------------------------------------------*/
225
226 inline static _upwind_weight_t *
_assign_weight_func(const cs_param_advection_scheme_t scheme)227 _assign_weight_func(const cs_param_advection_scheme_t scheme)
228 {
229 switch (scheme) {
230
231 case CS_PARAM_ADVECTION_SCHEME_UPWIND:
232 return _get_upwind_weight;
233
234 case CS_PARAM_ADVECTION_SCHEME_SAMARSKII:
235 return _get_samarskii_weight;
236
237 case CS_PARAM_ADVECTION_SCHEME_SG: /* Sharfetter-Gummel */
238 return _get_sg_weight;
239
240 default:
241 bft_error(__FILE__, __LINE__, 0,
242 " Incompatible type of algorithm to compute the weight of"
243 " upwind.");
244
245 } /* Switch on the type of function to return */
246
247 return NULL; /* Avoid warning */
248 }
249
250 /*----------------------------------------------------------------------------*/
251 /*!
252 * \brief Retrieve the face edge id from a given cell edge id.
253 * Compute the weights related to each vertex of face from pre-computed
254 * quantities
255 *
256 * \param[in] cm pointer to a cs_cell_mesh_t structure
257 * \param[in] e edge id in the cell numbering
258 * \param[in] f face id in the cell numbering
259 * \param[in] tef value of the area of the triangles tef
260 * \param[in, out] wvf weights of vertices for the current face
261 *
262 * \return the face edge id
263 */
264 /*----------------------------------------------------------------------------*/
265
266 inline static short int
_set_fquant(const cs_cell_mesh_t * cm,short int e,short int f,const double * tef,double * wvf)267 _set_fquant(const cs_cell_mesh_t *cm,
268 short int e,
269 short int f,
270 const double *tef,
271 double *wvf)
272 {
273 short int ef = -1;
274
275 for (short int v = 0; v < cm->n_vc; v++)
276 wvf[v] = 0;
277
278 for (short int ee = cm->f2e_idx[f]; ee < cm->f2e_idx[f+1]; ee++) {
279
280 short int eab = cm->f2e_ids[ee];
281 short int va = cm->e2v_ids[2*eab];
282 short int vb = cm->e2v_ids[2*eab+1];
283
284 wvf[va] += tef[ee];
285 wvf[vb] += tef[ee];
286
287 if (eab == e) ef = ee;
288
289 } /* Loop on edges of f1 */
290
291 const double invf = 0.5/cm->face[f].meas;
292 for (short int v = 0; v < cm->n_vc; v++)
293 wvf[v] *= invf;
294
295 assert(ef != -1);
296 return ef;
297 }
298
299 /*----------------------------------------------------------------------------*/
300 /*!
301 * \brief Assemble a face matrix into a cell matrix
302 *
303 * \param[in] cm pointer to a cs_cell_mesh_t structure
304 * \param[in] fm pointer to a cs_face_mesh_t structure
305 * \param[in] af pointer to a cs_sdm_t structure related to a face
306 * \param[in, out] ac pointer to a cs_sdm_t structure related to a cell
307 */
308 /*----------------------------------------------------------------------------*/
309
310 static void
_assemble_face(const cs_cell_mesh_t * cm,const cs_face_mesh_t * fm,const cs_sdm_t * af,cs_sdm_t * ac)311 _assemble_face(const cs_cell_mesh_t *cm,
312 const cs_face_mesh_t *fm,
313 const cs_sdm_t *af,
314 cs_sdm_t *ac)
315 {
316 /* Add the face matrix to the cell matrix */
317 for (short int vi = 0; vi < fm->n_vf; vi++) {
318
319 const double *afi = af->val + af->n_rows*vi;
320
321 double *aci = ac->val + ac->n_rows*fm->v_ids[vi];
322 for (short int vj = 0; vj < fm->n_vf; vj++)
323 aci[fm->v_ids[vj]] += afi[vj]; /* (i,j) face --> cell */
324 aci[cm->n_vc] += afi[fm->n_vf]; /* (i,c) face --> cell */
325
326 }
327
328 const double *afc = af->val + af->n_rows*fm->n_vf;
329
330 double *acc = ac->val + ac->n_rows*cm->n_vc;
331 for (short int vj = 0; vj < fm->n_vf; vj++)
332 acc[fm->v_ids[vj]] += afc[vj]; /* (c,j) face --> cell */
333 acc[cm->n_vc] += afc[fm->n_vf]; /* (c,c) */
334 }
335
336 /*----------------------------------------------------------------------------*/
337 /*!
338 * \brief Define the local convection operator between primal vertices and
339 * dual faces. (Conservative formulation and Upwind flux)
340 *
341 * \param[in] cm pointer to a cs_cell_mesh_t structure
342 * \param[in] get_weight pointer to a function for evaluating upw. weight
343 * \param[in] fluxes array of fluxes on dual faces
344 * \param[in] upwcoef array of coefficient to the weight the upwinding
345 * \param[in, out] adv pointer to a cs_sdm_t structure
346 */
347 /*----------------------------------------------------------------------------*/
348
349 static void
_build_cell_epcd_upw(const cs_cell_mesh_t * cm,_upwind_weight_t * get_weight,const cs_real_t fluxes[],const cs_real_t upwcoef[],cs_sdm_t * adv)350 _build_cell_epcd_upw(const cs_cell_mesh_t *cm,
351 _upwind_weight_t *get_weight,
352 const cs_real_t fluxes[],
353 const cs_real_t upwcoef[],
354 cs_sdm_t *adv)
355 {
356 /* Loop on cell edges */
357 for (short int e = 0; e < cm->n_ec; e++) {
358
359 const short int sgn_v1 = cm->e2v_sgn[e];
360 const cs_real_t beta_flx = fluxes[e] * sgn_v1;
361
362 if (fabs(beta_flx) > 0) {
363
364 /* Compute the upwind coefficient knowing that fd(e),cd(v) = -v,e */
365 const double wv1 = get_weight(-sgn_v1 * upwcoef[e]);
366 const double c1mw = beta_flx * (1 - wv1);
367 const double cw = beta_flx * wv1;
368
369 /* Update local convection matrix.
370 Remember that sgn_v2 = -sgn_v1 */
371 short int v1 = cm->e2v_ids[2*e];
372 short int v2 = cm->e2v_ids[2*e+1];
373 assert(v1 != -1 && v2 != -1); /* Sanity check */
374
375 /* Update for the vertex v1 */
376 double *m1 = adv->val + v1*adv->n_rows;
377 m1[v1] += c1mw; /* diagonal term */
378 m1[v2] = -c1mw; /* extra-diag term */
379
380
381 /* Update for the vertex v2 */
382 double *m2 = adv->val + v2*adv->n_rows;
383 m2[v2] += -cw; /* diagonal term */
384 m2[v1] = cw; /* extra-diag term */
385
386 } /* convective flux is greater than zero */
387
388 } /* Loop on cell edges */
389
390 }
391
392 /*----------------------------------------------------------------------------*/
393 /*!
394 * \brief Define the local convection operator between primal vertices and
395 * dual faces. (Non-conservative formulation and centered flux)
396 *
397 * \param[in] cm pointer to a cs_cell_mesh_t structure
398 * \param[in] fluxes array of fluxes on dual faces
399 * \param[in, out] adv pointer to a cs_sdm_t structure
400 */
401 /*----------------------------------------------------------------------------*/
402
403 static void
_build_cell_epcd_cen(const cs_cell_mesh_t * cm,const cs_real_t fluxes[],cs_sdm_t * adv)404 _build_cell_epcd_cen(const cs_cell_mesh_t *cm,
405 const cs_real_t fluxes[],
406 cs_sdm_t *adv)
407 {
408 /* Weight is always equal to 0.5 Loop on cell edges */
409 for (short int e = 0; e < cm->n_ec; e++) {
410
411 const cs_real_t wflx = 0.5 * fluxes[e] * cm->e2v_sgn[e]; /* sgn_v1 */
412
413 if (fabs(wflx) > 0) {
414
415 /* Update local convection matrix */
416 short int v1 = cm->e2v_ids[2*e];
417 short int v2 = cm->e2v_ids[2*e+1];
418 assert(v1 != -1 && v2 != -1); /* Sanity check */
419
420 /* Update for the vertex v1 */
421 double *adv1 = adv->val + v1*adv->n_rows;
422 adv1[v1] += wflx; /* diagonal term */
423 adv1[v2] = -wflx; /* extra-diag term */
424
425 /* Update for the vertex v2
426 Use the fact that fd(e),cd(v) = -v,e and sgn_v1 = -sgn_v2 */
427 double *adv2 = adv->val + v2*adv->n_rows;
428 adv2[v2] += -wflx; /* diagonal term */
429 adv2[v1] = wflx; /* extra-diag term */
430
431 } /* convective flux is greater than zero */
432
433 } /* Loop on cell edges */
434
435 }
436
437 /*----------------------------------------------------------------------------*/
438 /*!
439 * \brief Define the local convection operator between primal vertices and
440 * dual faces. (Conservative formulation and Upwind flux)
441 *
442 * \param[in] cm pointer to a cs_cell_mesh_t structure
443 * \param[in] get_weight pointer to a function for evaluating upw. weight
444 * \param[in] fluxes array of fluxes on dual faces
445 * \param[in] upwcoef array of coefficient to the weight the upwinding
446 * \param[in, out] adv pointer to a cs_sdm_t structure
447 */
448 /*----------------------------------------------------------------------------*/
449
450 static void
_build_cell_vpfd_upw(const cs_cell_mesh_t * cm,_upwind_weight_t * get_weight,const cs_real_t fluxes[],const cs_real_t upwcoef[],cs_sdm_t * adv)451 _build_cell_vpfd_upw(const cs_cell_mesh_t *cm,
452 _upwind_weight_t *get_weight,
453 const cs_real_t fluxes[],
454 const cs_real_t upwcoef[],
455 cs_sdm_t *adv)
456 {
457 /* Loop on cell edges */
458 for (short int e = 0; e < cm->n_ec; e++) {
459
460 const cs_real_t beta_flx = fluxes[e];
461
462 if (fabs(beta_flx) > 0) {
463
464 short int sgn_v1 = cm->e2v_sgn[e];
465
466 /* Compute the upwind coefficient knowing that fd(e),cd(v) = -v,e */
467 const double wv1 = get_weight(-sgn_v1 * upwcoef[e]); /* in [0,1] */
468 const double cw1 = sgn_v1 * beta_flx * wv1;
469 const double cw2 = sgn_v1 * beta_flx * (1 - wv1);
470
471 /* Update local convection matrix */
472 short int v1 = cm->e2v_ids[2*e];
473 short int v2 = cm->e2v_ids[2*e+1];
474 assert(v1 != -1 && v2 != -1); /* Sanity check */
475
476 /* Update for the vertex v1 */
477 double *m1 = adv->val + v1*adv->n_rows;
478 m1[v1] += -cw1; /* diagonal term */
479 m1[v2] = -cw2; /* extra-diag term */
480
481 /* Update for the vertex v2 */
482 double *m2 = adv->val + v2*adv->n_rows;
483 m2[v2] += cw2; /* diagonal term */
484 m2[v1] = cw1; /* extra-diag term */
485
486 } /* convective flux is greater than zero */
487
488 } /* Loop on cell edges */
489
490 }
491
492 /*----------------------------------------------------------------------------*/
493 /*!
494 * \brief Define the local convection operator between primal vertices and
495 * dual faces. (Conservative formulation and centered flux)
496 *
497 * \param[in] cm pointer to a cs_cell_mesh_t structure
498 * \param[in] fluxes array of fluxes on dual faces
499 * \param[in, out] adv pointer to a cs_sdm_t structure
500 */
501 /*----------------------------------------------------------------------------*/
502
503 static void
_build_cell_vpfd_cen(const cs_cell_mesh_t * cm,const cs_real_t fluxes[],cs_sdm_t * adv)504 _build_cell_vpfd_cen(const cs_cell_mesh_t *cm,
505 const cs_real_t fluxes[],
506 cs_sdm_t *adv)
507 {
508 /* Weight is always equal to 0.5. Loop on cell edges */
509 for (short int e = 0; e < cm->n_ec; e++) {
510
511 const cs_real_t wflx = 0.5*fluxes[e]*cm->e2v_sgn[e];
512
513 if (fabs(wflx) > 0) { /* Update local convection matrix */
514
515 short int v1 = cm->e2v_ids[2*e];
516 short int v2 = cm->e2v_ids[2*e+1];
517 assert(v1 != -1 && v2 != -1); /* Sanity check */
518
519 /* Update for the vertex v1 */
520 double *adv1 = adv->val + v1*adv->n_rows;
521 adv1[v1] += -wflx; /* diagonal term */
522 adv1[v2] = -wflx; /* extra-diag term */
523
524 /* Update for the vertex v2 (sgn_v2 = -sgn_v1) */
525 double *adv2 = adv->val + v2*adv->n_rows;
526 adv2[v2] += wflx; /* diagonal term */
527 adv2[v1] = wflx; /* extra-diag term */
528
529 } /* convective flux is greater than zero */
530
531 } /* Loop on cell edges */
532
533 }
534
535 /*----------------------------------------------------------------------------*/
536 /*!
537 * \brief Define the local convection operator between primal vertices and
538 * dual faces. (Conservative formulation and mixed centered/upwind)
539 *
540 * \param[in] cm pointer to a cs_cell_mesh_t structure
541 * \param[in] upwp portion of upwinding to apply
542 * \param[in] fluxes array of fluxes on dual faces
543 * \param[in, out] adv pointer to a cs_sdm_t structure
544 */
545 /*----------------------------------------------------------------------------*/
546
547 static void
_build_cell_vpfd_mcu(const cs_cell_mesh_t * cm,const cs_real_t upwp,const cs_real_t fluxes[],cs_sdm_t * adv)548 _build_cell_vpfd_mcu(const cs_cell_mesh_t *cm,
549 const cs_real_t upwp,
550 const cs_real_t fluxes[],
551 cs_sdm_t *adv)
552 {
553 /* Loop on cell edges */
554 for (short int e = 0; e < cm->n_ec; e++) {
555
556 const cs_real_t beta_flx = fluxes[e];
557
558 if (fabs(beta_flx) > 0) {
559
560 short int sgn_v1 = cm->e2v_sgn[e];
561
562 /* Compute the upwind coefficient knowing that fd(e),cd(v) = -v,e */
563 const double wup1 = _get_upwind_weight(-sgn_v1*beta_flx); /* 0 or 1 */
564 /* (1-upwp)*0.5 --> centered scheme contribution
565 upwp * wup1 --> upwind scheme contribution */
566 const double wv1 = wup1 * upwp + 0.5*(1 - upwp);
567 const double cw1 = sgn_v1 * beta_flx * wv1;
568 const double cw2 = sgn_v1 * beta_flx * (1 - wv1);
569
570 /* Update local convection matrix */
571 short int v1 = cm->e2v_ids[2*e];
572 short int v2 = cm->e2v_ids[2*e+1];
573 assert(v1 != -1 && v2 != -1); /* Sanity check */
574
575 /* Update for the vertex v1 */
576 double *m1 = adv->val + v1*adv->n_rows;
577 m1[v1] += -cw1; /* diagonal term */
578 m1[v2] = -cw2; /* extra-diag term */
579
580 /* Update for the vertex v2 */
581 double *m2 = adv->val + v2*adv->n_rows;
582 m2[v2] += cw2; /* diagonal term */
583 m2[v1] = cw1; /* extra-diag term */
584
585 } /* convective flux is greater than zero */
586
587 } /* Loop on cell edges */
588
589 }
590
591
592 /*----------------------------------------------------------------------------*/
593 /*!
594 * \brief Update the cellwise system with the contributions arising from
595 * the boundary conditions
596 *
597 * \param[in] beta_nf advection field coefficient
598 * \param[in] fm pointer to a cs_face_mesh_t structure
599 * \param[in] matf pointer to a local dense matrix related to a face
600 * \param[in, out] csys pointer to a cs_cell_sys_t structure
601 */
602 /*----------------------------------------------------------------------------*/
603
604 static void
_update_vcb_system_with_bc(const cs_real_t beta_nf,const cs_face_mesh_t * fm,const cs_sdm_t * matf,cs_cell_sys_t * csys)605 _update_vcb_system_with_bc(const cs_real_t beta_nf,
606 const cs_face_mesh_t *fm,
607 const cs_sdm_t *matf,
608 cs_cell_sys_t *csys)
609 {
610 double _dirf[10], _rhsf[10];
611 double *dirf = NULL, *rhsf = NULL;
612
613 if (csys->has_dirichlet) {
614
615 if (fm->n_vf > 10) {
616 BFT_MALLOC(dirf, 2*fm->n_vf, double);
617 rhsf = dirf + fm->n_vf;
618 }
619 else {
620 dirf = _dirf;
621 rhsf = _rhsf;
622 }
623
624 /* Compute the Dirichlet contribution to RHS */
625 for (short int vfi = 0; vfi < fm->n_vf; vfi++)
626 dirf[vfi] = beta_nf * csys->dir_values[fm->v_ids[vfi]];
627 cs_sdm_square_matvec(matf, dirf, rhsf);
628
629 /* Update RHS */
630 for (short int vfi = 0; vfi < fm->n_vf; vfi++)
631 csys->rhs[fm->v_ids[vfi]] += rhsf[vfi];
632
633 } /* There is at least one dirichlet */
634
635 /* Update cellwise matrix */
636 const int n_cell_dofs = csys->mat->n_rows;
637 double *matc = csys->mat->val;
638
639 for (short int vfi = 0; vfi < fm->n_vf; vfi++) {
640
641 const double *mfi = matf->val + vfi*fm->n_vf;
642 double *mi = matc + fm->v_ids[vfi]*n_cell_dofs;
643
644 for (short int vfj = 0; vfj < fm->n_vf; vfj++) {
645 const short int vj = fm->v_ids[vfj];
646 mi[vj] += beta_nf * mfi[vfj];
647 }
648
649 } /* Loop on face vertices */
650
651 if (dirf != _dirf)
652 BFT_FREE(dirf);
653 }
654
655 /*----------------------------------------------------------------------------*/
656 /*!
657 * \brief Compute the consistent part of the convection operator attached to
658 * a cell with a CDO vertex+cell-based scheme when the advection is
659 * considered as constant inside a cell
660 *
661 * \param[in] adv_cell constant vector field inside the current cell
662 * \param[in] cm pointer to a cs_cell_mesh_t structure
663 * \param[in] fm pointer to a cs_face_mesh_t structure
664 * \param[in, out] cb pointer to a cs_cell_builder_t structure
665 */
666 /*----------------------------------------------------------------------------*/
667
668 static void
_vcb_cellwise_consistent_part(const cs_nvec3_t adv_cell,const cs_cell_mesh_t * cm,const cs_face_mesh_t * fm,cs_cell_builder_t * cb)669 _vcb_cellwise_consistent_part(const cs_nvec3_t adv_cell,
670 const cs_cell_mesh_t *cm,
671 const cs_face_mesh_t *fm,
672 cs_cell_builder_t *cb)
673 {
674 cs_real_3_t grd_v1, grd_v2, grd_c;
675
676 const int n_sysf = fm->n_vf + 1;
677 const short int fshift = cm->f2e_idx[fm->f_id];
678
679 /* Temporary buffers storing useful quantities to build the consistent part */
680 cs_sdm_t *af = cb->aux;
681
682 /* tef_save is not set here but its length is 2*n_ec */
683 double *bgc_save = cb->values; /* size = n_fc */
684 double *l_vc = cb->values + cm->n_fc + 2*cm->n_ec; /* size = n_vc */
685
686 cs_real_3_t *bgvf = cb->vectors + fshift; /* size 2*n_ec */
687 cs_real_3_t *u_vc = cb->vectors + 2*cm->n_ec; /* size n_vc */
688
689 /* Useful quantities are stored in cb->values and cb->vectors */
690 const double hf_coef = cs_math_1ov3 * fm->hfc;
691
692 /* Set the consistent part for already known part.
693 Consistent part (c,c) contribution: sum_(f \in F_c) |pfc| bgc
694 Consistent part (i,c) contribution: sum_(f \in F_c) 0.75 * wif * |pfc| bgc
695 */
696
697 /* Compute the gradient of the Lagrange function related to xc which is
698 constant inside p_{f,c} */
699 cs_compute_grdfc_fw(fm, grd_c);
700
701 const double bgc = _dp3(grd_c, adv_cell.unitv);
702 const double pfc_bgc = adv_cell.meas * fm->pvol * bgc;
703
704 bgc_save[fm->f_id] = bgc; /* Store it for a future use */
705
706 af->val[n_sysf*fm->n_vf + fm->n_vf] = 0.25 * pfc_bgc; /* (c,c) */
707 for (short int v = 0; v < fm->n_vf; v++)
708 af->val[n_sysf*v + fm->n_vf] = 0.75 * fm->wvf[v] * pfc_bgc; /* (i,c) */
709
710 /* Compute xc --> xv length and unit vector for all face vertices */
711 for (short int v = 0; v < fm->n_vf; v++)
712 cs_math_3_length_unitv(fm->xc, fm->xv + 3*v, l_vc + v, u_vc[v]);
713
714 /* Compute the gradient of Lagrange functions related to the partition
715 into tetrahedron of p_(fc) */
716 for (short int e = 0; e < fm->n_ef; e++) {
717
718 const double pef_coef = 0.25 * hf_coef * fm->tef[e] * adv_cell.meas;
719 const short int v1 = fm->e2v_ids[2*e];
720 const short int v2 = fm->e2v_ids[2*e+1];
721
722 /* Gradient of the Lagrange function related to v1 and v2 */
723 cs_compute_grd_ve(v1, v2, fm->dedge, (const cs_real_t (*)[3])u_vc, l_vc,
724 grd_v1, grd_v2);
725
726 const double bgv1 = _dp3(grd_v1, adv_cell.unitv);
727 const double bgv2 = _dp3(grd_v2, adv_cell.unitv);
728
729 /* Gradient of the Lagrange function related to a face.
730 This formula is a consequence of the Partition of the Unity
731 grd_f = -( grd_c + grd_v1 + grd_v2) */
732 const double bgf = -(bgc + bgv1 + bgv2);
733
734 for (short int vi = 0; vi < fm->n_vf; vi++) {
735
736 double *afi = af->val + n_sysf*vi;
737 double lvci_part = fm->wvf[vi];
738 if (vi == v1 || vi == v2)
739 lvci_part += 1;
740 const double lvci = pef_coef * lvci_part;
741
742 for (short int vj = 0; vj < fm->n_vf; vj++) {
743
744 double glj = fm->wvf[vj]*bgf;
745 if (vj == v1)
746 glj += bgv1;
747 else if (vj == v2)
748 glj += bgv2;
749
750 afi[vj] += glj * lvci; /* consistent part (i,j) face mat. */
751
752 } /* Loop on vj in Vf */
753
754 } /* Loop on vi in Vf */
755
756 /* (c,j) entries */
757 double *afc = af->val + n_sysf*fm->n_vf;
758 for (short int vj = 0; vj < fm->n_vf; vj++) {
759
760 double glj = fm->wvf[vj]*bgf;
761 if (vj == v1)
762 glj += bgv1;
763 else if (vj == v2)
764 glj += bgv2;
765
766 afc[vj] += glj * pef_coef; /* consistent part (c,j) face mat. */
767
768 } /* Loop on vj in Vf */
769
770 /* Store bgv1, bgv2, bgf */
771 bgvf[e][0] = bgv1;
772 bgvf[e][1] = bgv2;
773 bgvf[e][2] = bgf;
774
775 } /* Loop on face edges */
776
777 }
778
779 /*----------------------------------------------------------------------------*/
780 /*!
781 * \brief Compute the consistent part of the convection operator attached to
782 * a cell with a CDO vertex+cell-based scheme when the advection field
783 * is not considered constant inside the current cell
784 *
785 * \param[in] adv_field pointer to a cs_adv_field_t structure
786 * \param[in] adv_cell constant vector field inside the current cell
787 * \param[in] cm pointer to a cs_cell_mesh_t structure
788 * \param[in] fm pointer to a cs_face_mesh_t structure
789 * \param[in] t_eval time at which one evaluates the advection field
790 * \param[in, out] cb pointer to a cs_cell_builder_t structure
791 */
792 /*----------------------------------------------------------------------------*/
793
794 static void
_vcb_consistent_part(const cs_adv_field_t * adv_field,const cs_nvec3_t adv_cell,const cs_cell_mesh_t * cm,const cs_face_mesh_t * fm,cs_real_t t_eval,cs_cell_builder_t * cb)795 _vcb_consistent_part(const cs_adv_field_t *adv_field,
796 const cs_nvec3_t adv_cell,
797 const cs_cell_mesh_t *cm,
798 const cs_face_mesh_t *fm,
799 cs_real_t t_eval,
800 cs_cell_builder_t *cb)
801 {
802 cs_real_3_t grd_v1, grd_v2, grd_c, xg;
803 cs_nvec3_t bnv;
804
805 const short int fshift = cm->f2e_idx[fm->f_id];
806 const int n_sysf = fm->n_vf + 1;
807 const cs_nvec3_t deq = fm->dedge;
808 const cs_quant_t pfq = fm->face;
809 const double hf_coef = cs_math_1ov3 * cm->hfc[fm->f_id];
810
811 /* Useful quantities are stored in cb->values and cb->vectors */
812 cs_sdm_t *af = cb->aux;
813 double *bgc_save = cb->values; /* size n_fc */
814 double *l_vc = cb->values + cm->n_fc + 2*cm->n_ec; /* size n_vc */
815 cs_real_3_t *bgvf = cb->vectors + fshift;
816 cs_real_3_t *u_vc = cb->vectors + 2*cm->n_ec;
817
818 /* Set the consistent part for already known part.
819 Consistent part (c,c) contribution: sum_(f \in F_c) |pfc| bgc
820 Consistent part (i,c) contribution: sum_(f \in F_c) 0.75 * wif * |pfc| bgc
821 */
822
823 /* Compute the gradient of the Lagrange function related to xc which is
824 constant inside p_{f,c} */
825 const cs_real_t ohf = -fm->f_sgn/cm->hfc[fm->f_id];
826 for (int k = 0; k < 3; k++) grd_c[k] = ohf * pfq.unitv[k];
827
828 const double bgcc = _dp3(grd_c, adv_cell.unitv);
829
830 bgc_save[fm->f_id] = bgcc; /* Store it for using it in the stabilization */
831
832 /* Compute xc --> xv length and unit vector for all face vertices */
833 for (short int v = 0; v < fm->n_vf; v++)
834 cs_math_3_length_unitv(fm->xc, fm->xv + 3*v, l_vc + v, u_vc[v]);
835
836 /* Compute the gradient of Lagrange functions related to the partition
837 into tetrahedron of p_(fc) */
838 for (short int e = 0; e < fm->n_ef; e++) {
839
840 const short int v1 = fm->e2v_ids[2*e];
841 const short int v2 = fm->e2v_ids[2*e+1];
842
843 for (int k = 0; k < 3; k++)
844 xg[k] = 0.25 * (fm->xv[3*v1+k]+fm->xv[3*v2+k]+fm->xc[k]+pfq.center[k]);
845 cs_advection_field_cw_eval_at_xyz(adv_field, cm, xg, t_eval, &bnv);
846
847 const double bgc = _dp3(grd_c, bnv.unitv);
848 const double pef_coef = 0.25 * bnv.meas * hf_coef * fm->tef[e];
849
850 /* Gradient of the Lagrange function related to v1 and v2 */
851 cs_compute_grd_ve(v1, v2, deq, (const cs_real_t (*)[3])u_vc, l_vc,
852 grd_v1, grd_v2);
853
854 const double bgv1 = _dp3(grd_v1, bnv.unitv);
855 const double bgv2 = _dp3(grd_v2, bnv.unitv);
856
857 /* Gradient of the Lagrange function related to a face.
858 This formula is a consequence of the Partition of the Unity
859 grd_f = -( grd_c + grd_v1 + grd_v2) */
860 const double bgf = -(bgc + bgv1 + bgv2);
861
862 for (short int vi = 0; vi < fm->n_vf; vi++) {
863
864 double *afi = af->val + n_sysf*vi;
865 double lvci_part = fm->wvf[vi];
866 if (vi == v1 || vi == v2)
867 lvci_part += 1;
868 const double lvci = pef_coef * lvci_part;
869
870 for (short int vj = 0; vj < fm->n_vf; vj++) {
871
872 double glj = fm->wvf[vj]*bgf;
873 if (vj == v1)
874 glj += bgv1;
875 else if (vj == v2)
876 glj += bgv2;
877
878 afi[vj] += glj * lvci; /* consistent part (i,j) face mat. */
879
880 } /* Loop on vj in Vf */
881
882 /* (i, c) entries */
883 afi[fm->n_vf] += lvci * bgc;
884
885 } /* Loop on vi in Vf */
886
887 /* (c,j) entries */
888 double *afc = af->val + n_sysf*fm->n_vf;
889 for (short int vj = 0; vj < fm->n_vf; vj++) {
890
891 double glj = fm->wvf[vj]*bgf;
892 if (vj == v1)
893 glj += bgv1;
894 else if (vj == v2)
895 glj += bgv2;
896
897 afc[vj] += glj * pef_coef; /* consistent part (c,j) face mat. */
898
899 } /* Loop on vj in Vf */
900
901 /* (c,c) entries */
902 afc[fm->n_vf] += pef_coef * bgc;
903
904 /* Store bgv1, bgv2, bgf */
905 bgvf[e][0] = _dp3(grd_v1, adv_cell.unitv);
906 bgvf[e][1] = _dp3(grd_v2, adv_cell.unitv);
907 /* This formula is a consequence of the Partition of the Unity
908 grd_f = -( grd_c + grd_v1 + grd_v2) */
909 bgvf[e][2] = -(bgcc + bgvf[e][0] + bgvf[e][1]);
910
911 } /* Loop on face edges */
912
913 }
914
915 /*----------------------------------------------------------------------------*/
916 /*!
917 * \brief Compute the stabilization part of the convection operator attached
918 * to a cell with a CDO vertex+cell-based scheme (inside pfc)
919 *
920 * \param[in] cm pointer to a cs_cell_mesh_t structure
921 * \param[in] fm pointer to a cs_face_mesh_t structure
922 * \param[in] stab_coef default value of the stabilization coefficient
923 * \param[in, out] cb pointer to a cs_cell_builder_t structure
924 */
925 /*----------------------------------------------------------------------------*/
926
927 static void
_vcb_stabilization_part1(const cs_cell_mesh_t * cm,const cs_face_mesh_t * fm,const double stab_coef,cs_cell_builder_t * cb)928 _vcb_stabilization_part1(const cs_cell_mesh_t *cm,
929 const cs_face_mesh_t *fm,
930 const double stab_coef,
931 cs_cell_builder_t *cb)
932 {
933 const short int n_sysf = fm->n_vf + 1;
934 const short int fshift = cm->f2e_idx[fm->f_id];
935
936 cs_real_3_t *bgvf = cb->vectors + fshift;
937 cs_sdm_t *af = cb->aux;
938
939 /* First part of the stabilization (inside each face) */
940 for (short int e1 = 0; e1 < fm->n_ef; e1++) {
941
942 short int e2 = e1 - 1;
943 if (e1 == 0)
944 e2 = fm->n_ef - 1; /* Last edge face */
945
946 short int v1e1 = fm->e2v_ids[2*e1];
947 short int v2e1 = fm->e2v_ids[2*e1+1];
948 short int v1e2 = fm->e2v_ids[2*e2]; /* v1_prev */
949 short int v2e2 = fm->e2v_ids[2*e2+1]; /* v2_prev */
950
951 const double jump_bgf = bgvf[e1][2] - bgvf[e2][2];
952
953 /* Area of the triangle which is the interface between the two pefc */
954 double svfc = stab_coef;
955 if (v1e1 == v1e2 || v1e1 == v2e2)
956 svfc *= cs_math_surftri(fm->xv + 3*v1e1, fm->face.center, fm->xc);
957 else
958 svfc *= cs_math_surftri(fm->xv + 3*v2e1, fm->face.center, fm->xc);
959
960 /* Update convection matrix related to the current face
961 (symmetric contributions) */
962 for (short int vi = 0; vi < fm->n_vf; vi++) {
963
964 double *afi = af->val + n_sysf*vi;
965 double jump_i = fm->wvf[vi] * jump_bgf;
966
967 if (vi == v1e1)
968 jump_i += bgvf[e1][0];
969 else if (vi == v2e1)
970 jump_i += bgvf[e1][1];
971
972 if (vi == v1e2)
973 jump_i -= bgvf[e2][0];
974 else if (vi == v2e2)
975 jump_i -= bgvf[e2][1];
976
977 const double coef_i = svfc * jump_i;
978
979 afi[vi] += coef_i * jump_i; /* Stab. (i,i) face mat. */
980
981 for (short int vj = vi + 1; vj < fm->n_vf; vj++) {
982
983 double jump_j = fm->wvf[vj] * jump_bgf;
984
985 if (vj == v1e1)
986 jump_j += bgvf[e1][0];
987 else if (vj == v2e1)
988 jump_j += bgvf[e1][1];
989
990 if (vj == v1e2)
991 jump_j -= bgvf[e2][0];
992 else if (vj == v2e2)
993 jump_j -= bgvf[e2][1];
994
995 const double coef_ij = coef_i * jump_j;
996
997 afi[vj] += coef_ij; /* Stab. (i,j) face mat. */
998 af->val[vj*n_sysf + vi] += coef_ij; /* Stab. (j,i) face mat. symm. */
999
1000 } /* Loop on vj in Vf */
1001 } /* Loop on vi in Vf */
1002
1003 } /* Loop on edge faces */
1004
1005 }
1006
1007 /*----------------------------------------------------------------------------*/
1008 /*!
1009 * \brief Compute the stabilization part of the convection operator attached
1010 * to a cell with a CDO vertex+cell-based scheme (between pfc)
1011 *
1012 * \param[in] cm pointer to a cs_cell_mesh_t structure
1013 * \param[in] stab_coef value of the stabilization coefficient
1014 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1015 */
1016 /*----------------------------------------------------------------------------*/
1017
1018 static void
_vcb_stabilization_part2(const cs_cell_mesh_t * cm,const double stab_coef,cs_cell_builder_t * cb)1019 _vcb_stabilization_part2(const cs_cell_mesh_t *cm,
1020 const double stab_coef,
1021 cs_cell_builder_t *cb)
1022 {
1023 const short int n_sysc = cm->n_vc + 1;
1024
1025 cs_sdm_t *a = cb->loc;
1026
1027 /* Temporary buffers used to store pre-computed data */
1028 double *bgc_save, *tef_save, *wvf1, *wvf2; /* scalar-valued buffers */
1029 int m_shft = 0;
1030
1031 bgc_save = cb->values, m_shft = cm->n_fc; /* size = n_fc */
1032 tef_save = cb->values + m_shft, m_shft += 2*cm->n_ec; /* size = 2*n_ec */
1033 wvf1 = cb->values + m_shft, m_shft += cm->n_vc; /* size = n_vc */
1034 wvf2 = cb->values + m_shft, m_shft += cm->n_vc; /* size = n_vc */
1035
1036 cs_real_3_t *bgvf_save = cb->vectors; /* size = 2*n_ec */
1037
1038 for (short int e = 0; e < cm->n_ec; e++) {
1039
1040 const double tec = stab_coef * cs_compute_area_from_quant(cm->edge[e],
1041 cm->xc);
1042 const short int eshift = 2*e;
1043 const short int v1 = cm->e2v_ids[eshift];
1044 const short int v2 = cm->e2v_ids[eshift+1];
1045 const short int _v = (v1 < v2) ? 0 : 1;
1046 const short int f1 = cm->e2f_ids[eshift];
1047 const short int f2 = cm->e2f_ids[eshift+1];
1048 const double jump_c = bgc_save[f2] - bgc_save[f1];
1049
1050 /* (c, c) contrib */
1051 double *ac = a->val + cm->n_vc*n_sysc;
1052 ac[cm->n_vc] += tec * jump_c * jump_c;
1053
1054 const short int ef1 = _set_fquant(cm, e, f1, tef_save, wvf1);
1055 const short int ef2 = _set_fquant(cm, e, f2, tef_save, wvf2);
1056
1057 for (short int vi = 0; vi < cm->n_vc; vi++) {
1058 if (wvf2[vi] + wvf1[vi] > 0) { /* vi belongs at least to f1 or f2 */
1059
1060 double *ai = a->val + vi*n_sysc;
1061 double jump_i =
1062 wvf2[vi]*bgvf_save[ef2][2] - wvf1[vi]*bgvf_save[ef1][2];
1063
1064 if (vi == v1)
1065 jump_i += bgvf_save[ef2][_v] - bgvf_save[ef1][_v];
1066 else if (vi == v2)
1067 jump_i += bgvf_save[ef2][1-_v] - bgvf_save[ef1][1-_v];
1068
1069 const double coef_i = tec * jump_i;
1070
1071 /* (i, i) contrib */
1072 ai[vi] += jump_i * coef_i;
1073
1074 for (short int vj = vi + 1; vj < cm->n_vc; vj++) {
1075 if (wvf2[vj] + wvf1[vj] > 0) { /* vj belongs at least to f1 or f2 */
1076
1077 double jump_j =
1078 wvf2[vj]*bgvf_save[ef2][2] - wvf1[vj]*bgvf_save[ef1][2];
1079 if (vj == v1)
1080 jump_j += bgvf_save[ef2][_v] - bgvf_save[ef1][_v];
1081 else if (vj == v2)
1082 jump_j += bgvf_save[ef2][1-_v] - bgvf_save[ef1][1-_v];
1083
1084 const double coef_ij = coef_i * jump_j;
1085 ai[vj] += coef_ij;
1086 a->val[vj*n_sysc + vi] += coef_ij; /* symmetric */
1087
1088 } /* vj belongs to f1 or f2 */
1089 } /* vj loop */
1090
1091 /* (i, c) contrib */
1092 const double coef_ic = coef_i * jump_c;
1093 ai[cm->n_vc] += coef_ic;
1094 ac[vi] += coef_ic; /* symmetric */
1095
1096 } /* vi belongs to f1 or f2 */
1097 } /* vi loop */
1098
1099 } /* End of loop on cell edges */
1100
1101 }
1102
1103 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1104
1105 /*============================================================================
1106 * Public function prototypes
1107 *============================================================================*/
1108
1109 /*----------------------------------------------------------------------------*/
1110 /*!
1111 * \brief Set the value of the stabilization coefficient used in CIP scheme
1112 *
1113 * \param[in] new_value value of the stabilization coefficient
1114 */
1115 /*----------------------------------------------------------------------------*/
1116
1117 void
cs_cdo_advection_set_cip_coef(double new_value)1118 cs_cdo_advection_set_cip_coef(double new_value)
1119 {
1120 cs_cip_stab_coef = new_value;
1121 }
1122
1123 /*----------------------------------------------------------------------------*/
1124 /*!
1125 * \brief Get the value of the stabilization coefficient used in CIP scheme
1126 *
1127 * \return the value the stabilization coefficient
1128 */
1129 /*----------------------------------------------------------------------------*/
1130
1131 double
cs_cdo_advection_get_cip_coef(void)1132 cs_cdo_advection_get_cip_coef(void)
1133 {
1134 return cs_cip_stab_coef;
1135 }
1136
1137 /*----------------------------------------------------------------------------*/
1138 /*!
1139 * \brief Perform preprocessing such as the computation of the advection flux
1140 * at the expected location in order to be able to build the advection
1141 * matrix. Follow the prototype given by cs_cdofb_adv_open_hook_t
1142 * Default case.
1143 *
1144 * \param[in] eqp pointer to a cs_equation_param_t structure
1145 * \param[in] cm pointer to a cs_cell_mesh_t structure
1146 * \param[in] csys pointer to a cs_cell_sys_t structure
1147 * \param[in, out] input NULL or pointer to a structure cast on-the-fly
1148 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1149 */
1150 /*----------------------------------------------------------------------------*/
1151
1152 void
cs_cdofb_advection_open_default(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,void * input,cs_cell_builder_t * cb)1153 cs_cdofb_advection_open_default(const cs_equation_param_t *eqp,
1154 const cs_cell_mesh_t *cm,
1155 const cs_cell_sys_t *csys,
1156 void *input,
1157 cs_cell_builder_t *cb)
1158 {
1159 CS_UNUSED(csys);
1160 CS_UNUSED(input);
1161
1162 assert(eqp->adv_extrapol == CS_PARAM_ADVECTION_EXTRAPOL_NONE);
1163
1164 /* Compute the flux across the primal faces. Store in cb->adv_fluxes */
1165 cs_advection_field_cw_face_flux(cm, eqp->adv_field, cb->t_bc_eval,
1166 cb->adv_fluxes);
1167
1168 if (eqp->adv_scaling_property != NULL) {
1169
1170 cs_real_t scaling = eqp->adv_scaling_property->ref_value;
1171 if (cs_property_is_uniform(eqp->adv_scaling_property) == false)
1172 scaling = cs_property_value_in_cell(cm,
1173 eqp->adv_scaling_property,
1174 cb->t_pty_eval);
1175
1176 for (int f = 0; f < cm->n_fc; f++)
1177 cb->adv_fluxes[f] *= scaling;
1178
1179 } /* Apply a scaling factor to the advective flux */
1180 }
1181
1182 /*----------------------------------------------------------------------------*/
1183 /*!
1184 * \brief Operation done after the matrix related to the advection term has
1185 * been defined.
1186 * Follow the prototype given by cs_cdofb_adv_close_hook_t
1187 * Default scalar-valued case.
1188 *
1189 * \param[in] eqp pointer to a cs_equation_param_t structure
1190 * \param[in] cm pointer to a cs_cell_mesh_t structure
1191 * \param[in, out] csys pointer to a cs_cell_sys_t structure
1192 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1193 * \param[in, out] adv pointer to the local advection matrix
1194 */
1195 /*----------------------------------------------------------------------------*/
1196
1197 void
cs_cdofb_advection_close_default_scal(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1198 cs_cdofb_advection_close_default_scal(const cs_equation_param_t *eqp,
1199 const cs_cell_mesh_t *cm,
1200 cs_cell_sys_t *csys,
1201 cs_cell_builder_t *cb,
1202 cs_sdm_t *adv)
1203 {
1204 CS_UNUSED(eqp);
1205 CS_UNUSED(cm);
1206 CS_UNUSED(adv);
1207
1208 cs_sdm_add(csys->mat, cb->loc);
1209 }
1210
1211 /*----------------------------------------------------------------------------*/
1212 /*!
1213 * \brief Operation done after the matrix related to the advection term has
1214 * been defined.
1215 * Follow the prototype given by cs_cdofb_adv_close_hook_t
1216 * Default vector-valued case.
1217 *
1218 * \param[in] eqp pointer to a cs_equation_param_t structure
1219 * \param[in] cm pointer to a cs_cell_mesh_t structure
1220 * \param[in, out] csys pointer to a cs_cell_sys_t structure
1221 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1222 * \param[in, out] adv pointer to the local advection matrix
1223 */
1224 /*----------------------------------------------------------------------------*/
1225
1226 void
cs_cdofb_advection_close_default_vect(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1227 cs_cdofb_advection_close_default_vect(const cs_equation_param_t *eqp,
1228 const cs_cell_mesh_t *cm,
1229 cs_cell_sys_t *csys,
1230 cs_cell_builder_t *cb,
1231 cs_sdm_t *adv)
1232 {
1233 CS_UNUSED(eqp);
1234 CS_UNUSED(cb);
1235
1236 /* Add the local scalar-valued advection operator to the local vector-valued
1237 system */
1238 const cs_real_t *sval = adv->val;
1239 for (int bi = 0; bi < cm->n_fc + 1; bi++) {
1240 for (int bj = 0; bj < cm->n_fc + 1; bj++) {
1241
1242 /* Retrieve the 3x3 matrix */
1243 cs_sdm_t *bij = cs_sdm_get_block(csys->mat, bi, bj);
1244 assert(bij->n_rows == bij->n_cols && bij->n_rows == 3);
1245
1246 const cs_real_t _val = sval[(cm->n_fc+1)*bi+bj];
1247 bij->val[0] += _val;
1248 bij->val[4] += _val;
1249 bij->val[8] += _val;
1250
1251 }
1252 }
1253
1254 }
1255
1256 /*----------------------------------------------------------------------------*/
1257 /*!
1258 * \brief Operation done after the matrix related to the advection term has
1259 * been defined.
1260 * Follow the prototype given by cs_cdofb_adv_close_hook_t
1261 * Explicit treatment without extrapolation for scalar-valued DoFs.
1262 *
1263 * \param[in] eqp pointer to a cs_equation_param_t structure
1264 * \param[in] cm pointer to a cs_cell_mesh_t structure
1265 * \param[in, out] csys pointer to a cs_cell_sys_t structure
1266 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1267 * \param[in, out] adv pointer to the local advection matrix
1268 */
1269 /*----------------------------------------------------------------------------*/
1270
1271 void
cs_cdofb_advection_close_exp_none_scal(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1272 cs_cdofb_advection_close_exp_none_scal(const cs_equation_param_t *eqp,
1273 const cs_cell_mesh_t *cm,
1274 cs_cell_sys_t *csys,
1275 cs_cell_builder_t *cb,
1276 cs_sdm_t *adv)
1277 {
1278 CS_UNUSED(eqp);
1279
1280 /* Update the RHS: u_n is the previous time step.
1281 * This is done before the static condensation. Thus, there is cell unknown
1282 */
1283 double *adv_u_n = cb->values;
1284
1285 cs_sdm_matvec(adv, csys->val_n, adv_u_n);
1286
1287 /* Update the RHS (interlaced values) */
1288 for (int i = 0; i < cm->n_fc + 1; i++)
1289 csys->rhs[i] -= adv_u_n[i];
1290 }
1291
1292 /*----------------------------------------------------------------------------*/
1293 /*!
1294 * \brief Operation done after the matrix related to the advection term has
1295 * been defined.
1296 * Follow the prototype given by cs_cdofb_adv_close_hook_t
1297 * Explicit treatment without extrapolation for vector-valued DoFs.
1298 *
1299 * \param[in] eqp pointer to a cs_equation_param_t structure
1300 * \param[in] cm pointer to a cs_cell_mesh_t structure
1301 * \param[in, out] csys pointer to a cs_cell_sys_t structure
1302 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1303 * \param[in, out] adv pointer to the local advection matrix
1304 */
1305 /*----------------------------------------------------------------------------*/
1306
1307 void
cs_cdofb_advection_close_exp_none_vect(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1308 cs_cdofb_advection_close_exp_none_vect(const cs_equation_param_t *eqp,
1309 const cs_cell_mesh_t *cm,
1310 cs_cell_sys_t *csys,
1311 cs_cell_builder_t *cb,
1312 cs_sdm_t *adv)
1313 {
1314 CS_UNUSED(eqp);
1315 assert(eqp->dim == 3);
1316
1317 /* Update the RHS component by component: u_n is the previous time step.
1318 * This is done before the static condensation. Thus, there is cell unknown
1319 */
1320 double *u_n = cb->values;
1321 double *adv_u_n = cb->values + cm->n_fc + 1;
1322
1323 for (int k = 0; k < 3; k++) {
1324
1325 /* De-interlace the local variable in order to perform a matrix-vector
1326 product */
1327 for (int i = 0; i < cm->n_fc + 1; i++)
1328 u_n[i] = csys->val_n[3*i+k];
1329
1330 cs_sdm_matvec(adv, u_n, adv_u_n);
1331
1332 /* Update the RHS (interlaced values) */
1333 for (int i = 0; i < cm->n_fc + 1; i++)
1334 csys->rhs[3*i+k] -= adv_u_n[i];
1335
1336 }
1337
1338 }
1339
1340 /*----------------------------------------------------------------------------*/
1341 /*!
1342 * \brief Build the cellwise advection operator for CDO-Fb schemes
1343 * The local matrix related to this operator is stored in cb->loc
1344 *
1345 * Case of an advection term without a diffusion operator. In this
1346 * situation, a numerical issue may arise if an internal or a border
1347 * face is such that there is no advective flux. A specil treatment
1348 * is performed to tackle this issue.
1349 *
1350 * \param[in] eqp pointer to a cs_equation_param_t structure
1351 * \param[in] cm pointer to a cs_cell_mesh_t structure
1352 * \param[in] csys pointer to a cs_cell_sys_t structure
1353 * \param[in] scheme_func pointer to the function building the system
1354 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1355 */
1356 /*----------------------------------------------------------------------------*/
1357
1358 void
cs_cdofb_advection_no_diffusion(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,cs_cdofb_adv_scheme_t * scheme_func,cs_cell_builder_t * cb)1359 cs_cdofb_advection_no_diffusion(const cs_equation_param_t *eqp,
1360 const cs_cell_mesh_t *cm,
1361 const cs_cell_sys_t *csys,
1362 cs_cdofb_adv_scheme_t *scheme_func,
1363 cs_cell_builder_t *cb)
1364 {
1365 /* Sanity checks */
1366 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOFB);
1367 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PF | CS_FLAG_COMP_PFQ));
1368
1369 /* Initialize the local matrix structure */
1370 cs_sdm_t *adv = cb->loc;
1371 cs_sdm_square_init(cm->n_fc + 1, adv);
1372
1373 if (cb->cell_flag & CS_FLAG_SOLID_CELL)
1374 return; /* Nothing to do. No advection in the current cell volume */
1375
1376 /* Define the local operator for advection. Boundary conditions are also
1377 treated here since there are always weakly enforced */
1378 scheme_func(eqp->dim, cm, csys, cb, adv);
1379
1380 /* Handle the specific case when there is no diffusion and no advection
1381 * flux. In this case, a zero row may appear leading to the divergence of
1382 * linear solver. To circumvent this issue, one set the boundary face value
1383 * to the cell value. This is equivalent to enforce a homogeneous Neumann
1384 * behavior. */
1385 assert(cs_equation_param_has_diffusion(eqp) == false);
1386
1387 cs_real_t max_abs_flux = 0.;
1388 for (int f = 0; f < cm->n_fc; f++)
1389 max_abs_flux = fmax(max_abs_flux, fabs(cb->adv_fluxes[f]));
1390
1391 const cs_real_t threshold = max_abs_flux * cs_math_epzero;
1392
1393 for (int f = 0; f < cm->n_fc; f++) {
1394
1395 if (fabs(cb->adv_fluxes[f]) < threshold) {
1396
1397 cs_real_t *f_row = adv->val + f*adv->n_rows;
1398
1399 f_row[cm->n_fc] += -1;
1400 f_row[f] += 1;
1401
1402 }
1403
1404 } /* Loop on cell faces */
1405
1406 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
1407 if (cs_dbg_cw_test(eqp, cm, csys)) {
1408 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection fluxes");
1409 cs_log_printf(CS_LOG_DEFAULT, "\n beta_fluxes>>");
1410 for (int f = 0; f < cm->n_fc; f++)
1411 cs_log_printf(CS_LOG_DEFAULT, "f%d;% -5.3e|",
1412 cm->f_ids[f], cb->adv_fluxes[f]);
1413 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
1414 cs_sdm_dump(cm->c_id, NULL, NULL, adv);
1415 }
1416 #endif
1417 }
1418
1419 /*----------------------------------------------------------------------------*/
1420 /*!
1421 * \brief Main function to build the cellwise advection operator for CDO
1422 * face-based schemes.
1423 * The local matrix related to this operator is stored in cb->loc
1424 *
1425 * One assumes that a diffusion term is present so that there is no
1426 * need to perform additional checkings on the well-posedness of the
1427 * operator.
1428 *
1429 * \param[in] eqp pointer to a cs_equation_param_t structure
1430 * \param[in] cm pointer to a cs_cell_mesh_t structure
1431 * \param[in] csys pointer to a cs_cell_sys_t structure
1432 * \param[in] scheme_func pointer to the function building the system
1433 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1434 */
1435 /*----------------------------------------------------------------------------*/
1436
1437 void
cs_cdofb_advection(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,cs_cdofb_adv_scheme_t * scheme_func,cs_cell_builder_t * cb)1438 cs_cdofb_advection(const cs_equation_param_t *eqp,
1439 const cs_cell_mesh_t *cm,
1440 const cs_cell_sys_t *csys,
1441 cs_cdofb_adv_scheme_t *scheme_func,
1442 cs_cell_builder_t *cb)
1443 {
1444 /* Sanity checks */
1445 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOFB);
1446 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PF | CS_FLAG_COMP_PFQ));
1447
1448 /* Initialize the local matrix structure */
1449 cs_sdm_t *adv = cb->loc;
1450 cs_sdm_square_init(cm->n_fc + 1, adv);
1451
1452 if (cb->cell_flag & CS_FLAG_SOLID_CELL)
1453 return; /* Nothing to do. No advection in the current cell volume */
1454
1455 /* Remark: The flux across the primal faces is stored in cb->adv_fluxes and
1456 should have been computed previously in a function compliant with the
1457 cs_cdofb_adv_open_hook_t prototype */
1458
1459 /* Define the local operator for advection. Boundary conditions are also
1460 treated here since there are always weakly enforced */
1461 scheme_func(eqp->dim, cm, csys, cb, adv);
1462
1463 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
1464 if (cs_dbg_cw_test(eqp, cm, csys)) {
1465 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection fluxes");
1466 cs_log_printf(CS_LOG_DEFAULT, "\n beta_fluxes>>");
1467 for (int f = 0; f < cm->n_fc; f++)
1468 cs_log_printf(CS_LOG_DEFAULT, "f%d;% -5.3e|",
1469 cm->f_ids[f], cb->adv_fluxes[f]);
1470 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
1471 cs_sdm_dump(cm->c_id, NULL, NULL, adv);
1472 }
1473 #endif
1474 }
1475
1476 /*----------------------------------------------------------------------------*/
1477 /*!
1478 * \brief Compute the convection operator attached to a cell with a CDO
1479 * face-based scheme
1480 * - non-conservative formulation beta \cdot \grad
1481 * - upwind scheme
1482 * Rely on the work performed during R. Milani's PhD
1483 *
1484 * A scalar-valued version is built. Only the enforcement of the
1485 * boundary condition depends on the variable dimension.
1486 * Remark: Usually the local matrix called hereafter adv is stored
1487 * in cb->loc
1488 *
1489 * \param[in] dim dimension of the variable (1 or 3)
1490 * \param[in] cm pointer to a cs_cell_mesh_t structure
1491 * \param[in] csys pointer to a cs_cell_sys_t structure
1492 * \param[in] cb pointer to a cs_cell_builder_t structure
1493 * \param[in, out] adv pointer to a local matrix to build
1494 */
1495 /*----------------------------------------------------------------------------*/
1496
1497 void
cs_cdofb_advection_upwnoc(int dim,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1498 cs_cdofb_advection_upwnoc(int dim,
1499 const cs_cell_mesh_t *cm,
1500 const cs_cell_sys_t *csys,
1501 cs_cell_builder_t *cb,
1502 cs_sdm_t *adv)
1503 {
1504 const cs_real_t *fluxes = cb->adv_fluxes;
1505
1506 /* Access the row containing current cell */
1507 const short int c = cm->n_fc; /* current cell's location in the matrix */
1508 double *c_row = adv->val + c*adv->n_rows;
1509
1510 if ((cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) && csys != NULL) {
1511
1512 /* There is at least one boundary face associated to this cell */
1513
1514 for (short int f = 0; f < cm->n_fc; f++) {
1515
1516 const cs_real_t beta_flx = cm->f_sgn[f]*fluxes[f];
1517 const cs_real_t beta_minus = 0.5*(fabs(beta_flx) - beta_flx);
1518 const cs_real_t beta_plus = 0.5*(fabs(beta_flx) + beta_flx);
1519
1520 /* access the row containing the current face */
1521 double *f_row = adv->val + f*adv->n_rows;
1522
1523 f_row[f] += beta_minus;
1524 f_row[c] -= beta_plus;
1525 c_row[f] -= beta_minus;
1526 c_row[c] += beta_minus;
1527
1528 if (csys->bf_ids[f] > -1) { /* This is a boundary face */
1529 if (csys->bf_flag[f] & CS_CDO_BC_DIRICHLET) {
1530
1531 /* Inward flux: add beta_minus = 0.5*(abs(flux) - flux) */
1532 f_row[f] += beta_minus;
1533
1534 /* Weak enforcement of the Dirichlet BCs.
1535 Update RHS for faces attached to a boundary face */
1536 for (int k = 0; k < dim; k++)
1537 csys->rhs[dim*f+k] += beta_minus * csys->dir_values[dim*f+k];
1538
1539 }
1540 }
1541
1542 } /* Loop on cell faces */
1543
1544 }
1545 else {
1546
1547 /* There is no boundary face associated to this cell */
1548
1549 for (short int f = 0; f < cm->n_fc; f++) {
1550
1551 const cs_real_t beta_flx = cm->f_sgn[f]*fluxes[f];
1552 const cs_real_t beta_minus = 0.5*(fabs(beta_flx) - beta_flx);
1553 const cs_real_t beta_plus = 0.5*(fabs(beta_flx) + beta_flx);
1554
1555 /* access the row containing the current face */
1556 double *f_row = adv->val + f*adv->n_rows;
1557
1558 f_row[f] += beta_minus;
1559 f_row[c] -= beta_plus;
1560 c_row[f] -= beta_minus;
1561 c_row[c] += beta_minus;
1562
1563 } /* Loop on cell faces */
1564
1565 }
1566
1567 }
1568
1569 /*----------------------------------------------------------------------------*/
1570 /*!
1571 * \brief Compute the convection operator attached to a cell with a CDO
1572 * face-based scheme
1573 * - conservative formulation div(\beta )
1574 * - upwind scheme
1575 * Rely on the work performed during R. Milani's PhD
1576 *
1577 * A scalar-valued version is built. Only the enforcement of the
1578 * boundary condition depends on the variable dimension.
1579 * Remark: Usually the local matrix called hereafter adv is stored
1580 * in cb->loc
1581 *
1582 * \param[in] dim dimension of the variable (1 or 3)
1583 * \param[in] cm pointer to a cs_cell_mesh_t structure
1584 * \param[in] csys pointer to a cs_cell_sys_t structure
1585 * \param[in] cb pointer to a cs_cell_builder_t structure
1586 * \param[in, out] adv pointer to a local matrix to build
1587 */
1588 /*----------------------------------------------------------------------------*/
1589
1590 void
cs_cdofb_advection_upwcsv(int dim,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1591 cs_cdofb_advection_upwcsv(int dim,
1592 const cs_cell_mesh_t *cm,
1593 const cs_cell_sys_t *csys,
1594 cs_cell_builder_t *cb,
1595 cs_sdm_t *adv)
1596 {
1597 const cs_real_t *fluxes = cb->adv_fluxes;
1598
1599 /* Access the row containing current cell */
1600 const short int c = cm->n_fc; /* current cell's location in the matrix */
1601 double *c_row = adv->val + c*adv->n_rows;
1602
1603 if ((cb->cell_flag & CS_FLAG_BOUNDARY_CELL_BY_FACE) && csys != NULL) {
1604
1605 /* There is at least one boundary face associated to this cell */
1606
1607 for (short int f = 0; f < cm->n_fc; f++) {
1608
1609 const cs_real_t beta_flx = cm->f_sgn[f]*fluxes[f];
1610 const cs_real_t beta_minus = 0.5*(fabs(beta_flx) - beta_flx);
1611 const cs_real_t beta_plus = 0.5*(fabs(beta_flx) + beta_flx);
1612
1613 /* access the row containing the current face */
1614 double *f_row = adv->val + f*adv->n_rows;
1615
1616 f_row[f] += beta_minus;
1617 f_row[c] -= beta_plus;
1618 c_row[f] -= beta_minus;
1619 c_row[c] += beta_plus;
1620
1621 if (csys->bf_ids[f] > -1) { /* This is a boundary face */
1622 if (csys->bf_flag[f] & CS_CDO_BC_DIRICHLET) {
1623
1624 /* Inward flux: add beta_minus = 0.5*(abs(flux) - flux) */
1625 f_row[f] += beta_minus;
1626
1627 /* Weak enforcement of the Dirichlet BCs.
1628 Update RHS for faces attached to a boundary face */
1629 for (int k = 0; k < dim; k++)
1630 csys->rhs[dim*f+k] += beta_minus * csys->dir_values[dim*f+k];
1631
1632 }
1633 }
1634
1635 } /* Loop on cell faces */
1636
1637 }
1638 else {
1639
1640 /* There is no boundary face associated to this cell */
1641
1642 for (short int f = 0; f < cm->n_fc; f++) {
1643
1644 const cs_real_t beta_flx = cm->f_sgn[f]*fluxes[f];
1645 const cs_real_t beta_minus = 0.5*(fabs(beta_flx) - beta_flx);
1646 const cs_real_t beta_plus = 0.5*(fabs(beta_flx) + beta_flx);
1647
1648 /* access the row containing the current face */
1649 double *f_row = adv->val + f*adv->n_rows;
1650
1651 f_row[f] += beta_minus;
1652 f_row[c] -= beta_plus;
1653 c_row[f] -= beta_minus;
1654 c_row[c] += beta_plus;
1655
1656 } /* Loop on cell faces */
1657
1658 }
1659
1660 }
1661
1662 /*----------------------------------------------------------------------------*/
1663 /*!
1664 * \brief Compute the convection operator attached to a cell with a CDO
1665 * face-based scheme
1666 * - non-conservative formulation beta.grad
1667 * - centered scheme
1668 * Rely on the work performed during R. Milani's PhD
1669 *
1670 * A scalar-valued version is built. Only the enforcement of the
1671 * boundary condition depends on the variable dimension.
1672 * Remark: Usually the local matrix called hereafter adv is stored
1673 * in cb->loc
1674 *
1675 * \param[in] dim dimension of the variable (1 or 3)
1676 * \param[in] cm pointer to a cs_cell_mesh_t structure
1677 * \param[in] csys pointer to a cs_cell_sys_t structure
1678 * \param[in] cb pointer to a cs_cell_builder_t structure
1679 * \param[in, out] adv pointer to a local matrix to build
1680 */
1681 /*----------------------------------------------------------------------------*/
1682
1683 void
cs_cdofb_advection_cennoc(int dim,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1684 cs_cdofb_advection_cennoc(int dim,
1685 const cs_cell_mesh_t *cm,
1686 const cs_cell_sys_t *csys,
1687 cs_cell_builder_t *cb,
1688 cs_sdm_t *adv)
1689 {
1690 const short int c = cm->n_fc; /* current cell's location in the matrix */
1691 const cs_real_t *fluxes = cb->adv_fluxes;
1692
1693 /* Access the row containing current cell */
1694 double *c_row = adv->val + c*adv->n_rows;
1695
1696 /* Loop on cell faces */
1697 for (short int f = 0; f < cm->n_fc; f++) {
1698
1699 const cs_real_t half_beta_flx = 0.5*cm->f_sgn[f]*fluxes[f];
1700
1701 /* access the row containing the current face */
1702 double *f_row = adv->val + f*adv->n_rows;
1703
1704 f_row[c] -= half_beta_flx;
1705 f_row[f] += half_beta_flx; /* Could be avoided:
1706 \sum_c(f) u_f v_f (\beta \cdot \nu_fc) = 0 */
1707 c_row[f] += half_beta_flx;
1708 c_row[c] -= half_beta_flx;
1709
1710 /* Apply boundary conditions */
1711 if (csys->bf_flag[f] & CS_CDO_BC_DIRICHLET) {
1712
1713 const cs_real_t beta_minus = 0.5*fabs(fluxes[f]) - half_beta_flx;
1714
1715 /* Inward flux: add beta_minus = 0.5*(abs(flux) - flux) */
1716 f_row[f] += beta_minus;
1717
1718 /* Weak enforcement of the Dirichlet BCs. Update RHS for faces attached to
1719 a boundary face */
1720 for (int k = 0; k < dim; k++)
1721 csys->rhs[dim*f+k] += beta_minus * csys->dir_values[dim*f+k];
1722
1723 }
1724
1725 } /* Loop on cell faces */
1726
1727 }
1728
1729 /*----------------------------------------------------------------------------*/
1730 /*!
1731 * \brief Compute the convection operator attached to a cell with a CDO
1732 * face-based scheme
1733 * - conservative formulation div(beta )
1734 * - centered scheme
1735 * Rely on the work performed during R. Milani's PhD
1736 *
1737 * A scalar-valued version is built. Only the enforcement of the
1738 * boundary condition depends on the variable dimension.
1739 * Remark: Usually the local matrix called hereafter adv is stored
1740 * in cb->loc
1741 *
1742 * \param[in] dim dimension of the variable (1 or 3)
1743 * \param[in] cm pointer to a cs_cell_mesh_t structure
1744 * \param[in] csys pointer to a cs_cell_sys_t structure
1745 * \param[in] cb pointer to a cs_cell_builder_t structure
1746 * \param[in, out] adv pointer to a local matrix to build
1747 */
1748 /*----------------------------------------------------------------------------*/
1749
1750 void
cs_cdofb_advection_cencsv(int dim,const cs_cell_mesh_t * cm,const cs_cell_sys_t * csys,cs_cell_builder_t * cb,cs_sdm_t * adv)1751 cs_cdofb_advection_cencsv(int dim,
1752 const cs_cell_mesh_t *cm,
1753 const cs_cell_sys_t *csys,
1754 cs_cell_builder_t *cb,
1755 cs_sdm_t *adv)
1756 {
1757 const short int c = cm->n_fc; /* current cell's location in the matrix */
1758 const cs_real_t *fluxes = cb->adv_fluxes;
1759
1760 cs_real_t div_c = 0;
1761
1762 /* Access the row containing current cell */
1763 double *c_row = adv->val + c*adv->n_rows;
1764
1765 /* Loop on cell faces */
1766 for (short int f = 0; f < cm->n_fc; f++) {
1767
1768 const cs_real_t half_beta_flx = 0.5*cm->f_sgn[f]*fluxes[f];
1769
1770 div_c += half_beta_flx;
1771
1772 /* access the row containing the current face */
1773 double *f_row = adv->val + f*adv->n_rows;
1774
1775 f_row[c] -= half_beta_flx;
1776 f_row[f] += half_beta_flx; /* Could be avoided:
1777 \sum_c(f) u_f v_f (\beta \cdot \nu_fc) = 0 */
1778 c_row[f] += half_beta_flx;
1779 c_row[c] -= half_beta_flx;
1780
1781 /* Apply boundary conditions */
1782 if (csys->bf_flag[f] & CS_CDO_BC_DIRICHLET) {
1783
1784 const cs_real_t beta_minus = 0.5*fabs(fluxes[f]) - half_beta_flx;
1785
1786 /* Inward flux: add beta_minus = 0.5*(abs(flux) - flux) */
1787 f_row[f] += beta_minus;
1788
1789 /* Weak enforcement of the Dirichlet BCs. Update RHS for faces attached to
1790 a boundary face */
1791 for (int k = 0; k < dim; k++)
1792 csys->rhs[dim*f+k] += beta_minus * csys->dir_values[dim*f+k];
1793
1794 }
1795
1796 } /* Loop on cell faces */
1797
1798 c_row[c] += 2*div_c; /* since half_beta_flx has been considered */
1799 }
1800
1801 /*----------------------------------------------------------------------------*/
1802 /*!
1803 * \brief Compute the convection operator attached to a cell with a CDO
1804 * vertex-based scheme with an upwind scheme and a conservative
1805 * formulation. The portion of upwinding relies on an evaluation
1806 * of the weigth associated to the given property
1807 * The local matrix related to this operator is stored in cb->loc
1808 * Predefined prototype to match the function pointer
1809 * cs_cdovb_advection_t
1810 *
1811 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
1812 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1813 * \param[in] diff_pty pointer to a \ref cs_property_data_t structure
1814 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1815 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1816 */
1817 /*----------------------------------------------------------------------------*/
1818
1819 void
cs_cdo_advection_vb_upwcsv_wpty(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)1820 cs_cdo_advection_vb_upwcsv_wpty(const cs_equation_param_t *eqp,
1821 const cs_cell_mesh_t *cm,
1822 const cs_property_data_t *diff_pty,
1823 cs_face_mesh_t *fm,
1824 cs_cell_builder_t *cb)
1825 {
1826 CS_UNUSED(fm);
1827
1828 /* Sanity checks */
1829 assert(diff_pty != NULL);
1830 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVB);
1831 assert(cs_eflag_test(cm->flag,
1832 CS_FLAG_COMP_PV | CS_FLAG_COMP_EV | CS_FLAG_COMP_PEQ |
1833 CS_FLAG_COMP_DFQ));
1834
1835 const cs_param_advection_scheme_t adv_scheme = eqp->adv_scheme;
1836
1837 /* Initialize the local matrix structure */
1838 cs_sdm_t *adv = cb->loc;
1839 cs_sdm_square_init(cm->n_vc, adv);
1840
1841 /* Compute the flux across the dual face attached to each edge of the cell */
1842 cs_real_t *fluxes = cb->values; /* size n_ec */
1843 cs_advection_field_cw_dface_flux(cm, eqp->adv_field, cb->t_bc_eval, fluxes);
1844
1845 /* Compute the criterion attached to each edge of the cell which is used
1846 to evaluate how to upwind */
1847 cs_real_t *upwcoef = cb->values + cm->n_ec;
1848
1849 for (short int e = 0; e < cm->n_ec; e++) {
1850
1851 const cs_nvec3_t dfq = cm->dface[e];
1852
1853 cs_real_t pty_contrib;
1854 if (diff_pty->is_iso)
1855 pty_contrib = diff_pty->value;
1856
1857 else { /* Property is considered as tensor-valued */
1858
1859 cs_real_t matnu[3];
1860 cs_math_33_3_product(diff_pty->tensor, dfq.unitv, matnu);
1861 pty_contrib = _dp3(dfq.unitv, matnu);
1862
1863 }
1864
1865 /* Define a coefficient close to a Peclet number */
1866 const double mean_flux = fluxes[e]/dfq.meas;
1867
1868 if (pty_contrib > cs_math_zero_threshold)
1869 upwcoef[e] = cm->edge[e].meas * mean_flux / pty_contrib;
1870 else
1871 upwcoef[e] = mean_flux * cs_math_big_r; /* dominated by convection */
1872
1873 } /* Loop on cell edges */
1874
1875 /* Set the function to compute the weight of upwinding */
1876 _upwind_weight_t *get_weight = _assign_weight_func(adv_scheme);
1877
1878 /* Define the local operator for advection */
1879 _build_cell_vpfd_upw(cm, get_weight, fluxes, upwcoef, adv);
1880
1881 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
1882 if (cs_dbg_cw_test(eqp, cm, NULL)) {
1883 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
1884 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
1885 }
1886 #endif
1887 }
1888
1889 /*----------------------------------------------------------------------------*/
1890 /*!
1891 * \brief Compute the convection operator attached to a cell with a CDO
1892 * vertex-based scheme without diffusion and an upwind scheme and a
1893 * conservative formulation is used.
1894 * The local matrix related to this operator is stored in cb->loc
1895 * Predefined prototype to match the function pointer
1896 * cs_cdovb_advection_t
1897 *
1898 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
1899 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1900 * \param[in] diff_pty pointer to a \ref cs_property_data_t structure
1901 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1902 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1903 */
1904 /*----------------------------------------------------------------------------*/
1905
1906 void
cs_cdo_advection_vb_upwcsv(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)1907 cs_cdo_advection_vb_upwcsv(const cs_equation_param_t *eqp,
1908 const cs_cell_mesh_t *cm,
1909 const cs_property_data_t *diff_pty,
1910 cs_face_mesh_t *fm,
1911 cs_cell_builder_t *cb)
1912 {
1913 CS_UNUSED(fm);
1914 CS_UNUSED(diff_pty);
1915
1916 /* Sanity checks */
1917 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVB);
1918 assert(cs_eflag_test(cm->flag,
1919 CS_FLAG_COMP_PV | CS_FLAG_COMP_EV | CS_FLAG_COMP_DFQ));
1920
1921 const cs_param_advection_scheme_t adv_scheme = eqp->adv_scheme;
1922
1923 /* Initialize the local matrix structure */
1924 cs_sdm_t *adv = cb->loc;
1925 cs_sdm_square_init(cm->n_vc, adv);
1926
1927 /* Compute the flux across the dual face attached to each edge of the cell */
1928 cs_real_t *fluxes = cb->values; /* size n_ec */
1929 cs_advection_field_cw_dface_flux(cm, eqp->adv_field, cb->t_bc_eval, fluxes);
1930
1931 /* Compute the criterion attached to each edge of the cell which is used
1932 to evaluate how to upwind */
1933 cs_real_t *upwcoef = cb->values + cm->n_ec;
1934 for (short int e = 0; e < cm->n_ec; e++)
1935 upwcoef[e] = fluxes[e]/cm->dface[e].meas;
1936
1937 /* Set the function to compute the weight of upwinding */
1938 _upwind_weight_t *get_weight = _assign_weight_func(adv_scheme);
1939
1940 /* Define the local operator for advection */
1941 _build_cell_vpfd_upw(cm, get_weight, fluxes, upwcoef, adv);
1942
1943 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
1944 if (cs_dbg_cw_test(eqp, cm, NULL)) {
1945 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
1946 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
1947 }
1948 #endif
1949 }
1950
1951 /*----------------------------------------------------------------------------*/
1952 /*!
1953 * \brief Compute the convection operator attached to a cell with a CDO
1954 * vertex-based scheme when a centered scheme and a conservative
1955 * formulation is used.
1956 * The local matrix related to this operator is stored in cb->loc
1957 * Predefined prototype to match the function pointer
1958 * cs_cdovb_advection_t
1959 *
1960 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
1961 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1962 * \param[in] diff_pty pointer to a \ref cs_property_data_t structure
1963 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
1964 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
1965 */
1966 /*----------------------------------------------------------------------------*/
1967
1968 void
cs_cdo_advection_vb_cencsv(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)1969 cs_cdo_advection_vb_cencsv(const cs_equation_param_t *eqp,
1970 const cs_cell_mesh_t *cm,
1971 const cs_property_data_t *diff_pty,
1972 cs_face_mesh_t *fm,
1973 cs_cell_builder_t *cb)
1974 {
1975 CS_UNUSED(fm);
1976 CS_UNUSED(diff_pty);
1977
1978 /* Sanity check */
1979 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVB); /* Sanity check */
1980 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PV | CS_FLAG_COMP_EV));
1981
1982 /* Initialize the local matrix structure */
1983 cs_sdm_t *adv = cb->loc;
1984 cs_sdm_square_init(cm->n_vc, adv);
1985
1986 /* Compute the flux across the dual face attached to each edge of the cell */
1987 cs_real_t *fluxes = cb->values; /* size n_ec */
1988 cs_advection_field_cw_dface_flux(cm, eqp->adv_field, cb->t_bc_eval, fluxes);
1989
1990 /* Define the local operator for advection */
1991 _build_cell_vpfd_cen(cm, fluxes, adv);
1992
1993 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
1994 if (cs_dbg_cw_test(eqp, cm, NULL)) {
1995 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
1996 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
1997 }
1998 #endif
1999 }
2000
2001 /*----------------------------------------------------------------------------*/
2002 /*!
2003 * \brief Compute the convection operator attached to a cell with a CDO
2004 * vertex-based scheme when a mixed centered/upwind scheme with
2005 * a conservative formulation is used.
2006 * The local matrix related to this operator is stored in cb->loc
2007 * Predefined prototype to match the function pointer
2008 * cs_cdovb_advection_t
2009 *
2010 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
2011 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2012 * \param[in] diff_pty pointer to a \ref cs_property_data_t structure
2013 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2014 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2015 */
2016 /*----------------------------------------------------------------------------*/
2017
2018 void
cs_cdo_advection_vb_mcucsv(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)2019 cs_cdo_advection_vb_mcucsv(const cs_equation_param_t *eqp,
2020 const cs_cell_mesh_t *cm,
2021 const cs_property_data_t *diff_pty,
2022 cs_face_mesh_t *fm,
2023 cs_cell_builder_t *cb)
2024 {
2025 CS_UNUSED(fm);
2026 CS_UNUSED(diff_pty);
2027
2028 /* Sanity check */
2029 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVB); /* Sanity check */
2030 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PV | CS_FLAG_COMP_EV));
2031
2032 /* Initialize the local matrix structure */
2033 cs_sdm_t *adv = cb->loc;
2034 cs_sdm_square_init(cm->n_vc, adv);
2035
2036 /* Compute the flux across the dual face attached to each edge of the cell */
2037 cs_real_t *fluxes = cb->values; /* size n_ec */
2038 cs_advection_field_cw_dface_flux(cm, eqp->adv_field, cb->t_bc_eval, fluxes);
2039
2040 /* Define the local operator for advection */
2041 _build_cell_vpfd_mcu(cm, eqp->upwind_portion, fluxes, adv);
2042
2043 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
2044 if (cs_dbg_cw_test(eqp, cm, NULL)) {
2045 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
2046 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
2047 }
2048 #endif
2049 }
2050
2051 /*----------------------------------------------------------------------------*/
2052 /*!
2053 * \brief Compute the convection operator attached to a cell with a CDO
2054 * vertex-based scheme with an upwind scheme and a conservative
2055 * formulation. The portion of upwinding relies on an evaluation
2056 * of the weigth associated to the given property.
2057 * The local matrix related to this operator is stored in cb->loc
2058 * Predefined prototype to match the function pointer
2059 * cs_cdovb_advection_t
2060 *
2061 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
2062 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2063 * \param[in] diff_pty pointer to a \ref cs_property_data_t structure
2064 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2065 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2066 */
2067 /*----------------------------------------------------------------------------*/
2068
2069 void
cs_cdo_advection_vb_upwnoc_wpty(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)2070 cs_cdo_advection_vb_upwnoc_wpty(const cs_equation_param_t *eqp,
2071 const cs_cell_mesh_t *cm,
2072 const cs_property_data_t *diff_pty,
2073 cs_face_mesh_t *fm,
2074 cs_cell_builder_t *cb)
2075 {
2076 CS_UNUSED(fm);
2077
2078 /* Sanity checks */
2079 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVB);
2080 assert(cs_eflag_test(cm->flag,
2081 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ));
2082
2083 const cs_param_advection_scheme_t adv_scheme = eqp->adv_scheme;
2084
2085 /* Initialize the local matrix structure */
2086 cs_sdm_t *adv = cb->loc;
2087 cs_sdm_square_init(cm->n_vc, adv);
2088
2089 /* Compute the flux across the dual face attached to each edge of the cell */
2090 cs_real_t *fluxes = cb->values; /* size n_ec */
2091 cs_advection_field_cw_dface_flux(cm, eqp->adv_field, cb->t_bc_eval, fluxes);
2092
2093 /* Compute the criterion attached to each edge of the cell which is used
2094 to evaluate how to upwind */
2095 cs_real_t *upwcoef = cb->values + cm->n_ec;
2096
2097 for (short int e = 0; e < cm->n_ec; e++) {
2098
2099 const cs_nvec3_t dfq = cm->dface[e];
2100
2101 cs_real_t pty_contrib;
2102 if (diff_pty->is_iso)
2103 pty_contrib = diff_pty->value;
2104
2105 else { /* Property is considered as tensor-valued */
2106
2107 cs_real_t matnu[3];
2108 cs_math_33_3_product(diff_pty->tensor, dfq.unitv, matnu);
2109 pty_contrib = _dp3(dfq.unitv, matnu);
2110
2111 }
2112
2113 /* Define a coefficient close to a Peclet number */
2114 const double mean_flux = fluxes[e]/dfq.meas;
2115
2116 if (pty_contrib > cs_math_zero_threshold)
2117 upwcoef[e] = cm->edge[e].meas * mean_flux / pty_contrib;
2118 else
2119 upwcoef[e] = mean_flux * cs_math_big_r; /* dominated by convection */
2120
2121 } /* Loop on cell edges */
2122
2123 /* Set the function to compute the weight of upwinding */
2124 _upwind_weight_t *get_weight = _assign_weight_func(adv_scheme);
2125
2126 /* Define the local operator for advection */
2127 _build_cell_epcd_upw(cm, get_weight, fluxes, upwcoef, adv);
2128
2129 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
2130 if (cs_dbg_cw_test(eqp, cm, NULL)) {
2131 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
2132 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
2133 }
2134 #endif
2135 }
2136
2137 /*----------------------------------------------------------------------------*/
2138 /*!
2139 * \brief Compute the convection operator attached to a cell with a CDO
2140 * vertex-based scheme without diffusion when an upwind scheme and a
2141 * conservative formulation is used.
2142 * The local matrix related to this operator is stored in cb->loc
2143 * Predefined prototype to match the function pointer
2144 * cs_cdovb_advection_t
2145 *
2146 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
2147 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2148 * \param[in] diff_pty pointer to a \ref cs_property_data_t structure
2149 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2150 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2151 */
2152 /*----------------------------------------------------------------------------*/
2153
2154 void
cs_cdo_advection_vb_upwnoc(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)2155 cs_cdo_advection_vb_upwnoc(const cs_equation_param_t *eqp,
2156 const cs_cell_mesh_t *cm,
2157 const cs_property_data_t *diff_pty,
2158 cs_face_mesh_t *fm,
2159 cs_cell_builder_t *cb)
2160 {
2161 CS_UNUSED(fm);
2162 CS_UNUSED(diff_pty);
2163
2164 /* Sanity checks */
2165 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVB); /* Sanity check */
2166 assert(cs_eflag_test(cm->flag,
2167 CS_FLAG_COMP_PV | CS_FLAG_COMP_DFQ | CS_FLAG_COMP_EV));
2168
2169 const cs_param_advection_scheme_t adv_scheme = eqp->adv_scheme;
2170
2171 /* Initialize the local matrix structure */
2172 cs_sdm_t *adv = cb->loc;
2173 cs_sdm_square_init(cm->n_vc, adv);
2174
2175 /* Compute the flux across the dual face attached to each edge of the cell */
2176 cs_real_t *fluxes = cb->values; /* size n_ec */
2177 cs_advection_field_cw_dface_flux(cm, eqp->adv_field, cb->t_bc_eval, fluxes);
2178
2179 /* Compute the criterion attached to each edge of the cell which is used
2180 to evaluate how to upwind */
2181 cs_real_t *upwcoef = cb->values + cm->n_ec;
2182 for (short int e = 0; e < cm->n_ec; e++)
2183 upwcoef[e] = fluxes[e]/cm->dface[e].meas;
2184
2185 /* Set the function to compute the weight of upwinding */
2186 _upwind_weight_t *get_weight = _assign_weight_func(adv_scheme);
2187
2188 /* Define the local operator for advection */
2189 _build_cell_epcd_upw(cm, get_weight, fluxes, upwcoef, adv);
2190
2191 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
2192 if (cs_dbg_cw_test(eqp, cm, NULL)) {
2193 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
2194 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
2195 }
2196 #endif
2197 }
2198
2199 /*----------------------------------------------------------------------------*/
2200 /*!
2201 * \brief Compute the convection operator attached to a cell with a CDO
2202 * vertex-based scheme when a centered scheme and a non-conservative
2203 * formulation is used.
2204 * The local matrix related to this operator is stored in cb->loc
2205 * Predefined prototype to match the function pointer
2206 * cs_cdovb_advection_t
2207 *
2208 * \param[in] eqp pointer to a \ref cs_equation_param_t structure
2209 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
2210 * \param[in] diff_pty pointer to a \ref cs_property_data_t structure
2211 * \param[in, out] fm pointer to a \ref cs_face_mesh_t structure
2212 * \param[in, out] cb pointer to a \ref cs_cell_builder_t structure
2213 */
2214 /*----------------------------------------------------------------------------*/
2215
2216 void
cs_cdo_advection_vb_cennoc(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)2217 cs_cdo_advection_vb_cennoc(const cs_equation_param_t *eqp,
2218 const cs_cell_mesh_t *cm,
2219 const cs_property_data_t *diff_pty,
2220 cs_face_mesh_t *fm,
2221 cs_cell_builder_t *cb)
2222 {
2223 CS_UNUSED(fm);
2224 CS_UNUSED(diff_pty);
2225
2226 /* Sanity checks */
2227 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVB);
2228 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PV | CS_FLAG_COMP_EV));
2229
2230 /* Initialize the local matrix structure */
2231 cs_sdm_t *adv = cb->loc;
2232 cs_sdm_square_init(cm->n_vc, adv);
2233
2234 /* Compute the flux across the dual face attached to each edge of the cell */
2235 cs_real_t *fluxes = cb->values; /* size n_ec */
2236 cs_advection_field_cw_dface_flux(cm, eqp->adv_field, cb->t_bc_eval, fluxes);
2237
2238 /* Define the local operator for advection */
2239 _build_cell_epcd_cen(cm, fluxes, adv);
2240
2241 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
2242 if (cs_dbg_cw_test(eqp, cm, NULL)) {
2243 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
2244 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
2245 }
2246 #endif
2247 }
2248
2249 /*----------------------------------------------------------------------------*/
2250 /*!
2251 * \brief Compute the convection operator attached to a cell with a CDO
2252 * vertex+cell-based scheme when the advection field is cellwise
2253 * constant
2254 *
2255 * \param[in] eqp pointer to a cs_equation_param_t structure
2256 * \param[in] cm pointer to a cs_cell_mesh_t structure
2257 * \param[in] diff_pty pointer to the property associated to diffusion
2258 * \param[in, out] fm pointer to a cs_face_mesh_t structure
2259 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2260 */
2261 /*----------------------------------------------------------------------------*/
2262
2263 void
cs_cdo_advection_vcb_cw_cst(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)2264 cs_cdo_advection_vcb_cw_cst(const cs_equation_param_t *eqp,
2265 const cs_cell_mesh_t *cm,
2266 const cs_property_data_t *diff_pty,
2267 cs_face_mesh_t *fm,
2268 cs_cell_builder_t *cb)
2269 {
2270 CS_UNUSED(diff_pty);
2271
2272 /* Sanity checks */
2273 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVCB);
2274 assert(eqp->adv_formulation == CS_PARAM_ADVECTION_FORM_NONCONS);
2275 assert(eqp->adv_scheme == CS_PARAM_ADVECTION_SCHEME_CIP_CW);
2276 assert(cs_eflag_test(cm->flag,
2277 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ |
2278 CS_FLAG_COMP_DEQ | CS_FLAG_COMP_EF | CS_FLAG_COMP_FEQ |
2279 CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ));
2280
2281 const int n_sysc = cm->n_vc + 1;
2282
2283 /* Initialize local matrix structure */
2284 cs_sdm_t *a = cb->loc;
2285 cs_sdm_square_init(n_sysc, a);
2286
2287 /* Use a cellwise constant approximation of the advection field */
2288 cs_nvec3_t adv_cell;
2289 cs_advection_field_get_cell_vector(cm->c_id, eqp->adv_field, &adv_cell);
2290
2291 if (adv_cell.meas < cs_math_get_machine_epsilon())
2292 return;
2293
2294 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 2
2295 const int n_vc = cm->n_vc;
2296 cs_sdm_t *cons = cs_sdm_square_create(n_sysc);
2297
2298 cons->n_rows = n_sysc;
2299 for (short int v = 0; v < n_sysc*n_sysc; v++) cons->val[v] = 0;
2300 #endif
2301
2302 /* Stabilization coefficent * |beta_c| */
2303 const double stab_coef = cs_cip_stab_coef * adv_cell.meas;
2304
2305 /* Temporary buffers
2306 bgc stored in cb->values (size n_fc)
2307 tef stored in cb->values + cm->n_fc (size 2*n_ec)
2308
2309 bgvf stored in cb->vectors (size: 2*n_ec)
2310 */
2311 cs_sdm_t *af = cb->aux;
2312 double *tef_save = cb->values + cm->n_fc; /* size = 2*n_ec */
2313
2314 for (short int f = 0; f < cm->n_fc; f++) { /* Loop on cell faces */
2315
2316 /* Build a facewise view of the mesh */
2317 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2318
2319 const int n_sysf = fm->n_vf + 1;
2320
2321 /* Initialize the local face matrix */
2322 af->n_rows = n_sysf;
2323 for (short int i = 0; i < n_sysf*n_sysf; i++) af->val[i] = 0.;
2324
2325 /* Store tef areas for a future usage (Second part of the stabilization) */
2326 const short int fshift = cm->f2e_idx[f];
2327 double *tef = tef_save + fshift;
2328 for (short int e = 0; e < fm->n_ef; e++) tef[e] = fm->tef[e];
2329
2330 /* Update the face matrix inside with the consistent */
2331 _vcb_cellwise_consistent_part(adv_cell, cm, fm, cb);
2332
2333 /* Build the first part of the stabilization. (inside pfc) */
2334 _vcb_stabilization_part1(cm, fm, stab_coef, cb);
2335
2336 /* Reorder v1 and v2 to insure coherency for edges shared between
2337 two faces */
2338 cs_real_3_t *bgvf = cb->vectors + fshift;
2339 for (short int e = 0; e < fm->n_ef; e++) {
2340 if (fm->v_ids[fm->e2v_ids[2*e]] > fm->v_ids[fm->e2v_ids[2*e+1]]) {
2341 double save = bgvf[e][0];
2342 bgvf[e][0] = bgvf[e][1];
2343 bgvf[e][1] = save;
2344 }
2345 }
2346
2347 /* Add the face matrix to the cell matrix */
2348 for (short int vi = 0; vi < fm->n_vf; vi++) {
2349
2350 double *aci = a->val + n_sysc*fm->v_ids[vi];
2351 const double *afi = af->val + n_sysf*vi;
2352 for (short int vj = 0; vj < fm->n_vf; vj++)
2353 aci[fm->v_ids[vj]] += afi[vj]; /* (i,j) face --> cell */
2354 aci[cm->n_vc] += afi[fm->n_vf]; /* (i,c) face --> cell */
2355
2356 }
2357
2358 double *acc = a->val + n_sysc*cm->n_vc;
2359 const double *afc = af->val + n_sysf*fm->n_vf;
2360 for (short int vj = 0; vj < fm->n_vf; vj++)
2361 acc[fm->v_ids[vj]] += afc[vj]; /* (c,j) face --> cell */
2362 acc[cm->n_vc] += afc[fm->n_vf];
2363
2364 } /* Loop on cell faces */
2365
2366 /* Build the second part of the stabilization. */
2367 _vcb_stabilization_part2(cm, stab_coef, cb);
2368
2369 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
2370 if (cs_dbg_cw_test(eqp, cm, NULL)) {
2371 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix (CW version)");
2372 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
2373 }
2374 #if CS_CDO_ADVECTION_DBG > 2
2375 cons = cs_sdm_free(cons);
2376 #endif
2377 #endif
2378 }
2379
2380 /*----------------------------------------------------------------------------*/
2381 /*!
2382 * \brief Compute the convection operator attached to a cell with a CDO
2383 * vertex+cell-based scheme
2384 *
2385 * \param[in] eqp pointer to a cs_equation_param_t structure
2386 * \param[in] cm pointer to a cs_cell_mesh_t structure
2387 * \param[in] diff_pty pointer to the property associated to diffusion
2388 * \param[in, out] fm pointer to a cs_face_mesh_t structure
2389 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2390 */
2391 /*----------------------------------------------------------------------------*/
2392
2393 void
cs_cdo_advection_vcb(const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,const cs_property_data_t * diff_pty,cs_face_mesh_t * fm,cs_cell_builder_t * cb)2394 cs_cdo_advection_vcb(const cs_equation_param_t *eqp,
2395 const cs_cell_mesh_t *cm,
2396 const cs_property_data_t *diff_pty,
2397 cs_face_mesh_t *fm,
2398 cs_cell_builder_t *cb)
2399 {
2400 CS_UNUSED(diff_pty);
2401
2402 /* Sanity checks */
2403 assert(eqp->space_scheme == CS_SPACE_SCHEME_CDOVCB);
2404 assert(eqp->adv_formulation == CS_PARAM_ADVECTION_FORM_NONCONS);
2405 assert(eqp->adv_scheme == CS_PARAM_ADVECTION_SCHEME_CIP);
2406 assert(cs_eflag_test(cm->flag,
2407 CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ |
2408 CS_FLAG_COMP_DEQ | CS_FLAG_COMP_EF | CS_FLAG_COMP_FEQ |
2409 CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ));
2410
2411 const int n_sysc = cm->n_vc + 1;
2412
2413 /* Initialize local matrix structure */
2414 cs_sdm_t *a = cb->loc;
2415 cs_sdm_square_init(n_sysc, a);
2416
2417 /* Use a cellwise constant approximation of the advection field */
2418 cs_nvec3_t adv_cell;
2419 cs_advection_field_get_cell_vector(cm->c_id, eqp->adv_field, &adv_cell);
2420
2421 if (adv_cell.meas < cs_math_get_machine_epsilon())
2422 return;
2423
2424 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 2
2425 const int n_vc = cm->n_vc;
2426 cs_sdm_t *cons = cs_sdm_square_create(n_sysc);
2427
2428 cons->n_rows = n_sysc;
2429 for (short int v = 0; v < n_sysc*n_sysc; v++) cons->val[v] = 0;
2430 #endif
2431
2432 /* Stabilization coefficent * |beta_c| */
2433 const double stab_coef = cs_cip_stab_coef * adv_cell.meas;
2434
2435 /* Temporary buffers:
2436 bgc stored in cb->values (size n_fc)
2437 tef stored in cb->values + cm->n_fc (size 2*n_ec)
2438
2439 bgvf stored in cb->vectors (size: 2*n_ec)
2440 */
2441 cs_sdm_t *af = cb->aux;
2442 double *tef_save = cb->values + cm->n_fc; /* size = 2*n_ec */
2443
2444 for (short int f = 0; f < cm->n_fc; f++) { /* Loop on cell faces */
2445
2446 /* Build a facewise view of the mesh */
2447 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2448
2449 const int n_sysf = fm->n_vf + 1;
2450
2451 /* Initialize the local face matrix */
2452 cs_sdm_square_init(n_sysf, af);
2453
2454 /* Store tef areas for a future usage (Second part of the stabilization) */
2455 const short int fshift = cm->f2e_idx[f];
2456 double *tef = tef_save + fshift;
2457 for (short int e = 0; e < fm->n_ef; e++) tef[e] = fm->tef[e];
2458
2459 /* Initialize and update the face matrix inside (build bgvf inside) */
2460 _vcb_consistent_part(eqp->adv_field, adv_cell, cm, fm, cb->t_bc_eval, cb);
2461
2462 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 2
2463 _assemble_face(cm, fm, af, cons);
2464 #endif
2465
2466 /* Build the first part of the stabilization: that inside pfc */
2467 _vcb_stabilization_part1(cm, fm, stab_coef, cb);
2468
2469 /* Reorder v1 and v2 to insure coherency for edges shared between
2470 two faces */
2471 cs_real_3_t *bgvf = cb->vectors + fshift;
2472 for (short int e = 0; e < fm->n_ef; e++) {
2473 if (fm->v_ids[fm->e2v_ids[2*e]] > fm->v_ids[fm->e2v_ids[2*e+1]]) {
2474 double save = bgvf[e][0];
2475 bgvf[e][0] = bgvf[e][1];
2476 bgvf[e][1] = save;
2477 }
2478 }
2479
2480 /* Add the face matrix to the cell matrix */
2481 _assemble_face(cm, fm, af, a);
2482
2483 } /* Loop on cell faces */
2484
2485 /* Build the second part of the stabilization. */
2486 _vcb_stabilization_part2(cm, stab_coef, cb);
2487
2488 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 0
2489 if (cs_dbg_cw_test(eqp, cm, NULL)) {
2490 cs_log_printf(CS_LOG_DEFAULT, "\n>> Cell advection matrix");
2491 cs_sdm_dump(cm->c_id, NULL, NULL, cb->loc);
2492 #if CS_CDO_ADVECTION_DBG > 2
2493 cs_log_printf(CS_LOG_DEFAULT, "\n>> Advection matrix (CONSISTENCY PART)");
2494 cs_sdm_dump(cm->c_id, NULL, NULL, cons);
2495 cons = cs_sdm_free(cons);
2496 #endif
2497 }
2498 #endif
2499 }
2500
2501 /*----------------------------------------------------------------------------*/
2502 /*!
2503 * \brief Compute the BC contribution for the convection operator
2504 *
2505 * \param[in] cm pointer to a cs_cell_mesh_t structure
2506 * \param[in] eqp pointer to a cs_equation_param_t structure
2507 * \param[in] t_eval time at which one evaluates the advection field
2508 * \param[in, out] fm pointer to a cs_face_mesh_t structure
2509 * \param[in, out] cb pointer to a convection builder structure
2510 * \param[in, out] csys cell-wise structure storing the local system
2511 */
2512 /*----------------------------------------------------------------------------*/
2513
2514 void
cs_cdo_advection_vb_bc(const cs_cell_mesh_t * cm,const cs_equation_param_t * eqp,cs_real_t t_eval,cs_face_mesh_t * fm,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2515 cs_cdo_advection_vb_bc(const cs_cell_mesh_t *cm,
2516 const cs_equation_param_t *eqp,
2517 cs_real_t t_eval,
2518 cs_face_mesh_t *fm,
2519 cs_cell_builder_t *cb,
2520 cs_cell_sys_t *csys)
2521 {
2522 CS_UNUSED(fm);
2523
2524 /* Sanity checks */
2525 assert(cs_eflag_test(cm->flag,
2526 CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FEQ |
2527 CS_FLAG_COMP_EV | CS_FLAG_COMP_FV));
2528
2529 cs_real_t *tmp_rhs = cb->values;
2530 cs_real_t *mat_diag = cb->values + cm->n_vc;
2531 cs_real_t *v_nflx = cb->values + 2*cm->n_vc;
2532
2533 const cs_adv_field_t *adv = eqp->adv_field;
2534
2535 /* Reset local temporary RHS and diagonal contributions */
2536 for (short int v = 0; v < cm->n_vc; v++) mat_diag[v] = tmp_rhs[v] = 0;
2537
2538 cs_real_t scaling = 1; /* By default no scaling */
2539 if (eqp->adv_scaling_property != NULL) {
2540 if (cs_property_is_uniform(eqp->adv_scaling_property))
2541 scaling = eqp->adv_scaling_property->ref_value;
2542 else
2543 scaling = cs_property_value_in_cell(cm,
2544 eqp->adv_scaling_property,
2545 cb->t_pty_eval);
2546 }
2547
2548 /* Add diagonal term for vertices attached to a boundary face where
2549 the advection field points inward. */
2550 for (short int i = 0; i < csys->n_bc_faces; i++) { /* Loop on border faces */
2551
2552 /* Get the boundary face in the cell numbering */
2553 const short int f = csys->_f_ids[i];
2554
2555 cs_advection_field_cw_boundary_f2v_flux(cm, adv, f, t_eval, v_nflx);
2556
2557 if (eqp->adv_scaling_property != NULL)
2558 for (short int v = 0; v < cm->n_vc; v++) v_nflx[v] *= scaling;
2559
2560 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 1
2561 if (cs_dbg_cw_test(eqp, cm, csys)) {
2562 cs_log_printf(CS_LOG_DEFAULT, " %s: f:%d --> bndy_flux:", __func__, f);
2563 for (short int v = 0; v < cm->n_vc; v++)
2564 cs_log_printf(CS_LOG_DEFAULT, " v%d:%e", v, v_nflx[v]);
2565 cs_log_printf(CS_LOG_DEFAULT, "\n");
2566 }
2567 #endif
2568
2569 if (eqp->adv_formulation == CS_PARAM_ADVECTION_FORM_CONSERV) {
2570
2571 for (int iv = cm->f2v_idx[f]; iv < cm->f2v_idx[f+1]; iv++) {
2572
2573 const short int v_id = cm->f2v_ids[iv];
2574
2575 if (v_nflx[v_id] < 0) {
2576 /* advection field is inward w.r.t. the face normal */
2577 if (csys->bf_flag[f] & CS_CDO_BC_DIRICHLET)
2578 /* Homogeneous Dirichlet don't contribute. Other Bcs are invalid */
2579 tmp_rhs[v_id] -= v_nflx[v_id] * csys->dir_values[v_id];
2580 }
2581 else /* advection is oriented outward */
2582 mat_diag[v_id] += v_nflx[v_id];
2583
2584 } /* Loop on face vertices */
2585
2586 }
2587 else { /* Non-conservative formulation */
2588
2589 for (int iv = cm->f2v_idx[f]; iv < cm->f2v_idx[f+1]; iv++) {
2590
2591 const short int v_id = cm->f2v_ids[iv];
2592
2593 if (v_nflx[v_id] < 0) {
2594
2595 /* advection field is inward w.r.t. the face normal */
2596 if (csys->bf_flag[f] & CS_CDO_BC_DIRICHLET)
2597 /* Homogeneous Dirichlet don't contribute. Other Bcs are invalid */
2598 tmp_rhs[v_id] -= v_nflx[v_id] * csys->dir_values[v_id];
2599
2600 mat_diag[v_id] -= v_nflx[v_id];
2601
2602 }
2603
2604 } /* Loop on face vertices */
2605
2606 } /* Type of formulation */
2607
2608 } /* Loop on border faces */
2609
2610 /* Update the diagonal and the RHS of the local system matrix */
2611 for (short int v = 0; v < cm->n_vc; v++) {
2612 csys->mat->val[v*cm->n_vc + v] += mat_diag[v];
2613 csys->rhs[v] += tmp_rhs[v];
2614 }
2615
2616 }
2617
2618 /*----------------------------------------------------------------------------*/
2619 /*!
2620 * \brief Compute the BC contribution for the convection operator with CDO
2621 * V+C schemes
2622 *
2623 * \param[in] cm pointer to a cs_cell_mesh_t structure
2624 * \param[in] eqp pointer to a cs_equation_param_t structure
2625 * \param[in] t_eval time at which one evaluates the advection field
2626 * \param[in, out] fm pointer to a cs_face_mesh_t structure
2627 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2628 * \param[in, out] csys cell-wise structure storing the local system
2629 */
2630 /*----------------------------------------------------------------------------*/
2631
2632 void
cs_cdo_advection_vcb_bc(const cs_cell_mesh_t * cm,const cs_equation_param_t * eqp,cs_real_t t_eval,cs_face_mesh_t * fm,cs_cell_builder_t * cb,cs_cell_sys_t * csys)2633 cs_cdo_advection_vcb_bc(const cs_cell_mesh_t *cm,
2634 const cs_equation_param_t *eqp,
2635 cs_real_t t_eval,
2636 cs_face_mesh_t *fm,
2637 cs_cell_builder_t *cb,
2638 cs_cell_sys_t *csys)
2639 {
2640 if (csys == NULL)
2641 return;
2642
2643 /* Sanity checks */
2644 assert(cm != NULL && eqp != NULL);
2645 assert(csys->mat->n_rows == cm->n_vc + 1);
2646 assert(eqp->adv_formulation == CS_PARAM_ADVECTION_FORM_NONCONS);
2647 assert(eqp->adv_scheme == CS_PARAM_ADVECTION_SCHEME_CIP ||
2648 eqp->adv_scheme == CS_PARAM_ADVECTION_SCHEME_CIP_CW);
2649 assert(cs_eflag_test(cm->flag,
2650 CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
2651 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
2652
2653 const cs_adv_field_t *adv = eqp->adv_field;
2654
2655 /* Loop on border faces */
2656 for (short int i = 0; i < csys->n_bc_faces; i++) {
2657
2658 /* Get the boundary face in the cell numbering */
2659 const short int f = csys->_f_ids[i];
2660 const cs_real_t nflx = cs_advection_field_cw_boundary_face_flux(t_eval, f,
2661 cm, adv);
2662 const cs_real_t beta_nflx = 0.5 * (fabs(nflx) - nflx);
2663
2664 #if defined(DEBUG) && !defined(NDEBUG) && CS_CDO_ADVECTION_DBG > 1
2665 if (cs_dbg_cw_test(eqp, cm, csys))
2666 cs_log_printf(CS_LOG_DEFAULT, " %s: f:%d --> bndy_flux: %e\n",
2667 __func__, f, nflx);
2668 #endif
2669
2670 if (beta_nflx > 0) {
2671
2672 cs_face_mesh_build_from_cell_mesh(cm, f, fm);
2673
2674 cs_hodge_compute_wbs_surfacic(fm, cb->aux);
2675
2676 _update_vcb_system_with_bc(beta_nflx/fm->face.meas, fm, cb->aux, csys);
2677
2678 } /* At least for one face vertex, beta_nf > 0 */
2679
2680 } /* Loop on border faces */
2681
2682 }
2683
2684 /*----------------------------------------------------------------------------*/
2685 /*!
2686 * \brief Compute the value of the upwinding coefficient in each cell knowing
2687 * the related Peclet number
2688 *
2689 * \param[in] cdoq pointer to the cdo quantities structure
2690 * \param[in] scheme type of scheme used for the advection term
2691 * \param[in, out] coefval pointer to the pointer of real numbers to fill
2692 * in: Peclet number in each cell
2693 * out: value of the upwind coefficient
2694 */
2695 /*----------------------------------------------------------------------------*/
2696
2697 void
cs_cdo_advection_cell_upwind_coef(const cs_cdo_quantities_t * cdoq,cs_param_advection_scheme_t scheme,cs_real_t coefval[])2698 cs_cdo_advection_cell_upwind_coef(const cs_cdo_quantities_t *cdoq,
2699 cs_param_advection_scheme_t scheme,
2700 cs_real_t coefval[])
2701 {
2702 /* Sanity check */
2703 assert(coefval != NULL);
2704
2705 /* Set the function to compute the weight of upwinding */
2706 _upwind_weight_t *get_weight = _assign_weight_func(scheme);
2707
2708 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++)
2709 coefval[c_id] = get_weight(coefval[c_id]);
2710 }
2711
2712 /*----------------------------------------------------------------------------*/
2713
2714 #undef _dp3
2715
2716 END_C_DECLS
2717