1 /*============================================================================
2 * Manage the definition/setting of advection fields
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 <assert.h>
34 #include <ctype.h>
35 #include <float.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include <bft_mem.h>
45
46 #include "cs_boundary_zone.h"
47 #include "cs_evaluate.h"
48 #include "cs_equation_common.h"
49 #include "cs_field.h"
50 #include "cs_log.h"
51 #include "cs_math.h"
52 #include "cs_mesh_location.h"
53 #include "cs_param_cdo.h"
54 #include "cs_reco.h"
55 #include "cs_volume_zone.h"
56 #include "cs_xdef.h"
57 #include "cs_zone.h"
58
59 /*----------------------------------------------------------------------------
60 * Header for the current file
61 *----------------------------------------------------------------------------*/
62
63 #include "cs_advection_field.h"
64
65 /*----------------------------------------------------------------------------*/
66
67 BEGIN_C_DECLS
68
69 /*=============================================================================
70 * Local Macro definitions and structure definitions
71 *============================================================================*/
72
73 /* Redefined names of function from cs_math to get shorter names */
74 #define _dp3 cs_math_3_dot_product
75
76 #define CS_ADVECTION_FIELD_DBG 0
77
78 #define CS_ADVECTION_FIELD_ID_NOT_SET -1
79 #define CS_ADVECTION_FIELD_ID_TO_BE_SET -2
80
81 /*============================================================================
82 * Private variables
83 *============================================================================*/
84
85 static const char _err_empty_adv[] =
86 " Stop setting an empty cs_adv_field_t structure.\n"
87 " Please check your settings.\n";
88
89 /* Pointer to shared structures (owned by a cs_domain_t structure) */
90 static const cs_cdo_quantities_t *cs_cdo_quant;
91 static const cs_cdo_connect_t *cs_cdo_connect;
92
93 /* Advection fields attached to the computational domain */
94 static int _n_adv_fields = 0;
95 static cs_adv_field_t **_adv_fields = NULL;
96
97 /*============================================================================
98 * Inline private function prototypes
99 *============================================================================*/
100
101 /*----------------------------------------------------------------------------*/
102 /*!
103 * \brief Return the required dimension for the definition of an advection
104 * field
105 *
106 * \param[in] adv pointer to an advection field structure
107 *
108 * \return the dimension
109 */
110 /*----------------------------------------------------------------------------*/
111
112 static inline int
_get_dim_def(const cs_adv_field_t * adv)113 _get_dim_def(const cs_adv_field_t *adv)
114 {
115 int dim = -1;
116
117 if (adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR)
118 dim = 3;
119 else if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX)
120 dim = 1;
121 else
122 bft_error(__FILE__, __LINE__, 0,
123 "%s: Invalid dimension for the advection field.", __func__);
124
125 return dim;
126 }
127
128 /*============================================================================
129 * Private function prototypes
130 *============================================================================*/
131
132 /*----------------------------------------------------------------------------*/
133 /*!
134 * \brief Update the contribution of the flux
135 *
136 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
137 * \param[in] f2e face --> edge connectivity
138 * \param[in] e2v edge --> vertices connectivity
139 * \param[in] bf_id boundary face id
140 * \param[in] face_flux value of the normal flux across the face
141 * \param[in, out] fluxes array of values to update
142 */
143 /*----------------------------------------------------------------------------*/
144
145 static void
_fill_uniform_boundary_flux(const cs_cdo_quantities_t * const cdoq,const cs_adjacency_t * const f2e,const cs_adjacency_t * const e2v,cs_lnum_t bf_id,cs_real_t face_flux,cs_real_t * fluxes)146 _fill_uniform_boundary_flux(const cs_cdo_quantities_t *const cdoq,
147 const cs_adjacency_t *const f2e,
148 const cs_adjacency_t *const e2v,
149 cs_lnum_t bf_id,
150 cs_real_t face_flux,
151 cs_real_t *fluxes)
152 {
153 const cs_real_t face_coef = 0.5*face_flux/cdoq->b_face_surf[bf_id];
154 const cs_real_t *xf = cdoq->b_face_center + 3*bf_id;
155 const cs_lnum_t f_id = cdoq->n_i_faces + bf_id;
156
157 for (cs_lnum_t i = f2e->idx[f_id]; i < f2e->idx[f_id+1]; i++) {
158
159 const cs_lnum_t eshift = 2*f2e->ids[i];
160 const cs_lnum_t v0 = e2v->ids[eshift];
161 const cs_lnum_t v1 = e2v->ids[eshift+1];
162 const double tef = cs_math_surftri(cdoq->vtx_coord + 3*v0,
163 cdoq->vtx_coord + 3*v1,
164 xf);
165 const double weighted_flux = tef * face_coef;
166
167 fluxes[v0] += weighted_flux;
168 fluxes[v1] += weighted_flux;
169
170 } /* Loop on face edges */
171
172 }
173
174 /*----------------------------------------------------------------------------*/
175 /*!
176 * \brief Update the contribution of the flux for each vertex
177 *
178 * \param[in] cm pointer to a cs_cell_mesh_t structure
179 * \param[in] f face id in the cellwise numbering
180 * \param[in] face_flux value of the normal flux across the face
181 * \param[in, out] fluxes normal boundary flux for each vertex of the face
182 */
183 /*----------------------------------------------------------------------------*/
184
185 static void
_cw_fill_uniform_boundary_flux(const cs_cell_mesh_t * cm,short int f,cs_real_t face_flux,cs_real_t * fluxes)186 _cw_fill_uniform_boundary_flux(const cs_cell_mesh_t *cm,
187 short int f,
188 cs_real_t face_flux,
189 cs_real_t *fluxes)
190 {
191 const cs_real_t face_coef = 0.5*face_flux/cm->face[f].meas;
192
193 for (short int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
194
195 const short int eshift = 2*cm->f2e_ids[i];
196 const short int v0 = cm->e2v_ids[eshift];
197 const short int v1 = cm->e2v_ids[eshift+1];
198 const double weighted_flux = face_coef * cm->tef[i];
199
200 fluxes[v0] += weighted_flux;
201 fluxes[v1] += weighted_flux;
202
203 } /* Loop on face edges */
204 }
205
206 /*----------------------------------------------------------------------------*/
207 /*!
208 * \brief Compute a vector-valued reconstruction of the advection field at
209 * vertices from the definition of the advection field
210 *
211 * \param[in] quant additional mesh quantities struct.
212 * \param[in] connect pointer to a cs_cdo_connect_t struct.
213 * \param[in] def pointer to a cs_xdef_t struct.
214 * \param[in, out] vtx_values values of the reconstruction
215 */
216 /*----------------------------------------------------------------------------*/
217
218 static void
_compute_adv_vector_at_vertices(const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,cs_xdef_t * def,cs_real_t * vtx_values)219 _compute_adv_vector_at_vertices(const cs_cdo_quantities_t *quant,
220 const cs_cdo_connect_t *connect,
221 cs_xdef_t *def,
222 cs_real_t *vtx_values)
223 {
224 memset(vtx_values, 0, 3*quant->n_vertices*sizeof(cs_real_t));
225
226 const cs_adjacency_t *c2v = connect->c2v;
227
228 switch (def->type) {
229
230 case CS_XDEF_BY_ARRAY:
231 {
232 cs_real_t cell_vector[3];
233 cs_xdef_array_context_t *ai = (cs_xdef_array_context_t *)def->context;
234
235 if (cs_flag_test(ai->loc, cs_flag_dual_face_byc)) {
236
237 assert(ai->index == connect->c2e->idx);
238
239 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
240
241 cs_reco_dfbyc_at_cell_center(c_id,
242 connect->c2e,
243 quant,
244 ai->values,
245 cell_vector);
246
247 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++) {
248
249 const cs_real_t vc_vol = quant->dcell_vol[j];
250 cs_real_t *_val = vtx_values + 3*c2v->ids[j];
251
252 _val[0] += vc_vol * cell_vector[0];
253 _val[1] += vc_vol * cell_vector[1];
254 _val[2] += vc_vol * cell_vector[2];
255
256 }
257
258 } /* Loop on cells */
259
260 }
261 else if (cs_flag_test(ai->loc, cs_flag_primal_face)) {
262
263 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
264
265 cs_reco_cell_vector_by_face_dofs(c_id,
266 connect->c2f,
267 quant,
268 ai->values,
269 false, /* local input ? */
270 cell_vector);
271
272 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++) {
273
274 const cs_real_t vc_vol = quant->dcell_vol[j];
275 cs_real_t *_val = vtx_values + 3*c2v->ids[j];
276
277 _val[0] += vc_vol * cell_vector[0];
278 _val[1] += vc_vol * cell_vector[1];
279 _val[2] += vc_vol * cell_vector[2];
280
281 }
282
283 } /* Loop on cells */
284
285 }
286 else
287 bft_error(__FILE__, __LINE__, 0,
288 " %s: Invalid location for array", __func__);
289
290 }
291 break; /* by array */
292
293 case CS_XDEF_BY_DOF_FUNCTION:
294 {
295 #if defined(HAVE_OPENMP) /* Determine the default number of OpenMP threads */
296 int t_id = omp_get_thread_num();
297 #else
298 int t_id = 0;
299 #endif
300 cs_real_t cell_vector[3];
301 cs_real_t *fluxes = cs_cdo_local_get_d_buffer(t_id);
302 cs_xdef_dof_context_t *cx = (cs_xdef_dof_context_t *)def->context;
303
304 if (cs_flag_test(cx->loc, cs_flag_primal_face) == false)
305 bft_error(__FILE__, __LINE__, 0,
306 "%s: Invalid location for definition by DoFs.\n", __func__);
307
308 const cs_lnum_t *idx = connect->c2f->idx;
309
310 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
311
312 cs_lnum_t n_fc = idx[c_id+1] - idx[c_id];
313 const cs_lnum_t *_ids = connect->c2f->ids + idx[c_id];
314
315 /* Values of the function are defined at the primal faces */
316 cx->func(n_fc, _ids, true, /* dense output */
317 cx->input,
318 fluxes);
319
320 cs_reco_cell_vector_by_face_dofs(c_id,
321 connect->c2f,
322 quant,
323 fluxes,
324 true, /* local input ? */
325 cell_vector);
326
327 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++) {
328
329 const cs_real_t vc_vol = quant->dcell_vol[j];
330 cs_real_t *_val = vtx_values + 3*c2v->ids[j];
331
332 _val[0] += vc_vol * cell_vector[0];
333 _val[1] += vc_vol * cell_vector[1];
334 _val[2] += vc_vol * cell_vector[2];
335
336 }
337
338 } /* Loop on cells */
339
340 }
341 break;
342
343 default:
344 bft_error(__FILE__, __LINE__, 0, " %s: Incompatible type of definition.",
345 __func__);
346 break;
347
348 } /* Type of definition */
349
350 }
351
352 /*============================================================================
353 * Public function prototypes
354 *============================================================================*/
355
356 /*----------------------------------------------------------------------------*/
357 /*!
358 * \brief Set shared pointers to main domain members
359 *
360 * \param[in] quant additional mesh quantities struct.
361 * \param[in] connect pointer to a cs_cdo_connect_t struct.
362 */
363 /*----------------------------------------------------------------------------*/
364
365 void
cs_advection_field_set_shared_pointers(const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect)366 cs_advection_field_set_shared_pointers(const cs_cdo_quantities_t *quant,
367 const cs_cdo_connect_t *connect)
368 {
369 /* Assign static const pointers */
370 cs_cdo_quant = quant;
371 cs_cdo_connect = connect;
372 }
373
374 /*----------------------------------------------------------------------------*/
375 /*!
376 * \brief Get the number of allocated cs_adv_field_t structures
377 *
378 * \return the number of advection fields
379 */
380 /*----------------------------------------------------------------------------*/
381
382 int
cs_advection_field_get_n_fields(void)383 cs_advection_field_get_n_fields(void)
384 {
385 return _n_adv_fields;
386 }
387
388 /*----------------------------------------------------------------------------*/
389 /*!
390 * \brief Search in the array of advection field structures which one has
391 * the name given in argument
392 *
393 * \param[in] name name of the advection field
394 *
395 * \return a pointer to a cs_adv_field_t structure or NULL if not found
396 */
397 /*----------------------------------------------------------------------------*/
398
399 cs_adv_field_t *
cs_advection_field_by_name(const char * name)400 cs_advection_field_by_name(const char *name)
401 {
402 if (_n_adv_fields <= 0)
403 return NULL;
404
405 for (int i = 0; i < _n_adv_fields; i++) {
406
407 cs_adv_field_t *adv = _adv_fields[i];
408 if (cs_advection_field_check_name(adv, name))
409 return adv;
410
411 }
412
413 return NULL;
414 }
415
416 /*----------------------------------------------------------------------------*/
417 /*!
418 * \brief Search in the array of advection field structures which one has
419 * the id given in argument
420 *
421 * \param[in] id identification number
422 *
423 * \return a pointer to a cs_adv_field_t structure or NULL if not found
424 */
425 /*----------------------------------------------------------------------------*/
426
427 cs_adv_field_t *
cs_advection_field_by_id(int id)428 cs_advection_field_by_id(int id)
429 {
430 if (_n_adv_fields <= 0)
431 return NULL;
432 if (id >= _n_adv_fields || id < 0)
433 return NULL;
434 if (_adv_fields == NULL)
435 return NULL;
436
437 return _adv_fields[id];
438 }
439
440 /*----------------------------------------------------------------------------*/
441 /*!
442 * \brief Add and initialize a new user-defined advection field structure
443 *
444 * \param[in] name name of the advection field
445 *
446 * \return a pointer to the new allocated cs_adv_field_t structure
447 */
448 /*----------------------------------------------------------------------------*/
449
450 cs_adv_field_t *
cs_advection_field_add_user(const char * name)451 cs_advection_field_add_user(const char *name)
452 {
453 return cs_advection_field_add(name, CS_ADVECTION_FIELD_USER);
454 }
455
456 /*----------------------------------------------------------------------------*/
457 /*!
458 * \brief Add and initialize a new advection field structure
459 *
460 * \param[in] name name of the advection field
461 * \param[in] status status of the advection field to create
462 *
463 * \return a pointer to the new allocated cs_adv_field_t structure
464 */
465 /*----------------------------------------------------------------------------*/
466
467 cs_adv_field_t *
cs_advection_field_add(const char * name,cs_advection_field_status_t status)468 cs_advection_field_add(const char *name,
469 cs_advection_field_status_t status)
470 {
471 if (name == NULL)
472 bft_error(__FILE__, __LINE__, 0,
473 "%s: A non-empty name is mandatory to add a new advection field",
474 __func__);
475
476 cs_adv_field_t *adv = cs_advection_field_by_name(name);
477 if (adv != NULL) {
478 cs_base_warn(__FILE__, __LINE__);
479 cs_log_printf(CS_LOG_DEFAULT,
480 _(" An existing advection field has already the name %s.\n"
481 " Stop adding this advection field.\n"), name);
482 return adv;
483 }
484
485 /* Sanity check */
486 if (!(status & CS_ADVECTION_FIELD_USER) &&
487 !(status & CS_ADVECTION_FIELD_GWF) &&
488 !(status & CS_ADVECTION_FIELD_NAVSTO))
489 bft_error(__FILE__, __LINE__, 0,
490 "%s: No category associated to the advection field %s.",
491 __func__, name);
492
493 int new_id = _n_adv_fields;
494 _n_adv_fields++;
495 BFT_REALLOC(_adv_fields, _n_adv_fields, cs_adv_field_t *);
496 _adv_fields[new_id] = NULL;
497
498 BFT_MALLOC(adv, 1, cs_adv_field_t);
499
500 /* Copy name */
501 size_t len = strlen(name);
502 BFT_MALLOC(adv->name, len + 1, char);
503 strncpy(adv->name, name, len);
504 adv->name[len] = '\0';
505
506 adv->id = new_id;
507 adv->status = status;
508 adv->post_flag = 0;
509
510 /* Default initialization */
511 adv->definition = NULL;
512
513 adv->n_bdy_flux_defs = 0;
514 adv->bdy_flux_defs = NULL;
515 adv->bdy_def_ids = NULL;
516
517 adv->cell_field_id = CS_ADVECTION_FIELD_ID_NOT_SET;
518 adv->int_field_id = CS_ADVECTION_FIELD_ID_NOT_SET;
519
520 adv->vtx_field_id = CS_ADVECTION_FIELD_ID_NOT_SET;
521 if (adv->status & CS_ADVECTION_FIELD_DEFINE_AT_VERTICES)
522 adv->vtx_field_id = CS_ADVECTION_FIELD_ID_TO_BE_SET;
523
524 adv->bdy_field_id = CS_ADVECTION_FIELD_ID_NOT_SET;
525 if (adv->status & CS_ADVECTION_FIELD_DEFINE_AT_BOUNDARY_FACES)
526 adv->bdy_field_id = CS_ADVECTION_FIELD_ID_TO_BE_SET;
527
528 /* Update the status to a scalar flux in case of an advection field coming
529 from the legacy FV part */
530 if ((status & CS_ADVECTION_FIELD_NAVSTO) &&
531 (status & CS_ADVECTION_FIELD_LEGACY_FV))
532 adv->status |= CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX;
533
534 /* If no type is set, set vector-valued by default */
535 if (!(status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX) &&
536 !(status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR))
537 adv->status |= CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR;
538
539 /* Store the new advection field */
540 _adv_fields[new_id] = adv;
541
542 return adv;
543 }
544
545 /*----------------------------------------------------------------------------*/
546 /*!
547 * \brief Free all allocated cs_adv_field_t structures and its related array
548 *
549 * \return a NULL pointer
550 */
551 /*----------------------------------------------------------------------------*/
552
553 void
cs_advection_field_destroy_all(void)554 cs_advection_field_destroy_all(void)
555 {
556 if (_adv_fields == NULL)
557 return;
558
559 for (int i = 0; i < _n_adv_fields; i++) {
560
561 cs_adv_field_t *adv = _adv_fields[i];
562
563 adv->definition = cs_xdef_free(adv->definition);
564
565 for (int id = 0; id < adv->n_bdy_flux_defs; id++)
566 adv->bdy_flux_defs[id] = cs_xdef_free(adv->bdy_flux_defs[id]);
567
568 if (adv->n_bdy_flux_defs > 0) BFT_FREE(adv->bdy_flux_defs);
569 if (adv->bdy_def_ids != NULL) BFT_FREE(adv->bdy_def_ids);
570
571 BFT_FREE(adv->name);
572 BFT_FREE(adv);
573
574 /* All other pointers are shared */
575
576 } /* Loop on advection fields */
577
578 BFT_FREE(_adv_fields);
579 _n_adv_fields = 0;
580 }
581
582 /*----------------------------------------------------------------------------*/
583 /*!
584 * \brief Check if the given advection field has the name ref_name
585 *
586 * \param[in] adv pointer to a cs_adv_field_t structure to test
587 * \param[in] ref_name name of the advection field to find
588 *
589 * \return true if the name of the advection field is ref_name otherwise false
590 */
591 /*----------------------------------------------------------------------------*/
592
593 bool
cs_advection_field_check_name(const cs_adv_field_t * adv,const char * ref_name)594 cs_advection_field_check_name(const cs_adv_field_t *adv,
595 const char *ref_name)
596 {
597 if (adv == NULL)
598 return false;
599
600 int reflen = strlen(ref_name);
601 int len = strlen(adv->name);
602
603 if (reflen == len) {
604 if (strcmp(ref_name, adv->name) == 0)
605 return true;
606 else
607 return false;
608 }
609 else
610 return false;
611 }
612
613 /*----------------------------------------------------------------------------*/
614 /*!
615 * \brief Print all setup information related to cs_adv_field_t structures
616 */
617 /*----------------------------------------------------------------------------*/
618
619 void
cs_advection_field_log_setup(void)620 cs_advection_field_log_setup(void)
621 {
622 if (_adv_fields == NULL)
623 return;
624
625 char prefix[256];
626
627 cs_log_printf(CS_LOG_SETUP, "\nSummary of the advection field\n");
628 cs_log_printf(CS_LOG_SETUP, "%s\n", cs_sep_h1);
629
630 for (int i = 0; i < _n_adv_fields; i++) {
631
632 const cs_adv_field_t *adv = _adv_fields[i];
633
634 if (adv == NULL)
635 continue;
636 assert(strlen(adv->name) < 200);
637
638 /* Category of advection field */
639 cs_log_printf(CS_LOG_SETUP, " * %s | Category: ", adv->name);
640 if (adv->status & CS_ADVECTION_FIELD_NAVSTO)
641 cs_log_printf(CS_LOG_SETUP, "Related to Navier-Stokes\n");
642 else if (adv->status & CS_ADVECTION_FIELD_GWF)
643 cs_log_printf(CS_LOG_SETUP,
644 "Related to the \"Groundwater Flow\" module\n");
645 else if (CS_ADVECTION_FIELD_USER)
646 cs_log_printf(CS_LOG_SETUP, "User-defined\n");
647
648 /* Type of advection field */
649 cs_log_printf(CS_LOG_SETUP, " * %s | Type: ", adv->name);
650 if (adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR)
651 cs_log_printf(CS_LOG_SETUP, "Velocity vector\n");
652 else if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX)
653 cs_log_printf(CS_LOG_SETUP, "Scalar flux\n");
654
655 /* Additional information */
656 if (adv->status & CS_ADVECTION_FIELD_LEGACY_FV)
657 cs_log_printf(CS_LOG_SETUP, " * %s | %s\n",
658 adv->name, "Related to Legacy FV schemes\n");
659 if (adv->status & CS_ADVECTION_FIELD_STEADY)
660 cs_log_printf(CS_LOG_SETUP,
661 " * %s | Time status: Steady-state\n", adv->name);
662 else
663 cs_log_printf(CS_LOG_SETUP,
664 " * %s | Time status: Unsteady\n", adv->name);
665
666 if (adv->post_flag & CS_ADVECTION_FIELD_POST_COURANT)
667 cs_log_printf(CS_LOG_SETUP,
668 " * %s | Postprocess the Courant number\n", adv->name);
669
670 /* Where fields are defined */
671 bool at_cells =
672 (adv->cell_field_id > CS_ADVECTION_FIELD_ID_NOT_SET) ? true : false;
673 bool at_vertices =
674 (adv->vtx_field_id > CS_ADVECTION_FIELD_ID_NOT_SET) ? true : false;
675 bool at_bfaces =
676 (adv->bdy_field_id > CS_ADVECTION_FIELD_ID_NOT_SET) ? true : false;
677 bool at_ifaces =
678 (adv->int_field_id > CS_ADVECTION_FIELD_ID_NOT_SET) ? true : false;
679
680 cs_log_printf(CS_LOG_SETUP,
681 " * %s | Fields defined at cells: %s; at vertices: %s\n",
682 adv->name, cs_base_strtf(at_cells),
683 cs_base_strtf(at_vertices));
684 cs_log_printf(CS_LOG_SETUP,
685 " * %s | Fields defined at boundary faces: %s;"
686 " at interior faces: %s\n\n",
687 adv->name, cs_base_strtf(at_bfaces),
688 cs_base_strtf(at_ifaces));
689
690 sprintf(prefix, " Definition");
691 cs_xdef_log(prefix, adv->definition);
692
693 /* Boundary flux definition */
694 cs_log_printf(CS_LOG_SETUP,
695 " * %s | Number of boundary flux definitions: %d\n",
696 adv->name, adv->n_bdy_flux_defs);
697 if (adv->n_bdy_flux_defs > 0)
698 cs_log_printf(CS_LOG_SETUP, "\n");
699 for (int ib = 0; ib < adv->n_bdy_flux_defs; ib++) {
700 sprintf(prefix, " Definition %2d", ib);
701 cs_xdef_log(prefix, adv->bdy_flux_defs[ib]);
702 }
703
704 } /* Loop on advection fields */
705
706 }
707
708 /*----------------------------------------------------------------------------*/
709 /*!
710 * \brief Set optional post-processings
711 *
712 * \param[in, out] adv pointer to a cs_adv_field_t structure
713 * \param[in] post_flag flag to set
714 */
715 /*----------------------------------------------------------------------------*/
716
717 void
cs_advection_field_set_postprocess(cs_adv_field_t * adv,cs_flag_t post_flag)718 cs_advection_field_set_postprocess(cs_adv_field_t *adv,
719 cs_flag_t post_flag)
720 {
721 if (adv == NULL)
722 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
723
724 adv->post_flag |= post_flag;
725 }
726
727 /*----------------------------------------------------------------------------*/
728 /*!
729 * \brief Define the value of a cs_adv_field_t structure
730 *
731 * \param[in, out] adv pointer to a cs_adv_field_t structure
732 * \param[in] values values to set
733 */
734 /*----------------------------------------------------------------------------*/
735
736 void
cs_advection_field_def_by_value(cs_adv_field_t * adv,cs_real_t * values)737 cs_advection_field_def_by_value(cs_adv_field_t *adv,
738 cs_real_t *values)
739 {
740 if (adv == NULL)
741 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
742
743 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM |
744 CS_FLAG_STATE_CELLWISE | CS_FLAG_STATE_STEADY;
745 cs_flag_t meta_flag = CS_FLAG_FULL_LOC;
746 int dim = _get_dim_def(adv);
747
748 adv->definition = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
749 dim,
750 0, /* zone_id = 0 => all cells */
751 state_flag,
752 meta_flag,
753 values);
754 }
755
756 /*----------------------------------------------------------------------------*/
757 /*!
758 * \brief Define a cs_adv_field_t structure thanks to an analytic function
759 *
760 * \param[in, out] adv pointer to a cs_adv_field_t structure
761 * \param[in] func pointer to a function
762 * \param[in] input NULL or pointer to a structure cast on-the-fly
763 */
764 /*----------------------------------------------------------------------------*/
765
766 void
cs_advection_field_def_by_analytic(cs_adv_field_t * adv,cs_analytic_func_t * func,void * input)767 cs_advection_field_def_by_analytic(cs_adv_field_t *adv,
768 cs_analytic_func_t *func,
769 void *input)
770 {
771 if (adv == NULL)
772 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
773
774 cs_flag_t state_flag = 0;
775 cs_flag_t meta_flag = CS_FLAG_FULL_LOC;
776 int dim = _get_dim_def(adv);
777 cs_xdef_analytic_context_t ac = { .z_id = 0, /* all cells */
778 .func = func,
779 .input = input,
780 .free_input = NULL };
781
782 adv->definition = cs_xdef_volume_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
783 dim,
784 0, /* zone_id = 0 => all cells */
785 state_flag,
786 meta_flag,
787 &ac);
788 }
789
790 /*----------------------------------------------------------------------------*/
791 /*!
792 * \brief Define a cs_adv_field_t structure using a cs_dof_func_t
793 *
794 * \param[in, out] adv pointer to a cs_adv_field_t structure
795 * \param[in] loc_flag location where values are computed
796 * \param[in] func pointer to a cs_dof_func_t function
797 * \param[in] input NULL or pointer to a structure cast on-the-fly
798 */
799 /*----------------------------------------------------------------------------*/
800
801 void
cs_advection_field_def_by_dof_func(cs_adv_field_t * adv,cs_flag_t loc_flag,cs_dof_func_t * func,void * input)802 cs_advection_field_def_by_dof_func(cs_adv_field_t *adv,
803 cs_flag_t loc_flag,
804 cs_dof_func_t *func,
805 void *input)
806 {
807 if (adv == NULL)
808 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
809
810 cs_flag_t state_flag = 0;
811 cs_flag_t meta_flag = CS_FLAG_FULL_LOC;
812 int dim = _get_dim_def(adv);
813 cs_xdef_dof_context_t cx = { .z_id = 0, /* all cells */
814 .loc = loc_flag,
815 .func = func,
816 .input = input,
817 .free_input = NULL };
818
819 adv->definition = cs_xdef_volume_create(CS_XDEF_BY_DOF_FUNCTION,
820 dim,
821 0, /* zone_id = 0 => all cells */
822 state_flag,
823 meta_flag,
824 &cx);
825 }
826
827 /*----------------------------------------------------------------------------*/
828 /*!
829 * \brief Define a cs_adv_field_t structure thanks to an array of values
830 *
831 * \param[in, out] adv pointer to a cs_adv_field_t structure
832 * \param[in] loc information to know where are located values
833 * \param[in] array pointer to an array
834 * \param[in] is_owner transfer the lifecycle to the cs_xdef_t structure
835 * (true or false)
836 * \param[in] index optional pointer to the array index
837 */
838 /*----------------------------------------------------------------------------*/
839
840 void
cs_advection_field_def_by_array(cs_adv_field_t * adv,cs_flag_t loc,cs_real_t * array,bool is_owner,cs_lnum_t * index)841 cs_advection_field_def_by_array(cs_adv_field_t *adv,
842 cs_flag_t loc,
843 cs_real_t *array,
844 bool is_owner,
845 cs_lnum_t *index)
846 {
847 if (adv == NULL)
848 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
849
850 cs_flag_t state_flag = 0; /* Will be updated during the creation */
851 cs_flag_t meta_flag = CS_FLAG_FULL_LOC;
852 cs_xdef_array_context_t context = {.z_id = 0, /* all cells */
853 .stride = 3, /* default initialization */
854 .loc = loc,
855 .values = array,
856 .is_owner = is_owner,
857 .index = index };
858
859 /* Set the stride accord to the status (flux or vector) */
860 context.stride = _get_dim_def(adv);
861
862 /* Checkings */
863 if ((loc & CS_FLAG_SCALAR) && context.stride == 3)
864 bft_error(__FILE__, __LINE__, 0,
865 "%s: Incompatible setting for advection field %s\n"
866 " Array is set as a flux while the advection field as a vector.",
867 __func__, adv->name);
868
869 if ((loc & CS_FLAG_VECTOR) && context.stride == 1)
870 bft_error(__FILE__, __LINE__, 0,
871 "%s: Incompatible setting for advection field %s\n"
872 " Array is set as a vector while the advection field as a flux.",
873 __func__, adv->name);
874
875 adv->definition = cs_xdef_volume_create(CS_XDEF_BY_ARRAY,
876 context.stride,
877 0, /* zone_id = all cells */
878 state_flag,
879 meta_flag,
880 (void *)&context);
881 }
882
883 /*----------------------------------------------------------------------------*/
884 /*!
885 * \brief Define a cs_adv_field_t structure thanks to a field structure
886 *
887 * \param[in, out] adv pointer to a cs_adv_field_t structure
888 * \param[in] field pointer to a cs_field_t structure
889 */
890 /*----------------------------------------------------------------------------*/
891
892 void
cs_advection_field_def_by_field(cs_adv_field_t * adv,cs_field_t * field)893 cs_advection_field_def_by_field(cs_adv_field_t *adv,
894 cs_field_t *field)
895 {
896 if (adv == NULL)
897 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
898
899 /* Flags will be updated during the creation */
900 cs_flag_t state_flag = 0;
901 cs_flag_t meta_flag = 0;
902
903 /* Set the stride accord to the status (flux or vector) */
904 int dim = _get_dim_def(adv);
905 if (field->dim != dim)
906 bft_error(__FILE__, __LINE__, 0,
907 " %s: Inconsistency found between the field dimension and the"
908 " definition of the advection field %s.\n",
909 __func__, adv->name);
910
911 adv->definition = cs_xdef_volume_create(CS_XDEF_BY_FIELD,
912 dim,
913 0, /* zone_id */
914 state_flag,
915 meta_flag,
916 field);
917 }
918
919 /*----------------------------------------------------------------------------*/
920 /*!
921 * \brief Define the value of the boundary normal flux for the given
922 * cs_adv_field_t structure
923 *
924 * \param[in, out] adv pointer to a cs_adv_field_t structure
925 * \param[in] zname name of the boundary zone to consider
926 * \param[in] normal_flux value to set
927 */
928 /*----------------------------------------------------------------------------*/
929
930 void
cs_advection_field_def_boundary_flux_by_value(cs_adv_field_t * adv,const char * zname,cs_real_t normal_flux)931 cs_advection_field_def_boundary_flux_by_value(cs_adv_field_t *adv,
932 const char *zname,
933 cs_real_t normal_flux)
934 {
935 if (adv == NULL)
936 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
937
938 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM | CS_FLAG_STATE_FACEWISE;
939 cs_flag_t meta_flag = 0;
940
941 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_VALUE,
942 1, /* dim. */
943 cs_get_bdy_zone_id(zname),
944 state_flag,
945 meta_flag,
946 (void *)&normal_flux);
947
948 int def_id = adv->n_bdy_flux_defs;
949 adv->n_bdy_flux_defs += 1;
950 BFT_REALLOC(adv->bdy_flux_defs, adv->n_bdy_flux_defs, cs_xdef_t *);
951 adv->bdy_flux_defs[def_id] = d;
952 }
953
954 /*----------------------------------------------------------------------------*/
955 /*!
956 * \brief Define the value of the boundary normal flux for the given
957 * cs_adv_field_t structure using an analytic function
958 *
959 * \param[in, out] adv pointer to a cs_adv_field_t structure
960 * \param[in] zname name of the boundary zone to consider
961 * \param[in] func pointer to a function
962 * \param[in] input NULL or pointer to a structure cast on-the-fly
963 */
964 /*----------------------------------------------------------------------------*/
965
966 void
cs_advection_field_def_boundary_flux_by_analytic(cs_adv_field_t * adv,const char * zname,cs_analytic_func_t * func,void * input)967 cs_advection_field_def_boundary_flux_by_analytic(cs_adv_field_t *adv,
968 const char *zname,
969 cs_analytic_func_t *func,
970 void *input)
971 {
972 if (adv == NULL)
973 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
974
975 cs_flag_t state_flag = 0, meta_flag = 0;
976 int z_id = cs_get_bdy_zone_id(zname);
977 cs_xdef_analytic_context_t ac = { .z_id = z_id,
978 .func = func,
979 .input = input,
980 .free_input = NULL };
981
982 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
983 1, /* dim. */
984 z_id,
985 state_flag,
986 meta_flag,
987 &ac);
988
989 int def_id = adv->n_bdy_flux_defs;
990 adv->n_bdy_flux_defs += 1;
991 BFT_REALLOC(adv->bdy_flux_defs, adv->n_bdy_flux_defs, cs_xdef_t *);
992 adv->bdy_flux_defs[def_id] = d;
993 }
994
995 /*----------------------------------------------------------------------------*/
996 /*!
997 * \brief Define the value of the boundary normal flux for the given
998 * cs_adv_field_t structure using an array of values
999 *
1000 * \param[in, out] adv pointer to a cs_adv_field_t structure
1001 * \param[in] zname name of the boundary zone to consider
1002 * \param[in] loc information to know where are located values
1003 * \param[in] array pointer to an array
1004 * \param[in] is_owner transfer the lifecycle to the cs_xdef_t structure
1005 * (true or false)
1006 * \param[in] index optional pointer to the array index
1007 */
1008 /*----------------------------------------------------------------------------*/
1009
1010 void
cs_advection_field_def_boundary_flux_by_array(cs_adv_field_t * adv,const char * zname,cs_flag_t loc,cs_real_t * array,bool is_owner,cs_lnum_t * index)1011 cs_advection_field_def_boundary_flux_by_array(cs_adv_field_t *adv,
1012 const char *zname,
1013 cs_flag_t loc,
1014 cs_real_t *array,
1015 bool is_owner,
1016 cs_lnum_t *index)
1017 {
1018 if (adv == NULL)
1019 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
1020
1021 if (loc & CS_FLAG_VECTOR)
1022 bft_error(__FILE__, __LINE__, 0,
1023 "%s: Advection field: %s\n"
1024 " The boundary flux is not compatible with a vector-valued"
1025 " definition.\n", __func__, adv->name);
1026
1027 cs_flag_t state_flag = 0;
1028 cs_flag_t meta_flag = 0;
1029
1030 int z_id = cs_get_bdy_zone_id(zname);
1031 if (z_id == 0)
1032 meta_flag |= CS_FLAG_FULL_LOC;
1033
1034 cs_xdef_array_context_t context = {.z_id = z_id,
1035 .stride = 1,
1036 .loc = loc,
1037 .values = array,
1038 .is_owner = is_owner,
1039 .index = index };
1040
1041 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_ARRAY,
1042 1, /* dim. */
1043 z_id,
1044 state_flag,
1045 meta_flag,
1046 (void *)&context);
1047
1048 int def_id = adv->n_bdy_flux_defs;
1049 adv->n_bdy_flux_defs += 1;
1050 BFT_REALLOC(adv->bdy_flux_defs, adv->n_bdy_flux_defs, cs_xdef_t *);
1051 adv->bdy_flux_defs[def_id] = d;
1052 }
1053
1054 /*----------------------------------------------------------------------------*/
1055 /*!
1056 * \brief Define the value of the boundary normal flux for the given
1057 * cs_adv_field_t structure using a field structure
1058 *
1059 * \param[in, out] adv pointer to a cs_adv_field_t structure
1060 * \param[in] field pointer to a cs_field_t structure
1061 */
1062 /*----------------------------------------------------------------------------*/
1063
1064 void
cs_advection_field_def_boundary_flux_by_field(cs_adv_field_t * adv,cs_field_t * field)1065 cs_advection_field_def_boundary_flux_by_field(cs_adv_field_t *adv,
1066 cs_field_t *field)
1067 {
1068 if (adv == NULL)
1069 bft_error(__FILE__, __LINE__, 0, _(_err_empty_adv));
1070
1071 /* Flags will be updated during the creation */
1072 cs_flag_t state_flag = 0;
1073 cs_flag_t meta_flag = CS_FLAG_FULL_LOC; /* all boundary faces are handled */
1074
1075 /* Set the stride accord to the status (flux or vector) */
1076 if (field->dim != 1)
1077 bft_error(__FILE__, __LINE__, 0,
1078 " %s: Inconsistency found in the field dimension.\n"
1079 " A flux is requested (dim = 1) for advection field %s\n",
1080 __func__, adv->name);
1081
1082 cs_xdef_t *d = cs_xdef_boundary_create(CS_XDEF_BY_FIELD,
1083 1, /* dim. */
1084 0, /* all boundary faces */
1085 state_flag,
1086 meta_flag,
1087 field);
1088
1089 int def_id = adv->n_bdy_flux_defs;
1090 assert(def_id == 0);
1091 adv->n_bdy_flux_defs += 1;
1092 BFT_REALLOC(adv->bdy_flux_defs, adv->n_bdy_flux_defs, cs_xdef_t *);
1093 adv->bdy_flux_defs[def_id] = d;
1094 }
1095
1096 /*----------------------------------------------------------------------------*/
1097 /*!
1098 * \brief Create all needed cs_field_t structures related to an advection
1099 * field
1100 */
1101 /*----------------------------------------------------------------------------*/
1102
1103 void
cs_advection_field_create_fields(void)1104 cs_advection_field_create_fields(void)
1105 {
1106 int len;
1107 char *field_name = NULL;
1108
1109 for (int i = 0; i < _n_adv_fields; i++) {
1110
1111 cs_adv_field_t *adv = _adv_fields[i];
1112 assert(adv != NULL);
1113
1114 bool has_previous = (adv->status & CS_ADVECTION_FIELD_STEADY) ? true:false;
1115 int field_mask = CS_FIELD_PROPERTY | CS_FIELD_CDO;
1116
1117 /* Always add a field attached to cells (it may be used to define the
1118 numerical flux for advection, to compute adimensional numbers or to
1119 postprocess the advection field in case of definition by flux */
1120 if (adv->cell_field_id < 0) {
1121
1122 if (cs_flag_test(adv->status, CS_ADVECTION_FIELD_NAVSTO)) {
1123
1124 /* If this is the advection field related to the Navier-Stokes equations
1125 then there is no need to allocate a new buffer. One links to the
1126 existing velocity field */
1127 adv->cell_field_id = cs_field_id_by_name("velocity");
1128 assert(adv->cell_field_id > -1);
1129
1130 }
1131 else {
1132
1133 /* Define the name of the new field related the cell array */
1134 len = strlen(adv->name) + strlen("_cells") + 1;
1135 BFT_MALLOC(field_name, len, char);
1136 sprintf(field_name, "%s_cells", adv->name);
1137
1138 cs_field_t *fld = cs_field_create(field_name,
1139 field_mask,
1140 CS_MESH_LOCATION_CELLS,
1141 3, /* always a vector-valued field */
1142 has_previous);
1143
1144 cs_field_set_key_int(fld, cs_field_key_id("log"), 1);
1145 cs_field_set_key_int(fld, cs_field_key_id("post_vis"), 1);
1146
1147 adv->cell_field_id = cs_field_id_by_name(field_name);
1148
1149 BFT_FREE(field_name);
1150
1151 }
1152
1153 } /* Add a field at cells */
1154
1155 if (adv->vtx_field_id == CS_ADVECTION_FIELD_ID_TO_BE_SET) {
1156
1157 /* The creation of a field to store the value at vertices has been
1158 requested. Add this field */
1159 len = strlen(adv->name) + strlen("_vertices") + 1;
1160 BFT_MALLOC(field_name, len, char);
1161 sprintf(field_name, "%s_vertices", adv->name);
1162
1163 cs_field_t *fld = cs_field_create(field_name,
1164 field_mask,
1165 CS_MESH_LOCATION_VERTICES,
1166 3, /* always a vector-valued field */
1167 has_previous);
1168
1169 cs_field_set_key_int(fld, cs_field_key_id("log"), 1);
1170 cs_field_set_key_int(fld, cs_field_key_id("post_vis"), 1);
1171
1172 adv->vtx_field_id = cs_field_id_by_name(field_name);
1173
1174 BFT_FREE(field_name);
1175
1176 } /* Add a field attached to vertices */
1177
1178 if (adv->bdy_field_id == CS_ADVECTION_FIELD_ID_TO_BE_SET) {
1179
1180 /* Create a new field at the boundary faces for taking
1181 into account the normal flux used in the treatment of the boundary
1182 conditions for instance */
1183 len = strlen(adv->name) + strlen("_boundary_flux") + 1;
1184 BFT_MALLOC(field_name, len, char);
1185 sprintf(field_name, "%s_boundary_flux", adv->name);
1186
1187 cs_field_t *fld = cs_field_create(field_name,
1188 field_mask,
1189 CS_MESH_LOCATION_BOUNDARY_FACES,
1190 1, /* always a scalar-valued field */
1191 has_previous);
1192
1193 cs_field_set_key_int(fld, cs_field_key_id("log"), 1);
1194 cs_field_set_key_int(fld, cs_field_key_id("post_vis"), 1);
1195
1196 adv->bdy_field_id = cs_field_id_by_name(field_name);
1197
1198 BFT_FREE(field_name);
1199
1200 } /* Add a field attached to boundary faces */
1201
1202 } /* Loop on advection fields */
1203
1204 }
1205
1206 /*----------------------------------------------------------------------------*/
1207 /*!
1208 * \brief Last stage of the definition of an advection field based on several
1209 * definitions (i.e. definition by subdomains on the boundary)
1210 */
1211 /*----------------------------------------------------------------------------*/
1212
1213 void
cs_advection_field_finalize_setup(void)1214 cs_advection_field_finalize_setup(void)
1215 {
1216 if (_n_adv_fields == 0)
1217 return;
1218
1219 for (int i = 0; i < _n_adv_fields; i++) {
1220
1221 cs_adv_field_t *adv = _adv_fields[i];
1222 assert(adv != NULL);
1223
1224 /* Case of an advection field defined as the mass flux from the legacy FV
1225 part */
1226 if ((adv->status & CS_ADVECTION_FIELD_NAVSTO) &&
1227 (adv->status & CS_ADVECTION_FIELD_LEGACY_FV)) {
1228
1229 /* Automatic definition and checkings */
1230 cs_field_t *fld = cs_field_by_name("boundary_mass_flux");
1231 assert(fld != NULL);
1232 adv->bdy_field_id = fld->id;
1233
1234 if (adv->bdy_flux_defs == NULL)
1235 cs_advection_field_def_boundary_flux_by_field(adv, fld);
1236
1237 else { /* Sanity check */
1238 if (adv->n_bdy_flux_defs > 1 ||
1239 adv->bdy_flux_defs[0]->type != CS_XDEF_BY_FIELD)
1240 bft_error(__FILE__, __LINE__, 0,
1241 "%s: Invalid setting found for the advection field %s\n"
1242 " No need to perform additional setting when the legacy"
1243 " FV mass flux is used.\n", __func__, adv->name);
1244 }
1245
1246 fld = cs_field_by_name("inner_mass_flux");
1247 assert(fld != NULL);
1248 cs_advection_field_def_by_field(adv, fld);
1249 adv->int_field_id = fld->id;
1250
1251 if (adv->definition == NULL)
1252 cs_advection_field_def_by_field(adv, fld);
1253
1254 else { /* Sanity check */
1255 if (adv->definition->type != CS_XDEF_BY_FIELD)
1256 bft_error(__FILE__, __LINE__, 0,
1257 "%s: Invalid setting found for the advection field %s\n"
1258 " No need to perform additional setting when the legacy"
1259 " FV mass flux is used.\n", __func__, adv->name);
1260 }
1261
1262 }
1263
1264 if (adv->n_bdy_flux_defs > 1) {
1265
1266 const cs_lnum_t n_b_faces = cs_cdo_quant->n_b_faces;
1267
1268 BFT_MALLOC(adv->bdy_def_ids, n_b_faces, short int);
1269 # pragma omp parallel for if (n_b_faces > CS_THR_MIN)
1270 for (cs_lnum_t j = 0; j < n_b_faces; j++)
1271 adv->bdy_def_ids[j] = -1; /* Unset by default */
1272
1273 for (short int def_id = 0; def_id < adv->n_bdy_flux_defs; def_id++) {
1274
1275 const cs_xdef_t *def = adv->bdy_flux_defs[def_id];
1276 assert(def->z_id > 0);
1277 assert(def->support == CS_XDEF_SUPPORT_BOUNDARY);
1278
1279 const cs_zone_t *z = cs_boundary_zone_by_id(def->z_id);
1280 assert(z != NULL);
1281
1282 # pragma omp parallel for if (z->n_elts > CS_THR_MIN)
1283 for (cs_lnum_t j = 0; j < z->n_elts; j++)
1284 adv->bdy_def_ids[z->elt_ids[j]] = def_id;
1285
1286 }
1287
1288 } /* More than one definition */
1289
1290 } /* Loop on advection fields */
1291
1292 }
1293
1294 /*----------------------------------------------------------------------------*/
1295 /*!
1296 * \brief Compute the value of the advection field at the cell center
1297 *
1298 * \param[in] c_id id of the current cell
1299 * \param[in] adv pointer to a cs_adv_field_t structure
1300 * \param[in, out] vect pointer to a cs_nvec3_t structure (meas + unitv)
1301 */
1302 /*----------------------------------------------------------------------------*/
1303
1304 void
cs_advection_field_get_cell_vector(cs_lnum_t c_id,const cs_adv_field_t * adv,cs_nvec3_t * vect)1305 cs_advection_field_get_cell_vector(cs_lnum_t c_id,
1306 const cs_adv_field_t *adv,
1307 cs_nvec3_t *vect)
1308 {
1309 /* Initialize the vector */
1310 vect->meas = 0.;
1311 for (int k = 0; k < 3; k++)
1312 vect->unitv[k] = 0;
1313
1314 if (adv == NULL)
1315 return;
1316
1317 cs_field_t *f = cs_advection_field_get_field(adv, CS_MESH_LOCATION_CELLS);
1318
1319 cs_nvec3(f->val + 3*c_id, vect);
1320 }
1321
1322 /*----------------------------------------------------------------------------*/
1323 /*!
1324 * \brief Compute the vector-valued interpolation of the advection field at
1325 * a given location inside a cell
1326 *
1327 * \param[in] adv pointer to a cs_adv_field_t structure
1328 * \param[in] cm pointer to a cs_cell_mesh_t structure
1329 * \param[in] xyz location where to perform the evaluation
1330 * \param[in] time_eval physical time at which one evaluates the term
1331 * \param[in, out] eval pointer to a cs_nvec3_t
1332 */
1333 /*----------------------------------------------------------------------------*/
1334
1335 void
cs_advection_field_cw_eval_at_xyz(const cs_adv_field_t * adv,const cs_cell_mesh_t * cm,const cs_real_3_t xyz,cs_real_t time_eval,cs_nvec3_t * eval)1336 cs_advection_field_cw_eval_at_xyz(const cs_adv_field_t *adv,
1337 const cs_cell_mesh_t *cm,
1338 const cs_real_3_t xyz,
1339 cs_real_t time_eval,
1340 cs_nvec3_t *eval)
1341 {
1342 if (adv == NULL)
1343 return;
1344
1345 cs_xdef_t *def = adv->definition;
1346 cs_real_3_t vector_values = {0, 0, 0};
1347
1348 if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX) {
1349
1350 assert(def->dim == 1);
1351 switch (def->type) {
1352
1353 case CS_XDEF_BY_ARRAY: /* P0 approximation */
1354 {
1355 cs_xdef_array_context_t *ai = (cs_xdef_array_context_t *)def->context;
1356
1357 if (cs_flag_test(ai->loc, cs_flag_dual_face_byc)) {
1358
1359 assert(ai->index == cs_cdo_connect->c2e->idx);
1360 cs_reco_dfbyc_in_cell(cm,
1361 ai->values + ai->index[cm->c_id],
1362 vector_values);
1363
1364 }
1365 else if (cs_flag_test(ai->loc, cs_flag_primal_face)) {
1366
1367 const cs_real_t *i_flux = ai->values;
1368 const cs_real_t *b_flux = ai->values + cs_cdo_quant->n_i_faces;
1369
1370 cs_reco_cw_cell_vect_from_face_dofs(cm,
1371 i_flux,
1372 b_flux,
1373 vector_values);
1374 }
1375 else
1376 bft_error(__FILE__, __LINE__, 0,
1377 " %s: Invalid location for array", __func__);
1378
1379 cs_nvec3(vector_values, eval);
1380 }
1381 break;
1382
1383 case CS_XDEF_BY_FIELD: /* P0 approximation */
1384 {
1385 const cs_real_t *i_flux = (cs_field_by_id(adv->int_field_id))->val;
1386 const cs_real_t *b_flux = (cs_field_by_id(adv->bdy_field_id))->val;
1387
1388 cs_reco_cw_cell_vect_from_face_dofs(cm,
1389 i_flux,
1390 b_flux,
1391 vector_values);
1392 cs_nvec3(vector_values, eval);
1393 }
1394 break;
1395
1396 case CS_XDEF_BY_DOF_FUNCTION: /* P0 approximation */
1397 {
1398 #if defined(HAVE_OPENMP) /* Determine the default number of OpenMP threads */
1399 int t_id = omp_get_thread_num();
1400 #else
1401 int t_id = 0;
1402 #endif
1403 cs_real_t *fluxes = cs_cdo_local_get_d_buffer(t_id);
1404 cs_xdef_dof_context_t *cx = (cs_xdef_dof_context_t *)def->context;
1405
1406 /* Values of the function are defined at the primal faces */
1407 if (cs_flag_test(cx->loc, cs_flag_primal_face))
1408 cx->func(cm->n_fc, cm->f_ids, true, /* dense output */
1409 cx->input, fluxes);
1410 else
1411 bft_error(__FILE__, __LINE__, 0,
1412 "%s: Invalid location for definition by DoFs.\n", __func__);
1413
1414 cs_reco_cw_cell_vect_from_flux(cm, fluxes, vector_values);
1415 cs_nvec3(vector_values, eval);
1416 }
1417 break;
1418
1419 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1420 case CS_XDEF_BY_VALUE:
1421 bft_error(__FILE__, __LINE__, 0,
1422 " %s: Type of definition not implemented yet.", __func__);
1423 break;
1424
1425 default:
1426 bft_error(__FILE__, __LINE__, 0,
1427 " %s: Incompatible type of definition.", __func__);
1428 break;
1429
1430 } /* Type of definition */
1431
1432 }
1433 else {
1434
1435 assert(adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR);
1436 assert(def->dim == 3);
1437
1438 switch (def->type) {
1439
1440 case CS_XDEF_BY_VALUE:
1441 {
1442 const cs_real_t *constant_val = (cs_real_t *)def->context;
1443
1444 cs_nvec3(constant_val, eval);
1445 }
1446 break; /* definition by value */
1447
1448 case CS_XDEF_BY_ARRAY: /* P0 approximation */
1449 cs_xdef_cw_eval_vector_at_xyz_by_array(cm, 1, xyz, time_eval,
1450 def->context,
1451 vector_values);
1452 cs_nvec3(vector_values, eval);
1453 break;
1454
1455 case CS_XDEF_BY_FIELD:
1456 if (adv->vtx_field_id < 0 && adv->cell_field_id < 0)
1457 bft_error(__FILE__, __LINE__, 0,
1458 "%s: This support is not available for this function.\n",
1459 __func__);
1460
1461 cs_xdef_cw_eval_vector_at_xyz_by_field(cm, 1, xyz, time_eval,
1462 def->context,
1463 vector_values);
1464 cs_nvec3(vector_values, eval);
1465 break;
1466
1467 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1468 cs_xdef_cw_eval_at_xyz_by_analytic(cm, 1, xyz, time_eval, def->context,
1469 vector_values);
1470 cs_nvec3(vector_values, eval);
1471 break;
1472
1473 default:
1474 bft_error(__FILE__, __LINE__, 0, " %s: Incompatible type of definition.",
1475 __func__);
1476 break;
1477
1478 } /* Type of definition */
1479
1480 } /* Velocity vector */
1481 }
1482
1483 /*----------------------------------------------------------------------------*/
1484 /*!
1485 * \brief Compute the mean-value of the vector-valued field related to the
1486 * advection field inside each cell
1487 *
1488 * \param[in] adv pointer to a cs_adv_field_t structure
1489 * \param[in] time_eval physical time at which one evaluates the term
1490 * \param[in, out] cell_values array of values at cell centers
1491 */
1492 /*----------------------------------------------------------------------------*/
1493
1494 void
cs_advection_field_in_cells(const cs_adv_field_t * adv,cs_real_t time_eval,cs_real_t * cell_values)1495 cs_advection_field_in_cells(const cs_adv_field_t *adv,
1496 cs_real_t time_eval,
1497 cs_real_t *cell_values)
1498 {
1499 if (adv == NULL)
1500 return;
1501
1502 const cs_cdo_quantities_t *cdoq = cs_cdo_quant;
1503
1504 cs_xdef_t *def = adv->definition;
1505 assert(def != NULL); /* Sanity checks */
1506
1507 if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX) {
1508
1509 switch (def->type) {
1510
1511 case CS_XDEF_BY_ARRAY:
1512 {
1513 cs_xdef_array_context_t *ai = (cs_xdef_array_context_t *)def->context;
1514
1515 if (cs_flag_test(ai->loc, cs_flag_dual_face_byc)) {
1516
1517 assert(ai->index == cs_cdo_connect->c2e->idx);
1518
1519 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
1520 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++)
1521 cs_reco_dfbyc_at_cell_center(c_id,
1522 cs_cdo_connect->c2e,
1523 cdoq,
1524 ai->values,
1525 cell_values + 3*c_id);
1526
1527 }
1528 else if (cs_flag_test(ai->loc, cs_flag_primal_face)) {
1529
1530 assert(adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX);
1531
1532 cs_reco_cell_vectors_by_face_dofs(cs_cdo_connect->c2f,
1533 cdoq,
1534 ai->values,
1535 cell_values);
1536
1537 }
1538 else
1539 bft_error(__FILE__, __LINE__, 0,
1540 " %s: Invalid location for array", __func__);
1541
1542 }
1543 break;
1544
1545 case CS_XDEF_BY_FIELD:
1546 {
1547 const cs_field_t *field = (cs_field_t *)def->context;
1548 assert(field != NULL);
1549 assert(field->dim == 1);
1550 const cs_mesh_location_type_t loc_type =
1551 cs_mesh_location_get_type(field->location_id);
1552
1553 if (loc_type == CS_MESH_LOCATION_INTERIOR_FACES) {
1554
1555 /* Sanity checks */
1556 assert(adv->bdy_field_id > -1);
1557 cs_field_t *b_mflx_fld = cs_field_by_id(adv->bdy_field_id);
1558 assert(b_mflx_fld != NULL);
1559
1560 cs_reco_cell_vectors_by_ib_face_dofs(cs_cdo_connect->c2f,
1561 cdoq,
1562 field->val,
1563 b_mflx_fld->val,
1564 cell_values);
1565 }
1566 else
1567 bft_error(__FILE__, __LINE__, 0,
1568 " %s: Invalid support for the input field", __func__);
1569 }
1570 break;
1571
1572 case CS_XDEF_BY_DOF_FUNCTION: /* P0 approximation in each cell */
1573 {
1574 cs_real_t *fluxes = cs_equation_get_tmpbuf();
1575 cs_xdef_dof_context_t *cx = (cs_xdef_dof_context_t *)def->context;
1576
1577 if (cs_flag_test(cx->loc, cs_flag_primal_face) == false)
1578 bft_error(__FILE__, __LINE__, 0,
1579 "%s: Invalid location for definition by DoFs.\n", __func__);
1580 assert(cs_equation_get_tmpbuf_size() <= (size_t)cdoq->n_faces);
1581
1582 /* Values of the function are defined at the primal faces */
1583 cx->func(cdoq->n_faces, NULL, true, /* dense output */
1584 cx->input, fluxes);
1585
1586 cs_reco_cell_vectors_by_face_dofs(cs_cdo_connect->c2f,
1587 cdoq,
1588 fluxes,
1589 cell_values);
1590 }
1591 break;
1592
1593 default:
1594 bft_error(__FILE__, __LINE__, 0, "%s: Incompatible type of definition.",
1595 __func__);
1596 break;
1597
1598 } /* Type of definition */
1599
1600 }
1601 else {
1602
1603 assert(adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR);
1604
1605 switch (def->type) {
1606
1607 case CS_XDEF_BY_VALUE:
1608 {
1609 const cs_real_t *constant_val = (cs_real_t *)def->context;
1610
1611 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
1612 for (cs_lnum_t i = 0; i < cdoq->n_cells; i++) {
1613 cell_values[3*i ] = constant_val[0];
1614 cell_values[3*i+1] = constant_val[1];
1615 cell_values[3*i+2] = constant_val[2];
1616 }
1617 }
1618 break; /* definition by value */
1619
1620 case CS_XDEF_BY_ARRAY:
1621 {
1622 cs_xdef_array_context_t *ai = (cs_xdef_array_context_t *)def->context;
1623
1624 assert(ai->stride == 3);
1625 if (!cs_flag_test(ai->loc, cs_flag_primal_cell))
1626 bft_error(__FILE__, __LINE__, 0,
1627 " %s: Invalid location for array", __func__);
1628
1629 memcpy(cell_values, ai->values,
1630 ai->stride*cdoq->n_cells*sizeof(cs_real_t));
1631 }
1632 break; /* definition by array */
1633
1634 case CS_XDEF_BY_FIELD:
1635 {
1636 const cs_field_t *field = (cs_field_t *)def->context;
1637 assert(field != NULL);
1638 const cs_mesh_location_type_t loc_type =
1639 cs_mesh_location_get_type(field->location_id);
1640
1641 switch(loc_type) {
1642
1643 case CS_MESH_LOCATION_CELLS:
1644 assert(field->dim == 3);
1645 /* If not the field associated to the advection field */
1646 if (field->id != adv->cell_field_id)
1647 memcpy(cell_values, field->val, 3*cdoq->n_cells*sizeof(cs_real_t));
1648 break;
1649
1650 case CS_MESH_LOCATION_VERTICES:
1651 assert(field->dim == 3);
1652 cs_reco_vect_pv_at_cell_centers(cs_cdo_connect->c2v,
1653 cdoq,
1654 field->val,
1655 cell_values);
1656 break;
1657
1658 default:
1659 bft_error(__FILE__, __LINE__, 0,
1660 " %s: Invalid support for the input field", __func__);
1661 break;
1662
1663 } /* Switch on location type */
1664
1665 }
1666 break; /* definition by field */
1667
1668 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1669 cs_evaluate_average_on_cells_by_analytic(def, time_eval, cell_values);
1670 break; /* definition by analytic */
1671
1672 default:
1673 bft_error(__FILE__, __LINE__, 0, " %s: Incompatible type of definition.",
1674 __func__);
1675 break;
1676
1677 } /* Type of definition */
1678
1679 } /* Velocity vector */
1680 }
1681
1682 /*----------------------------------------------------------------------------*/
1683 /*!
1684 * \brief Compute the value of the advection field at vertices
1685 *
1686 * \param[in] adv pointer to a cs_adv_field_t structure
1687 * \param[in] time_eval physical time at which one evaluates the term
1688 * \param[in, out] vtx_values array storing the results
1689 */
1690 /*----------------------------------------------------------------------------*/
1691
1692 void
cs_advection_field_at_vertices(const cs_adv_field_t * adv,cs_real_t time_eval,cs_real_t * vtx_values)1693 cs_advection_field_at_vertices(const cs_adv_field_t *adv,
1694 cs_real_t time_eval,
1695 cs_real_t *vtx_values)
1696 {
1697 if (adv == NULL)
1698 return;
1699
1700 const cs_cdo_quantities_t *cdoq = cs_cdo_quant;
1701 const cs_cdo_connect_t *connect = cs_cdo_connect;
1702
1703 cs_xdef_t *def = adv->definition;
1704
1705 if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX) {
1706
1707 /* Compute the vector-valued vector at each vertex */
1708 _compute_adv_vector_at_vertices(cdoq, connect, def, vtx_values);
1709
1710 /* Synchronization of values at vertices */
1711 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL)
1712 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
1713 cdoq->n_vertices,
1714 3, /* stride */
1715 true, /* = interlace */
1716 CS_REAL_TYPE,
1717 vtx_values);
1718
1719 cs_real_t *dual_vol = NULL;
1720 BFT_MALLOC(dual_vol, cdoq->n_vertices, cs_real_t);
1721 cs_cdo_quantities_compute_dual_volumes(cdoq, connect->c2v, dual_vol);
1722
1723 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL)
1724 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
1725 cdoq->n_vertices,
1726 1, /* stride */
1727 true, /* = interlace */
1728 CS_REAL_TYPE,
1729 dual_vol);
1730
1731 # pragma omp parallel for if (cdoq->n_vertices > CS_THR_MIN)
1732 for (cs_lnum_t v_id = 0; v_id < cdoq->n_vertices; v_id++) {
1733 const cs_real_t invvol = 1./dual_vol[v_id];
1734 for (int k = 0; k < 3; k++)
1735 vtx_values[3*v_id+k] *= invvol;
1736 }
1737
1738 BFT_FREE(dual_vol);
1739
1740 }
1741 else {
1742
1743 assert(adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR);
1744
1745 switch (def->type) {
1746
1747 case CS_XDEF_BY_VALUE:
1748 {
1749 const cs_real_t *constant_val = (cs_real_t *)def->context;
1750
1751 # pragma omp parallel for if (cdoq->n_vertices > CS_THR_MIN)
1752 for (cs_lnum_t i = 0; i < cdoq->n_vertices; i++) {
1753 vtx_values[3*i ] = constant_val[0];
1754 vtx_values[3*i+1] = constant_val[1];
1755 vtx_values[3*i+2] = constant_val[2];
1756 }
1757
1758 }
1759 break; /* definition by value */
1760
1761 case CS_XDEF_BY_ARRAY:
1762 {
1763 cs_xdef_array_context_t *ai =
1764 (cs_xdef_array_context_t *)def->context;
1765 assert(ai->stride == 3);
1766
1767 if (cs_flag_test(ai->loc, cs_flag_primal_vtx))
1768 memcpy(vtx_values, ai->values, 3*cdoq->n_vertices*sizeof(cs_real_t));
1769
1770 else if (cs_flag_test(ai->loc, cs_flag_primal_cell))
1771 cs_reco_vect_pv_from_pc(connect->c2v, cdoq, ai->values, vtx_values);
1772
1773 else
1774 bft_error(__FILE__, __LINE__, 0,
1775 " %s: Invalid location for array", __func__);
1776 }
1777 break; /* definition by array */
1778
1779 case CS_XDEF_BY_FIELD:
1780 {
1781 cs_field_t *field = (cs_field_t *)def->context;
1782 assert(field != NULL);
1783 assert(field->dim == 3);
1784
1785 switch(cs_mesh_location_get_type(field->location_id)) {
1786
1787 case CS_MESH_LOCATION_CELLS:
1788 cs_reco_vect_pv_from_pc(cs_cdo_connect->c2v,
1789 cdoq,
1790 field->val,
1791 vtx_values);
1792 break;
1793
1794 case CS_MESH_LOCATION_VERTICES:
1795 /* If not the field associated to the advection field */
1796 if (field->id != adv->vtx_field_id)
1797 memcpy(vtx_values, field->val,
1798 3*cdoq->n_vertices*sizeof(cs_real_t));
1799 break;
1800
1801 default:
1802 bft_error(__FILE__, __LINE__, 0,
1803 " %s: Invalid case for the input field", __func__);
1804 break;
1805
1806 }
1807
1808 }
1809 break; /* definition by field */
1810
1811 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1812 cs_evaluate_potential_at_vertices_by_analytic(def,
1813 time_eval,
1814 cdoq->n_vertices,
1815 NULL,
1816 vtx_values);
1817 break; /* definition by analytic */
1818
1819 default:
1820 bft_error(__FILE__, __LINE__, 0, " %s: Incompatible type of definition.",
1821 __func__);
1822 break;
1823
1824 } /* Type of definition */
1825
1826 } /* velocity vector */
1827
1828 }
1829
1830 /*----------------------------------------------------------------------------*/
1831 /*!
1832 * \brief Compute the value of the normal flux of the advection field
1833 * across the boundary faces
1834 *
1835 * \param[in] adv pointer to a cs_adv_field_t structure
1836 * \param[in] time_eval physical time at which one evaluates the term
1837 * \param[in, out] flx_values array storing the results
1838 */
1839 /*----------------------------------------------------------------------------*/
1840
1841 void
cs_advection_field_across_boundary(const cs_adv_field_t * adv,cs_real_t time_eval,cs_real_t * flx_values)1842 cs_advection_field_across_boundary(const cs_adv_field_t *adv,
1843 cs_real_t time_eval,
1844 cs_real_t *flx_values)
1845 {
1846 if (adv == NULL)
1847 return;
1848 if (adv->status & CS_ADVECTION_FIELD_LEGACY_FV)
1849 return;
1850
1851 const cs_cdo_quantities_t *cdoq = cs_cdo_quant;
1852 const cs_lnum_t n_b_faces = cdoq->n_b_faces;
1853 const cs_lnum_t n_i_faces = cdoq->n_i_faces;
1854 const cs_xdef_t *def = adv->definition;
1855
1856 if (adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR) {
1857
1858 assert(def->dim == 3);
1859 if (adv->n_bdy_flux_defs > 0)
1860 bft_error(__FILE__, __LINE__, 0,
1861 "%s: Case where one expects a volume definition only.",
1862 __func__);
1863
1864 /* No specific definition of the boundary flux should be set.
1865 * So, one uses the definition related to the volume */
1866
1867 switch (def->type) {
1868
1869 case CS_XDEF_BY_VALUE:
1870 {
1871 const cs_real_t *constant_val = (cs_real_t *)def->context;
1872
1873 # pragma omp parallel for if (n_b_faces > CS_THR_MIN)
1874 for (cs_lnum_t i = 0; i < n_b_faces; i++)
1875 flx_values[i] = _dp3(cdoq->b_face_normal + 3*i, constant_val);
1876
1877 }
1878 break;
1879
1880 case CS_XDEF_BY_ARRAY:
1881 {
1882 cs_xdef_array_context_t *ai = (cs_xdef_array_context_t *)def->context;
1883 assert(ai->stride == 3);
1884
1885 if (cs_flag_test(ai->loc, cs_flag_primal_face)) {
1886
1887 /* Boundary faces come after interior faces */
1888 const cs_real_t *const bface_vel = ai->values + 3*n_i_faces;
1889
1890 # pragma omp parallel for if (n_b_faces > CS_THR_MIN)
1891 for (cs_lnum_t i = 0; i < n_b_faces; i++)
1892 flx_values[i] = _dp3(cdoq->b_face_normal + 3*i, bface_vel + 3*i);
1893
1894 }
1895 else if (cs_flag_test(ai->loc, cs_flag_primal_cell)) {
1896
1897 const cs_adjacency_t *f2c = cs_cdo_connect->f2c;
1898 const cs_lnum_t *const bf2c_ids = f2c->ids + 2*n_i_faces;
1899 const cs_real_t *const cell_vel = ai->values;
1900
1901 # pragma omp parallel for if (n_b_faces > CS_THR_MIN)
1902 for (cs_lnum_t i = 0; i < n_b_faces; i++) {
1903 const cs_lnum_t c_id = bf2c_ids[i];
1904 flx_values[i] = _dp3(cdoq->b_face_normal + 3*i, cell_vel + 3*c_id);
1905 }
1906
1907 }
1908 else
1909 bft_error(__FILE__, __LINE__, 0,
1910 " %s: Incompatible location when defined by array.",
1911 __func__);
1912
1913 }
1914 break;
1915
1916 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1917 {
1918 const cs_adjacency_t *f2e = cs_cdo_connect->f2e;
1919 const cs_adjacency_t *e2v = cs_cdo_connect->e2v;
1920 const cs_real_t *xv = cdoq->vtx_coord;
1921
1922 cs_quadrature_tria_integral_t
1923 *compute_integral = cs_quadrature_get_tria_integral(def->dim,
1924 def->qtype);
1925 cs_xdef_analytic_context_t *ac =
1926 (cs_xdef_analytic_context_t *)def->context;
1927
1928 for (cs_lnum_t i = 0; i < n_b_faces; i++) {
1929
1930 const cs_lnum_t f_id = n_i_faces + i;
1931 const cs_quant_t pfq = cs_quant_set_face(f_id, cdoq);
1932 const cs_lnum_t start_idx = f2e->idx[f_id];
1933 const cs_lnum_t end_idx = f2e->idx[f_id+1];
1934
1935 cs_real_t val[3] = {0, 0, 0};
1936
1937 switch (end_idx - start_idx) {
1938
1939 case CS_TRIANGLE_CASE: /* Triangle: one-shot computation */
1940 {
1941 cs_lnum_t v1, v2, v3;
1942
1943 cs_connect_get_next_3_vertices(f2e->ids, e2v->ids, start_idx,
1944 &v1, &v2, &v3);
1945
1946 compute_integral(time_eval,
1947 xv + 3*v1, xv + 3*v2, xv + 3*v3, pfq.meas,
1948 ac->func, ac->input, val);
1949 }
1950 break;
1951
1952 default:
1953 for (cs_lnum_t j = start_idx; j < end_idx; j++) {
1954
1955 const cs_lnum_t _2e = 2*f2e->ids[j];
1956 const cs_real_t *xv1 = xv + 3*e2v->ids[_2e];
1957 const cs_real_t *xv2 = xv + 3*e2v->ids[_2e+1];
1958 const cs_real_t tef_meas = cs_math_surftri(xv1, xv2, pfq.center);
1959
1960 /* val is updated (+=) */
1961 compute_integral(time_eval, xv1, xv2, pfq.center, tef_meas,
1962 ac->func, ac->input, val);
1963
1964 } /* Loop on face edges */
1965
1966 } /* End of switch */
1967
1968 flx_values[i] = _dp3(pfq.unitv, val);
1969
1970 } /* Loop on border faces */
1971
1972 }
1973 break; /* definition by analytic */
1974
1975 default:
1976 bft_error(__FILE__, __LINE__, 0, " %s: Incompatible type of definition.",
1977 __func__);
1978 break;
1979
1980 } /* Type of definition */
1981
1982 } /* Velocity vector */
1983 else {
1984
1985 assert(adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX);
1986
1987 if (adv->n_bdy_flux_defs == 0) {
1988
1989 switch (def->type) {
1990
1991 case CS_XDEF_BY_ARRAY:
1992 {
1993 const cs_xdef_array_context_t *ac =
1994 (cs_xdef_array_context_t *)def->context;
1995
1996 assert(ac->stride == 1);
1997 assert(def->meta & CS_FLAG_FULL_LOC); /* over the volume hence the
1998 shift with n_i_faces */
1999
2000 if (cs_flag_test(ac->loc, cs_flag_primal_face))
2001 memcpy(flx_values, ac->values + n_i_faces,
2002 sizeof(cs_real_t)*n_b_faces);
2003 }
2004 break;
2005
2006 case CS_XDEF_BY_DOF_FUNCTION:
2007 {
2008 cs_xdef_dof_context_t *cx = (cs_xdef_dof_context_t *)def->context;
2009
2010 if (cs_flag_test(cx->loc, cs_flag_primal_face) == false)
2011 bft_error(__FILE__, __LINE__, 0,
2012 "%s: Invalid location for definition by DoFs.\n",
2013 __func__);
2014
2015 cs_lnum_t *bface_ids = NULL;
2016 BFT_MALLOC(bface_ids, n_b_faces, cs_lnum_t);
2017 for (cs_lnum_t i = 0; i < n_b_faces; i++)
2018 bface_ids[i] = n_i_faces + i;
2019
2020 cx->func(n_b_faces, bface_ids, true, /* dense output */
2021 cx->input,
2022 flx_values);
2023
2024 BFT_FREE(bface_ids);
2025 }
2026 break;
2027
2028 default:
2029 bft_error(__FILE__, __LINE__, 0,
2030 "%s: Invalid type of definition in case of flux defined\n"
2031 "on the whole domain. It should a definition by array.",
2032 __func__);
2033 }
2034
2035 }
2036 else {
2037
2038 /* Consider the definition of the boundary flux for updating the field
2039 values */
2040 for (int def_id = 0; def_id < adv->n_bdy_flux_defs; def_id++) {
2041
2042 const cs_xdef_t *bdef = adv->bdy_flux_defs[def_id];
2043 const cs_zone_t *z = cs_boundary_zone_by_id(bdef->z_id);
2044
2045 assert(bdef->support == CS_XDEF_SUPPORT_BOUNDARY);
2046 assert(bdef->dim == 1);
2047
2048 switch (bdef->type) {
2049
2050 case CS_XDEF_BY_VALUE:
2051 {
2052 const cs_real_t *constant_val = (cs_real_t *)bdef->context;
2053
2054 if (z->elt_ids == NULL) {
2055 assert(z->n_elts == n_b_faces);
2056 # pragma omp parallel for if (n_b_faces > CS_THR_MIN)
2057 for (cs_lnum_t i = 0; i < n_b_faces; i++)
2058 flx_values[i] = constant_val[0];
2059 }
2060 else {
2061 # pragma omp parallel for if (z->n_elts > CS_THR_MIN)
2062 for (cs_lnum_t i = 0; i < z->n_elts; i++)
2063 flx_values[z->elt_ids[i]] = constant_val[0];
2064 }
2065
2066 }
2067 break;
2068
2069 case CS_XDEF_BY_ARRAY:
2070 {
2071 const cs_xdef_array_context_t *ac =
2072 (cs_xdef_array_context_t *)bdef->context;
2073 const cs_real_t *val = ac->values;
2074
2075 assert(ac->stride == 1);
2076 assert(bdef->meta & CS_FLAG_FULL_LOC || z->elt_ids == NULL);
2077
2078 if (cs_flag_test(ac->loc, cs_flag_primal_face))
2079 memcpy(flx_values, val, sizeof(cs_real_t)*n_b_faces);
2080
2081 else if (cs_flag_test(ac->loc, cs_flag_dual_closure_byf)) {
2082
2083 for (cs_lnum_t bf_id = 0; bf_id < n_b_faces; bf_id++) {
2084
2085 flx_values[bf_id] = 0;
2086 for (cs_lnum_t i = ac->index[bf_id]; i < ac->index[bf_id+1];
2087 i++)
2088 flx_values[bf_id] += val[i];
2089
2090 } /* Loop on border faces */
2091
2092 }
2093 else
2094 bft_error(__FILE__, __LINE__, 0, "%s: Invalid case.", __func__);
2095
2096 }
2097 break; /* definition by array */
2098
2099 case CS_XDEF_BY_FIELD:
2100 {
2101 cs_field_t *field = (cs_field_t *)bdef->context;
2102
2103 assert(field->dim == 1);
2104 assert(bdef->meta & CS_FLAG_FULL_LOC || z->elt_ids == NULL);
2105
2106 if (field->location_id ==
2107 cs_mesh_location_get_id_by_name(N_("boundary faces")))
2108 memcpy(flx_values, field->val, sizeof(cs_real_t)*n_b_faces);
2109 else
2110 bft_error(__FILE__, __LINE__, 0, " %s: Invalid case.", __func__);
2111
2112 }
2113 break; /* definition by field */
2114
2115 case CS_XDEF_BY_ANALYTIC_FUNCTION:
2116 {
2117 cs_xdef_analytic_context_t *ac =
2118 (cs_xdef_analytic_context_t *)bdef->context;
2119
2120 ac->func(time_eval,
2121 z->n_elts, z->elt_ids, cdoq->b_face_center,
2122 false, /* compacted output ? */
2123 ac->input,
2124 flx_values);
2125 }
2126 break; /* definition by analytic */
2127
2128 default:
2129 bft_error(__FILE__, __LINE__, 0,
2130 " %s: Incompatible type of definition.", __func__);
2131 break;
2132
2133 } /* Type of definition */
2134
2135 } /* Loop on boundary definitions */
2136
2137 } /* There is at least one boundary definition */
2138
2139 } /* Scalar-valued flux */
2140
2141 }
2142
2143 /*----------------------------------------------------------------------------*/
2144 /*!
2145 * \brief Compute the value of the normal flux of the advection field across
2146 * a boundary face f (cellwise version)
2147 *
2148 * \param[in] time_eval physical time at which one evaluates the term
2149 * \param[in] f face id in the cellwise numbering
2150 * \param[in] cm pointer to a cs_cell_mesh_t structure
2151 * \param[in] adv pointer to a cs_adv_field_t structure
2152 *
2153 * \return the normal boundary flux for the face f
2154 */
2155 /*----------------------------------------------------------------------------*/
2156
2157 cs_real_t
cs_advection_field_cw_boundary_face_flux(const cs_real_t time_eval,const short int f,const cs_cell_mesh_t * cm,const cs_adv_field_t * adv)2158 cs_advection_field_cw_boundary_face_flux(const cs_real_t time_eval,
2159 const short int f,
2160 const cs_cell_mesh_t *cm,
2161 const cs_adv_field_t *adv)
2162 {
2163 cs_real_t f_flux = 0.;
2164
2165 if (adv == NULL)
2166 return f_flux;
2167
2168 const cs_quant_t pfq = cm->face[f];
2169 const cs_lnum_t bf_id = cm->f_ids[f] - cm->bface_shift;
2170
2171 assert(bf_id > -1);
2172 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ));
2173
2174 if (adv->bdy_field_id > -1) { /* Use the current computed advection field
2175 across the boundary */
2176
2177 cs_field_t *fld = cs_field_by_id(adv->bdy_field_id);
2178
2179 return fld->val[bf_id];
2180
2181 }
2182
2183 const cs_xdef_t *def = adv->definition;
2184
2185 if (adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR) {
2186
2187 /* No specific definition of the boundary flux. Use the definition related
2188 to the volume */
2189 assert(adv->n_bdy_flux_defs == 0);
2190 assert(def->dim == 3);
2191
2192 switch (def->type) {
2193
2194 case CS_XDEF_BY_VALUE:
2195 {
2196 const cs_real_t *constant_val = (cs_real_t *)def->context;
2197
2198 f_flux = pfq.meas * _dp3(pfq.unitv, constant_val);
2199 }
2200 break;
2201
2202 case CS_XDEF_BY_ANALYTIC_FUNCTION:
2203 {
2204 cs_real_t adv_val[3] = {0, 0, 0};
2205 cs_quadrature_tria_integral_t *compute_integral =
2206 cs_quadrature_get_tria_integral(def->dim, def->qtype);
2207 cs_xdef_analytic_context_t *ac =
2208 (cs_xdef_analytic_context_t *)def->context;
2209
2210 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PV | CS_FLAG_COMP_FEQ |
2211 CS_FLAG_COMP_EV));
2212
2213 const cs_lnum_t start_idx = cm->f2e_idx[f];
2214 const cs_lnum_t end_idx = cm->f2e_idx[f+1];
2215 const short int n_vf = end_idx - start_idx; /* #vertices (=#edges) */
2216 const short int *f2e_ids = cm->f2e_ids + start_idx;
2217
2218 switch (n_vf) {
2219
2220 case CS_TRIANGLE_CASE: /* Triangle: one-shot computation */
2221 {
2222 short int v1, v2, v3;
2223 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids,
2224 &v1, &v2, &v3);
2225
2226 compute_integral(time_eval,
2227 cm->xv + 3*v1, cm->xv + 3*v2, cm->xv + 3*v3,
2228 pfq.meas, ac->func, ac->input,
2229 adv_val);
2230 }
2231 break;
2232
2233 default:
2234 {
2235 const double *tef = cm->tef + start_idx;
2236
2237 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
2238
2239 const short int _2e = 2*f2e_ids[e];
2240 const cs_real_t *xv1 = cm->xv + 3*cm->e2v_ids[_2e];
2241 const cs_real_t *xv2 = cm->xv + 3*cm->e2v_ids[_2e+1];
2242
2243 /* adv_val is updated (+=) */
2244 compute_integral(time_eval, xv1, xv2, pfq.center, tef[e],
2245 ac->func, ac->input,
2246 adv_val);
2247 } /* Loop on face edges */
2248
2249 }
2250
2251 } /* End of switch (n_vf) */
2252
2253 f_flux = _dp3(pfq.unitv, adv_val);
2254
2255 }
2256 break; /* definition by analytic */
2257
2258 default:
2259 bft_error(__FILE__, __LINE__, 0, " %s: Incompatible type of definition.",
2260 __func__);
2261 break;
2262
2263 } /* Type of definition */
2264
2265 }
2266 else {
2267
2268 assert(adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX);
2269 if (adv->n_bdy_flux_defs == 0) {
2270
2271 switch (def->type) {
2272
2273 case CS_XDEF_BY_ARRAY:
2274 {
2275 const cs_xdef_array_context_t *ac =
2276 (cs_xdef_array_context_t *)def->context;
2277
2278 assert(ac->stride == 1);
2279 if (cs_flag_test(ac->loc, cs_flag_primal_face))
2280 f_flux = ac->values[cm->f_ids[f]];
2281 else
2282 bft_error(__FILE__, __LINE__, 0,
2283 "%s: Invalid array support.", __func__);
2284 }
2285 break;
2286
2287 case CS_XDEF_BY_DOF_FUNCTION:
2288 {
2289 cs_xdef_dof_context_t *cx = (cs_xdef_dof_context_t *)def->context;
2290
2291 if (cs_flag_test(cx->loc, cs_flag_primal_face) == false)
2292 bft_error(__FILE__, __LINE__, 0,
2293 "%s: Invalid location for definition by DoFs.\n",
2294 __func__);
2295
2296 /* Values of the function are defined at the primal faces */
2297 cx->func(1, cm->f_ids + f, true, /* dense output */
2298 cx->input,
2299 &f_flux);
2300 }
2301 break;
2302
2303 default:
2304 bft_error(__FILE__, __LINE__, 0,
2305 "%s: Invalid type of definition.", __func__);
2306
2307 } /* def->type */
2308
2309 }
2310 else { /* n_bdy_flux_defs > 0 */
2311
2312 const cs_xdef_t *bdef = (adv->bdy_def_ids == NULL) ?
2313 adv->bdy_flux_defs[0] : adv->bdy_flux_defs[adv->bdy_def_ids[bf_id]];
2314
2315 assert(bdef != NULL);
2316 assert(bdef->support == CS_XDEF_SUPPORT_BOUNDARY);
2317
2318 switch (bdef->type) {
2319
2320 case CS_XDEF_BY_VALUE:
2321 {
2322 const cs_real_t *constant_val = (cs_real_t *)bdef->context;
2323 f_flux = constant_val[0];
2324 }
2325 break;
2326
2327 case CS_XDEF_BY_ARRAY:
2328 {
2329 const cs_xdef_array_context_t *ac =
2330 (cs_xdef_array_context_t *)bdef->context;
2331
2332 assert(ac->stride == 1);
2333 if (cs_flag_test(ac->loc, cs_flag_primal_face))
2334 f_flux = ac->values[bf_id];
2335
2336 else if (cs_flag_test(ac->loc, cs_flag_dual_closure_byf)) {
2337
2338 const cs_adjacency_t *bf2v = cs_cdo_connect->bf2v;
2339 for (cs_lnum_t i = bf2v->idx[bf_id]; i < bf2v->idx[bf_id+1]; i++)
2340 f_flux += ac->values[i];
2341
2342 }
2343 else
2344 bft_error(__FILE__, __LINE__, 0, " %s: Invalid support", __func__);
2345
2346 }
2347 break; /* definition by array */
2348
2349 case CS_XDEF_BY_FIELD:
2350 {
2351 cs_field_t *field = (cs_field_t *)bdef->context;
2352 assert(field->dim == 1);
2353
2354 if (field->location_id ==
2355 cs_mesh_location_get_id_by_name(N_("boundary faces")))
2356 f_flux = field->val[bf_id];
2357 else
2358 bft_error(__FILE__, __LINE__, 0, " %s: Invalid case.", __func__);
2359
2360 }
2361 break; /* definition by field */
2362
2363 case CS_XDEF_BY_ANALYTIC_FUNCTION:
2364 {
2365 cs_xdef_analytic_context_t *ac =
2366 (cs_xdef_analytic_context_t *)def->context;
2367
2368 ac->func(time_eval,
2369 1, NULL, pfq.center,
2370 true, /* compacted output ? */
2371 ac->input,
2372 &f_flux);
2373 }
2374 break; /* definition by analytic */
2375
2376 default:
2377 bft_error(__FILE__, __LINE__, 0,
2378 " %s: Incompatible type of definition.", __func__);
2379 break;
2380
2381 } /* Type of definition */
2382
2383 } /* There is at least one definition of the normal boundary flux */
2384
2385 } /* Scalar-valued flux */
2386
2387 return f_flux;
2388 }
2389
2390 /*----------------------------------------------------------------------------*/
2391 /*!
2392 * \brief Compute the value of the normal flux of the advection field across
2393 * the closure of the dual cell related to each vertex belonging to the
2394 * boundary face f
2395 *
2396 * \param[in] cm pointer to a cs_cell_mesh_t structure
2397 * \param[in] adv pointer to a cs_adv_field_t structure
2398 * \param[in] f face id in the cellwise numbering
2399 * \param[in] time_eval physical time at which one evaluates the term
2400 * \param[in, out] fluxes normal boundary flux for each vertex of the face
2401 */
2402 /*----------------------------------------------------------------------------*/
2403
2404 void
cs_advection_field_cw_boundary_f2v_flux(const cs_cell_mesh_t * cm,const cs_adv_field_t * adv,short int f,cs_real_t time_eval,cs_real_t * fluxes)2405 cs_advection_field_cw_boundary_f2v_flux(const cs_cell_mesh_t *cm,
2406 const cs_adv_field_t *adv,
2407 short int f,
2408 cs_real_t time_eval,
2409 cs_real_t *fluxes)
2410 {
2411 if (adv == NULL)
2412 return;
2413
2414 if (fluxes == NULL)
2415 bft_error(__FILE__, __LINE__, 0,
2416 " %s: Array of fluxes should be allocated before the call.",
2417 __func__);
2418
2419 /* Reset flux values */
2420 for (short int v = 0; v < cm->n_vc; v++) fluxes[v] = 0;
2421
2422 const cs_quant_t pfq = cm->face[f];
2423 const cs_lnum_t bf_id = cm->f_ids[f] - cm->bface_shift;
2424
2425 assert(bf_id > -1);
2426
2427 /* If the boundary flux has been computed, used this array expected in case of
2428 analytical definition on the boundary */
2429 if (adv->bdy_field_id > -1 && adv->n_bdy_flux_defs == 0) {
2430
2431 cs_field_t *fld = cs_field_by_id(adv->bdy_field_id);
2432
2433 _cw_fill_uniform_boundary_flux(cm, f, fld->val[bf_id], fluxes);
2434
2435 }
2436 else if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX) {
2437
2438 if (adv->n_bdy_flux_defs == 0) {
2439
2440 const cs_xdef_t *def = adv->definition;
2441
2442 if (def->type != CS_XDEF_BY_ARRAY)
2443 bft_error(__FILE__, __LINE__, 0,
2444 "%s: Invalid type of definition", __func__);
2445
2446 const cs_xdef_array_context_t *ac =
2447 (cs_xdef_array_context_t *)def->context;
2448
2449 assert(cs_flag_test(ac->loc, cs_flag_primal_face));
2450 _cw_fill_uniform_boundary_flux(cm, f, ac->values[bf_id], fluxes);
2451
2452 }
2453 else {
2454
2455 assert(cs_eflag_test(cm->flag,
2456 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FEQ |
2457 CS_FLAG_COMP_EV | CS_FLAG_COMP_FE));
2458
2459 /* Consider the definition of the boundary flux to set the values */
2460 const cs_xdef_t *def = (adv->bdy_def_ids == NULL) ?
2461 adv->bdy_flux_defs[0] : adv->bdy_flux_defs[adv->bdy_def_ids[bf_id]];
2462
2463 switch (def->type) {
2464
2465 case CS_XDEF_BY_VALUE:
2466 {
2467 const cs_real_t *constant_val = (cs_real_t *)def->context;
2468
2469 _cw_fill_uniform_boundary_flux(cm, f, constant_val[0], fluxes);
2470 }
2471 break; /* by value */
2472
2473 case CS_XDEF_BY_ARRAY:
2474 {
2475 const cs_xdef_array_context_t *ac
2476 = (cs_xdef_array_context_t *)def->context;
2477
2478 assert(ac->stride == 1);
2479 if (cs_flag_test(ac->loc, cs_flag_primal_face))
2480 _cw_fill_uniform_boundary_flux(cm, f, ac->values[bf_id], fluxes);
2481
2482 else if (cs_flag_test(ac->loc, cs_flag_dual_closure_byf)) {
2483
2484 const cs_adjacency_t *bf2v = cs_cdo_connect->bf2v;
2485 for (cs_lnum_t i = bf2v->idx[bf_id]; i < bf2v->idx[bf_id+1]; i++) {
2486 const short int v = cs_cell_mesh_get_v(bf2v->ids[i], cm);
2487 assert(v > -1 && v != cm->n_vc);
2488 fluxes[v] += ac->values[i];
2489 }
2490
2491 }
2492 else
2493 bft_error(__FILE__, __LINE__, 0, " %s: Invalid support", __func__);
2494 }
2495 break; /* by_array */
2496
2497 case CS_XDEF_BY_FIELD:
2498 {
2499 cs_field_t *field = (cs_field_t *)def->context;
2500 assert(field->dim == 1);
2501
2502 if (cs_mesh_location_get_type(field->location_id) ==
2503 CS_MESH_LOCATION_BOUNDARY_FACES)
2504 _cw_fill_uniform_boundary_flux(cm, f, field->val[bf_id], fluxes);
2505 else
2506 bft_error(__FILE__, __LINE__, 0, " %s: Invalid support", __func__);
2507 }
2508 break; /* definition by field */
2509
2510 case CS_XDEF_BY_ANALYTIC_FUNCTION:
2511 {
2512 cs_real_t f_flux = 0;
2513 cs_xdef_analytic_context_t *ac =
2514 (cs_xdef_analytic_context_t *)def->context;
2515
2516 ac->func(time_eval,
2517 1, NULL, pfq.center,
2518 true, /* compacted output ? */
2519 ac->input,
2520 &f_flux);
2521
2522 _cw_fill_uniform_boundary_flux(cm, f, f_flux, fluxes);
2523 }
2524 break; /* definition by analytic */
2525
2526 default:
2527 bft_error(__FILE__, __LINE__, 0,
2528 " %s: Invalid type of definition", __func__);
2529 break;
2530
2531 } /* End of switch on the type of definition */
2532
2533 }
2534
2535 }
2536 else {
2537
2538 assert(adv->n_bdy_flux_defs == 0);
2539 assert(adv->bdy_def_ids == NULL);
2540 assert(adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR);
2541
2542 /* Implicit definition of the boundary flux from the definition of the
2543 advection field inside the computational domain */
2544 cs_xdef_t *def = adv->definition;
2545
2546 assert(cs_eflag_test(cm->flag,
2547 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FEQ |
2548 CS_FLAG_COMP_EV | CS_FLAG_COMP_FE));
2549
2550 switch (def->type) {
2551
2552 case CS_XDEF_BY_ANALYTIC_FUNCTION:
2553 {
2554 cs_quadrature_tria_integral_t
2555 *compute_integral = cs_quadrature_get_tria_integral(def->dim,
2556 def->qtype);
2557 cs_xdef_analytic_context_t *ac =
2558 (cs_xdef_analytic_context_t *)def->context;
2559
2560 const short int end_idx = cm->f2e_idx[f+1];
2561 const short int start_idx = cm->f2e_idx[f];
2562 const short int n_ef = end_idx - start_idx; /* #vertices (= #edges) */
2563 const short int *f2e_ids = cm->f2e_ids + start_idx;
2564 const double *tef = cm->tef + start_idx;
2565
2566 for (short int e = 0; e < n_ef; e++) { /* Loop on face edges */
2567
2568 const short int _2e = 2*f2e_ids[e];
2569 const short int v1 = cm->e2v_ids[_2e];
2570 const short int v2 = cm->e2v_ids[_2e+1];
2571
2572 cs_real_t val[3] = {0, 0, 0};
2573
2574 /* val is updated (+=) */
2575 compute_integral(time_eval,
2576 cm->xv + 3*v1, cm->xv + 3*v2, pfq.center, tef[e],
2577 ac->func, ac->input, val);
2578
2579 const double e_contrib = 0.5*_dp3(pfq.unitv, val);
2580
2581 fluxes[v1] += e_contrib;
2582 fluxes[v2] += e_contrib;
2583
2584 } /* Loop on face edges */
2585 }
2586 break;
2587
2588 case CS_XDEF_BY_VALUE:
2589 {
2590 const cs_real_t *constant_val = (cs_real_t *)def->context;
2591 const cs_real_t half_face_density = 0.5*_dp3(pfq.unitv, constant_val);
2592
2593 for (short int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
2594
2595 const short int eshift = 2*cm->f2e_ids[i];
2596 const short int v0 = cm->e2v_ids[eshift];
2597 const short int v1 = cm->e2v_ids[eshift+1];
2598 const double weighted_flux = half_face_density * cm->tef[i];
2599
2600 fluxes[v0] += weighted_flux;
2601 fluxes[v1] += weighted_flux;
2602
2603 } /* Loop on face edges */
2604
2605 }
2606 break;
2607
2608 default:
2609 bft_error(__FILE__, __LINE__, 0,
2610 "%s: Invalid type of definition", __func__);
2611 break;
2612
2613 } /* End of switch */
2614
2615 } /* n_bdy_flux_defs == 0 && CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR */
2616
2617 }
2618
2619 /*----------------------------------------------------------------------------*/
2620 /*!
2621 * \brief Compute the value of the flux of the advection field across the
2622 * the (primal) faces of a cell
2623 *
2624 * \param[in] cm pointer to a cs_cell_mesh_t structure
2625 * \param[in] adv pointer to a cs_adv_field_t structure
2626 * \param[in] time_eval physical time at which one evaluates the term
2627 * \param[in, out] fluxes array of values attached to primal faces of a cell
2628 */
2629 /*----------------------------------------------------------------------------*/
2630
2631 void
cs_advection_field_cw_face_flux(const cs_cell_mesh_t * cm,const cs_adv_field_t * adv,cs_real_t time_eval,cs_real_t * fluxes)2632 cs_advection_field_cw_face_flux(const cs_cell_mesh_t *cm,
2633 const cs_adv_field_t *adv,
2634 cs_real_t time_eval,
2635 cs_real_t *fluxes)
2636 {
2637 if (adv == NULL)
2638 return;
2639
2640 if (fluxes == NULL) {
2641 bft_error(__FILE__, __LINE__, 0,
2642 " %s: The array of local fluxes should be already allocated.",
2643 __func__);
2644 return;
2645 }
2646
2647 cs_xdef_t *def = adv->definition;
2648 assert(def != NULL); /* Sanity check */
2649
2650 if (adv->status & CS_ADVECTION_FIELD_TYPE_SCALAR_FLUX) {
2651
2652 switch (def->type) {
2653
2654 case CS_XDEF_BY_ARRAY:
2655 {
2656 const cs_real_t *xdef_fluxes = cs_xdef_get_array(def);
2657 assert(fluxes != NULL);
2658 for (short int f = 0; f < cm->n_fc; f++)
2659 fluxes[f] = xdef_fluxes[cm->f_ids[f]];
2660
2661 }
2662 break;
2663
2664 case CS_XDEF_BY_FIELD:
2665 {
2666 cs_field_t *fld = (cs_field_t *)def->context;
2667
2668 /* Sanity check */
2669 assert(cs_flag_test(cm->flag, CS_FLAG_COMP_PFQ));
2670 assert(cs_mesh_location_get_type(fld->location_id)
2671 == CS_MESH_LOCATION_INTERIOR_FACES);
2672 assert(adv->bdy_field_id > -1);
2673
2674 cs_field_t *b_mflx_fld = cs_field_by_id(adv->bdy_field_id);
2675 assert(b_mflx_fld != NULL);
2676
2677 const cs_real_t *b_flux = b_mflx_fld->val;
2678 const cs_real_t *i_flux = fld->val;
2679
2680 /* Loop on cell faces */
2681 for (short int f = 0; f < cm->n_fc; f++) {
2682
2683 const cs_lnum_t f_id = cm->f_ids[f];
2684 if (cm->f_ids[f] < cm->bface_shift) /* Interior face */
2685 fluxes[f] = i_flux[f_id];
2686 else
2687 fluxes[f] = b_flux[f_id - cm->bface_shift];
2688
2689 } /* Loop on cell faces */
2690 }
2691 break;
2692
2693 case CS_XDEF_BY_DOF_FUNCTION:
2694 {
2695 cs_xdef_dof_context_t *cx = (cs_xdef_dof_context_t *)def->context;
2696
2697 assert(cs_flag_test(cx->loc, cs_flag_primal_face));
2698 cx->func(cm->n_fc, cm->f_ids, true, /* dense_output */
2699 cx->input,
2700 fluxes);
2701 }
2702 break;
2703
2704 default:
2705 bft_error(__FILE__, __LINE__, 0,
2706 "%s: Invalid type of definition", __func__);
2707
2708 } /* End of switch def->type */
2709
2710 }
2711 else {
2712
2713 assert(adv->status & CS_ADVECTION_FIELD_TYPE_VELOCITY_VECTOR);
2714 assert(def->dim == 3);
2715
2716 switch (def->type) {
2717
2718 case CS_XDEF_BY_VALUE:
2719 {
2720 /* Loop on cell faces */
2721 for (short int f = 0; f < cm->n_fc; f++)
2722 cs_xdef_cw_eval_flux_by_val(cm, f, time_eval, def->context, fluxes);
2723 }
2724 break;
2725
2726 case CS_XDEF_BY_ANALYTIC_FUNCTION:
2727 {
2728 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_FEQ | CS_FLAG_COMP_PFQ));
2729
2730 /* Loop on cell faces */
2731 for (short int f = 0; f < cm->n_fc; f++)
2732 cs_xdef_cw_eval_flux_by_analytic(cm, f, time_eval,
2733 def->context, def->qtype,
2734 fluxes);
2735 }
2736 break;
2737
2738 case CS_XDEF_BY_ARRAY:
2739 {
2740 cs_nvec3_t nv;
2741 cs_xdef_array_context_t *ac = (cs_xdef_array_context_t *)def->context;
2742
2743 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ));
2744 assert(ac->values != NULL);
2745
2746 /* Test if location has at least the pattern of the reference support */
2747 if (cs_flag_test(ac->loc, cs_flag_primal_face)) {
2748
2749 for (short int f = 0; f < cm->n_fc; f++) {
2750
2751 /* Retrieve the advection field:
2752 Switch to a cs_nvec3_t representation */
2753 cs_nvec3(ac->values + 3*cm->f_ids[f], &nv);
2754
2755 fluxes[f] = nv.meas*cm->face[f].meas * _dp3(nv.unitv,
2756 cm->face[f].unitv);
2757
2758 } /* Loop on cell faces */
2759
2760 }
2761 else if (cs_flag_test(ac->loc, cs_flag_primal_cell)) {
2762
2763 /* Retrieve the advection field:
2764 Switch to a cs_nvec3_t representation */
2765 cs_nvec3(ac->values + 3*cm->c_id, &nv);
2766
2767 /* Loop on cell faces */
2768 for (short int f = 0; f < cm->n_fc; f++)
2769 fluxes[f] = nv.meas*cm->face[f].meas * _dp3(nv.unitv,
2770 cm->face[f].unitv);
2771
2772 }
2773 else
2774 bft_error(__FILE__, __LINE__, 0,
2775 " %s: Invalid support for evaluating the advection field %s"
2776 " at the cell center of cell %ld.",
2777 __func__, adv->name, (long)cm->c_id);
2778 }
2779 break; /* Definition by array */
2780
2781 case CS_XDEF_BY_FIELD:
2782 {
2783 cs_field_t *fld = (cs_field_t *)def->context;
2784
2785 /* Sanity check */
2786 assert(cs_flag_test(cm->flag, CS_FLAG_COMP_PFQ));
2787 assert(cs_mesh_location_get_type(fld->location_id)
2788 == CS_MESH_LOCATION_CELLS);
2789
2790 /* Retrieve the advection field:
2791 Switch to a cs_nvec3_t representation */
2792 cs_nvec3_t nv;
2793 cs_nvec3(fld->val + 3*cm->c_id, &nv);
2794
2795 /* Loop on cell faces */
2796 for (short int f = 0; f < cm->n_fc; f++)
2797 fluxes[f] = nv.meas*cm->face[f].meas * _dp3(nv.unitv,
2798 cm->face[f].unitv);
2799 }
2800 break; /* Definition by field */
2801
2802 default:
2803 bft_error(__FILE__, __LINE__, 0,
2804 "%s: Incompatible type of definition.", __func__);
2805 break;
2806
2807 } /* def_type */
2808
2809 } /* Velocity vector */
2810
2811 }
2812
2813 /*----------------------------------------------------------------------------*/
2814 /*!
2815 * \brief Compute the value of the flux of the advection field across the
2816 * the dual faces of a cell
2817 *
2818 * \param[in] cm pointer to a cs_cell_mesh_t structure
2819 * \param[in] adv pointer to a cs_adv_field_t structure
2820 * \param[in] time_eval physical time at which one evaluates the term
2821 * \param[in, out] fluxes array of values attached to dual faces of a cell
2822 */
2823 /*----------------------------------------------------------------------------*/
2824
2825 void
cs_advection_field_cw_dface_flux(const cs_cell_mesh_t * cm,const cs_adv_field_t * adv,cs_real_t time_eval,cs_real_t * fluxes)2826 cs_advection_field_cw_dface_flux(const cs_cell_mesh_t *cm,
2827 const cs_adv_field_t *adv,
2828 cs_real_t time_eval,
2829 cs_real_t *fluxes)
2830 {
2831 if (adv == NULL)
2832 return;
2833
2834 if (fluxes == NULL) {
2835 bft_error(__FILE__, __LINE__, 0,
2836 " fluxes array should be allocated before the call.");
2837 return;
2838 }
2839
2840 cs_xdef_t *def = adv->definition;
2841
2842 /* Sanity checks */
2843 assert(def != NULL);
2844
2845 switch (def->type) {
2846
2847 case CS_XDEF_BY_VALUE:
2848 {
2849 /* Retrieve the advection field: Switch to a cs_nvec3_t representation */
2850 cs_real_3_t cell_vector;
2851 cs_xdef_cw_eval_vector_by_val(cm, time_eval, def->context, cell_vector);
2852 cs_nvec3_t nv;
2853 cs_nvec3(cell_vector, &nv);
2854
2855 /* Sanity check */
2856 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_DFQ | CS_FLAG_COMP_PEQ));
2857
2858 /* Loop on cell edges */
2859 for (short int e = 0; e < cm->n_ec; e++)
2860 fluxes[e] = nv.meas * cm->dface[e].meas * _dp3(nv.unitv,
2861 cm->dface[e].unitv);
2862 }
2863 break;
2864
2865 case CS_XDEF_BY_ANALYTIC_FUNCTION:
2866 {
2867 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PEQ | CS_FLAG_COMP_SEF |
2868 CS_FLAG_COMP_PFQ));
2869
2870 cs_quadrature_type_t qtype = cs_xdef_get_quadrature(def);
2871
2872 /* Reset fluxes across dual faces */
2873 memset(fluxes, 0, cm->n_ec*sizeof(cs_real_t));
2874
2875 switch (qtype) {
2876
2877 case CS_QUADRATURE_NONE:
2878 case CS_QUADRATURE_BARY:
2879 case CS_QUADRATURE_BARY_SUBDIV:
2880 {
2881 cs_real_3_t xg, adv_xg;
2882
2883 /* Loop on cell faces */
2884 for (short int f = 0; f < cm->n_fc; f++) {
2885
2886 const cs_quant_t face = cm->face[f];
2887 const cs_real_3_t xfc = { cm->xc[0] + face.center[0],
2888 cm->xc[1] + face.center[1],
2889 cm->xc[2] + face.center[2] };
2890
2891 for (short int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
2892
2893 const short int e = cm->f2e_ids[i];
2894 const cs_quant_t edge = cm->edge[e];
2895 const cs_nvec3_t sefc = cm->sefc[i];
2896
2897 for (int k = 0; k < 3; k++)
2898 xg[k] = cs_math_1ov3 * (xfc[k] + edge.center[k]);
2899
2900 cs_xdef_cw_eval_at_xyz_by_analytic(cm,
2901 1, (const cs_real_t *)xg,
2902 time_eval,
2903 def->context,
2904 (cs_real_t *)adv_xg);
2905
2906 fluxes[e] += sefc.meas * _dp3(adv_xg, sefc.unitv);
2907
2908 } /* Loop on face edges */
2909
2910 } /* Loop on cell faces */
2911 }
2912 break; /* quadrature NONE, BARY or BARY_SUBDIV */
2913
2914 case CS_QUADRATURE_HIGHER:
2915 {
2916 /* 4 Gauss points by triangle spanned by x_c, x_f, x_e */
2917 cs_real_t w[4], eval[12];
2918 cs_real_3_t gpts[4];
2919
2920 /* Loop on cell faces */
2921 for (short int f = 0; f < cm->n_fc; f++) {
2922
2923 const cs_quant_t face = cm->face[f];
2924
2925 for (short int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
2926
2927 const short int e = cm->f2e_ids[i];
2928 const cs_quant_t edge = cm->edge[e];
2929 const cs_nvec3_t sefc = cm->sefc[i];
2930
2931 cs_quadrature_tria_4pts(edge.center, face.center, cm->xc,
2932 sefc.meas,
2933 gpts, w);
2934
2935 cs_xdef_cw_eval_at_xyz_by_analytic(cm,
2936 4, (const cs_real_t *)gpts,
2937 time_eval,
2938 def->context,
2939 eval);
2940
2941 cs_real_t add = 0;
2942 for (int p = 0; p < 4; p++)
2943 add += w[p] * _dp3(eval + 3*p, sefc.unitv);
2944 fluxes[e] += add;
2945
2946 } /* Loop on face edges */
2947
2948 } /* Loop on cell faces */
2949 }
2950 break; /* CS_QUADRATURE_HIGHER */
2951
2952 case CS_QUADRATURE_HIGHEST:
2953 {
2954 /* 7 Gauss points by triangle spanned by x_c, x_f, x_e */
2955 cs_real_t w[7], eval[21];
2956 cs_real_3_t gpts[7];
2957
2958 /* Loop on cell faces */
2959 for (short int f = 0; f < cm->n_fc; f++) {
2960
2961 const cs_quant_t face = cm->face[f];
2962
2963 for (short int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
2964
2965 const short int e = cm->f2e_ids[i];
2966 const cs_quant_t edge = cm->edge[e];
2967 const cs_nvec3_t sefc = cm->sefc[i];
2968
2969 cs_quadrature_tria_7pts(edge.center, face.center, cm->xc,
2970 sefc.meas,
2971 gpts, w);
2972
2973 cs_xdef_cw_eval_at_xyz_by_analytic(cm,
2974 7, (const cs_real_t *)gpts,
2975 time_eval,
2976 def->context,
2977 eval);
2978
2979 cs_real_t add = 0;
2980 for (int p = 0; p < 7; p++)
2981 add += w[p] * _dp3(eval + 3*p, sefc.unitv);
2982 fluxes[e] += add;
2983
2984 } /* Loop on face edges */
2985
2986 } /* Loop on cell faces */
2987 }
2988 break; /* CS_QUADRATURE_HIGHEST */
2989
2990 default:
2991 bft_error(__FILE__, __LINE__, 0,
2992 " %s: Invalid type of quadrature.", __func__);
2993 break;
2994
2995 } /* Switch type of quadrature */
2996
2997 } /* Definition by analytic function */
2998 break;
2999
3000 case CS_XDEF_BY_ARRAY:
3001 {
3002 cs_xdef_array_context_t *ac = (cs_xdef_array_context_t *)def->context;
3003
3004 /* Test if location has at least the pattern of the reference support */
3005 if (cs_flag_test(ac->loc, cs_flag_dual_face_byc)) {
3006
3007 /* Sanity check */
3008 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PE));
3009
3010 const cs_real_t *flux_array =
3011 ac->values + cs_cdo_connect->c2e->idx[cm->c_id];
3012 for (short int e = 0; e < cm->n_ec; e++)
3013 fluxes[e] = flux_array[e];
3014
3015 }
3016 else if (cs_flag_test(ac->loc, cs_flag_primal_face)) {
3017
3018 const cs_real_t *i_mflx = ac->values;
3019 const cs_real_t *b_mflx = ac->values + cm->bface_shift;
3020
3021 /* Sanity check */
3022 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
3023 CS_FLAG_COMP_DFQ | CS_FLAG_COMP_PEQ));
3024
3025 cs_real_t cell_mflx[3] = {0., 0., 0.};
3026 cs_reco_cw_cell_vect_from_face_dofs(cm, i_mflx, b_mflx, cell_mflx);
3027
3028 /* Loop on cell edges */
3029 for (short int e = 0; e < cm->n_ec; e++)
3030 fluxes[e] = cm->dface[e].meas * _dp3(cell_mflx, cm->dface[e].unitv);
3031
3032 }
3033 else if (cs_flag_test(ac->loc, cs_flag_primal_cell)) {
3034
3035 /* Retrieve the advection field:
3036 Switch to a cs_nvec3_t representation */
3037 cs_nvec3_t nv;
3038 cs_nvec3(ac->values + 3*cm->c_id, &nv);
3039
3040 /* Sanity check */
3041 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_DFQ | CS_FLAG_COMP_PEQ));
3042
3043 /* Loop on cell edges */
3044 for (short int e = 0; e < cm->n_ec; e++)
3045 fluxes[e] = nv.meas*cm->dface[e].meas*_dp3(nv.unitv,
3046 cm->dface[e].unitv);
3047
3048 }
3049 else
3050 bft_error(__FILE__, __LINE__, 0,
3051 " Invalid support for evaluating the advection field %s"
3052 " at the cell center of cell %ld.", adv->name,
3053 (long)cm->c_id);
3054 }
3055 break;
3056
3057 case CS_XDEF_BY_FIELD:
3058 {
3059 cs_field_t *f = (cs_field_t *)def->context;
3060 const cs_mesh_location_type_t loc_type =
3061 cs_mesh_location_get_type(f->location_id);
3062
3063 switch(loc_type) {
3064
3065 case CS_MESH_LOCATION_CELLS:
3066 {
3067 /* Retrieve the advection field:
3068 Switch to a cs_nvec3_t representation */
3069 cs_nvec3_t nv;
3070 cs_nvec3(f->val + 3*cm->c_id, &nv);
3071
3072 /* Sanity check */
3073 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_DFQ | CS_FLAG_COMP_PEQ));
3074
3075 /* Loop on cell edges */
3076 for (short int e = 0; e < cm->n_ec; e++)
3077 fluxes[e] = nv.meas*cm->dface[e].meas*_dp3(nv.unitv,
3078 cm->dface[e].unitv);
3079 }
3080 break;
3081
3082 case CS_MESH_LOCATION_INTERIOR_FACES:
3083 {
3084 assert(adv->bdy_field_id > -1);
3085 cs_field_t *b_mflx_fld = cs_field_by_id(adv->bdy_field_id);
3086 assert(b_mflx_fld != NULL);
3087
3088 const cs_real_t *b_mflx = b_mflx_fld->val;
3089 const cs_real_t *i_mflx = f->val;
3090
3091 /* Sanity check */
3092 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
3093 CS_FLAG_COMP_DFQ | CS_FLAG_COMP_PEQ));
3094
3095 cs_real_t cell_mflx[3] = {0., 0., 0.};
3096 cs_reco_cw_cell_vect_from_face_dofs(cm, i_mflx, b_mflx, cell_mflx);
3097
3098 /* Loop on cell edges */
3099 for (short int e = 0; e < cm->n_ec; e++)
3100 fluxes[e] = cm->dface[e].meas * _dp3(cell_mflx, cm->dface[e].unitv);
3101 }
3102 break;
3103
3104 default:
3105 bft_error(__FILE__, __LINE__, 0, "%s: Invalid location type.\nTODO.",
3106 __func__);
3107
3108 }
3109
3110 }
3111 break;
3112
3113 default:
3114 bft_error(__FILE__, __LINE__, 0, "%s: Incompatible type of definition.",
3115 __func__);
3116 break;
3117
3118 } /* def_type */
3119
3120 }
3121
3122 /*----------------------------------------------------------------------------*/
3123 /*!
3124 * \brief For each cs_adv_field_t structures, update the values of the
3125 * related field(s)
3126 *
3127 * \param[in] t_eval physical time at which one evaluates the term
3128 * \param[in] cur2prev true or false
3129 */
3130 /*----------------------------------------------------------------------------*/
3131
3132 void
cs_advection_field_update(cs_real_t t_eval,bool cur2prev)3133 cs_advection_field_update(cs_real_t t_eval,
3134 bool cur2prev)
3135 {
3136 for (int i = 0; i < _n_adv_fields; i++) {
3137
3138 cs_adv_field_t *adv = _adv_fields[i];
3139
3140 /* Sanity checks */
3141 assert(adv != NULL);
3142
3143 if (t_eval > 0 && (adv->status & CS_ADVECTION_FIELD_STEADY))
3144 continue;
3145
3146 /* GWF and NAVSTO categories of advection fields are updated elsewhere
3147 except if there is a field defined at vertices */
3148
3149 if ((adv->status & CS_ADVECTION_FIELD_USER) ||
3150 (adv->status & CS_ADVECTION_FIELD_LEGACY_FV)) {
3151
3152 /* Field stored at cell centers */
3153 assert(adv->cell_field_id > -1);
3154 cs_field_t *cfld = cs_field_by_id(adv->cell_field_id);
3155
3156 /* Copy current field values to previous values */
3157 if (cur2prev)
3158 cs_field_current_to_previous(cfld);
3159
3160 /* Set the new values */
3161 cs_advection_field_in_cells(adv, t_eval, cfld->val);
3162
3163 /* Set the new values (only in the case of a user-defined advection
3164 field */
3165 if ((adv->status & CS_ADVECTION_FIELD_USER) && adv->bdy_field_id > -1) {
3166
3167 /* Field storing the boundary normal flux */
3168 cs_field_t *bfld = cs_field_by_id(adv->bdy_field_id);
3169
3170 /* Copy current field values to previous values */
3171 if (cur2prev)
3172 cs_field_current_to_previous(bfld);
3173
3174 cs_advection_field_across_boundary(adv, t_eval, bfld->val);
3175
3176 }
3177
3178 }
3179
3180 if (adv->vtx_field_id > -1) { /* Field stored at vertices */
3181
3182 cs_field_t *vfld = cs_field_by_id(adv->vtx_field_id);
3183
3184 /* Copy current field values to previous values */
3185 if (cur2prev)
3186 cs_field_current_to_previous(vfld);
3187
3188 /* Set the new values */
3189 cs_advection_field_at_vertices(adv, t_eval, vfld->val);
3190
3191 }
3192
3193 } /* Loop on advection fields */
3194
3195 }
3196
3197 /*----------------------------------------------------------------------------*/
3198 /*!
3199 * \brief Compute the Peclet number in each cell
3200 *
3201 * \param[in] adv pointer to the advection field struct.
3202 * \param[in] diff pointer to the diffusion property struct.
3203 * \param[in] t_eval time at which one evaluates the advection field
3204 * \param[in, out] peclet pointer to an array storing the Peclet number
3205 */
3206 /*----------------------------------------------------------------------------*/
3207
3208 void
cs_advection_get_peclet(const cs_adv_field_t * adv,const cs_property_t * diff,cs_real_t t_eval,cs_real_t peclet[])3209 cs_advection_get_peclet(const cs_adv_field_t *adv,
3210 const cs_property_t *diff,
3211 cs_real_t t_eval,
3212 cs_real_t peclet[])
3213 {
3214 /* Sanity checks */
3215 assert(adv != NULL);
3216 assert(diff != NULL);
3217 assert(peclet != NULL);
3218
3219 cs_real_t ptymat[3][3];
3220 cs_real_3_t ptydir;
3221 cs_nvec3_t adv_c;
3222
3223 const bool pty_uniform = cs_property_is_uniform(diff);
3224 const cs_cdo_quantities_t *cdoq = cs_cdo_quant;
3225
3226 /* Get the value of the material property at the first cell center */
3227 if (pty_uniform)
3228 cs_property_get_cell_tensor(0, t_eval, diff, false, ptymat);
3229
3230 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
3231
3232 /* Get the value of the material property at the cell center */
3233 if (!pty_uniform)
3234 cs_property_get_cell_tensor(c_id, t_eval, diff, false, ptymat);
3235
3236 const cs_real_t hc = cbrt(cdoq->cell_vol[c_id]);
3237
3238 cs_advection_field_get_cell_vector(c_id, adv, &adv_c);
3239 cs_math_33_3_product((const cs_real_t (*)[3])ptymat, adv_c.unitv, ptydir);
3240 peclet[c_id] = hc * adv_c.meas / _dp3(adv_c.unitv, ptydir);
3241
3242 } /* Loop on cells */
3243
3244 }
3245
3246 /*----------------------------------------------------------------------------*/
3247 /*!
3248 * \brief Compute the Courant number in each cell
3249 *
3250 * \param[in] adv pointer to the advection field struct.
3251 * \param[in] dt_cur current time step
3252 * \param[in, out] courant pointer to an array storing the Courant number
3253 */
3254 /*----------------------------------------------------------------------------*/
3255
3256 void
cs_advection_get_courant(const cs_adv_field_t * adv,cs_real_t dt_cur,cs_real_t courant[])3257 cs_advection_get_courant(const cs_adv_field_t *adv,
3258 cs_real_t dt_cur,
3259 cs_real_t courant[])
3260 {
3261 assert(courant != NULL); /* Sanity check */
3262 const cs_cdo_quantities_t *cdoq = cs_cdo_quant;
3263 const cs_adjacency_t *c2f = cs_cdo_connect->c2f;
3264
3265 assert(adv->cell_field_id > -1); /* field should be defined at cell centers */
3266
3267 const cs_field_t *fld = cs_field_by_id(adv->cell_field_id);
3268
3269 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
3270 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
3271
3272 const cs_real_t *vel_c = fld->val + 3*c_id;
3273 const cs_real_t ovol_c = 1./cdoq->cell_vol[c_id];
3274
3275 cs_real_t _courant = 0.;
3276 for (cs_lnum_t i = c2f->idx[c_id]; i < c2f->idx[c_id+1]; i++) {
3277
3278 const cs_real_t *f_area = cs_quant_get_face_vector_area(c2f->ids[i],
3279 cdoq);
3280 _courant = fmax(_courant, fabs(_dp3(f_area, vel_c)) * ovol_c);
3281
3282 }
3283 courant[c_id] = _courant * dt_cur;
3284
3285 } /* Loop on cells */
3286
3287 }
3288
3289 /*----------------------------------------------------------------------------*/
3290 /*!
3291 * \brief Compute the divergence of the advection field at vertices
3292 * Useful for CDO Vertex-based schemes
3293 *
3294 * \param[in] adv pointer to the advection field struct.
3295 * \param[in] t_eval time at which one evaluates the advection field
3296 *
3297 * \return a pointer to an array storing the result
3298 */
3299 /*----------------------------------------------------------------------------*/
3300
3301 cs_real_t *
cs_advection_field_divergence_at_vertices(const cs_adv_field_t * adv,cs_real_t t_eval)3302 cs_advection_field_divergence_at_vertices(const cs_adv_field_t *adv,
3303 cs_real_t t_eval)
3304 {
3305 CS_UNUSED(t_eval); /* Useful in case of analytic definition */
3306
3307 cs_real_t *divergence = NULL;
3308
3309 if (adv == NULL)
3310 return divergence;
3311
3312 const cs_cdo_quantities_t *cdoq = cs_cdo_quant;
3313 const cs_cdo_connect_t *connect = cs_cdo_connect;
3314 const cs_adjacency_t *f2e = connect->f2e;
3315 const cs_adjacency_t *e2v = connect->e2v;
3316
3317 BFT_MALLOC(divergence, cdoq->n_vertices, cs_real_t);
3318 memset(divergence, 0, sizeof(cs_real_t)*cdoq->n_vertices);
3319
3320 { /* Volume part */
3321 const cs_xdef_t *def = adv->definition;
3322
3323 switch (def->type) {
3324
3325 case CS_XDEF_BY_ARRAY:
3326 {
3327 cs_xdef_array_context_t *ac = (cs_xdef_array_context_t *)def->context;
3328
3329 if (cs_flag_test(ac->loc, cs_flag_dual_face_byc)) {
3330
3331 const cs_adjacency_t *c2e = connect->c2e;
3332
3333 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
3334 for (cs_lnum_t j = c2e->idx[c_id]; j < c2e->idx[c_id+1]; j++) {
3335
3336 const cs_lnum_t e_id = c2e->ids[j];
3337 const cs_real_t flx = ac->values[j];
3338 const cs_lnum_t eshift = 2*e_id;
3339 const cs_lnum_t v0 = e2v->ids[eshift];
3340 const cs_lnum_t v1 = e2v->ids[eshift+1];
3341 const short int sgn = e2v->sgn[eshift]; /* Sign for v0 */
3342
3343 divergence[v0] += -sgn*flx;
3344 divergence[v1] += sgn*flx;
3345
3346 } /* Loop on cell edges */
3347 } /* Loop on cells */
3348
3349 }
3350 else
3351 bft_error(__FILE__, __LINE__, 0,
3352 " %s: Invalid location for the array.", __func__);
3353
3354 }
3355 break;
3356
3357 default:
3358 bft_error(__FILE__, __LINE__, 0, " %s: Invalid case.", __func__);
3359
3360 } /* End of switch */
3361
3362 } /* Volume part */
3363
3364 /* Boundary part */
3365 if (adv->n_bdy_flux_defs > 0) {
3366
3367 for (int def_id = 0; def_id < adv->n_bdy_flux_defs; def_id++) {
3368
3369 const cs_xdef_t *def = adv->bdy_flux_defs[def_id];
3370 const cs_zone_t *z = cs_boundary_zone_by_id(def->z_id);
3371 assert(def->support == CS_XDEF_SUPPORT_BOUNDARY);
3372
3373 switch (def->type) {
3374
3375 case CS_XDEF_BY_VALUE:
3376 {
3377 const cs_real_t *constant_val = (cs_real_t *)def->context;
3378
3379 for (cs_lnum_t id = 0; id < z->n_elts; id++) {
3380
3381 const cs_lnum_t bf_id = (z->elt_ids == NULL) ? id: z->elt_ids[id];
3382
3383 _fill_uniform_boundary_flux(cdoq, f2e, e2v, bf_id, constant_val[0],
3384 divergence);
3385
3386 } /* Loop on boundary faces */
3387 }
3388 break; /* by value */
3389
3390 case CS_XDEF_BY_ARRAY:
3391 {
3392 const cs_xdef_array_context_t *ac =
3393 (cs_xdef_array_context_t *)def->context;
3394 const cs_real_t *val = ac->values;
3395
3396 assert(ac->stride == 1);
3397 assert(z->id == 0);
3398
3399 if (cs_flag_test(ac->loc, cs_flag_primal_face)) {
3400
3401 for (cs_lnum_t bf_id = 0; bf_id < cdoq->n_b_faces; bf_id++)
3402 _fill_uniform_boundary_flux(cdoq, f2e, e2v, bf_id, val[bf_id],
3403 divergence);
3404
3405 }
3406 else if (cs_flag_test(ac->loc, cs_flag_dual_closure_byf)) {
3407
3408 const cs_adjacency_t *bf2v = connect->bf2v;
3409 const cs_lnum_t *idx = bf2v->idx;
3410 assert(idx == ac->index);
3411
3412 for (cs_lnum_t bf_id = 0; bf_id < cdoq->n_b_faces; bf_id++) {
3413 for (cs_lnum_t i = idx[bf_id]; i < idx[bf_id+1]; i++) {
3414 const cs_lnum_t v_id = bf2v->ids[i];
3415 divergence[v_id] += val[i];
3416 }
3417 }
3418
3419 }
3420 else
3421 bft_error(__FILE__, __LINE__, 0, " %s: Invalid case.", __func__);
3422
3423 }
3424 break; /* by_array */
3425
3426 default:
3427 bft_error(__FILE__, __LINE__, 0, " %s: Invalid case", __func__);
3428 break;
3429
3430 } /* End of switch on the type of definition */
3431 } /* Loop on definitions */
3432
3433 }
3434 else {
3435
3436 /* Handle the boundary flux */
3437 const cs_field_t *bflx =
3438 cs_advection_field_get_field(adv, CS_MESH_LOCATION_BOUNDARY_FACES);
3439
3440 for (cs_lnum_t bf_id = 0; bf_id < cdoq->n_b_faces; bf_id++) {
3441
3442 const cs_real_t face_flx = bflx->val[bf_id];
3443 const cs_real_t invsurf = 1./cdoq->b_face_surf[bf_id];
3444 const cs_lnum_t f_id = cdoq->n_i_faces + bf_id;
3445
3446 for (cs_lnum_t i = f2e->idx[f_id]; i < f2e->idx[f_id+1]; i++) {
3447
3448 const cs_lnum_t e_id = f2e->ids[i];
3449 const cs_lnum_t eshift = 2*e_id;
3450 const cs_lnum_t v0 = e2v->ids[eshift];
3451 const cs_lnum_t v1 = e2v->ids[eshift+1];
3452 const double tef = cs_math_surftri(cdoq->vtx_coord + 3*v0,
3453 cdoq->vtx_coord + 3*v1,
3454 cdoq->b_face_center + 3*bf_id);
3455
3456 const double weighted_flux = 0.5 * tef * invsurf * face_flx;
3457
3458 divergence[v0] += weighted_flux;
3459 divergence[v1] += weighted_flux;
3460
3461 } /* Loop on face edges */
3462
3463 } /* Loop on boundary faces */
3464
3465 } /* Boundary part */
3466
3467 #if defined(HAVE_MPI) /* Synchronization if needed */
3468 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL)
3469 cs_interface_set_sum(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
3470 cdoq->n_vertices,
3471 1, /* stride */
3472 false, /* interlace (not useful here) */
3473 CS_REAL_TYPE,
3474 divergence);
3475 #endif /* MPI */
3476
3477 return divergence;
3478 }
3479
3480 /*----------------------------------------------------------------------------*/
3481
3482 #undef _dp3
3483
3484 END_C_DECLS
3485