1 /*============================================================================
2 * Field based algebraic operators.
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------*/
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *----------------------------------------------------------------------------*/
34
35 #include <assert.h>
36 #include <math.h>
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <string.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include "bft_mem.h"
46 #include "bft_error.h"
47 #include "bft_printf.h"
48
49 #include "cs_field.h"
50 #include "cs_field_default.h"
51 #include "cs_gradient.h"
52 #include "cs_halo.h"
53 #include "cs_halo_perio.h"
54 #include "cs_mesh.h"
55 #include "cs_log.h"
56 #include "cs_map.h"
57 #include "cs_parameters.h"
58 #include "cs_parall.h"
59 #include "cs_mesh.h"
60 #include "cs_mesh_adjacencies.h"
61 #include "cs_mesh_location.h"
62 #include "cs_mesh_quantities.h"
63 #include "cs_internal_coupling.h"
64
65 /*----------------------------------------------------------------------------
66 * Header for the current file
67 *----------------------------------------------------------------------------*/
68
69 #include "cs_field_operator.h"
70
71 /*----------------------------------------------------------------------------*/
72
73 BEGIN_C_DECLS
74
75 /*=============================================================================
76 * Additional doxygen documentation
77 *============================================================================*/
78
79 /*!
80 \file cs_field_operator.c
81 Field based algebraic operators.
82
83 \enum cs_field_interpolate_t
84
85 \brief Field interpolation modes
86
87 \var CS_FIELD_INTERPOLATE_MEAN
88 Mean element value (P0 interpolation)
89 \var CS_FIELD_INTERPOLATE_GRADIENT
90 Mean element value + gradient correction (pseudo-P1)
91 */
92
93 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
94
95 /*=============================================================================
96 * Macro definitions
97 *============================================================================*/
98
99 /*============================================================================
100 * Type definitions
101 *============================================================================*/
102
103 /*============================================================================
104 * Static global variables
105 *============================================================================*/
106
107 /*============================================================================
108 * Prototypes for functions intended for use only by Fortran wrappers.
109 * (descriptions follow, with function bodies).
110 *============================================================================*/
111
112 void
113 cs_f_field_gradient_scalar(int f_id,
114 int use_previous_t,
115 int imrgr0,
116 int inc,
117 int recompute_cocg,
118 cs_real_3_t *restrict grad);
119
120 void
121 cs_f_field_gradient_potential(int f_id,
122 int use_previous_t,
123 int imrgr0,
124 int inc,
125 int recompute_cocg,
126 int hyd_p_flag,
127 cs_real_3_t f_ext[],
128 cs_real_3_t *restrict grad);
129
130 void
131 cs_f_field_gradient_vector(int f_id,
132 int use_previous_t,
133 int imrgr0,
134 int inc,
135 cs_real_33_t *restrict grad);
136
137 void
138 cs_f_field_gradient_tensor(int f_id,
139 int use_previous_t,
140 int imrgr0,
141 int inc,
142 cs_real_63_t *restrict grad);
143
144 void
145 cs_f_field_set_volume_average(int f_id,
146 cs_real_t va);
147
148 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
149
150 /*============================================================================
151 * Private function definitions
152 *============================================================================*/
153
154 /*----------------------------------------------------------------------------
155 * Interpolate field values at a given set of points using P0 interpolation.
156 *
157 * parameters:
158 * f <-- pointer to field
159 * n_points <-- number of points at which interpolation
160 * is required
161 * point_location <-- location of points in mesh elements
162 * (based on the field location)
163 * val --> interpolated values
164 *----------------------------------------------------------------------------*/
165
166 static void
_field_interpolate_by_mean(const cs_field_t * f,cs_lnum_t n_points,const cs_lnum_t point_location[],cs_real_t * val)167 _field_interpolate_by_mean(const cs_field_t *f,
168 cs_lnum_t n_points,
169 const cs_lnum_t point_location[],
170 cs_real_t *val)
171 {
172 for (cs_lnum_t i = 0; i < n_points; i++) {
173
174 cs_lnum_t cell_id = point_location[i];
175
176 for (cs_lnum_t j = 0; j < f->dim; j++)
177 val[i*f->dim + j] = f->val[cell_id*f->dim + j];
178
179 }
180 }
181
182 /*----------------------------------------------------------------------------
183 * Interpolate field values at a given set of points using gradient-corrected
184 * interpolation.
185 *
186 * parameters:
187 * f <-- pointer to field
188 * n_points <-- number of points at which interpolation
189 * is required
190 * point_location <-- location of points in mesh elements
191 * (based on the field location)
192 * point_coords <-- point coordinates
193 * val --> interpolated values
194 *----------------------------------------------------------------------------*/
195
196 static void
_field_interpolate_by_gradient(const cs_field_t * f,cs_lnum_t n_points,const cs_lnum_t point_location[],const cs_real_3_t point_coords[],cs_real_t * val)197 _field_interpolate_by_gradient(const cs_field_t *f,
198 cs_lnum_t n_points,
199 const cs_lnum_t point_location[],
200 const cs_real_3_t point_coords[],
201 cs_real_t *val)
202 {
203 const cs_lnum_t dim = f->dim;
204 const cs_lnum_t n_cells_ext = cs_glob_mesh->n_cells_with_ghosts;
205 const cs_real_3_t *cell_cen
206 = (const cs_real_3_t *)(cs_glob_mesh_quantities->cell_cen);
207
208 /* Currently possible only for fields on cell location */
209
210 if (f->location_id != CS_MESH_LOCATION_CELLS)
211 bft_error(__FILE__, __LINE__, 0,
212 _("Field gradient interpolation for field %s :\n"
213 " not implemented for fields on location %s."),
214 f->name, cs_mesh_location_type_name[f->location_id]);
215
216 /* Compute field cell gradient */
217
218 cs_real_t *grad;
219 BFT_MALLOC(grad, 3*dim*n_cells_ext, cs_real_t);
220
221 if (dim == 1)
222 cs_field_gradient_scalar(f,
223 true, /* use_previous_t */
224 1, /* inc */
225 true, /* recompute_cocg */
226 (cs_real_3_t *)grad);
227 else if (dim == 3)
228 cs_field_gradient_vector(f,
229 true, /* use_previous_t */
230 1, /* inc */
231 (cs_real_33_t *)grad);
232
233 else
234 bft_error(__FILE__, __LINE__, 0,
235 _("Field gradient interpolation for field %s of dimension %d:\n"
236 " not implemented."),
237 f->name, (int)dim);
238
239 /* Now interpolated values */
240
241 for (cs_lnum_t i = 0; i < n_points; i++) {
242
243 cs_lnum_t cell_id = point_location[i];
244
245 cs_real_3_t d = {point_coords[i][0] - cell_cen[cell_id][0],
246 point_coords[i][1] - cell_cen[cell_id][1],
247 point_coords[i][2] - cell_cen[cell_id][2]};
248
249 for (cs_lnum_t j = 0; j < f->dim; j++) {
250 cs_lnum_t k = (cell_id*dim + j)*3;
251 val[i*dim + j] = f->val[cell_id*dim + j]
252 + d[0] * grad[k]
253 + d[1] * grad[k+1]
254 + d[2] * grad[k+2];
255 }
256
257 }
258
259 BFT_FREE(grad);
260 }
261
262 /*----------------------------------------------------------------------------
263 * For each mesh cell this function finds the local extrema of a
264 * scalar field.
265 *
266 * parameters:
267 * pvar <-- scalar values
268 * halo_type <-- halo type
269 * local_max --> local maximum value
270 * local_min --> local minimum value
271 *----------------------------------------------------------------------------*/
272
273 static void
_local_extrema_scalar(const cs_real_t * restrict pvar,cs_halo_type_t halo_type,cs_real_t * local_max,cs_real_t * local_min)274 _local_extrema_scalar(const cs_real_t *restrict pvar,
275 cs_halo_type_t halo_type,
276 cs_real_t *local_max,
277 cs_real_t *local_min)
278 {
279 const cs_mesh_t *m = cs_glob_mesh;
280 const cs_lnum_t n_cells = m->n_cells;
281 const cs_lnum_t n_vertices = m->n_vertices;
282
283 const cs_adjacency_t *c2v = cs_mesh_adjacencies_cell_vertices();
284 const cs_lnum_t *c2v_idx = c2v->idx;
285 const cs_lnum_t *c2v_ids = c2v->ids;
286
287 # pragma omp parallel for if (n_cells > CS_THR_MIN)
288 for (cs_lnum_t ii = 0; ii < n_cells; ii++) {
289 local_max[ii] = pvar[ii];
290 local_min[ii] = pvar[ii];
291 }
292
293 cs_real_t *v_min, *v_max;
294 BFT_MALLOC(v_min, n_vertices, cs_real_t);
295 BFT_MALLOC(v_max, n_vertices, cs_real_t);
296
297 # pragma omp parallel for if (n_vertices > CS_THR_MIN)
298 for (cs_lnum_t ii = 0; ii < n_vertices; ii++) {
299 v_max[ii] = -HUGE_VAL;
300 v_min[ii] = HUGE_VAL;
301 }
302
303 /* Scatter min/max values to vertices */
304
305 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
306 cs_lnum_t s_id = c2v_idx[c_id];
307 cs_lnum_t e_id = c2v_idx[c_id+1];
308 cs_real_t _c_var = pvar[c_id];
309 for (cs_lnum_t j = s_id; j < e_id; j++) {
310 cs_lnum_t v_id = c2v_ids[j];
311 if (_c_var < v_min[v_id])
312 v_min[v_id] = _c_var;
313 if (_c_var > v_max[v_id])
314 v_max[v_id] = _c_var;
315 }
316 }
317
318 if (m->vtx_interfaces != NULL) {
319 cs_interface_set_min(m->vtx_interfaces,
320 m->n_vertices,
321 1,
322 true,
323 CS_REAL_TYPE,
324 v_min);
325 cs_interface_set_max(m->vtx_interfaces,
326 m->n_vertices,
327 1,
328 true,
329 CS_REAL_TYPE,
330 v_max);
331 }
332
333 /* Gather min/max values from vertices */
334
335 # pragma omp parallel for if (n_cells > CS_THR_MIN)
336 for (cs_lnum_t c_id = 0; c_id < n_cells; c_id++) {
337 cs_lnum_t s_id = c2v_idx[c_id];
338 cs_lnum_t e_id = c2v_idx[c_id+1];
339 for (cs_lnum_t j = s_id; j < e_id; j++) {
340 cs_lnum_t v_id = c2v_ids[j];
341 if (v_min[v_id] < local_min[c_id])
342 local_min[c_id] = v_min[v_id];
343 if (v_max[v_id] > local_max[c_id])
344 local_max[c_id] = v_max[v_id];
345 }
346 }
347
348 /* Free memory */
349 BFT_FREE(v_min);
350 BFT_FREE(v_max);
351
352 /* Synchronisation */
353 if (m->halo != NULL) {
354 cs_halo_sync_var(m->halo, halo_type, local_min);
355 cs_halo_sync_var(m->halo, halo_type, local_max);
356 }
357 }
358
359 /*============================================================================
360 * Fortran wrapper function definitions
361 *============================================================================*/
362
363 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
364
365 /*----------------------------------------------------------------------------
366 * Compute cell gradient of scalar field or component of vector or
367 * tensor field.
368 *
369 * parameters:
370 * f_id <-- field id
371 * use_previous_t <-- should we use values from the previous time step ?
372 * imrgr0 <-- ignored
373 * inc <-- if 0, solve on increment; 1 otherwise
374 * recompute_cocg <-- should COCG FV quantities be recomputed ?
375 * grad --> gradient
376 *----------------------------------------------------------------------------*/
377
378 void
cs_f_field_gradient_scalar(int f_id,int use_previous_t,int imrgr0,int inc,int recompute_cocg,cs_real_3_t * restrict grad)379 cs_f_field_gradient_scalar(int f_id,
380 int use_previous_t,
381 int imrgr0,
382 int inc,
383 int recompute_cocg,
384 cs_real_3_t *restrict grad)
385 {
386 CS_UNUSED(imrgr0);
387
388 bool _use_previous_t = use_previous_t ? true : false;
389 bool _recompute_cocg = recompute_cocg ? true : false;
390
391 const cs_field_t *f = cs_field_by_id(f_id);
392
393 cs_field_gradient_scalar(f,
394 _use_previous_t,
395 inc,
396 _recompute_cocg,
397 grad);
398 }
399
400 /*----------------------------------------------------------------------------
401 * Compute cell gradient of scalar field or component of vector or
402 * tensor field.
403 *
404 * parameters:
405 * f_id <-- field id
406 * use_previous_t <-- should we use values from the previous time step ?
407 * imrgr0 <-- ignored
408 * halo_type <-- halo type
409 * inc <-- if 0, solve on increment; 1 otherwise
410 * recompute_cocg <-- should COCG FV quantities be recomputed ?
411 * hyd_p_flag <-- flag for hydrostatic pressure
412 * f_ext <-- exterior force generating the hydrostatic pressure
413 * grad --> gradient
414 *----------------------------------------------------------------------------*/
415
416 void
cs_f_field_gradient_potential(int f_id,int use_previous_t,int imrgr0,int inc,int recompute_cocg,int hyd_p_flag,cs_real_3_t f_ext[],cs_real_3_t * restrict grad)417 cs_f_field_gradient_potential(int f_id,
418 int use_previous_t,
419 int imrgr0,
420 int inc,
421 int recompute_cocg,
422 int hyd_p_flag,
423 cs_real_3_t f_ext[],
424 cs_real_3_t *restrict grad)
425 {
426 CS_UNUSED(imrgr0);
427
428 bool _use_previous_t = use_previous_t ? true : false;
429 bool _recompute_cocg = recompute_cocg ? true : false;
430
431 const cs_field_t *f = cs_field_by_id(f_id);
432
433 cs_field_gradient_potential(f,
434 _use_previous_t,
435 inc,
436 _recompute_cocg,
437 hyd_p_flag,
438 f_ext,
439 grad);
440 }
441
442 /*----------------------------------------------------------------------------
443 * Compute cell gradient of scalar field or component of vector or
444 * tensor field.
445 *
446 * parameters:
447 * f_id <-- field id
448 * use_previous_t <-- should we use values from the previous time step ?
449 * imrgr0 <-- ignored
450 * inc <-- if 0, solve on increment; 1 otherwise
451 * recompute_cocg <-- should COCG FV quantities be recomputed ?
452 * grad --> gradient
453 *----------------------------------------------------------------------------*/
454
455 void
cs_f_field_gradient_vector(int f_id,int use_previous_t,int imrgr0,int inc,cs_real_33_t * restrict grad)456 cs_f_field_gradient_vector(int f_id,
457 int use_previous_t,
458 int imrgr0,
459 int inc,
460 cs_real_33_t *restrict grad)
461 {
462 CS_UNUSED(imrgr0);
463
464 bool _use_previous_t = use_previous_t ? true : false;
465
466 const cs_field_t *f = cs_field_by_id(f_id);
467
468 cs_field_gradient_vector(f,
469 _use_previous_t,
470 inc,
471 grad);
472 }
473
474 /*----------------------------------------------------------------------------
475 * Compute cell gradient of scalar field or component of vector or
476 * tensor field.
477 *
478 * parameters:
479 * f_id <-- field id
480 * use_previous_t <-- should we use values from the previous time step ?
481 * imrgr0 <-- ignored
482 * inc <-- if 0, solve on increment; 1 otherwise
483 * recompute_cocg <-- should COCG FV quantities be recomputed ?
484 * grad --> gradient
485 *----------------------------------------------------------------------------*/
486
487 void
cs_f_field_gradient_tensor(int f_id,int use_previous_t,int imrgr0,int inc,cs_real_63_t * restrict grad)488 cs_f_field_gradient_tensor(int f_id,
489 int use_previous_t,
490 int imrgr0,
491 int inc,
492 cs_real_63_t *restrict grad)
493 {
494 CS_UNUSED(imrgr0);
495
496 bool _use_previous_t = use_previous_t ? true : false;
497
498 const cs_field_t *f = cs_field_by_id(f_id);
499
500 cs_field_gradient_tensor(f,
501 _use_previous_t,
502 inc,
503 grad);
504 }
505
506 /*----------------------------------------------------------------------------
507 * Shift field values in order to set its spatial average to a given value.
508 *
509 * parameters:
510 * f_id <-- field id
511 * va <-- real value of volume average to be set
512 *----------------------------------------------------------------------------*/
513
514 void
cs_f_field_set_volume_average(int f_id,cs_real_t va)515 cs_f_field_set_volume_average(int f_id,
516 cs_real_t va)
517 {
518 cs_field_t *f = cs_field_by_id(f_id);
519
520 cs_field_set_volume_average(f,
521 va);
522 }
523
524 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
525
526 /*=============================================================================
527 * Public function definitions
528 *============================================================================*/
529
530 /*----------------------------------------------------------------------------*/
531 /*!
532 * \brief Compute cell gradient of scalar field or component of vector or
533 * tensor field.
534 *
535 * \param[in] f pointer to field
536 * \param[in] use_previous_t should we use values from the previous
537 * time step ?
538 * \param[in] inc if 0, solve on increment; 1 otherwise
539 * \param[in] recompute_cocg should COCG FV quantities be recomputed ?
540 * \param[out] grad gradient
541 */
542 /*----------------------------------------------------------------------------*/
543
544 void
cs_field_gradient_scalar(const cs_field_t * f,bool use_previous_t,int inc,bool recompute_cocg,cs_real_3_t * restrict grad)545 cs_field_gradient_scalar(const cs_field_t *f,
546 bool use_previous_t,
547 int inc,
548 bool recompute_cocg,
549 cs_real_3_t *restrict grad)
550 {
551 cs_halo_type_t halo_type = CS_HALO_STANDARD;
552 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
553
554 /* Does the field have a parent (variable) ?
555 Field is its own parent if not parent is specified */
556
557 const cs_field_t *parent_f = f;
558
559 const int f_parent_id
560 = cs_field_get_key_int(f, cs_field_key_id("parent_field_id"));
561 if (f_parent_id > -1)
562 parent_f = cs_field_by_id(f_parent_id);
563
564 int imrgra = cs_glob_space_disc->imrgra;
565 cs_var_cal_opt_t eqp_default = cs_parameters_var_cal_opt_default();
566
567 /* Get the calculation option from the field */
568 const cs_equation_param_t
569 *eqp = cs_field_get_equation_param_const(parent_f);
570
571 if (eqp != NULL)
572 imrgra = eqp->imrgra;
573 else
574 eqp = &eqp_default;
575
576 cs_gradient_type_by_imrgra(imrgra,
577 &gradient_type,
578 &halo_type);
579
580 int w_stride = 1;
581 cs_real_t *c_weight = NULL;
582 cs_internal_coupling_t *cpl = NULL;
583
584 if (parent_f->type & CS_FIELD_VARIABLE && eqp->idiff > 0) {
585
586 if (eqp->iwgrec == 1) {
587 /* Weighted gradient coefficients */
588 int key_id = cs_field_key_id("gradient_weighting_id");
589 int diff_id = cs_field_get_key_int(parent_f, key_id);
590 if (diff_id > -1) {
591 cs_field_t *f_weight = cs_field_by_id(diff_id);
592 c_weight = f_weight->val;
593 w_stride = f_weight->dim;
594 }
595 }
596
597 /* Internal coupling structure */
598 int key_id = cs_field_key_id_try("coupling_entity");
599 if (key_id > -1) {
600 int coupl_id = cs_field_get_key_int(parent_f, key_id);
601 if (coupl_id > -1)
602 cpl = cs_internal_coupling_by_id(coupl_id);
603 }
604
605 }
606
607 if (f->n_time_vals < 2 && use_previous_t)
608 bft_error(__FILE__, __LINE__, 0,
609 _("%s: field %s does not maintain previous time step values\n"
610 "so \"use_previous_t\" can not be handled."),
611 __func__, f->name);
612
613 cs_real_t *var = (use_previous_t) ? f->val_pre : f->val;
614
615 const cs_real_t *bc_coeff_a = NULL, *bc_coeff_b = NULL;
616 if (f->bc_coeffs != NULL) {
617 bc_coeff_a = f->bc_coeffs->a;
618 bc_coeff_b = f->bc_coeffs->b;
619 }
620
621 cs_gradient_scalar(f->name,
622 gradient_type,
623 halo_type,
624 inc,
625 recompute_cocg,
626 eqp->nswrgr,
627 0, /* ignored */
628 0, /* hyd_p_flag */
629 w_stride,
630 eqp->verbosity,
631 eqp->imligr,
632 eqp->epsrgr,
633 eqp->climgr,
634 NULL, /* f_ext */
635 bc_coeff_a,
636 bc_coeff_b,
637 var,
638 c_weight,
639 cpl, /* internal coupling */
640 grad);
641 }
642
643 /*----------------------------------------------------------------------------*/
644 /*!
645 * \brief Compute cell gradient of scalar field or component of vector or
646 * tensor field.
647 *
648 * \param[in] f pointer to field
649 * \param[in] use_previous_t should we use values from the previous
650 * time step ?
651 * \param[in] inc if 0, solve on increment; 1 otherwise
652 * \param[in] recompute_cocg should COCG FV quantities be recomputed ?
653 * \param[in] hyd_p_flag flag for hydrostatic pressure
654 * \param[in] f_ext exterior force generating
655 * the hydrostatic pressure
656 * \param[out] grad gradient
657 */
658 /*----------------------------------------------------------------------------*/
659
660 void
cs_field_gradient_potential(const cs_field_t * f,bool use_previous_t,int inc,bool recompute_cocg,int hyd_p_flag,cs_real_3_t f_ext[],cs_real_3_t * restrict grad)661 cs_field_gradient_potential(const cs_field_t *f,
662 bool use_previous_t,
663 int inc,
664 bool recompute_cocg,
665 int hyd_p_flag,
666 cs_real_3_t f_ext[],
667 cs_real_3_t *restrict grad)
668 {
669 cs_halo_type_t halo_type = CS_HALO_STANDARD;
670 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
671
672 /* Does the field have a parent (variable) ?
673 Field is its own parent if not parent is specified */
674
675 const cs_field_t *parent_f = f;
676
677 const int f_parent_id
678 = cs_field_get_key_int(f, cs_field_key_id("parent_field_id"));
679 if (f_parent_id > -1)
680 parent_f = cs_field_by_id(f_parent_id);
681
682 int imrgra = cs_glob_space_disc->imrgra;
683 cs_var_cal_opt_t eqp_default = cs_parameters_var_cal_opt_default();
684
685 /* Get the calculation option from the field */
686 const cs_equation_param_t
687 *eqp = cs_field_get_equation_param_const(parent_f);
688
689 if (eqp != NULL)
690 imrgra = eqp->imrgra;
691 else
692 eqp = &eqp_default;
693
694 cs_gradient_type_by_imrgra(imrgra,
695 &gradient_type,
696 &halo_type);
697
698 if (f->n_time_vals < 2 && use_previous_t)
699 bft_error(__FILE__, __LINE__, 0,
700 _("%s: field %s does not maintain previous time step values\n"
701 "so \"use_previous_t\" can not be handled."),
702 __func__, f->name);
703
704 int w_stride = 1;
705 cs_real_t *var = (use_previous_t) ? f->val_pre : f->val;
706
707 cs_real_t *c_weight = NULL;
708 cs_internal_coupling_t *cpl = NULL;
709
710 if (parent_f->type & CS_FIELD_VARIABLE && eqp->idiff > 0) {
711
712 if (eqp->iwgrec == 1) {
713 int key_id = cs_field_key_id("gradient_weighting_id");
714 int diff_id = cs_field_get_key_int(parent_f, key_id);
715 if (diff_id > -1) {
716 cs_field_t *f_weight = cs_field_by_id(diff_id);
717 c_weight = f_weight->val;
718 w_stride = f_weight->dim;
719 }
720 }
721
722 /* Internal coupling structure */
723 int key_id = cs_field_key_id_try("coupling_entity");
724 if (key_id > -1) {
725 int coupl_id = cs_field_get_key_int(parent_f, key_id);
726 if (coupl_id > -1)
727 cpl = cs_internal_coupling_by_id(coupl_id);
728 }
729
730 }
731
732 const cs_real_t *bc_coeff_a = NULL, *bc_coeff_b = NULL;
733 if (f->bc_coeffs != NULL) {
734 bc_coeff_a = f->bc_coeffs->a;
735 bc_coeff_b = f->bc_coeffs->b;
736 }
737
738 cs_gradient_scalar(f->name,
739 gradient_type,
740 halo_type,
741 inc,
742 recompute_cocg,
743 eqp->nswrgr,
744 0, /* tr_dim */
745 hyd_p_flag,
746 w_stride,
747 eqp->verbosity,
748 eqp->imligr,
749 eqp->epsrgr,
750 eqp->climgr,
751 f_ext,
752 bc_coeff_a,
753 bc_coeff_b,
754 var,
755 c_weight,
756 cpl, /* internal coupling */
757 grad);
758 }
759
760 /*----------------------------------------------------------------------------*/
761 /*!
762 * \brief Compute cell gradient of vector field.
763 *
764 * \param[in] f pointer to field
765 * \param[in] use_previous_t should we use values from the previous
766 * time step ?
767 * \param[in] inc if 0, solve on increment; 1 otherwise
768 * \param[out] grad gradient
769 */
770 /*----------------------------------------------------------------------------*/
771
772 void
cs_field_gradient_vector(const cs_field_t * f,bool use_previous_t,int inc,cs_real_33_t * restrict grad)773 cs_field_gradient_vector(const cs_field_t *f,
774 bool use_previous_t,
775 int inc,
776 cs_real_33_t *restrict grad)
777 {
778 cs_halo_type_t halo_type = CS_HALO_STANDARD;
779 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
780
781 int imrgra = cs_glob_space_disc->imrgra;
782 cs_var_cal_opt_t eqp_default = cs_parameters_var_cal_opt_default();
783
784 /* Get the calculation option from the field */
785 const cs_equation_param_t *eqp = cs_field_get_equation_param_const(f);
786
787 if (eqp != NULL)
788 imrgra = eqp->imrgra;
789 else
790 eqp = &eqp_default;
791
792 cs_gradient_type_by_imrgra(imrgra,
793 &gradient_type,
794 &halo_type);
795
796 cs_real_t *c_weight = NULL;
797 cs_internal_coupling_t *cpl = NULL;
798
799 if (f->type & CS_FIELD_VARIABLE && eqp->idiff > 0) {
800
801 if (eqp->iwgrec == 1) {
802 /* Weighted gradient coefficients */
803 int key_id = cs_field_key_id("gradient_weighting_id");
804 int diff_id = cs_field_get_key_int(f, key_id);
805 if (diff_id > -1) {
806 cs_field_t *f_weight = cs_field_by_id(diff_id);
807 c_weight = f_weight->val;
808 }
809 }
810
811 int key_id = cs_field_key_id_try("coupling_entity");
812 if (key_id > -1) {
813 int coupl_id = cs_field_get_key_int(f, key_id);
814 if (coupl_id > -1)
815 cpl = cs_internal_coupling_by_id(coupl_id);
816 }
817
818 }
819
820 if (f->n_time_vals < 2 && use_previous_t)
821 bft_error(__FILE__, __LINE__, 0,
822 _("%s: field %s does not maintain previous time step values\n"
823 "so \"use_previous_t\" can not be handled."),
824 __func__, f->name);
825
826 cs_real_3_t *var = (use_previous_t) ? (cs_real_3_t *)(f->val_pre)
827 : (cs_real_3_t *)(f->val);
828
829 const cs_real_3_t *bc_coeff_a = NULL;
830 const cs_real_33_t *bc_coeff_b = NULL;
831
832 if (f->bc_coeffs != NULL) {
833 int coupled_key_id = cs_field_key_id_try("coupled");
834 if (coupled_key_id > 1) {
835 if (cs_field_get_key_int(f, coupled_key_id) > 0) {
836 bc_coeff_a = (const cs_real_3_t *)f->bc_coeffs->a;
837 bc_coeff_b = (const cs_real_33_t *)f->bc_coeffs->b;
838 }
839 }
840 }
841
842 cs_gradient_vector(f->name,
843 gradient_type,
844 halo_type,
845 inc,
846 eqp->nswrgr,
847 eqp->verbosity,
848 eqp->imligr,
849 eqp->epsrgr,
850 eqp->climgr,
851 bc_coeff_a,
852 bc_coeff_b,
853 var,
854 c_weight,
855 cpl,
856 grad);
857 }
858
859 /*----------------------------------------------------------------------------*/
860 /*!
861 * \brief Compute cell gradient of tensor field.
862 *
863 * \param[in] f pointer to field
864 * \param[in] use_previous_t should we use values from the previous
865 * time step ?
866 * \param[in] inc if 0, solve on increment; 1 otherwise
867 * \param[out] grad gradient
868 */
869 /*----------------------------------------------------------------------------*/
870
871 void
cs_field_gradient_tensor(const cs_field_t * f,bool use_previous_t,int inc,cs_real_63_t * restrict grad)872 cs_field_gradient_tensor(const cs_field_t *f,
873 bool use_previous_t,
874 int inc,
875 cs_real_63_t *restrict grad)
876 {
877 cs_halo_type_t halo_type = CS_HALO_STANDARD;
878 cs_gradient_type_t gradient_type = CS_GRADIENT_GREEN_ITER;
879
880 int imrgra = cs_glob_space_disc->imrgra;
881 cs_var_cal_opt_t eqp_default = cs_parameters_var_cal_opt_default();
882
883 /* Get the calculation option from the field */
884 const cs_equation_param_t *eqp = cs_field_get_equation_param_const(f);
885
886 if (eqp != NULL)
887 imrgra = eqp->imrgra;
888 else
889 eqp = &eqp_default;
890
891 cs_gradient_type_by_imrgra(imrgra,
892 &gradient_type,
893 &halo_type);
894
895 if (f->n_time_vals < 2 && use_previous_t)
896 bft_error(__FILE__, __LINE__, 0,
897 _("%s: field %s does not maintain previous time step values\n"
898 "so \"use_previous_t\" can not be handled."),
899 __func__, f->name);
900
901 cs_real_6_t *var = (use_previous_t) ? (cs_real_6_t *)(f->val_pre)
902 : (cs_real_6_t *)(f->val);
903
904 const cs_real_6_t *bc_coeff_a = NULL;
905 const cs_real_66_t *bc_coeff_b = NULL;
906 if (f->bc_coeffs != NULL) {
907 int coupled_key_id = cs_field_key_id_try("coupled");
908 if (coupled_key_id > 1) {
909 if (cs_field_get_key_int(f, coupled_key_id) > 0) {
910 bc_coeff_a = (const cs_real_6_t *)f->bc_coeffs->a;
911 bc_coeff_b = (const cs_real_66_t *)f->bc_coeffs->b;
912 }
913 }
914 }
915
916 cs_gradient_tensor(f->name,
917 gradient_type,
918 halo_type,
919 inc,
920 eqp->nswrgr,
921 eqp->verbosity,
922 eqp->imligr,
923 eqp->epsrgr,
924 eqp->climgr,
925 bc_coeff_a,
926 bc_coeff_b,
927 var,
928 grad);
929 }
930
931 /*----------------------------------------------------------------------------*/
932 /*!
933 * \brief Interpolate field values at a given set of points.
934 *
935 * \param[in] f pointer to field
936 * \param[in] interpolation_type interpolation type
937 * \param[in] n_points number of points at which interpolation
938 * is required
939 * \param[in] point_location location of points in mesh elements
940 * (based on the field location)
941 * \param[in] point_coords point coordinates
942 * \param[out] val interpolated values
943 */
944 /*----------------------------------------------------------------------------*/
945
946 void
cs_field_interpolate(cs_field_t * f,cs_field_interpolate_t interpolation_type,cs_lnum_t n_points,const cs_lnum_t point_location[],const cs_real_3_t point_coords[],cs_real_t * val)947 cs_field_interpolate(cs_field_t *f,
948 cs_field_interpolate_t interpolation_type,
949 cs_lnum_t n_points,
950 const cs_lnum_t point_location[],
951 const cs_real_3_t point_coords[],
952 cs_real_t *val)
953 {
954 switch (interpolation_type) {
955
956 case CS_FIELD_INTERPOLATE_MEAN:
957 _field_interpolate_by_mean(f,
958 n_points,
959 point_location,
960 val);
961 break;
962
963 case CS_FIELD_INTERPOLATE_GRADIENT:
964 _field_interpolate_by_gradient(f,
965 n_points,
966 point_location,
967 point_coords,
968 val);
969 break;
970
971 default:
972 assert(0);
973 break;
974 }
975 }
976
977 /*----------------------------------------------------------------------------*/
978 /*!
979 * \brief Find local extrema of a given scalar field at each cell
980 *
981 * This assumes the field values have been synchronized.
982 *
983 * \param[in] f_id scalar field id
984 * \param[in] halo_type halo type
985 * \param[in, out] local_max local maximum value
986 * \param[in, out] local_min local minimum value
987 */
988 /*----------------------------------------------------------------------------*/
989
990 void
cs_field_local_extrema_scalar(int f_id,cs_halo_type_t halo_type,cs_real_t * local_max,cs_real_t * local_min)991 cs_field_local_extrema_scalar(int f_id,
992 cs_halo_type_t halo_type,
993 cs_real_t *local_max,
994 cs_real_t *local_min)
995 {
996 const cs_mesh_t *m = cs_glob_mesh;
997 const cs_lnum_t n_cells_ext = m->n_cells_with_ghosts;
998
999 cs_field_t *f = cs_field_by_id(f_id);
1000 const cs_real_t *restrict pvar = f->val;
1001
1002 _local_extrema_scalar(pvar,
1003 halo_type,
1004 local_max,
1005 local_min);
1006
1007 /* Initialisation of local extrema */
1008
1009 int key_scamax_id = cs_field_key_id("max_scalar");
1010 int key_scamin_id = cs_field_key_id("min_scalar");
1011
1012 cs_real_t scalar_max = cs_field_get_key_double(f, key_scamax_id);
1013 cs_real_t scalar_min = cs_field_get_key_double(f, key_scamin_id);
1014
1015 /* Bounded by the global extrema */
1016 # pragma omp parallel for
1017 for (cs_lnum_t ii = 0; ii < n_cells_ext; ii++) {
1018 local_max[ii] = CS_MIN(local_max[ii], scalar_max);
1019 local_min[ii] = CS_MAX(local_min[ii], scalar_min);
1020 }
1021 }
1022
1023 /*----------------------------------------------------------------------------*/
1024 /*!
1025 * \brief Shift field values in order to set its spatial average (fluid volume
1026 * average) to a given value.
1027 *
1028 * \param[in] f pointer to field
1029 * \param[in] va real value of fluid volume average to be set
1030 */
1031 /*----------------------------------------------------------------------------*/
1032
1033 void
cs_field_set_volume_average(cs_field_t * f,const cs_real_t va)1034 cs_field_set_volume_average(cs_field_t *f,
1035 const cs_real_t va)
1036 {
1037 const cs_lnum_t n_cells = cs_glob_mesh->n_cells;
1038
1039 cs_mesh_quantities_t *mq = cs_glob_mesh_quantities;
1040 const cs_real_t *cell_f_vol = mq->cell_f_vol;
1041 const cs_real_t tot_f_vol = mq->tot_f_vol;
1042
1043 cs_real_t *restrict val = f->val;
1044 cs_real_t p_va = 0.;
1045
1046 # pragma omp parallel for reduction(+:p_va)
1047 for (cs_lnum_t c_id = 0 ; c_id < n_cells ; c_id++) {
1048 p_va += cell_f_vol[c_id]*val[c_id];
1049 }
1050
1051 cs_parall_sum(1, CS_DOUBLE, &p_va);
1052 p_va = p_va / tot_f_vol;
1053
1054 cs_real_t shift = va - p_va;
1055
1056 # pragma omp parallel for
1057 for (cs_lnum_t c_id = 0 ; c_id < n_cells ; c_id++) {
1058 val[c_id] += shift;
1059 }
1060 }
1061
1062 /*----------------------------------------------------------------------------*/
1063 /*!
1064 * \brief Synchronize current parallel and periodic field values.
1065 *
1066 * This function currently only upates fields based on CS_MESH_LOCATION_CELLS.
1067 *
1068 * \param[in, out] f pointer to field
1069 * \param[in] halo_type halo type
1070 */
1071 /*----------------------------------------------------------------------------*/
1072
1073 void
cs_field_synchronize(cs_field_t * f,cs_halo_type_t halo_type)1074 cs_field_synchronize(cs_field_t *f,
1075 cs_halo_type_t halo_type)
1076 {
1077 if (f->location_id == CS_MESH_LOCATION_CELLS) {
1078
1079 const cs_halo_t *halo = cs_glob_mesh->halo;
1080
1081 if (halo != NULL) {
1082
1083 if (f->dim == 1)
1084 cs_halo_sync_var(halo, halo_type, f->val);
1085
1086 else {
1087
1088 cs_halo_sync_var_strided(halo, halo_type, f->val, f->dim);
1089
1090 if (cs_glob_mesh->n_init_perio > 0) {
1091 switch(f->dim) {
1092 case 9:
1093 cs_halo_perio_sync_var_tens(halo, halo_type, f->val);
1094 break;
1095 case 6:
1096 cs_halo_perio_sync_var_sym_tens(halo, halo_type, f->val);
1097 break;
1098 case 3:
1099 cs_halo_perio_sync_var_vect(halo, halo_type, f->val, 3);
1100 break;
1101 default:
1102 break;
1103 }
1104
1105 }
1106
1107 }
1108
1109 }
1110
1111 }
1112 }
1113
1114 /*----------------------------------------------------------------------------*/
1115
1116 END_C_DECLS
1117